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