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-2024 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  SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors);
493 
494  /* free memory */
495  SCIPfreeBufferArray(scip, &intersection);
496  SCIPfreeBufferArray(scip, &rowvars);
497  SCIPhashmapFree(&rowmap);
498  SCIPfreeExpriter(&it);
499 
500 #ifdef SCIP_STATISTIC
501  totaltime += SCIPgetTotalTime(scip);
502  SCIPstatisticMessage("MINOR DETECT %s %f %d %d\n", SCIPgetProbName(scip), totaltime, sepadata->nminors, maxminors);
503 #endif
504 
505  return SCIP_OKAY;
506 }
507 
508 /** constructs map between lp position of a basic variable and its row in the tableau */
509 static
511  SCIP* scip, /**< SCIP data structure */
512  int* map /**< buffer to store the map */
513  )
514 {
515  int* basisind;
516  int nrows;
517  int i;
518 
519  nrows = SCIPgetNLPRows(scip);
520  SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
521 
522  SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
523  for( i = 0; i < nrows; ++i )
524  {
525  if( basisind[i] >= 0 )
526  map[basisind[i]] = i;
527  }
528 
529  SCIPfreeBufferArray(scip, &basisind);
530 
531  return SCIP_OKAY;
532 }
533 
534 /** The restriction of the function representing the maximal S-free set to zlp + t * ray has the form
535  * sqrt(A t^2 + B t + C) - (D t + E).
536  * This function computes the coefficients A, B, C, D, E for the given ray.
537  */
538 static
540  SCIP* scip, /**< SCIP data structure */
541  SCIP_Real* ray, /**< coefficients of ray */
542  SCIP_VAR** vars, /**< variables */
543  SCIP_Real* coefs, /**< buffer to store A, B, C, D, and E of cases 1, 2, 3, or 4a*/
544  SCIP_Real* coefs4b, /**< buffer to store A, B, C, D, and E of case 4b */
545  SCIP_Real* coefscondition, /**< buffer to store coefs for checking whether we are in case 4a or 4b */
546  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
547  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
548  SCIP_Bool* success /**< FALSE if we need to abort generation because of numerics */
549  )
550 {
551  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};
552  SCIP_Real eigenvalues[4] = {0.5, 0.5, -0.5, -0.5};
553  SCIP_Real eigencoef = 0.7071067811865475244008443621048490;
554  SCIP_Real* a;
555  SCIP_Real* b;
556  SCIP_Real* c;
557  SCIP_Real* d;
558  SCIP_Real* e;
559  SCIP_Real min;
560  SCIP_Real max;
561  SCIP_Real norm1;
562  SCIP_Real norm2;
563  int negidx;
564  int posidx;
565  int i;
566 
567  *success = TRUE;
568 
569  /* set all coefficients to zero */
570  memset(coefs, 0, 5 * sizeof(SCIP_Real));
571  memset(coefs4b, 0, 5 * sizeof(SCIP_Real));
572  norm1 = 0.0;
573  norm2 = 0.0;
574 
575  a = coefs;
576  b = coefs + 1;
577  c = coefs + 2;
578  d = coefs + 3;
579  e = coefs + 4;
580 
581  negidx = 2;
582  posidx = 0;
583  for( i = 0; i < 4; ++i )
584  {
585  int j;
586  SCIP_Real vzlp;
587  SCIP_Real vdotray;
588 
589  vzlp = 0;
590  vdotray = 0;
591 
592  /* compute eigenvec * ray and eigenvec * solution */
593  for( j = 0; j < 4; ++j )
594  {
595  vdotray += eigencoef * eigenvectors[4 * i + j] * ray[j];
596  vzlp += eigencoef * eigenvectors[4 * i + j] * SCIPvarGetLPSol(vars[j]);
597  }
598 
599  if( eigenvalues[i] > 0 )
600  {
601  /* positive eigenvalue: compute D and E */
602  *d += eigenvalues[i] * vzlp * vdotray;
603  *e += eigenvalues[i] * SQR( vzlp );
604 
605  if( usebounds )
606  {
607  norm1 += eigenvalues[i] * (1 - SQR( ad[posidx] )) * SQR( vzlp );
608  norm2 += sqrt( eigenvalues[i] ) * ad[posidx] * vzlp;
609  ++posidx;
610  }
611  }
612  else
613  {
614  /* negative eigenvalue: compute A, B, and C */
615  *a -= eigenvalues[i] * SQR( vdotray );
616  *b -= 2.0 * eigenvalues[i] * vzlp * vdotray;
617  *c -= eigenvalues[i] * SQR( vzlp );
618 
619  if( usebounds )
620  {
621  coefs4b[0] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vdotray );
622  coefs4b[1] -= 2.0 * eigenvalues[i] * (1 - SQR( ad[negidx] )) * vzlp * vdotray;
623  coefs4b[2] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vzlp );
624  coefs4b[3] += sqrt( -eigenvalues[i] ) * ad[negidx] * vdotray;
625  coefs4b[4] += sqrt( -eigenvalues[i] ) * ad[negidx] * vzlp;
626  ++negidx;
627  }
628  }
629  }
630 
631  assert(*e > 0);
632 
633  if( sqrt( *c ) - sqrt( *e ) >= 0.0 )
634  {
635  assert(sqrt( *c ) - sqrt( *e ) < 1e-6);
636  *success = FALSE;
637  return SCIP_OKAY;
638  }
639 
640  /* finish computation of coefficients when using bounds */
641  if( usebounds )
642  {
643  coefscondition[0] = norm2 / sqrt( *e );
644  coefscondition[1] = coefs4b[3];
645  coefscondition[2] = coefs4b[4];
646 
647  coefs4b[0] *= norm1 / *e;
648  coefs4b[1] *= norm1 / *e;
649  coefs4b[2] *= norm1 / *e;
650  coefs4b[3] *= norm2 / sqrt( *e );
651  coefs4b[4] *= norm2 / sqrt( *e );
652 
653  coefs4b[3] += *d / sqrt( *e );
654  coefs4b[4] += sqrt( *e );
655 
656  assert( sqrt( coefs4b[2] ) - coefs4b[4] < 0.0 );
657  }
658 
659  /* finish computation of D and E */
660  *e = sqrt( *e );
661  *d /= *e;
662 
663  /* maybe we want to avoid a large dynamism between A, B and C */
664  max = 0.0;
665  min = SCIPinfinity(scip);
666  for( i = 0; i < 3; ++i )
667  {
668  SCIP_Real absval;
669 
670  absval = ABS(coefs[i]);
671  if( max < absval )
672  max = absval;
673  if( absval != 0.0 && absval < min )
674  min = absval;
675  }
676 
677  if( SCIPisHugeValue(scip, max / min) )
678  {
679 #ifdef DEBUG_INTERSECTIONCUT
680  printf("Bad numerics: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min);
681 #endif
682  *success = FALSE;
683  return SCIP_OKAY;
684  }
685 
686  /* some sanity checks */
687  assert(*c >= 0); /* radicand at zero */
688  assert(sqrt( *c ) - *e < 0); /* the function at 0 must be negative */
689  assert(*a >= 0); /* the function inside the root is convex */
690 
691 #ifdef DEBUG_INTERSECTIONCUT
692  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]);
693 #endif
694 
695  return SCIP_OKAY;
696 }
697 
698 /** returns phi(zlp + t * ray) = sqrt(A t^2 + B t + C) - (D t + E) */ /*lint -e{715}*/
699 static
701  SCIP* scip, /**< SCIP data structure */
702  SCIP_Real t, /**< argument of phi restricted to ray */
703  SCIP_Real a, /**< value of A */
704  SCIP_Real b, /**< value of B */
705  SCIP_Real c, /**< value of C */
706  SCIP_Real d, /**< value of D */
707  SCIP_Real e /**< value of E */
708  )
709 {
710 #ifdef INTERCUTS_DBLDBL
711  SCIP_Real QUAD(lin);
712  SCIP_Real QUAD(disc);
713  SCIP_Real QUAD(tmp);
714  SCIP_Real QUAD(root);
715 
716  /* d * t + e */
717  SCIPquadprecProdDD(lin, d, t);
718  SCIPquadprecSumQD(lin, lin, e);
719 
720  /* a * t * t */
721  SCIPquadprecSquareD(disc, t);
722  SCIPquadprecProdQD(disc, disc, a);
723 
724  /* b * t */
725  SCIPquadprecProdDD(tmp, b, t);
726 
727  /* a * t * t + b * t */
728  SCIPquadprecSumQQ(disc, disc, tmp);
729 
730  /* a * t * t + b * t + c */
731  SCIPquadprecSumQD(disc, disc, c);
732 
733  /* sqrt(above): can't take sqrt of 0! */
734  if( QUAD_TO_DBL(disc) == 0 )
735  {
736  QUAD_ASSIGN(root, 0.0);
737  }
738  else
739  {
740  SCIPquadprecSqrtQ(root, disc);
741  }
742 
743  /* final result */
744  QUAD_SCALE(lin, -1.0);
745  SCIPquadprecSumQQ(tmp, root, lin);
746 
747  assert(!SCIPisInfinity(scip, t) || QUAD_TO_DBL(tmp) <= 0);
748 
749  return QUAD_TO_DBL(tmp);
750 #else
751  return sqrt( a * t * t + b * t + c ) - ( d * t + e );
752 #endif
753 }
754 
755 /** helper function of computeRoot: we want phi to be <= 0 */
756 static
758  SCIP* scip, /**< SCIP data structure */
759  SCIP_Real a, /**< value of A */
760  SCIP_Real b, /**< value of B */
761  SCIP_Real c, /**< value of C */
762  SCIP_Real d, /**< value of D */
763  SCIP_Real e, /**< value of E */
764  SCIP_Real* sol /**< buffer to store solution; also gives initial point */
765  )
766 {
767  SCIP_Real lb = 0.0;
768  SCIP_Real ub = *sol;
769  SCIP_Real curr;
770  int i;
771 
772  for( i = 0; i < BINSEARCH_MAXITERS; ++i )
773  {
774  SCIP_Real phival;
775 
776  curr = (lb + ub) / 2.0;
777  phival = evalPhiAtRay(scip, curr, a, b, c, d, e);
778 #ifdef INTERCUT_MOREDEBUG
779  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));
780 #endif
781 
782  if( phival <= 0.0 )
783  {
784  lb = curr;
785  if( SCIPisFeasZero(scip, phival) || SCIPisFeasEQ(scip, ub, lb) )
786  break;
787  }
788  else
789  ub = curr;
790  }
791 
792  *sol = lb;
793 }
794 
795 /** checks if we are in case 4a, i.e., if
796  * (num(xhat_{r+1}(zlp)) / E) * sqrt(A * tsol^2 + B * tsol + C) + w(ray) * tsol + num(yhat_{s+1}(zlp)) <= 0
797  */
798 static
800  SCIP_Real tsol, /**< t in the above formula */
801  SCIP_Real* coefs, /**< coefficients A, B, C, D, and E of case 4a */
802  SCIP_Real* coefscondition /**< extra coefficients needed for the evaluation of the condition:
803  * num(xhat_{r+1}(zlp)) / E; w(ray); num(yhat_{s+1}(zlp)) */
804  )
805 {
806  return (coefscondition[0] * sqrt( coefs[0] * SQR( tsol ) + coefs[1] * tsol + coefs[2] ) + coefscondition[1] *
807  tsol + coefscondition[2]) <= 0.0;
808 }
809 
810 /** finds smallest positive root phi by finding the smallest positive root of
811  * (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) = 0
812  *
813  * However, we are conservative and want a solution such that phi is negative, but close to 0;
814  * thus we correct the result with a binary search
815  */
816 static
818  SCIP* scip, /**< SCIP data structure */
819  SCIP_Real* coefs /**< value of A */
820  )
821 {
822  SCIP_Real sol;
823  SCIP_INTERVAL bounds;
824  SCIP_INTERVAL result;
825  SCIP_Real a = coefs[0];
826  SCIP_Real b = coefs[1];
827  SCIP_Real c = coefs[2];
828  SCIP_Real d = coefs[3];
829  SCIP_Real e = coefs[4];
830 
831  /* there is an intersection point if and only if sqrt(A) > D: here we are beliving in math, this might cause
832  * numerical issues
833  */
834  if( sqrt( a ) <= d )
835  {
836  sol = SCIPinfinity(scip);
837 
838  return sol;
839  }
840 
841  SCIPintervalSetBounds(&bounds, 0.0, SCIPinfinity(scip));
842 
843  /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds.
844  * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus
845  * tsol is the smallest t such that (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) >= 0
846  */
848  e, -(c - e * e), bounds);
849 
850  /* it can still be empty because of our infinity, I guess... */
852 
853  /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary
854  * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the
855  * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g.,
856  * ex8_3_1, bchoco05, etc
857  */
858  if( evalPhiAtRay(scip, sol, a, b, c, d, e) <= 1e-10 )
859  {
860 #ifdef INTERCUT_MOREDEBUG
861  printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e));
862  printf("don't do bin search\n");
863 #endif
864 
865  return sol;
866  }
867  else
868  {
869  /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */
870 #ifdef INTERCUT_MOREDEBUG
871  printf("do bin search because phival is %g\n", evalPhiAtRay(scip, sol, a, b, c, d, e));
872 #endif
873  doBinarySearch(scip, a, b, c, d, e, &sol);
874  }
875 
876  return sol;
877 }
878 
879 /** The maximal S-free set is gamma(z) <= 0; we find the intersection point of the ray `ray` starting from zlp with the
880  * boundary of the S-free set.
881  * That is, we find t >= 0 such that gamma(zlp + t * ray) = 0.
882  *
883  * In cases 1,2, and 3, gamma is of the form
884  * gamma(zlp + t * ray) = sqrt(A t^2 + B t + C) - (D t + E)
885  *
886  * In the case 4 gamma is of the form
887  * gamma(zlp + t * ray) = sqrt(A t^2 + B t + C) - (D t + E) if some condition holds
888  * sqrt(A' t^2 + B' t + C') - (D' t + E') otherwise
889  *
890  * It can be shown (given the special properties of gamma) that the smallest positive root of each function of the form
891  * sqrt(a t^2 + b t + c) - (d t + e)
892  * is the same as the smallest positive root of the quadratic equation:
893  * (sqrt(a t^2 + b t + c) - (d t + e)) * (sqrt(a t^2 + b t + c) + (d t + e)) = 0
894  * <==> (a - d^2) t^2 + (b - 2 d*e) t + (c - e^2) = 0
895  *
896  * So, in cases 1, 2, and 3, this function just returns the solution of the above equation.
897  * In case 4, it first solves the equation assuming we are in the first piece.
898  * If there is no solution, then the second piece can't have a solution (first piece >= second piece for all t)
899  * Then we check if the solution satisfies the condition.
900  * If it doesn't then we solve the equation for the second piece.
901  * If it has a solution, then it _has_ to be the solution.
902  */
903 static
905  SCIP* scip, /**< SCIP data structure */
906  SCIP_Bool usebounds, /**< whether we are in case 4 or not */
907  SCIP_Real* coefs, /**< values of A, B, C, D, and E of cases 1, 2, 3, or 4a */
908  SCIP_Real* coefs4b, /**< values of A, B, C, D, and E of case 4b */
909  SCIP_Real* coefscondition /**< values needed to evaluate condition of case 4 */
910  )
911 {
912  SCIP_Real sol;
913  SCIP_Real sol4b;
914 
915  assert(coefs != NULL);
916 
917  if( ! usebounds )
918  return computeRoot(scip, coefs);
919 
920  assert(coefs4b != NULL);
921  assert(coefscondition != NULL);
922 
923  /* compute solution of first piece */
924  sol = computeRoot(scip, coefs);
925 
926  /* if there is no solution --> second piece doesn't have solution */
927  if( SCIPisInfinity(scip, sol) )
928  {
929  /* this assert fails on multiplants_mtg5 the problem is that sqrt(A) <= D in 4a but not in 4b,
930  * now, this is impossible since the phi4a >= phi4b, so actually sqrt(A) is 10e-15 away from
931  * D in 4b
932  */
933  /* assert(SCIPisInfinity(scip, computeRoot(scip, coefs4b))); */
934  return sol;
935  }
936 
937  /* if solution of 4a is in 4a, then return */
938  if( isCase4a(sol, coefs, coefscondition) )
939  return sol;
940 
941  /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should
942  * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function
943  * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the
944  * binary search) we can find a slightly smaller solution; thus, we just keep the larger one
945  */
946  sol4b = computeRoot(scip, coefs4b);
947 
948  return MAX(sol, sol4b);
949 }
950 
951 /** adds cutcoef * (col - col*) to rowprep */
952 static
954  SCIP* scip, /**< SCIP data structure */
955  SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */
956  SCIP_Real cutcoef, /**< cut coefficient */
957  SCIP_COL* col /**< column to add to rowprep */
958  )
959 {
960  assert(col != NULL);
961 
962 #ifdef DEBUG_INTERCUTS_NUMERICS
963  SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)),
965  SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" :
966  "upper bound" , SCIPcolGetPrimsol(col));
967 #endif
968 
969  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(col), cutcoef) );
970  SCIProwprepAddConstant(rowprep, -cutcoef * SCIPcolGetPrimsol(col) );
971 
972  return SCIP_OKAY;
973 }
974 
975 /** adds cutcoef * (slack - slack*) to rowprep
976  *
977  * row is lhs <= <coefs, vars> + constant <= rhs, thus slack is defined by
978  * slack + <coefs, vars> + constant = side
979  * If row (slack) is at upper, it means that <coefs,vars*> + constant = rhs, and so
980  * slack* = side - rhs --> slack - slack* = rhs - <coefs, vars> - constant.
981  * If row (slack) is at lower, then <coefs,vars*> + constant = lhs, and so
982  * slack* = side - lhs --> slack - slack* = lhs - <coefs, vars> - constant.
983  */
984 static
986  SCIP* scip, /**< SCIP data structure */
987  SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */
988  SCIP_Real cutcoef, /**< cut coefficient */
989  SCIP_ROW* row, /**< row, whose slack we are adding to rowprep */
990  SCIP_Bool* success /**< buffer to store whether the row is nonbasic enough */
991  )
992 {
993  int i;
994  SCIP_COL** rowcols;
995  SCIP_Real* rowcoefs;
996  int nnonz;
997 
998  assert(row != NULL);
999 
1000  rowcols = SCIProwGetCols(row);
1001  rowcoefs = SCIProwGetVals(row);
1002  nnonz = SCIProwGetNLPNonz(row);
1003 
1004 #ifdef DEBUG_INTERCUTS_NUMERICS
1005  SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row));
1006  SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" :
1008  SCIPgetRowActivity(scip, row));
1009 #endif
1010 
1012  {
1013  assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row)));
1014  if( ! SCIPisFeasEQ(scip, SCIProwGetLhs(row), SCIPgetRowActivity(scip, row)) )
1015  {
1016  *success = FALSE;
1017  return SCIP_OKAY;
1018  }
1019 
1020  SCIProwprepAddConstant(rowprep, SCIProwGetLhs(row) * cutcoef);
1021  }
1022  else
1023  {
1024  assert(!SCIPisInfinity(scip, SCIProwGetRhs(row)));
1025  if( ! SCIPisFeasEQ(scip, SCIProwGetRhs(row), SCIPgetRowActivity(scip, row)) )
1026  {
1027  *success = FALSE;
1028  return SCIP_OKAY;
1029  }
1030 
1031  SCIProwprepAddConstant(rowprep, SCIProwGetRhs(row) * cutcoef);
1032  }
1033 
1034  for( i = 0; i < nnonz; i++ )
1035  {
1036  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) );
1037  }
1038 
1039  SCIProwprepAddConstant(rowprep, -SCIProwGetConstant(row) * cutcoef);
1040 
1041  return SCIP_OKAY;
1042 }
1043 
1044 /** get the tableau rows of the variables in vars */
1045 static
1047  SCIP* scip, /**< SCIP data structure */
1048  SCIP_VAR** vars, /**< variables in the minor */
1049  int* basicvarpos2tableaurow,/**< map between basic var and its tableau row */
1050  SCIP_HASHMAP* tableau, /**< map between var an its tableau row */
1051  SCIP_Real** tableaurows, /**< buffer to store tableau row */
1052  SCIP_Bool* success /**< set to TRUE if no variable had basisstat = ZERO */
1053  )
1054 {
1055  int v;
1056  int nrows;
1057  int ncols;
1058 
1059  *success = TRUE;
1060 
1061  nrows = SCIPgetNLPRows(scip);
1062  ncols = SCIPgetNLPCols(scip);
1063 
1064  /* check if we have the tableau row of the variable and if not compute it */
1065  for( v = 0; v < 4; ++v )
1066  {
1067  if( ! SCIPhashmapExists(tableau, (void*)vars[v]) )
1068  {
1069  SCIP_COL* col;
1070 
1071  /* get column of variable */
1072  col = SCIPvarGetCol(vars[v]);
1073 
1074  /* if variable is basic, then get its tableau row and insert it in the hashmap */
1076  {
1077  int lppos;
1078  SCIP_Real* densetableaurow;
1079 
1080  lppos = SCIPcolGetLPPos(col);
1081  SCIP_CALL( SCIPallocBufferArray(scip, &densetableaurow, ncols + nrows) );
1082 
1083  SCIP_CALL( SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], NULL, NULL) );
1084  SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], densetableaurow, NULL, NULL) );
1085 
1086  /* insert tableau row in hashmap*/
1087  SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)densetableaurow) );
1088  }
1089  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1090  {
1091  *success = FALSE;
1092  return SCIP_OKAY; /* don't even bother */
1093  }
1094  else
1095  {
1096  SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)NULL) );
1097  }
1098  }
1099 
1100  /* get tableau row of var */
1101  tableaurows[v] = (SCIP_Real *)SCIPhashmapGetImage(tableau, (void*)vars[v]);
1102  }
1103  return SCIP_OKAY;
1104 }
1105 
1106 /** computes the cut coefs of the non-basic (non-slack) variables (correspond to cols) and adds them to the
1107  * intersection cut
1108  */
1109 static
1111  SCIP* scip, /**< SCIP data structure */
1112  SCIP_VAR** vars, /**< variables */
1113  SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */
1114  SCIP_ROWPREP* rowprep, /**< store cut */
1115  SCIP_Real* rays, /**< buffer to store rays */
1116  int* nrays, /**< pointer to store number of nonzero rays */
1117  int* rayslppos, /**< buffer to store lppos of nonzero rays */
1118  SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */
1119  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1120  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1121  SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */
1122  )
1123 {
1124  int i;
1125  int ncols;
1126  SCIP_COL** cols;
1127 
1128  *success = TRUE;
1129 
1130  /* loop over non-basic (non-slack) variables */
1131  cols = SCIPgetLPCols(scip);
1132  ncols = SCIPgetNLPCols(scip);
1133  for( i = 0; i < ncols; ++i )
1134  {
1135  SCIP_COL* col;
1136  SCIP_Real coefs[5];
1137  SCIP_Real coefs4b[5];
1138  SCIP_Real coefscondition[3];
1139  SCIP_Real factor;
1140  SCIP_Bool israynonzero;
1141  SCIP_Real cutcoef;
1142  SCIP_Real interpoint;
1143  int v;
1144 
1145  col = cols[i];
1146 
1147  /* set factor to store entries of ray as = [-BinvL, BinvU] */
1149  factor = -1.0;
1150  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
1151  factor = 1.0;
1152  else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO )
1153  {
1154  *success = FALSE;
1155  return SCIP_OKAY;
1156  }
1157  else
1158  continue;
1159 
1160  /* build the ray */
1161  israynonzero = FALSE;
1162  for( v = 0; v < 4; ++v )
1163  {
1164  if( tableaurows[v] != NULL )
1165  rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][i]) ? 0.0 : tableaurows[v][i]);
1166  else
1167  {
1168  if( col == SCIPvarGetCol(vars[v]) )
1169  rays[(*nrays) * 4 + v] = -factor;
1170  else
1171  rays[(*nrays) * 4 + v] = 0.0;
1172  }
1173 
1174  israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1175  }
1176 
1177  /* do nothing if ray is 0 */
1178  if( ! israynonzero )
1179  continue;
1180 
1181  /* compute the cut */
1182  SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1183  ad, success) );
1184 
1185  if( *success == FALSE )
1186  return SCIP_OKAY;
1187 
1188  /* compute intersection point */
1189  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1190 
1191  /* store intersection points */
1192  interpoints[*nrays] = interpoint;
1193 
1194  /* remember lppos */
1195  rayslppos[*nrays] = i;
1196 
1197  /* count nonzero rays */
1198  *nrays += 1;
1199 
1200  /* compute cut coef */
1201  cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1202 
1203  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1205  SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ? -cutcoef :
1206  cutcoef, col) );
1207  }
1208 
1209  return SCIP_OKAY;
1210 }
1211 
1212 /** computes the cut coefs of the non-basic slack variables (correspond to rows) and adds them to the
1213  * intersection cut
1214  */
1215 static
1217  SCIP* scip, /**< SCIP data structure */
1218  SCIP_VAR** vars, /**< variables */
1219  SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */
1220  SCIP_ROWPREP* rowprep, /**< store cut */
1221  SCIP_Real* rays, /**< buffer to store rays */
1222  int* nrays, /**< pointer to store number of nonzero rays */
1223  int* rayslppos, /**< buffer to store lppos of nonzero rays */
1224  SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */
1225  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1226  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1227  SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */
1228  )
1229 {
1230  int i;
1231  int nrows;
1232  int ncols;
1233  SCIP_ROW** rows;
1234 
1235  nrows = SCIPgetNLPRows(scip);
1236  ncols = SCIPgetNLPCols(scip);
1237 
1238  *success = TRUE;
1239 
1240  /* loop over non-basic slack variables */
1241  rows = SCIPgetLPRows(scip);
1242  for( i = 0; i < nrows; ++i )
1243  {
1244  SCIP_ROW* row;
1245  SCIP_Real coefs[5];
1246  SCIP_Real coefs4b[5];
1247  SCIP_Real coefscondition[3];
1248  SCIP_Real factor;
1249  SCIP_Bool israynonzero;
1250  SCIP_Real cutcoef;
1251  SCIP_Real interpoint;
1252  int v;
1253 
1254  row = rows[i];
1255 
1256  /* set factor to store entries of ray as = [BinvL, -BinvU] */
1258  factor = 1.0;
1259  else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER )
1260  factor = -1.0;
1261  else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_ZERO )
1262  {
1263  *success = FALSE;
1264  return SCIP_OKAY;
1265  }
1266  else
1267  continue;
1268 
1269  /* build the ray */
1270  israynonzero = FALSE;
1271  for( v = 0; v < 4; ++v )
1272  {
1273  int idx;
1274 
1275  idx = ncols + i;
1276 
1277  if( tableaurows[v] != NULL )
1278  rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][idx]) ? 0.0 : tableaurows[v][idx]);
1279  else
1280  {
1281  /* TODO: We assume that slack variables can never occure in the minor. This is correct, right? */
1282  rays[(*nrays) * 4 + v] = 0.0;
1283  }
1284 
1285  israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0);
1286  }
1287 
1288  /* do nothing if ray is 0 */
1289  if( ! israynonzero )
1290  continue;
1291 
1292  /* compute the cut */
1293  SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds,
1294  ad, success) );
1295 
1296  if( *success == FALSE )
1297  return SCIP_OKAY;
1298 
1299  /* compute intersection point */
1300  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1301 
1302  /* store intersection points */
1303  interpoints[*nrays] = interpoint;
1304 
1305  /* store lppos of ray, make it negative so we can differentiate between cols and rows */
1306  rayslppos[*nrays] = -i - 1;
1307 
1308  /* count nonzero rays */
1309  *nrays += 1;
1310 
1311  /* compute cut coef */
1312  cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint;
1313 
1314  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1316 
1317  SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ? cutcoef :
1318  -cutcoef, row, success) ); /* rows have flipper base status! */
1319  }
1320 
1321  return SCIP_OKAY;
1322 }
1323 
1324 /* checks if two rays are linearly dependent */
1325 static
1327  SCIP* scip, /**< SCIP data structure */
1328  SCIP_Real* ray1, /**< coefficients of ray 1 */
1329  SCIP_Real* ray2, /**< coefficients of ray 2 */
1330  SCIP_Real* coef /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are dependent */
1331  )
1332 {
1333  int i;
1334 
1335  *coef = 0.0;
1336 
1337  for( i = 0; i < 4; ++i )
1338  {
1339  /* rays cannot be dependent if one ray has zero entry and the other one doesn't */
1340  if( (SCIPisZero(scip, ray1[i]) && ! SCIPisZero(scip, ray2[i])) ||
1341  (! SCIPisZero(scip, ray1[i]) && SCIPisZero(scip, ray2[i])) )
1342  {
1343  return FALSE;
1344  }
1345 
1346  if( *coef != 0.0 )
1347  {
1348  /* cannot be dependent if the coefs aren't equal for all entries */
1349  if( ! SCIPisFeasEQ(scip, *coef, ray1[i] / ray2[i]) )
1350  return FALSE;
1351  }
1352  else
1353  *coef = ray1[i] / ray2[i];
1354  }
1355 
1356  return TRUE;
1357 }
1358 
1359 /** finds the smallest negative steplength for the current ray r_idx such that the combination
1360  * of r_idx with all rays not in the recession cone is in the recession cone
1361  */
1362 static
1364  SCIP* scip, /**< SCIP data structure */
1365  SCIP_Real* rays, /**< rays */
1366  int nrays, /**< number of nonzero rays */
1367  int idx, /**< index of current ray we want to find rho for */
1368  SCIP_Real* interpoints, /**< intersection points of nonzero rays */
1369  SCIP_VAR** vars, /**< variables */
1370  SCIP_Real* rho, /**< pointer to store the optimal rho */
1371  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1372  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1373  SCIP_Bool* success /**< TRUE if computation of rho was successful */
1374  )
1375 {
1376  int i;
1377 
1378  *success = TRUE;
1379 
1380  /* go through all rays not in the recession cone and compute the largest negative steplength possible. The
1381  * smallest of them is then the steplength rho we use for the current ray */
1382  *rho = 0;
1383  for( i = 0; i < nrays; ++i )
1384  {
1385  SCIP_Real currentrho;
1386  SCIP_Real coef;
1387 
1388  if( SCIPisInfinity(scip, interpoints[i]) )
1389  continue;
1390 
1391  /* if the rays are linearly independent, we don't need to search for rho */
1392  if( raysAreDependent(scip, &rays[4 * i], &rays[4 * idx], &coef) )
1393  currentrho = coef * interpoints[i];
1394  else
1395  {
1396  SCIP_Real lb;
1397  SCIP_Real ub;
1398  SCIP_Real alpha;
1399  int j;
1400 
1401  /* do binary search by lookig at the convex combinations of r_i and r_j */
1402  lb = 0.0;
1403  ub = 1.0;
1404 
1405  for( j = 0; j < BINSEARCH_MAXITERS; ++j )
1406  {
1407  SCIP_Real coefs[5];
1408  SCIP_Real coefs4b[5];
1409  SCIP_Real coefscondition[3];
1410  SCIP_Real newray[4];
1411  SCIP_Real interpoint;
1412  int k;
1413 
1414  alpha = (lb + ub) / 2.0;
1415 
1416  /* build the ray alpha * ray_i + (1 - alpha) * ray_idx */
1417  for( k = 0; k < 4; ++k )
1418  newray[k] = alpha * rays[4 * i + k] - (1 - alpha) * rays[4 * idx + k];
1419 
1420  /* restrict phi to the "new" ray */
1421  SCIP_CALL( computeRestrictionToRay(scip, newray, vars, coefs, coefs4b, coefscondition, usebounds,
1422  ad, success) );
1423 
1424  if( ! *success )
1425  return SCIP_OKAY;
1426 
1427  /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is
1428  * positive
1429  */
1430 
1431  /* compute intersection point */
1432  interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition);
1433 
1434  /* no root exists */
1435  if( SCIPisInfinity(scip, interpoint) )
1436  {
1437  lb = alpha;
1438  if( SCIPisEQ(scip, ub, lb) )
1439  break;
1440  }
1441  else
1442  ub = alpha;
1443  }
1444 
1445  /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we
1446  * cannot move the ray in the recession cone, i.e. strengthening is not possible */
1447  if( SCIPisZero(scip, alpha) )
1448  {
1449  *rho = -SCIPinfinity(scip);
1450  return SCIP_OKAY;
1451  }
1452  else
1453  currentrho = (alpha - 1) * interpoints[i] / alpha;
1454  }
1455 
1456  if( currentrho < *rho )
1457  *rho = currentrho;
1458  }
1459 
1460  return SCIP_OKAY;
1461 }
1462 
1463 /** computes negative steplengths for the rays that are in the recession cone of the S-free set, i.e.,
1464  * which have an infinite intersection point.
1465  */
1466 static
1468  SCIP* scip, /**< SCIP data structure */
1469  SCIP_VAR** vars, /**< variables */
1470  SCIP_Real* rays, /**< rays */
1471  int nrays, /**< number of nonzero rays */
1472  int* rayslppos, /**< lppos of nonzero rays */
1473  SCIP_Real* interpoints, /**< intersection points */
1474  SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */
1475  SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */
1476  SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */
1477  SCIP_Bool* success /**< if a cut candidate could be computed */
1478  )
1479 {
1480  SCIP_COL** cols;
1481  SCIP_ROW** rows;
1482  int i;
1483 
1484  *success = TRUE;
1485 
1486  cols = SCIPgetLPCols(scip);
1487  rows = SCIPgetLPRows(scip);
1488 
1489  /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the
1490  * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */
1491  for( i = 0; i < nrays ; ++i )
1492  {
1493  SCIP_Real rho;
1494  SCIP_Real cutcoef;
1495  int lppos;
1496 
1497  if( !SCIPisInfinity(scip, interpoints[i]) )
1498  continue;
1499 
1500  /* compute the smallest rho */
1501  SCIP_CALL( findRho(scip, rays, nrays, i, interpoints, vars, &rho, usebounds, ad, success) );
1502 
1503  if( ! *success )
1504  continue;
1505 
1506  /* compute cut coef */
1507  cutcoef = SCIPisInfinity(scip, -rho) ? 0.0 : 1.0 / rho;
1508 
1509  /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */
1510  lppos = rayslppos[i];
1511  if( lppos < 0 )
1512  {
1513  lppos = -lppos - 1;
1514 
1515  assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) ==
1517 
1518  SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef :
1519  -cutcoef, rows[lppos], success) ); /* rows have flipped base status! */
1520 
1521  if( ! *success )
1522  return SCIP_OKAY;
1523  }
1524  else
1525  {
1526  assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) ==
1528  SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef :
1529  cutcoef, cols[lppos]) );
1530  }
1531  }
1532 
1533  return SCIP_OKAY;
1534 }
1535 
1536 /** separates cuts for stored principal minors */
1537 static
1539  SCIP* scip, /**< SCIP data structure */
1540  SCIP_SEPA* sepa, /**< separator */
1541  SCIP_SEPADATA* sepadata, /**< separator data */
1542  SCIP_VAR* xik, /**< variable X_ik = x_i * x_k */
1543  SCIP_VAR* xil, /**< variable X_il = x_i * x_l */
1544  SCIP_VAR* xjk, /**< variable X_jk = x_j * x_k */
1545  SCIP_VAR* xjl, /**< variable X_jl = x_j * x_l */
1546  SCIP_Bool* isxikdiag, /**< is X_ik diagonal? (i.e. i = k) */
1547  SCIP_Bool* isxildiag, /**< is X_il diagonal? (i.e. i = l) */
1548  SCIP_Bool* isxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */
1549  SCIP_Bool* isxjldiag, /**< is X_jl diagonal? (i.e. j = l) */
1550  int* basicvarpos2tableaurow,/**< map from basic var to its tableau row */
1551  SCIP_HASHMAP* tableau, /**< map from var to its tableau row */
1552  SCIP_RESULT* result /**< pointer to store the result of the separation call */
1553  )
1554 {
1555  SCIP_ROWPREP* rowprep;
1556  SCIP_VAR* vars[4] = {xik, xjl, xil, xjk};
1557  SCIP_Real* tableaurows[4];
1558  SCIP_Real* interpoints;
1559  SCIP_Real* rays;
1560  int nrays;
1561  int* rayslppos;
1562  int ncols;
1563  int nrows;
1564  SCIP_Bool success;
1565  SCIP_Real ad[4] = {0.0, 0.0, 0.0, 0.0};
1566  SCIP_Real solxik;
1567  SCIP_Real solxil;
1568  SCIP_Real solxjk;
1569  SCIP_Real solxjl;
1570 
1571  ncols = SCIPgetNLPCols(scip);
1572  nrows = SCIPgetNLPRows(scip);
1573 
1574  /* allocate memory for intersection points */
1575  SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, ncols + nrows) );
1576 
1577  /* allocate memory for rays */
1578  SCIP_CALL( SCIPallocBufferArray(scip, &rays, 4 * (ncols + nrows)) );
1579  SCIP_CALL( SCIPallocBufferArray(scip, &rayslppos, ncols + nrows) );
1580 
1581  /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */
1582  SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, TRUE) );
1583  SCIProwprepAddSide(rowprep, 1.0);
1584 
1585  /* check if we have the tableau row of the variable and if not compute it */
1586  SCIP_CALL( getTableauRows(scip, vars, basicvarpos2tableaurow, tableau, tableaurows, &success) );
1587 
1588  if( ! success )
1589  goto CLEANUP;
1590 
1591  /* if we want to enforce bounds, set the right a and d to enforce aTx + dTy <= 0 */
1592  if( sepadata->usebounds )
1593  {
1594  solxik = SCIPvarGetLPSol(xik);
1595  solxil = SCIPvarGetLPSol(xil);
1596  solxjk = SCIPvarGetLPSol(xjk);
1597  solxjl = SCIPvarGetLPSol(xjl);
1598 
1599  if( isxikdiag && SCIPisFeasNegative(scip, solxik) )
1600  {
1601  ad[0] = -1.0;
1602  ad[2] = 1.0;
1603  }
1604  else if( isxjldiag && SCIPisFeasNegative(scip, solxjl) )
1605  {
1606  ad[0] = -1.0;
1607  ad[2] = -1.0;
1608  }
1609  else if( isxildiag && SCIPisFeasNegative(scip, solxil) )
1610  {
1611  ad[1] = 1.0;
1612  ad[3] = -1.0;
1613  }
1614  else if( isxjkdiag && SCIPisFeasNegative(scip, solxjk) )
1615  {
1616  ad[1] = -1.0;
1617  ad[3] = -1.0;
1618  }
1619  }
1620 
1621  nrays = 0;
1622  /* loop over each non-basic var; get the ray; compute cut coefficient */
1623  SCIP_CALL( addCols(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1624 
1625  if( ! success )
1626  goto CLEANUP;
1627 
1628  /* loop over non-basic slack variables */
1629  SCIP_CALL( addRows(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) );
1630 
1631  if( ! success )
1632  goto CLEANUP;
1633 
1634  /* do strengthening */
1635  if( sepadata->usestrengthening )
1636  {
1637  SCIP_CALL( computeNegCutcoefs(scip, vars, rays, nrays, rayslppos, interpoints, rowprep, sepadata->usebounds, ad, &success) );
1638 
1639  if( ! success )
1640  goto CLEANUP;
1641  }
1642 
1643  /* merge coefficients that belong to same variable */
1644  SCIPmergeRowprepTerms(scip, rowprep);
1645 
1646  SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, sepadata->mincutviol, NULL, &success) );
1647 
1648  /* if cleanup was successfull, create row out of rowprep and add it */
1649  if( success )
1650  {
1651  SCIP_ROW* row;
1652  SCIP_Bool infeasible;
1653 
1654  /* create row */
1655  SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) );
1656 
1657  assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0);
1658 
1659  /* add row */
1660  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
1661 
1662  if( infeasible )
1663  *result = SCIP_CUTOFF;
1664  else
1665  *result = SCIP_SEPARATED;
1666 
1667  SCIP_CALL( SCIPreleaseRow(scip, &row) );
1668  }
1669 
1670 CLEANUP:
1671  SCIPfreeRowprep(scip, &rowprep);
1672  SCIPfreeBuffer(scip, &rayslppos);
1673  SCIPfreeBuffer(scip, &rays);
1674  SCIPfreeBuffer(scip, &interpoints);
1675 
1676  return SCIP_OKAY;
1677 }
1678 
1679 
1680 /** separates cuts for stored principal minors */
1681 static
1683  SCIP* scip, /**< SCIP data structure */
1684  SCIP_SEPA* sepa, /**< separator */
1685  SCIP_RESULT* result /**< pointer to store the result of the separation call */
1686  )
1687 {
1688  SCIP_SEPADATA* sepadata;
1689  SCIP_HASHMAP* tableau = NULL;
1690  int* basicvarpos2tableaurow = NULL; /* map between basic var and its tableau row */
1691  int i;
1692 
1693  assert(sepa != NULL);
1694  assert(result != NULL);
1695 
1696  *result = SCIP_DIDNOTRUN;
1697 
1698  sepadata = SCIPsepaGetData(sepa);
1699  assert(sepadata != NULL);
1700 
1701  /* check whether there are some minors available */
1702  if( sepadata->nminors == 0 )
1703  return SCIP_OKAY;
1704 
1705  *result = SCIP_DIDNOTFIND;
1706 
1707  /* loop over the minors and if they are violated build cut */
1708  for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i )
1709  {
1710  SCIP_VAR* auxvarxik;
1711  SCIP_VAR* auxvarxil;
1712  SCIP_VAR* auxvarxjk;
1713  SCIP_VAR* auxvarxjl;
1714  SCIP_Bool isauxvarxikdiag;
1715  SCIP_Bool isauxvarxildiag;
1716  SCIP_Bool isauxvarxjkdiag;
1717  SCIP_Bool isauxvarxjldiag;
1718  SCIP_Real solxik;
1719  SCIP_Real solxil;
1720  SCIP_Real solxjk;
1721  SCIP_Real solxjl;
1722  SCIP_Real det;
1723 
1724  /* get variables of the i-th minor */
1725  SCIP_CALL( getMinorVars(sepadata, i, &auxvarxik, &auxvarxil, &auxvarxjk, &auxvarxjl, &isauxvarxikdiag,
1726  &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag) );
1727 
1728  /* get current solution values */
1729  solxik = SCIPvarGetLPSol(auxvarxik);
1730  solxil = SCIPvarGetLPSol(auxvarxil);
1731  solxjk = SCIPvarGetLPSol(auxvarxjk);
1732  solxjl = SCIPvarGetLPSol(auxvarxjl);
1733 
1734  det = solxik * solxjl - solxil * solxjk;
1735 
1736  if( SCIPisFeasZero(scip, det) )
1737  continue;
1738 
1739  if( basicvarpos2tableaurow == NULL )
1740  {
1741  /* allocate memory */
1742  SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, SCIPgetNLPCols(scip)) );
1743  SCIP_CALL( SCIPhashmapCreate(&tableau, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1744 
1745  /* construct basicvar to tableau row map */
1746  SCIP_CALL( constructBasicVars2TableauRowMap(scip, basicvarpos2tableaurow) );
1747  }
1748  assert(tableau != NULL);
1749 
1750  if( SCIPisFeasPositive(scip, det) )
1751  {
1752  SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxik, auxvarxil, auxvarxjk, auxvarxjl, &isauxvarxikdiag,
1753  &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag, basicvarpos2tableaurow, tableau, result) );
1754  }
1755  else
1756  {
1757  assert(SCIPisFeasNegative(scip, det));
1758  SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxil, auxvarxik, auxvarxjl, auxvarxjk, &isauxvarxildiag,
1759  &isauxvarxikdiag, &isauxvarxjldiag, &isauxvarxjkdiag, basicvarpos2tableaurow, tableau, result) );
1760  }
1761  }
1762 
1763  /* all minors were feasible, so no memory to free */
1764  if( basicvarpos2tableaurow == NULL )
1765  return SCIP_OKAY;
1766 
1767  /* free memory */
1768  for( i = 0; i < SCIPhashmapGetNEntries(tableau); ++i )
1769  {
1770  SCIP_HASHMAPENTRY* entry;
1771 
1772  entry = SCIPhashmapGetEntry(tableau, i);
1773 
1774  if( entry != NULL )
1775  {
1776  SCIP_Real* tableaurow;
1777 
1778  tableaurow = (SCIP_Real *) SCIPhashmapEntryGetImage(entry);
1779 
1780  SCIPfreeBufferArrayNull(scip, &tableaurow);
1781  }
1782  }
1783  SCIPhashmapFree(&tableau);
1784  SCIPfreeBufferArray(scip, &basicvarpos2tableaurow);
1785 
1786  return SCIP_OKAY;
1787 }
1788 
1789 /*
1790  * Callback methods of separator
1791  */
1792 
1793 /** copy method for separator plugins (called when SCIP copies plugins) */
1794 static
1795 SCIP_DECL_SEPACOPY(sepaCopyMinor)
1796 { /*lint --e{715}*/
1797  assert(scip != NULL);
1798  assert(sepa != NULL);
1799  assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
1800 
1801  /* call inclusion method of constraint handler */
1803 
1804  return SCIP_OKAY;
1805 }
1806 
1807 
1808 /** destructor of separator to free user data (called when SCIP is exiting) */
1809 static
1810 SCIP_DECL_SEPAFREE(sepaFreeMinor)
1811 { /*lint --e{715}*/
1812  SCIP_SEPADATA* sepadata;
1813 
1814  sepadata = SCIPsepaGetData(sepa);
1815  assert(sepadata != NULL);
1816  assert(sepadata->minors == NULL);
1817  assert(sepadata->nminors == 0);
1818  assert(sepadata->minorssize == 0);
1819 
1820  /* free separator data */
1821  SCIPfreeBlockMemory(scip, &sepadata);
1822  SCIPsepaSetData(sepa, NULL);
1823 
1824  return SCIP_OKAY;
1825 }
1826 
1827 
1828 /** initialization method of separator (called after problem was transformed) */
1829 static
1830 SCIP_DECL_SEPAINIT(sepaInitMinor)
1831 { /*lint --e{715}*/
1832  SCIP_SEPADATA* sepadata;
1833 
1834  /* get separator data */
1835  sepadata = SCIPsepaGetData(sepa);
1836  assert(sepadata != NULL);
1837  assert(sepadata->randnumgen == NULL);
1838 
1839  /* create random number generator */
1840  SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
1841 
1842  return SCIP_OKAY;
1843 }
1844 
1845 
1846 /** deinitialization method of separator (called before transformed problem is freed) */
1847 static
1848 SCIP_DECL_SEPAEXIT(sepaExitMinor)
1849 { /*lint --e{715}*/
1850  SCIP_SEPADATA* sepadata;
1851 
1852  /* get separator data */
1853  sepadata = SCIPsepaGetData(sepa);
1854  assert(sepadata != NULL);
1855  assert(sepadata->randnumgen != NULL);
1856 
1857  /* free random number generator */
1858  SCIPfreeRandom(scip, &sepadata->randnumgen);
1859 
1860  return SCIP_OKAY;
1861 }
1862 
1863 
1864 /** solving process initialization method of separator (called when branch and bound process is about to begin) */
1865 static
1866 SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
1867 { /*lint --e{715}*/
1868  return SCIP_OKAY;
1869 }
1870 
1871 
1872 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */
1873 static
1874 SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
1875 { /*lint --e{715}*/
1876  SCIP_SEPADATA* sepadata;
1877 
1878  sepadata = SCIPsepaGetData(sepa);
1879  assert(sepadata != NULL);
1880 
1881  /* clear separation data */
1882  SCIP_CALL( sepadataClear(scip, sepadata) );
1883 
1884  return SCIP_OKAY;
1885 }
1886 
1887 
1888 /** LP solution separation method of separator */
1889 static
1890 SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
1891 { /*lint --e{715}*/
1892  SCIP_SEPADATA* sepadata;
1893  int ncalls;
1894  int currentdepth;
1895 
1896  /* need routine to compute eigenvalues/eigenvectors */
1897  if( !SCIPisIpoptAvailableIpopt() )
1898  return SCIP_OKAY;
1899 
1900  sepadata = SCIPsepaGetData(sepa);
1901  assert(sepadata != NULL);
1902  currentdepth = SCIPgetDepth(scip);
1903  ncalls = SCIPsepaGetNCallsAtNode(sepa);
1904 
1905  /* only call the separator a given number of times at each node */
1906  if( (currentdepth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
1907  || (currentdepth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
1908  {
1909  SCIPdebugMsg(scip, "reached round limit for node\n");
1910  return SCIP_OKAY;
1911  }
1912 
1913  /* try to detect minors */
1914  SCIP_CALL( detectMinors(scip, sepadata) );
1915 
1916  /* call separation method */
1917  SCIP_CALL( separatePoint(scip, sepa, result) );
1918 
1919  return SCIP_OKAY;
1920 }
1921 
1922 
1923 /*
1924  * separator specific interface methods
1925  */
1926 
1927 /** creates the minor separator and includes it in SCIP */
1929  SCIP* scip /**< SCIP data structure */
1930  )
1931 {
1932  SCIP_SEPADATA* sepadata = NULL;
1933  SCIP_SEPA* sepa = NULL;
1934 
1935  /* create minor separator data */
1936  SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
1937  BMSclearMemory(sepadata);
1938 
1939  /* include separator */
1942  sepaExeclpMinor, NULL,
1943  sepadata) );
1944 
1945  assert(sepa != NULL);
1946 
1947  /* set non fundamental callbacks via setter functions */
1948  SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) );
1949  SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) );
1950  SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) );
1951  SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) );
1952  SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) );
1953  SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) );
1954 
1955  /* add minor separator parameters */
1957  "separating/" SEPA_NAME "/usestrengthening",
1958  "whether to use strengthened intersection cuts to separate minors",
1959  &sepadata->usestrengthening, FALSE, DEFAULT_USESTRENGTHENING, NULL, NULL) );
1960 
1962  "separating/" SEPA_NAME "/usebounds",
1963  "whether to also enforce nonegativity bounds of principle minors",
1964  &sepadata->usebounds, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) );
1965 
1967  "separating/" SEPA_NAME "/mincutviol",
1968  "minimum required violation of a cut",
1969  &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
1970 
1971  SCIP_CALL( SCIPaddIntParam(scip,
1972  "separating/" SEPA_NAME "/maxrounds",
1973  "maximal number of separation rounds per node (-1: unlimited)",
1974  &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
1975 
1976  SCIP_CALL( SCIPaddIntParam(scip,
1977  "separating/" SEPA_NAME "/maxroundsroot",
1978  "maximal number of separation rounds in the root node (-1: unlimited)",
1979  &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
1980 
1981  return SCIP_OKAY;
1982 }
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:3570
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)
#define NULL
Definition: def.h:267
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:501
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:941
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3854
#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:17031
int SCIPcalcMemGrowSize(SCIP *scip, int num)
Definition: scip_mem.c:139
#define SQR(x)
Definition: def.h:214
SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:18135
static SCIP_DECL_SEPACOPY(sepaCopyMinor)
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1250
#define DEFAULT_RANDSEED
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
int SCIProwGetNLPNonz(SCIP_ROW *row)
Definition: lp.c:17227
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4595
private functions to work with algebraic expressions
SCIP_Real SCIProwGetLhs(SCIP_ROW *row)
Definition: lp.c:17292
#define FALSE
Definition: def.h:94
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition: misc.c:3074
SCIP_RETCODE SCIPgetRowprepRowSepa(SCIP *scip, SCIP_ROW **row, SCIP_ROWPREP *rowprep, SCIP_SEPA *sepa)
SCIP_BASESTAT SCIProwGetBasisStatus(SCIP_ROW *row)
Definition: lp.c:17340
static SCIP_Real computeRoot(SCIP *scip, SCIP_Real *coefs)
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3457
#define TRUE
Definition: def.h:93
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:17769
#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:3261
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:1464
#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:1067
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3423
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3864
#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:16996
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:3541
SCIP_HASHMAPENTRY * SCIPhashmapGetEntry(SCIP_HASHMAP *hashmap, int entryidx)
Definition: misc.c:3549
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:8216
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17420
static SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
void SCIPcomputeArraysIntersectionInt(int *array1, int narray1, int *array2, int narray2, int *intersectarray, int *nintersectarray)
Definition: misc.c:10559
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition: misc.c:3108
void SCIPsepaSetData(SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata)
Definition: sepa.c:643
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
int * vals
power and signed power expression handlers
SCIP_Real SCIPvarGetLPSol(SCIP_VAR *var)
Definition: var.c:18453
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:583
SCIP_HASHMAP * auxvars
#define SCIP_CALL(x)
Definition: def.h:380
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
Definition: misc_rowprep.c:746
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:17302
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:195
SCIP_COL ** SCIProwGetCols(SCIP_ROW *row)
Definition: lp.c:17238
int SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4638
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2337
#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:17248
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:91
SCIP_EXPR * SCIPexpriterRestartDFS(SCIP_EXPRITER *iterator, SCIP_EXPR *expr)
Definition: expriter.c:630
void SCIProwprepAddConstant(SCIP_ROWPREP *rowprep, SCIP_Real constant)
Definition: misc_rowprep.c:760
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
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:858
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:17790
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:1475
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
#define BMSclearMemory(ptr)
Definition: memory.h:129
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:67
static SCIP_DECL_SEPAEXECLP(sepaExeclpMinor)
int SCIPgetNVars(SCIP *scip)
Definition: scip_prob.c:1992
#define SCIP_REAL_MAX
Definition: def.h:174
void SCIPfreeExpriter(SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2351
product expression handler
SCIP_Real SCIProwGetConstant(SCIP_ROW *row)
Definition: lp.c:17258
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
#define MAX(x, y)
Definition: def.h:239
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:17042
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:1947
int SCIProwGetLPPos(SCIP_ROW *row)
Definition: lp.c:17501
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
SCIP_RETCODE SCIPcaptureVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_var.c:1216
#define SCIP_Real
Definition: def.h:173
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:563
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:18145
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:913
#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:3156
SCIP_Bool SCIPexpriterIsEnd(SCIP_EXPRITER *iterator)
Definition: expriter.c:969
static SCIP_RETCODE constructBasicVars2TableauRowMap(SCIP *scip, int *map)
static SCIP_DECL_SEPAEXIT(sepaExitMinor)
int SCIPcolGetLPPos(SCIP_COL *col)
Definition: lp.c:17093
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
#define ABS(x)
Definition: def.h:235
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:2104
#define SCIPreallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:128