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_strongcg.c
17 * @ingroup DEFPLUGINS_SEPA
18 * @brief Strong CG Cuts (Letchford & Lodi)
19 * @author Kati Wolter
20 * @author Tobias Achterberg
21 */
22
23 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
24
25 #include "blockmemshell/memory.h"
26 #include "scip/cuts.h"
27 #include "scip/pub_lp.h"
28 #include "scip/pub_message.h"
29 #include "scip/pub_misc.h"
30 #include "scip/pub_misc_sort.h"
31 #include "scip/pub_sepa.h"
32 #include "scip/pub_var.h"
33 #include "scip/scip_branch.h"
34 #include "scip/scip_cut.h"
35 #include "scip/scip_general.h"
36 #include "scip/scip_lp.h"
37 #include "scip/scip_mem.h"
38 #include "scip/scip_message.h"
39 #include "scip/scip_numerics.h"
40 #include "scip/scip_param.h"
41 #include "scip/scip_prob.h"
42 #include "scip/scip_randnumgen.h"
43 #include "scip/scip_sepa.h"
44 #include "scip/scip_solvingstats.h"
45 #include "scip/scip_tree.h"
46 #include "scip/sepa_strongcg.h"
47 #include <string.h>
48
49
50 #define SEPA_NAME "strongcg"
51 #define SEPA_DESC "Strong CG cuts separator (Letchford and Lodi)"
52 #define SEPA_PRIORITY -2000
53 #define SEPA_FREQ 10
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_MAXROUNDS 5 /**< maximal number of strong CG separation rounds per node (-1: unlimited) */
59 #define DEFAULT_MAXROUNDSROOT 20 /**< maximal number of strong CG separation rounds in the root node (-1: unlimited) */
60 #define DEFAULT_MAXSEPACUTS 20 /**< maximal number of strong CG cuts separated per separation round */
61 #define DEFAULT_MAXSEPACUTSROOT 500 /**< maximal number of strong CG cuts separated per separation round in root node */
62 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
63 #define DEFAULT_RANDSEED 54 /**< initial random seed */
64
65 #define SEPARATEROWS /* separate rows with integral slack */
66
67 #define BOUNDSWITCH 0.9999
68 #define POSTPROCESS TRUE
69 #define USEVBDS TRUE
70 #define MINFRAC 0.05
71 #define MAXFRAC 0.95
72
73 #define MAXAGGRLEN(nvars) (0.1*(nvars)+1000) /**< maximal length of base inequality */
74
75
76 /** separator data */
77 struct SCIP_SepaData
78 {
79 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
80 int maxrounds; /**< maximal number of strong CG separation rounds per node (-1: unlimited) */
81 int maxroundsroot; /**< maximal number of strong CG separation rounds in the root node (-1: unlimited) */
82 int maxsepacuts; /**< maximal number of strong CG cuts separated per separation round */
83 int maxsepacutsroot; /**< maximal number of strong CG cuts separated per separation round in root node */
84 int lastncutsfound; /**< total number of cuts found after last call of separator */
85 SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */
86 };
87
88 /*
89 * Callback methods
90 */
91
92 /** copy method for separator plugins (called when SCIP copies plugins) */
93 static
SCIP_DECL_SEPACOPY(sepaCopyStrongcg)94 SCIP_DECL_SEPACOPY(sepaCopyStrongcg)
95 { /*lint --e{715}*/
96 assert(scip != NULL);
97 assert(sepa != NULL);
98 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
99
100 /* call inclusion method of constraint handler */
101 SCIP_CALL( SCIPincludeSepaStrongcg(scip) );
102
103 return SCIP_OKAY;
104 }
105
106 /** destructor of separator to free user data (called when SCIP is exiting) */
107 static
SCIP_DECL_SEPAFREE(sepaFreeStrongcg)108 SCIP_DECL_SEPAFREE(sepaFreeStrongcg)
109 { /*lint --e{715}*/
110 SCIP_SEPADATA* sepadata;
111
112 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
113
114 /* free separator data */
115 sepadata = SCIPsepaGetData(sepa);
116 assert(sepadata != NULL);
117
118 SCIPfreeBlockMemory(scip, &sepadata);
119
120 SCIPsepaSetData(sepa, NULL);
121
122 return SCIP_OKAY;
123 }
124
125 /** initialization method of separator (called after problem was transformed) */
126 static
SCIP_DECL_SEPAINIT(sepaInitStrongcg)127 SCIP_DECL_SEPAINIT(sepaInitStrongcg)
128 {
129 SCIP_SEPADATA* sepadata;
130
131 sepadata = SCIPsepaGetData(sepa);
132 assert(sepadata != NULL);
133
134 /* create and initialize random number generator */
135 SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) );
136
137 return SCIP_OKAY;
138 }
139
140 /** deinitialization method of separator (called before transformed problem is freed) */
141 static
SCIP_DECL_SEPAEXIT(sepaExitStrongcg)142 SCIP_DECL_SEPAEXIT(sepaExitStrongcg)
143 { /*lint --e{715}*/
144 SCIP_SEPADATA* sepadata;
145
146 sepadata = SCIPsepaGetData(sepa);
147 assert(sepadata != NULL);
148
149 SCIPfreeRandom(scip, &sepadata->randnumgen);
150
151 return SCIP_OKAY;
152 }
153
154 /** LP solution separation method of separator */
155 static
SCIP_DECL_SEPAEXECLP(sepaExeclpStrongcg)156 SCIP_DECL_SEPAEXECLP(sepaExeclpStrongcg)
157 { /*lint --e{715}*/
158 SCIP_SEPADATA* sepadata;
159 SCIP_VAR** vars;
160 SCIP_COL** cols;
161 SCIP_ROW** rows;
162 SCIP_AGGRROW* aggrrow;
163 SCIP_Real* varsolvals;
164 SCIP_Real* binvrow;
165 SCIP_Real* cutcoefs;
166 SCIP_Real* basisfrac;
167 SCIP_Real cutrhs;
168 int* basisind;
169 int* basisperm;
170 int* inds;
171 int* cutinds;
172 int cutnnz;
173 int ninds;
174 int nvars;
175 int ncols;
176 int nrows;
177 int ncalls;
178 int depth;
179 int maxsepacuts;
180 int ncuts;
181 int c;
182 int i;
183 int cutrank;
184 SCIP_Bool success;
185 SCIP_Bool cutislocal;
186
187 assert(sepa != NULL);
188 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
189 assert(scip != NULL);
190 assert(result != NULL);
191
192 *result = SCIP_DIDNOTRUN;
193
194 sepadata = SCIPsepaGetData(sepa);
195 assert(sepadata != NULL);
196
197 depth = SCIPgetDepth(scip);
198 ncalls = SCIPsepaGetNCallsAtNode(sepa);
199
200 /* only call separator, if we are not close to terminating */
201 if( SCIPisStopped(scip) || !allowlocal )
202 return SCIP_OKAY;
203
204 /* only call the strong CG cut separator a given number of times at each node */
205 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
206 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
207 return SCIP_OKAY;
208
209 /* only call separator, if an optimal LP solution is at hand */
210 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
211 return SCIP_OKAY;
212
213 /* only call separator, if the LP solution is basic */
214 if( !SCIPisLPSolBasic(scip) )
215 return SCIP_OKAY;
216
217 /* only call separator, if there are fractional variables */
218 if( SCIPgetNLPBranchCands(scip) == 0 )
219 return SCIP_OKAY;
220
221 /* get variables data */
222 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
223
224 /* get LP data */
225 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
226 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
227 if( ncols == 0 || nrows == 0 )
228 return SCIP_OKAY;
229
230 #if 0 /* if too many columns, separator is usually very slow: delay it until no other cuts have been found */
231 if( ncols >= 50*nrows )
232 return SCIP_OKAY;
233 if( ncols >= 5*nrows )
234 {
235 int ncutsfound;
236
237 ncutsfound = SCIPgetNCutsFound(scip);
238 if( ncutsfound > sepadata->lastncutsfound || !SCIPsepaWasLPDelayed(sepa) )
239 {
240 sepadata->lastncutsfound = ncutsfound;
241 *result = SCIP_DELAYED;
242 return SCIP_OKAY;
243 }
244 }
245 #endif
246
247 *result = SCIP_DIDNOTFIND;
248
249 /* allocate temporary memory */
250 SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
251 SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, nvars) );
252 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) );
253 SCIP_CALL( SCIPallocBufferArray(scip, &basisperm, nrows) );
254 SCIP_CALL( SCIPallocBufferArray(scip, &basisfrac, nrows) );
255 SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) );
256 SCIP_CALL( SCIPallocBufferArray(scip, &inds, nrows) );
257 SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );
258
259 varsolvals = NULL; /* allocate this later, if needed */
260
261 /* get basis indices */
262 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) );
263
264 for( i = 0; i < nrows; ++i )
265 {
266 SCIP_Real frac = 0.0;
267
268 c = basisind[i];
269
270 basisperm[i] = i;
271
272 if( c >= 0 )
273 {
274 SCIP_VAR* var;
275
276 assert(c < ncols);
277 var = SCIPcolGetVar(cols[c]);
278 if( SCIPvarGetType(var) != SCIP_VARTYPE_CONTINUOUS )
279 {
280 frac = SCIPfeasFrac(scip, SCIPcolGetPrimsol(cols[c]));
281 frac = MIN(frac, 1.0 - frac);
282 }
283 }
284 #ifdef SEPARATEROWS
285 else
286 {
287 SCIP_ROW* row;
288
289 assert(0 <= -c-1 && -c-1 < nrows);
290 row = rows[-c-1];
291 if( SCIProwIsIntegral(row) && !SCIProwIsModifiable(row) )
292 {
293 frac = SCIPfeasFrac(scip, SCIPgetRowActivity(scip, row));
294 frac = MIN(frac, 1.0 - frac);
295 }
296 }
297 #endif
298
299 if( frac >= MINFRAC )
300 {
301 /* slightly change fractionality to have random order for equal fractions */
302 basisfrac[i] = frac + SCIPrandomGetReal(sepadata->randnumgen, -1e-6, 1e-6);
303 }
304 else
305 {
306 basisfrac[i] = 0.0;
307 }
308 }
309
310 /* sort basis indices by fractionality */
311 SCIPsortDownRealInt(basisfrac, basisperm, nrows);
312
313 /* get the maximal number of cuts allowed in a separation round */
314 if( depth == 0 )
315 maxsepacuts = sepadata->maxsepacutsroot;
316 else
317 maxsepacuts = sepadata->maxsepacuts;
318
319 SCIPdebugMsg(scip, "searching strong CG cuts: %d cols, %d rows, maxcuts=%d\n",
320 ncols, nrows, maxsepacuts);
321
322 /* for all basic columns belonging to integer variables, try to generate a strong CG cut */
323 ncuts = 0;
324 for( i = 0; i < nrows && ncuts < maxsepacuts && !SCIPisStopped(scip) && *result != SCIP_CUTOFF; ++i )
325 {
326 int j;
327 SCIP_Real cutefficacy;
328
329 if( basisfrac[i] == 0.0 )
330 break;
331
332 j = basisperm[i];
333 c = basisind[j];
334
335 /* get the row of B^-1 for this basic integer variable with fractional solution value */
336 SCIP_CALL( SCIPgetLPBInvRow(scip, j, binvrow, inds, &ninds) );
337
338 #ifdef SCIP_DEBUG
339 /* initialize variables, that might not have been initialized in SCIPcalcMIR if success == FALSE */
340 cutefficacy = 0.0;
341 cutrhs = SCIPinfinity(scip);
342 #endif
343
344 /* create the aggregation row using the B^-1 row as weights */
345 SCIP_CALL( SCIPaggrRowSumRows(scip, aggrrow, binvrow, inds, ninds,
346 FALSE, allowlocal, 1, (int) MAXAGGRLEN(nvars), &success) );
347
348 if( !success )
349 continue;
350
351 /* create a strong CG cut out of the aggregation row */
352 SCIP_CALL( SCIPcalcStrongCG(scip, NULL, POSTPROCESS, BOUNDSWITCH, USEVBDS, allowlocal, MINFRAC, MAXFRAC,
353 1.0, aggrrow, cutcoefs, &cutrhs, cutinds, &cutnnz, &cutefficacy, &cutrank, &cutislocal, &success) );
354
355 assert(allowlocal || !cutislocal);
356 SCIPdebugMsg(scip, " -> success=%u: rhs: %g, efficacy: %g\n", success, cutrhs, cutefficacy);
357
358 if( !success )
359 continue;
360
361 /* if successful, convert dense cut into sparse row, and add the row as a cut */
362 if( SCIPisEfficacious(scip, cutefficacy) )
363 {
364 SCIP_ROW* cut;
365 char cutname[SCIP_MAXSTRLEN];
366 int v;
367
368 SCIPdebugMsg(scip, " -> strong CG cut for <%s>: act=%f, rhs=%f, norm=%f, eff=%f, rank=%d\n",
369 c >= 0 ? SCIPvarGetName(SCIPcolGetVar(cols[c])) : SCIProwGetName(rows[-c-1]),
370 cutefficacy * SCIPgetVectorEfficacyNorm(scip, cutcoefs, cutnnz) + cutrhs, cutrhs, SCIPgetVectorEfficacyNorm(scip, cutcoefs, cutnnz), cutefficacy, cutrank);
371
372 /* create the cut */
373 if( c >= 0 )
374 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "scg%" SCIP_LONGINT_FORMAT "_x%d", SCIPgetNLPs(scip), c);
375 else
376 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "scg%" SCIP_LONGINT_FORMAT "_s%d", SCIPgetNLPs(scip), -c-1);
377
378 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, sepadata->dynamiccuts) );
379
380 /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, cut, NULL)) );*/
381 SCIProwChgRank(cut, cutrank);
382
383 /* cache the row extension and only flush them if the cut gets added */
384 SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
385
386 /* collect all non-zero coefficients */
387 for( v = 0; v < cutnnz; ++v )
388 {
389 SCIP_CALL( SCIPaddVarToRow(scip, cut, vars[cutinds[v]], cutcoefs[v]) );
390 }
391
392 assert(success);
393
394 if( !SCIPisCutEfficacious(scip, NULL, cut) )
395 {
396 SCIPdebugMsg(scip, " -> strong CG cut <%s> no longer efficacious: act=%f, rhs=%f, norm=%f, eff=%f\n",
397 cutname, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
398 SCIPgetCutEfficacy(scip, NULL, cut));
399 /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, cut, NULL)) );*/
400 success = FALSE;
401 }
402 else
403 {
404 SCIP_Bool infeasible = FALSE;
405
406 /* flush all changes before adding the cut */
407 SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
408
409 SCIPdebugMsg(scip, " -> found strong CG cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
410 cutname, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
411 SCIPgetCutEfficacy(scip, NULL, cut),
412 SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
413 SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
414 /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, cut, NULL)) );*/
415
416 if( SCIPisCutNew(scip, cut) )
417 {
418 if( !cutislocal )
419 {
420 SCIP_CALL( SCIPaddPoolCut(scip, cut) );
421 }
422 else
423 {
424 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, &infeasible) );
425 }
426
427 ncuts++;
428
429 if( infeasible )
430 {
431 *result = SCIP_CUTOFF;
432 }
433 else
434 {
435 *result = SCIP_SEPARATED;
436 }
437 }
438 }
439
440 /* release the row */
441 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
442 }
443 }
444
445 /* free temporary memory */
446 SCIPaggrRowFree(scip, &aggrrow);
447 SCIPfreeBufferArrayNull(scip, &varsolvals);
448 SCIPfreeBufferArray(scip, &inds);
449 SCIPfreeBufferArray(scip, &binvrow);
450 SCIPfreeBufferArray(scip, &basisfrac);
451 SCIPfreeBufferArray(scip, &basisperm);
452 SCIPfreeBufferArray(scip, &basisind);
453 SCIPfreeBufferArray(scip, &cutinds);
454 SCIPfreeBufferArray(scip, &cutcoefs);
455
456 SCIPdebugMsg(scip, "end searching strong CG cuts: found %d cuts\n", ncuts);
457
458 sepadata->lastncutsfound = SCIPgetNCutsFound(scip);
459
460 return SCIP_OKAY;
461 }
462
463
464 /*
465 * separator specific interface methods
466 */
467
468 /** creates the Strong CG cut separator and includes it in SCIP */
SCIPincludeSepaStrongcg(SCIP * scip)469 SCIP_RETCODE SCIPincludeSepaStrongcg(
470 SCIP* scip /**< SCIP data structure */
471 )
472 {
473 SCIP_SEPADATA* sepadata;
474 SCIP_SEPA* sepa;
475
476 /* create separator data */
477 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) );
478 sepadata->lastncutsfound = 0;
479
480 /* include separator */
481 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
482 SEPA_USESSUBSCIP, SEPA_DELAY,
483 sepaExeclpStrongcg, NULL,
484 sepadata) );
485
486 assert(sepa != NULL);
487
488 /* set non-NULL pointers to callback methods */
489 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyStrongcg) );
490 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeStrongcg) );
491 SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitStrongcg) );
492 SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitStrongcg) );
493
494 /* add separator parameters */
495 SCIP_CALL( SCIPaddIntParam(scip,
496 "separating/strongcg/maxrounds",
497 "maximal number of strong CG separation rounds per node (-1: unlimited)",
498 &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
499 SCIP_CALL( SCIPaddIntParam(scip,
500 "separating/strongcg/maxroundsroot",
501 "maximal number of strong CG separation rounds in the root node (-1: unlimited)",
502 &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
503 SCIP_CALL( SCIPaddIntParam(scip,
504 "separating/strongcg/maxsepacuts",
505 "maximal number of strong CG cuts separated per separation round",
506 &sepadata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
507 SCIP_CALL( SCIPaddIntParam(scip,
508 "separating/strongcg/maxsepacutsroot",
509 "maximal number of strong CG cuts separated per separation round in the root node",
510 &sepadata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
511 SCIP_CALL( SCIPaddBoolParam(scip,
512 "separating/strongcg/dynamiccuts",
513 "should generated cuts be removed from the LP if they are no longer tight?",
514 &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
515
516 return SCIP_OKAY;
517 }
518
519