Scippy

SCIP

Solving Constraint Integer Programs

sepa_interminor.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file sepa_interminor.c
26  * @ingroup DEFPLUGINS_SEPA
27  * @brief minor separator with intersection cuts
28  * @author Felipe Serrano
29  * @author Antonia Chmiela
30  *
31  * Let X be the matrix of auxiliary variables added for bilinear terms, X_{ij} = x_ix_j.
32  * The separator enforces quadratic constraints det(2x2 minor of X) = 0 via intersection cuts.
33  */
34 
35 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
36 
37 #include <assert.h>
38 #include <string.h>
39 
40 #include "scip/sepa_interminor.h"
41 #include "scip/expr.h"
42 #include "scip/expr_var.h"
43 #include "scip/expr_pow.h"
44 #include "scip/expr_product.h"
45 #include "scip/nlpi_ipopt.h"
46 #include "scip/cons_nonlinear.h"
47 
48 
49 
50 #define SEPA_NAME "interminor"
51 #define SEPA_DESC "intersection cuts separator to ensure that 2x2 minors of X (= xx') have determinant 0"
52 #define SEPA_PRIORITY 0
53 #define SEPA_FREQ -1
54 #define SEPA_MAXBOUNDDIST 1.0
55 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
56 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
57 
58 #define DEFAULT_MINCUTVIOL 1e-4 /**< default minimum required violation of a cut */
59 #define DEFAULT_RANDSEED 157 /**< default random seed */
60 #define DEFAULT_MAXROUNDS 10 /**< maximal number of separation rounds per node (-1: unlimited) */
61 #define DEFAULT_MAXROUNDSROOT -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
62 #define BINSEARCH_MAXITERS 120 /**< default iteration limit for binary search */
63 #define DEFAULT_USESTRENGTHENING FALSE /**< default for using strengthend intersection cuts to separate */
64 #define DEFAULT_USEBOUNDS FALSE /**< default for using nonnegativity bounds when separating */
65 
66 /*
67  * Data structures
68  */
69 
70 /** separator data */
71 struct SCIP_SepaData
72 {
73  SCIP_VAR** minors; /**< variables of 2x2 minors; each minor is stored like (auxvar_x^2,auxvar_y^2,auxvar_xy) */
74  SCIP_Bool* isdiagonal; /**< bool array determining if the variables appearing in the minor are diagonal */
75  int nminors; /**< total number of minors */
76  int minorssize; /**< size of minors array */
77  int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */
78  int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */
79  SCIP_Bool detectedminors; /**< has minor detection be called? */
80  SCIP_Real mincutviol; /**< minimum required violation of a cut */
81  SCIP_RANDNUMGEN* randnumgen; /**< random number generation */
82  SCIP_Bool usestrengthening; /**< whether to use strengthened intersection cuts to separate minors */
83  SCIP_Bool usebounds; /**< whether to also enforce nonegativity bounds of principle minors */
84 };
85 
86 /* these represent a row */
87 struct rowdata
88 {
89  int* vals; /**< index of the column */
90  int rowidx; /**< index corresponding to variable of that row */
91  int nvals; /**< number of nonzero entries in column */
92  int valssize; /**< size of the array that is currently allocated */
93  SCIP_HASHMAP* auxvars; /**< entry of the matrix */
94 };
95 
96 /*
97  * Local methods
98  */
99 
100 /** helper method to store a 2x2 minor in the separation data */
101 static
103  SCIP* scip, /**< SCIP data structure */
104  SCIP_SEPADATA* sepadata, /**< separator data */
105  SCIP_VAR* auxvarxik, /**< auxiliary variable X_ik = x_i * x_k */
106  SCIP_VAR* auxvarxil, /**< auxiliary variable X_il = x_i * x_l */
107  SCIP_VAR* auxvarxjk, /**< auxiliary variable X_jk = x_j * x_k */
108  SCIP_VAR* auxvarxjl, /**< auxiliary variable X_jl = x_j * x_l */
109  SCIP_Bool isauxvarxikdiag, /**< is X_ik diagonal? (i.e. i = k) */
110  SCIP_Bool isauxvarxildiag, /**< is X_il diagonal? (i.e. i = l) */
111  SCIP_Bool isauxvarxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */
112  SCIP_Bool isauxvarxjldiag /**< is X_jl diagonal? (i.e. j = l) */
113  )
114 {
115  assert(sepadata != NULL);
116  assert(auxvarxik != NULL);
117  assert(auxvarxil != NULL);
118  assert(auxvarxjk != NULL);
119  assert(auxvarxjl != NULL);
120  assert(auxvarxik != auxvarxil);
121  assert(auxvarxjk != auxvarxjl);
122 
123  SCIPdebugMsg(scip, "store 2x2 minor: [%s %s, %s %s]\n", SCIPvarGetName(auxvarxik), SCIPvarGetName(auxvarxil),
124  SCIPvarGetName(auxvarxjk), SCIPvarGetName(auxvarxjl));
125 
126  /* reallocate if necessary */
127  if( sepadata->minorssize < 4 * (sepadata->nminors + 1) )
128  {
129  int newsize = SCIPcalcMemGrowSize(scip, 4 * (sepadata->nminors + 1));
130  assert(newsize >= 4 * (sepadata->nminors + 1));
131 
132  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->minors), sepadata->minorssize, newsize) );
133  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->isdiagonal), sepadata->minorssize, newsize) );
134  sepadata->minorssize = newsize;
135  }
136 
137  /* store minor */
138  sepadata->minors[4 * sepadata->nminors] = auxvarxik;
139  sepadata->minors[4 * sepadata->nminors + 1] = auxvarxil;
140  sepadata->minors[4 * sepadata->nminors + 2] = auxvarxjk;
141  sepadata->minors[4 * sepadata->nminors + 3] = auxvarxjl;
142  sepadata->isdiagonal[4 * sepadata->nminors] = isauxvarxikdiag;
143  sepadata->isdiagonal[4 * sepadata->nminors + 1] = isauxvarxildiag;
144  sepadata->isdiagonal[4 * sepadata->nminors + 2] = isauxvarxjkdiag;
145  sepadata->isdiagonal[4 * sepadata->nminors + 3] = isauxvarxjldiag;
146  ++(sepadata->nminors);
147 
148  /* capture variables */
149  SCIP_CALL( SCIPcaptureVar(scip, auxvarxik) );
150  SCIP_CALL( SCIPcaptureVar(scip, auxvarxil) );
151  SCIP_CALL( SCIPcaptureVar(scip, auxvarxjk) );
152  SCIP_CALL( SCIPcaptureVar(scip, auxvarxjl) );
153 
154  return SCIP_OKAY;
155 }
156 
157 /** helper method to clear separation data */
158 static
160  SCIP* scip, /**< SCIP data structure */
161  SCIP_SEPADATA* sepadata /**< separator data */
162  )
163 {
164  int i;
165 
166  assert(sepadata != NULL);
167 
168  SCIPdebugMsg(scip, "clear separation data\n");
169 
170  /* release captured variables */
171  for( i = 0; i < 4 * sepadata->nminors; ++i )
172  {
173  assert(sepadata->minors[i] != NULL);
174  SCIP_CALL( SCIPreleaseVar(scip, &sepadata->minors[i]) );
175  }
176 
177  /* free memory */
178  SCIPfreeBlockMemoryArrayNull(scip, &sepadata->minors, sepadata->minorssize);
179  SCIPfreeBlockMemoryArrayNull(scip, &sepadata->isdiagonal, sepadata->minorssize);
180 
181  /* reset counters */
182  sepadata->nminors = 0;
183  sepadata->minorssize = 0;
184 
185  return SCIP_OKAY;
186 }
187 
188 /** helper method to get the variables associated to a minor */
189 static
191  SCIP_SEPADATA* sepadata, /**< separator data */
192  int idx, /**< index of the stored minor */
193  SCIP_VAR** auxvarxik, /**< auxiliary variable X_ik = x_i * x_k */
194  SCIP_VAR** auxvarxil, /**< auxiliary variable X_il = x_i * x_l */
195  SCIP_VAR** auxvarxjk, /**< auxiliary variable X_jk = x_j * x_k */
196  SCIP_VAR** auxvarxjl, /**< auxiliary variable X_jl = x_j * x_l */
197  SCIP_Bool* isauxvarxikdiag, /**< is X_ik diagonal? (i.e. i = k) */
198  SCIP_Bool* isauxvarxildiag, /**< is X_il diagonal? (i.e. i = l) */
199  SCIP_Bool* isauxvarxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */
200  SCIP_Bool* isauxvarxjldiag /**< is X_jl diagonal? (i.e. j = l) */
201  )
202 {
203  assert(auxvarxik != NULL);
204  assert(auxvarxil != NULL);
205  assert(auxvarxjk != NULL);
206  assert(auxvarxjl != NULL);
207 
208  *auxvarxik = sepadata->minors[4 * idx];
209  *auxvarxil = sepadata->minors[4 * idx + 1];
210  *auxvarxjk = sepadata->minors[4 * idx + 2];
211  *auxvarxjl = sepadata->minors[4 * idx + 3];
212 
213  *isauxvarxikdiag = sepadata->isdiagonal[4 * idx];
214  *isauxvarxildiag = sepadata->isdiagonal[4 * idx + 1];
215  *isauxvarxjkdiag = sepadata->isdiagonal[4 * idx + 2];
216  *isauxvarxjldiag = sepadata->isdiagonal[4 * idx + 3];
217 
218  return SCIP_OKAY;
219 }
220 
221 
222 /** adds a new entry (i.e., auxvar) of in (row, col) of matrix M.
223  *
224  * we have a matrix, M, indexed by the variables
225  * M(xi, xk) is the auxiliary variable of xi * xk if it exists
226  * We store, for each row of the matrix, the indices of the nonzero column entries (assoc with the given row) and the auxiliary variable for xi * xk
227  * The nonzero column entries are stored as an array (struct rowdata)
228  * So we have a hashmap mapping each variable (row of the matrix) with its array representing the nonzero entries of the row.
229  */
230 static
232  SCIP* scip, /**< SCIP data structure */
233  SCIP_HASHMAP* rowmap, /**< hashmap of the rows of the matrix */
234  SCIP_VAR* row, /**< variable corresponding to row of new entry */
235  SCIP_VAR* col, /**< variable corresponding to column of new entry */
236  SCIP_VAR* auxvar, /**< auxvar to insert into the matrix */
237  int* rowindices, /**< array of indices of all variables corresponding to a row */
238  int* nrows /**< number of rows */
239  )
240 {
241  SCIPdebugMsg(scip, "inserting %s in row %s and col %s \n", SCIPvarGetName(auxvar), SCIPvarGetName(row), SCIPvarGetName(col));
242 
243  /* check whether variable has an array associated to it */
244  if( SCIPhashmapExists(rowmap, (void*)row) )
245  {
246  struct rowdata* arr;
247 
248  arr = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)row);
249 
250  /* reallocate if necessary */
251  if( arr->valssize < arr->nvals + 1 )
252  {
253  int newsize = SCIPcalcMemGrowSize(scip, arr->nvals + 1);
254  assert(newsize > arr->nvals + 1);
255 
256  SCIP_CALL( SCIPreallocBufferArray(scip, &(arr->vals), newsize) );
257  arr->valssize = newsize;
258  }
259 
260  /* insert */
261  arr->vals[arr->nvals] = SCIPvarGetProbindex(col);
262  SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) );
263  arr->nvals += 1;
264  }
265  else
266  {
267  struct rowdata* arr;
268 
269  /* create index array */
270  SCIP_CALL( SCIPallocBuffer(scip, &arr) );
271  arr->valssize = 10;
272  arr->nvals = 0;
273  SCIP_CALL( SCIPallocBufferArray(scip, &arr->vals, arr->valssize) );
274  SCIP_CALL( SCIPhashmapCreate(&arr->auxvars, SCIPblkmem(scip), arr->valssize) );
275 
276  /* insert */
277  arr->rowidx = SCIPvarGetProbindex(row);
278  arr->vals[arr->nvals] = SCIPvarGetProbindex(col);
279  SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) );
280  arr->nvals += 1;
281 
282  /* store in hashmap */
283  SCIP_CALL( SCIPhashmapInsert(rowmap, (void*)row, (void *)arr) );
284 
285  /* remember the new row */
286  rowindices[*nrows] = SCIPvarGetProbindex(row);
287  *nrows += 1;
288  }
289 
290  return SCIP_OKAY;
291 }
292 
293 /** method to detect and store principal minors */
294 static
296  SCIP* scip, /**< SCIP data structure */
297  SCIP_SEPADATA* sepadata /**< separator data */
298  )
299 {
300  SCIP_CONSHDLR* conshdlr;
301  SCIP_EXPRITER* it;
302  SCIP_HASHMAP* rowmap;
303  int* rowvars = NULL;
304  int* intersection;
305  int nrowvars = 0;
306  int c;
307  int i;
308 
309 #ifdef SCIP_STATISTIC
310  SCIP_Real totaltime = -SCIPgetTotalTime(scip);
311 #endif
312 
313  assert(sepadata != NULL);
314 
315  /* check whether minor detection has been called already */
316  if( sepadata->detectedminors )
317  return SCIP_OKAY;
318 
319  assert(sepadata->minors == NULL);
320  assert(sepadata->nminors == 0);
321 
322  /* we assume that the auxiliary variables in the nonlinear constraint handler have been already generated */
323  sepadata->detectedminors = TRUE;
324 
325  /* check whether there are nonlinear constraints available */
326  conshdlr = SCIPfindConshdlr(scip, "nonlinear");
327  if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 )
328  return SCIP_OKAY;
329 
330  SCIPdebugMsg(scip, "call detectMinors()\n");
331 
332  /* allocate memory */
333  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
334  SCIP_CALL( SCIPhashmapCreate(&rowmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
335  SCIP_CALL( SCIPallocBufferArray(scip, &rowvars, SCIPgetNVars(scip)) );
336  SCIP_CALL( SCIPallocBufferArray(scip, &intersection, SCIPgetNVars(scip)) );
337 
338  /* initialize iterator */
340 
341  for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
342  {
343  SCIP_CONS* cons;
344  SCIP_EXPR* expr;
345  SCIP_EXPR* root;
346 
347  cons = SCIPconshdlrGetConss(conshdlr)[c];
348  assert(cons != NULL);
349  root = SCIPgetExprNonlinear(cons);
350  assert(root != NULL);
351 
352  for( expr = SCIPexpriterRestartDFS(it, root); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*//*lint !e440*/
353  {
354  SCIP_EXPR** children;
355  SCIP_VAR* auxvar;
356 
357  SCIPdebugMsg(scip, "visit expression %p in constraint %s\n", (void*)expr, SCIPconsGetName(cons));
358 
359  /* check whether the expression has an auxiliary variable */
360  auxvar = SCIPgetExprAuxVarNonlinear(expr);
361  if( auxvar == NULL )
362  {
363  SCIPdebugMsg(scip, "expression has no auxiliary variable -> skip\n");
364  continue;
365  }
366 
367  children = SCIPexprGetChildren(expr);
368 
369  /* check for expr = (x)^2 */
370  if( SCIPexprGetNChildren(expr) == 1 && SCIPisExprPower(scip, expr)
371  && SCIPgetExponentExprPow(expr) == 2.0
372  && SCIPgetExprAuxVarNonlinear(children[0]) != NULL )
373  {
374  SCIP_VAR* quadvar;
375 
376  assert(children[0] != NULL);
377 
378  quadvar = SCIPgetExprAuxVarNonlinear(children[0]);
379  assert(quadvar != NULL);
380 
381  SCIP_CALL( insertIndex(scip, rowmap, quadvar, quadvar, auxvar, rowvars, &nrowvars) );
382  }
383  /* check for expr = x_i * x_k */
384  else if( SCIPexprGetNChildren(expr) == 2 && SCIPisExprProduct(scip, expr)
385  && SCIPgetExprAuxVarNonlinear(children[0]) != NULL && SCIPgetExprAuxVarNonlinear(children[1]) != NULL )
386  {
387  SCIP_VAR* xi;
388  SCIP_VAR* xk;
389 
390  assert(children[0] != NULL);
391  assert(children[1] != NULL);
392 
393  xi = SCIPgetExprAuxVarNonlinear(children[0]);
394  xk = SCIPgetExprAuxVarNonlinear(children[1]);
395 
396  SCIP_CALL( insertIndex(scip, rowmap, xk, xi, auxvar, rowvars, &nrowvars) );
397  SCIP_CALL( insertIndex(scip, rowmap, xi, xk, auxvar, rowvars, &nrowvars) );
398  }
399  }
400  }
401 
402  /* sort the column entries */
403  for( i = 0; i < nrowvars; ++i )
404  {
405  struct rowdata* row;
406 
407  row = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]);
408  SCIPsortInt(row->vals, row->nvals);
409  }
410 
411  /* store 2x2 minors */
412  /* TODO: we might store some minors twice since the matrix is symmetric. Handle that! (see unit test for example) */
413  for( i = 0; i < nrowvars; ++i )
414  {
415  int j;
416  struct rowdata* rowi;
417 
418  rowi = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]);
419 
420  for( j = i + 1; j < nrowvars; ++j )
421  {
422  struct rowdata* rowj;
423  int ninter;
424 
425  rowj = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[j]]);
426 
427  SCIPcomputeArraysIntersectionInt(rowi->vals, rowi->nvals, rowj->vals, rowj->nvals, intersection, &ninter);
428 
429  if( ninter > 1)
430  {
431  int p;
432 
433  for( p = 0; p < ninter - 1; ++p )
434  {
435  int q;
436 
437  for( q = p + 1; q < ninter; ++q )
438  {
439  SCIP_HASHMAP* rowicols;
440  SCIP_HASHMAP* rowjcols;
441  SCIP_VAR* colk;
442  SCIP_VAR* coll;
443  SCIP_VAR* auxvarik;
444  SCIP_VAR* auxvaril;
445  SCIP_VAR* auxvarjk;
446  SCIP_VAR* auxvarjl;
447  int ii;
448  int jj;
449  int k;
450  int l;
451  SCIP_Bool isauxvarikdiag = FALSE;
452  SCIP_Bool isauxvarildiag = FALSE;
453  SCIP_Bool isauxvarjkdiag = FALSE;
454  SCIP_Bool isauxvarjldiag = FALSE;
455 
456  ii = rowi->rowidx;
457  jj = rowj->rowidx;
458  k = intersection[p];
459  l = intersection[q];
460 
461  rowicols = rowi->auxvars;
462  rowjcols = rowj->auxvars;
463 
464  colk = SCIPgetVars(scip)[k];
465  coll = SCIPgetVars(scip)[l];
466 
467  auxvarik = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, colk);
468  auxvaril = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, coll);
469  auxvarjk = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, colk);
470  auxvarjl = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, coll);
471 
472  if( ii == k )
473  isauxvarikdiag = TRUE;
474  else if( ii == l )
475  isauxvarildiag = TRUE;
476  if( jj == k )
477  isauxvarjkdiag = TRUE;
478  else if( jj == l )
479  isauxvarjldiag = TRUE;
480 
481  SCIP_CALL( sepadataAddMinor(scip, sepadata, auxvarik, auxvaril, auxvarjk, auxvarjl,
482  isauxvarikdiag, isauxvarildiag, isauxvarjkdiag, isauxvarjldiag) );
483  }
484  }
485  }
486  }
487  SCIPfreeBufferArrayNull(scip, &rowi->vals);
488  SCIPhashmapFree(&rowi->auxvars);
489  SCIPfreeBufferArrayNull(scip, &rowi);
490  }
491 
492 
493  SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors);
494 
495  /* free memory */
496  SCIPfreeBufferArray(scip, &intersection);
497  SCIPfreeBufferArray(scip, &rowvars);
498  SCIPhashmapFree(&rowmap);
499  SCIPfreeExpriter(&it);
500 
501 #ifdef SCIP_STATISTIC
502  totaltime += SCIPgetTotalTime(scip);
503  SCIPstatisticMessage("MINOR DETECT %s %f %d %d\n", SCIPgetProbName(scip), totaltime, sepadata->nminors, maxminors);
504 #endif
505 
506  return SCIP_OKAY;
507 }
508 
509 /** constructs map between lp position of a basic variable and its row in the tableau */
510 static
512  SCIP* scip, /**< SCIP data structure */
513  int* map /**< buffer to store the map */
514  )
515 {
516  int* basisind;
517  int nrows;
518  int i;
519 
520  nrows = SCIPgetNLPRows(scip);
521  SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
522 
523  SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
524  for( i = 0; i < nrows; ++i )
525  {
526  if( basisind[i] >= 0 )
527  map[basisind[i]] = i;
528  }
529 
530  SCIPfreeBufferArray(scip, &basisind);
531 
532  return SCIP_OKAY;
533 }
534 
535 /** The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
536  * SQRT(A t^2 + B t + C) - (D t + E).
537  * This function computes the coefficients A, B, C, D, E for the given ray.
538  */
539 static
541  SCIP* scip, /**< SCIP data structure */
542  SCIP_Real* ray, /**< coefficients of ray */
543  SCIP_VAR** vars, /**< variables */
544  SCIP_Real* coefs, /**< buffer to store A, B, C, D, and E of cases 1, 2, 3, or 4a*/
545  SCIP_Real* coefs4b, /**< buffer to store A, B, C, D, and E of case 4b */
546  SCIP_Real* coefscondition, /**< buffer to store coefs for checking whether we are in case 4a or 4b */
547  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
548  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
549  SCIP_Bool* success /**< FALSE if we need to abort generation because of numerics */
550  )
551 {
552  SCIP_Real eigenvectors[16] = {1.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0};
553  SCIP_Real eigenvalues[4] = {0.5, 0.5, -0.5, -0.5};
554  SCIP_Real eigencoef = 0.7071067811865475244008443621048490;
555  SCIP_Real* a;
556  SCIP_Real* b;
557  SCIP_Real* c;
558  SCIP_Real* d;
559  SCIP_Real* e;
560  SCIP_Real min;
561  SCIP_Real max;
562  SCIP_Real norm1;
563  SCIP_Real norm2;
564  int negidx;
565  int posidx;
566  int i;
567 
568  *success = TRUE;
569 
570  /* set all coefficients to zero */
571  memset(coefs, 0, 5 * sizeof(SCIP_Real));
572  memset(coefs4b, 0, 5 * sizeof(SCIP_Real));
573  norm1 = 0.0;
574  norm2 = 0.0;
575 
576  a = coefs;
577  b = coefs + 1;
578  c = coefs + 2;
579  d = coefs + 3;
580  e = coefs + 4;
581 
582  negidx = 2;
583  posidx = 0;
584  for( i = 0; i < 4; ++i )
585  {
586  int j;
587  SCIP_Real vzlp;
588  SCIP_Real vdotray;
589 
590  vzlp = 0;
591  vdotray = 0;
592 
593  /* compute eigenvec * ray and eigenvec * solution */
594  for( j = 0; j < 4; ++j )
595  {
596  vdotray += eigencoef * eigenvectors[4 * i + j] * ray[j];
597  vzlp += eigencoef * eigenvectors[4 * i + j] * SCIPvarGetLPSol(vars[j]);
598  }
599 
600  if( eigenvalues[i] > 0 )
601  {
602  /* positive eigenvalue: compute D and E */
603  *d += eigenvalues[i] * vzlp * vdotray;
604  *e += eigenvalues[i] * SQR( vzlp );
605 
606  if( usebounds )
607  {
608  norm1 += eigenvalues[i] * (1 - SQR( ad[posidx] )) * SQR( vzlp );
609  norm2 += SQRT( eigenvalues[i] ) * ad[posidx] * vzlp;
610  ++posidx;
611  }
612 
613  }
614  else
615  {
616  /* negative eigenvalue: compute A, B, and C */
617  *a -= eigenvalues[i] * SQR( vdotray );
618  *b -= 2.0 * eigenvalues[i] * vzlp * vdotray;
619  *c -= eigenvalues[i] * SQR( vzlp );
620 
621  if( usebounds )
622  {
623  coefs4b[0] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vdotray );
624  coefs4b[1] -= 2.0 * eigenvalues[i] * (1 - SQR( ad[negidx] )) * vzlp * vdotray;
625  coefs4b[2] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vzlp );
626  coefs4b[3] += SQRT( -eigenvalues[i] ) * ad[negidx] * vdotray;
627  coefs4b[4] += SQRT( -eigenvalues[i] ) * ad[negidx] * vzlp;
628  ++negidx;
629  }
630  }
631  }
632 
633  assert(*e > 0);
634 
635  if( SQRT( *c ) - SQRT( *e ) >= 0.0 )
636  {
637  assert(SQRT( *c ) - SQRT( *e ) < 1e-6);
638  *success = FALSE;
639  return SCIP_OKAY;
640  }
641 
642  /* finish computation of coefficients when using bounds */
643  if( usebounds )
644  {
645  coefscondition[0] = norm2 / SQRT( *e );
646  coefscondition[1] = coefs4b[3];
647  coefscondition[2] = coefs4b[4];
648 
649  coefs4b[0] *= norm1 / *e;
650  coefs4b[1] *= norm1 / *e;
651  coefs4b[2] *= norm1 / *e;
652  coefs4b[3] *= norm2 / SQRT( *e );
653  coefs4b[4] *= norm2 / SQRT( *e );
654 
655  coefs4b[3] += *d / SQRT( *e );
656  coefs4b[4] += SQRT( *e );
657 
658  assert( SQRT( coefs4b[2] ) - coefs4b[4] < 0.0 );
659  }
660 
661  /* finish computation of D and E */
662  *e = SQRT( *e );
663  *d /= *e;
664 
665  /* maybe we want to avoid a large dynamism between A, B and C */
666  max = 0.0;
667  min = SCIPinfinity(scip);
668  for( i = 0; i < 3; ++i )
669  {
670  SCIP_Real absval;
671 
672  absval = ABS(coefs[i]);
673  if( max < absval )
674  max = absval;
675  if( absval != 0.0 && absval < min )
676  min = absval;
677  }
678 
679  if( SCIPisHugeValue(scip, max / min) )
680  {
681 #ifdef DEBUG_INTERSECTIONCUT
682  printf("Bad numerics: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min);
683 #endif
684  *success = FALSE;
685  return SCIP_OKAY;
686  }
687 
688  /* some sanity checks */
689  assert(*c >= 0); /* radicand at zero */
690  assert(SQRT( *c ) - *e < 0); /* the function at 0 must be negative */
691  assert(*a >= 0); /* the function inside the root is convex */
692 
693 #ifdef DEBUG_INTERSECTIONCUT
694  SCIPinfoMessage(scip, NULL, "Restriction yields: a,b,c,d,e %g %g %g %g %g\n", coefs[0], coefs[1], coefs[2], coefs[3], coefs[4]);
695 #endif
696 
697  return SCIP_OKAY;
698 }
699 
700 /** returns phi(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) */ /*lint -e{715}*/
701 static
703  SCIP* scip, /**< SCIP data structure */
704  SCIP_Real t, /**< argument of phi restricted to ray */
705  SCIP_Real a, /**< value of A */
706  SCIP_Real b, /**< value of B */
707  SCIP_Real c, /**< value of C */
708  SCIP_Real d, /**< value of D */
709  SCIP_Real e /**< value of E */
710  )
711 {
712 #ifdef INTERCUTS_DBLDBL
713  SCIP_Real QUAD(lin);
714  SCIP_Real QUAD(disc);
715  SCIP_Real QUAD(tmp);
716  SCIP_Real QUAD(root);
717 
718  /* d * t + e */
719  SCIPquadprecProdDD(lin, d, t);
720  SCIPquadprecSumQD(lin, lin, e);
721 
722  /* a * t * t */
723  SCIPquadprecSquareD(disc, t);
724  SCIPquadprecProdQD(disc, disc, a);
725 
726  /* b * t */
727  SCIPquadprecProdDD(tmp, b, t);
728 
729  /* a * t * t + b * t */
730  SCIPquadprecSumQQ(disc, disc, tmp);
731 
732  /* a * t * t + b * t + c */
733  SCIPquadprecSumQD(disc, disc, c);
734 
735  /* sqrt(above): can't take sqrt of 0! */
736  if( QUAD_TO_DBL(disc) == 0 )
737  {
738  QUAD_ASSIGN(root, 0.0);
739  }
740  else
741  {
742  SCIPquadprecSqrtQ(root, disc);
743  }
744 
745  /* final result */
746  QUAD_SCALE(lin, -1.0);
747  SCIPquadprecSumQQ(tmp, root, lin);
748 
749  assert(!SCIPisInfinity(scip, t) || QUAD_TO_DBL(tmp) <= 0);
750 
751  return QUAD_TO_DBL(tmp);
752 #else
753  return SQRT( a * t * t + b * t + c ) - ( d * t + e );
754 #endif
755 }
756 
757 /** helper function of computeRoot: we want phi to be <= 0 */
758 static
760  SCIP* scip, /**< SCIP data structure */
761  SCIP_Real a, /**< value of A */
762  SCIP_Real b, /**< value of B */
763  SCIP_Real c, /**< value of C */
764  SCIP_Real d, /**< value of D */
765  SCIP_Real e, /**< value of E */
766  SCIP_Real* sol /**< buffer to store solution; also gives initial point */
767  )
768 {
769  SCIP_Real lb = 0.0;
770  SCIP_Real ub = *sol;
771  SCIP_Real curr;
772  int i;
773 
774  for( i = 0; i < BINSEARCH_MAXITERS; ++i )
775  {
776  SCIP_Real phival;
777 
778  curr = (lb + ub) / 2.0;
779  phival = evalPhiAtRay(scip, curr, a, b, c, d, e);
780 #ifdef INTERCUT_MOREDEBUG
781  printf("%d: lb,ub %.10f, %.10f. curr = %g -> phi at curr %g -> phi at lb %g \n", i, lb, ub, curr, phival, evalPhiAtRay(scip, lb, a, b, c, d, e));
782 #endif
783 
784  if( phival <= 0.0 )
785  {
786  lb = curr;
787  if( SCIPisFeasZero(scip, phival) || SCIPisFeasEQ(scip, ub, lb) )
788  break;
789  }
790  else
791  ub = curr;
792  }
793 
794  *sol = lb;
795 
796 }
797 
798 /** checks if we are in case 4a, i.e., if
799  * (num(xhat_{r+1}(zlp)) / E) * SQRT(A * tsol^2 + B * tsol + C) + w(ray) * tsol + num(yhat_{s+1}(zlp)) <= 0
800  */
801 static
803  SCIP_Real tsol, /**< t in the above formula */
804  SCIP_Real* coefs, /**< coefficients A, B, C, D, and E of case 4a */
805  SCIP_Real* coefscondition /**< extra coefficients needed for the evaluation of the condition:
806  num(xhat_{r+1}(zlp)) / E; w(ray); num(yhat_{s+1}(zlp)) */
807  )
808 {
809  return (coefscondition[0] * SQRT( coefs[0] * SQR( tsol ) + coefs[1] * tsol + coefs[2] ) + coefscondition[1] *
810  tsol + coefscondition[2]) <= 0.0;
811 }
812 
813 /** finds smallest positive root phi by finding the smallest positive root of
814  * (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) = 0
815  *
816  * However, we are conservative and want a solution such that phi is negative, but close to 0;
817  * thus we correct the result with a binary search
818  */
819 static
821  SCIP* scip, /**< SCIP data structure */
822  SCIP_Real* coefs /**< value of A */
823  )
824 {
825  SCIP_Real sol;
826  SCIP_INTERVAL bounds;
827  SCIP_INTERVAL result;
828  SCIP_Real a = coefs[0];
829  SCIP_Real b = coefs[1];
830  SCIP_Real c = coefs[2];
831  SCIP_Real d = coefs[3];
832  SCIP_Real e = coefs[4];
833 
834  /* there is an intersection point if and only if SQRT(A) > D: here we are beliving in math, this might cause
835  * numerical issues
836  */
837  if( SQRT( a ) <= d )
838  {
839  sol = SCIPinfinity(scip);
840 
841  return sol;
842  }
843 
844  SCIPintervalSetBounds(&bounds, 0.0, SCIPinfinity(scip));
845 
846  /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds.
847  * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus
848  * tsol is the smallest t such that (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) >= 0
849  */
851  e, -(c - e * e), bounds);
852 
853  /* it can still be empty because of our infinity, I guess... */
855 
856  /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary
857  * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the
858  * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g.,
859  * ex8_3_1, bchoco05, etc
860  */
861  if( evalPhiAtRay(scip, sol, a, b, c, d, e) <= 1e-10 )
862  {
863 #ifdef INTERCUT_MOREDEBUG
864  printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e));
865  printf("don't do bin search\n");
866 #endif
867 
868  return sol;
869  }
870  else
871  {
872  /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */
873 #ifdef INTERCUT_MOREDEBUG
874  printf("do bin search because phival is %g\n", evalPhiAtRay(scip, sol, a, b, c, d, e));
875 #endif
876  doBinarySearch(scip, a, b, c, d, e, &sol);
877  }
878 
879  return sol;
880 }
881 
882 /** The maximal S-free set is gamma(z) <= 0; we find the intersection point of the ray `ray` starting from zlp with the
883  * boundary of the S-free set.
884  * That is, we find t >= 0 such that gamma(zlp + t * ray) = 0.
885  *
886  * In cases 1,2, and 3, gamma is of the form
887  * gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E)
888  *
889  * In the case 4 gamma is of the form
890  * gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) if some condition holds
891  * SQRT(A' t^2 + B' t + C') - (D' t + E') otherwise
892  *
893  * It can be shown (given the special properties of gamma) that the smallest positive root of each function of the form
894  * SQRT(a t^2 + b t + c) - (d t + e)
895  * is the same as the smallest positive root of the quadratic equation:
896  * (SQRT(a t^2 + b t + c) - (d t + e)) * (SQRT(a t^2 + b t + c) + (d t + e)) = 0
897  * <==> (a - d^2) t^2 + (b - 2 d*e) t + (c - e^2) = 0
898  *
899  * So, in cases 1, 2, and 3, this function just returns the solution of the above equation.
900  * In case 4, it first solves the equation assuming we are in the first piece.
901  * If there is no solution, then the second piece can't have a solution (first piece >= second piece for all t)
902  * Then we check if the solution satisfies the condition.
903  * If it doesn't then we solve the equation for the second piece.
904  * If it has a solution, then it _has_ to be the solution.
905  */
906 static
908  SCIP* scip, /**< SCIP data structure */
909  SCIP_Bool usebounds, /**< whether we are in case 4 or not */
910  SCIP_Real* coefs, /**< values of A, B, C, D, and E of cases 1, 2, 3, or 4a */
911  SCIP_Real* coefs4b, /**< values of A, B, C, D, and E of case 4b */
912  SCIP_Real* coefscondition /**< values needed to evaluate condition of case 4 */
913  )
914 {
915  SCIP_Real sol;
916  SCIP_Real sol4b;
917 
918  assert(coefs != NULL);
919 
920  if( ! usebounds )
921  return computeRoot(scip, coefs);
922 
923  assert(coefs4b != NULL);
924  assert(coefscondition != NULL);
925 
926  /* compute solution of first piece */
927  sol = computeRoot(scip, coefs);
928 
929  /* if there is no solution --> second piece doesn't have solution */
930  if( SCIPisInfinity(scip, sol) )
931  {
932  /* this assert fails on multiplants_mtg5 the problem is that sqrt(A) <= D in 4a but not in 4b,
933  * now, this is impossible since the phi4a >= phi4b, so actually sqrt(A) is 10e-15 away from
934  * D in 4b
935  */
936  /* assert(SCIPisInfinity(scip, computeRoot(scip, coefs4b))); */
937  return sol;
938  }
939 
940  /* if solution of 4a is in 4a, then return */
941  if( isCase4a(sol, coefs, coefscondition) )
942  return sol;
943 
944  /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should
945  * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function
946  * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the
947  * binary search) we can find a slightly smaller solution; thus, we just keep the larger one
948  */
949  sol4b = computeRoot(scip, coefs4b);
950 
951  return MAX(sol, sol4b);
952 }
953 
954 /** adds cutcoef * (col - col*) to rowprep */
955 static
957  SCIP* scip, /**< SCIP data structure */
958  SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */
959  SCIP_Real cutcoef, /**< cut coefficient */
960  SCIP_COL* col /**< column to add to rowprep */
961  )
962 {
963  assert(col != NULL);
964 
965 #ifdef DEBUG_INTERCUTS_NUMERICS
966  SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)),
968  SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" :
969  "upper bound" , SCIPcolGetPrimsol(col));
970 #endif
971 
972  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(col), cutcoef) );
973  SCIProwprepAddConstant(rowprep, -cutcoef * SCIPcolGetPrimsol(col) );
974 
975  return SCIP_OKAY;
976 }
977 
978 /** adds cutcoef * (slack - slack*) to rowprep
979  *
980  * row is lhs <= <coefs, vars> + constant <= rhs, thus slack is defined by
981  * slack + <coefs, vars> + constant = side
982  * If row (slack) is at upper, it means that <coefs,vars*> + constant = rhs, and so
983  * slack* = side - rhs --> slack - slack* = rhs - <coefs, vars> - constant.
984  * If row (slack) is at lower, then <coefs,vars*> + constant = lhs, and so
985  * slack* = side - lhs --> slack - slack* = lhs - <coefs, vars> - constant.
986  */
987 static
989  SCIP* scip, /**< SCIP data structure */
990  SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */
991  SCIP_Real cutcoef, /**< cut coefficient */
992  SCIP_ROW* row, /**< row, whose slack we are adding to rowprep */
993  SCIP_Bool* success /**< buffer to store whether the row is nonbasic enough */
994  )
995 {
996  int i;
997  SCIP_COL** rowcols;
998  SCIP_Real* rowcoefs;
999  int nnonz;
1000 
1001  assert(row != NULL);
1002 
1003  rowcols = SCIProwGetCols(row);
1004  rowcoefs = SCIProwGetVals(row);
1005  nnonz = SCIProwGetNLPNonz(row);
1006 
1007 #ifdef DEBUG_INTERCUTS_NUMERICS
1008  SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row));
1009  SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" :
1011  SCIPgetRowActivity(scip, row));
1012 #endif
1013 
1015  {
1016  assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row)));
1017  if( ! SCIPisFeasEQ(scip, SCIProwGetLhs(row), SCIPgetRowActivity(scip, row)) )
1018  {
1019  *success = FALSE;
1020  return SCIP_OKAY;
1021  }
1022 
1023  SCIProwprepAddConstant(rowprep, SCIProwGetLhs(row) * cutcoef);
1024  }
1025  else
1026  {
1027  assert(!SCIPisInfinity(scip, SCIProwGetRhs(row)));
1028  if( ! SCIPisFeasEQ(scip, SCIProwGetRhs(row), SCIPgetRowActivity(scip, row)) )
1029  {
1030  *success = FALSE;
1031  return SCIP_OKAY;
1032  }
1033 
1034  SCIProwprepAddConstant(rowprep, SCIProwGetRhs(row) * cutcoef);
1035  }
1036 
1037  for( i = 0; i < nnonz; i++ )
1038  {
1039  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) );
1040  }
1041 
1042  SCIProwprepAddConstant(rowprep, -SCIProwGetConstant(row) * cutcoef);
1043 
1044  return SCIP_OKAY;
1045 }
1046 
1047 /** get the tableau rows of the variables in vars */
1048 static
1050  SCIP* scip, /**< SCIP data structure */
1051  SCIP_VAR** vars, /**< variables in the minor */
1052  int* basicvarpos2tableaurow,/**< map between basic var and its tableau row */
1053  SCIP_HASHMAP* tableau, /**< map between var an its tableau row */
1054  SCIP_Real** tableaurows, /**< buffer to store tableau row */
1055  SCIP_Bool* success /**< set to TRUE if no variable had basisstat = ZERO */
1056  )
1057 {
1058  int v;
1059  int nrows;
1060  int ncols;
1061 
1062  *success = TRUE;
1063 
1064  nrows = SCIPgetNLPRows(scip);
1065  ncols = SCIPgetNLPCols(scip);
1066 
1067  /* check if we have the tableau row of the variable and if not compute it */
1068  for( v = 0; v < 4; ++v )
1069  {
1070  if( ! SCIPhashmapExists(tableau, (void*)vars[v]) )
1071  {
1072  SCIP_COL* col;
1073 
1074  /* get column of variable */
1075  col = SCIPvarGetCol(vars[v]);
1076 
1077  /* if variable is basic, then get its tableau row and insert it in the hashmap */
1079  {
1080  int lppos;
1081  SCIP_Real* densetableaurow;
1082 
1083  lppos = SCIPcolGetLPPos(col);
1084  SCIP_CALL( SCIPallocBufferArray(scip, &densetableaurow, ncols + nrows) );
1085 
1086  SCIP_CALL( SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], NULL, NULL) );
1087  SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], densetableaurow, NULL, NULL) );
1088 
1089  /* insert tableau row in hashmap*/
1090  SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)densetableaurow) );
1091 
1092  }
1093  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1094  {
1095  *success = FALSE;
1096  return SCIP_OKAY; /* don't even bother */
1097  }
1098  else
1099  {
1100  SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)NULL) );
1101  }
1102 
1103  }
1104 
1105  /* get tableau row of var */
1106  tableaurows[v] = (SCIP_Real *)SCIPhashmapGetImage(tableau, (void*)vars[v]);
1107  }
1108  return SCIP_OKAY;
1109 }
1110 
1111 /** computes the cut coefs of the non-basic (non-slack) variables (correspond to cols) and adds them to the
1112  * intersection cut
1113  */
1114 static
1116  SCIP* scip, /**< SCIP data structure */
1117  SCIP_VAR** vars, /**< variables */
1118  SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */
1119  SCIP_ROWPREP* rowprep, /**< store cut */
1120  SCIP_Real* rays, /**< buffer to store rays */
1121  int* nrays, /**< pointer to store number of nonzero rays */
1122  int* rayslppos, /**< buffer to store lppos of nonzero rays */
1123  SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */
1124  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1125  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1126  SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */
1127  )
1128 {
1129  int i;
1130  int ncols;
1131  SCIP_COL** cols;
1132 
1133  *success = TRUE;
1134 
1135  /* loop over non-basic (non-slack) variables */
1136  cols = SCIPgetLPCols(scip);
1137  ncols = SCIPgetNLPCols(scip);
1138  for( i = 0; i < ncols; ++i )
1139  {
1140  SCIP_COL* col;
1141  SCIP_Real coefs[5];
1142  SCIP_Real coefs4b[5];
1143  SCIP_Real coefscondition[3];
1144  SCIP_Real factor;
1145  SCIP_Bool israynonzero;
1146  SCIP_Real cutcoef;
1147  SCIP_Real interpoint;
1148  int v;
1149 
1150  col = cols[i];
1151 
1152  /* set factor to store entries of ray as = [-BinvL, BinvU] */
1154  factor = -1.0;
1155  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
1156  factor = 1.0;
1157  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1158  {
1159  *success = FALSE;
1160  return SCIP_OKAY;
1161  }
1162  else
1163  continue;
1164 
1165  /* build the ray */
1166  israynonzero = FALSE;
1167  for( v = 0; v < 4; ++v )
1168  {
1169  if( tableaurows[v] != NULL )
1170  rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][i]) ? 0.0 : tableaurows[v][i]);
1171  else
1172  {
1173  if( col == SCIPvarGetCol(vars[v]) )
1174  rays[(*nrays) * 4 + v] = -factor;
1175  else
1176  rays[(*nrays) * 4 + v] = 0.0;
1177  }
1178 
1179  israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1180  }
1181 
1182  /* do nothing if ray is 0 */
1183  if( ! israynonzero )
1184  continue;
1185 
1186  /* compute the cut */
1187  SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1188  ad, success) );
1189 
1190  if( *success == FALSE )
1191  return SCIP_OKAY;
1192 
1193  /* compute intersection point */
1194  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1195 
1196  /* store intersection points */
1197  interpoints[*nrays] = interpoint;
1198 
1199  /* remember lppos */
1200  rayslppos[*nrays] = i;
1201 
1202  /* count nonzero rays */
1203  *nrays += 1;
1204 
1205  /* compute cut coef */
1206  cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1207 
1208  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1210  SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ? -cutcoef :
1211  cutcoef, col) );
1212  }
1213 
1214  return SCIP_OKAY;
1215 }
1216 
1217 /** computes the cut coefs of the non-basic slack variables (correspond to rows) and adds them to the
1218  * intersection cut
1219  */
1220 static
1222  SCIP* scip, /**< SCIP data structure */
1223  SCIP_VAR** vars, /**< variables */
1224  SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */
1225  SCIP_ROWPREP* rowprep, /**< store cut */
1226  SCIP_Real* rays, /**< buffer to store rays */
1227  int* nrays, /**< pointer to store number of nonzero rays */
1228  int* rayslppos, /**< buffer to store lppos of nonzero rays */
1229  SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */
1230  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1231  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1232  SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */
1233  )
1234 {
1235  int i;
1236  int nrows;
1237  int ncols;
1238  SCIP_ROW** rows;
1239 
1240  nrows = SCIPgetNLPRows(scip);
1241  ncols = SCIPgetNLPCols(scip);
1242 
1243  *success = TRUE;
1244 
1245  /* loop over non-basic slack variables */
1246  rows = SCIPgetLPRows(scip);
1247  for( i = 0; i < nrows; ++i )
1248  {
1249  SCIP_ROW* row;
1250  SCIP_Real coefs[5];
1251  SCIP_Real coefs4b[5];
1252  SCIP_Real coefscondition[3];
1253  SCIP_Real factor;
1254  SCIP_Bool israynonzero;
1255  SCIP_Real cutcoef;
1256  SCIP_Real interpoint;
1257  int v;
1258 
1259  row = rows[i];
1260 
1261  /* set factor to store entries of ray as = [BinvL, -BinvU] */
1263  factor = 1.0;
1264  else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER )
1265  factor = -1.0;
1266  else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_ZERO )
1267  {
1268  *success = FALSE;
1269  return SCIP_OKAY;
1270  }
1271  else
1272  continue;
1273 
1274  /* build the ray */
1275  israynonzero = FALSE;
1276  for( v = 0; v < 4; ++v )
1277  {
1278  int idx;
1279 
1280  idx = ncols + i;
1281 
1282  if( tableaurows[v] != NULL )
1283  rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][idx]) ? 0.0 : tableaurows[v][idx]);
1284  else
1285  {
1286  /* TODO: We assume that slack variables can never occure in the minor. This is correct, right? */
1287  rays[(*nrays) * 4 + v] = 0.0;
1288  }
1289 
1290  israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1291  }
1292 
1293  /* do nothing if ray is 0 */
1294  if( ! israynonzero )
1295  continue;
1296 
1297  /* compute the cut */
1298  SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1299  ad, success) );
1300 
1301  if( *success == FALSE )
1302  return SCIP_OKAY;
1303 
1304  /* compute intersection point */
1305  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1306 
1307  /* store intersection points */
1308  interpoints[*nrays] = interpoint;
1309 
1310  /* store lppos of ray, make it negative so we can differentiate between cols and rows */
1311  rayslppos[*nrays] = -i - 1;
1312 
1313  /* count nonzero rays */
1314  *nrays += 1;
1315 
1316  /* compute cut coef */
1317  cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1318 
1319  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1321 
1322  SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ? cutcoef :
1323  -cutcoef, row, success) ); /* rows have flipper base status! */
1324  }
1325 
1326  return SCIP_OKAY;
1327 }
1328 
1329 /* checks if two rays are linearly dependent */
1330 static
1332  SCIP* scip, /**< SCIP data structure */
1333  SCIP_Real* ray1, /**< coefficients of ray 1 */
1334  SCIP_Real* ray2, /**< coefficients of ray 2 */
1335  SCIP_Real* coef /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are
1336  dependent */
1337  )
1338 {
1339  int i;
1340 
1341  *coef = 0.0;
1342 
1343  for( i = 0; i < 4; ++i )
1344  {
1345  /* rays cannot be dependent if one ray has zero entry and the other one doesn't */
1346  if( (SCIPisZero(scip, ray1[i]) && ! SCIPisZero(scip, ray2[i])) ||
1347  (! SCIPisZero(scip, ray1[i]) && SCIPisZero(scip, ray2[i])) )
1348  {
1349  return FALSE;
1350  }
1351 
1352  if( *coef != 0.0 )
1353  {
1354  /* cannot be dependent if the coefs aren't equal for all entries */
1355  if( ! SCIPisFeasEQ(scip, *coef, ray1[i] / ray2[i]) )
1356  return FALSE;
1357  }
1358  else
1359  *coef = ray1[i] / ray2[i];
1360  }
1361 
1362  return TRUE;
1363 }
1364 
1365 /** finds the smallest negative steplength for the current ray r_idx such that the combination
1366  * of r_idx with all rays not in the recession cone is in the recession cone
1367  */
1368 static
1370  SCIP* scip, /**< SCIP data structure */
1371  SCIP_Real* rays, /**< rays */
1372  int nrays, /**< number of nonzero rays */
1373  int idx, /**< index of current ray we want to find rho for */
1374  SCIP_Real* interpoints, /**< intersection points of nonzero rays */
1375  SCIP_VAR** vars, /**< variables */
1376  SCIP_Real* rho, /**< pointer to store the optimal rho */
1377  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1378  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1379  SCIP_Bool* success /**< TRUE if computation of rho was successful */
1380  )
1381 {
1382  int i;
1383 
1384  *success = TRUE;
1385 
1386  /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
1387  * smallest of them is then the steplength rho we use for the current ray */
1388  *rho = 0;
1389  for( i = 0; i < nrays; ++i )
1390  {
1391  SCIP_Real currentrho;
1392  SCIP_Real coef;
1393 
1394  if( SCIPisInfinity(scip, interpoints[i]) )
1395  continue;
1396 
1397  /* if the rays are linearly independent, we don't need to search for rho */
1398  if( raysAreDependent(scip, &rays[4 * i], &rays[4 * idx], &coef) )
1399  currentrho = coef * interpoints[i];
1400  else
1401  {
1402  SCIP_Real lb;
1403  SCIP_Real ub;
1404  SCIP_Real alpha;
1405  int j;
1406 
1407  /* do binary search by lookig at the convex combinations of r_i and r_j */
1408  lb = 0.0;
1409  ub = 1.0;
1410 
1411  for( j = 0; j < BINSEARCH_MAXITERS; ++j )
1412  {
1413  SCIP_Real coefs[5];
1414  SCIP_Real coefs4b[5];
1415  SCIP_Real coefscondition[3];
1416  SCIP_Real newray[4];
1417  SCIP_Real interpoint;
1418  int k;
1419 
1420  alpha = (lb + ub) / 2.0;
1421 
1422  /* build the ray alpha * ray_i + (1 - alpha) * ray_idx */
1423  for( k = 0; k < 4; ++k )
1424  newray[k] = alpha * rays[4 * i + k] - (1 - alpha) * rays[4 * idx + k];
1425 
1426  /* restrict phi to the "new" ray */
1427  SCIP_CALL( computeRestrictionToRay(scip, newray, vars, coefs, coefs4b, coefscondition, usebounds,
1428  ad, success) );
1429 
1430  if( ! *success )
1431  return SCIP_OKAY;
1432 
1433  /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
1434  * positive
1435  */
1436 
1437  /* compute intersection point */
1438  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1439 
1440  /* no root exists */
1441  if( SCIPisInfinity(scip, interpoint) )
1442  {
1443  lb = alpha;
1444  if( SCIPisEQ(scip, ub, lb) )
1445  break;
1446  }
1447  else
1448  ub = alpha;
1449  }
1450 
1451  /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
1452  * cannot move the ray in the recession cone, i.e. strengthening is not possible */
1453  if( SCIPisZero(scip, alpha) )
1454  {
1455  *rho = -SCIPinfinity(scip);
1456  return SCIP_OKAY;
1457  }
1458  else
1459  currentrho = (alpha - 1) * interpoints[i] / alpha;
1460  }
1461 
1462  if( currentrho < *rho )
1463  *rho = currentrho;
1464  }
1465 
1466  return SCIP_OKAY;
1467 }
1468 
1469 /** computes negative steplengths for the rays that are in the recession cone of the S-free set, i.e.,
1470  * which have an infinite intersection point.
1471  */
1472 static
1474  SCIP* scip, /**< SCIP data structure */
1475  SCIP_VAR** vars, /**< variables */
1476  SCIP_Real* rays, /**< rays */
1477  int nrays, /**< number of nonzero rays */
1478  int* rayslppos, /**< lppos of nonzero rays */
1479  SCIP_Real* interpoints, /**< intersection points */
1480  SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */
1481  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1482  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1483  SCIP_Bool* success /**< if a cut candidate could be computed */
1484  )
1485 {
1486  SCIP_COL** cols;
1487  SCIP_ROW** rows;
1488  int i;
1489 
1490  *success = TRUE;
1491 
1492  cols = SCIPgetLPCols(scip);
1493  rows = SCIPgetLPRows(scip);
1494 
1495  /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
1496  * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
1497  for( i = 0; i < nrays ; ++i )
1498  {
1499  SCIP_Real rho;
1500  SCIP_Real cutcoef;
1501  int lppos;
1502 
1503  if( !SCIPisInfinity(scip, interpoints[i]) )
1504  continue;
1505 
1506  /* compute the smallest rho */
1507  SCIP_CALL( findRho(scip, rays, nrays, i, interpoints, vars, &rho, usebounds, ad, success) );
1508 
1509  if( ! *success )
1510  continue;
1511 
1512  /* compute cut coef */
1513  cutcoef = SCIPisInfinity(scip, -rho) ? 0.0 : 1.0 / rho;
1514 
1515  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1516  lppos = rayslppos[i];
1517  if( lppos < 0 )
1518  {
1519  lppos = -lppos - 1;
1520 
1521  assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
1523 
1524  SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
1525  -cutcoef, rows[lppos], success) ); /* rows have flipped base status! */
1526 
1527  if( ! *success )
1528  return SCIP_OKAY;
1529  }
1530  else
1531  {
1532  assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
1534  SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
1535  cutcoef, cols[lppos]) );
1536  }
1537  }
1538 
1539  return SCIP_OKAY;
1540 }
1541 
1542 /** separates cuts for stored principal minors */
1543 static
1545  SCIP* scip, /**< SCIP data structure */
1546  SCIP_SEPA* sepa, /**< separator */
1547  SCIP_SEPADATA* sepadata, /**< separator data */
1548  SCIP_VAR* xik, /**< variable X_ik = x_i * x_k */
1549  SCIP_VAR* xil, /**< variable X_il = x_i * x_l */
1550  SCIP_VAR* xjk, /**< variable X_jk = x_j * x_k */
1551  SCIP_VAR* xjl, /**< variable X_jl = x_j * x_l */
1552  SCIP_Bool* isxikdiag, /**< is X_ik diagonal? (i.e. i = k) */
1553  SCIP_Bool* isxildiag, /**< is X_il diagonal? (i.e. i = l) */
1554  SCIP_Bool* isxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */
1555  SCIP_Bool* isxjldiag, /**< is X_jl diagonal? (i.e. j = l) */
1556  int* basicvarpos2tableaurow,/**< map from basic var to its tableau row */
1557  SCIP_HASHMAP* tableau, /**< map from var to its tableau row */
1558  SCIP_RESULT* result /**< pointer to store the result of the separation call */
1559  )
1560 {
1561  SCIP_ROWPREP* rowprep;
1562  SCIP_VAR* vars[4] = {xik, xjl, xil, xjk};
1563  SCIP_Real* tableaurows[4];
1564  SCIP_Real* interpoints;
1565  SCIP_Real* rays;
1566  int nrays;
1567  int* rayslppos;
1568  int ncols;
1569  int nrows;
1570  SCIP_Bool success;
1571  SCIP_Real ad[4] = {0.0, 0.0, 0.0, 0.0};
1572  SCIP_Real solxik;
1573  SCIP_Real solxil;
1574  SCIP_Real solxjk;
1575  SCIP_Real solxjl;
1576 
1577  ncols = SCIPgetNLPCols(scip);
1578  nrows = SCIPgetNLPRows(scip);
1579 
1580  /* allocate memory for intersection points */
1581  SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, ncols + nrows) );
1582 
1583  /* allocate memory for rays */
1584  SCIP_CALL( SCIPallocBufferArray(scip, &rays, 4 * (ncols + nrows)) );
1585  SCIP_CALL( SCIPallocBufferArray(scip, &rayslppos, ncols + nrows) );
1586 
1587  /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */
1588  SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, TRUE) );
1589  SCIProwprepAddSide(rowprep, 1.0);
1590 
1591  /* check if we have the tableau row of the variable and if not compute it */
1592  SCIP_CALL( getTableauRows(scip, vars, basicvarpos2tableaurow, tableau, tableaurows, &success) );
1593 
1594  if( ! success )
1595  goto CLEANUP;
1596 
1597  /* if we want to enforce bounds, set the right a and d to enforce aTx + dTy <= 0 */
1598  if( sepadata->usebounds )
1599  {
1600  solxik = SCIPvarGetLPSol(xik);
1601  solxil = SCIPvarGetLPSol(xil);
1602  solxjk = SCIPvarGetLPSol(xjk);
1603  solxjl = SCIPvarGetLPSol(xjl);
1604 
1605  if( isxikdiag && SCIPisFeasNegative(scip, solxik) )
1606  {
1607  ad[0] = -1.0;
1608  ad[2] = 1.0;
1609  }
1610  else if( isxjldiag && SCIPisFeasNegative(scip, solxjl) )
1611  {
1612  ad[0] = -1.0;
1613  ad[2] = -1.0;
1614  }
1615  else if( isxildiag && SCIPisFeasNegative(scip, solxil) )
1616  {
1617  ad[1] = 1.0;
1618  ad[3] = -1.0;
1619  }
1620  else if( isxjkdiag && SCIPisFeasNegative(scip, solxjk) )
1621  {
1622  ad[1] = -1.0;
1623  ad[3] = -1.0;
1624  }
1625  }
1626 
1627  nrays = 0;
1628  /* loop over each non-basic var; get the ray; compute cut coefficient */
1629  SCIP_CALL( addCols(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1630 
1631  if( ! success )
1632  goto CLEANUP;
1633 
1634  /* loop over non-basic slack variables */
1635  SCIP_CALL( addRows(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1636 
1637  if( ! success )
1638  goto CLEANUP;
1639 
1640  /* do strengthening */
1641  if( sepadata->usestrengthening )
1642  {
1643  SCIP_CALL( computeNegCutcoefs(scip, vars, rays, nrays, rayslppos, interpoints, rowprep, sepadata->usebounds, ad, &success) );
1644 
1645  if( ! success )
1646  goto CLEANUP;
1647  }
1648 
1649  /* merge coefficients that belong to same variable */
1650  SCIPmergeRowprepTerms(scip, rowprep);
1651 
1652  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, sepadata->mincutviol, NULL, &success) );
1653 
1654  /* if cleanup was successfull, create row out of rowprep and add it */
1655  if( success )
1656  {
1657  SCIP_ROW* row;
1658  SCIP_Bool infeasible;
1659 
1660  /* create row */
1661  SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) );
1662 
1663  assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0);
1664 
1665  /* add row */
1666  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
1667 
1668  if( infeasible )
1669  *result = SCIP_CUTOFF;
1670  else
1671  *result = SCIP_SEPARATED;
1672 
1673  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1674  }
1675 
1676 CLEANUP:
1677  SCIPfreeRowprep(scip, &rowprep);
1678  SCIPfreeBuffer(scip, &rayslppos);
1679  SCIPfreeBuffer(scip, &rays);
1680  SCIPfreeBuffer(scip, &interpoints);
1681 
1682  return SCIP_OKAY;
1683 }
1684 
1685 
1686 /** separates cuts for stored principal minors */
1687 static
1689  SCIP* scip, /**< SCIP data structure */
1690  SCIP_SEPA* sepa, /**< separator */
1691  SCIP_RESULT* result /**< pointer to store the result of the separation call */
1692  )
1693 {
1694  SCIP_SEPADATA* sepadata;
1695  SCIP_HASHMAP* tableau = NULL;
1696  int* basicvarpos2tableaurow = NULL; /* map between basic var and its tableau row */
1697  int i;
1698 
1699  assert(sepa != NULL);
1700  assert(result != NULL);
1701 
1702  *result = SCIP_DIDNOTRUN;
1703 
1704  sepadata = SCIPsepaGetData(sepa);
1705  assert(sepadata != NULL);
1706 
1707  /* check whether there are some minors available */
1708  if( sepadata->nminors == 0 )
1709  return SCIP_OKAY;
1710 
1711  *result = SCIP_DIDNOTFIND;
1712 
1713  /* loop over the minors and if they are violated build cut */
1714  for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i )
1715  {
1716  SCIP_VAR* auxvarxik;
1717  SCIP_VAR* auxvarxil;
1718  SCIP_VAR* auxvarxjk;
1719  SCIP_VAR* auxvarxjl;
1720  SCIP_Bool isauxvarxikdiag;
1721  SCIP_Bool isauxvarxildiag;
1722  SCIP_Bool isauxvarxjkdiag;
1723  SCIP_Bool isauxvarxjldiag;
1724  SCIP_Real solxik;
1725  SCIP_Real solxil;
1726  SCIP_Real solxjk;
1727  SCIP_Real solxjl;
1728  SCIP_Real det;
1729 
1730  /* get variables of the i-th minor */
1731  SCIP_CALL( getMinorVars(sepadata, i, &auxvarxik, &auxvarxil, &auxvarxjk, &auxvarxjl, &isauxvarxikdiag,
1732  &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag) );
1733 
1734  /* get current solution values */
1735  solxik = SCIPvarGetLPSol(auxvarxik);
1736  solxil = SCIPvarGetLPSol(auxvarxil);
1737  solxjk = SCIPvarGetLPSol(auxvarxjk);
1738  solxjl = SCIPvarGetLPSol(auxvarxjl);
1739 
1740  det = solxik * solxjl - solxil * solxjk;
1741 
1742  if( SCIPisFeasZero(scip, det) )
1743  continue;
1744 
1745  if( basicvarpos2tableaurow == NULL )
1746  {
1747  /* allocate memory */
1748  SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, SCIPgetNLPCols(scip)) );
1749  SCIP_CALL( SCIPhashmapCreate(&tableau, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1750 
1751  /* construct basicvar to tableau row map */
1752  SCIP_CALL( constructBasicVars2TableauRowMap(scip, basicvarpos2tableaurow) );
1753  }
1754  assert(tableau != NULL);
1755 
1756  if( SCIPisFeasPositive(scip, det) )
1757  {
1758  SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxik, auxvarxil, auxvarxjk, auxvarxjl, &isauxvarxikdiag,
1759  &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag, basicvarpos2tableaurow, tableau, result) );
1760  }
1761  else
1762  {
1763  assert(SCIPisFeasNegative(scip, det));
1764  SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxil, auxvarxik, auxvarxjl, auxvarxjk, &isauxvarxildiag,
1765  &isauxvarxikdiag, &isauxvarxjldiag, &isauxvarxjkdiag, basicvarpos2tableaurow, tableau, result) );
1766  }
1767  }
1768 
1769  /* all minors were feasible, so no memory to free */
1770  if( basicvarpos2tableaurow == NULL )
1771  return SCIP_OKAY;
1772 
1773  /* free memory */
1774  for( i = 0; i < SCIPhashmapGetNEntries(tableau); ++i )
1775  {
1776  SCIP_HASHMAPENTRY* entry;
1777 
1778  entry = SCIPhashmapGetEntry(tableau, i);
1779 
1780  if( entry != NULL )
1781  {
1782  SCIP_Real* tableaurow;
1783 
1784  tableaurow = (SCIP_Real *) SCIPhashmapEntryGetImage(entry);
1785 
1786  SCIPfreeBufferArrayNull(scip, &tableaurow);
1787  }
1788  }
1789  SCIPhashmapFree(&tableau);
1790  SCIPfreeBufferArray(scip, &basicvarpos2tableaurow);
1791 
1792  return SCIP_OKAY;
1793 }
1794 
1795 /*
1796  * Callback methods of separator
1797  */
1798 
1799 /** copy method for separator plugins (called when SCIP copies plugins) */
1800 static
1801 SCIP_DECL_SEPACOPY(sepaCopyMinor)
1802 { /*lint --e{715}*/
1803  assert(scip != NULL);
1804  assert(sepa != NULL);
1805  assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
1806 
1807  /* call inclusion method of constraint handler */
1809 
1810  return SCIP_OKAY;
1811 }
1812 
1813 
1814 /** destructor of separator to free user data (called when SCIP is exiting) */
1815 static
1816 SCIP_DECL_SEPAFREE(sepaFreeMinor)
1817 { /*lint --e{715}*/
1818  SCIP_SEPADATA* sepadata;
1819 
1820  sepadata = SCIPsepaGetData(sepa);
1821  assert(sepadata != NULL);
1822  assert(sepadata->minors == NULL);
1823  assert(sepadata->nminors == 0);
1824  assert(sepadata->minorssize == 0);
1825 
1826  /* free separator data */
1827  SCIPfreeBlockMemory(scip, &sepadata);
1828  SCIPsepaSetData(sepa, NULL);
1829 
1830  return SCIP_OKAY;
1831 }
1832 
1833 
1834 /** initialization method of separator (called after problem was transformed) */
1835 static
1836 SCIP_DECL_SEPAINIT(sepaInitMinor)
1837 { /*lint --e{715}*/
1838  SCIP_SEPADATA* sepadata;
1839 
1840  /* get separator data */
1841  sepadata = SCIPsepaGetData(sepa);
1842  assert(sepadata != NULL);
1843  assert(sepadata->randnumgen == NULL);
1844 
1845  /* create random number generator */
1846  SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
1847 
1848  return SCIP_OKAY;
1849 }
1850 
1851 
1852 /** deinitialization method of separator (called before transformed problem is freed) */
1853 static
1854 SCIP_DECL_SEPAEXIT(sepaExitMinor)
1855 { /*lint --e{715}*/
1856  SCIP_SEPADATA* sepadata;
1857 
1858  /* get separator data */
1859  sepadata = SCIPsepaGetData(sepa);
1860  assert(sepadata != NULL);
1861  assert(sepadata->randnumgen != NULL);
1862 
1863  /* free random number generator */
1864  SCIPfreeRandom(scip, &sepadata->randnumgen);
1865 
1866  return SCIP_OKAY;
1867 }
1868 
1869 
1870 /** solving process initialization method of separator (called when branch and bound process is about to begin) */
1871 static
1872 SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
1873 { /*lint --e{715}*/
1874  return SCIP_OKAY;
1875 }
1876 
1877 
1878 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */
1879 static
1880 SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
1881 { /*lint --e{715}*/
1882  SCIP_SEPADATA* sepadata;
1883 
1884  sepadata = SCIPsepaGetData(sepa);
1885  assert(sepadata != NULL);
1886 
1887  /* clear separation data */
1888  SCIP_CALL( sepadataClear(scip, sepadata) );
1889 
1890  return SCIP_OKAY;
1891 }
1892 
1893 
1894 /** LP solution separation method of separator */
1895 static
1896 SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
1897 { /*lint --e{715}*/
1898  SCIP_SEPADATA* sepadata;
1899  int ncalls;
1900  int currentdepth;
1901 
1902  /* need routine to compute eigenvalues/eigenvectors */
1903  if( !SCIPisIpoptAvailableIpopt() )
1904  return SCIP_OKAY;
1905 
1906  sepadata = SCIPsepaGetData(sepa);
1907  assert(sepadata != NULL);
1908  currentdepth = SCIPgetDepth(scip);
1909  ncalls = SCIPsepaGetNCallsAtNode(sepa);
1910 
1911  /* only call the separator a given number of times at each node */
1912  if( (currentdepth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
1913  || (currentdepth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
1914  {
1915  SCIPdebugMsg(scip, "reached round limit for node\n");
1916  return SCIP_OKAY;
1917  }
1918 
1919  /* try to detect minors */
1920  SCIP_CALL( detectMinors(scip, sepadata) );
1921 
1922  /* call separation method */
1923  SCIP_CALL( separatePoint(scip, sepa, result) );
1924 
1925  return SCIP_OKAY;
1926 }
1927 
1928 
1929 /*
1930  * separator specific interface methods
1931  */
1932 
1933 /** creates the minor separator and includes it in SCIP */
1935  SCIP* scip /**< SCIP data structure */
1936  )
1937 {
1938  SCIP_SEPADATA* sepadata = NULL;
1939  SCIP_SEPA* sepa = NULL;
1940 
1941  /* create minor separator data */
1942  SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
1943  BMSclearMemory(sepadata);
1944 
1945  /* include separator */
1948  sepaExeclpMinor, NULL,
1949  sepadata) );
1950 
1951  assert(sepa != NULL);
1952 
1953  /* set non fundamental callbacks via setter functions */
1954  SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) );
1955  SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) );
1956  SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) );
1957  SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) );
1958  SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) );
1959  SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) );
1960 
1961  /* add minor separator parameters */
1963  "separating/" SEPA_NAME "/usestrengthening",
1964  "whether to use strengthened intersection cuts to separate minors",
1965  &sepadata->usestrengthening, FALSE, DEFAULT_USESTRENGTHENING, NULL, NULL) );
1966 
1968  "separating/" SEPA_NAME "/usebounds",
1969  "whether to also enforce nonegativity bounds of principle minors",
1970  &sepadata->usebounds, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) );
1971 
1973  "separating/" SEPA_NAME "/mincutviol",
1974  "minimum required violation of a cut",
1975  &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1976 
1977  SCIP_CALL( SCIPaddIntParam(scip,
1978  "separating/" SEPA_NAME "/maxrounds",
1979  "maximal number of separation rounds per node (-1: unlimited)",
1980  &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
1981 
1982  SCIP_CALL( SCIPaddIntParam(scip,
1983  "separating/" SEPA_NAME "/maxroundsroot",
1984  "maximal number of separation rounds in the root node (-1: unlimited)",
1985  &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
1986 
1987  return SCIP_OKAY;
1988 }
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:61
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
SCIP_ROW ** SCIPgetLPRows(SCIP *scip)
Definition: scip_lp.c:605
SCIP_Bool SCIPisFeasZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
#define SCIPreallocBlockMemoryArray(scip, ptr, oldnum, newnum)
Definition: scip_mem.h:99
#define SEPA_USESSUBSCIP
void * SCIPhashmapEntryGetImage(SCIP_HASHMAPENTRY *entry)
Definition: misc.c:3520
SCIP_RETCODE SCIPgetLPBInvRow(SCIP *scip, int r, SCIP_Real *coefs, int *inds, int *ninds)
Definition: scip_lp.c:714
static SCIP_RETCODE separateDeterminant(SCIP *scip, SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata, SCIP_VAR *xik, SCIP_VAR *xil, SCIP_VAR *xjk, SCIP_VAR *xjl, SCIP_Bool *isxikdiag, SCIP_Bool *isxildiag, SCIP_Bool *isxjkdiag, SCIP_Bool *isxjldiag, int *basicvarpos2tableaurow, SCIP_HASHMAP *tableau, SCIP_RESULT *result)
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:500
static SCIP_RETCODE detectMinors(SCIP *scip, SCIP_SEPADATA *sepadata)
static SCIP_Real isCase4a(SCIP_Real tsol, SCIP_Real *coefs, SCIP_Real *coefscondition)
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define DEFAULT_USEBOUNDS
#define SEPA_MAXBOUNDDIST
#define SCIPquadprecSumQD(r, a, b)
Definition: dbldblarith.h:62
SCIP_CONSHDLR * SCIPfindConshdlr(SCIP *scip, const char *name)
Definition: scip_cons.c:886
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3808
#define DEFAULT_MINCUTVIOL
static SCIP_RETCODE getTableauRows(SCIP *scip, SCIP_VAR **vars, int *basicvarpos2tableaurow, SCIP_HASHMAP *tableau, SCIP_Real **tableaurows, SCIP_Bool *success)
SCIP_BASESTAT SCIPcolGetBasisStatus(SCIP_COL *col)
Definition: lp.c:17034
int SCIPcalcMemGrowSize(SCIP *scip, int num)
Definition: scip_mem.c:139
SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:17979
static SCIP_DECL_SEPACOPY(sepaCopyMinor)
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1248
#define DEFAULT_RANDSEED
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
int SCIProwGetNLPNonz(SCIP_ROW *row)
Definition: lp.c:17230
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4554
private functions to work with algebraic expressions
SCIP_Real SCIProwGetLhs(SCIP_ROW *row)
Definition: lp.c:17295
#define FALSE
Definition: def.h:96
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition: misc.c:3024
SCIP_RETCODE SCIPgetRowprepRowSepa(SCIP *scip, SCIP_ROW **row, SCIP_ROWPREP *rowprep, SCIP_SEPA *sepa)
SCIP_BASESTAT SCIProwGetBasisStatus(SCIP_ROW *row)
Definition: lp.c:17343
static SCIP_Real computeRoot(SCIP *scip, SCIP_Real *coefs)
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3352
#define TRUE
Definition: def.h:95
const char * SCIPsepaGetName(SCIP_SEPA *sepa)
Definition: sepa.c:743
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
#define SCIPstatisticMessage
Definition: pub_message.h:123
#define BINSEARCH_MAXITERS
void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_Real lincoeff, SCIP_Real rhs, SCIP_INTERVAL xbnds)
int SCIPvarGetProbindex(SCIP_VAR *var)
Definition: var.c:17613
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:108
#define QUAD_SCALE(x, a)
Definition: dbldblarith.h:50
void * SCIPhashmapGetImage(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3211
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:136
#define QUAD_ASSIGN(a, constant)
Definition: dbldblarith.h:51
static SCIP_RETCODE computeRestrictionToRay(SCIP *scip, SCIP_Real *ray, SCIP_VAR **vars, SCIP_Real *coefs, SCIP_Real *coefs4b, SCIP_Real *coefscondition, SCIP_Bool usebounds, SCIP_Real *ad, SCIP_Bool *success)
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:89
variable expression handler
SCIP_RETCODE SCIPsetSepaCopy(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPACOPY((*sepacopy)))
Definition: scip_sepa.c:151
#define SEPA_DELAY
#define SCIPdebugMsg
Definition: scip_message.h:78
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:83
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
#define SEPA_NAME
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1454
#define QUAD_TO_DBL(x)
Definition: dbldblarith.h:49
SCIP_SEPADATA * SCIPsepaGetData(SCIP_SEPA *sepa)
Definition: sepa.c:633
static SCIP_RETCODE insertIndex(SCIP *scip, SCIP_HASHMAP *rowmap, SCIP_VAR *row, SCIP_VAR *col, SCIP_VAR *auxvar, int *rowindices, int *nrows)
static SCIP_Real computeIntersectionPoint(SCIP *scip, SCIP_Bool usebounds, SCIP_Real *coefs, SCIP_Real *coefs4b, SCIP_Real *coefscondition)
const char * SCIPgetProbName(SCIP *scip)
Definition: scip_prob.c:1075
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3373
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3818
#define DEFAULT_USESTRENGTHENING
static SCIP_RETCODE addRows(SCIP *scip, SCIP_VAR **vars, SCIP_Real **tableaurows, SCIP_ROWPREP *rowprep, SCIP_Real *rays, int *nrays, int *rayslppos, SCIP_Real *interpoints, SCIP_Bool usebounds, SCIP_Real *ad, SCIP_Bool *success)
SCIP_Real SCIPcolGetPrimsol(SCIP_COL *col)
Definition: lp.c:16999
static SCIP_RETCODE computeNegCutcoefs(SCIP *scip, SCIP_VAR **vars, SCIP_Real *rays, int nrays, int *rayslppos, SCIP_Real *interpoints, SCIP_ROWPREP *rowprep, SCIP_Bool usebounds, SCIP_Real *ad, SCIP_Bool *success)
int SCIPhashmapGetNEntries(SCIP_HASHMAP *hashmap)
Definition: misc.c:3491
SCIP_HASHMAPENTRY * SCIPhashmapGetEntry(SCIP_HASHMAP *hashmap, int entryidx)
Definition: misc.c:3499
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
SCIP_RETCODE SCIPsetSepaExit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAEXIT((*sepaexit)))
Definition: scip_sepa.c:199
#define SCIPallocBuffer(scip, ptr)
Definition: scip_mem.h:122
static void doBinarySearch(SCIP *scip, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e, SCIP_Real *sol)
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:137
int SCIPsepaGetNCallsAtNode(SCIP_SEPA *sepa)
Definition: sepa.c:880
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:57
#define QUAD(x)
Definition: dbldblarith.h:47
const char * SCIPconsGetName(SCIP_CONS *cons)
Definition: cons.c:8090
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17264
static SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
void SCIPcomputeArraysIntersectionInt(int *array1, int narray1, int *array2, int narray2, int *intersectarray, int *nintersectarray)
Definition: misc.c:10470
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition: misc.c:3058
void SCIPsepaSetData(SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata)
Definition: sepa.c:643
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
#define NULL
Definition: lpi_spx1.cpp:164
int * vals
power and signed power expression handlers
SCIP_Real SCIPvarGetLPSol(SCIP_VAR *var)
Definition: var.c:18297
SCIP_RETCODE SCIPsetSepaInitsol(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAINITSOL((*sepainitsol)))
Definition: scip_sepa.c:215
int SCIPgetNLPRows(SCIP *scip)
Definition: scip_lp.c:626
void SCIPfreeRowprep(SCIP *scip, SCIP_ROWPREP **rowprep)
Definition: misc_rowprep.c:581
SCIP_HASHMAP * auxvars
#define SCIP_CALL(x)
Definition: def.h:394
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
Definition: misc_rowprep.c:737
static SCIP_RETCODE findRho(SCIP *scip, SCIP_Real *rays, int nrays, int idx, SCIP_Real *interpoints, SCIP_VAR **vars, SCIP_Real *rho, SCIP_Bool usebounds, SCIP_Real *ad, SCIP_Bool *success)
SCIP_Real SCIProwGetRhs(SCIP_ROW *row)
Definition: lp.c:17305
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:250
#define SCIPquadprecProdDD(r, a, b)
Definition: dbldblarith.h:58
#define SCIP_INTERVAL_INFINITY
Definition: def.h:208
SCIP_COL ** SCIProwGetCols(SCIP_ROW *row)
Definition: lp.c:17241
int SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4597
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2302
#define SCIPquadprecSqrtQ(r, a)
Definition: dbldblarith.h:71
SCIP_RETCODE SCIPincludeSepaBasic(SCIP *scip, SCIP_SEPA **sepa, const char *name, const char *desc, int priority, int freq, SCIP_Real maxbounddist, SCIP_Bool usessubscip, SCIP_Bool delay, SCIP_DECL_SEPAEXECLP((*sepaexeclp)), SCIP_DECL_SEPAEXECSOL((*sepaexecsol)), SCIP_SEPADATA *sepadata)
Definition: scip_sepa.c:109
SCIP_Bool SCIPisHugeValue(SCIP *scip, SCIP_Real val)
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
Ipopt NLP interface.
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:124
SCIP_RETCODE SCIPgetLPBInvARow(SCIP *scip, int r, SCIP_Real *binvrow, SCIP_Real *coefs, int *inds, int *ninds)
Definition: scip_lp.c:785
SCIP_Real * SCIProwGetVals(SCIP_ROW *row)
Definition: lp.c:17251
SCIP_RETCODE SCIPsetSepaExitsol(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAEXITSOL((*sepaexitsol)))
Definition: scip_sepa.c:231
static SCIP_DECL_SEPAINIT(sepaInitMinor)
static SCIP_RETCODE sepadataAddMinor(SCIP *scip, SCIP_SEPADATA *sepadata, SCIP_VAR *auxvarxik, SCIP_VAR *auxvarxil, SCIP_VAR *auxvarxjk, SCIP_VAR *auxvarxjl, SCIP_Bool isauxvarxikdiag, SCIP_Bool isauxvarxildiag, SCIP_Bool isauxvarxjkdiag, SCIP_Bool isauxvarxjldiag)
#define SCIP_Bool
Definition: def.h:93
SCIP_EXPR * SCIPexpriterRestartDFS(SCIP_EXPRITER *iterator, SCIP_EXPR *expr)
Definition: expriter.c:629
void SCIProwprepAddConstant(SCIP_ROWPREP *rowprep, SCIP_Real constant)
Definition: misc_rowprep.c:751
static SCIP_RETCODE addRowToCut(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_Real cutcoef, SCIP_ROW *row, SCIP_Bool *success)
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:670
constraint handler for nonlinear constraints specified by algebraic expressions
void SCIPmergeRowprepTerms(SCIP *scip, SCIP_ROWPREP *rowprep)
#define SCIPquadprecSquareD(r, a)
Definition: dbldblarith.h:59
#define MAX(x, y)
Definition: tclique_def.h:92
static SCIP_RETCODE sepadataClear(SCIP *scip, SCIP_SEPADATA *sepadata)
SCIP_Real SCIPgetCutEfficacy(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:94
SCIP_EXPR * SCIPexpriterGetNext(SCIP_EXPRITER *iterator)
Definition: expriter.c:857
static SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
#define SCIPquadprecProdQD(r, a, b)
Definition: dbldblarith.h:63
static SCIP_RETCODE separatePoint(SCIP *scip, SCIP_SEPA *sepa, SCIP_RESULT *result)
#define SEPA_FREQ
SCIP_COL * SCIPvarGetCol(SCIP_VAR *var)
Definition: var.c:17634
SCIP_RETCODE SCIPgetLPBasisInd(SCIP *scip, int *basisind)
Definition: scip_lp.c:686
SCIP_COL ** SCIPgetLPCols(SCIP *scip)
Definition: scip_lp.c:506
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1465
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
#define BMSclearMemory(ptr)
Definition: memory.h:131
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:67
static SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
int SCIPgetNVars(SCIP *scip)
Definition: scip_prob.c:2000
#define SCIP_REAL_MAX
Definition: def.h:187
void SCIPfreeExpriter(SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2316
product expression handler
SCIP_Real SCIProwGetConstant(SCIP_ROW *row)
Definition: lp.c:17261
SCIP_VAR ** b
Definition: circlepacking.c:65
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1562
static SCIP_DECL_SEPAFREE(sepaFreeMinor)
#define DEFAULT_MAXROUNDSROOT
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:134
SCIP_RETCODE SCIPsetSepaFree(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAFREE((*sepafree)))
Definition: scip_sepa.c:167
SCIP_RETCODE SCIPsetSepaInit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAINIT((*sepainit)))
Definition: scip_sepa.c:183
#define SEPA_PRIORITY
SCIP_VAR * SCIPcolGetVar(SCIP_COL *col)
Definition: lp.c:17045
static SCIP_Real evalPhiAtRay(SCIP *scip, SCIP_Real t, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e)
SCIP_VAR * a
Definition: circlepacking.c:66
SCIP_Bool SCIPisFeasPositive(SCIP *scip, SCIP_Real val)
SCIP_VAR ** SCIPgetVars(SCIP *scip)
Definition: scip_prob.c:1955
int SCIProwGetLPPos(SCIP_ROW *row)
Definition: lp.c:17504
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
SCIP_RETCODE SCIPcaptureVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_var.c:1214
#define SCIP_Real
Definition: def.h:186
static SCIP_RETCODE getMinorVars(SCIP_SEPADATA *sepadata, int idx, SCIP_VAR **auxvarxik, SCIP_VAR **auxvarxil, SCIP_VAR **auxvarxjk, SCIP_VAR **auxvarxjl, SCIP_Bool *isauxvarxikdiag, SCIP_Bool *isauxvarxildiag, SCIP_Bool *isauxvarxjkdiag, SCIP_Bool *isauxvarxjldiag)
static SCIP_RETCODE addCols(SCIP *scip, SCIP_VAR **vars, SCIP_Real **tableaurows, SCIP_ROWPREP *rowprep, SCIP_Real *rays, int *nrays, int *rayslppos, SCIP_Real *interpoints, SCIP_Bool usebounds, SCIP_Real *ad, SCIP_Bool *success)
#define SEPA_DESC
static SCIP_Bool raysAreDependent(SCIP *scip, SCIP_Real *ray1, SCIP_Real *ray2, SCIP_Real *coef)
void SCIPsortInt(int *intarray, int len)
#define DEFAULT_MAXROUNDS
SCIP_Real SCIPgetTotalTime(SCIP *scip)
Definition: scip_timing.c:351
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
Definition: misc_rowprep.c:561
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_RETCODE SCIPincludeSepaInterminor(SCIP *scip)
int SCIPgetNLPCols(SCIP *scip)
Definition: scip_lp.c:527
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17989
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:904
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:111
SCIP_RETCODE SCIPcleanupRowprep(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Real minviol, SCIP_Real *viol, SCIP_Bool *success)
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
SCIP_RETCODE SCIPhashmapInsert(SCIP_HASHMAP *hashmap, void *origin, void *image)
Definition: misc.c:3106
SCIP_Bool SCIPexpriterIsEnd(SCIP_EXPRITER *iterator)
Definition: expriter.c:968
static SCIP_RETCODE constructBasicVars2TableauRowMap(SCIP *scip, int *map)
static SCIP_DECL_SEPAEXIT(sepaExitMinor)
int SCIPcolGetLPPos(SCIP_COL *col)
Definition: lp.c:17096
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:139
struct SCIP_SepaData SCIP_SEPADATA
Definition: type_sepa.h:52
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:57
static SCIP_RETCODE addColToCut(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_Real cutcoef, SCIP_COL *col)
SCIP_Real SCIPgetRowActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:2101
#define SCIPreallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:128