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
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 */
71struct 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 */
87struct 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 */
101static
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 */
158static
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 */
189static
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 */
230static
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 */
271 arr->valssize = 10;
272 arr->nvals = 0;
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 */
294static
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 */
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 }
488 SCIPhashmapFree(&rowi->auxvars);
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 */
509static
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 */
538static
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}*/
699static
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 */
756static
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 */
798static
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 */
816static
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
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 */
903static
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 */
952static
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 */
984static
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)));
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)));
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 */
1045static
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 */
1109static
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;
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 */
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 */
1215static
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;
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
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 */
1325static
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 */
1362static
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 */
1466static
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 */
1537static
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 */
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
1670CLEANUP:
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 */
1681static
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)) );
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) */
1794static
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) */
1809static
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) */
1829static
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) */
1847static
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) */
1865static
1866SCIP_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) */
1873static
1874SCIP_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 */
1889static
1890SCIP_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 */
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
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
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}
SCIP_VAR * a
Definition: circlepacking.c:66
SCIP_VAR ** b
Definition: circlepacking.c:65
constraint handler for nonlinear constraints specified by algebraic expressions
#define SCIPquadprecSqrtQ(r, a)
Definition: dbldblarith.h:71
#define SCIPquadprecProdDD(r, a, b)
Definition: dbldblarith.h:58
#define SCIPquadprecProdQD(r, a, b)
Definition: dbldblarith.h:63
#define QUAD_SCALE(x, a)
Definition: dbldblarith.h:50
#define SCIPquadprecSumQD(r, a, b)
Definition: dbldblarith.h:62
#define QUAD_ASSIGN(a, constant)
Definition: dbldblarith.h:51
#define QUAD(x)
Definition: dbldblarith.h:47
#define SCIPquadprecSquareD(r, a)
Definition: dbldblarith.h:59
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:67
#define QUAD_TO_DBL(x)
Definition: dbldblarith.h:49
#define NULL
Definition: def.h:267
#define SCIP_REAL_MAX
Definition: def.h:174
#define SCIP_INTERVAL_INFINITY
Definition: def.h:195
#define SCIP_Bool
Definition: def.h:91
#define SCIP_Real
Definition: def.h:173
#define ABS(x)
Definition: def.h:235
#define SQR(x)
Definition: def.h:214
#define TRUE
Definition: def.h:93
#define FALSE
Definition: def.h:94
#define MAX(x, y)
Definition: def.h:239
#define SCIP_CALL(x)
Definition: def.h:374
private functions to work with algebraic expressions
power and signed power expression handlers
product expression handler
variable expression handler
void SCIPcomputeArraysIntersectionInt(int *array1, int narray1, int *array2, int narray2, int *intersectarray, int *nintersectarray)
Definition: misc.c:10559
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
const char * SCIPgetProbName(SCIP *scip)
Definition: scip_prob.c:1067
int SCIPgetNVars(SCIP *scip)
Definition: scip_prob.c:1992
SCIP_VAR ** SCIPgetVars(SCIP *scip)
Definition: scip_prob.c:1947
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition: misc.c:3108
void * SCIPhashmapEntryGetImage(SCIP_HASHMAPENTRY *entry)
Definition: misc.c:3570
void * SCIPhashmapGetImage(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3261
SCIP_RETCODE SCIPhashmapInsert(SCIP_HASHMAP *hashmap, void *origin, void *image)
Definition: misc.c:3156
int SCIPhashmapGetNEntries(SCIP_HASHMAP *hashmap)
Definition: misc.c:3541
SCIP_HASHMAPENTRY * SCIPhashmapGetEntry(SCIP_HASHMAP *hashmap, int entryidx)
Definition: misc.c:3549
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition: misc.c:3074
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3423
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
#define SCIPdebugMsg
Definition: scip_message.h:78
SCIP_Bool SCIPisIpoptAvailableIpopt(void)
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
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
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
int SCIPcolGetLPPos(SCIP_COL *col)
Definition: lp.c:17093
SCIP_VAR * SCIPcolGetVar(SCIP_COL *col)
Definition: lp.c:17042
SCIP_Real SCIPcolGetPrimsol(SCIP_COL *col)
Definition: lp.c:16996
SCIP_BASESTAT SCIPcolGetBasisStatus(SCIP_COL *col)
Definition: lp.c:17031
int SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4636
SCIP_CONSHDLR * SCIPfindConshdlr(SCIP *scip, const char *name)
Definition: scip_cons.c:941
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4593
const char * SCIPconsGetName(SCIP_CONS *cons)
Definition: cons.c:8214
SCIP_Real SCIPgetCutEfficacy(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:94
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:250
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3860
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3456
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1464
SCIP_Bool SCIPexpriterIsEnd(SCIP_EXPRITER *iterator)
Definition: expriter.c:969
SCIP_EXPR * SCIPexpriterRestartDFS(SCIP_EXPRITER *iterator, SCIP_EXPR *expr)
Definition: expriter.c:630
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2337
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1475
SCIP_EXPR * SCIPexpriterGetNext(SCIP_EXPRITER *iterator)
Definition: expriter.c:858
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3870
void SCIPfreeExpriter(SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2351
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:501
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_Real sqrcoeff, SCIP_Real lincoeff, SCIP_Real rhs, SCIP_INTERVAL xbnds)
SCIP_RETCODE SCIPgetLPBasisInd(SCIP *scip, int *basisind)
Definition: scip_lp.c:686
SCIP_ROW ** SCIPgetLPRows(SCIP *scip)
Definition: scip_lp.c:605
int SCIPgetNLPRows(SCIP *scip)
Definition: scip_lp.c:626
SCIP_RETCODE SCIPgetLPBInvARow(SCIP *scip, int r, SCIP_Real *binvrow, SCIP_Real *coefs, int *inds, int *ninds)
Definition: scip_lp.c:785
SCIP_COL ** SCIPgetLPCols(SCIP *scip)
Definition: scip_lp.c:506
int SCIPgetNLPCols(SCIP *scip)
Definition: scip_lp.c:527
SCIP_RETCODE SCIPgetLPBInvRow(SCIP *scip, int r, SCIP_Real *coefs, int *inds, int *ninds)
Definition: scip_lp.c:714
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:134
int SCIPcalcMemGrowSize(SCIP *scip, int num)
Definition: scip_mem.c:139
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:124
#define SCIPreallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:128
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:136
#define SCIPallocBuffer(scip, ptr)
Definition: scip_mem.h:122
#define SCIPreallocBlockMemoryArray(scip, ptr, oldnum, newnum)
Definition: scip_mem.h:99
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:108
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:111
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:137
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:89
SCIP_Real SCIProwGetLhs(SCIP_ROW *row)
Definition: lp.c:17292
SCIP_COL ** SCIProwGetCols(SCIP_ROW *row)
Definition: lp.c:17238
SCIP_Real SCIProwGetRhs(SCIP_ROW *row)
Definition: lp.c:17302
int SCIProwGetNLPNonz(SCIP_ROW *row)
Definition: lp.c:17227
int SCIProwGetLPPos(SCIP_ROW *row)
Definition: lp.c:17501
SCIP_Real SCIPgetRowActivity(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:2104
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1562
SCIP_Real SCIProwGetConstant(SCIP_ROW *row)
Definition: lp.c:17258
SCIP_Real * SCIProwGetVals(SCIP_ROW *row)
Definition: lp.c:17248
SCIP_BASESTAT SCIProwGetBasisStatus(SCIP_ROW *row)
Definition: lp.c:17340
SCIP_RETCODE SCIPsetSepaExit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAEXIT((*sepaexit)))
Definition: scip_sepa.c:199
SCIP_RETCODE SCIPsetSepaInitsol(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAINITSOL((*sepainitsol)))
Definition: scip_sepa.c:215
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
const char * SCIPsepaGetName(SCIP_SEPA *sepa)
Definition: sepa.c:743
int SCIPsepaGetNCallsAtNode(SCIP_SEPA *sepa)
Definition: sepa.c:880
SCIP_RETCODE SCIPsetSepaFree(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAFREE((*sepafree)))
Definition: scip_sepa.c:167
SCIP_RETCODE SCIPsetSepaExitsol(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAEXITSOL((*sepaexitsol)))
Definition: scip_sepa.c:231
SCIP_SEPADATA * SCIPsepaGetData(SCIP_SEPA *sepa)
Definition: sepa.c:633
void SCIPsepaSetData(SCIP_SEPA *sepa, SCIP_SEPADATA *sepadata)
Definition: sepa.c:643
SCIP_RETCODE SCIPsetSepaCopy(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPACOPY((*sepacopy)))
Definition: scip_sepa.c:151
SCIP_RETCODE SCIPsetSepaInit(SCIP *scip, SCIP_SEPA *sepa, SCIP_DECL_SEPAINIT((*sepainit)))
Definition: scip_sepa.c:183
SCIP_Real SCIPgetTotalTime(SCIP *scip)
Definition: scip_timing.c:351
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisFeasZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisHugeValue(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisFeasPositive(SCIP *scip, SCIP_Real val)
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:670
SCIP_COL * SCIPvarGetCol(SCIP_VAR *var)
Definition: var.c:17789
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:18144
int SCIPvarGetProbindex(SCIP_VAR *var)
Definition: var.c:17768
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17419
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1248
SCIP_Real SCIPvarGetLPSol(SCIP_VAR *var)
Definition: var.c:18452
SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:18134
SCIP_RETCODE SCIPcaptureVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_var.c:1214
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
void SCIPmergeRowprepTerms(SCIP *scip, SCIP_ROWPREP *rowprep)
void SCIProwprepAddConstant(SCIP_ROWPREP *rowprep, SCIP_Real constant)
Definition: misc_rowprep.c:760
SCIP_RETCODE SCIPgetRowprepRowSepa(SCIP *scip, SCIP_ROW **row, SCIP_ROWPREP *rowprep, SCIP_SEPA *sepa)
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:913
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
Definition: misc_rowprep.c:563
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
Definition: misc_rowprep.c:746
SCIP_RETCODE SCIPcleanupRowprep(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_SOL *sol, SCIP_Real minviol, SCIP_Real *viol, SCIP_Bool *success)
void SCIPfreeRowprep(SCIP *scip, SCIP_ROWPREP **rowprep)
Definition: misc_rowprep.c:583
SCIP_RETCODE SCIPincludeSepaInterminor(SCIP *scip)
void SCIPsortInt(int *intarray, int len)
#define BMSclearMemory(ptr)
Definition: memory.h:129
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:57
Ipopt NLP interface.
#define SCIPstatisticMessage
Definition: pub_message.h:123
static SCIP_RETCODE insertIndex(SCIP *scip, SCIP_HASHMAP *rowmap, SCIP_VAR *row, SCIP_VAR *col, SCIP_VAR *auxvar, int *rowindices, int *nrows)
#define SEPA_PRIORITY
static SCIP_Real computeIntersectionPoint(SCIP *scip, SCIP_Bool usebounds, SCIP_Real *coefs, SCIP_Real *coefs4b, SCIP_Real *coefscondition)
#define DEFAULT_USEBOUNDS
static SCIP_DECL_SEPAINIT(sepaInitMinor)
#define DEFAULT_USESTRENGTHENING
#define SEPA_DELAY
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)
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_DECL_SEPAEXECLP(sepaExeclpMinor)
static SCIP_RETCODE sepadataClear(SCIP *scip, SCIP_SEPADATA *sepadata)
#define SEPA_DESC
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)
#define DEFAULT_MAXROUNDSROOT
#define SEPA_USESSUBSCIP
static SCIP_RETCODE getTableauRows(SCIP *scip, SCIP_VAR **vars, int *basicvarpos2tableaurow, SCIP_HASHMAP *tableau, SCIP_Real **tableaurows, SCIP_Bool *success)
static SCIP_DECL_SEPAFREE(sepaFreeMinor)
static SCIP_DECL_SEPACOPY(sepaCopyMinor)
static SCIP_Real isCase4a(SCIP_Real tsol, SCIP_Real *coefs, SCIP_Real *coefscondition)
static SCIP_RETCODE addRowToCut(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_Real cutcoef, SCIP_ROW *row, SCIP_Bool *success)
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)
static SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor)
static SCIP_Real computeRoot(SCIP *scip, SCIP_Real *coefs)
static SCIP_Real evalPhiAtRay(SCIP *scip, SCIP_Real t, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e)
#define DEFAULT_MINCUTVIOL
static SCIP_RETCODE detectMinors(SCIP *scip, SCIP_SEPADATA *sepadata)
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 SEPA_MAXBOUNDDIST
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)
static SCIP_DECL_SEPAEXIT(sepaExitMinor)
#define SEPA_FREQ
#define DEFAULT_RANDSEED
#define SEPA_NAME
static SCIP_DECL_SEPAINITSOL(sepaInitsolMinor)
#define BINSEARCH_MAXITERS
#define DEFAULT_MAXROUNDS
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)
static void doBinarySearch(SCIP *scip, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e, SCIP_Real *sol)
static SCIP_RETCODE addColToCut(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_Real cutcoef, SCIP_COL *col)
static SCIP_RETCODE constructBasicVars2TableauRowMap(SCIP *scip, int *map)
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)
static SCIP_Bool raysAreDependent(SCIP *scip, SCIP_Real *ray1, SCIP_Real *ray2, SCIP_Real *coef)
static SCIP_RETCODE separatePoint(SCIP *scip, SCIP_SEPA *sepa, SCIP_RESULT *result)
int * vals
SCIP_HASHMAP * auxvars
@ SCIP_EXPRITER_DFS
Definition: type_expr.h:716
@ SCIP_SIDETYPE_LEFT
Definition: type_lp.h:64
@ SCIP_BASESTAT_BASIC
Definition: type_lpi.h:92
@ SCIP_BASESTAT_UPPER
Definition: type_lpi.h:93
@ SCIP_BASESTAT_LOWER
Definition: type_lpi.h:91
@ SCIP_BASESTAT_ZERO
Definition: type_lpi.h:94
@ SCIP_DIDNOTRUN
Definition: type_result.h:42
@ SCIP_CUTOFF
Definition: type_result.h:48
@ SCIP_DIDNOTFIND
Definition: type_result.h:44
@ SCIP_SEPARATED
Definition: type_result.h:49
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:61
@ SCIP_OKAY
Definition: type_retcode.h:42
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
struct SCIP_SepaData SCIP_SEPADATA
Definition: type_sepa.h:52