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_minor.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief principal minor separator 28 * @author Benjamin Mueller 29 * 30 * @todo detect non-principal minors and use them to derive split cuts 31 */ 32 33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 34 35 #include <assert.h> 36 #include <string.h> 37 38 #include "scip/sepa_minor.h" 39 #include "scip/cons_nonlinear.h" 40 #include "scip/lapack_calls.h" 41 42 #define SEPA_NAME "minor" 43 #define SEPA_DESC "separator to ensure that 2x2 principal minors of X - xx' are positive semi-definite" 44 #define SEPA_PRIORITY 0 45 #define SEPA_FREQ 10 46 #define SEPA_MAXBOUNDDIST 1.0 47 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 48 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */ 49 50 #define DEFAULT_MAXMINORSCONST 3000 /**< default constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */ 51 #define DEFAULT_MAXMINORSFAC 10.0 /**< default factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */ 52 #define DEFAULT_MINCUTVIOL 1e-4 /**< default minimum required violation of a cut */ 53 #define DEFAULT_RANDSEED 157 /**< default random seed */ 54 #define DEFAULT_MAXROUNDS 10 /**< maximal number of separation rounds per node (-1: unlimited) */ 55 #define DEFAULT_MAXROUNDSROOT -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */ 56 #define DEFAULT_IGNOREPACKINGCONSS TRUE /**< default for ignoring circle packing constraints during minor detection */ 57 58 /* 59 * Data structures 60 */ 61 62 /** separator data */ 63 struct SCIP_SepaData 64 { 65 SCIP_VAR** minors; /**< variables of 2x2 minors; each minor is stored like (auxvar_x^2,auxvar_y^2,auxvar_xy) */ 66 int nminors; /**< total number of minors */ 67 int minorssize; /**< size of minors array */ 68 int maxminorsconst; /**< constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */ 69 SCIP_Real maxminorsfac; /**< factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms) */ 70 int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */ 71 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */ 72 SCIP_Bool detectedminors; /**< has minor detection be called? */ 73 SCIP_Real mincutviol; /**< minimum required violation of a cut */ 74 SCIP_RANDNUMGEN* randnumgen; /**< random number generation */ 75 SCIP_Bool ignorepackingconss; /**< whether to ignore circle packing constraints during minor detection */ 76 }; 77 78 /* 79 * Local methods 80 */ 81 82 /** helper method to store a 2x2 minor in the separation data */ 83 static 84 SCIP_RETCODE sepadataAddMinor( 85 SCIP* scip, /**< SCIP data structure */ 86 SCIP_SEPADATA* sepadata, /**< separator data */ 87 SCIP_VAR* x, /**< x variable */ 88 SCIP_VAR* y, /**< y variable */ 89 SCIP_VAR* auxvarxx, /**< auxiliary variable for x*x */ 90 SCIP_VAR* auxvaryy, /**< auxiliary variable for y*y */ 91 SCIP_VAR* auxvarxy /**< auxiliary variable for x*y */ 92 ) 93 { 94 assert(sepadata != NULL); 95 assert(x != NULL); 96 assert(y != NULL); 97 assert(x != y); 98 assert(auxvarxx != NULL); 99 assert(auxvaryy != NULL); 100 assert(auxvarxy != NULL); 101 assert(auxvarxx != auxvaryy); 102 assert(auxvarxx != auxvarxy); 103 assert(auxvaryy != auxvarxy); 104 105 SCIPdebugMsg(scip, "store 2x2 minor: %s %s %s for x=%s y=%s\n", SCIPvarGetName(auxvarxx), SCIPvarGetName(auxvaryy), 106 SCIPvarGetName(auxvarxy), SCIPvarGetName(x), SCIPvarGetName(y)); 107 108 /* reallocate if necessary */ 109 if( sepadata->minorssize < 5 * (sepadata->nminors + 1) ) 110 { 111 int newsize = SCIPcalcMemGrowSize(scip, 5 * (sepadata->nminors + 1)); 112 assert(newsize > 5 * (sepadata->nminors + 1)); 113 114 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->minors), sepadata->minorssize, newsize) ); 115 sepadata->minorssize = newsize; 116 } 117 118 /* store minor */ 119 sepadata->minors[5 * sepadata->nminors] = x; 120 sepadata->minors[5 * sepadata->nminors + 1] = y; 121 sepadata->minors[5 * sepadata->nminors + 2] = auxvarxx; 122 sepadata->minors[5 * sepadata->nminors + 3] = auxvaryy; 123 sepadata->minors[5 * sepadata->nminors + 4] = auxvarxy; 124 ++(sepadata->nminors); 125 126 /* capture variables */ 127 SCIP_CALL( SCIPcaptureVar(scip, x) ); 128 SCIP_CALL( SCIPcaptureVar(scip, y) ); 129 SCIP_CALL( SCIPcaptureVar(scip, auxvarxx) ); 130 SCIP_CALL( SCIPcaptureVar(scip, auxvaryy) ); 131 SCIP_CALL( SCIPcaptureVar(scip, auxvarxy) ); 132 133 return SCIP_OKAY; 134 } 135 136 /** helper method to clear separation data */ 137 static 138 SCIP_RETCODE sepadataClear( 139 SCIP* scip, /**< SCIP data structure */ 140 SCIP_SEPADATA* sepadata /**< separator data */ 141 ) 142 { 143 int i; 144 145 assert(sepadata != NULL); 146 147 SCIPdebugMsg(scip, "clear separation data\n"); 148 149 /* release captured variables */ 150 for( i = 0; i < 5 * sepadata->nminors; ++i ) 151 { 152 assert(sepadata->minors[i] != NULL); 153 SCIP_CALL( SCIPreleaseVar(scip, &sepadata->minors[i]) ); 154 } 155 156 /* free memory */ 157 SCIPfreeBlockMemoryArrayNull(scip, &sepadata->minors, sepadata->minorssize); 158 159 /* reset counters */ 160 sepadata->nminors = 0; 161 sepadata->minorssize = 0; 162 163 return SCIP_OKAY; 164 } 165 166 /** helper method to identify non-overlapping constraints in circle packing */ 167 static 168 SCIP_Bool isPackingCons( 169 SCIP* scip, /**< SCIP data structure */ 170 SCIP_CONS* cons /**< nonlinear constraint */ 171 ) 172 { 173 SCIP_EXPR* root; 174 SCIP_VAR* quadvars[4] = {NULL, NULL, NULL, NULL}; 175 SCIP_VAR* bilinvars[4] = {NULL, NULL, NULL, NULL}; 176 int nbilinvars = 0; 177 int nquadvars = 0; 178 int nchildren; 179 int i; 180 181 assert(scip != NULL); 182 assert(cons != NULL); 183 184 root = SCIPgetExprNonlinear(cons); 185 assert(root != NULL); 186 nchildren = SCIPexprGetNChildren(root); 187 188 /* non-overlapping constraint has 6 terms (2 bilinear + 4 quadratic) */ 189 if( nchildren != 6 || !SCIPisExprSum(scip, root) ) 190 return FALSE; 191 192 for( i = 0; i < nchildren; ++i ) 193 { 194 SCIP_EXPR* expr; 195 SCIP_EXPR** children; 196 197 /* get child */ 198 expr = SCIPexprGetChildren(root)[i]; 199 assert(expr != NULL); 200 children = SCIPexprGetChildren(expr); 201 202 /* case: expr = x^2; x is no auxiliary variable */ 203 if( SCIPisExprPower(scip, expr) && SCIPgetExponentExprPow(expr) == 2.0 204 && SCIPisExprVar(scip, children[0]) ) 205 { 206 SCIP_VAR* x; 207 208 /* too many quadratic variables -> stop */ 209 if( nquadvars > 3 ) 210 return FALSE; 211 212 x = SCIPgetVarExprVar(children[0]); 213 assert(x != NULL); 214 215 quadvars[nquadvars++] = x; 216 } 217 /* case: expr = x * y; x and y are no auxiliary variables */ 218 else if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) == 2 219 && SCIPisExprVar(scip, children[0]) && SCIPisExprVar(scip, children[1]) ) 220 { 221 SCIP_VAR* x; 222 SCIP_VAR* y; 223 224 /* too many bilinear variables -> stop */ 225 if( nbilinvars > 2 ) 226 return FALSE; 227 228 x = SCIPgetVarExprVar(children[0]); 229 assert(x != NULL); 230 y = SCIPgetVarExprVar(children[1]); 231 assert(y != NULL); 232 assert(x != y); 233 234 bilinvars[nbilinvars++] = x; 235 bilinvars[nbilinvars++] = y; 236 } 237 else 238 { 239 return FALSE; 240 } 241 } 242 243 /* number of bilinear and quadratic terms do not fit */ 244 if( nbilinvars != 4 || nquadvars != 4 ) 245 return FALSE; 246 247 /* each quadratic variable has to appear in exactly one bilinear terms */ 248 for( i = 0; i < nquadvars; ++i ) 249 { 250 int counter = 0; 251 int j; 252 253 for( j = 0; j < nbilinvars; ++j ) 254 { 255 if( quadvars[i] == bilinvars[j] ) 256 ++counter; 257 } 258 259 if( counter != 1 ) 260 return FALSE; 261 } 262 263 return TRUE; 264 } 265 266 /** helper method to get the variables associated to a minor */ 267 static 268 SCIP_RETCODE getMinorVars( 269 SCIP_SEPADATA* sepadata, /**< separator data */ 270 int idx, /**< index of the stored minor */ 271 SCIP_VAR** x, /**< pointer to store x variable */ 272 SCIP_VAR** y, /**< pointer to store x variable */ 273 SCIP_VAR** auxvarxx, /**< pointer to store auxiliary variable for x*x */ 274 SCIP_VAR** auxvaryy, /**< pointer to store auxiliary variable for y*y */ 275 SCIP_VAR** auxvarxy /**< pointer to store auxiliary variable for x*y */ 276 ) 277 { 278 assert(sepadata != NULL); 279 assert(idx >= 0 && idx < sepadata->nminors); 280 assert(auxvarxx != NULL); 281 assert(auxvaryy != NULL); 282 assert(auxvarxy != NULL); 283 284 *x = sepadata->minors[5 * idx]; 285 *y = sepadata->minors[5 * idx + 1]; 286 *auxvarxx = sepadata->minors[5 * idx + 2]; 287 *auxvaryy = sepadata->minors[5 * idx + 3]; 288 *auxvarxy = sepadata->minors[5 * idx + 4]; 289 290 return SCIP_OKAY; 291 } 292 293 /** method to detect and store principal minors */ 294 static 295 SCIP_RETCODE detectMinors( 296 SCIP* scip, /**< SCIP data structure */ 297 SCIP_SEPADATA* sepadata /**< separator data */ 298 ) 299 { 300 SCIP_CONSHDLR* conshdlr; 301 SCIP_EXPRITER* it; 302 SCIP_HASHMAP* quadmap; 303 SCIP_VAR** xs; 304 SCIP_VAR** ys; 305 SCIP_VAR** auxvars; 306 int* perm = NULL; 307 int nbilinterms = 0; 308 int nquadterms = 0; 309 int maxminors; 310 int c; 311 int i; 312 313 #ifdef SCIP_STATISTIC 314 SCIP_Real totaltime = -SCIPgetTotalTime(scip); 315 #endif 316 317 assert(sepadata != NULL); 318 319 /* check whether minor detection has been called already */ 320 if( sepadata->detectedminors ) 321 return SCIP_OKAY; 322 323 assert(sepadata->minors == NULL); 324 assert(sepadata->nminors == 0); 325 326 /* we assume that the auxiliary variables in the nonlinear constraint handler have been already generated */ 327 sepadata->detectedminors = TRUE; 328 329 /* check whether there are nonlinear constraints available */ 330 conshdlr = SCIPfindConshdlr(scip, "nonlinear"); 331 if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 ) 332 return SCIP_OKAY; 333 334 SCIPdebugMsg(scip, "call detectMinors()\n"); 335 336 /* allocate memory */ 337 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 338 SCIP_CALL( SCIPhashmapCreate(&quadmap, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 339 SCIP_CALL( SCIPallocBufferArray(scip, &xs, SCIPgetNVars(scip)) ); 340 SCIP_CALL( SCIPallocBufferArray(scip, &ys, SCIPgetNVars(scip)) ); 341 SCIP_CALL( SCIPallocBufferArray(scip, &auxvars, SCIPgetNVars(scip)) ); 342 343 /* initialize iterator */ 344 SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) ); 345 SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_ENTEREXPR); 346 347 for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c ) 348 { 349 SCIP_CONS* cons; 350 SCIP_EXPR* expr; 351 SCIP_EXPR* root; 352 353 cons = SCIPconshdlrGetConss(conshdlr)[c]; 354 assert(cons != NULL); 355 root = SCIPgetExprNonlinear(cons); 356 assert(root != NULL); 357 358 /* ignore circle packing constraints; the motivation for this is that in circle packing instance not only the SDP 359 * relaxation is weak (see "Packing circles in a square: a theoretical comparison of various convexification 360 * techniques", http://www.optimization-online.org/DB_HTML/2017/03/5911.html), but it also hurts performance 361 */ 362 if( sepadata->ignorepackingconss && isPackingCons(scip, cons) ) 363 { 364 SCIPdebugMsg(scip, "ignore packing constraints %s\n", SCIPconsGetName(cons)); 365 continue; 366 } 367 368 for( expr = SCIPexpriterRestartDFS(it, root); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*/ /*lint !e440*/ 369 { 370 SCIP_EXPR** children; 371 SCIP_VAR* auxvar; 372 373 SCIPdebugMsg(scip, "visit expression %p in constraint %s\n", (void*)expr, SCIPconsGetName(cons)); 374 375 /* check whether the expression has an auxiliary variable */ 376 auxvar = SCIPgetExprAuxVarNonlinear(expr); 377 if( auxvar == NULL ) 378 { 379 SCIPdebugMsg(scip, "expression has no auxiliary variable -> skip\n"); 380 continue; 381 } 382 383 children = SCIPexprGetChildren(expr); 384 385 /* check for expr = (x)^2 */ 386 if( SCIPexprGetNChildren(expr) == 1 && SCIPisExprPower(scip, expr) 387 && SCIPgetExponentExprPow(expr) == 2.0 388 && SCIPgetExprAuxVarNonlinear(children[0]) != NULL ) 389 { 390 SCIP_VAR* quadvar; 391 392 assert(children[0] != NULL); 393 394 quadvar = SCIPgetExprAuxVarNonlinear(children[0]); 395 assert(quadvar != NULL); 396 assert(!SCIPhashmapExists(quadmap, (void*)quadvar)); 397 SCIPdebugMsg(scip, "found %s = (%s)^2\n", SCIPvarGetName(auxvar), SCIPvarGetName(quadvar)); 398 399 /* hash the quadratic variable to its corresponding auxiliary variable */ 400 SCIP_CALL( SCIPhashmapInsert(quadmap, (void*)quadvar, auxvar) ); 401 ++nquadterms; 402 } 403 /* check for expr = x * y */ 404 else if( SCIPexprGetNChildren(expr) == 2 && SCIPisExprProduct(scip, expr) 405 && SCIPgetExprAuxVarNonlinear(children[0]) != NULL && SCIPgetExprAuxVarNonlinear(children[1]) != NULL ) 406 { 407 SCIP_VAR* x; 408 SCIP_VAR* y; 409 410 assert(children[0] != NULL); 411 assert(children[1] != NULL); 412 413 x = SCIPgetExprAuxVarNonlinear(children[0]); 414 y = SCIPgetExprAuxVarNonlinear(children[1]); 415 416 /* ignore binary variables */ 417 if( !SCIPvarIsBinary(x) && !SCIPvarIsBinary(y) ) 418 { 419 xs[nbilinterms] = SCIPgetExprAuxVarNonlinear(children[0]); 420 ys[nbilinterms] = SCIPgetExprAuxVarNonlinear(children[1]); 421 auxvars[nbilinterms] = auxvar; 422 SCIPdebugMsg(scip, "found %s = %s * %s\n", SCIPvarGetName(auxvar), SCIPvarGetName(xs[nbilinterms]), SCIPvarGetName(ys[nbilinterms])); 423 ++nbilinterms; 424 } 425 } 426 } 427 } 428 assert(nbilinterms < SCIPgetNVars(scip)); 429 SCIPdebugMsg(scip, "stored %d bilinear terms in total\n", nbilinterms); 430 431 /* use max(maxminorsconst, maxminorsfac * # quadratic terms) as a limit for the maximum number of minors */ 432 maxminors = (int) MAX(sepadata->maxminorsconst, sepadata->maxminorsfac * nquadterms); 433 SCIPdebugMsg(scip, "maximum number of minors = %d\n", maxminors); 434 435 /* permute bilinear terms if there are too many of them; the motivation for this is that we don't want to 436 * prioritize variables because of the order in the bilinear terms where they appear; however, variables that 437 * appear more often in bilinear terms might be more important than others so the corresponding bilinear terms 438 * are more likely to be chosen 439 */ 440 if( maxminors < nbilinterms && maxminors < SQR(nquadterms) ) 441 { 442 SCIP_CALL( SCIPallocBufferArray(scip, &perm, nbilinterms) ); 443 444 for( i = 0; i < nbilinterms; ++i ) 445 perm[i] = i; 446 447 /* permute array */ 448 SCIPrandomPermuteIntArray(sepadata->randnumgen, perm, 0, nbilinterms); 449 } 450 451 /* store 2x2 principal minors */ 452 for( i = 0; i < nbilinterms && sepadata->nminors < maxminors; ++i ) 453 { 454 SCIP_VAR* x; 455 SCIP_VAR* y; 456 SCIP_VAR* auxvarxy; 457 458 if( perm == NULL ) 459 { 460 x = xs[i]; 461 y = ys[i]; 462 auxvarxy = auxvars[i]; 463 } 464 else 465 { 466 x = xs[perm[i]]; 467 y = ys[perm[i]]; 468 auxvarxy = auxvars[perm[i]]; 469 } 470 471 assert(x != NULL); 472 assert(y != NULL); 473 assert(auxvarxy != NULL); 474 assert(x != y); 475 476 if( SCIPhashmapExists(quadmap, (void*)x) && SCIPhashmapExists(quadmap, (void*)y) ) 477 { 478 SCIP_VAR* auxvarxx; 479 SCIP_VAR* auxvaryy; 480 481 auxvarxx = (SCIP_VAR*)SCIPhashmapGetImage(quadmap, (void*)x); 482 assert(auxvarxx != NULL); 483 auxvaryy = (SCIP_VAR*)SCIPhashmapGetImage(quadmap, (void*)y); 484 assert(auxvaryy != NULL); 485 486 /* store minor into the separation data */ 487 SCIP_CALL( sepadataAddMinor(scip, sepadata, x, y, auxvarxx, auxvaryy, auxvarxy) ); 488 } 489 } 490 SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors); 491 492 /* free memory */ 493 SCIPfreeBufferArrayNull(scip, &perm); 494 SCIPfreeBufferArray(scip, &auxvars); 495 SCIPfreeBufferArray(scip, &ys); 496 SCIPfreeBufferArray(scip, &xs); 497 SCIPhashmapFree(&quadmap); 498 SCIPfreeExpriter(&it); 499 500 #ifdef SCIP_STATISTIC 501 totaltime += SCIPgetTotalTime(scip); 502 SCIPstatisticMessage("MINOR DETECT %s %f %d %d\n", SCIPgetProbName(scip), totaltime, sepadata->nminors, maxminors); 503 #endif 504 505 return SCIP_OKAY; 506 } 507 508 /** helper method to compute eigenvectors and eigenvalues */ 509 static 510 SCIP_RETCODE getEigenValues( 511 SCIP* scip, /**< SCIP data structure */ 512 SCIP_Real x, /**< solution value of x */ 513 SCIP_Real y, /**< solution value of y */ 514 SCIP_Real xx, /**< solution value of x*x */ 515 SCIP_Real yy, /**< solution value of y*y */ 516 SCIP_Real xy, /**< solution value of x*y */ 517 SCIP_Real* eigenvals, /**< array to store eigenvalues (at least of size 3) */ 518 SCIP_Real* eigenvecs, /**< array to store eigenvalues (at least of size 9) */ 519 SCIP_Bool* success /**< pointer to store whether eigenvalue computation was successful */ 520 ) 521 { 522 assert(eigenvals != NULL); 523 assert(eigenvecs != NULL); 524 assert(success != NULL); 525 526 *success = TRUE; 527 528 /* construct matrix */ 529 eigenvecs[0] = 1.0; 530 eigenvecs[1] = x; 531 eigenvecs[2] = y; 532 eigenvecs[3] = x; 533 eigenvecs[4] = xx; 534 eigenvecs[5] = xy; 535 eigenvecs[6] = y; 536 eigenvecs[7] = xy; 537 eigenvecs[8] = yy; 538 539 /* use LAPACK to compute the eigenvalues and eigenvectors */ 540 if( SCIPlapackComputeEigenvalues(SCIPbuffer(scip), TRUE, 3, eigenvecs, eigenvals) != SCIP_OKAY ) 541 { 542 SCIPdebugMsg(scip, "Failed to compute eigenvalues and eigenvectors of augmented quadratic form matrix.\n"); 543 *success = FALSE; 544 } 545 546 return SCIP_OKAY; 547 } 548 549 /** generate and add a cut */ 550 static 551 SCIP_RETCODE addCut( 552 SCIP* scip, /**< SCIP data structure */ 553 SCIP_SEPA* sepa, /**< separator */ 554 SCIP_SOL* sol, /**< solution to separate (might be NULL) */ 555 SCIP_VAR* x, /**< x variable */ 556 SCIP_VAR* y, /**< y variable */ 557 SCIP_VAR* xx, /**< auxiliary variable for x*x */ 558 SCIP_VAR* yy, /**< auxiliary variable for y*y */ 559 SCIP_VAR* xy, /**< auxiliary variable for x*y */ 560 SCIP_Real* eigenvec, /**< array containing an eigenvector */ 561 SCIP_Real eigenval, /**< eigenvalue */ 562 SCIP_Real mincutviol, /**< minimal required violation */ 563 SCIP_RESULT* result /**< pointer to update the result */ 564 ) 565 { 566 SCIP_VAR* vars[5] = {x, y, xx, yy, xy}; 567 SCIP_Real coefs[5]; 568 SCIP_Real constant; 569 SCIP_ROWPREP* rowprep; 570 SCIP_Bool success; 571 572 assert(x != NULL); 573 assert(y != NULL); 574 assert(xx != NULL); 575 assert(yy != NULL); 576 assert(xy != NULL); 577 assert(eigenvec != NULL); 578 assert(mincutviol >= 0.0); 579 assert(result != NULL); 580 581 /* check whether the resulting cut is violated enough */ 582 if( !SCIPisFeasLT(scip, eigenval, -mincutviol) ) 583 return SCIP_OKAY; 584 585 /* the resulting cut reads as 586 * (1 x y ) (v0) 587 * (v0 v1 v2) (x xx xy) (v1) >= 0 588 * (y xy yy) (v2) 589 * where v is the eigenvector corresponding to a negative eigenvalue 590 * that is, 591 * v0^2 + 2 v0 v1 * x + 2 v0 v2 * y + v1^2 * xx + v2^2 * yy + 2 v1 v2 * xy >= 0 592 */ 593 constant = SQR(eigenvec[0]); 594 coefs[0] = 2.0 * eigenvec[0] * eigenvec[1]; 595 coefs[1] = 2.0 * eigenvec[0] * eigenvec[2]; 596 coefs[2] = SQR(eigenvec[1]); 597 coefs[3] = SQR(eigenvec[2]); 598 coefs[4] = 2.0 * eigenvec[1] * eigenvec[2]; 599 600 /* create rowprep */ 601 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, FALSE) ); 602 SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, 5, vars, coefs) ); 603 SCIProwprepAddConstant(rowprep, constant); 604 SCIPdebug( SCIPprintRowprep(scip, rowprep, NULL) ); 605 SCIPdebugMsg(scip, "cut violation %g mincutviol = %g\n", SCIPgetRowprepViolation(scip, rowprep, sol, NULL), mincutviol); 606 607 /* cleanup coefficient and side, esp treat epsilon to integral values; don't consider scaling up here */ 608 SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, 0.0, NULL, &success) ); 609 610 /* check cut violation */ 611 if( success && SCIPgetRowprepViolation(scip, rowprep, sol, NULL) > mincutviol ) 612 { 613 SCIP_ROW* row; 614 SCIP_Bool infeasible; 615 616 /* set name of rowprep */ 617 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "minor_%s_%s_%s_%lld", SCIPvarGetName(xx), SCIPvarGetName(yy), 618 SCIPvarGetName(xy), SCIPgetNLPs(scip)); 619 620 /* create, add, and release row */ 621 SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) ); 622 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) ); 623 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 624 625 /* update result pointer */ 626 *result = infeasible ? SCIP_CUTOFF : SCIP_SEPARATED; 627 } 628 629 /* free rowprep */ 630 SCIPfreeRowprep(scip, &rowprep); 631 632 return SCIP_OKAY; 633 } 634 635 /** separates cuts for stored principal minors */ 636 static 637 SCIP_RETCODE separatePoint( 638 SCIP* scip, /**< SCIP data structure */ 639 SCIP_SEPA* sepa, /**< separator */ 640 SCIP_SOL* sol, /**< primal solution that should be separated, or NULL for LP solution */ 641 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 642 ) 643 { 644 SCIP_SEPADATA* sepadata; 645 int i; 646 647 assert(sepa != NULL); 648 assert(result != NULL); 649 650 *result = SCIP_DIDNOTRUN; 651 652 sepadata = SCIPsepaGetData(sepa); 653 assert(sepadata != NULL); 654 655 /* check whether there are some minors available */ 656 if( sepadata->nminors == 0 ) 657 return SCIP_OKAY; 658 659 *result = SCIP_DIDNOTFIND; 660 661 for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i ) 662 { 663 SCIP_Real eigenvals[3]; 664 SCIP_Real eigenvecs[9]; 665 SCIP_VAR* x; 666 SCIP_VAR* y; 667 SCIP_VAR* xx; 668 SCIP_VAR* yy; 669 SCIP_VAR* xy; 670 SCIP_Real solx; 671 SCIP_Real soly; 672 SCIP_Real solxx; 673 SCIP_Real solyy; 674 SCIP_Real solxy; 675 SCIP_Bool success; 676 int k; 677 678 /* get variables of the i-th minor */ 679 SCIP_CALL( getMinorVars(sepadata, i, &x, &y, &xx, &yy, &xy) ); 680 assert(x != NULL); 681 assert(y != NULL); 682 assert(xx != NULL); 683 assert(yy != NULL); 684 assert(xy != NULL); 685 686 /* get current solution values */ 687 solx = SCIPgetSolVal(scip, sol, x); 688 soly = SCIPgetSolVal(scip, sol, y); 689 solxx = SCIPgetSolVal(scip, sol, xx); 690 solyy = SCIPgetSolVal(scip, sol, yy); 691 solxy = SCIPgetSolVal(scip, sol, xy); 692 SCIPdebugMsg(scip, "solution values (x,y,xx,yy,xy)=(%g,%g,%g,%g,%g)\n", solx, soly, solxx, solyy, solxy); 693 694 /* compute eigenvalues and eigenvectors */ 695 SCIP_CALL( getEigenValues(scip, solx, soly, solxx, solyy, solxy, eigenvals, eigenvecs, &success) ); 696 if( !success ) 697 continue; 698 699 /* try to generate a cut for each negative eigenvalue */ 700 for( k = 0; k < 3 && (*result != SCIP_CUTOFF); ++k ) 701 { 702 SCIPdebugMsg(scip, "eigenvalue = %g eigenvector = (%g,%g,%g)\n", eigenvals[k], eigenvecs[3*k], eigenvecs[3*k + 1], eigenvecs[3*k + 2]); 703 SCIP_CALL( addCut(scip, sepa, sol, x, y, xx, yy, xy, &eigenvecs[3*k], eigenvals[k], sepadata->mincutviol, result) ); 704 SCIPdebugMsg(scip, "result: %d\n", *result); 705 } 706 } 707 708 return SCIP_OKAY; 709 } 710 711 /* 712 * Callback methods of separator 713 */ 714 715 /** copy method for separator plugins (called when SCIP copies plugins) */ 716 static 717 SCIP_DECL_SEPACOPY(sepaCopyMinor) 718 { /*lint --e{715}*/ 719 assert(scip != NULL); 720 assert(sepa != NULL); 721 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 722 723 /* call inclusion method of constraint handler */ 724 SCIP_CALL( SCIPincludeSepaMinor(scip) ); 725 726 return SCIP_OKAY; 727 } 728 729 730 /** destructor of separator to free user data (called when SCIP is exiting) */ 731 static 732 SCIP_DECL_SEPAFREE(sepaFreeMinor) 733 { /*lint --e{715}*/ 734 SCIP_SEPADATA* sepadata; 735 736 sepadata = SCIPsepaGetData(sepa); 737 assert(sepadata != NULL); 738 assert(sepadata->minors == NULL); 739 assert(sepadata->nminors == 0); 740 assert(sepadata->minorssize == 0); 741 742 /* free separator data */ 743 SCIPfreeBlockMemory(scip, &sepadata); 744 SCIPsepaSetData(sepa, NULL); 745 746 return SCIP_OKAY; 747 } 748 749 750 /** initialization method of separator (called after problem was transformed) */ 751 static 752 SCIP_DECL_SEPAINIT(sepaInitMinor) 753 { /*lint --e{715}*/ 754 SCIP_SEPADATA* sepadata; 755 756 /* get separator data */ 757 sepadata = SCIPsepaGetData(sepa); 758 assert(sepadata != NULL); 759 assert(sepadata->randnumgen == NULL); 760 761 /* create random number generator */ 762 SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) ); 763 764 return SCIP_OKAY; 765 } 766 767 768 /** deinitialization method of separator (called before transformed problem is freed) */ 769 static 770 SCIP_DECL_SEPAEXIT(sepaExitMinor) 771 { /*lint --e{715}*/ 772 SCIP_SEPADATA* sepadata; 773 774 /* get separator data */ 775 sepadata = SCIPsepaGetData(sepa); 776 assert(sepadata != NULL); 777 assert(sepadata->randnumgen != NULL); 778 779 /* free random number generator */ 780 SCIPfreeRandom(scip, &sepadata->randnumgen); 781 782 return SCIP_OKAY; 783 } 784 785 786 /** solving process initialization method of separator (called when branch and bound process is about to begin) */ 787 static 788 SCIP_DECL_SEPAINITSOL(sepaInitsolMinor) 789 { /*lint --e{715}*/ 790 return SCIP_OKAY; 791 } 792 793 794 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */ 795 static 796 SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor) 797 { /*lint --e{715}*/ 798 SCIP_SEPADATA* sepadata; 799 800 sepadata = SCIPsepaGetData(sepa); 801 assert(sepadata != NULL); 802 803 /* clear separation data */ 804 SCIP_CALL( sepadataClear(scip, sepadata) ); 805 806 return SCIP_OKAY; 807 } 808 809 810 /** LP solution separation method of separator */ 811 static 812 SCIP_DECL_SEPAEXECLP(sepaExeclpMinor) 813 { /*lint --e{715}*/ 814 SCIP_SEPADATA* sepadata; 815 int ncalls; 816 817 /* need routine to compute eigenvalues/eigenvectors */ 818 if( ! SCIPlapackIsAvailable() ) 819 return SCIP_OKAY; 820 821 sepadata = SCIPsepaGetData(sepa); 822 assert(sepadata != NULL); 823 ncalls = SCIPsepaGetNCallsAtNode(sepa); 824 825 /* only call the separator a given number of times at each node */ 826 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 827 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 828 { 829 SCIPdebugMsg(scip, "reached round limit for node\n"); 830 return SCIP_OKAY; 831 } 832 833 /* try to detect minors */ 834 SCIP_CALL( detectMinors(scip, sepadata) ); 835 836 /* call separation method */ 837 SCIP_CALL( separatePoint(scip, sepa, NULL, result) ); 838 839 return SCIP_OKAY; 840 } 841 842 843 /** arbitrary primal solution separation method of separator */ 844 static 845 SCIP_DECL_SEPAEXECSOL(sepaExecsolMinor) 846 { /*lint --e{715}*/ 847 SCIP_SEPADATA* sepadata; 848 int ncalls; 849 850 /* need routine to compute eigenvalues/eigenvectors */ 851 if( ! SCIPlapackIsAvailable() ) 852 return SCIP_OKAY; 853 854 sepadata = SCIPsepaGetData(sepa); 855 assert(sepadata != NULL); 856 ncalls = SCIPsepaGetNCallsAtNode(sepa); 857 858 /* only call the separator a given number of times at each node */ 859 if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 860 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 861 { 862 SCIPdebugMsg(scip, "reached round limit for node\n"); 863 return SCIP_OKAY; 864 } 865 866 /* try to detect minors */ 867 SCIP_CALL( detectMinors(scip, SCIPsepaGetData(sepa)) ); 868 869 /* call separation method */ 870 SCIP_CALL( separatePoint(scip, sepa, sol, result) ); 871 872 return SCIP_OKAY; 873 } 874 875 /* 876 * separator specific interface methods 877 */ 878 879 /** creates the minor separator and includes it in SCIP */ 880 SCIP_RETCODE SCIPincludeSepaMinor( 881 SCIP* scip /**< SCIP data structure */ 882 ) 883 { 884 SCIP_SEPADATA* sepadata = NULL; 885 SCIP_SEPA* sepa = NULL; 886 887 /* create minor separator data */ 888 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) ); 889 BMSclearMemory(sepadata); 890 891 /* include separator */ 892 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 893 SEPA_USESSUBSCIP, SEPA_DELAY, 894 sepaExeclpMinor, sepaExecsolMinor, 895 sepadata) ); 896 897 assert(sepa != NULL); 898 899 /* set non fundamental callbacks via setter functions */ 900 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) ); 901 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) ); 902 SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) ); 903 SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) ); 904 SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) ); 905 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) ); 906 907 /* add minor separator parameters */ 908 SCIP_CALL( SCIPaddIntParam(scip, 909 "separating/" SEPA_NAME "/maxminorsconst", 910 "constant for the maximum number of minors, i.e., max(const, fac * # quadratic terms)", 911 &sepadata->maxminorsconst, FALSE, DEFAULT_MAXMINORSCONST, 0, INT_MAX, NULL, NULL) ); 912 913 SCIP_CALL( SCIPaddRealParam(scip, 914 "separating/" SEPA_NAME "/maxminorsfac", 915 "factor for the maximum number of minors, i.e., max(const, fac * # quadratic terms)", 916 &sepadata->maxminorsfac, FALSE, DEFAULT_MAXMINORSFAC, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 917 918 SCIP_CALL( SCIPaddRealParam(scip, 919 "separating/" SEPA_NAME "/mincutviol", 920 "minimum required violation of a cut", 921 &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 922 923 SCIP_CALL( SCIPaddIntParam(scip, 924 "separating/" SEPA_NAME "/maxrounds", 925 "maximal number of separation rounds per node (-1: unlimited)", 926 &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) ); 927 928 SCIP_CALL( SCIPaddIntParam(scip, 929 "separating/" SEPA_NAME "/maxroundsroot", 930 "maximal number of separation rounds in the root node (-1: unlimited)", 931 &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) ); 932 933 SCIP_CALL( SCIPaddBoolParam(scip, 934 "separating/" SEPA_NAME "/ignorepackingconss", 935 "whether to ignore circle packing constraints during minor detection", 936 &sepadata->ignorepackingconss, FALSE, DEFAULT_IGNOREPACKINGCONSS, NULL, NULL) ); 937 938 return SCIP_OKAY; 939 } 940