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