1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file sepa_disjunctive.c
17 * @ingroup DEFPLUGINS_SEPA
18 * @brief disjunctive cut separator
19 * @author Tobias Fischer
20 * @author Marc Pfetsch
21 *
22 * We separate disjunctive cuts for two term disjunctions of the form \f$x_1 = 0 \vee x_2 = 0\f$. They can be generated
23 * directly from the simplex tableau. For further information, we refer to@n
24 * "A complementarity-based partitioning and disjunctive cut algorithm for mathematical programming problems with
25 * equilibrium constraints"@n
26 * Júdice, J.J., Sherali, H.D., Ribeiro, I.M., Faustino, A.M., Journal of Global Optimization 36(1), 89–114 (2006)
27 *
28 * Cut coefficients belonging to integer variables can be strengthened by the 'monoidal cut strengthening' procedure, see@n
29 * "Strengthening cuts for mixed integer programs"@n
30 * Egon Balas, Robert G. Jeroslow, European Journal of Operational Research, Volume 4, Issue 4, 1980, Pages 224-234
31 */
32
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34
35 #include "blockmemshell/memory.h"
36 #include "scip/cons_sos1.h"
37 #include "scip/pub_cons.h"
38 #include "scip/pub_lp.h"
39 #include "scip/pub_message.h"
40 #include "scip/pub_misc.h"
41 #include "scip/pub_misc_sort.h"
42 #include "scip/pub_sepa.h"
43 #include "scip/pub_var.h"
44 #include "scip/scip_cons.h"
45 #include "scip/scip_cut.h"
46 #include "scip/scip_general.h"
47 #include "scip/scip_lp.h"
48 #include "scip/scip_mem.h"
49 #include "scip/scip_message.h"
50 #include "scip/scip_numerics.h"
51 #include "scip/scip_param.h"
52 #include "scip/scip_sepa.h"
53 #include "scip/scip_sol.h"
54 #include "scip/scip_solvingstats.h"
55 #include "scip/scip_tree.h"
56 #include "scip/sepa_disjunctive.h"
57 #include <string.h>
58
59
60 #define SEPA_NAME "disjunctive"
61 #define SEPA_DESC "disjunctive cut separator"
62 #define SEPA_PRIORITY 10 /**< priority for separation */
63 #define SEPA_FREQ 0 /**< frequency for separating cuts; zero means to separate only in the root node */
64 #define SEPA_MAXBOUNDDIST 0.0 /**< maximal relative distance from the current node's dual bound to primal bound
65 * compared to best node's dual bound for applying separation.*/
66 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
67 #define SEPA_DELAY TRUE /**< should separation method be delayed, if other separators found cuts? */
68
69 #define DEFAULT_MAXRANK 20 /**< maximal rank of a cut that could not be scaled to integral coefficients (-1: unlimited) */
70 #define DEFAULT_MAXRANKINTEGRAL -1 /**< maximal rank of a cut that could be scaled to integral coefficients (-1: unlimited) */
71 #define DEFAULT_MAXWEIGHTRANGE 1e+03 /**< maximal valid range max(|weights|)/min(|weights|) of row weights */
72 #define DEFAULT_STRENGTHEN TRUE /**< strengthen cut if integer variables are present */
73
74 #define DEFAULT_MAXDEPTH -1 /**< node depth of separating cuts (-1: no limit) */
75 #define DEFAULT_MAXROUNDS 25 /**< maximal number of separation rounds in a branching node (-1: no limit) */
76 #define DEFAULT_MAXROUNDSROOT 100 /**< maximal number of separation rounds in the root node (-1: no limit) */
77 #define DEFAULT_MAXINVCUTS 50 /**< maximal number of cuts investigated per iteration in a branching node */
78 #define DEFAULT_MAXINVCUTSROOT 250 /**< maximal number of cuts investigated per iteration in the root node */
79 #define DEFAULT_MAXCONFSDELAY 100000 /**< delay separation if number of conflict graph edges is larger than predefined value (-1: no limit) */
80
81 #define MAKECONTINTEGRAL FALSE /**< convert continuous variable to integral variables in SCIPmakeRowIntegral() */
82
83
84 /** separator data */
85 struct SCIP_SepaData
86 {
87 SCIP_Bool strengthen; /**< strengthen cut if integer variables are present */
88 SCIP_CONSHDLR* conshdlr; /**< SOS1 constraint handler */
89 SCIP_Real maxweightrange; /**< maximal valid range max(|weights|)/min(|weights|) of row weights */
90 int maxrank; /**< maximal rank of a cut that could not be scaled to integral coefficients (-1: unlimited) */
91 int maxrankintegral; /**< maximal rank of a cut that could be scaled to integral coefficients (-1: unlimited) */
92 int maxdepth; /**< node depth of separating cuts (-1: no limit) */
93 int maxrounds; /**< maximal number of separation rounds in a branching node (-1: no limit) */
94 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: no limit) */
95 int maxinvcuts; /**< maximal number of cuts separated per iteration in a branching node */
96 int maxinvcutsroot; /**< maximal number of cuts separated per iteration in the root node */
97 int maxconfsdelay; /**< delay separation if number of conflict graph edges is larger than predefined value (-1: no limit) */
98 int lastncutsfound; /**< total number of cuts found after last call of separator */
99 };
100
101
102 /** gets rank of variable corresponding to row of \f$B^{-1}\f$ */
103 static
getVarRank(SCIP * scip,SCIP_Real * binvrow,SCIP_Real * rowsmaxval,SCIP_Real maxweightrange,SCIP_ROW ** rows,int nrows)104 int getVarRank(
105 SCIP* scip, /**< SCIP pointer */
106 SCIP_Real* binvrow, /**< row of \f$B^{-1}\f$ */
107 SCIP_Real* rowsmaxval, /**< maximal row multiplicator from nonbasic matrix A_N */
108 SCIP_Real maxweightrange, /**< maximal valid range max(|weights|)/min(|weights|) of row weights */
109 SCIP_ROW** rows, /**< rows of LP relaxation of scip */
110 int nrows /**< number of rows */
111 )
112 {
113 SCIP_Real maxweight = 0.0;
114 int maxrank = 0;
115 int r;
116
117 assert( scip != NULL );
118 assert( binvrow != NULL || nrows == 0 );
119 assert( rowsmaxval != NULL || nrows == 0 );
120 assert( rows != NULL || nrows == 0 );
121
122 /* compute maximum row weights resulting from multiplication */
123 for (r = 0; r < nrows; ++r)
124 {
125 SCIP_Real val;
126
127 val = REALABS(binvrow[r] * rowsmaxval[r]);/*lint !e613*/
128 if ( SCIPisGT(scip, val, maxweight) )
129 maxweight = val;
130 }
131
132 /* compute rank */
133 for (r = 0; r < nrows; ++r)
134 {
135 SCIP_Real val;
136 int rank;
137
138 val = REALABS(binvrow[r] * rowsmaxval[r]);/*lint !e613*/
139 rank = SCIProwGetRank(rows[r]);/*lint !e613*/
140 if ( rank > maxrank && SCIPisGT(scip, val * maxweightrange, maxweight) )
141 maxrank = rank;
142 }
143
144 return maxrank;
145 }
146
147
148 /** gets the nonbasic coefficients of a simplex row */
149 static
getSimplexCoefficients(SCIP * scip,SCIP_ROW ** rows,int nrows,SCIP_COL ** cols,int ncols,SCIP_Real * coef,SCIP_Real * binvrow,SCIP_Real * simplexcoefs,int * nonbasicnumber)150 SCIP_RETCODE getSimplexCoefficients(
151 SCIP* scip, /**< SCIP pointer */
152 SCIP_ROW** rows, /**< LP rows */
153 int nrows, /**< number LP rows */
154 SCIP_COL** cols, /**< LP columns */
155 int ncols, /**< number of LP columns */
156 SCIP_Real* coef, /**< row of \f$B^{-1} \cdot A\f$ */
157 SCIP_Real* binvrow, /**< row of \f$B^{-1}\f$ */
158 SCIP_Real* simplexcoefs, /**< pointer to store the nonbasic simplex-coefficients */
159 int* nonbasicnumber /**< pointer to store the number of nonbasic simplex-coefficients */
160 )
161 {
162 int r;
163 int c;
164
165 assert( scip != NULL );
166 assert( rows != NULL );
167 assert( nonbasicnumber != NULL );
168 assert( simplexcoefs != NULL );
169 assert( cols != NULL );
170
171 *nonbasicnumber = 0;
172
173 /* note: the non-slack variables have to be added first (see the function generateDisjCutSOS1()) */
174
175 /* get simplex-coefficients of the non-basic non-slack variables */
176 for (c = 0; c < ncols; ++c)
177 {
178 SCIP_COL* col;
179
180 col = cols[c];
181 assert( col != NULL );
182 if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER || SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
183 simplexcoefs[(*nonbasicnumber)++] = coef[c];
184 }
185
186 /* get simplex-coefficients of the non-basic slack variables */
187 for (r = 0; r < nrows; ++r)
188 {
189 SCIP_ROW* row;
190 row = rows[r];
191 assert( row != NULL );
192
193 if ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER )
194 {
195 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, row) - SCIProwGetRhs(row)) || SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, row) - SCIProwGetLhs(row)) );
196
197 simplexcoefs[(*nonbasicnumber)++] = binvrow[r];
198 }
199 }
200
201 return SCIP_OKAY;
202 }
203
204
205 /** computes a disjunctive cut inequality based on two simplex taubleau rows */
206 static
generateDisjCutSOS1(SCIP * scip,SCIP_SEPA * sepa,SCIP_ROW ** rows,int nrows,SCIP_COL ** cols,int ncols,int ndisjcuts,SCIP_Bool scale,SCIP_Bool strengthen,SCIP_Real cutlhs1,SCIP_Real cutlhs2,SCIP_Real bound1,SCIP_Real bound2,SCIP_Real * simplexcoefs1,SCIP_Real * simplexcoefs2,SCIP_Real * cutcoefs,SCIP_ROW ** row,SCIP_Bool * madeintegral)207 SCIP_RETCODE generateDisjCutSOS1(
208 SCIP* scip, /**< SCIP pointer */
209 SCIP_SEPA* sepa, /**< separator */
210 SCIP_ROW** rows, /**< LP rows */
211 int nrows, /**< number of LP rows */
212 SCIP_COL** cols, /**< LP columns */
213 int ncols, /**< number of LP columns */
214 int ndisjcuts, /**< number of disjunctive cuts found so far */
215 SCIP_Bool scale, /**< should cut be scaled */
216 SCIP_Bool strengthen, /**< should cut be strengthened if integer variables are present */
217 SCIP_Real cutlhs1, /**< left hand side of the first simplex row */
218 SCIP_Real cutlhs2, /**< left hand side of the second simplex row */
219 SCIP_Real bound1, /**< bound of first simplex row */
220 SCIP_Real bound2, /**< bound of second simplex row */
221 SCIP_Real* simplexcoefs1, /**< simplex coefficients of first row */
222 SCIP_Real* simplexcoefs2, /**< simplex coefficients of second row */
223 SCIP_Real* cutcoefs, /**< pointer to store cut coefficients (length: nscipvars) */
224 SCIP_ROW** row, /**< pointer to store disjunctive cut inequality */
225 SCIP_Bool* madeintegral /**< pointer to store whether cut has been scaled to integral values */
226 )
227 {
228 char cutname[SCIP_MAXSTRLEN];
229 SCIP_COL** rowcols;
230 SCIP_COL* col;
231 SCIP_Real* rowvals;
232 SCIP_Real lhsrow;
233 SCIP_Real rhsrow;
234 SCIP_Real cutlhs;
235 SCIP_Real sgn;
236 SCIP_Real lb;
237 SCIP_Real ub;
238 int nonbasicnumber = 0;
239 int rownnonz;
240 int ind;
241 int r;
242 int c;
243
244 assert( scip != NULL );
245 assert( row != NULL );
246 assert( rows != NULL );
247 assert( cols != NULL );
248 assert( simplexcoefs1 != NULL );
249 assert( simplexcoefs2 != NULL );
250 assert( cutcoefs != NULL );
251 assert( sepa != NULL );
252 assert( madeintegral != NULL );
253
254 *madeintegral = FALSE;
255
256 /* check signs */
257 if ( SCIPisFeasPositive(scip, cutlhs1) == SCIPisFeasPositive(scip, cutlhs2) )
258 sgn = 1.0;
259 else
260 sgn = -1.0;
261
262 /* check bounds */
263 if ( SCIPisInfinity(scip, REALABS(bound1)) || SCIPisInfinity(scip, REALABS(bound2)) )
264 strengthen = FALSE;
265
266 /* compute left hand side of row (a later update is possible, see below) */
267 cutlhs = sgn * cutlhs1 * cutlhs2;
268
269 /* add cut-coefficients of the non-basic non-slack variables */
270 for (c = 0; c < ncols; ++c)
271 {
272 col = cols[c];
273 assert( col != NULL );
274 ind = SCIPcolGetLPPos(col);
275 assert( ind >= 0 );
276
277 if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER )
278 {
279 lb = SCIPcolGetLb(col);
280
281 /* for integer variables we may obtain stronger coefficients */
282 if ( strengthen && SCIPcolIsIntegral(col) )
283 {
284 SCIP_Real mval;
285 SCIP_Real mvalfloor;
286 SCIP_Real mvalceil;
287
288 mval = (cutlhs2 * simplexcoefs1[nonbasicnumber] - cutlhs1 * simplexcoefs2[nonbasicnumber]) / (cutlhs2 * bound1 + cutlhs1 * bound2);
289 mvalfloor = SCIPfloor(scip, mval);
290 mvalceil = SCIPceil(scip, mval);
291
292 cutcoefs[ind] = MIN(sgn * cutlhs2 * (simplexcoefs1[nonbasicnumber] - mvalfloor * bound1), sgn * cutlhs1 * (simplexcoefs2[nonbasicnumber] + mvalceil * bound2));
293 assert( SCIPisFeasLE(scip, cutcoefs[ind], MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber])) );
294 }
295 else
296 cutcoefs[ind] = MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]);
297
298 cutlhs += cutcoefs[ind] * lb;
299 ++nonbasicnumber;
300 }
301 else if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER )
302 {
303 ub = SCIPcolGetUb(col);
304
305 /* for integer variables we may obtain stronger coefficients */
306 if ( strengthen && SCIPcolIsIntegral(col) )
307 {
308 SCIP_Real mval;
309 SCIP_Real mvalfloor;
310 SCIP_Real mvalceil;
311
312 mval = (cutlhs2 * simplexcoefs1[nonbasicnumber] - cutlhs1 * simplexcoefs2[nonbasicnumber]) / (cutlhs2 * bound1 + cutlhs1 * bound2);
313 mvalfloor = SCIPfloor(scip, -mval);
314 mvalceil = SCIPceil(scip, -mval);
315
316 cutcoefs[ind] = MAX(sgn * cutlhs2 * (simplexcoefs1[nonbasicnumber] + mvalfloor * bound1), sgn * cutlhs1 * (simplexcoefs2[nonbasicnumber] - mvalceil * bound2));
317 assert( SCIPisFeasLE(scip, -cutcoefs[ind], -MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber])) );
318 }
319 else
320 cutcoefs[ind] = MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]);
321
322 cutlhs += cutcoefs[ind] * ub;
323 ++nonbasicnumber;
324 }
325 else
326 {
327 assert( SCIPcolGetBasisStatus(col) != SCIP_BASESTAT_ZERO );
328 cutcoefs[ind] = 0.0;
329 }
330 }
331
332 /* add cut-coefficients of the non-basic slack variables */
333 for (r = 0; r < nrows; ++r)
334 {
335 rhsrow = SCIProwGetRhs(rows[r]) - SCIProwGetConstant(rows[r]);
336 lhsrow = SCIProwGetLhs(rows[r]) - SCIProwGetConstant(rows[r]);
337
338 assert( SCIProwGetBasisStatus(rows[r]) != SCIP_BASESTAT_ZERO );
339 assert( SCIPisFeasZero(scip, lhsrow - rhsrow) || SCIPisNegative(scip, lhsrow - rhsrow) );
340 assert( SCIProwIsInLP(rows[r]) );
341
342 if ( SCIProwGetBasisStatus(rows[r]) != SCIP_BASESTAT_BASIC )
343 {
344 SCIP_Real cutcoef;
345
346 if ( SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_UPPER )
347 {
348 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, rows[r]) - SCIProwGetRhs(rows[r])) );
349
350 cutcoef = MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]);
351 cutlhs -= cutcoef * rhsrow;
352 ++nonbasicnumber;
353 }
354 else /* SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_LOWER */
355 {
356 assert( SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_LOWER );
357 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, rows[r]) - SCIProwGetLhs(rows[r])) );
358
359 cutcoef = MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]);
360 cutlhs -= cutcoef * lhsrow;
361 ++nonbasicnumber;
362 }
363
364 rownnonz = SCIProwGetNNonz(rows[r]);
365 rowvals = SCIProwGetVals(rows[r]);
366 rowcols = SCIProwGetCols(rows[r]);
367
368 for (c = 0; c < rownnonz; ++c)
369 {
370 ind = SCIPcolGetLPPos(rowcols[c]);
371
372 /* if column is not in LP, then return without generating cut */
373 if ( ind < 0 )
374 {
375 *row = NULL;
376 return SCIP_OKAY;
377 }
378
379 cutcoefs[ind] -= cutcoef * rowvals[c];
380 }
381 }
382 }
383
384 /* create cut */
385 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d", SCIPsepaGetName(sepa), SCIPgetNLPs(scip), ndisjcuts);
386 if ( SCIPgetDepth(scip) == 0 )
387 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, cutname, cutlhs, SCIPinfinity(scip), FALSE, FALSE, TRUE) );
388 else
389 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, cutname, cutlhs, SCIPinfinity(scip), TRUE, FALSE, TRUE) );
390
391 SCIP_CALL( SCIPcacheRowExtensions(scip, *row) );
392 for (c = 0; c < ncols; ++c)
393 {
394 ind = SCIPcolGetLPPos(cols[c]);
395 assert( ind >= 0 );
396 if ( ! SCIPisFeasZero(scip, cutcoefs[ind]) )
397 {
398 SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPcolGetVar(cols[c]), cutcoefs[ind] ) );
399 }
400 }
401 SCIP_CALL( SCIPflushRowExtensions(scip, *row) );
402
403 /* try to scale the cut to integral values
404 * @todo find better but still stable disjunctive cut settings
405 */
406 if ( scale )
407 {
408 SCIP_Longint maxdnom;
409 SCIP_Real maxscale;
410
411 assert( SCIPgetDepth(scip) >= 0 );
412 if( SCIPgetDepth(scip) == 0 )
413 {
414 maxdnom = 100;
415 maxscale = 100.0;
416 }
417 else
418 {
419 maxdnom = 10;
420 maxscale = 10.0;
421 }
422
423 SCIP_CALL( SCIPmakeRowIntegral(scip, *row, -SCIPepsilon(scip), SCIPsumepsilon(scip), maxdnom, maxscale, TRUE, madeintegral) );
424 }
425
426 return SCIP_OKAY;
427 }
428
429
430 /*
431 * Callback methods
432 */
433
434 /** copy method for separator plugins (called when SCIP copies plugins) */
435 static
SCIP_DECL_SEPACOPY(sepaCopyDisjunctive)436 SCIP_DECL_SEPACOPY(sepaCopyDisjunctive)
437 {
438 assert( scip != NULL );
439 assert( sepa != NULL );
440 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
441
442 /* call inclusion method of constraint handler */
443 SCIP_CALL( SCIPincludeSepaDisjunctive(scip) );
444
445 return SCIP_OKAY;
446 }
447
448
449 /** destructor of separator to free user data (called when SCIP is exiting) */
450 static
SCIP_DECL_SEPAFREE(sepaFreeDisjunctive)451 SCIP_DECL_SEPAFREE(sepaFreeDisjunctive)/*lint --e{715}*/
452 {
453 SCIP_SEPADATA* sepadata;
454
455 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
456
457 /* free separator data */
458 sepadata = SCIPsepaGetData(sepa);
459 assert( sepadata != NULL );
460
461 SCIPfreeBlockMemory(scip, &sepadata);
462
463 SCIPsepaSetData(sepa, NULL);
464
465 return SCIP_OKAY;
466 }
467
468
469 /** solving process initialization method of separator */
470 static
SCIP_DECL_SEPAEXITSOL(sepaInitsolDisjunctive)471 SCIP_DECL_SEPAEXITSOL(sepaInitsolDisjunctive)
472 { /*lint --e{715}*/
473 SCIP_SEPADATA* sepadata;
474
475 sepadata = SCIPsepaGetData(sepa);
476 assert(sepadata != NULL);
477
478 sepadata->conshdlr = SCIPfindConshdlr(scip, "SOS1");
479
480 return SCIP_OKAY;
481 }
482
483
484 /** LP solution separation method for disjunctive cuts */
485 static
SCIP_DECL_SEPAEXECLP(sepaExeclpDisjunctive)486 SCIP_DECL_SEPAEXECLP(sepaExeclpDisjunctive)
487 {
488 SCIP_SEPADATA* sepadata;
489 SCIP_CONSHDLR* conshdlr;
490 SCIP_DIGRAPH* conflictgraph;
491 SCIP_ROW** rows;
492 SCIP_COL** cols;
493 SCIP_Real* cutcoefs = NULL;
494 SCIP_Real* simplexcoefs1 = NULL;
495 SCIP_Real* simplexcoefs2 = NULL;
496 SCIP_Real* coef = NULL;
497 SCIP_Real* binvrow = NULL;
498 SCIP_Real* rowsmaxval = NULL;
499 SCIP_Real* violationarray = NULL;
500 int* fixings1 = NULL;
501 int* fixings2 = NULL;
502 int* basisind = NULL;
503 int* basisrow = NULL;
504 int* varrank = NULL;
505 int* edgearray = NULL;
506 int nedges;
507 int ndisjcuts;
508 int nrelevantedges;
509 int nsos1vars;
510 int nconss;
511 int maxcuts;
512 int ncalls;
513 int depth;
514 int ncols;
515 int nrows;
516 int ind;
517 int j;
518 int i;
519
520 assert( sepa != NULL );
521 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
522 assert( scip != NULL );
523 assert( result != NULL );
524
525 *result = SCIP_DIDNOTRUN;
526
527 if( !allowlocal )
528 return SCIP_OKAY;
529
530 /* only generate disjunctive cuts if we are not close to terminating */
531 if ( SCIPisStopped(scip) )
532 return SCIP_OKAY;
533
534 /* only generate disjunctive cuts if an optimal LP solution is at hand */
535 if ( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
536 return SCIP_OKAY;
537
538 /* only generate disjunctive cuts if the LP solution is basic */
539 if ( ! SCIPisLPSolBasic(scip) )
540 return SCIP_OKAY;
541
542 /* get LP data */
543 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
544 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
545
546 /* return if LP has no columns or no rows */
547 if ( ncols == 0 || nrows == 0 )
548 return SCIP_OKAY;
549
550 assert( cols != NULL );
551 assert( rows != NULL );
552
553 /* get sepa data */
554 sepadata = SCIPsepaGetData(sepa);
555 assert( sepadata != NULL );
556
557 /* get constraint handler */
558 conshdlr = sepadata->conshdlr;
559 if ( conshdlr == NULL )
560 return SCIP_OKAY;
561
562 /* get number of constraints */
563 nconss = SCIPconshdlrGetNConss(conshdlr);
564 if ( nconss == 0 )
565 return SCIP_OKAY;
566
567 /* check for maxdepth < depth, maxinvcutsroot = 0 and maxinvcuts = 0 */
568 depth = SCIPgetDepth(scip);
569 if ( ( sepadata->maxdepth >= 0 && sepadata->maxdepth < depth )
570 || ( depth == 0 && sepadata->maxinvcutsroot == 0 )
571 || ( depth > 0 && sepadata->maxinvcuts == 0 ) )
572 return SCIP_OKAY;
573
574 /* only call the cut separator a given number of times at each node */
575 ncalls = SCIPsepaGetNCallsAtNode(sepa);
576 if ( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
577 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
578 return SCIP_OKAY;
579
580 /* get conflict graph and number of conflict graph edges (note that the digraph arcs were added in both directions) */
581 conflictgraph = SCIPgetConflictgraphSOS1(conshdlr);
582 if( conflictgraph == NULL )
583 return SCIP_OKAY;
584
585 nedges = (int)SCIPceil(scip, (SCIP_Real)SCIPdigraphGetNArcs(conflictgraph)/2);
586
587 /* if too many conflict graph edges, the separator can be slow: delay it until no other cuts have been found */
588 if ( sepadata->maxconfsdelay >= 0 && nedges >= sepadata->maxconfsdelay )
589 {
590 int ncutsfound;
591
592 ncutsfound = SCIPgetNCutsFound(scip);
593 if ( ncutsfound > sepadata->lastncutsfound || ! SCIPsepaWasLPDelayed(sepa) )
594 {
595 sepadata->lastncutsfound = ncutsfound;
596 *result = SCIP_DELAYED;
597 return SCIP_OKAY;
598 }
599 }
600
601 /* check basis status */
602 for (j = 0; j < ncols; ++j)
603 {
604 if ( SCIPcolGetBasisStatus(cols[j]) == SCIP_BASESTAT_ZERO )
605 return SCIP_OKAY;
606 }
607
608 /* check accuracy of rows */
609 for (j = 0; j < nrows; ++j)
610 {
611 SCIP_ROW* row;
612
613 row = rows[j];
614
615 if ( ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER && ! SCIPisEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetRhs(row)) )
616 || ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER && ! SCIPisEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row)) ) )
617 return SCIP_OKAY;
618 }
619
620 /* get number of SOS1 variables */
621 nsos1vars = SCIPgetNSOS1Vars(conshdlr);
622
623 /* allocate buffer arrays */
624 SCIP_CALL( SCIPallocBufferArray(scip, &edgearray, nedges) );
625 SCIP_CALL( SCIPallocBufferArray(scip, &fixings1, nedges) );
626 SCIP_CALL( SCIPallocBufferArray(scip, &fixings2, nedges) );
627 SCIP_CALL( SCIPallocBufferArray(scip, &violationarray, nedges) );
628
629 /* get all violated conflicts {i, j} in the conflict graph and sort them based on the degree of a violation value */
630 nrelevantedges = 0;
631 for (j = 0; j < nsos1vars; ++j)
632 {
633 SCIP_VAR* var;
634
635 var = SCIPnodeGetVarSOS1(conflictgraph, j);
636
637 if ( SCIPvarIsActive(var) && ! SCIPisFeasZero(scip, SCIPcolGetPrimsol(SCIPvarGetCol(var))) && SCIPcolGetBasisStatus(SCIPvarGetCol(var)) == SCIP_BASESTAT_BASIC )
638 {
639 int* succ;
640 int nsucc;
641
642 /* get successors and number of successors */
643 nsucc = SCIPdigraphGetNSuccessors(conflictgraph, j);
644 succ = SCIPdigraphGetSuccessors(conflictgraph, j);
645
646 for (i = 0; i < nsucc; ++i)
647 {
648 SCIP_VAR* varsucc;
649 int succind;
650
651 succind = succ[i];
652 varsucc = SCIPnodeGetVarSOS1(conflictgraph, succind);
653 if ( SCIPvarIsActive(varsucc) && succind < j && ! SCIPisFeasZero(scip, SCIPgetSolVal(scip, NULL, varsucc) ) &&
654 SCIPcolGetBasisStatus(SCIPvarGetCol(varsucc)) == SCIP_BASESTAT_BASIC )
655 {
656 fixings1[nrelevantedges] = j;
657 fixings2[nrelevantedges] = succind;
658 edgearray[nrelevantedges] = nrelevantedges;
659 violationarray[nrelevantedges++] = SCIPgetSolVal(scip, NULL, var) * SCIPgetSolVal(scip, NULL, varsucc);
660 }
661 }
662 }
663 }
664
665 /* sort violation score values */
666 if ( nrelevantedges > 0)
667 SCIPsortDownRealInt(violationarray, edgearray, nrelevantedges);
668 else
669 {
670 SCIPfreeBufferArrayNull(scip, &violationarray);
671 SCIPfreeBufferArrayNull(scip, &fixings2);
672 SCIPfreeBufferArrayNull(scip, &fixings1);
673 SCIPfreeBufferArrayNull(scip, &edgearray);
674
675 return SCIP_OKAY;
676 }
677 SCIPfreeBufferArrayNull(scip, &violationarray);
678
679 /* compute maximal number of cuts */
680 if ( SCIPgetDepth(scip) == 0 )
681 maxcuts = MIN(sepadata->maxinvcutsroot, nrelevantedges);
682 else
683 maxcuts = MIN(sepadata->maxinvcuts, nrelevantedges);
684 assert( maxcuts > 0 );
685
686 /* allocate buffer arrays */
687 SCIP_CALL( SCIPallocBufferArray(scip, &varrank, ncols) );
688 SCIP_CALL( SCIPallocBufferArray(scip, &rowsmaxval, nrows) );
689 SCIP_CALL( SCIPallocBufferArray(scip, &basisrow, ncols) );
690 SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
691 SCIP_CALL( SCIPallocBufferArray(scip, &coef, ncols) );
692 SCIP_CALL( SCIPallocBufferArray(scip, &simplexcoefs1, ncols) );
693 SCIP_CALL( SCIPallocBufferArray(scip, &simplexcoefs2, ncols) );
694 SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) );
695 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
696
697 /* get basis indices */
698 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
699
700 /* create vector "basisrow" with basisrow[column of non-slack basis variable] = corresponding row of B^-1;
701 * compute maximum absolute value of nonbasic row coefficients */
702 for (j = 0; j < nrows; ++j)
703 {
704 SCIP_COL** rowcols;
705 SCIP_Real* rowvals;
706 SCIP_ROW* row;
707 SCIP_Real val;
708 SCIP_Real max = 0.0;
709 int nnonz;
710
711 /* fill basisrow vector */
712 ind = basisind[j];
713 if ( ind >= 0 )
714 basisrow[ind] = j;
715
716 /* compute maximum absolute value of nonbasic row coefficients */
717 row = rows[j];
718 assert( row != NULL );
719 rowvals = SCIProwGetVals(row);
720 nnonz = SCIProwGetNNonz(row);
721 rowcols = SCIProwGetCols(row);
722
723 for (i = 0; i < nnonz; ++i)
724 {
725 if ( SCIPcolGetBasisStatus(rowcols[i]) == SCIP_BASESTAT_LOWER || SCIPcolGetBasisStatus(rowcols[i]) == SCIP_BASESTAT_UPPER )
726 {
727 val = REALABS(rowvals[i]);
728 if ( SCIPisFeasGT(scip, val, max) )
729 max = REALABS(val);
730 }
731 }
732
733 /* handle slack variable coefficient and save maximum value */
734 rowsmaxval[j] = MAX(max, 1.0);
735 }
736
737 /* initialize variable ranks with -1 */
738 for (j = 0; j < ncols; ++j)
739 varrank[j] = -1;
740
741 /* free buffer array */
742 SCIPfreeBufferArrayNull(scip, &basisind);
743
744 /* for the most promising disjunctions: try to generate disjunctive cuts */
745 ndisjcuts = 0;
746 for (i = 0; i < maxcuts; ++i)
747 {
748 SCIP_Bool madeintegral;
749 SCIP_Real cutlhs1;
750 SCIP_Real cutlhs2;
751 SCIP_Real bound1;
752 SCIP_Real bound2;
753 SCIP_ROW* row = NULL;
754 SCIP_VAR* var;
755 SCIP_COL* col;
756
757 int nonbasicnumber;
758 int cutrank = 0;
759 int edgenumber;
760 int rownnonz;
761
762 edgenumber = edgearray[i];
763
764 /* determine first simplex row */
765 var = SCIPnodeGetVarSOS1(conflictgraph, fixings1[edgenumber]);
766 col = SCIPvarGetCol(var);
767 ind = SCIPcolGetLPPos(col);
768 assert( ind >= 0 );
769 assert( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC );
770
771 /* get the 'ind'th row of B^-1 and B^-1 \cdot A */
772 SCIP_CALL( SCIPgetLPBInvRow(scip, basisrow[ind], binvrow, NULL, NULL) );
773 SCIP_CALL( SCIPgetLPBInvARow(scip, basisrow[ind], binvrow, coef, NULL, NULL) );
774
775 /* get the simplex-coefficients of the non-basic variables */
776 SCIP_CALL( getSimplexCoefficients(scip, rows, nrows, cols, ncols, coef, binvrow, simplexcoefs1, &nonbasicnumber) );
777
778 /* get rank of variable if not known already */
779 if ( varrank[ind] < 0 )
780 varrank[ind] = getVarRank(scip, binvrow, rowsmaxval, sepadata->maxweightrange, rows, nrows);
781 cutrank = MAX(cutrank, varrank[ind]);
782
783 /* get right hand side and bound of simplex talbeau row */
784 cutlhs1 = SCIPcolGetPrimsol(col);
785 if ( SCIPisFeasPositive(scip, cutlhs1) )
786 bound1 = SCIPcolGetUb(col);
787 else
788 bound1 = SCIPcolGetLb(col);
789
790 /* determine second simplex row */
791 var = SCIPnodeGetVarSOS1(conflictgraph, fixings2[edgenumber]);
792 col = SCIPvarGetCol(var);
793 ind = SCIPcolGetLPPos(col);
794 assert( ind >= 0 );
795 assert( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC );
796
797 /* get the 'ind'th row of B^-1 and B^-1 \cdot A */
798 SCIP_CALL( SCIPgetLPBInvRow(scip, basisrow[ind], binvrow, NULL, NULL) );
799 SCIP_CALL( SCIPgetLPBInvARow(scip, basisrow[ind], binvrow, coef, NULL, NULL) );
800
801 /* get the simplex-coefficients of the non-basic variables */
802 SCIP_CALL( getSimplexCoefficients(scip, rows, nrows, cols, ncols, coef, binvrow, simplexcoefs2, &nonbasicnumber) );
803
804 /* get rank of variable if not known already */
805 if ( varrank[ind] < 0 )
806 varrank[ind] = getVarRank(scip, binvrow, rowsmaxval, sepadata->maxweightrange, rows, nrows);
807 cutrank = MAX(cutrank, varrank[ind]);
808
809 /* get right hand side and bound of simplex talbeau row */
810 cutlhs2 = SCIPcolGetPrimsol(col);
811 if ( SCIPisFeasPositive(scip, cutlhs2) )
812 bound2 = SCIPcolGetUb(col);
813 else
814 bound2 = SCIPcolGetLb(col);
815
816 /* add coefficients to cut */
817 SCIP_CALL( generateDisjCutSOS1(scip, sepa, rows, nrows, cols, ncols, ndisjcuts, MAKECONTINTEGRAL, sepadata->strengthen, cutlhs1, cutlhs2, bound1, bound2, simplexcoefs1, simplexcoefs2, cutcoefs, &row, &madeintegral) );
818 if ( row == NULL )
819 continue;
820
821 /* raise cutrank for present cut */
822 ++cutrank;
823
824 /* check if there are numerical evidences */
825 if ( ( madeintegral && ( sepadata->maxrankintegral == -1 || cutrank <= sepadata->maxrankintegral ) )
826 || ( ! madeintegral && ( sepadata->maxrank == -1 || cutrank <= sepadata->maxrank ) ) )
827 {
828 /* possibly add cut to LP if it is useful; in case the lhs of the cut is minus infinity (due to scaling) the cut is useless */
829 rownnonz = SCIProwGetNNonz(row);
830 if ( rownnonz > 0 && ! SCIPisInfinity(scip, -SCIProwGetLhs(row)) && ! SCIProwIsInLP(row) && SCIPisCutEfficacious(scip, NULL, row) )
831 {
832 SCIP_Bool infeasible;
833
834 /* set cut rank */
835 SCIProwChgRank(row, cutrank);
836
837 /* add cut */
838 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
839 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
840 if ( infeasible )
841 {
842 *result = SCIP_CUTOFF;
843 break;
844 }
845 ++ndisjcuts;
846 }
847 }
848
849 /* release row */
850 SCIP_CALL( SCIPreleaseRow(scip, &row) );
851 }
852
853 /* save total number of cuts found so far */
854 sepadata->lastncutsfound = SCIPgetNCutsFound(scip);
855
856 /* evaluate the result of the separation */
857 if ( *result != SCIP_CUTOFF )
858 {
859 if ( ndisjcuts > 0 )
860 *result = SCIP_SEPARATED;
861 else
862 *result = SCIP_DIDNOTFIND;
863 }
864
865 SCIPdebugMsg(scip, "Number of found disjunctive cuts: %d.\n", ndisjcuts);
866
867 /* free buffer arrays */
868 SCIPfreeBufferArrayNull(scip, &cutcoefs);
869 SCIPfreeBufferArrayNull(scip, &simplexcoefs2);
870 SCIPfreeBufferArrayNull(scip, &simplexcoefs1);
871 SCIPfreeBufferArrayNull(scip, &coef);
872 SCIPfreeBufferArrayNull(scip, &binvrow);
873 SCIPfreeBufferArrayNull(scip, &basisrow);
874 SCIPfreeBufferArrayNull(scip, &rowsmaxval);
875 SCIPfreeBufferArrayNull(scip, &varrank);
876 SCIPfreeBufferArrayNull(scip, &fixings2);
877 SCIPfreeBufferArrayNull(scip, &fixings1);
878 SCIPfreeBufferArrayNull(scip, &edgearray);
879
880 return SCIP_OKAY;
881 }
882
883
884 /** creates the disjunctive cut separator and includes it in SCIP */
SCIPincludeSepaDisjunctive(SCIP * scip)885 SCIP_RETCODE SCIPincludeSepaDisjunctive(
886 SCIP* scip /**< SCIP data structure */
887 )
888 {
889 SCIP_SEPADATA* sepadata = NULL;
890 SCIP_SEPA* sepa = NULL;
891
892 /* create separator data */
893 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
894 sepadata->conshdlr = NULL;
895 sepadata->lastncutsfound = 0;
896
897 /* include separator */
898 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
899 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpDisjunctive, NULL, sepadata) );
900
901 assert( sepa != NULL );
902
903 /* set non fundamental callbacks via setter functions */
904 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyDisjunctive) );
905 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeDisjunctive) );
906 SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolDisjunctive) );
907
908 /* add separator parameters */
909 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/strengthen",
910 "strengthen cut if integer variables are present.",
911 &sepadata->strengthen, TRUE, DEFAULT_STRENGTHEN, NULL, NULL) );
912
913 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxdepth",
914 "node depth of separating bipartite disjunctive cuts (-1: no limit)",
915 &sepadata->maxdepth, TRUE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) );
916
917 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds",
918 "maximal number of separation rounds per iteration in a branching node (-1: no limit)",
919 &sepadata->maxrounds, TRUE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
920
921 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot",
922 "maximal number of separation rounds in the root node (-1: no limit)",
923 &sepadata->maxroundsroot, TRUE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
924
925 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxinvcuts",
926 "maximal number of cuts investigated per iteration in a branching node",
927 &sepadata->maxinvcuts, TRUE, DEFAULT_MAXINVCUTS, 0, INT_MAX, NULL, NULL) );
928
929 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxinvcutsroot",
930 "maximal number of cuts investigated per iteration in the root node",
931 &sepadata->maxinvcutsroot, TRUE, DEFAULT_MAXINVCUTSROOT, 0, INT_MAX, NULL, NULL) );
932
933 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxconfsdelay",
934 "delay separation if number of conflict graph edges is larger than predefined value (-1: no limit)",
935 &sepadata->maxconfsdelay, TRUE, DEFAULT_MAXCONFSDELAY, -1, INT_MAX, NULL, NULL) );
936
937 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrank",
938 "maximal rank of a disj. cut that could not be scaled to integral coefficients (-1: unlimited)",
939 &sepadata->maxrank, FALSE, DEFAULT_MAXRANK, -1, INT_MAX, NULL, NULL) );
940
941 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrankintegral",
942 "maximal rank of a disj. cut that could be scaled to integral coefficients (-1: unlimited)",
943 &sepadata->maxrankintegral, FALSE, DEFAULT_MAXRANKINTEGRAL, -1, INT_MAX, NULL, NULL) );
944
945 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/maxweightrange",
946 "maximal valid range max(|weights|)/min(|weights|) of row weights",
947 &sepadata->maxweightrange, TRUE, DEFAULT_MAXWEIGHTRANGE, 1.0, SCIP_REAL_MAX, NULL, NULL) );
948
949 return SCIP_OKAY;
950 }
951