1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file sepa_rlt.c
17 * @ingroup DEFPLUGINS_SEPA
18 * @brief separator for cuts generated by Reformulation-Linearization-Technique (RLT)
19 * @author Fabian Wegscheider
20 * @author Ksenia Bestuzheva
21 *
22 * @todo implement the possibility to add extra auxiliary variables for RLT (like in DOI 10.1080/10556788.2014.916287)
23 * @todo add RLT cuts for the product of equality constraints
24 * @todo implement dynamic addition of RLT cuts during branching (see DOI 10.1007/s10898-012-9874-7)
25 * @todo use SCIPvarIsBinary instead of SCIPvarGetType() == SCIP_VARTYPE_BINARY ?
26 * @todo parameter maxusedvars seems arbitrary (too large for small problems; too small for large problems); something more adaptive we can do? (e.g., all variables with priority >= x% of highest prio)
27 */
28
29 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
30
31 #include <assert.h>
32 #include <string.h>
33
34 #include "scip/sepa_rlt.h"
35 #include "scip/cons_nonlinear.h"
36 #include "scip/pub_lp.h"
37 #include "scip/expr_pow.h"
38 #include "scip/nlhdlr_bilinear.h"
39 #include "scip/cutsel_hybrid.h"
40
41
42 #define SEPA_NAME "rlt"
43 #define SEPA_DESC "reformulation-linearization-technique separator"
44 #define SEPA_PRIORITY 10 /**< priority for separation */
45 #define SEPA_FREQ 0 /**< frequency for separating cuts; zero means to separate only in the root node */
46 #define SEPA_MAXBOUNDDIST 1.0 /**< maximal relative distance from the current node's dual bound to primal bound
47 * compared to best node's dual bound for applying separation.*/
48 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
49 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
50
51 #define DEFAULT_MAXUNKNOWNTERMS 0 /**< maximum number of unknown bilinear terms a row can have to be used */
52 #define DEFAULT_MAXUSEDVARS 100 /**< maximum number of variables that will be used to compute rlt cuts */
53 #define DEFAULT_MAXNCUTS -1 /**< maximum number of cuts that will be added per round */
54 #define DEFAULT_MAXROUNDS 1 /**< maximum number of separation rounds per node (-1: unlimited) */
55 #define DEFAULT_MAXROUNDSROOT 10 /**< maximum number of separation rounds in the root node (-1: unlimited) */
56 #define DEFAULT_ONLYEQROWS FALSE /**< whether only equality rows should be used for rlt cuts */
57 #define DEFAULT_ONLYCONTROWS FALSE /**< whether only continuous rows should be used for rlt cuts */
58 #define DEFAULT_ONLYORIGINAL TRUE /**< whether only original variables and rows should be used for rlt cuts */
59 #define DEFAULT_USEINSUBSCIP FALSE /**< whether the separator should also be used in sub-scips */
60 #define DEFAULT_USEPROJECTION FALSE /**< whether the separator should first check projected rows */
61 #define DEFAULT_DETECTHIDDEN FALSE /**< whether implicit products should be detected and separated by McCormick */
62 #define DEFAULT_HIDDENRLT FALSE /**< whether RLT cuts should be added for hidden products */
63 #define DEFAULT_ADDTOPOOL TRUE /**< whether globally valid RLT cuts are added to the global cut pool */
64
65 #define DEFAULT_GOODSCORE 1.0 /**< threshold for score of cut relative to best score to be considered good,
66 * so that less strict filtering is applied */
67 #define DEFAULT_BADSCORE 0.5 /**< threshold for score of cut relative to best score to be discarded */
68 #define DEFAULT_OBJPARALWEIGHT 0.0 /**< weight of objective parallelism in cut score calculation */
69 #define DEFAULT_EFFICACYWEIGHT 1.0 /**< weight of efficacy in cut score calculation */
70 #define DEFAULT_DIRCUTOFFDISTWEIGHT 0.0 /**< weight of directed cutoff distance in cut score calculation */
71 #define DEFAULT_GOODMAXPARALL 0.1 /**< maximum parallelism for good cuts */
72 #define DEFAULT_MAXPARALL 0.1 /**< maximum parallelism for non-good cuts */
73
74 #define MAXVARBOUND 1e+5 /**< maximum allowed variable bound for computing an RLT-cut */
75
76 /*
77 * Data structures
78 */
79
80 /** data object for pairs and triples of variables */
81 struct HashData
82 {
83 SCIP_VAR* vars[3]; /**< variables in the pair or triple, used for hash comparison */
84 int nvars; /**< number of variables */
85 int nrows; /**< number of rows */
86 int firstrow; /**< beginning of the corresponding row linked list */
87 };
88 typedef struct HashData HASHDATA;
89
90 /** data structure representing an array of variables together with number of elements and size;
91 * used for storing variables that are in some sense adjacent to a given variable
92 */
93 struct AdjacentVarData
94 {
95 SCIP_VAR** adjacentvars; /**< adjacent vars */
96 int nadjacentvars; /**< number of vars in adjacentvars */
97 int sadjacentvars; /**< size of adjacentvars */
98 };
99 typedef struct AdjacentVarData ADJACENTVARDATA;
100
101 /** separator data */
102 struct SCIP_SepaData
103 {
104 SCIP_CONSHDLR* conshdlr; /**< nonlinear constraint handler */
105 SCIP_Bool iscreated; /**< indicates whether the sepadata has been initialized yet */
106 SCIP_Bool isinitialround; /**< indicates that this is the first round and original rows are used */
107
108 /* bilinear variables */
109 SCIP_VAR** varssorted; /**< variables that occur in bilinear terms sorted by priority */
110 SCIP_HASHMAP* bilinvardatamap; /**< maps each bilinear var to ADJACENTVARDATA containing vars appearing
111 together with it in bilinear products */
112 int* varpriorities; /**< priorities of variables */
113 int nbilinvars; /**< total number of variables occurring in bilinear terms */
114 int sbilinvars; /**< size of arrays for variables occurring in bilinear terms */
115
116 /* information about bilinear terms */
117 int* eqauxexpr; /**< position of the auxexpr that is equal to the product (-1 if none) */
118 int nbilinterms; /**< total number of bilinear terms */
119
120 /* parameters */
121 int maxunknownterms; /**< maximum number of unknown bilinear terms a row can have to be used (-1: unlimited) */
122 int maxusedvars; /**< maximum number of variables that will be used to compute rlt cuts (-1: unlimited) */
123 int maxncuts; /**< maximum number of cuts that will be added per round (-1: unlimited) */
124 int maxrounds; /**< maximum number of separation rounds per node (-1: unlimited) */
125 int maxroundsroot; /**< maximum number of separation rounds in the root node (-1: unlimited) */
126 SCIP_Bool onlyeqrows; /**< whether only equality rows should be used for rlt cuts */
127 SCIP_Bool onlycontrows; /**< whether only continuous rows should be used for rlt cuts */
128 SCIP_Bool onlyoriginal; /**< whether only original rows and variables should be used for rlt cuts */
129 SCIP_Bool useinsubscip; /**< whether the separator should also be used in sub-scips */
130 SCIP_Bool useprojection; /**< whether the separator should first check projected rows */
131 SCIP_Bool detecthidden; /**< whether implicit products should be detected and separated by McCormick */
132 SCIP_Bool hiddenrlt; /**< whether RLT cuts should be added for hidden products */
133 SCIP_Bool addtopool; /**< whether globally valid RLT cuts are added to the global cut pool */
134
135 /* cut selection parameters */
136 SCIP_Real goodscore; /**< threshold for score of cut relative to best score to be considered good,
137 * so that less strict filtering is applied */
138 SCIP_Real badscore; /**< threshold for score of cut relative to best score to be discarded */
139 SCIP_Real objparalweight; /**< weight of objective parallelism in cut score calculation */
140 SCIP_Real efficacyweight; /**< weight of efficacy in cut score calculation */
141 SCIP_Real dircutoffdistweight;/**< weight of directed cutoff distance in cut score calculation */
142 SCIP_Real goodmaxparall; /**< maximum parallelism for good cuts */
143 SCIP_Real maxparall; /**< maximum parallelism for non-good cuts */
144 };
145
146 /* a simplified representation of an LP row */
147 struct RLT_SimpleRow
148 {
149 const char* name; /**< name of the row */
150 SCIP_Real* coefs; /**< coefficients */
151 SCIP_VAR** vars; /**< variables */
152 SCIP_Real rhs; /**< right hand side */
153 SCIP_Real lhs; /**< left hand side */
154 SCIP_Real cst; /**< constant */
155 int nnonz; /**< number of nonzeroes */
156 int size; /**< size of the coefs and vars arrays */
157 };
158 typedef struct RLT_SimpleRow RLT_SIMPLEROW;
159
160 /*
161 * Local methods
162 */
163
164 /** returns TRUE iff both keys are equal
165 *
166 * two variable pairs/triples are equal if the variables are equal
167 */
168 static
169 SCIP_DECL_HASHKEYEQ(hashdataKeyEqConss)
170 { /*lint --e{715}*/
171 HASHDATA* hashdata1;
172 HASHDATA* hashdata2;
173 int v;
174
175 hashdata1 = (HASHDATA*)key1;
176 hashdata2 = (HASHDATA*)key2;
177
178 /* check data structure */
179 assert(hashdata1->nvars == hashdata2->nvars);
180 assert(hashdata1->firstrow != -1 || hashdata2->firstrow != -1);
181
182 for( v = hashdata1->nvars-1; v >= 0; --v )
183 {
184 /* tests if variables are equal */
185 if( hashdata1->vars[v] != hashdata2->vars[v] )
186 return FALSE;
187
188 assert(SCIPvarCompare(hashdata1->vars[v], hashdata2->vars[v]) == 0);
189 }
190
191 /* if two hashdata objects have the same variables, then either one of them doesn't have a row list yet
192 * (firstrow == -1) or they both point to the same row list
193 */
194 assert(hashdata1->firstrow == -1 || hashdata2->firstrow == -1 || hashdata1->firstrow == hashdata2->firstrow);
195
196 return TRUE;
197 }
198
199 /** returns the hash value of the key */
200 static
201 SCIP_DECL_HASHKEYVAL(hashdataKeyValConss)
202 { /*lint --e{715}*/
203 HASHDATA* hashdata;
204 int minidx;
205 int mididx;
206 int maxidx;
207 int idx[3];
208
209 hashdata = (HASHDATA*)key;
210 assert(hashdata != NULL);
211 assert(hashdata->nvars == 3 || hashdata->nvars == 2);
212
213 idx[0] = SCIPvarGetIndex(hashdata->vars[0]);
214 idx[1] = SCIPvarGetIndex(hashdata->vars[1]);
215 idx[2] = SCIPvarGetIndex(hashdata->vars[hashdata->nvars - 1]);
216
217 minidx = MIN3(idx[0], idx[1], idx[2]);
218 maxidx = MAX3(idx[0], idx[1], idx[2]);
219 if( idx[0] == maxidx )
220 mididx = MAX(idx[1], idx[2]);
221 else
222 mididx = MAX(idx[0], MIN(idx[1], idx[2]));
223
224 /* vars should already be sorted by index */
225 assert(minidx <= mididx && mididx <= maxidx);
226
227 return SCIPhashFour(hashdata->nvars, minidx, mididx, maxidx);
228 }
229
230 /** store a pair of adjacent variables */
231 static
232 SCIP_RETCODE addAdjacentVars(
233 SCIP* scip, /**< SCIP data structure */
234 SCIP_HASHMAP* adjvarmap, /**< hashmap mapping variables to their ADJACENTVARDATAs */
235 SCIP_VAR** vars /**< variable pair to be stored */
236 )
237 {
238 int v1;
239 int v2;
240 SCIP_Bool found;
241 int pos2;
242 int i;
243 ADJACENTVARDATA* adjacentvardata;
244
245 assert(adjvarmap != NULL);
246
247 /* repeat for each variable of the new pair */
248 for( v1 = 0; v1 < 2; ++v1 )
249 {
250 v2 = 1 - v1;
251
252 /* look for data of the first variable */
253 adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapGetImage(adjvarmap, (void*)(size_t) SCIPvarGetIndex(vars[v1]));
254
255 /* if the first variable has not been added to adjvarmap yet, add it here */
256 if( adjacentvardata == NULL )
257 {
258 SCIP_CALL( SCIPallocClearBlockMemory(scip, &adjacentvardata) );
259 SCIP_CALL( SCIPhashmapInsert(adjvarmap, (void*)(size_t) SCIPvarGetIndex(vars[v1]), adjacentvardata) );
260 }
261
262 assert(adjacentvardata != NULL);
263
264 /* look for second variable in adjacentvars */
265 if( adjacentvardata->adjacentvars == NULL )
266 {
267 found = FALSE;
268 pos2 = 0;
269 }
270 else
271 {
272 found = SCIPsortedvecFindPtr((void**) adjacentvardata->adjacentvars, SCIPvarComp, vars[v2],
273 adjacentvardata->nadjacentvars, &pos2);
274 }
275
276 /* add second var to adjacentvardata->adjacentvars, if not already added */
277 if( !found )
278 {
279 /* ensure size of adjacentvardata->adjacentvars */
280 SCIP_CALL( SCIPensureBlockMemoryArray(scip, &adjacentvardata->adjacentvars, &adjacentvardata->sadjacentvars,
281 adjacentvardata->nadjacentvars + 1) );
282
283 /* insert second var at the correct position */
284 for( i = adjacentvardata->nadjacentvars; i > pos2; --i )
285 adjacentvardata->adjacentvars[i] = adjacentvardata->adjacentvars[i-1];
286 adjacentvardata->adjacentvars[pos2] = vars[v2];
287 ++adjacentvardata->nadjacentvars;
288 }
289
290 /* if this is a self-adjacent var, only need to add the connection once */
291 if( vars[v1] == vars[v2] )
292 break;
293 }
294
295 return SCIP_OKAY;
296 }
297
298 /** returns the array of adjacent variables for a given variable */
299 static
300 SCIP_VAR** getAdjacentVars(
301 SCIP_HASHMAP* adjvarmap, /**< hashmap mapping variables to their ADJACENTVARDATAs */
302 SCIP_VAR* var, /**< variable */
303 int* nadjacentvars /**< buffer to store the number of variables in the returned array */
304 )
305 {
306 ADJACENTVARDATA* adjacentvardata;
307
308 assert(adjvarmap != NULL);
309
310 *nadjacentvars = 0;
311 adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapGetImage(adjvarmap, (void*)(size_t) SCIPvarGetIndex(var));
312
313 if( adjacentvardata == NULL )
314 return NULL;
315
316 *nadjacentvars = adjacentvardata->nadjacentvars;
317
318 return adjacentvardata->adjacentvars;
319 }
320
321 /** frees all ADJACENTVARDATAs stored in a hashmap */
322 static
323 void clearVarAdjacency(
324 SCIP* scip, /**< SCIP data structure */
325 SCIP_HASHMAP* adjvarmap /**< hashmap mapping variables to their ADJACENTVARDATAs */
326 )
327 {
328 int i;
329 SCIP_HASHMAPENTRY* entry;
330 ADJACENTVARDATA* adjacentvardata;
331
332 assert(adjvarmap != NULL);
333
334 for( i = 0; i < SCIPhashmapGetNEntries(adjvarmap); ++i )
335 {
336 entry = SCIPhashmapGetEntry(adjvarmap, i);
337
338 if( entry == NULL )
339 continue;
340
341 adjacentvardata = (ADJACENTVARDATA*) SCIPhashmapEntryGetImage(entry);
342
343 /* if adjacentvardata has been added to the hashmap, it can't be empty */
344 assert(adjacentvardata->adjacentvars != NULL);
345
346 SCIPfreeBlockMemoryArray(scip, &adjacentvardata->adjacentvars, adjacentvardata->sadjacentvars);
347 SCIPfreeBlockMemory(scip, &adjacentvardata);
348 }
349 }
350
351 /** free separator data */
352 static
353 SCIP_RETCODE freeSepaData(
354 SCIP* scip, /**< SCIP data structure */
355 SCIP_SEPADATA* sepadata /**< separation data */
356 )
357 { /*lint --e{715}*/
358 int i;
359
360 assert(sepadata->iscreated);
361
362 if( sepadata->nbilinvars != 0 )
363 {
364 /* release bilinvars that were captured for rlt and free all related arrays */
365
366 /* if there are bilinear vars, some of them must also participate in the same product */
367 assert(sepadata->bilinvardatamap != NULL);
368
369 clearVarAdjacency(scip, sepadata->bilinvardatamap);
370
371 for( i = 0; i < sepadata->nbilinvars; ++i )
372 {
373 assert(sepadata->varssorted[i] != NULL);
374 SCIP_CALL( SCIPreleaseVar(scip, &(sepadata->varssorted[i])) );
375 }
376
377 SCIPhashmapFree(&sepadata->bilinvardatamap);
378 SCIPfreeBlockMemoryArray(scip, &sepadata->varssorted, sepadata->sbilinvars);
379 SCIPfreeBlockMemoryArray(scip, &sepadata->varpriorities, sepadata->sbilinvars);
380 sepadata->nbilinvars = 0;
381 sepadata->sbilinvars = 0;
382 }
383
384 /* free the remaining array */
385 if( sepadata->nbilinterms > 0 )
386 {
387 SCIPfreeBlockMemoryArray(scip, &sepadata->eqauxexpr, sepadata->nbilinterms);
388 }
389
390 sepadata->iscreated = FALSE;
391
392 return SCIP_OKAY;
393 }
394
395 /** creates and returns rows of original linear constraints */
396 static
397 SCIP_RETCODE getOriginalRows(
398 SCIP* scip, /**< SCIP data structure */
399 SCIP_ROW*** rows, /**< buffer to store the rows */
400 int* nrows /**< buffer to store the number of linear rows */
401 )
402 {
403 SCIP_CONS** conss;
404 int nconss;
405 int i;
406
407 assert(rows != NULL);
408 assert(nrows != NULL);
409
410 conss = SCIPgetConss(scip);
411 nconss = SCIPgetNConss(scip);
412 *nrows = 0;
413
414 SCIP_CALL( SCIPallocBufferArray(scip, rows, nconss) );
415
416 for( i = 0; i < nconss; ++i )
417 {
418 SCIP_ROW *row;
419
420 row = SCIPconsGetRow(scip, conss[i]);
421
422 if( row != NULL )
423 {
424 (*rows)[*nrows] = row;
425 ++*nrows;
426 }
427 }
428
429 return SCIP_OKAY;
430 }
431
432 /** fills an array of rows suitable for RLT cut generation */
433 static
434 SCIP_RETCODE storeSuitableRows(
435 SCIP* scip, /**< SCIP data structure */
436 SCIP_SEPA* sepa, /**< separator */
437 SCIP_SEPADATA* sepadata, /**< separator data */
438 SCIP_ROW** prob_rows, /**< problem rows */
439 SCIP_ROW** rows, /**< an array to be filled with suitable rows */
440 int* nrows, /**< buffer to store the number of suitable rows */
441 SCIP_HASHMAP* row_to_pos, /**< hashmap linking row indices to positions in rows */
442 SCIP_Bool allowlocal /**< are local rows allowed? */
443 )
444 {
445 int new_nrows;
446 int r;
447 int j;
448 SCIP_Bool iseqrow;
449 SCIP_COL** cols;
450 SCIP_Bool iscontrow;
451
452 new_nrows = 0;
453
454 for( r = 0; r < *nrows; ++r )
455 {
456 iseqrow = SCIPisEQ(scip, SCIProwGetLhs(prob_rows[r]), SCIProwGetRhs(prob_rows[r]));
457
458 /* if equality rows are requested, only those can be used */
459 if( sepadata->onlyeqrows && !iseqrow )
460 continue;
461
462 /* if global cuts are requested, only globally valid rows can be used */
463 if( !allowlocal && SCIProwIsLocal(prob_rows[r]) )
464 continue;
465
466 /* if continuous rows are requested, only those can be used */
467 if( sepadata->onlycontrows )
468 {
469 cols = SCIProwGetCols(prob_rows[r]);
470 iscontrow = TRUE;
471
472 /* check row for integral variables */
473 for( j = 0; j < SCIProwGetNNonz(prob_rows[r]); ++j )
474 {
475 if( SCIPcolIsIntegral(cols[j]) )
476 {
477 iscontrow = FALSE;
478 break;
479 }
480 }
481
482 if( !iscontrow )
483 continue;
484 }
485
486 /* don't try to use rows that have been generated by the RLT separator */
487 if( SCIProwGetOriginSepa(prob_rows[r]) == sepa )
488 continue;
489
490 /* if we are here, the row has passed all checks and should be added to rows */
491 rows[new_nrows] = prob_rows[r];
492 SCIP_CALL( SCIPhashmapSetImageInt(row_to_pos, (void*)(size_t)SCIProwGetIndex(prob_rows[r]), new_nrows) ); /*lint !e571 */
493 ++new_nrows;
494 }
495
496 *nrows = new_nrows;
497
498 return SCIP_OKAY;
499 }
500
501 /** make sure that the arrays in sepadata are large enough to store information on n variables */
502 static
503 SCIP_RETCODE ensureVarsSize(
504 SCIP* scip, /**< SCIP data structure */
505 SCIP_SEPADATA* sepadata, /**< separator data */
506 int n /**< number of variables that we need to store */
507 )
508 {
509 int newsize;
510
511 /* check whether array is large enough */
512 if( n <= sepadata->sbilinvars )
513 return SCIP_OKAY;
514
515 /* compute new size */
516 newsize = SCIPcalcMemGrowSize(scip, n);
517 assert(n <= newsize);
518
519 /* realloc arrays */
520 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->varssorted, sepadata->sbilinvars, newsize) );
521 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &sepadata->varpriorities, sepadata->sbilinvars, newsize) );
522
523 sepadata->sbilinvars = newsize;
524
525 return SCIP_OKAY;
526 }
527
528 /** saves variables x and y to separator data and stores information about their connection
529 *
530 * variables must be captured separately
531 */
532 static
533 SCIP_RETCODE addProductVars(
534 SCIP* scip, /**< SCIP data structure */
535 SCIP_SEPADATA* sepadata, /**< separator data */
536 SCIP_VAR* x, /**< x variable */
537 SCIP_VAR* y, /**< y variable */
538 SCIP_HASHMAP* varmap, /**< hashmap linking var index to position */
539 int nlocks /**< number of locks */
540 )
541 {
542 int xpos;
543 int ypos;
544 int xidx;
545 int yidx;
546 SCIP_VAR* vars[2];
547
548 if( sepadata->bilinvardatamap == NULL )
549 {
550 int varmapsize;
551 int nvars;
552
553 /* the number of variables participating in bilinear products cannot exceed twice the number of bilinear terms;
554 * however, if we detect hidden products, the number of terms is yet unknown, so use the number of variables
555 */
556 nvars = SCIPgetNVars(scip);
557 varmapsize = sepadata->detecthidden ? nvars : MIN(nvars, sepadata->nbilinterms * 2);
558
559 SCIP_CALL( SCIPhashmapCreate(&sepadata->bilinvardatamap, SCIPblkmem(scip), varmapsize) );
560 }
561
562 xidx = SCIPvarGetIndex(x);
563 yidx = SCIPvarGetIndex(y);
564
565 xpos = SCIPhashmapGetImageInt(varmap, (void*)(size_t) xidx); /*lint !e571 */
566
567 if( xpos == INT_MAX )
568 {
569 /* add x to sepadata and initialise its priority */
570 SCIP_CALL( SCIPhashmapInsertInt(varmap, (void*)(size_t) xidx, sepadata->nbilinvars) ); /*lint !e571*/
571 SCIP_CALL( ensureVarsSize(scip, sepadata, sepadata->nbilinvars + 1) );
572 sepadata->varssorted[sepadata->nbilinvars] = x;
573 sepadata->varpriorities[sepadata->nbilinvars] = 0;
574 xpos = sepadata->nbilinvars;
575 ++sepadata->nbilinvars;
576 }
577
578 assert(xpos >= 0 && xpos < sepadata->nbilinvars );
579 assert(xpos == SCIPhashmapGetImageInt(varmap, (void*)(size_t) xidx)); /*lint !e571 */
580
581 /* add locks to priority of x */
582 sepadata->varpriorities[xpos] += nlocks;
583
584 if( xidx != yidx )
585 {
586 ypos = SCIPhashmapGetImageInt(varmap, (void*)(size_t) yidx); /*lint !e571 */
587
588 if( ypos == INT_MAX )
589 {
590 /* add y to sepadata and initialise its priority */
591 SCIP_CALL( SCIPhashmapInsertInt(varmap, (void*)(size_t) yidx, sepadata->nbilinvars) ); /*lint !e571*/
592 SCIP_CALL( ensureVarsSize(scip, sepadata, sepadata->nbilinvars + 1) );
593 sepadata->varssorted[sepadata->nbilinvars] = y;
594 sepadata->varpriorities[sepadata->nbilinvars] = 0;
595 ypos = sepadata->nbilinvars;
596 ++sepadata->nbilinvars;
597 }
598
599 assert(ypos >= 0 && ypos < sepadata->nbilinvars);
600 assert(ypos == SCIPhashmapGetImageInt(varmap, (void*)(size_t) yidx)); /*lint !e571 */
601
602 /* add locks to priority of y */
603 sepadata->varpriorities[ypos] += nlocks;
604 }
605
606 /* remember the connection between x and y */
607 vars[0] = x;
608 vars[1] = y;
609 SCIP_CALL( addAdjacentVars(scip, sepadata->bilinvardatamap, vars) );
610
611 return SCIP_OKAY;
612 }
613
614 /** extract a bilinear product from two linear relations, if possible
615 *
616 * First, the two given rows are brought to the form:
617 * \f[
618 * a_1x + b_1w + c_1y \leq/\geq d_1,\\
619 * a_2x + b_2w + c_2y \leq/\geq d_2,
620 * \f]
621 * where \f$ a_1a_2 \leq 0 \f$ and the first implied relation is enabled when \f$ x = 1 \f$
622 * and the second when \f$ x = 0 \f$, and \f$ b_1, b_2 > 0 \f$, the product relation can be written as:
623 * \f[
624 * \frac{b_1b_2w + (b_2(a_1 - d_1) + b_1d_2)x + b_1c_2y - b_1d_2}{b_1c_2 - c_1b_2} \leq/\geq xy.
625 * \f]
626 * The inequality sign in the product relation is similar to that in the given linear relations if
627 * \f$ b_1c_2 - c_1b_2 > 0 \f$ and opposite if \f$ b_1c_2 - c_1b_2 > 0 \f$.
628 *
629 * To obtain this formula, the given relations are first multiplied by scaling factors \f$ \alpha \f$
630 * and \f$ \beta \f$, which is necessary in order for the solution to always exist, and written as
631 * implications:
632 * \f{align}{
633 * x = 1 & ~\Rightarrow~ \alpha b_1w + \alpha c_1y \leq/\geq \alpha(d_1 - a_1), \\
634 * x = 0 & ~\Rightarrow~ \beta b_2w + \beta c_2y \leq/\geq \beta d_2.
635 * \f}
636 * Then a linear system is solved which ensures that the coefficients of the two implications of the product
637 * relation are equal to the corresponding coefficients in the linear relations.
638 * If the product relation is written as:
639 * \f[
640 * Ax + Bw + Cy + D \leq/\geq xy,
641 * \f]
642 * then the system is
643 * \f[
644 * B = \alpha b_1, ~C - 1 = \alpha c_1, ~D+A = \alpha(a_1-d_1),\\
645 * B = \beta b_2, ~C = \beta c_2, ~D = -\beta d_2.
646 * \f]
647 */
648 static
649 SCIP_RETCODE extractProducts(
650 SCIP* scip, /**< SCIP data structure */
651 SCIP_SEPADATA* sepadata, /**< separator data */
652 SCIP_VAR** vars_xwy, /**< 3 variables involved in the inequalities in the order x,w,y */
653 SCIP_Real* coefs1, /**< coefficients of the first inequality (always implied, i.e. has x) */
654 SCIP_Real* coefs2, /**< coefficients of the second inequality (can be unconditional) */
655 SCIP_Real d1, /**< side of the first inequality */
656 SCIP_Real d2, /**< side of the second inequality */
657 SCIP_SIDETYPE sidetype1, /**< side type (lhs or rls) in the first inequality */
658 SCIP_SIDETYPE sidetype2, /**< side type (lhs or rhs) in the second inequality */
659 SCIP_HASHMAP* varmap, /**< variable map */
660 SCIP_Bool f /**< the first relation is an implication x == f */
661 )
662 {
663 SCIP_Real mult;
664
665 /* coefficients and constant of the auxexpr */
666 SCIP_Real A; /* coefficient of x */
667 SCIP_Real B; /* coefficient of w */
668 SCIP_Real C; /* coefficient of y */
669 SCIP_Real D; /* constant */
670
671 /* variables */
672 SCIP_VAR* w;
673 SCIP_VAR* x;
674 SCIP_VAR* y;
675
676 /* does auxexpr overestimate the product? */
677 SCIP_Bool overestimate;
678
679 /* coefficients in given relations: a for x, b for w, c for y; 1 and 2 for 1st and 2nd relation, respectively */
680 SCIP_Real a1 = coefs1[0];
681 SCIP_Real b1 = coefs1[1];
682 SCIP_Real c1 = coefs1[2];
683 SCIP_Real a2 = coefs2[0];
684 SCIP_Real b2 = coefs2[1];
685 SCIP_Real c2 = coefs2[2];
686
687 x = vars_xwy[0];
688 w = vars_xwy[1];
689 y = vars_xwy[2];
690
691 /* check given linear relations and decide if to continue */
692
693 assert(SCIPvarGetType(x) == SCIP_VARTYPE_BINARY); /* x must be binary */
694 assert(a1 != 0.0); /* the first relation is always conditional */
695 assert(b1 != 0.0 || b2 != 0.0); /* at least one w coefficient must be nonzero */
696
697 SCIPdebugMsg(scip, "Extracting product from two implied relations:\n");
698 SCIPdebugMsg(scip, "Relation 1: <%s> == %u => %g<%s> + %g<%s> %s %g\n", SCIPvarGetName(x), f, b1,
699 SCIPvarGetName(w), c1, SCIPvarGetName(y), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
700 f ? d1 - a1 : d1);
701 SCIPdebugMsg(scip, "Relation 2: <%s> == %d => %g<%s> + %g<%s> %s %g\n", SCIPvarGetName(x), !f, b2,
702 SCIPvarGetName(w), c2, SCIPvarGetName(y), sidetype2 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
703 f ? d2 : d2 - a2);
704
705 /* cannot use a global bound on x to detect a product */
706 if( (b1 == 0.0 && c1 == 0.0) || (b2 == 0.0 && c2 == 0.0) )
707 return SCIP_OKAY;
708
709 /* cannot use a global bound on y to detect a non-redundant product relation */
710 if( a2 == 0.0 && b2 == 0.0 ) /* only check the 2nd relation because the 1st at least has x */
711 {
712 SCIPdebugMsg(scip, "Ignoring a global bound on y\n");
713 return SCIP_OKAY;
714 }
715
716 SCIPdebugMsg(scip, "binary var = <%s>, product of its coefs: %g\n", SCIPvarGetName(x), a1*a2);
717
718 /* rewrite the linear relations in a standard form:
719 * a1x + b1w + c1y <=/>= d1,
720 * a2x + b2w + c2y <=/>= d2,
721 * where b1 > 0, b2 > 0 and first implied relation is activated when x == 1
722 */
723
724 /* if needed, multiply the rows by -1 so that coefs of w are positive */
725 if( b1 < 0 )
726 {
727 a1 *= -1.0;
728 b1 *= -1.0;
729 c1 *= -1.0;
730 d1 *= -1.0;
731 sidetype1 = sidetype1 == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
732 }
733 if( b2 < 0 )
734 {
735 a2 *= -1.0;
736 b2 *= -1.0;
737 c2 *= -1.0;
738 d2 *= -1.0;
739 sidetype2 = sidetype2 == SCIP_SIDETYPE_LEFT ? SCIP_SIDETYPE_RIGHT : SCIP_SIDETYPE_LEFT;
740 }
741
742 /* the linear relations imply a product only if the inequality signs are similar */
743 if( sidetype1 != sidetype2 )
744 return SCIP_OKAY;
745
746 /* when b1c2 = b2c1, the linear relations do not imply a product relation */
747 if( SCIPisRelEQ(scip, b2*c1, c2*b1) )
748 {
749 SCIPdebugMsg(scip, "Ignoring a pair of linear relations because b1c2 = b2c1\n");
750 return SCIP_OKAY;
751 }
752
753 if( !f )
754 {
755 /* swap the linear relations so that the relation implied by x == TRUE goes first */
756 SCIPswapReals(&a1, &a2);
757 SCIPswapReals(&b1, &b2);
758 SCIPswapReals(&c1, &c2);
759 SCIPswapReals(&d1, &d2);
760 }
761
762 /* all conditions satisfied, we can extract the product and write it as:
763 * (1/(b1c2 - c1b2))*(b1b2w + (b2(a1 - d1) + b1d2)x + b1c2y - b1d2) >=/<= xy,
764 * where the inequality sign in the product relation is similar to that in the given linear relations
765 * if b1c2 - c1b2 > 0 and opposite if b1c2 - c1b2 > 0
766 */
767
768 /* compute the multiplier */
769 mult = 1/(b1*c2 - c1*b2);
770
771 /* determine the inequality sign; only check sidetype1 because sidetype2 is equal to it */
772 overestimate = (sidetype1 == SCIP_SIDETYPE_LEFT && mult > 0.0) || (sidetype1 == SCIP_SIDETYPE_RIGHT && mult < 0.0);
773
774 SCIPdebugMsg(scip, "found suitable implied rels (w,x,y): %g<%s> + %g<%s> + %g<%s> <= %g\n", a1,
775 SCIPvarGetName(x), b1, SCIPvarGetName(w), c1, SCIPvarGetName(y), d1);
776 SCIPdebugMsg(scip, " and %g<%s> + %g<%s> + %g<%s> <= %g\n", a2, SCIPvarGetName(x),
777 b2, SCIPvarGetName(w), c2, SCIPvarGetName(y), d2);
778
779 /* compute the coefficients for x, w and y and the constant in auxexpr */
780 A = (b2*a1 - d1*b2 + d2*b1)*mult;
781 B = b1*b2*mult;
782 C = b1*c2*mult;
783 D = -b1*d2*mult;
784
785 SCIPdebugMsg(scip, "product: <%s><%s> %s %g<%s> + %g<%s> + %g<%s> + %g\n", SCIPvarGetName(x), SCIPvarGetName(y),
786 overestimate ? "<=" : ">=", A, SCIPvarGetName(x), B, SCIPvarGetName(w), C, SCIPvarGetName(y), D);
787
788 SCIP_CALL( addProductVars(scip, sepadata, x, y, varmap, 1) );
789 SCIP_CALL( SCIPinsertBilinearTermImplicitNonlinear(scip, sepadata->conshdlr, x, y, w, A, C, B, D, overestimate) );
790
791 return SCIP_OKAY;
792 }
793
794 /** convert an implied bound: `binvar` = `binval` ⇒ `implvar` ≤/≥ `implbnd` into a big-M constraint */
795 static
796 void implBndToBigM(
797 SCIP* scip, /**< SCIP data structure */
798 SCIP_VAR** vars_xwy, /**< variables in order x,w,y */
799 int binvarpos, /**< position of binvar in vars_xwy */
800 int implvarpos, /**< position of implvar in vars_xwy */
801 SCIP_BOUNDTYPE bndtype, /**< type of implied bound */
802 SCIP_Bool binval, /**< value of binvar which implies the bound */
803 SCIP_Real implbnd, /**< value of the implied bound */
804 SCIP_Real* coefs, /**< coefficients of the big-M constraint */
805 SCIP_Real* side /**< side of the big-M constraint */
806 )
807 {
808 SCIP_VAR* implvar;
809 SCIP_Real globbnd;
810
811 assert(vars_xwy != NULL);
812 assert(coefs != NULL);
813 assert(side != NULL);
814 assert(binvarpos != implvarpos);
815
816 implvar = vars_xwy[implvarpos];
817 globbnd = bndtype == SCIP_BOUNDTYPE_LOWER ? SCIPvarGetLbGlobal(implvar) : SCIPvarGetUbGlobal(implvar);
818
819 /* Depending on the bound type and binval, there are four possibilities:
820 * binvar == 1 => implvar >= implbnd <=> (implvar^l - implbnd)binvar + implvar >= implvar^l;
821 * binvar == 0 => implvar >= implbnd <=> (implbnd - implvar^l)binvar + implvar >= implbnd;
822 * binvar == 1 => implvar <= implbnd <=> (implvar^u - implbnd)binvar + implvar <= implvar^u;
823 * binvar == 0 => implvar <= implbnd <=> (implbnd - implvar^u)binvar + implvar <= implbnd.
824 */
825
826 coefs[0] = 0.0;
827 coefs[1] = 0.0;
828 coefs[2] = 0.0;
829 coefs[binvarpos] = binval ? globbnd - implbnd : implbnd - globbnd;
830 coefs[implvarpos] = 1.0;
831 *side = binval ? globbnd : implbnd;
832
833 SCIPdebugMsg(scip, "Got an implied relation with binpos = %d, implpos = %d, implbnd = %g, "
834 "bnd type = %s, binval = %u, glbbnd = %g\n", binvarpos, implvarpos, implbnd,
835 bndtype == SCIP_BOUNDTYPE_LOWER ? "lower" : "upper", binval, globbnd);
836 SCIPdebugMsg(scip, "Constructed big-M: %g*bvar + implvar %s %g\n", coefs[binvarpos],
837 bndtype == SCIP_BOUNDTYPE_LOWER ? ">=" : "<=", *side);
838 }
839
840 /** extract products from a relation given by coefs1, vars, side1 and sidetype1 and
841 * implied bounds of the form `binvar` = `!f` ⇒ `implvar` ≥/≤ `implbnd`
842 */
843 static
844 SCIP_RETCODE detectProductsImplbnd(
845 SCIP* scip, /**< SCIP data structure */
846 SCIP_SEPADATA* sepadata, /**< separator data */
847 SCIP_Real* coefs1, /**< coefficients of the first linear relation */
848 SCIP_VAR** vars_xwy, /**< variables in the order x, w, y */
849 SCIP_Real side1, /**< side of the first relation */
850 SCIP_SIDETYPE sidetype1, /**< is the left or right hand side given for the first relation? */
851 int binvarpos, /**< position of the indicator variable in the vars_xwy array */
852 int implvarpos, /**< position of the variable that is bounded */
853 SCIP_HASHMAP* varmap, /**< variable map */
854 SCIP_Bool f /**< the value of x that activates the first relation */
855 )
856 {
857 SCIP_Real coefs2[3] = { 0., 0., 0. };
858 SCIP_Real impllb;
859 SCIP_Real implub;
860 SCIP_VAR* binvar;
861 SCIP_VAR* implvar;
862 SCIP_Real side2;
863 int i;
864 SCIP_Bool binvals[2] = {!f, f};
865
866 assert(binvarpos != implvarpos);
867 assert(implvarpos != 0); /* implied variable must be continuous, therefore it can't be x */
868
869 binvar = vars_xwy[binvarpos];
870 implvar = vars_xwy[implvarpos];
871
872 assert(SCIPvarGetType(binvar) == SCIP_VARTYPE_BINARY);
873 assert(SCIPvarGetType(implvar) != SCIP_VARTYPE_BINARY);
874
875 /* loop over binvals; if binvar is x (case binvarpos == 0), then we want to use only implications from
876 * binvar == !f (which is the option complementing the first relation, which is implied from f); if
877 * binvar is not x, this doesn't matter since the implbnd doesn't depend on x, therefore try both !f and f
878 */
879 for( i = 0; i < (binvarpos == 0 ? 1 : 2); ++i )
880 {
881 /* get implications binvar == binval => implvar <=/>= implbnd */
882 SCIPvarGetImplicVarBounds(binvar, binvals[i], implvar, &impllb, &implub);
883
884 if( impllb != SCIP_INVALID ) /*lint !e777*/
885 {
886 /* write the implied bound as a big-M constraint */
887 implBndToBigM(scip, vars_xwy, binvarpos, implvarpos, SCIP_BOUNDTYPE_LOWER, binvals[i], impllb, coefs2, &side2);
888
889 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
890 SCIP_SIDETYPE_LEFT, varmap, f) );
891 }
892
893 if( implub != SCIP_INVALID ) /*lint !e777*/
894 {
895 /* write the implied bound as a big-M constraint */
896 implBndToBigM(scip, vars_xwy, binvarpos, implvarpos, SCIP_BOUNDTYPE_UPPER, binvals[i], implub, coefs2, &side2);
897
898 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
899 SCIP_SIDETYPE_RIGHT, varmap, f) );
900 }
901 }
902
903 return SCIP_OKAY;
904 }
905
906 /** extract products from a relation given by `coefs1`, `vars_xwy`, `side1` and `sidetype1` and
907 * cliques containing `vars_xwy[varpos1]` and `vars_xwy[varpos2]`
908 */
909 static
910 SCIP_RETCODE detectProductsClique(
911 SCIP* scip, /**< SCIP data structure */
912 SCIP_SEPADATA* sepadata, /**< separator data */
913 SCIP_Real* coefs1, /**< coefficients of the first linear relation */
914 SCIP_VAR** vars_xwy, /**< variables of the first relation in the order x, w, y */
915 SCIP_Real side1, /**< side of the first relation */
916 SCIP_SIDETYPE sidetype1, /**< is the left or right hand side given for the first relation? */
917 int varpos1, /**< position of the first variable in the vars_xwy array */
918 int varpos2, /**< position of the second variable in the vars_xwy array */
919 SCIP_HASHMAP* varmap, /**< variable map */
920 SCIP_Bool f /**< the value of x that activates the first relation */
921 )
922 {
923 SCIP_Real coefs2[3] = { 0., 0., 0. };
924 SCIP_VAR* var1;
925 SCIP_VAR* var2;
926 SCIP_Real side2;
927 int i;
928 int imax;
929 SCIP_Bool binvals[2] = {!f, f};
930
931 var1 = vars_xwy[varpos1];
932 var2 = vars_xwy[varpos2];
933
934 /* this decides whether we do one or two iterations of the loop for binvals: if var1
935 * or var2 is x, we only want cliques with x = !f (which is the option complementing
936 * the first relation, which is implied from f); otherwise this doesn't matter since
937 * the clique doesn't depend on x, therefore try both !f and f
938 */
939 imax = (varpos1 == 0 || varpos2 == 0) ? 1 : 2;
940
941 assert(SCIPvarGetType(var1) == SCIP_VARTYPE_BINARY);
942 assert(SCIPvarGetType(var2) == SCIP_VARTYPE_BINARY);
943
944 for( i = 0; i < imax; ++i )
945 {
946 /* if var1=TRUE and var2=TRUE are in a clique (binvals[i] == TRUE), the relation var1 + var2 <= 1 is implied
947 * if var1=FALSE and var2=TRUE are in a clique (binvals[i] == FALSE), the relation (1 - var1) + var2 <= 1 is implied
948 */
949 if( SCIPvarsHaveCommonClique(var1, binvals[i], var2, TRUE, TRUE) )
950 {
951 SCIPdebugMsg(scip, "vars %s<%s> and <%s> are in a clique\n", binvals[i] ? "" : "!", SCIPvarGetName(var1), SCIPvarGetName(var2));
952 coefs2[varpos1] = binvals[i] ? 1.0 : -1.0;
953 coefs2[varpos2] = 1.0;
954 side2 = binvals[i] ? 1.0 : 0.0;
955
956 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
957 SCIP_SIDETYPE_RIGHT, varmap, f) );
958 }
959
960 /* if var1=TRUE and var2=FALSE are in the same clique, the relation var1 + (1-var2) <= 1 is implied
961 * if var1=FALSE and var2=FALSE are in the same clique, the relation (1-var1) + (1-var2) <= 1 is implied
962 */
963 if( SCIPvarsHaveCommonClique(var1, binvals[i], var2, FALSE, TRUE) )
964 {
965 SCIPdebugMsg(scip, "vars %s<%s> and !<%s> are in a clique\n", binvals[i] ? "" : "!", SCIPvarGetName(var1), SCIPvarGetName(var2));
966 coefs2[varpos1] = binvals[i] ? 1.0 : -1.0;
967 coefs2[varpos2] = -1.0;
968 side2 = binvals[i] ? 0.0 : -1.0;
969
970 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
971 SCIP_SIDETYPE_RIGHT, varmap, f) );
972 }
973 }
974
975 return SCIP_OKAY;
976 }
977
978
979 /** extract products from a relation given by `coefs1`, `vars`, `side1` and `sidetype1` and unconditional relations
980 * (inequalities with 2 nonzeros) containing `vars[varpos1]` and `vars[varpos2]`
981 */
982 static
983 SCIP_RETCODE detectProductsUnconditional(
984 SCIP* scip, /**< SCIP data structure */
985 SCIP_SEPADATA* sepadata, /**< separator data */
986 SCIP_ROW** rows, /**< problem rows */
987 int* row_list, /**< linked list of rows corresponding to 2 or 3 var sets */
988 SCIP_HASHTABLE* hashtable, /**< hashtable storing unconditional relations */
989 SCIP_Real* coefs1, /**< coefficients of the first linear relation */
990 SCIP_VAR** vars_xwy, /**< variables of the first relation in the order x, w, y */
991 SCIP_Real side1, /**< side of the first relation */
992 SCIP_SIDETYPE sidetype1, /**< is the left or right hand side given for the first relation? */
993 int varpos1, /**< position of the first unconditional variable in the vars_xwy array */
994 int varpos2, /**< position of the second unconditional variable in the vars_xwy array */
995 SCIP_HASHMAP* varmap, /**< variable map */
996 SCIP_Bool f /**< the value of x that activates the first relation */
997 )
998 {
999 HASHDATA hashdata;
1000 HASHDATA* foundhashdata;
1001 SCIP_ROW* row2;
1002 int r2;
1003 int pos1;
1004 int pos2;
1005 SCIP_Real coefs2[3] = { 0., 0., 0. };
1006 SCIP_VAR* var1;
1007 SCIP_VAR* var2;
1008
1009 /* always unconditional, therefore x must not be one of the two variables */
1010 assert(varpos1 != 0);
1011 assert(varpos2 != 0);
1012
1013 var1 = vars_xwy[varpos1];
1014 var2 = vars_xwy[varpos2];
1015
1016 hashdata.nvars = 2;
1017 hashdata.firstrow = -1;
1018 if( SCIPvarGetIndex(var1) < SCIPvarGetIndex(var2) )
1019 {
1020 pos1 = 0;
1021 pos2 = 1;
1022 }
1023 else
1024 {
1025 pos1 = 1;
1026 pos2 = 0;
1027 }
1028
1029 hashdata.vars[pos1] = var1;
1030 hashdata.vars[pos2] = var2;
1031
1032 foundhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable, &hashdata);
1033
1034 if( foundhashdata != NULL )
1035 {
1036 /* if the var pair exists, use all corresponding rows */
1037 r2 = foundhashdata->firstrow;
1038
1039 while( r2 != -1 )
1040 {
1041 row2 = rows[r2];
1042 assert(SCIProwGetNNonz(row2) == 2);
1043 assert(var1 == SCIPcolGetVar(SCIProwGetCols(row2)[pos1]));
1044 assert(var2 == SCIPcolGetVar(SCIProwGetCols(row2)[pos2]));
1045
1046 coefs2[varpos1] = SCIProwGetVals(row2)[pos1];
1047 coefs2[varpos2] = SCIProwGetVals(row2)[pos2];
1048
1049 SCIPdebugMsg(scip, "Unconditional:\n");
1050 if( !SCIPisInfinity(scip, -SCIProwGetLhs(row2)) )
1051 {
1052 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1053 SCIProwGetLhs(row2) - SCIProwGetConstant(row2), sidetype1, SCIP_SIDETYPE_LEFT, varmap, f) );
1054 }
1055 if( !SCIPisInfinity(scip, SCIProwGetRhs(row2)) )
1056 {
1057 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1058 SCIProwGetRhs(row2) - SCIProwGetConstant(row2), sidetype1, SCIP_SIDETYPE_RIGHT, varmap, f) );
1059 }
1060
1061 r2 = row_list[r2];
1062 }
1063 }
1064
1065 return SCIP_OKAY;
1066 }
1067
1068 /** finds and stores implied relations (x = f ⇒ ay + bw ≤ c, f can be 0 or 1) and 2-variable relations
1069 *
1070 * Fills the following:
1071 *
1072 * - An array of variables that participate in two variable relations; for each such variable, ADJACENTVARDATA
1073 * containing an array of variables that participate in two variable relations together with it; and a hashmap
1074 * mapping variables to ADJACENTVARDATAs.
1075 *
1076 * - Hashtables storing hashdata objects with the two or three variables and the position of the first row in the
1077 * `prob_rows` array, which in combination with the linked list (described below) will allow access to all rows that
1078 * depend only on the corresponding variables.
1079 *
1080 * - Linked lists of row indices. Each list corresponds to a pair or triple of variables and contains positions of rows
1081 * which depend only on those variables. All lists are stored in `row_list`, an array of length `nrows`, which is
1082 * possible because each row belongs to at most one list. The array indices of `row_list` represent the positions of
1083 * rows in `prob_rows`, and a value in the `row_list` array represents the next index in the list (-1 if there is no next
1084 * list element). The first index of each list is stored in one of the hashdata objects as firstrow.
1085 */
1086 static
1087 SCIP_RETCODE fillRelationTables(
1088 SCIP* scip, /**< SCIP data structure */
1089 SCIP_ROW** prob_rows, /**< linear rows of the problem */
1090 int nrows, /**< number of rows */
1091 SCIP_HASHTABLE* hashtable2, /**< hashtable to store 2-variable relations */
1092 SCIP_HASHTABLE* hashtable3, /**< hashtable to store implied relations */
1093 SCIP_HASHMAP* vars_in_2rels, /**< connections between variables that appear in 2-variable relations */
1094 int* row_list /**< linked lists of row positions for each 2 or 3 variable set */
1095 )
1096 {
1097 int r;
1098 SCIP_COL** cols;
1099 HASHDATA searchhashdata;
1100 HASHDATA* elementhashdata;
1101
1102 assert(prob_rows != NULL);
1103 assert(nrows > 0);
1104 assert(hashtable2 != NULL);
1105 assert(hashtable3 != NULL);
1106 assert(vars_in_2rels != NULL);
1107 assert(row_list != NULL);
1108
1109 for( r = 0; r < nrows; ++r )
1110 {
1111 assert(prob_rows[r] != NULL);
1112
1113 cols = SCIProwGetCols(prob_rows[r]);
1114 assert(cols != NULL);
1115
1116 /* initialise with the "end of list" value */
1117 row_list[r] = -1;
1118
1119 /* look for unconditional relations with 2 variables */
1120 if( SCIProwGetNNonz(prob_rows[r]) == 2 )
1121 {
1122 /* if at least one of the variables is binary, this is either an implied bound
1123 * or a clique; these are covered separately */
1124 if( SCIPvarGetType(SCIPcolGetVar(cols[0])) == SCIP_VARTYPE_BINARY ||
1125 SCIPvarGetType(SCIPcolGetVar(cols[1])) == SCIP_VARTYPE_BINARY )
1126 {
1127 SCIPdebugMsg(scip, "ignoring relation <%s> because a var is binary\n", SCIProwGetName(prob_rows[r]));
1128 continue;
1129 }
1130
1131 /* fill in searchhashdata so that to search for the two variables in hashtable2 */
1132 searchhashdata.nvars = 2;
1133 searchhashdata.firstrow = -1;
1134 searchhashdata.vars[0] = SCIPcolGetVar(cols[0]);
1135 searchhashdata.vars[1] = SCIPcolGetVar(cols[1]);
1136
1137 /* get the element corresponding to the two variables */
1138 elementhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable2, &searchhashdata);
1139
1140 if( elementhashdata != NULL )
1141 {
1142 /* if element exists, update it by adding the row */
1143 row_list[r] = elementhashdata->firstrow;
1144 elementhashdata->firstrow = r;
1145 ++elementhashdata->nrows;
1146 }
1147 else
1148 {
1149 /* create an element for the combination of two variables */
1150 SCIP_CALL( SCIPallocBuffer(scip, &elementhashdata) );
1151
1152 elementhashdata->nvars = 2;
1153 elementhashdata->nrows = 1;
1154 elementhashdata->vars[0] = searchhashdata.vars[0];
1155 elementhashdata->vars[1] = searchhashdata.vars[1];
1156 elementhashdata->firstrow = r;
1157
1158 SCIP_CALL( SCIPhashtableInsert(hashtable2, (void*)elementhashdata) );
1159
1160 /* hashdata.vars are two variables participating together in a two variable relation, therefore update
1161 * these variables' adjacency data
1162 */
1163 SCIP_CALL( addAdjacentVars(scip, vars_in_2rels, searchhashdata.vars) );
1164 }
1165 }
1166
1167 /* look for implied relations (three variables, at least one binary variable) */
1168 if( SCIProwGetNNonz(prob_rows[r]) == 3 )
1169 {
1170 /* an implied relation contains at least one binary variable */
1171 if( SCIPvarGetType(SCIPcolGetVar(cols[0])) != SCIP_VARTYPE_BINARY &&
1172 SCIPvarGetType(SCIPcolGetVar(cols[1])) != SCIP_VARTYPE_BINARY &&
1173 SCIPvarGetType(SCIPcolGetVar(cols[2])) != SCIP_VARTYPE_BINARY )
1174 continue;
1175
1176 /* fill in hashdata so that to search for the three variables in hashtable3 */
1177 searchhashdata.nvars = 3;
1178 searchhashdata.firstrow = -1;
1179 searchhashdata.vars[0] = SCIPcolGetVar(cols[0]);
1180 searchhashdata.vars[1] = SCIPcolGetVar(cols[1]);
1181 searchhashdata.vars[2] = SCIPcolGetVar(cols[2]);
1182
1183 /* get the element corresponding to the three variables */
1184 elementhashdata = (HASHDATA*)SCIPhashtableRetrieve(hashtable3, &searchhashdata);
1185
1186 if( elementhashdata != NULL )
1187 {
1188 /* if element exists, update it by adding the row */
1189 row_list[r] = elementhashdata->firstrow;
1190 elementhashdata->firstrow = r;
1191 ++elementhashdata->nrows;
1192 }
1193 else
1194 {
1195 /* create an element for the combination of three variables */
1196 SCIP_CALL( SCIPallocBuffer(scip, &elementhashdata) );
1197
1198 elementhashdata->nvars = 3;
1199 elementhashdata->nrows = 1;
1200 elementhashdata->vars[0] = searchhashdata.vars[0];
1201 elementhashdata->vars[1] = searchhashdata.vars[1];
1202 elementhashdata->vars[2] = searchhashdata.vars[2];
1203 elementhashdata->firstrow = r;
1204
1205 SCIP_CALL( SCIPhashtableInsert(hashtable3, (void*)elementhashdata) );
1206 }
1207 }
1208 }
1209
1210 return SCIP_OKAY;
1211 }
1212
1213 /** detect bilinear products encoded in linear constraints */
1214 static
1215 SCIP_RETCODE detectHiddenProducts(
1216 SCIP* scip, /**< SCIP data structure */
1217 SCIP_SEPADATA* sepadata, /**< separation data */
1218 SCIP_HASHMAP* varmap /**< variable map */
1219 )
1220 {
1221 int r1; /* first relation index */
1222 int r2; /* second relation index */
1223 int i; /* outer loop counter */
1224 int permwy; /* index for permuting w and y */
1225 int nrows;
1226 SCIP_ROW** prob_rows;
1227 SCIP_HASHTABLE* hashtable3;
1228 SCIP_HASHTABLE* hashtable2;
1229 HASHDATA* foundhashdata;
1230 SCIP_VAR* vars_xwy[3];
1231 SCIP_Real coefs1[3];
1232 SCIP_Real coefs2[3];
1233 SCIP_ROW* row1;
1234 SCIP_ROW* row2;
1235 int xpos;
1236 int ypos;
1237 int wpos;
1238 int f; /* value of the binary variable */
1239 SCIP_VAR** relatedvars;
1240 int nrelatedvars;
1241 SCIP_Bool xfixing;
1242 SCIP_SIDETYPE sidetype1;
1243 SCIP_SIDETYPE sidetype2;
1244 SCIP_Real side1;
1245 SCIP_Real side2;
1246 int* row_list;
1247 SCIP_HASHMAP* vars_in_2rels;
1248 int nvars;
1249
1250 /* get the (original) rows */
1251 SCIP_CALL( getOriginalRows(scip, &prob_rows, &nrows) );
1252
1253 if( nrows == 0 )
1254 {
1255 SCIPfreeBufferArray(scip, &prob_rows);
1256 return SCIP_OKAY;
1257 }
1258
1259 /* create tables of implied and unconditional relations */
1260 SCIP_CALL( SCIPhashtableCreate(&hashtable3, SCIPblkmem(scip), nrows, SCIPhashGetKeyStandard,
1261 hashdataKeyEqConss, hashdataKeyValConss, NULL) );
1262 SCIP_CALL( SCIPhashtableCreate(&hashtable2, SCIPblkmem(scip), nrows, SCIPhashGetKeyStandard,
1263 hashdataKeyEqConss, hashdataKeyValConss, NULL) );
1264 SCIP_CALL( SCIPallocBufferArray(scip, &row_list, nrows) );
1265
1266 /* allocate the adjacency data map for variables that appear in 2-var relations */
1267 nvars = SCIPgetNVars(scip);
1268 SCIP_CALL( SCIPhashmapCreate(&vars_in_2rels, SCIPblkmem(scip), MIN(nvars, nrows * 2)) );
1269
1270 /* fill the data structures that will be used for product detection: hashtables and linked lists allowing to access
1271 * two and three variable relations by the variables; and the hashmap for accessing variables participating in two
1272 * variable relations with each given variable */
1273 SCIP_CALL( fillRelationTables(scip, prob_rows, nrows, hashtable2, hashtable3, vars_in_2rels, row_list) );
1274
1275 /* start actually looking for products */
1276 /* go through all sets of three variables */
1277 for( i = 0; i < SCIPhashtableGetNEntries(hashtable3); ++i )
1278 {
1279 foundhashdata = (HASHDATA*)SCIPhashtableGetEntry(hashtable3, i);
1280 if( foundhashdata == NULL )
1281 continue;
1282
1283 SCIPdebugMsg(scip, "(<%s>, <%s>, <%s>): ", SCIPvarGetName(foundhashdata->vars[0]),
1284 SCIPvarGetName(foundhashdata->vars[1]), SCIPvarGetName(foundhashdata->vars[2]));
1285
1286 /* An implied relation has the form: x == f => l(w,y) <=/>= side (f is 0 or 1, l is a linear function). Given
1287 * a linear relation with three variables, any binary var can be x: we try them all here because this can
1288 * produce different products.
1289 */
1290 for( xpos = 0; xpos < 3; ++xpos )
1291 {
1292 /* in vars_xwy, the order of variables is always as in the name: x, w, y */
1293 vars_xwy[0] = foundhashdata->vars[xpos];
1294
1295 /* x must be binary */
1296 if( SCIPvarGetType(vars_xwy[0]) != SCIP_VARTYPE_BINARY )
1297 continue;
1298
1299 /* the first row might be an implication from f == 0 or f == 1: try both */
1300 for( f = 0; f <= 1; ++f )
1301 {
1302 xfixing = f == 1;
1303
1304 /* go through implied relations for the corresponding three variables */
1305 for( r1 = foundhashdata->firstrow; r1 != -1; r1 = row_list[r1] )
1306 {
1307 /* get the implied relation */
1308 row1 = prob_rows[r1];
1309
1310 assert(SCIProwGetNNonz(row1) == 3);
1311 /* the order of variables in all rows should be the same, and similar to the order in hashdata->vars,
1312 * therefore the x variable from vars_xwy should be similar to the column variable at xpos
1313 */
1314 assert(vars_xwy[0] == SCIPcolGetVar(SCIProwGetCols(row1)[xpos]));
1315
1316 coefs1[0] = SCIProwGetVals(row1)[xpos];
1317
1318 /* use the side for which the inequality becomes tighter when x == xfixing than when x == !xfixing */
1319 if( (!xfixing && coefs1[0] > 0.0) || (xfixing && coefs1[0] < 0.0) )
1320 {
1321 sidetype1 = SCIP_SIDETYPE_LEFT;
1322 side1 = SCIProwGetLhs(row1);
1323 }
1324 else
1325 {
1326 sidetype1 = SCIP_SIDETYPE_RIGHT;
1327 side1 = SCIProwGetRhs(row1);
1328 }
1329
1330 if( SCIPisInfinity(scip, REALABS(side1)) )
1331 continue;
1332
1333 side1 -= SCIProwGetConstant(row1);
1334
1335 /* permute w and y */
1336 for( permwy = 1; permwy <= 2; ++permwy )
1337 {
1338 wpos = (xpos + permwy) % 3;
1339 ypos = (xpos - permwy + 3) % 3;
1340 vars_xwy[1] = foundhashdata->vars[wpos];
1341 vars_xwy[2] = foundhashdata->vars[ypos];
1342
1343 assert(vars_xwy[1] == SCIPcolGetVar(SCIProwGetCols(row1)[wpos]));
1344 assert(vars_xwy[2] == SCIPcolGetVar(SCIProwGetCols(row1)[ypos]));
1345
1346 coefs1[1] = SCIProwGetVals(row1)[wpos];
1347 coefs1[2] = SCIProwGetVals(row1)[ypos];
1348
1349 /* look for the second relation: it should be tighter when x == !xfixing than when x == xfixing
1350 * and can be either another implied relation or one of several types of two and one variable
1351 * relations
1352 */
1353
1354 /* go through the remaining rows (implied relations) for these three variables */
1355 for( r2 = row_list[r1]; r2 != -1; r2 = row_list[r2] )
1356 {
1357 /* get the second implied relation */
1358 row2 = prob_rows[r2];
1359
1360 assert(SCIProwGetNNonz(row2) == 3);
1361 assert(vars_xwy[0] == SCIPcolGetVar(SCIProwGetCols(row2)[xpos]));
1362 assert(vars_xwy[1] == SCIPcolGetVar(SCIProwGetCols(row2)[wpos]));
1363 assert(vars_xwy[2] == SCIPcolGetVar(SCIProwGetCols(row2)[ypos]));
1364
1365 coefs2[0] = SCIProwGetVals(row2)[xpos];
1366 coefs2[1] = SCIProwGetVals(row2)[wpos];
1367 coefs2[2] = SCIProwGetVals(row2)[ypos];
1368
1369 /* use the side for which the inequality becomes tighter when x == !xfixing than when x == xfixing */
1370 if( (!xfixing && coefs2[0] > 0.0) || (xfixing && coefs2[0] < 0.0) )
1371 {
1372 sidetype2 = SCIP_SIDETYPE_RIGHT;
1373 side2 = SCIProwGetRhs(row2);
1374 }
1375 else
1376 {
1377 sidetype2 = SCIP_SIDETYPE_LEFT;
1378 side2 = SCIProwGetLhs(row2);
1379 }
1380
1381 if( SCIPisInfinity(scip, REALABS(side2)) )
1382 continue;
1383
1384 side2 -= SCIProwGetConstant(row2);
1385
1386 SCIPdebugMsg(scip, "Two implied relations:\n");
1387 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1, side2, sidetype1,
1388 sidetype2, varmap, xfixing) );
1389 }
1390
1391 /* use global bounds on w */
1392 coefs2[0] = 0.0;
1393 coefs2[1] = 1.0;
1394 coefs2[2] = 0.0;
1395 SCIPdebugMsg(scip, "w global bounds:\n");
1396 if( !SCIPisInfinity(scip, -SCIPvarGetLbGlobal(vars_xwy[1])) )
1397 {
1398 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1399 SCIPvarGetLbGlobal(vars_xwy[1]), sidetype1, SCIP_SIDETYPE_LEFT, varmap, xfixing) );
1400 }
1401
1402 if( !SCIPisInfinity(scip, SCIPvarGetUbGlobal(vars_xwy[1])) )
1403 {
1404 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1405 SCIPvarGetUbGlobal(vars_xwy[1]), sidetype1, SCIP_SIDETYPE_RIGHT, varmap, xfixing) );
1406 }
1407
1408 /* use implied bounds and cliques with w */
1409 if( SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY )
1410 {
1411 /* w is non-binary - look for implied bounds x == !f => w >=/<= bound */
1412 SCIPdebugMsg(scip, "Implied relation + implied bounds on w:\n");
1413 SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 0, 1,
1414 varmap, xfixing) );
1415 }
1416 else
1417 {
1418 /* w is binary - look for cliques containing x and w */
1419 SCIPdebugMsg(scip, "Implied relation + cliques with x and w:\n");
1420 SCIP_CALL( detectProductsClique(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 0, 1,
1421 varmap, xfixing) );
1422 }
1423
1424 /* use unconditional relations (i.e. relations of w and y) */
1425
1426 /* implied bound w == 0/1 => y >=/<= bound */
1427 if( SCIPvarGetType(vars_xwy[1]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) != SCIP_VARTYPE_BINARY )
1428 {
1429 SCIPdebugMsg(scip, "Implied relation + implied bounds with w and y:\n");
1430 SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
1431 }
1432
1433 /* implied bound y == 0/1 => w >=/<= bound */
1434 if( SCIPvarGetType(vars_xwy[2]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY )
1435 {
1436 SCIPdebugMsg(scip, "Implied relation + implied bounds with y and w:\n");
1437 SCIP_CALL( detectProductsImplbnd(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 2, 1, varmap, xfixing) );
1438 }
1439
1440 /* cliques containing w and y */
1441 if( SCIPvarGetType(vars_xwy[1]) == SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) == SCIP_VARTYPE_BINARY )
1442 {
1443 SCIPdebugMsg(scip, "Implied relation + cliques with w and y:\n");
1444 SCIP_CALL( detectProductsClique(scip, sepadata, coefs1, vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
1445 }
1446
1447 /* inequalities containing w and y */
1448 if( SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY && SCIPvarGetType(vars_xwy[2]) != SCIP_VARTYPE_BINARY )
1449 {
1450 SCIPdebugMsg(scip, "Implied relation + unconditional with w and y:\n");
1451 SCIP_CALL( detectProductsUnconditional(scip, sepadata, prob_rows, row_list, hashtable2, coefs1,
1452 vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
1453 }
1454 }
1455 }
1456 }
1457 }
1458 SCIPfreeBuffer(scip, &foundhashdata);
1459 }
1460
1461 /* also loop through implied bounds to look for products */
1462 for( i = 0; i < SCIPgetNBinVars(scip); ++i )
1463 {
1464 /* first choose the x variable: it can be any binary variable in the problem */
1465 vars_xwy[0] = SCIPgetVars(scip)[i];
1466
1467 assert(SCIPvarGetType(vars_xwy[0]) == SCIP_VARTYPE_BINARY);
1468
1469 /* consider both possible values of x */
1470 for( f = 0; f <= 1; ++f )
1471 {
1472 xfixing = f == 1;
1473
1474 /* go through implications of x */
1475 for( r1 = 0; r1 < SCIPvarGetNImpls(vars_xwy[0], xfixing); ++r1 )
1476 {
1477 /* w is the implication var */
1478 vars_xwy[1] = SCIPvarGetImplVars(vars_xwy[0], xfixing)[r1];
1479 assert(SCIPvarGetType(vars_xwy[1]) != SCIP_VARTYPE_BINARY);
1480
1481 /* write the implication as a big-M constraint */
1482 implBndToBigM(scip, vars_xwy, 0, 1, SCIPvarGetImplTypes(vars_xwy[0], xfixing)[r1], xfixing,
1483 SCIPvarGetImplBounds(vars_xwy[0], xfixing)[r1], coefs1, &side1);
1484 sidetype1 = SCIPvarGetImplTypes(vars_xwy[0], xfixing)[r1] == SCIP_BOUNDTYPE_LOWER ?
1485 SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT;
1486
1487 /* if the global bound is equal to the implied bound, there is nothing to do */
1488 if( SCIPisZero(scip, coefs1[0]) )
1489 continue;
1490
1491 SCIPdebugMsg(scip, "Implication %s == %u => %s %s %g\n", SCIPvarGetName(vars_xwy[0]), xfixing,
1492 SCIPvarGetName(vars_xwy[1]), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=",
1493 SCIPvarGetImplBounds(vars_xwy[0], xfixing)[r1]);
1494 SCIPdebugMsg(scip, "Written as big-M: %g%s + %s %s %g\n", coefs1[0], SCIPvarGetName(vars_xwy[0]),
1495 SCIPvarGetName(vars_xwy[1]), sidetype1 == SCIP_SIDETYPE_LEFT ? ">=" : "<=", side1);
1496
1497 /* the second relation is in w and y (y could be anything, but must be in relation with w) */
1498
1499 /* x does not participate in the second relation, so we immediately set its coefficient to 0.0 */
1500 coefs2[0] = 0.0;
1501
1502 SCIPdebugMsg(scip, "Implic of x = <%s> + implied lb on w = <%s>:\n", SCIPvarGetName(vars_xwy[0]), SCIPvarGetName(vars_xwy[1]));
1503
1504 /* use implied lower bounds on w: w >= b*y + d */
1505 for( r2 = 0; r2 < SCIPvarGetNVlbs(vars_xwy[1]); ++r2 )
1506 {
1507 vars_xwy[2] = SCIPvarGetVlbVars(vars_xwy[1])[r2];
1508 if( vars_xwy[2] == vars_xwy[0] )
1509 continue;
1510
1511 coefs2[1] = 1.0;
1512 coefs2[2] = -SCIPvarGetVlbCoefs(vars_xwy[1])[r2];
1513
1514 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1515 SCIPvarGetVlbConstants(vars_xwy[1])[r2], sidetype1, SCIP_SIDETYPE_LEFT, varmap, xfixing) );
1516 }
1517
1518 SCIPdebugMsg(scip, "Implic of x = <%s> + implied ub on w = <%s>:\n", SCIPvarGetName(vars_xwy[0]), SCIPvarGetName(vars_xwy[1]));
1519
1520 /* use implied upper bounds on w: w <= b*y + d */
1521 for( r2 = 0; r2 < SCIPvarGetNVubs(vars_xwy[1]); ++r2 )
1522 {
1523 vars_xwy[2] = SCIPvarGetVubVars(vars_xwy[1])[r2];
1524 if( vars_xwy[2] == vars_xwy[0] )
1525 continue;
1526
1527 coefs2[1] = 1.0;
1528 coefs2[2] = -SCIPvarGetVubCoefs(vars_xwy[1])[r2];
1529
1530 SCIP_CALL( extractProducts(scip, sepadata, vars_xwy, coefs1, coefs2, side1,
1531 SCIPvarGetVubConstants(vars_xwy[1])[r2], sidetype1, SCIP_SIDETYPE_RIGHT, varmap, xfixing) );
1532 }
1533
1534 /* use unconditional relations containing w */
1535 relatedvars = getAdjacentVars(vars_in_2rels, vars_xwy[1], &nrelatedvars);
1536 if( relatedvars == NULL )
1537 continue;
1538
1539 for( r2 = 0; r2 < nrelatedvars; ++r2 )
1540 {
1541 vars_xwy[2] = relatedvars[r2];
1542 SCIPdebugMsg(scip, "Implied bound + unconditional with w and y:\n");
1543 SCIP_CALL( detectProductsUnconditional(scip, sepadata, prob_rows, row_list, hashtable2, coefs1,
1544 vars_xwy, side1, sidetype1, 1, 2, varmap, xfixing) );
1545 }
1546 }
1547 }
1548 }
1549
1550 /* free memory */
1551 clearVarAdjacency(scip, vars_in_2rels);
1552 SCIPhashmapFree(&vars_in_2rels);
1553
1554 SCIPdebugMsg(scip, "Unconditional relations table:\n");
1555 for( i = 0; i < SCIPhashtableGetNEntries(hashtable2); ++i )
1556 {
1557 foundhashdata = (HASHDATA*)SCIPhashtableGetEntry(hashtable2, i);
1558 if( foundhashdata == NULL )
1559 continue;
1560
1561 SCIPdebugMsg(scip, "(%s, %s): ", SCIPvarGetName(foundhashdata->vars[0]),
1562 SCIPvarGetName(foundhashdata->vars[1]));
1563
1564 SCIPfreeBuffer(scip, &foundhashdata);
1565 }
1566
1567 SCIPfreeBufferArray(scip, &row_list);
1568
1569 SCIPhashtableFree(&hashtable2);
1570 SCIPhashtableFree(&hashtable3);
1571
1572 SCIPfreeBufferArray(scip, &prob_rows);
1573
1574 return SCIP_OKAY;
1575 }
1576
1577 /** helper method to create separation data */
1578 static
1579 SCIP_RETCODE createSepaData(
1580 SCIP* scip, /**< SCIP data structure */
1581 SCIP_SEPADATA* sepadata /**< separation data */
1582 )
1583 {
1584 SCIP_HASHMAP* varmap;
1585 int i;
1586 SCIP_CONSNONLINEAR_BILINTERM* bilinterms;
1587 int varmapsize;
1588 int nvars;
1589
1590 assert(sepadata != NULL);
1591
1592 /* initialize some fields of sepadata */
1593 sepadata->varssorted = NULL;
1594 sepadata->varpriorities = NULL;
1595 sepadata->bilinvardatamap = NULL;
1596 sepadata->eqauxexpr = NULL;
1597 sepadata->nbilinvars = 0;
1598 sepadata->sbilinvars = 0;
1599
1600 /* get total number of bilinear terms */
1601 sepadata->nbilinterms = SCIPgetNBilinTermsNonlinear(sepadata->conshdlr);
1602
1603 /* skip if there are no bilinear terms and implicit product detection is off */
1604 if( sepadata->nbilinterms == 0 && !sepadata->detecthidden )
1605 return SCIP_OKAY;
1606
1607 /* the number of variables participating in bilinear products cannot exceed twice the number of bilinear terms;
1608 * however, if we detect hidden products, the number of terms is yet unknown, so use the number of variables
1609 */
1610 nvars = SCIPgetNVars(scip);
1611 varmapsize = sepadata->detecthidden ? nvars : MIN(nvars, sepadata->nbilinterms * 2);
1612
1613 /* create variable map */
1614 SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(scip), varmapsize) );
1615
1616 /* get all bilinear terms from the nonlinear constraint handler */
1617 bilinterms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
1618
1619 /* store the information of all variables that appear bilinearly */
1620 for( i = 0; i < sepadata->nbilinterms; ++i )
1621 {
1622 assert(bilinterms[i].x != NULL);
1623 assert(bilinterms[i].y != NULL);
1624 assert(bilinterms[i].nlockspos + bilinterms[i].nlocksneg > 0);
1625
1626 /* skip bilinear term if it does not have an auxiliary variable */
1627 if( bilinterms[i].aux.var == NULL )
1628 continue;
1629
1630 /* if only original variables should be used, skip products that contain at least one auxiliary variable */
1631 if( sepadata->onlyoriginal && (SCIPvarIsRelaxationOnly(bilinterms[i].x) ||
1632 SCIPvarIsRelaxationOnly(bilinterms[i].y)) )
1633 continue;
1634
1635 SCIP_CALL( addProductVars(scip, sepadata, bilinterms[i].x, bilinterms[i].y, varmap,
1636 bilinterms[i].nlockspos + bilinterms[i].nlocksneg) );
1637 }
1638
1639 if( sepadata->detecthidden )
1640 {
1641 int oldnterms = sepadata->nbilinterms;
1642
1643 SCIP_CALL( detectHiddenProducts(scip, sepadata, varmap) );
1644
1645 /* update nbilinterms and bilinterms, as detectHiddenProducts might have found new terms */
1646 sepadata->nbilinterms = SCIPgetNBilinTermsNonlinear(sepadata->conshdlr);
1647 bilinterms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
1648
1649 if( sepadata->nbilinterms > oldnterms )
1650 {
1651 SCIPstatisticMessage(" Number of hidden products: %d\n", sepadata->nbilinterms - oldnterms);
1652 }
1653 }
1654
1655 SCIPhashmapFree(&varmap);
1656
1657 if( sepadata->nbilinterms == 0 )
1658 {
1659 return SCIP_OKAY;
1660 }
1661
1662 /* mark positions of aux.exprs that must be equal to the product */
1663 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &sepadata->eqauxexpr, sepadata->nbilinterms) );
1664
1665 for( i = 0; i < sepadata->nbilinterms; ++i )
1666 {
1667 int j;
1668
1669 sepadata->eqauxexpr[i] = -1;
1670 for( j = 0; j < bilinterms[i].nauxexprs; ++j )
1671 {
1672 assert(bilinterms[i].aux.exprs[j] != NULL);
1673
1674 if( bilinterms[i].aux.exprs[j]->underestimate && bilinterms[i].aux.exprs[j]->overestimate )
1675 {
1676 sepadata->eqauxexpr[i] = j;
1677 break;
1678 }
1679 }
1680 }
1681
1682 /* find maxnumber of variables that occur most often and sort them by number of occurrences
1683 * (same as normal sort, except that entries at positions maxusedvars..nbilinvars may be unsorted at end)
1684 */
1685 SCIPselectDownIntPtr(sepadata->varpriorities, (void**) sepadata->varssorted, MIN(sepadata->maxusedvars,sepadata->nbilinvars-1),
1686 sepadata->nbilinvars);
1687
1688 /* capture all variables */
1689 for( i = 0; i < sepadata->nbilinvars; ++i )
1690 {
1691 assert(sepadata->varssorted[i] != NULL);
1692 SCIP_CALL( SCIPcaptureVar(scip, sepadata->varssorted[i]) );
1693 }
1694
1695 /* mark that separation data has been created */
1696 sepadata->iscreated = TRUE;
1697 sepadata->isinitialround = TRUE;
1698
1699 if( SCIPgetNBilinTermsNonlinear(sepadata->conshdlr) > 0 )
1700 SCIPstatisticMessage(" Found bilinear terms\n");
1701 else
1702 SCIPstatisticMessage(" No bilinear terms\n");
1703
1704 return SCIP_OKAY;
1705 }
1706
1707 /** get the positions of the most violated auxiliary under- and overestimators for each product
1708 *
1709 * -1 means no relation with given product is violated
1710 */
1711 static
1712 void getBestEstimators(
1713 SCIP* scip, /**< SCIP data structure */
1714 SCIP_SEPADATA* sepadata, /**< separator data */
1715 SCIP_SOL* sol, /**< solution at which to evaluate the expressions */
1716 int* bestunderestimators,/**< array of indices of best underestimators for each term */
1717 int* bestoverestimators /**< array of indices of best overestimators for each term */
1718 )
1719 {
1720 SCIP_Real prodval;
1721 SCIP_Real auxval;
1722 SCIP_Real prodviol;
1723 SCIP_Real viol_below;
1724 SCIP_Real viol_above;
1725 int i;
1726 int j;
1727 SCIP_CONSNONLINEAR_BILINTERM* terms;
1728
1729 assert(bestunderestimators != NULL);
1730 assert(bestoverestimators != NULL);
1731
1732 terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
1733
1734 for( j = 0; j < SCIPgetNBilinTermsNonlinear(sepadata->conshdlr); ++j )
1735 {
1736 viol_below = 0.0;
1737 viol_above = 0.0;
1738
1739 /* evaluate the product expression */
1740 prodval = SCIPgetSolVal(scip, sol, terms[j].x) * SCIPgetSolVal(scip, sol, terms[j].y);
1741
1742 bestunderestimators[j] = -1;
1743 bestoverestimators[j] = -1;
1744
1745 /* if there are any auxexprs, look there */
1746 for( i = 0; i < terms[j].nauxexprs; ++i )
1747 {
1748 auxval = SCIPevalBilinAuxExprNonlinear(scip, terms[j].x, terms[j].y, terms[j].aux.exprs[i], sol);
1749 prodviol = auxval - prodval;
1750
1751 if( terms[j].aux.exprs[i]->underestimate && SCIPisFeasGT(scip, auxval, prodval) && prodviol > viol_below )
1752 {
1753 viol_below = prodviol;
1754 bestunderestimators[j] = i;
1755 }
1756 if( terms[j].aux.exprs[i]->overestimate && SCIPisFeasGT(scip, prodval, auxval) && -prodviol > viol_above )
1757 {
1758 viol_above = -prodviol;
1759 bestoverestimators[j] = i;
1760 }
1761 }
1762
1763 /* if the term has a plain auxvar, it will be treated differently - do nothing here */
1764 }
1765 }
1766
1767 /** tests if a row contains too many unknown bilinear terms w.r.t. the parameters */
1768 static
1769 SCIP_RETCODE isAcceptableRow(
1770 SCIP_SEPADATA* sepadata, /**< separation data */
1771 SCIP_ROW* row, /**< the row to be tested */
1772 SCIP_VAR* var, /**< the variable that is to be multiplied with row */
1773 int* currentnunknown, /**< buffer to store number of unknown terms in current row if acceptable */
1774 SCIP_Bool* acceptable /**< buffer to store the result */
1775 )
1776 {
1777 int i;
1778 int idx;
1779 SCIP_CONSNONLINEAR_BILINTERM* terms;
1780
1781 assert(row != NULL);
1782 assert(var != NULL);
1783
1784 *currentnunknown = 0;
1785 terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
1786
1787 for( i = 0; (i < SCIProwGetNNonz(row)) && (sepadata->maxunknownterms < 0 || *currentnunknown <= sepadata->maxunknownterms); ++i )
1788 {
1789 idx = SCIPgetBilinTermIdxNonlinear(sepadata->conshdlr, var, SCIPcolGetVar(SCIProwGetCols(row)[i]));
1790
1791 /* if the product hasn't been found, no auxiliary expressions for it are known */
1792 if( idx < 0 )
1793 {
1794 ++(*currentnunknown);
1795 continue;
1796 }
1797
1798 /* known terms are only those that have an aux.var or equality estimators */
1799 if( sepadata->eqauxexpr[idx] == -1 && !(terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL) )
1800 {
1801 ++(*currentnunknown);
1802 }
1803 }
1804
1805 *acceptable = sepadata->maxunknownterms < 0 || *currentnunknown <= sepadata->maxunknownterms;
1806
1807 return SCIP_OKAY;
1808 }
1809
1810 /** adds coefficients and constant of an auxiliary expression
1811 *
1812 * the variables the pointers are pointing to must already be initialized
1813 */
1814 static
1815 void addAuxexprCoefs(
1816 SCIP_VAR* var1, /**< first product variable */
1817 SCIP_VAR* var2, /**< second product variable */
1818 SCIP_CONSNONLINEAR_AUXEXPR* auxexpr, /**< auxiliary expression to be added */
1819 SCIP_Real coef, /**< coefficient of the auxiliary expression */
1820 SCIP_Real* coefaux, /**< pointer to add the coefficient of the auxiliary variable */
1821 SCIP_Real* coef1, /**< pointer to add the coefficient of the first variable */
1822 SCIP_Real* coef2, /**< pointer to add the coefficient of the second variable */
1823 SCIP_Real* cst /**< pointer to add the constant */
1824 )
1825 {
1826 assert(auxexpr != NULL);
1827 assert(auxexpr->auxvar != NULL);
1828 assert(coefaux != NULL);
1829 assert(coef1 != NULL);
1830 assert(coef2 != NULL);
1831 assert(cst != NULL);
1832
1833 *coefaux += auxexpr->coefs[0] * coef;
1834
1835 /* in auxexpr, x goes before y and has the smaller index,
1836 * so compare vars to figure out which one is x and which is y
1837 */
1838 if( SCIPvarCompare(var1, var2) < 1 )
1839 {
1840 *coef1 += auxexpr->coefs[1] * coef;
1841 *coef2 += auxexpr->coefs[2] * coef;
1842 }
1843 else
1844 {
1845 *coef1 += auxexpr->coefs[2] * coef;
1846 *coef2 += auxexpr->coefs[1] * coef;
1847 }
1848 *cst += coef * auxexpr->cst;
1849 }
1850
1851 /** add a linear term `coef`*`colvar` multiplied by a bound factor (var - lb(var)) or (ub(var) - var)
1852 *
1853 * adds the linear term with `colvar` to `cut` and updates `coefvar` and `cst`
1854 */
1855 static
1856 SCIP_RETCODE addRltTerm(
1857 SCIP* scip, /**< SCIP data structure */
1858 SCIP_SEPADATA* sepadata, /**< separator data */
1859 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
1860 int* bestunderest, /**< positions of most violated underestimators for each product term */
1861 int* bestoverest, /**< positions of most violated overestimators for each product term */
1862 SCIP_ROW* cut, /**< cut to which the term is to be added */
1863 SCIP_VAR* var, /**< multiplier variable */
1864 SCIP_VAR* colvar, /**< row variable to be multiplied */
1865 SCIP_Real coef, /**< coefficient of the bilinear term */
1866 SCIP_Bool uselb, /**< whether we multiply with (var - lb) or (ub - var) */
1867 SCIP_Bool uselhs, /**< whether to create a cut for the lhs or rhs */
1868 SCIP_Bool local, /**< whether local or global cuts should be computed */
1869 SCIP_Bool computeEqCut, /**< whether conditions are fulfilled to compute equality cuts */
1870 SCIP_Real* coefvar, /**< coefficient of var */
1871 SCIP_Real* cst, /**< buffer to store the constant part of the cut */
1872 SCIP_Bool* success /**< buffer to store whether cut was updated successfully */
1873 )
1874 {
1875 SCIP_Real lbvar;
1876 SCIP_Real ubvar;
1877 SCIP_Real refpointvar;
1878 SCIP_Real signfactor;
1879 SCIP_Real boundfactor;
1880 SCIP_Real coefauxvar;
1881 SCIP_Real coefcolvar;
1882 SCIP_Real coefterm;
1883 int auxpos;
1884 int idx;
1885 SCIP_CONSNONLINEAR_BILINTERM* terms;
1886 SCIP_VAR* auxvar;
1887
1888 terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
1889
1890 if( computeEqCut )
1891 {
1892 lbvar = 0.0;
1893 ubvar = 0.0;
1894 }
1895 else
1896 {
1897 lbvar = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
1898 ubvar = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
1899 }
1900
1901 refpointvar = MAX(lbvar, MIN(ubvar, SCIPgetSolVal(scip, sol, var))); /*lint !e666*/
1902
1903 signfactor = (uselb ? 1.0 : -1.0);
1904 boundfactor = (uselb ? -lbvar : ubvar);
1905
1906 coefterm = coef * signfactor; /* coefficient of the bilinear term */
1907 coefcolvar = coef * boundfactor; /* coefficient of the linear term */
1908 coefauxvar = 0.0; /* coefficient of the auxiliary variable corresponding to the bilinear term */
1909 auxvar = NULL;
1910
1911 assert(!SCIPisInfinity(scip, REALABS(coefterm)));
1912
1913 /* first, add the linearisation of the bilinear term */
1914
1915 idx = SCIPgetBilinTermIdxNonlinear(sepadata->conshdlr, var, colvar);
1916 auxpos = -1;
1917
1918 /* for an implicit term, get the position of the best estimator */
1919 if( idx >= 0 && terms[idx].nauxexprs > 0 )
1920 {
1921 if( computeEqCut )
1922 {
1923 /* use an equality auxiliary expression (which should exist for computeEqCut to be TRUE) */
1924 assert(sepadata->eqauxexpr[idx] >= 0);
1925 auxpos = sepadata->eqauxexpr[idx];
1926 }
1927 else if( (uselhs && coefterm > 0.0) || (!uselhs && coefterm < 0.0) )
1928 {
1929 /* use an overestimator */
1930 auxpos = bestoverest[idx];
1931 }
1932 else
1933 {
1934 /* use an underestimator */
1935 auxpos = bestunderest[idx];
1936 }
1937 }
1938
1939 /* if the term is implicit and a suitable auxiliary expression for var*colvar exists, add the coefficients
1940 * of the auxiliary expression for coefterm*var*colvar to coefauxvar, coefcolvar, coefvar and cst
1941 */
1942 if( auxpos >= 0 )
1943 {
1944 SCIPdebugMsg(scip, "auxiliary expression for <%s> and <%s> found, will be added to cut:\n",
1945 SCIPvarGetName(colvar), SCIPvarGetName(var));
1946 addAuxexprCoefs(var, colvar, terms[idx].aux.exprs[auxpos], coefterm, &coefauxvar, coefvar, &coefcolvar, cst);
1947 auxvar = terms[idx].aux.exprs[auxpos]->auxvar;
1948 }
1949 /* for an existing term, use the auxvar if there is one */
1950 else if( idx >= 0 && terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
1951 {
1952 SCIPdebugMsg(scip, "auxvar for <%s> and <%s> found, will be added to cut:\n",
1953 SCIPvarGetName(colvar), SCIPvarGetName(var));
1954 coefauxvar += coefterm;
1955 auxvar = terms[idx].aux.var;
1956 }
1957
1958 /* otherwise, use clique information or the McCormick estimator in place of the bilinear term */
1959 else if( colvar != var )
1960 {
1961 SCIP_Bool found_clique = FALSE;
1962 SCIP_Real lbcolvar = local ? SCIPvarGetLbLocal(colvar) : SCIPvarGetLbGlobal(colvar);
1963 SCIP_Real ubcolvar = local ? SCIPvarGetUbLocal(colvar) : SCIPvarGetUbGlobal(colvar);
1964 SCIP_Real refpointcolvar = MAX(lbcolvar, MIN(ubcolvar, SCIPgetSolVal(scip, sol, colvar))); /*lint !e666*/
1965
1966 assert(!computeEqCut);
1967
1968 if( REALABS(lbcolvar) > MAXVARBOUND || REALABS(ubcolvar) > MAXVARBOUND )
1969 {
1970 *success = FALSE;
1971 return SCIP_OKAY;
1972 }
1973
1974 SCIPdebugMsg(scip, "auxvar for <%s> and <%s> not found, will linearize the product\n", SCIPvarGetName(colvar), SCIPvarGetName(var));
1975
1976 /* if both variables are binary, check if they are contained together in some clique */
1977 if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY && SCIPvarGetType(colvar) == SCIP_VARTYPE_BINARY )
1978 {
1979 int c;
1980 SCIP_CLIQUE** varcliques;
1981
1982 varcliques = SCIPvarGetCliques(var, TRUE);
1983
1984 /* look through cliques containing var */
1985 for( c = 0; c < SCIPvarGetNCliques(var, TRUE); ++c )
1986 {
1987 if( SCIPcliqueHasVar(varcliques[c], colvar, TRUE) ) /* var + colvar <= 1 => var*colvar = 0 */
1988 {
1989 /* product is zero, add nothing */
1990 found_clique = TRUE;
1991 break;
1992 }
1993
1994 if( SCIPcliqueHasVar(varcliques[c], colvar, FALSE) ) /* var + (1-colvar) <= 1 => var*colvar = var */
1995 {
1996 *coefvar += coefterm;
1997 found_clique = TRUE;
1998 break;
1999 }
2000 }
2001
2002 if( !found_clique )
2003 {
2004 varcliques = SCIPvarGetCliques(var, FALSE);
2005
2006 /* look through cliques containing complement of var */
2007 for( c = 0; c < SCIPvarGetNCliques(var, FALSE); ++c )
2008 {
2009 if( SCIPcliqueHasVar(varcliques[c], colvar, TRUE) ) /* (1-var) + colvar <= 1 => var*colvar = colvar */
2010 {
2011 coefcolvar += coefterm;
2012 found_clique = TRUE;
2013 break;
2014 }
2015
2016 if( SCIPcliqueHasVar(varcliques[c], colvar, FALSE) ) /* (1-var) + (1-colvar) <= 1 => var*colvar = var + colvar - 1 */
2017 {
2018 *coefvar += coefterm;
2019 coefcolvar += coefterm;
2020 *cst -= coefterm;
2021 found_clique = TRUE;
2022 break;
2023 }
2024 }
2025 }
2026 }
2027
2028 if( !found_clique )
2029 {
2030 SCIPdebugMsg(scip, "clique for <%s> and <%s> not found or at least one of them is not binary, will use McCormick\n", SCIPvarGetName(colvar), SCIPvarGetName(var));
2031 SCIPaddBilinMcCormick(scip, coefterm, lbvar, ubvar, refpointvar, lbcolvar,
2032 ubcolvar, refpointcolvar, uselhs, coefvar, &coefcolvar, cst, success);
2033 if( !*success )
2034 return SCIP_OKAY;
2035 }
2036 }
2037
2038 /* or, if it's a quadratic term, use a secant for overestimation and a gradient for underestimation */
2039 else
2040 {
2041 SCIPdebugMsg(scip, "auxvar for <%s>^2 not found, will use gradient and secant estimators\n", SCIPvarGetName(colvar));
2042
2043 assert(!computeEqCut);
2044
2045 /* for a binary var, var^2 = var */
2046 if( SCIPvarGetType(var) == SCIP_VARTYPE_BINARY )
2047 {
2048 *coefvar += coefterm;
2049 }
2050 else
2051 {
2052 /* depending on over-/underestimation and the sign of the column variable, compute secant or tangent */
2053 if( (uselhs && coefterm > 0.0) || (!uselhs && coefterm < 0.0) )
2054 SCIPaddSquareSecant(scip, coefterm, lbvar, ubvar, coefvar, cst, success);
2055 else
2056 SCIPaddSquareLinearization(scip, coefterm, refpointvar, SCIPvarIsIntegral(var), coefvar, cst, success);
2057
2058 if( !*success )
2059 return SCIP_OKAY;
2060 }
2061 }
2062
2063 /* add the auxiliary variable if its coefficient is nonzero */
2064 if( !SCIPisZero(scip, coefauxvar) )
2065 {
2066 assert(auxvar != NULL);
2067 SCIP_CALL( SCIPaddVarToRow(scip, cut, auxvar, coefauxvar) );
2068 }
2069
2070 /* we are done with the product linearisation, now add the term which comes from multiplying
2071 * coef*colvar by the constant part of the bound factor
2072 */
2073
2074 if( colvar != var )
2075 {
2076 assert(!SCIPisInfinity(scip, REALABS(coefcolvar)));
2077 SCIP_CALL( SCIPaddVarToRow(scip, cut, colvar, coefcolvar) );
2078 }
2079 else
2080 *coefvar += coefcolvar;
2081
2082 return SCIP_OKAY;
2083 }
2084
2085 /** creates the RLT cut formed by multiplying a given row with (x - lb) or (ub - x)
2086 *
2087 * In detail:
2088 * - The row is multiplied either with (x - lb(x)) or with (ub(x) - x), depending on parameter `uselb`, or by x if
2089 * this is an equality cut
2090 * - The (inequality) cut is computed either for lhs or rhs, depending on parameter `uselhs`.
2091 * - Terms for which no auxiliary variable and no clique relation exists are replaced by either McCormick, secants,
2092 * or gradient linearization cuts.
2093 */
2094 static
2095 SCIP_RETCODE computeRltCut(
2096 SCIP* scip, /**< SCIP data structure */
2097 SCIP_SEPA* sepa, /**< separator */
2098 SCIP_SEPADATA* sepadata, /**< separation data */
2099 SCIP_ROW** cut, /**< buffer to store the cut */
2100 SCIP_ROW* row, /**< the row that is used for the rlt cut (NULL if using projected row) */
2101 RLT_SIMPLEROW* projrow, /**< projected row that is used for the rlt cut (NULL if using row) */
2102 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
2103 int* bestunderest, /**< positions of most violated underestimators for each product term */
2104 int* bestoverest, /**< positions of most violated overestimators for each product term */
2105 SCIP_VAR* var, /**< the variable that is used for the rlt cuts */
2106 SCIP_Bool* success, /**< buffer to store whether cut was created successfully */
2107 SCIP_Bool uselb, /**< whether we multiply with (var - lb) or (ub - var) */
2108 SCIP_Bool uselhs, /**< whether to create a cut for the lhs or rhs */
2109 SCIP_Bool local, /**< whether local or global cuts should be computed */
2110 SCIP_Bool computeEqCut, /**< whether conditions are fulfilled to compute equality cuts */
2111 SCIP_Bool useprojrow /**< whether to use projected row instead of normal row */
2112 )
2113 { /*lint --e{413}*/
2114 SCIP_Real signfactor;
2115 SCIP_Real boundfactor;
2116 SCIP_Real lbvar;
2117 SCIP_Real ubvar;
2118 SCIP_Real coefvar;
2119 SCIP_Real consside;
2120 SCIP_Real finalside;
2121 SCIP_Real cstterm;
2122 SCIP_Real lhs;
2123 SCIP_Real rhs;
2124 SCIP_Real rowcst;
2125 int i;
2126 const char* rowname;
2127 char cutname[SCIP_MAXSTRLEN];
2128
2129 assert(sepadata != NULL);
2130 assert(cut != NULL);
2131 assert(useprojrow || row != NULL);
2132 assert(!useprojrow || projrow != NULL);
2133 assert(var != NULL);
2134 assert(success != NULL);
2135
2136 lhs = useprojrow ? projrow->lhs : SCIProwGetLhs(row);
2137 rhs = useprojrow ? projrow->rhs : SCIProwGetRhs(row);
2138 rowname = useprojrow ? projrow->name : SCIProwGetName(row);
2139 rowcst = useprojrow ? projrow ->cst : SCIProwGetConstant(row);
2140
2141 assert(!computeEqCut || SCIPisEQ(scip, lhs, rhs));
2142
2143 *cut = NULL;
2144
2145 /* get data for given variable */
2146 if( computeEqCut )
2147 {
2148 lbvar = 0.0;
2149 ubvar = 0.0;
2150 }
2151 else
2152 {
2153 lbvar = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
2154 ubvar = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
2155 }
2156
2157 /* get row side */
2158 consside = uselhs ? lhs : rhs;
2159
2160 /* if the bounds are too large or the respective side is infinity, skip this cut */
2161 if( (uselb && REALABS(lbvar) > MAXVARBOUND) || (!uselb && REALABS(ubvar) > MAXVARBOUND)
2162 || SCIPisInfinity(scip, REALABS(consside)) )
2163 {
2164 SCIPdebugMsg(scip, "cut generation for %srow <%s>, %s, and variable <%s> with its %s %g not possible\n",
2165 useprojrow ? "projected " : "", rowname, uselhs ? "lhs" : "rhs", SCIPvarGetName(var),
2166 uselb ? "lower bound" : "upper bound", uselb ? lbvar : ubvar);
2167
2168 if( REALABS(lbvar) > MAXVARBOUND )
2169 SCIPdebugMsg(scip, " because of lower bound\n");
2170 if( REALABS(ubvar) > MAXVARBOUND )
2171 SCIPdebugMsg(scip, " because of upper bound\n");
2172 if( SCIPisInfinity(scip, REALABS(consside)) )
2173 SCIPdebugMsg(scip, " because of side %g\n", consside);
2174
2175 *success = FALSE;
2176 return SCIP_OKAY;
2177 }
2178
2179 /* initialize some factors needed for computation */
2180 coefvar = 0.0;
2181 cstterm = 0.0;
2182 signfactor = (uselb ? 1.0 : -1.0);
2183 boundfactor = (uselb ? -lbvar : ubvar);
2184 *success = TRUE;
2185
2186 /* create an empty row which we then fill with variables step by step */
2187 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "rlt_%scut_%s_%s_%s_%s_%d", useprojrow ? "proj" : "", rowname,
2188 uselhs ? "lhs" : "rhs", SCIPvarGetName(var), uselb ? "lb" : "ub", SCIPgetNLPs(scip));
2189 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, cut, sepa, cutname, -SCIPinfinity(scip), SCIPinfinity(scip),
2190 SCIPgetDepth(scip) > 0 && local, FALSE, FALSE) ); /* TODO SCIPgetDepth() should be replaced by depth that is passed on to the SEPAEXEC calls (?) */
2191
2192 SCIP_CALL( SCIPcacheRowExtensions(scip, *cut) );
2193
2194 /* iterate over all variables in the row and add the corresponding terms coef*colvar*(bound factor) to the cuts */
2195 for( i = 0; i < (useprojrow ? projrow->nnonz : SCIProwGetNNonz(row)); ++i )
2196 {
2197 SCIP_VAR* colvar;
2198
2199 colvar = useprojrow ? projrow->vars[i] : SCIPcolGetVar(SCIProwGetCols(row)[i]);
2200 SCIP_CALL( addRltTerm(scip, sepadata, sol, bestunderest, bestoverest, *cut, var, colvar,
2201 useprojrow ? projrow->coefs[i] : SCIProwGetVals(row)[i], uselb, uselhs, local, computeEqCut,
2202 &coefvar, &cstterm, success) );
2203 }
2204
2205 if( REALABS(cstterm) > MAXVARBOUND )
2206 {
2207 *success = FALSE;
2208 return SCIP_OKAY;
2209 }
2210
2211 /* multiply (x-lb) or (ub -x) with the lhs and rhs of the row */
2212 coefvar += signfactor * (rowcst - consside);
2213 finalside = boundfactor * (consside - rowcst) - cstterm;
2214
2215 assert(!SCIPisInfinity(scip, REALABS(coefvar)));
2216 assert(!SCIPisInfinity(scip, REALABS(finalside)));
2217
2218 /* set the coefficient of var and update the side */
2219 SCIP_CALL( SCIPaddVarToRow(scip, *cut, var, coefvar) );
2220 SCIP_CALL( SCIPflushRowExtensions(scip, *cut) );
2221 if( uselhs || computeEqCut )
2222 {
2223 SCIP_CALL( SCIPchgRowLhs(scip, *cut, finalside) );
2224 }
2225 if( !uselhs || computeEqCut )
2226 {
2227 SCIP_CALL( SCIPchgRowRhs(scip, *cut, finalside) );
2228 }
2229
2230 SCIPdebugMsg(scip, "%scut was generated successfully:\n", useprojrow ? "projected " : "");
2231 #ifdef SCIP_DEBUG
2232 SCIP_CALL( SCIPprintRow(scip, *cut, NULL) );
2233 #endif
2234
2235 return SCIP_OKAY;
2236 }
2237
2238 /** store a row projected by fixing all variables that are at bound at sol; the result is a simplified row */
2239 static
2240 SCIP_RETCODE createProjRow(
2241 SCIP* scip, /**< SCIP data structure */
2242 RLT_SIMPLEROW* simplerow, /**< pointer to the simplified row */
2243 SCIP_ROW* row, /**< row to be projected */
2244 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
2245 SCIP_Bool local /**< whether local bounds should be checked */
2246 )
2247 {
2248 int i;
2249 SCIP_VAR* var;
2250 SCIP_Real val;
2251 SCIP_Real vlb;
2252 SCIP_Real vub;
2253
2254 assert(simplerow != NULL);
2255
2256 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(simplerow->name), SCIProwGetName(row),
2257 strlen(SCIProwGetName(row))+1) ); /*lint !e666*/
2258 simplerow->nnonz = 0;
2259 simplerow->size = 0;
2260 simplerow->vars = NULL;
2261 simplerow->coefs = NULL;
2262 simplerow->lhs = SCIProwGetLhs(row);
2263 simplerow->rhs = SCIProwGetRhs(row);
2264 simplerow->cst = SCIProwGetConstant(row);
2265
2266 for( i = 0; i < SCIProwGetNNonz(row); ++i )
2267 {
2268 var = SCIPcolGetVar(SCIProwGetCols(row)[i]);
2269 val = SCIPgetSolVal(scip, sol, var);
2270 vlb = local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
2271 vub = local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
2272 if( SCIPisFeasEQ(scip, vlb, val) || SCIPisFeasEQ(scip, vub, val) )
2273 {
2274 /* if we are projecting and the var is at bound, add var as a constant to simplerow */
2275 if( !SCIPisInfinity(scip, -simplerow->lhs) )
2276 simplerow->lhs -= SCIProwGetVals(row)[i]*val;
2277 if( !SCIPisInfinity(scip, simplerow->rhs) )
2278 simplerow->rhs -= SCIProwGetVals(row)[i]*val;
2279 }
2280 else
2281 {
2282 if( simplerow->nnonz + 1 > simplerow->size )
2283 {
2284 int newsize;
2285
2286 newsize = SCIPcalcMemGrowSize(scip, simplerow->nnonz + 1);
2287 SCIP_CALL( SCIPreallocBufferArray(scip, &simplerow->coefs, newsize) );
2288 SCIP_CALL( SCIPreallocBufferArray(scip, &simplerow->vars, newsize) );
2289 simplerow->size = newsize;
2290 }
2291
2292 /* add the term to simplerow */
2293 simplerow->vars[simplerow->nnonz] = var;
2294 simplerow->coefs[simplerow->nnonz] = SCIProwGetVals(row)[i];
2295 ++(simplerow->nnonz);
2296 }
2297 }
2298
2299 return SCIP_OKAY;
2300 }
2301
2302 /** free the projected row */
2303 static
2304 void freeProjRow(
2305 SCIP* scip, /**< SCIP data structure */
2306 RLT_SIMPLEROW* simplerow /**< simplified row to be freed */
2307 )
2308 {
2309 assert(simplerow != NULL);
2310
2311 if( simplerow->size > 0 )
2312 {
2313 assert(simplerow->vars != NULL);
2314 assert(simplerow->coefs != NULL);
2315
2316 SCIPfreeBufferArray(scip, &simplerow->vars);
2317 SCIPfreeBufferArray(scip, &simplerow->coefs);
2318 }
2319 SCIPfreeBlockMemoryArray(scip, &simplerow->name, strlen(simplerow->name)+1);
2320 }
2321
2322 /** creates the projected problem
2323 *
2324 * All variables that are at their bounds at the current solution are added
2325 * to left and/or right hand sides as constant values.
2326 */
2327 static
2328 SCIP_RETCODE createProjRows(
2329 SCIP* scip, /**< SCIP data structure */
2330 SCIP_ROW** rows, /**< problem rows */
2331 int nrows, /**< number of rows */
2332 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
2333 RLT_SIMPLEROW** projrows, /**< the projected rows to be filled */
2334 SCIP_Bool local, /**< are local cuts allowed? */
2335 SCIP_Bool* allcst /**< buffer to store whether all projected rows have only constants */
2336 )
2337 {
2338 int i;
2339
2340 assert(scip != NULL);
2341 assert(rows != NULL);
2342 assert(projrows != NULL);
2343 assert(allcst != NULL);
2344
2345 *allcst = TRUE;
2346 SCIP_CALL( SCIPallocBufferArray(scip, projrows, nrows) );
2347
2348 for( i = 0; i < nrows; ++i )
2349 {
2350 /* get a simplified and projected row */
2351 SCIP_CALL( createProjRow(scip, &(*projrows)[i], rows[i], sol, local) );
2352 if( (*projrows)[i].nnonz > 0 )
2353 *allcst = FALSE;
2354 }
2355
2356 return SCIP_OKAY;
2357 }
2358
2359 #ifdef SCIP_DEBUG
2360 /* prints the projected LP */
2361 static
2362 void printProjRows(
2363 SCIP* scip, /**< SCIP data structure */
2364 RLT_SIMPLEROW* projrows, /**< the projected rows */
2365 int nrows, /**< number of projected rows */
2366 FILE* file /**< output file (or NULL for standard output) */
2367 )
2368 {
2369 int i;
2370 int j;
2371
2372 assert(projrows != NULL);
2373
2374 for( i = 0; i < nrows; ++i )
2375 {
2376 SCIPinfoMessage(scip, file, "\nproj_row[%d]: ", i);
2377 if( !SCIPisInfinity(scip, -projrows[i].lhs) )
2378 SCIPinfoMessage(scip, file, "%.15g <= ", projrows[i].lhs);
2379 for( j = 0; j < projrows[i].nnonz; ++j )
2380 {
2381 if( j == 0 )
2382 {
2383 if( projrows[i].coefs[j] < 0 )
2384 SCIPinfoMessage(scip, file, "-");
2385 }
2386 else
2387 {
2388 if( projrows[i].coefs[j] < 0 )
2389 SCIPinfoMessage(scip, file, " - ");
2390 else
2391 SCIPinfoMessage(scip, file, " + ");
2392 }
2393
2394 if( projrows[i].coefs[j] != 1.0 )
2395 SCIPinfoMessage(scip, file, "%.15g*", REALABS(projrows[i].coefs[j]));
2396 SCIPinfoMessage(scip, file, "<%s>", SCIPvarGetName(projrows[i].vars[j]));
2397 }
2398 if( projrows[i].cst > 0 )
2399 SCIPinfoMessage(scip, file, " + %.15g", projrows[i].cst);
2400 else if( projrows[i].cst < 0 )
2401 SCIPinfoMessage(scip, file, " - %.15g", REALABS(projrows[i].cst));
2402
2403 if( !SCIPisInfinity(scip, projrows[i].rhs) )
2404 SCIPinfoMessage(scip, file, " <= %.15g", projrows[i].rhs);
2405 }
2406 SCIPinfoMessage(scip, file, "\n");
2407 }
2408 #endif
2409
2410 /** frees the projected rows */
2411 static
2412 void freeProjRows(
2413 SCIP* scip, /**< SCIP data structure */
2414 RLT_SIMPLEROW** projrows, /**< the projected LP */
2415 int nrows /**< number of rows in projrows */
2416 )
2417 {
2418 int i;
2419
2420 for( i = 0; i < nrows; ++i )
2421 freeProjRow(scip, &(*projrows)[i]);
2422
2423 SCIPfreeBufferArray(scip, projrows);
2424 }
2425
2426 /** mark a row for rlt cut selection
2427 *
2428 * depending on the sign of the coefficient and violation, set or update mark which cut is required:
2429 * - 1 - cuts for axy < aw case,
2430 * - 2 - cuts for axy > aw case,
2431 * - 3 - cuts for both cases
2432 */
2433 static
2434 void addRowMark(
2435 int ridx, /**< row index */
2436 SCIP_Real a, /**< coefficient of x in the row */
2437 SCIP_Bool violatedbelow, /**< whether the relation auxexpr <= xy is violated */
2438 SCIP_Bool violatedabove, /**< whether the relation xy <= auxexpr is violated */
2439 int* row_idcs, /**< sparse array with indices of marked rows */
2440 unsigned int* row_marks, /**< sparse array to store the marks */
2441 int* nmarked /**< number of marked rows */
2442 )
2443 {
2444 unsigned int newmark;
2445 int pos;
2446 SCIP_Bool exists;
2447
2448 assert(a != 0.0);
2449
2450 if( (a > 0.0 && violatedbelow) || (a < 0.0 && violatedabove) )
2451 newmark = 1; /* axy < aw case */
2452 else
2453 newmark = 2; /* axy > aw case */
2454
2455 /* find row idx in row_idcs */
2456 exists = SCIPsortedvecFindInt(row_idcs, ridx, *nmarked, &pos);
2457
2458 if( exists )
2459 {
2460 /* we found the row index: update the mark at pos */
2461 row_marks[pos] |= newmark;
2462 }
2463 else /* the given row index does not yet exist in row_idcs */
2464 {
2465 int i;
2466
2467 /* insert row index at the correct position */
2468 for( i = *nmarked; i > pos; --i )
2469 {
2470 row_idcs[i] = row_idcs[i-1];
2471 row_marks[i] = row_marks[i-1];
2472 }
2473 row_idcs[pos] = ridx;
2474 row_marks[pos] = newmark;
2475 (*nmarked)++;
2476 }
2477 }
2478
2479 /** mark all rows that should be multiplied by xj */
2480 static
2481 SCIP_RETCODE markRowsXj(
2482 SCIP* scip, /**< SCIP data structure */
2483 SCIP_SEPADATA* sepadata, /**< separator data */
2484 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
2485 SCIP_SOL* sol, /**< point to be separated (can be NULL) */
2486 int j, /**< index of the multiplier variable in sepadata */
2487 SCIP_Bool local, /**< are local cuts allowed? */
2488 SCIP_HASHMAP* row_to_pos, /**< hashmap linking row indices to positions in array */
2489 int* bestunderest, /**< positions of most violated underestimators for each product term */
2490 int* bestoverest, /**< positions of most violated overestimators for each product term */
2491 unsigned int* row_marks, /**< sparse array storing the row marks */
2492 int* row_idcs, /**< sparse array storing the marked row positions */
2493 int* nmarked /**< number of marked rows */
2494 )
2495 {
2496 int i;
2497 int idx;
2498 int ncolrows;
2499 int r;
2500 int ridx;
2501 SCIP_VAR* xi;
2502 SCIP_VAR* xj;
2503 SCIP_Real vlb;
2504 SCIP_Real vub;
2505 SCIP_Real vali;
2506 SCIP_Real valj;
2507 SCIP_Real a;
2508 SCIP_COL* coli;
2509 SCIP_Real* colvals;
2510 SCIP_ROW** colrows;
2511 SCIP_CONSNONLINEAR_BILINTERM* terms;
2512 SCIP_Bool violatedbelow;
2513 SCIP_Bool violatedabove;
2514 SCIP_VAR** bilinadjvars;
2515 int nbilinadjvars;
2516
2517 *nmarked = 0;
2518
2519 xj = sepadata->varssorted[j];
2520 assert(xj != NULL);
2521
2522 valj = SCIPgetSolVal(scip, sol, xj);
2523 vlb = local ? SCIPvarGetLbLocal(xj) : SCIPvarGetLbGlobal(xj);
2524 vub = local ? SCIPvarGetUbLocal(xj) : SCIPvarGetUbGlobal(xj);
2525
2526 if( sepadata->useprojection && (SCIPisFeasEQ(scip, vlb, valj) || SCIPisFeasEQ(scip, vub, valj)) )
2527 {
2528 /* we don't want to multiply by variables that are at bound */
2529 SCIPdebugMsg(scip, "Rejected multiplier <%s> in [%g,%g] because it is at bound (current value %g)\n", SCIPvarGetName(xj), vlb, vub, valj);
2530 return SCIP_OKAY;
2531 }
2532
2533 terms = SCIPgetBilinTermsNonlinear(conshdlr);
2534 bilinadjvars = getAdjacentVars(sepadata->bilinvardatamap, xj, &nbilinadjvars);
2535 assert(bilinadjvars != NULL);
2536
2537 /* for each var which appears in a bilinear product together with xj, mark rows */
2538 for( i = 0; i < nbilinadjvars; ++i )
2539 {
2540 xi = bilinadjvars[i];
2541
2542 if( SCIPvarGetStatus(xi) != SCIP_VARSTATUS_COLUMN )
2543 continue;
2544
2545 vali = SCIPgetSolVal(scip, sol, xi);
2546 vlb = local ? SCIPvarGetLbLocal(xi) : SCIPvarGetLbGlobal(xi);
2547 vub = local ? SCIPvarGetUbLocal(xi) : SCIPvarGetUbGlobal(xi);
2548
2549 /* if we use projection, we aren't interested in products with variables that are at bound */
2550 if( sepadata->useprojection && (SCIPisFeasEQ(scip, vlb, vali) || SCIPisFeasEQ(scip, vub, vali)) )
2551 continue;
2552
2553 /* get the index of the bilinear product */
2554 idx = SCIPgetBilinTermIdxNonlinear(conshdlr, xj, xi);
2555 assert(idx >= 0 && idx < SCIPgetNBilinTermsNonlinear(conshdlr));
2556
2557 /* skip implicit products if we don't want to add RLT cuts for them */
2558 if( !sepadata->hiddenrlt && !terms[idx].existing )
2559 continue;
2560
2561 /* use the most violated under- and overestimators for this product;
2562 * if equality cuts are computed, we might end up using a different auxiliary expression;
2563 * so this is an optimistic (i.e. taking the largest possible violation) estimation
2564 */
2565 if( bestunderest == NULL || bestunderest[idx] == -1 )
2566 { /* no violated implicit underestimation relations -> either use auxvar or set violatedbelow to FALSE */
2567 if( terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
2568 {
2569 assert(terms[idx].existing);
2570 violatedbelow = SCIPisFeasPositive(scip, SCIPgetSolVal(scip, sol, terms[idx].aux.var) - valj * vali);
2571 }
2572 else
2573 {
2574 assert(bestunderest != NULL);
2575 violatedbelow = FALSE;
2576 }
2577 }
2578 else
2579 {
2580 assert(bestunderest[idx] >= 0 && bestunderest[idx] < terms[idx].nauxexprs);
2581
2582 /* if we are here, the relation with the best underestimator must be violated */
2583 assert(SCIPisFeasPositive(scip, SCIPevalBilinAuxExprNonlinear(scip, terms[idx].x, terms[idx].y,
2584 terms[idx].aux.exprs[bestunderest[idx]], sol) - valj * vali));
2585 violatedbelow = TRUE;
2586 }
2587
2588 if( bestoverest == NULL || bestoverest[idx] == -1 )
2589 { /* no violated implicit overestimation relations -> either use auxvar or set violatedabove to FALSE */
2590 if( terms[idx].nauxexprs == 0 && terms[idx].aux.var != NULL )
2591 {
2592 assert(terms[idx].existing);
2593 violatedabove = SCIPisFeasPositive(scip, valj * vali - SCIPgetSolVal(scip, sol, terms[idx].aux.var));
2594 }
2595 else
2596 {
2597 assert(bestoverest != NULL);
2598 violatedabove = FALSE;
2599 }
2600 }
2601 else
2602 {
2603 assert(bestoverest[idx] >= 0 && bestoverest[idx] < terms[idx].nauxexprs);
2604
2605 /* if we are here, the relation with the best overestimator must be violated */
2606 assert(SCIPisFeasPositive(scip, valj * vali - SCIPevalBilinAuxExprNonlinear(scip, terms[idx].x, terms[idx].y,
2607 terms[idx].aux.exprs[bestoverest[idx]], sol)));
2608 violatedabove = TRUE;
2609 }
2610
2611 /* only violated products contribute to row marks */
2612 if( !violatedbelow && !violatedabove )
2613 {
2614 SCIPdebugMsg(scip, "the product for vars <%s> and <%s> is not violated\n", SCIPvarGetName(xj), SCIPvarGetName(xi));
2615 continue;
2616 }
2617
2618 /* get the column of xi */
2619 coli = SCIPvarGetCol(xi);
2620 colvals = SCIPcolGetVals(coli);
2621 ncolrows = SCIPcolGetNNonz(coli);
2622 colrows = SCIPcolGetRows(coli);
2623
2624 SCIPdebugMsg(scip, "marking rows for xj = <%s>, xi = <%s>\n", SCIPvarGetName(xj), SCIPvarGetName(xi));
2625
2626 /* mark the rows */
2627 for( r = 0; r < ncolrows; ++r )
2628 {
2629 ridx = SCIProwGetIndex(colrows[r]);
2630
2631 if( !SCIPhashmapExists(row_to_pos, (void*)(size_t)ridx) )
2632 continue; /* if row index is not in row_to_pos, it means that storeSuitableRows decided to ignore this row */
2633
2634 a = colvals[r];
2635 if( a == 0.0 )
2636 continue;
2637
2638 SCIPdebugMsg(scip, "Marking row %d\n", ridx);
2639 addRowMark(ridx, a, violatedbelow, violatedabove, row_idcs, row_marks, nmarked);
2640 }
2641 }
2642
2643 return SCIP_OKAY;
2644 }
2645
2646 /** adds McCormick inequalities for implicit products */
2647 static
2648 SCIP_RETCODE separateMcCormickImplicit(
2649 SCIP* scip, /**< SCIP data structure */
2650 SCIP_SEPA* sepa, /**< separator */
2651 SCIP_SEPADATA* sepadata, /**< separator data */
2652 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
2653 int* bestunderestimators,/**< indices of auxiliary underestimators with largest violation in sol */
2654 int* bestoverestimators, /**< indices of auxiliary overestimators with largest violation in sol */
2655 SCIP_RESULT* result /**< pointer to store the result */
2656 )
2657 {
2658 int i;
2659 int j;
2660 SCIP_CONSNONLINEAR_BILINTERM* terms;
2661 SCIP_ROW* cut;
2662 char name[SCIP_MAXSTRLEN];
2663 SCIP_Bool underestimate;
2664 SCIP_Real xcoef;
2665 SCIP_Real ycoef;
2666 SCIP_Real auxcoef;
2667 SCIP_Real constant;
2668 SCIP_Bool success;
2669 SCIP_CONSNONLINEAR_AUXEXPR* auxexpr;
2670 SCIP_Bool cutoff;
2671 SCIP_Real refpointx;
2672 SCIP_Real refpointy;
2673 SCIP_INTERVAL bndx;
2674 SCIP_INTERVAL bndy;
2675 #ifndef NDEBUG
2676 SCIP_Real productval;
2677 SCIP_Real auxval;
2678 #endif
2679
2680 assert(sepadata->nbilinterms == SCIPgetNBilinTermsNonlinear(sepadata->conshdlr));
2681 assert(bestunderestimators != NULL && bestoverestimators != NULL);
2682
2683 cutoff = FALSE;
2684 terms = SCIPgetBilinTermsNonlinear(sepadata->conshdlr);
2685
2686 for( i = 0; i < sepadata->nbilinterms; ++i )
2687 {
2688 if( terms[i].existing )
2689 continue;
2690
2691 assert(terms[i].nauxexprs > 0);
2692
2693 bndx.inf = SCIPvarGetLbLocal(terms[i].x);
2694 bndx.sup = SCIPvarGetUbLocal(terms[i].x);
2695 bndy.inf = SCIPvarGetLbLocal(terms[i].y);
2696 bndy.sup = SCIPvarGetUbLocal(terms[i].y);
2697 refpointx = SCIPgetSolVal(scip, sol, terms[i].x);
2698 refpointy = SCIPgetSolVal(scip, sol, terms[i].y);
2699
2700 /* adjust the reference points */
2701 refpointx = MIN(MAX(refpointx, bndx.inf), bndx.sup); /*lint !e666*/
2702 refpointy = MIN(MAX(refpointy, bndy.inf), bndy.sup); /*lint !e666*/
2703
2704 /* one iteration for underestimation and one for overestimation */
2705 for( j = 0; j < 2; ++j )
2706 {
2707 /* if underestimate, separate xy <= auxexpr; if !underestimate, separate xy >= auxexpr;
2708 * the cuts will be:
2709 * if underestimate: McCormick_under(xy) - auxexpr <= 0,
2710 * if !underestimate: McCormick_over(xy) - auxexpr >= 0
2711 */
2712 underestimate = j == 0;
2713 if( underestimate && bestoverestimators[i] != -1 )
2714 auxexpr = terms[i].aux.exprs[bestoverestimators[i]];
2715 else if( !underestimate && bestunderestimators[i] != -1 )
2716 auxexpr = terms[i].aux.exprs[bestunderestimators[i]];
2717 else
2718 continue;
2719 assert(!underestimate || auxexpr->overestimate);
2720 assert(underestimate || auxexpr->underestimate);
2721
2722 #ifndef NDEBUG
2723 /* make sure that the term is violated */
2724 productval = SCIPgetSolVal(scip, sol, terms[i].x) * SCIPgetSolVal(scip, sol, terms[i].y);
2725 auxval = SCIPevalBilinAuxExprNonlinear(scip, terms[i].x, terms[i].y, auxexpr, sol);
2726
2727 /* if underestimate, then xy <= aux must be violated; otherwise aux <= xy must be violated */
2728 assert((underestimate && SCIPisFeasLT(scip, auxval, productval)) ||
2729 (!underestimate && SCIPisFeasLT(scip, productval, auxval)));
2730 #endif
2731
2732 /* create an empty row */
2733 (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "mccormick_%sestimate_implicit_%s*%s_%d",
2734 underestimate ? "under" : "over", SCIPvarGetName(terms[i].x), SCIPvarGetName(terms[i].y),
(1) Event invalid_type: |
Argument "SCIPgetNLPs(scip)" to format specifier "%d" was expected to have type "int" but has type "long long". [details] |
2735 SCIPgetNLPs(scip));
2736
2737 SCIP_CALL(SCIPcreateEmptyRowSepa(scip, &cut, sepa, name, -SCIPinfinity(scip), SCIPinfinity(scip), TRUE,
2738 FALSE, FALSE) );
2739
2740 xcoef = 0.0;
2741 ycoef = 0.0;
2742 auxcoef = 0.0;
2743 constant = 0.0;
2744 success = TRUE;
2745
2746 /* subtract auxexpr from the cut */
2747 addAuxexprCoefs(terms[i].x, terms[i].y, auxexpr, -1.0, &auxcoef, &xcoef, &ycoef, &constant);
2748
2749 /* add McCormick terms: ask for an underestimator if relation is xy <= auxexpr, and vice versa */
2750 SCIPaddBilinMcCormick(scip, 1.0, bndx.inf, bndx.sup, refpointx, bndy.inf, bndy.sup, refpointy, !underestimate,
2751 &xcoef, &ycoef, &constant, &success);
2752
2753 if( REALABS(constant) > MAXVARBOUND )
2754 success = FALSE;
2755
2756 if( success )
2757 {
2758 assert(!SCIPisInfinity(scip, REALABS(xcoef)));
2759 assert(!SCIPisInfinity(scip, REALABS(ycoef)));
2760 assert(!SCIPisInfinity(scip, REALABS(constant)));
2761
2762 SCIP_CALL( SCIPaddVarToRow(scip, cut, terms[i].x, xcoef) );
2763 SCIP_CALL( SCIPaddVarToRow(scip, cut, terms[i].y, ycoef) );
2764 SCIP_CALL( SCIPaddVarToRow(scip, cut, auxexpr->auxvar, auxcoef) );
2765
2766 /* set side */
2767 if( underestimate )
2768 SCIP_CALL( SCIPchgRowRhs(scip, cut, -constant) );
2769 else
2770 SCIP_CALL( SCIPchgRowLhs(scip, cut, -constant) );
2771
2772 /* if the cut is violated, add it to SCIP */
2773 if( SCIPisFeasNegative(scip, SCIPgetRowFeasibility(scip, cut)) )
2774 {
2775 SCIP_CALL( SCIPaddRow(scip, cut, FALSE, &cutoff) );
2776 *result = SCIP_SEPARATED;
2777 }
2778 else
2779 {
2780 SCIPdebugMsg(scip, "\nMcCormick cut for hidden product <%s>*<%s> was created successfully, but is not violated",
2781 SCIPvarGetName(terms[i].x), SCIPvarGetName(terms[i].y));
2782 }
2783 }
2784
2785 /* release the cut */
2786 if( cut != NULL )
2787 {
2788 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2789 }
2790
2791 if( cutoff )
2792 {
2793 *result = SCIP_CUTOFF;
2794 SCIPdebugMsg(scip, "exit separator because we found a cutoff -> skip\n");
2795 return SCIP_OKAY;
2796 }
2797 }
2798 }
2799
2800 return SCIP_OKAY;
2801 }
2802
2803 /** builds and adds the RLT cuts */
2804 static
2805 SCIP_RETCODE separateRltCuts(
2806 SCIP* scip, /**< SCIP data structure */
2807 SCIP_SEPA* sepa, /**< separator */
2808 SCIP_SEPADATA* sepadata, /**< separator data */
2809 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */
2810 SCIP_SOL* sol, /**< the point to be separated (can be NULL) */
2811 SCIP_HASHMAP* row_to_pos, /**< hashmap linking row indices to positions in array */
2812 RLT_SIMPLEROW* projrows, /**< projected rows */
2813 SCIP_ROW** rows, /**< problem rows */
2814 int nrows, /**< number of problem rows */
2815 SCIP_Bool allowlocal, /**< are local cuts allowed? */
2816 int* bestunderestimators,/**< indices of auxiliary underestimators with largest violation in sol */
2817 int* bestoverestimators, /**< indices of auxiliary overestimators with largest violation in sol */
2818 SCIP_RESULT* result /**< buffer to store whether separation was successful */
2819 )
2820 {
2821 int j;
2822 int r;
2823 int k;
2824 int nmarked;
2825 int cutssize;
2826 int ncuts;
2827 SCIP_VAR* xj;
2828 unsigned int* row_marks;
2829 int* row_idcs;
2830 SCIP_ROW* cut;
2831 SCIP_ROW** cuts;
2832 SCIP_Bool uselb[4] = {TRUE, TRUE, FALSE, FALSE};
2833 SCIP_Bool uselhs[4] = {TRUE, FALSE, TRUE, FALSE};
2834 SCIP_Bool success;
2835 SCIP_Bool infeasible;
2836 SCIP_Bool accepted;
2837 SCIP_Bool buildeqcut;
2838 SCIP_Bool iseqrow;
2839
2840 assert(!sepadata->useprojection || projrows != NULL);
2841 assert(!sepadata->detecthidden || (bestunderestimators != NULL && bestoverestimators != NULL));
2842
2843 ncuts = 0;
2844 cutssize = 0;
2845 cuts = NULL;
2846 *result = SCIP_DIDNOTFIND;
2847
2848 SCIP_CALL( SCIPallocCleanBufferArray(scip, &row_marks, nrows) );
2849 SCIP_CALL( SCIPallocBufferArray(scip, &row_idcs, nrows) );
2850
2851 /* loop through all variables that appear in bilinear products */
2852 for( j = 0; j < sepadata->nbilinvars && (sepadata->maxusedvars < 0 || j < sepadata->maxusedvars); ++j )
2853 {
2854 xj = sepadata->varssorted[j];
2855
2856 /* mark all rows for multiplier xj */
2857 SCIP_CALL( markRowsXj(scip, sepadata, conshdlr, sol, j, allowlocal, row_to_pos, bestunderestimators,
2858 bestoverestimators, row_marks, row_idcs, &nmarked) );
2859
2860 assert(nmarked <= nrows);
2861
2862 /* generate the projected cut and if it is violated, generate the actual cut */
2863 for( r = 0; r < nmarked; ++r )
2864 {
2865 int pos;
2866 int currentnunknown;
2867 SCIP_ROW* row;
2868
2869 assert(row_marks[r] != 0);
2870 assert(SCIPhashmapExists(row_to_pos, (void*)(size_t) row_idcs[r])); /*lint !e571 */
2871
2872 pos = SCIPhashmapGetImageInt(row_to_pos, (void*)(size_t) row_idcs[r]); /*lint !e571 */
2873 row = rows[pos];
2874 assert(SCIProwGetIndex(row) == row_idcs[r]);
2875
2876 /* check whether this row and var fulfill the conditions */
2877 SCIP_CALL( isAcceptableRow(sepadata, row, xj, ¤tnunknown, &accepted) );
2878 if( !accepted )
2879 {
2880 SCIPdebugMsg(scip, "rejected row <%s> for variable <%s> (introduces too many new products)\n", SCIProwGetName(row), SCIPvarGetName(xj));
2881 row_marks[r] = 0;
2882 continue;
2883 }
2884
2885 SCIPdebugMsg(scip, "accepted row <%s> for variable <%s>\n", SCIProwGetName(rows[r]), SCIPvarGetName(xj));
2886 #ifdef SCIP_DEBUG
2887 SCIP_CALL( SCIPprintRow(scip, rows[r], NULL) );
2888 #endif
2889 iseqrow = SCIPisEQ(scip, SCIProwGetLhs(row), SCIProwGetRhs(row));
2890
2891 /* if all terms are known and it is an equality row, compute equality cut, that is, multiply row with (x-lb) and/or (ub-x) (but see also @todo at top)
2892 * otherwise, multiply row w.r.t. lhs and/or rhs with (x-lb) and/or (ub-x) and estimate product terms that have no aux.var or aux.expr
2893 */
2894 buildeqcut = (currentnunknown == 0 && iseqrow);
2895
2896 /* go over all suitable combinations of sides and bounds and compute the respective cuts */
2897 for( k = 0; k < 4; ++k )
2898 {
2899 /* if equality cuts are possible, lhs and rhs cuts are equal so skip rhs */
2900 if( buildeqcut )
2901 {
2902 if( k != 1 )
2903 continue;
2904 }
2905 /* otherwise which cuts are generated depends on the marks */
2906 else
2907 {
2908 if( row_marks[r] == 1 && uselb[k] == uselhs[k] )
2909 continue;
2910
2911 if( row_marks[r] == 2 && uselb[k] != uselhs[k] )
2912 continue;
2913 }
2914
2915 success = TRUE;
2916 cut = NULL;
2917
2918 SCIPdebugMsg(scip, "row <%s>, uselb = %u, uselhs = %u\n", SCIProwGetName(row), uselb[k], uselhs[k]);
2919
2920 if( sepadata->useprojection )
2921 {
2922 /* if no variables are left in the projected row, the RLT cut will not be violated */
2923 if( projrows[pos].nnonz == 0 )
2924 continue;
2925
2926 /* compute the rlt cut for a projected row first */
2927 SCIP_CALL( computeRltCut(scip, sepa, sepadata, &cut, NULL, &(projrows[pos]), sol, bestunderestimators,
2928 bestoverestimators, xj, &success, uselb[k], uselhs[k], allowlocal, buildeqcut, TRUE) );
2929
2930 /* if the projected cut is not violated, set success to FALSE */
2931 if( cut != NULL )
2932 {
2933 SCIPdebugMsg(scip, "proj cut viol = %g\n", -SCIPgetRowFeasibility(scip, cut));
2934 }
2935 if( cut != NULL && !SCIPisFeasPositive(scip, -SCIPgetRowFeasibility(scip, cut)) )
2936 {
2937 SCIPdebugMsg(scip, "projected cut is not violated, feasibility = %g\n", SCIPgetRowFeasibility(scip, cut));
2938 success = FALSE;
2939 }
2940
2941 /* release the projected cut */
2942 if( cut != NULL )
2943 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2944 }
2945
2946 /* if we don't use projection or if the projected cut was generated successfully and is violated,
2947 * generate the actual cut */
2948 if( success )
2949 {
2950 SCIP_CALL( computeRltCut(scip, sepa, sepadata, &cut, row, NULL, sol, bestunderestimators,
2951 bestoverestimators, xj, &success, uselb[k], uselhs[k], allowlocal, buildeqcut, FALSE) );
2952 }
2953
2954 if( success )
2955 {
2956 success = SCIPisFeasNegative(scip, SCIPgetRowFeasibility(scip, cut)) || (sepadata->addtopool &&
2957 !SCIProwIsLocal(cut));
2958 }
2959
2960 /* if the cut was created successfully and is violated or (if addtopool == TRUE) globally valid,
2961 * it is added to the cuts array */
2962 if( success )
2963 {
2964 if( ncuts + 1 > cutssize )
2965 {
2966 int newsize;
2967
2968 newsize = SCIPcalcMemGrowSize(scip, ncuts + 1);
2969 SCIP_CALL( SCIPreallocBufferArray(scip, &cuts, newsize) );
2970 cutssize = newsize;
2971 }
2972 cuts[ncuts] = cut;
2973 (ncuts)++;
2974 }
2975 else
2976 {
2977 SCIPdebugMsg(scip, "the generation of the cut failed or cut not violated and not added to cutpool\n");
2978 /* release the cut */
2979 if( cut != NULL )
2980 {
2981 SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2982 }
2983 }
2984 }
2985
2986 /* clear row_marks[r] since it will be used for the next multiplier */
2987 row_marks[r] = 0;
2988 }
2989 }
2990
2991 /* if cuts were found, we apply an additional filtering procedure, which is similar to sepastore */
2992 if( ncuts > 0 )
2993 {
2994 int nselectedcuts;
2995 int i;
2996
2997 assert(cuts != NULL);
2998
2999 SCIP_CALL( SCIPselectCutsHybrid(scip, cuts, NULL, NULL, sepadata->goodscore, sepadata->badscore, sepadata->goodmaxparall,
3000 sepadata->maxparall, sepadata->dircutoffdistweight, sepadata->efficacyweight, sepadata->objparalweight,
3001 0.0, ncuts, 0, sepadata->maxncuts == -1 ? ncuts : sepadata->maxncuts, &nselectedcuts) );
3002
3003 for( i = 0; i < ncuts; ++i )
3004 {
3005 assert(cuts[i] != NULL);
3006
3007 if( i < nselectedcuts )
3008 {
3009 /* if selected, add global cuts to the pool and local cuts to the sepastore */
3010 if( SCIProwIsLocal(cuts[i]) || !sepadata->addtopool )
3011 {
3012 SCIP_CALL( SCIPaddRow(scip, cuts[i], FALSE, &infeasible) );
3013
3014 if( infeasible )
3015 {
3016 SCIPdebugMsg(scip, "CUTOFF! The cut <%s> revealed infeasibility\n", SCIProwGetName(cuts[i]));
3017 *result = SCIP_CUTOFF;
3018 }
3019 else
3020 {
3021 SCIPdebugMsg(scip, "SEPARATED: added cut to scip\n");
3022 *result = SCIP_SEPARATED;
3023 }
3024 }
3025 else
3026 {
3027 SCIP_CALL( SCIPaddPoolCut(scip, cuts[i]) );
3028 }
3029 }
3030
3031 /* release current cut */
3032 SCIP_CALL( SCIPreleaseRow(scip, &cuts[i]) );
3033 }
3034 }
3035
3036 SCIPdebugMsg(scip, "exit separator because cut calculation is finished\n");
3037
3038 SCIPfreeBufferArrayNull(scip, &cuts);
3039 SCIPfreeBufferArray(scip, &row_idcs);
3040 SCIPfreeCleanBufferArray(scip, &row_marks);
3041
3042 return SCIP_OKAY;
3043 }
3044
3045 /*
3046 * Callback methods of separator
3047 */
3048
3049 /** copy method for separator plugins (called when SCIP copies plugins) */
3050 static
3051 SCIP_DECL_SEPACOPY(sepaCopyRlt)
3052 { /*lint --e{715}*/
3053 assert(scip != NULL);
3054 assert(sepa != NULL);
3055 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3056
3057 /* call inclusion method of separator */
3058 SCIP_CALL( SCIPincludeSepaRlt(scip) );
3059
3060 return SCIP_OKAY;
3061 }
3062
3063 /** destructor of separator to free user data (called when SCIP is exiting) */
3064 static
3065 SCIP_DECL_SEPAFREE(sepaFreeRlt)
3066 { /*lint --e{715}*/
3067 SCIP_SEPADATA* sepadata;
3068
3069 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3070
3071 sepadata = SCIPsepaGetData(sepa);
3072 assert(sepadata != NULL);
3073
3074 /* free separator data */
3075 SCIPfreeBlockMemory(scip, &sepadata);
3076
3077 SCIPsepaSetData(sepa, NULL);
3078
3079 return SCIP_OKAY;
3080 }
3081
3082 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */
3083 static
3084 SCIP_DECL_SEPAEXITSOL(sepaExitsolRlt)
3085 { /*lint --e{715}*/
3086 SCIP_SEPADATA* sepadata;
3087
3088 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3089
3090 sepadata = SCIPsepaGetData(sepa);
3091 assert(sepadata != NULL);
3092
3093 if( sepadata->iscreated )
3094 {
3095 SCIP_CALL( freeSepaData(scip, sepadata) );
3096 }
3097
3098 return SCIP_OKAY;
3099 }
3100
3101 /** LP solution separation method of separator */
3102 static
3103 SCIP_DECL_SEPAEXECLP(sepaExeclpRlt)
3104 { /*lint --e{715}*/
3105 SCIP_ROW** prob_rows;
3106 SCIP_ROW** rows;
3107 SCIP_SEPADATA* sepadata;
3108 int ncalls;
3109 int nrows;
3110 SCIP_HASHMAP* row_to_pos;
3111 RLT_SIMPLEROW* projrows;
3112
3113 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
3114
3115 sepadata = SCIPsepaGetData(sepa);
3116 *result = SCIP_DIDNOTRUN;
3117
3118 if( sepadata->maxncuts == 0 )
3119 {
3120 SCIPdebugMsg(scip, "exit separator because maxncuts is set to 0\n");
3121 return SCIP_OKAY;
3122 }
3123
3124 /* don't run in a sub-SCIP or in probing */
3125 if( SCIPgetSubscipDepth(scip) > 0 && !sepadata->useinsubscip )
3126 {
3127 SCIPdebugMsg(scip, "exit separator because in sub-SCIP\n");
3128 return SCIP_OKAY;
3129 }
3130
3131 /* don't run in probing */
3132 if( SCIPinProbing(scip) )
3133 {
3134 SCIPdebugMsg(scip, "exit separator because in probing\n");
3135 return SCIP_OKAY;
3136 }
3137
3138 /* only call separator a given number of times at each node */
3139 ncalls = SCIPsepaGetNCallsAtNode(sepa);
3140 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
3141 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
3142 {
3143 SCIPdebugMsg(scip, "exit separator because round limit for this node is reached\n");
3144 return SCIP_OKAY;
3145 }
3146
3147 /* if this is called for the first time, create the sepadata and start the initial separation round */
3148 if( !sepadata->iscreated )
3149 {
3150 *result = SCIP_DIDNOTFIND;
3151 SCIP_CALL( createSepaData(scip, sepadata) );
3152 }
3153 assert(sepadata->iscreated || (sepadata->nbilinvars == 0 && sepadata->nbilinterms == 0));
3154 assert(sepadata->nbilinterms == SCIPgetNBilinTermsNonlinear(sepadata->conshdlr));
3155
3156 /* no bilinear terms available -> skip */
3157 if( sepadata->nbilinvars == 0 )
3158 {
3159 SCIPdebugMsg(scip, "exit separator because there are no known bilinear terms\n");
3160 return SCIP_OKAY;
3161 }
3162
3163 /* only call separator, if we are not close to terminating */
3164 if( SCIPisStopped(scip) )
3165 {
3166 SCIPdebugMsg(scip, "exit separator because we are too close to terminating\n");
3167 return SCIP_OKAY;
3168 }
3169
3170 /* only call separator, if an optimal LP solution is at hand */
3171 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
3172 {
3173 SCIPdebugMsg(scip, "exit separator because there is no LP solution at hand\n");
3174 return SCIP_OKAY;
3175 }
3176
3177 /* get the rows, depending on settings */
3178 if( sepadata->isinitialround || sepadata->onlyoriginal )
3179 {
3180 SCIP_CALL( getOriginalRows(scip, &prob_rows, &nrows) );
3181 }
3182 else
3183 {
3184 SCIP_CALL( SCIPgetLPRowsData(scip, &prob_rows, &nrows) );
3185 }
3186
3187 /* save the suitable rows */
3188 SCIP_CALL( SCIPallocBufferArray(scip, &rows, nrows) );
3189 SCIP_CALL( SCIPhashmapCreate(&row_to_pos, SCIPblkmem(scip), nrows) );
3190
3191 SCIP_CALL( storeSuitableRows(scip, sepa, sepadata, prob_rows, rows, &nrows, row_to_pos, allowlocal) );
3192
3193 if( nrows == 0 ) /* no suitable rows found, free memory and exit */
3194 {
3195 SCIPhashmapFree(&row_to_pos);
3196 SCIPfreeBufferArray(scip, &rows);
3197 if( sepadata->isinitialround || sepadata->onlyoriginal )
3198 {
3199 SCIPfreeBufferArray(scip, &prob_rows);
3200 sepadata->isinitialround = FALSE;
3201 }
3202 return SCIP_OKAY;
3203 }
3204
3205 /* create the projected problem */
3206 if( sepadata->useprojection )
3207 {
3208 SCIP_Bool allcst;
3209
3210 SCIP_CALL( createProjRows(scip, rows, nrows, NULL, &projrows, allowlocal, &allcst) );
3211
3212 /* if all projected rows have only constants left, quit */
3213 if( allcst )
3214 goto TERMINATE;
3215
3216 #ifdef SCIP_DEBUG
3217 printProjRows(scip, projrows, nrows, NULL);
3218 #endif
3219 }
3220 else
3221 {
3222 projrows = NULL;
3223 }
3224
3225 /* separate the cuts */
3226 if( sepadata->detecthidden )
3227 {
3228 int* bestunderestimators;
3229 int* bestoverestimators;
3230
3231 /* if we detect implicit products, a term might have more than one estimator in each direction;
3232 * save the indices of the most violated estimators
3233 */
3234 SCIP_CALL( SCIPallocBufferArray(scip, &bestunderestimators, sepadata->nbilinterms) );
3235 SCIP_CALL( SCIPallocBufferArray(scip, &bestoverestimators, sepadata->nbilinterms) );
3236 getBestEstimators(scip, sepadata, NULL, bestunderestimators, bestoverestimators);
3237
3238 /* also separate McCormick cuts for implicit products */
3239 SCIP_CALL( separateMcCormickImplicit(scip, sepa, sepadata, NULL, bestunderestimators, bestoverestimators,
3240 result) );
3241
3242 if( *result != SCIP_CUTOFF )
3243 {
3244 SCIP_CALL( separateRltCuts(scip, sepa, sepadata, sepadata->conshdlr, NULL, row_to_pos, projrows, rows, nrows,
3245 allowlocal, bestunderestimators, bestoverestimators, result) );
3246 }
3247
3248 SCIPfreeBufferArray(scip, &bestoverestimators);
3249 SCIPfreeBufferArray(scip, &bestunderestimators);
3250 }
3251 else
3252 {
3253 SCIP_CALL( separateRltCuts(scip, sepa, sepadata, sepadata->conshdlr, NULL, row_to_pos, projrows, rows, nrows,
3254 allowlocal, NULL, NULL, result) );
3255 }
3256
3257 TERMINATE:
3258 /* free the projected problem */
3259 if( sepadata->useprojection )
3260 {
3261 freeProjRows(scip, &projrows, nrows);
3262 }
3263
3264 SCIPhashmapFree(&row_to_pos);
3265 SCIPfreeBufferArray(scip, &rows);
3266
3267 if( sepadata->isinitialround || sepadata->onlyoriginal )
3268 {
3269 SCIPfreeBufferArray(scip, &prob_rows);
3270 sepadata->isinitialround = FALSE;
3271 }
3272
3273 return SCIP_OKAY;
3274 }
3275
3276 /*
3277 * separator specific interface methods
3278 */
3279
3280 /** creates the RLT separator and includes it in SCIP */
3281 SCIP_RETCODE SCIPincludeSepaRlt(
3282 SCIP* scip /**< SCIP data structure */
3283 )
3284 {
3285 SCIP_SEPADATA* sepadata;
3286 SCIP_SEPA* sepa;
3287
3288 /* create RLT separator data */
3289 SCIP_CALL( SCIPallocClearBlockMemory(scip, &sepadata) );
3290 sepadata->conshdlr = SCIPfindConshdlr(scip, "nonlinear");
3291 assert(sepadata->conshdlr != NULL);
3292
3293 /* include separator */
3294 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
3295 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpRlt, NULL, sepadata) );
3296
3297 /* set non fundamental callbacks via setter functions */
3298 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyRlt) );
3299 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeRlt) );
3300 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolRlt) );
3301
3302 /* add RLT separator parameters */
3303 SCIP_CALL( SCIPaddIntParam(scip,
3304 "separating/" SEPA_NAME "/maxncuts",
3305 "maximal number of rlt-cuts that are added per round (-1: unlimited)",
3306 &sepadata->maxncuts, FALSE, DEFAULT_MAXNCUTS, -1, INT_MAX, NULL, NULL) );
3307
3308 SCIP_CALL( SCIPaddIntParam(scip,
3309 "separating/" SEPA_NAME "/maxunknownterms",
3310 "maximal number of unknown bilinear terms a row is still used with (-1: unlimited)",
3311 &sepadata->maxunknownterms, FALSE, DEFAULT_MAXUNKNOWNTERMS, -1, INT_MAX, NULL, NULL) );
3312
3313 SCIP_CALL( SCIPaddIntParam(scip,
3314 "separating/" SEPA_NAME "/maxusedvars",
3315 "maximal number of variables used to compute rlt cuts (-1: unlimited)",
3316 &sepadata->maxusedvars, FALSE, DEFAULT_MAXUSEDVARS, -1, INT_MAX, NULL, NULL) );
3317
3318 SCIP_CALL( SCIPaddIntParam(scip,
3319 "separating/" SEPA_NAME "/maxrounds",
3320 "maximal number of separation rounds per node (-1: unlimited)",
3321 &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
3322
3323 SCIP_CALL( SCIPaddIntParam(scip,
3324 "separating/" SEPA_NAME "/maxroundsroot",
3325 "maximal number of separation rounds in the root node (-1: unlimited)",
3326 &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
3327
3328 SCIP_CALL( SCIPaddBoolParam(scip,
3329 "separating/" SEPA_NAME "/onlyeqrows",
3330 "if set to true, only equality rows are used for rlt cuts",
3331 &sepadata->onlyeqrows, FALSE, DEFAULT_ONLYEQROWS, NULL, NULL) );
3332
3333 SCIP_CALL( SCIPaddBoolParam(scip,
3334 "separating/" SEPA_NAME "/onlycontrows",
3335 "if set to true, only continuous rows are used for rlt cuts",
3336 &sepadata->onlycontrows, FALSE, DEFAULT_ONLYCONTROWS, NULL, NULL) );
3337
3338 SCIP_CALL( SCIPaddBoolParam(scip,
3339 "separating/" SEPA_NAME "/onlyoriginal",
3340 "if set to true, only original rows and variables are used",
3341 &sepadata->onlyoriginal, FALSE, DEFAULT_ONLYORIGINAL, NULL, NULL) );
3342
3343 SCIP_CALL( SCIPaddBoolParam(scip,
3344 "separating/" SEPA_NAME "/useinsubscip",
3345 "if set to true, rlt is also used in sub-scips",
3346 &sepadata->useinsubscip, FALSE, DEFAULT_USEINSUBSCIP, NULL, NULL) );
3347
3348 SCIP_CALL( SCIPaddBoolParam(scip,
3349 "separating/" SEPA_NAME "/useprojection",
3350 "if set to true, projected rows are checked first",
3351 &sepadata->useprojection, FALSE, DEFAULT_USEPROJECTION, NULL, NULL) );
3352
3353 SCIP_CALL( SCIPaddBoolParam(scip,
3354 "separating/" SEPA_NAME "/detecthidden",
3355 "if set to true, hidden products are detected and separated by McCormick cuts",
3356 &sepadata->detecthidden, FALSE, DEFAULT_DETECTHIDDEN, NULL, NULL) );
3357
3358 SCIP_CALL( SCIPaddBoolParam(scip,
3359 "separating/" SEPA_NAME "/hiddenrlt",
3360 "whether RLT cuts (TRUE) or only McCormick inequalities (FALSE) should be added for hidden products",
3361 &sepadata->hiddenrlt, FALSE, DEFAULT_HIDDENRLT, NULL, NULL) );
3362
3363 SCIP_CALL( SCIPaddBoolParam(scip,
3364 "separating/" SEPA_NAME "/addtopool",
3365 "if set to true, globally valid RLT cuts are added to the global cut pool",
3366 &sepadata->addtopool, FALSE, DEFAULT_ADDTOPOOL, NULL, NULL) );
3367
3368 SCIP_CALL( SCIPaddRealParam(scip,
3369 "separating/" SEPA_NAME "/goodscore",
3370 "threshold for score of cut relative to best score to be considered good, so that less strict filtering is applied",
3371 &sepadata->goodscore, TRUE, DEFAULT_GOODSCORE, 0.0, 1.0, NULL, NULL) );
3372
3373 SCIP_CALL( SCIPaddRealParam(scip,
3374 "separating/" SEPA_NAME "/badscore",
3375 "threshold for score of cut relative to best score to be discarded",
3376 &sepadata->badscore, TRUE, DEFAULT_BADSCORE, 0.0, 1.0, NULL, NULL) );
3377
3378 SCIP_CALL( SCIPaddRealParam(scip,
3379 "separating/" SEPA_NAME "/objparalweight",
3380 "weight of objective parallelism in cut score calculation",
3381 &sepadata->objparalweight, TRUE, DEFAULT_OBJPARALWEIGHT, 0.0, 1.0, NULL, NULL) );
3382
3383 SCIP_CALL( SCIPaddRealParam(scip,
3384 "separating/" SEPA_NAME "/efficacyweight",
3385 "weight of efficacy in cut score calculation",
3386 &sepadata->efficacyweight, TRUE, DEFAULT_EFFICACYWEIGHT, 0.0, 1.0, NULL, NULL) );
3387
3388 SCIP_CALL( SCIPaddRealParam(scip,
3389 "separating/" SEPA_NAME "/dircutoffdistweight",
3390 "weight of directed cutoff distance in cut score calculation",
3391 &sepadata->dircutoffdistweight, TRUE, DEFAULT_DIRCUTOFFDISTWEIGHT, 0.0, 1.0, NULL, NULL) );
3392
3393 SCIP_CALL( SCIPaddRealParam(scip,
3394 "separating/" SEPA_NAME "/goodmaxparall",
3395 "maximum parallelism for good cuts",
3396 &sepadata->goodmaxparall, TRUE, DEFAULT_GOODMAXPARALL, 0.0, 1.0, NULL, NULL) );
3397
3398 SCIP_CALL( SCIPaddRealParam(scip,
3399 "separating/" SEPA_NAME "/maxparall",
3400 "maximum parallelism for non-good cuts",
3401 &sepadata->maxparall, TRUE, DEFAULT_MAXPARALL, 0.0, 1.0, NULL, NULL) );
3402
3403 return SCIP_OKAY;
3404 }
3405