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_interminor.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief minor separator with intersection cuts 28 * @author Felipe Serrano 29 * @author Antonia Chmiela 30 * 31 * Let X be the matrix of auxiliary variables added for bilinear terms, X_{ij} = x_ix_j. 32 * The separator enforces quadratic constraints det(2x2 minor of X) = 0 via intersection cuts. 33 */ 34 35 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 36 37 #include <assert.h> 38 #include <string.h> 39 40 #include "scip/sepa_interminor.h" 41 #include "scip/expr.h" 42 #include "scip/expr_var.h" 43 #include "scip/expr_pow.h" 44 #include "scip/expr_product.h" 45 #include "scip/nlpi_ipopt.h" 46 #include "scip/cons_nonlinear.h" 47 48 49 50 #define SEPA_NAME "interminor" 51 #define SEPA_DESC "intersection cuts separator to ensure that 2x2 minors of X (= xx') have determinant 0" 52 #define SEPA_PRIORITY 0 53 #define SEPA_FREQ -1 54 #define SEPA_MAXBOUNDDIST 1.0 55 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 56 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */ 57 58 #define DEFAULT_MINCUTVIOL 1e-4 /**< default minimum required violation of a cut */ 59 #define DEFAULT_RANDSEED 157 /**< default random seed */ 60 #define DEFAULT_MAXROUNDS 10 /**< maximal number of separation rounds per node (-1: unlimited) */ 61 #define DEFAULT_MAXROUNDSROOT -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */ 62 #define BINSEARCH_MAXITERS 120 /**< default iteration limit for binary search */ 63 #define DEFAULT_USESTRENGTHENING FALSE /**< default for using strengthend intersection cuts to separate */ 64 #define DEFAULT_USEBOUNDS FALSE /**< default for using nonnegativity bounds when separating */ 65 66 /* 67 * Data structures 68 */ 69 70 /** separator data */ 71 struct SCIP_SepaData 72 { 73 SCIP_VAR** minors; /**< variables of 2x2 minors; each minor is stored like (auxvar_x^2,auxvar_y^2,auxvar_xy) */ 74 SCIP_Bool* isdiagonal; /**< bool array determining if the variables appearing in the minor are diagonal */ 75 int nminors; /**< total number of minors */ 76 int minorssize; /**< size of minors array */ 77 int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */ 78 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */ 79 SCIP_Bool detectedminors; /**< has minor detection be called? */ 80 SCIP_Real mincutviol; /**< minimum required violation of a cut */ 81 SCIP_RANDNUMGEN* randnumgen; /**< random number generation */ 82 SCIP_Bool usestrengthening; /**< whether to use strengthened intersection cuts to separate minors */ 83 SCIP_Bool usebounds; /**< whether to also enforce nonegativity bounds of principle minors */ 84 }; 85 86 /* these represent a row */ 87 struct rowdata 88 { 89 int* vals; /**< index of the column */ 90 int rowidx; /**< index corresponding to variable of that row */ 91 int nvals; /**< number of nonzero entries in column */ 92 int valssize; /**< size of the array that is currently allocated */ 93 SCIP_HASHMAP* auxvars; /**< entry of the matrix */ 94 }; 95 96 /* 97 * Local methods 98 */ 99 100 /** helper method to store a 2x2 minor in the separation data */ 101 static 102 SCIP_RETCODE sepadataAddMinor( 103 SCIP* scip, /**< SCIP data structure */ 104 SCIP_SEPADATA* sepadata, /**< separator data */ 105 SCIP_VAR* auxvarxik, /**< auxiliary variable X_ik = x_i * x_k */ 106 SCIP_VAR* auxvarxil, /**< auxiliary variable X_il = x_i * x_l */ 107 SCIP_VAR* auxvarxjk, /**< auxiliary variable X_jk = x_j * x_k */ 108 SCIP_VAR* auxvarxjl, /**< auxiliary variable X_jl = x_j * x_l */ 109 SCIP_Bool isauxvarxikdiag, /**< is X_ik diagonal? (i.e. i = k) */ 110 SCIP_Bool isauxvarxildiag, /**< is X_il diagonal? (i.e. i = l) */ 111 SCIP_Bool isauxvarxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */ 112 SCIP_Bool isauxvarxjldiag /**< is X_jl diagonal? (i.e. j = l) */ 113 ) 114 { 115 assert(sepadata != NULL); 116 assert(auxvarxik != NULL); 117 assert(auxvarxil != NULL); 118 assert(auxvarxjk != NULL); 119 assert(auxvarxjl != NULL); 120 assert(auxvarxik != auxvarxil); 121 assert(auxvarxjk != auxvarxjl); 122 123 SCIPdebugMsg(scip, "store 2x2 minor: [%s %s, %s %s]\n", SCIPvarGetName(auxvarxik), SCIPvarGetName(auxvarxil), 124 SCIPvarGetName(auxvarxjk), SCIPvarGetName(auxvarxjl)); 125 126 /* reallocate if necessary */ 127 if( sepadata->minorssize < 4 * (sepadata->nminors + 1) ) 128 { 129 int newsize = SCIPcalcMemGrowSize(scip, 4 * (sepadata->nminors + 1)); 130 assert(newsize >= 4 * (sepadata->nminors + 1)); 131 132 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->minors), sepadata->minorssize, newsize) ); 133 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &(sepadata->isdiagonal), sepadata->minorssize, newsize) ); 134 sepadata->minorssize = newsize; 135 } 136 137 /* store minor */ 138 sepadata->minors[4 * sepadata->nminors] = auxvarxik; 139 sepadata->minors[4 * sepadata->nminors + 1] = auxvarxil; 140 sepadata->minors[4 * sepadata->nminors + 2] = auxvarxjk; 141 sepadata->minors[4 * sepadata->nminors + 3] = auxvarxjl; 142 sepadata->isdiagonal[4 * sepadata->nminors] = isauxvarxikdiag; 143 sepadata->isdiagonal[4 * sepadata->nminors + 1] = isauxvarxildiag; 144 sepadata->isdiagonal[4 * sepadata->nminors + 2] = isauxvarxjkdiag; 145 sepadata->isdiagonal[4 * sepadata->nminors + 3] = isauxvarxjldiag; 146 ++(sepadata->nminors); 147 148 /* capture variables */ 149 SCIP_CALL( SCIPcaptureVar(scip, auxvarxik) ); 150 SCIP_CALL( SCIPcaptureVar(scip, auxvarxil) ); 151 SCIP_CALL( SCIPcaptureVar(scip, auxvarxjk) ); 152 SCIP_CALL( SCIPcaptureVar(scip, auxvarxjl) ); 153 154 return SCIP_OKAY; 155 } 156 157 /** helper method to clear separation data */ 158 static 159 SCIP_RETCODE sepadataClear( 160 SCIP* scip, /**< SCIP data structure */ 161 SCIP_SEPADATA* sepadata /**< separator data */ 162 ) 163 { 164 int i; 165 166 assert(sepadata != NULL); 167 168 SCIPdebugMsg(scip, "clear separation data\n"); 169 170 /* release captured variables */ 171 for( i = 0; i < 4 * sepadata->nminors; ++i ) 172 { 173 assert(sepadata->minors[i] != NULL); 174 SCIP_CALL( SCIPreleaseVar(scip, &sepadata->minors[i]) ); 175 } 176 177 /* free memory */ 178 SCIPfreeBlockMemoryArrayNull(scip, &sepadata->minors, sepadata->minorssize); 179 SCIPfreeBlockMemoryArrayNull(scip, &sepadata->isdiagonal, sepadata->minorssize); 180 181 /* reset counters */ 182 sepadata->nminors = 0; 183 sepadata->minorssize = 0; 184 185 return SCIP_OKAY; 186 } 187 188 /** helper method to get the variables associated to a minor */ 189 static 190 SCIP_RETCODE getMinorVars( 191 SCIP_SEPADATA* sepadata, /**< separator data */ 192 int idx, /**< index of the stored minor */ 193 SCIP_VAR** auxvarxik, /**< auxiliary variable X_ik = x_i * x_k */ 194 SCIP_VAR** auxvarxil, /**< auxiliary variable X_il = x_i * x_l */ 195 SCIP_VAR** auxvarxjk, /**< auxiliary variable X_jk = x_j * x_k */ 196 SCIP_VAR** auxvarxjl, /**< auxiliary variable X_jl = x_j * x_l */ 197 SCIP_Bool* isauxvarxikdiag, /**< is X_ik diagonal? (i.e. i = k) */ 198 SCIP_Bool* isauxvarxildiag, /**< is X_il diagonal? (i.e. i = l) */ 199 SCIP_Bool* isauxvarxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */ 200 SCIP_Bool* isauxvarxjldiag /**< is X_jl diagonal? (i.e. j = l) */ 201 ) 202 { 203 assert(auxvarxik != NULL); 204 assert(auxvarxil != NULL); 205 assert(auxvarxjk != NULL); 206 assert(auxvarxjl != NULL); 207 208 *auxvarxik = sepadata->minors[4 * idx]; 209 *auxvarxil = sepadata->minors[4 * idx + 1]; 210 *auxvarxjk = sepadata->minors[4 * idx + 2]; 211 *auxvarxjl = sepadata->minors[4 * idx + 3]; 212 213 *isauxvarxikdiag = sepadata->isdiagonal[4 * idx]; 214 *isauxvarxildiag = sepadata->isdiagonal[4 * idx + 1]; 215 *isauxvarxjkdiag = sepadata->isdiagonal[4 * idx + 2]; 216 *isauxvarxjldiag = sepadata->isdiagonal[4 * idx + 3]; 217 218 return SCIP_OKAY; 219 } 220 221 222 /** adds a new entry (i.e., auxvar) of in (row, col) of matrix M. 223 * 224 * we have a matrix, M, indexed by the variables 225 * M(xi, xk) is the auxiliary variable of xi * xk if it exists 226 * We store, for each row of the matrix, the indices of the nonzero column entries (assoc with the given row) and the auxiliary variable for xi * xk 227 * The nonzero column entries are stored as an array (struct rowdata) 228 * So we have a hashmap mapping each variable (row of the matrix) with its array representing the nonzero entries of the row. 229 */ 230 static 231 SCIP_RETCODE insertIndex( 232 SCIP* scip, /**< SCIP data structure */ 233 SCIP_HASHMAP* rowmap, /**< hashmap of the rows of the matrix */ 234 SCIP_VAR* row, /**< variable corresponding to row of new entry */ 235 SCIP_VAR* col, /**< variable corresponding to column of new entry */ 236 SCIP_VAR* auxvar, /**< auxvar to insert into the matrix */ 237 int* rowindices, /**< array of indices of all variables corresponding to a row */ 238 int* nrows /**< number of rows */ 239 ) 240 { 241 SCIPdebugMsg(scip, "inserting %s in row %s and col %s \n", SCIPvarGetName(auxvar), SCIPvarGetName(row), SCIPvarGetName(col)); 242 243 /* check whether variable has an array associated to it */ 244 if( SCIPhashmapExists(rowmap, (void*)row) ) 245 { 246 struct rowdata* arr; 247 248 arr = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)row); 249 250 /* reallocate if necessary */ 251 if( arr->valssize < arr->nvals + 1 ) 252 { 253 int newsize = SCIPcalcMemGrowSize(scip, arr->nvals + 1); 254 assert(newsize > arr->nvals + 1); 255 256 SCIP_CALL( SCIPreallocBufferArray(scip, &(arr->vals), newsize) ); 257 arr->valssize = newsize; 258 } 259 260 /* insert */ 261 arr->vals[arr->nvals] = SCIPvarGetProbindex(col); 262 SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) ); 263 arr->nvals += 1; 264 } 265 else 266 { 267 struct rowdata* arr; 268 269 /* create index array */ 270 SCIP_CALL( SCIPallocBuffer(scip, &arr) ); 271 arr->valssize = 10; 272 arr->nvals = 0; 273 SCIP_CALL( SCIPallocBufferArray(scip, &arr->vals, arr->valssize) ); 274 SCIP_CALL( SCIPhashmapCreate(&arr->auxvars, SCIPblkmem(scip), arr->valssize) ); 275 276 /* insert */ 277 arr->rowidx = SCIPvarGetProbindex(row); 278 arr->vals[arr->nvals] = SCIPvarGetProbindex(col); 279 SCIP_CALL( SCIPhashmapInsert(arr->auxvars, (void*)col, (void *)auxvar) ); 280 arr->nvals += 1; 281 282 /* store in hashmap */ 283 SCIP_CALL( SCIPhashmapInsert(rowmap, (void*)row, (void *)arr) ); 284 285 /* remember the new row */ 286 rowindices[*nrows] = SCIPvarGetProbindex(row); 287 *nrows += 1; 288 } 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* rowmap; 303 int* rowvars = NULL; 304 int* intersection; 305 int nrowvars = 0; 306 int c; 307 int i; 308 309 #ifdef SCIP_STATISTIC 310 SCIP_Real totaltime = -SCIPgetTotalTime(scip); 311 #endif 312 313 assert(sepadata != NULL); 314 315 /* check whether minor detection has been called already */ 316 if( sepadata->detectedminors ) 317 return SCIP_OKAY; 318 319 assert(sepadata->minors == NULL); 320 assert(sepadata->nminors == 0); 321 322 /* we assume that the auxiliary variables in the nonlinear constraint handler have been already generated */ 323 sepadata->detectedminors = TRUE; 324 325 /* check whether there are nonlinear constraints available */ 326 conshdlr = SCIPfindConshdlr(scip, "nonlinear"); 327 if( conshdlr == NULL || SCIPconshdlrGetNConss(conshdlr) == 0 ) 328 return SCIP_OKAY; 329 330 SCIPdebugMsg(scip, "call detectMinors()\n"); 331 332 /* allocate memory */ 333 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 334 SCIP_CALL( SCIPhashmapCreate(&rowmap, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 335 SCIP_CALL( SCIPallocBufferArray(scip, &rowvars, SCIPgetNVars(scip)) ); 336 SCIP_CALL( SCIPallocBufferArray(scip, &intersection, SCIPgetNVars(scip)) ); 337 338 /* initialize iterator */ 339 SCIP_CALL( SCIPexpriterInit(it, NULL, SCIP_EXPRITER_DFS, FALSE) ); 340 341 for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c ) 342 { 343 SCIP_CONS* cons; 344 SCIP_EXPR* expr; 345 SCIP_EXPR* root; 346 347 cons = SCIPconshdlrGetConss(conshdlr)[c]; 348 assert(cons != NULL); 349 root = SCIPgetExprNonlinear(cons); 350 assert(root != NULL); 351 352 for( expr = SCIPexpriterRestartDFS(it, root); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*//*lint !e440*/ 353 { 354 SCIP_EXPR** children; 355 SCIP_VAR* auxvar; 356 357 SCIPdebugMsg(scip, "visit expression %p in constraint %s\n", (void*)expr, SCIPconsGetName(cons)); 358 359 /* check whether the expression has an auxiliary variable */ 360 auxvar = SCIPgetExprAuxVarNonlinear(expr); 361 if( auxvar == NULL ) 362 { 363 SCIPdebugMsg(scip, "expression has no auxiliary variable -> skip\n"); 364 continue; 365 } 366 367 children = SCIPexprGetChildren(expr); 368 369 /* check for expr = (x)^2 */ 370 if( SCIPexprGetNChildren(expr) == 1 && SCIPisExprPower(scip, expr) 371 && SCIPgetExponentExprPow(expr) == 2.0 372 && SCIPgetExprAuxVarNonlinear(children[0]) != NULL ) 373 { 374 SCIP_VAR* quadvar; 375 376 assert(children[0] != NULL); 377 378 quadvar = SCIPgetExprAuxVarNonlinear(children[0]); 379 assert(quadvar != NULL); 380 381 SCIP_CALL( insertIndex(scip, rowmap, quadvar, quadvar, auxvar, rowvars, &nrowvars) ); 382 } 383 /* check for expr = x_i * x_k */ 384 else if( SCIPexprGetNChildren(expr) == 2 && SCIPisExprProduct(scip, expr) 385 && SCIPgetExprAuxVarNonlinear(children[0]) != NULL && SCIPgetExprAuxVarNonlinear(children[1]) != NULL ) 386 { 387 SCIP_VAR* xi; 388 SCIP_VAR* xk; 389 390 assert(children[0] != NULL); 391 assert(children[1] != NULL); 392 393 xi = SCIPgetExprAuxVarNonlinear(children[0]); 394 xk = SCIPgetExprAuxVarNonlinear(children[1]); 395 396 SCIP_CALL( insertIndex(scip, rowmap, xk, xi, auxvar, rowvars, &nrowvars) ); 397 SCIP_CALL( insertIndex(scip, rowmap, xi, xk, auxvar, rowvars, &nrowvars) ); 398 } 399 } 400 } 401 402 /* sort the column entries */ 403 for( i = 0; i < nrowvars; ++i ) 404 { 405 struct rowdata* row; 406 407 row = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]); 408 SCIPsortInt(row->vals, row->nvals); 409 } 410 411 /* store 2x2 minors */ 412 /* TODO: we might store some minors twice since the matrix is symmetric. Handle that! (see unit test for example) */ 413 for( i = 0; i < nrowvars; ++i ) 414 { 415 int j; 416 struct rowdata* rowi; 417 418 rowi = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[i]]); 419 420 for( j = i + 1; j < nrowvars; ++j ) 421 { 422 struct rowdata* rowj; 423 int ninter; 424 425 rowj = (struct rowdata*)SCIPhashmapGetImage(rowmap, (void *)SCIPgetVars(scip)[rowvars[j]]); 426 427 SCIPcomputeArraysIntersectionInt(rowi->vals, rowi->nvals, rowj->vals, rowj->nvals, intersection, &ninter); 428 429 if( ninter > 1) 430 { 431 int p; 432 433 for( p = 0; p < ninter - 1; ++p ) 434 { 435 int q; 436 437 for( q = p + 1; q < ninter; ++q ) 438 { 439 SCIP_HASHMAP* rowicols; 440 SCIP_HASHMAP* rowjcols; 441 SCIP_VAR* colk; 442 SCIP_VAR* coll; 443 SCIP_VAR* auxvarik; 444 SCIP_VAR* auxvaril; 445 SCIP_VAR* auxvarjk; 446 SCIP_VAR* auxvarjl; 447 int ii; 448 int jj; 449 int k; 450 int l; 451 SCIP_Bool isauxvarikdiag = FALSE; 452 SCIP_Bool isauxvarildiag = FALSE; 453 SCIP_Bool isauxvarjkdiag = FALSE; 454 SCIP_Bool isauxvarjldiag = FALSE; 455 456 ii = rowi->rowidx; 457 jj = rowj->rowidx; 458 k = intersection[p]; 459 l = intersection[q]; 460 461 rowicols = rowi->auxvars; 462 rowjcols = rowj->auxvars; 463 464 colk = SCIPgetVars(scip)[k]; 465 coll = SCIPgetVars(scip)[l]; 466 467 auxvarik = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, colk); 468 auxvaril = (SCIP_VAR*) SCIPhashmapGetImage(rowicols, coll); 469 auxvarjk = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, colk); 470 auxvarjl = (SCIP_VAR*) SCIPhashmapGetImage(rowjcols, coll); 471 472 if( ii == k ) 473 isauxvarikdiag = TRUE; 474 else if( ii == l ) 475 isauxvarildiag = TRUE; 476 if( jj == k ) 477 isauxvarjkdiag = TRUE; 478 else if( jj == l ) 479 isauxvarjldiag = TRUE; 480 481 SCIP_CALL( sepadataAddMinor(scip, sepadata, auxvarik, auxvaril, auxvarjk, auxvarjl, 482 isauxvarikdiag, isauxvarildiag, isauxvarjkdiag, isauxvarjldiag) ); 483 } 484 } 485 } 486 } 487 SCIPfreeBufferArrayNull(scip, &rowi->vals); 488 SCIPhashmapFree(&rowi->auxvars); 489 SCIPfreeBufferArrayNull(scip, &rowi); 490 } 491 492 SCIPdebugMsg(scip, "found %d principal minors in total\n", sepadata->nminors); 493 494 /* free memory */ 495 SCIPfreeBufferArray(scip, &intersection); 496 SCIPfreeBufferArray(scip, &rowvars); 497 SCIPhashmapFree(&rowmap); 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 /** constructs map between lp position of a basic variable and its row in the tableau */ 509 static 510 SCIP_RETCODE constructBasicVars2TableauRowMap( 511 SCIP* scip, /**< SCIP data structure */ 512 int* map /**< buffer to store the map */ 513 ) 514 { 515 int* basisind; 516 int nrows; 517 int i; 518 519 nrows = SCIPgetNLPRows(scip); 520 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) ); 521 522 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) ); 523 for( i = 0; i < nrows; ++i ) 524 { 525 if( basisind[i] >= 0 ) 526 map[basisind[i]] = i; 527 } 528 529 SCIPfreeBufferArray(scip, &basisind); 530 531 return SCIP_OKAY; 532 } 533 534 /** The restriction of the function representing the maximal S-free set to zlp + t * ray has the form 535 * SQRT(A t^2 + B t + C) - (D t + E). 536 * This function computes the coefficients A, B, C, D, E for the given ray. 537 */ 538 static 539 SCIP_RETCODE computeRestrictionToRay( 540 SCIP* scip, /**< SCIP data structure */ 541 SCIP_Real* ray, /**< coefficients of ray */ 542 SCIP_VAR** vars, /**< variables */ 543 SCIP_Real* coefs, /**< buffer to store A, B, C, D, and E of cases 1, 2, 3, or 4a*/ 544 SCIP_Real* coefs4b, /**< buffer to store A, B, C, D, and E of case 4b */ 545 SCIP_Real* coefscondition, /**< buffer to store coefs for checking whether we are in case 4a or 4b */ 546 SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */ 547 SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */ 548 SCIP_Bool* success /**< FALSE if we need to abort generation because of numerics */ 549 ) 550 { 551 SCIP_Real eigenvectors[16] = {1.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 1.0, -1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0}; 552 SCIP_Real eigenvalues[4] = {0.5, 0.5, -0.5, -0.5}; 553 SCIP_Real eigencoef = 0.7071067811865475244008443621048490; 554 SCIP_Real* a; 555 SCIP_Real* b; 556 SCIP_Real* c; 557 SCIP_Real* d; 558 SCIP_Real* e; 559 SCIP_Real min; 560 SCIP_Real max; 561 SCIP_Real norm1; 562 SCIP_Real norm2; 563 int negidx; 564 int posidx; 565 int i; 566 567 *success = TRUE; 568 569 /* set all coefficients to zero */ 570 memset(coefs, 0, 5 * sizeof(SCIP_Real)); 571 memset(coefs4b, 0, 5 * sizeof(SCIP_Real)); 572 norm1 = 0.0; 573 norm2 = 0.0; 574 575 a = coefs; 576 b = coefs + 1; 577 c = coefs + 2; 578 d = coefs + 3; 579 e = coefs + 4; 580 581 negidx = 2; 582 posidx = 0; 583 for( i = 0; i < 4; ++i ) 584 { 585 int j; 586 SCIP_Real vzlp; 587 SCIP_Real vdotray; 588 589 vzlp = 0; 590 vdotray = 0; 591 592 /* compute eigenvec * ray and eigenvec * solution */ 593 for( j = 0; j < 4; ++j ) 594 { 595 vdotray += eigencoef * eigenvectors[4 * i + j] * ray[j]; 596 vzlp += eigencoef * eigenvectors[4 * i + j] * SCIPvarGetLPSol(vars[j]); 597 } 598 599 if( eigenvalues[i] > 0 ) 600 { 601 /* positive eigenvalue: compute D and E */ 602 *d += eigenvalues[i] * vzlp * vdotray; 603 *e += eigenvalues[i] * SQR( vzlp ); 604 605 if( usebounds ) 606 { 607 norm1 += eigenvalues[i] * (1 - SQR( ad[posidx] )) * SQR( vzlp ); 608 norm2 += SQRT( eigenvalues[i] ) * ad[posidx] * vzlp; 609 ++posidx; 610 } 611 } 612 else 613 { 614 /* negative eigenvalue: compute A, B, and C */ 615 *a -= eigenvalues[i] * SQR( vdotray ); 616 *b -= 2.0 * eigenvalues[i] * vzlp * vdotray; 617 *c -= eigenvalues[i] * SQR( vzlp ); 618 619 if( usebounds ) 620 { 621 coefs4b[0] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vdotray ); 622 coefs4b[1] -= 2.0 * eigenvalues[i] * (1 - SQR( ad[negidx] )) * vzlp * vdotray; 623 coefs4b[2] -= eigenvalues[i] * (1 - SQR( ad[negidx] )) * SQR( vzlp ); 624 coefs4b[3] += SQRT( -eigenvalues[i] ) * ad[negidx] * vdotray; 625 coefs4b[4] += SQRT( -eigenvalues[i] ) * ad[negidx] * vzlp; 626 ++negidx; 627 } 628 } 629 } 630 631 assert(*e > 0); 632 633 if( SQRT( *c ) - SQRT( *e ) >= 0.0 ) 634 { 635 assert(SQRT( *c ) - SQRT( *e ) < 1e-6); 636 *success = FALSE; 637 return SCIP_OKAY; 638 } 639 640 /* finish computation of coefficients when using bounds */ 641 if( usebounds ) 642 { 643 coefscondition[0] = norm2 / SQRT( *e ); 644 coefscondition[1] = coefs4b[3]; 645 coefscondition[2] = coefs4b[4]; 646 647 coefs4b[0] *= norm1 / *e; 648 coefs4b[1] *= norm1 / *e; 649 coefs4b[2] *= norm1 / *e; 650 coefs4b[3] *= norm2 / SQRT( *e ); 651 coefs4b[4] *= norm2 / SQRT( *e ); 652 653 coefs4b[3] += *d / SQRT( *e ); 654 coefs4b[4] += SQRT( *e ); 655 656 assert( SQRT( coefs4b[2] ) - coefs4b[4] < 0.0 ); 657 } 658 659 /* finish computation of D and E */ 660 *e = SQRT( *e ); 661 *d /= *e; 662 663 /* maybe we want to avoid a large dynamism between A, B and C */ 664 max = 0.0; 665 min = SCIPinfinity(scip); 666 for( i = 0; i < 3; ++i ) 667 { 668 SCIP_Real absval; 669 670 absval = ABS(coefs[i]); 671 if( max < absval ) 672 max = absval; 673 if( absval != 0.0 && absval < min ) 674 min = absval; 675 } 676 677 if( SCIPisHugeValue(scip, max / min) ) 678 { 679 #ifdef DEBUG_INTERSECTIONCUT 680 printf("Bad numerics: max(A,B,C)/min(A,B,C) is too large (%g)\n", max / min); 681 #endif 682 *success = FALSE; 683 return SCIP_OKAY; 684 } 685 686 /* some sanity checks */ 687 assert(*c >= 0); /* radicand at zero */ 688 assert(SQRT( *c ) - *e < 0); /* the function at 0 must be negative */ 689 assert(*a >= 0); /* the function inside the root is convex */ 690 691 #ifdef DEBUG_INTERSECTIONCUT 692 SCIPinfoMessage(scip, NULL, "Restriction yields: a,b,c,d,e %g %g %g %g %g\n", coefs[0], coefs[1], coefs[2], coefs[3], coefs[4]); 693 #endif 694 695 return SCIP_OKAY; 696 } 697 698 /** returns phi(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) */ /*lint -e{715}*/ 699 static 700 SCIP_Real evalPhiAtRay( 701 SCIP* scip, /**< SCIP data structure */ 702 SCIP_Real t, /**< argument of phi restricted to ray */ 703 SCIP_Real a, /**< value of A */ 704 SCIP_Real b, /**< value of B */ 705 SCIP_Real c, /**< value of C */ 706 SCIP_Real d, /**< value of D */ 707 SCIP_Real e /**< value of E */ 708 ) 709 { 710 #ifdef INTERCUTS_DBLDBL 711 SCIP_Real QUAD(lin); 712 SCIP_Real QUAD(disc); 713 SCIP_Real QUAD(tmp); 714 SCIP_Real QUAD(root); 715 716 /* d * t + e */ 717 SCIPquadprecProdDD(lin, d, t); 718 SCIPquadprecSumQD(lin, lin, e); 719 720 /* a * t * t */ 721 SCIPquadprecSquareD(disc, t); 722 SCIPquadprecProdQD(disc, disc, a); 723 724 /* b * t */ 725 SCIPquadprecProdDD(tmp, b, t); 726 727 /* a * t * t + b * t */ 728 SCIPquadprecSumQQ(disc, disc, tmp); 729 730 /* a * t * t + b * t + c */ 731 SCIPquadprecSumQD(disc, disc, c); 732 733 /* sqrt(above): can't take sqrt of 0! */ 734 if( QUAD_TO_DBL(disc) == 0 ) 735 { 736 QUAD_ASSIGN(root, 0.0); 737 } 738 else 739 { 740 SCIPquadprecSqrtQ(root, disc); 741 } 742 743 /* final result */ 744 QUAD_SCALE(lin, -1.0); 745 SCIPquadprecSumQQ(tmp, root, lin); 746 747 assert(!SCIPisInfinity(scip, t) || QUAD_TO_DBL(tmp) <= 0); 748 749 return QUAD_TO_DBL(tmp); 750 #else 751 return SQRT( a * t * t + b * t + c ) - ( d * t + e ); 752 #endif 753 } 754 755 /** helper function of computeRoot: we want phi to be <= 0 */ 756 static 757 void doBinarySearch( 758 SCIP* scip, /**< SCIP data structure */ 759 SCIP_Real a, /**< value of A */ 760 SCIP_Real b, /**< value of B */ 761 SCIP_Real c, /**< value of C */ 762 SCIP_Real d, /**< value of D */ 763 SCIP_Real e, /**< value of E */ 764 SCIP_Real* sol /**< buffer to store solution; also gives initial point */ 765 ) 766 { 767 SCIP_Real lb = 0.0; 768 SCIP_Real ub = *sol; 769 SCIP_Real curr; 770 int i; 771 772 for( i = 0; i < BINSEARCH_MAXITERS; ++i ) 773 { 774 SCIP_Real phival; 775 776 curr = (lb + ub) / 2.0; 777 phival = evalPhiAtRay(scip, curr, a, b, c, d, e); 778 #ifdef INTERCUT_MOREDEBUG 779 printf("%d: lb,ub %.10f, %.10f. curr = %g -> phi at curr %g -> phi at lb %g \n", i, lb, ub, curr, phival, evalPhiAtRay(scip, lb, a, b, c, d, e)); 780 #endif 781 782 if( phival <= 0.0 ) 783 { 784 lb = curr; 785 if( SCIPisFeasZero(scip, phival) || SCIPisFeasEQ(scip, ub, lb) ) 786 break; 787 } 788 else 789 ub = curr; 790 } 791 792 *sol = lb; 793 } 794 795 /** checks if we are in case 4a, i.e., if 796 * (num(xhat_{r+1}(zlp)) / E) * SQRT(A * tsol^2 + B * tsol + C) + w(ray) * tsol + num(yhat_{s+1}(zlp)) <= 0 797 */ 798 static 799 SCIP_Real isCase4a( 800 SCIP_Real tsol, /**< t in the above formula */ 801 SCIP_Real* coefs, /**< coefficients A, B, C, D, and E of case 4a */ 802 SCIP_Real* coefscondition /**< extra coefficients needed for the evaluation of the condition: 803 * num(xhat_{r+1}(zlp)) / E; w(ray); num(yhat_{s+1}(zlp)) */ 804 ) 805 { 806 return (coefscondition[0] * SQRT( coefs[0] * SQR( tsol ) + coefs[1] * tsol + coefs[2] ) + coefscondition[1] * 807 tsol + coefscondition[2]) <= 0.0; 808 } 809 810 /** finds smallest positive root phi by finding the smallest positive root of 811 * (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) = 0 812 * 813 * However, we are conservative and want a solution such that phi is negative, but close to 0; 814 * thus we correct the result with a binary search 815 */ 816 static 817 SCIP_Real computeRoot( 818 SCIP* scip, /**< SCIP data structure */ 819 SCIP_Real* coefs /**< value of A */ 820 ) 821 { 822 SCIP_Real sol; 823 SCIP_INTERVAL bounds; 824 SCIP_INTERVAL result; 825 SCIP_Real a = coefs[0]; 826 SCIP_Real b = coefs[1]; 827 SCIP_Real c = coefs[2]; 828 SCIP_Real d = coefs[3]; 829 SCIP_Real e = coefs[4]; 830 831 /* there is an intersection point if and only if SQRT(A) > D: here we are beliving in math, this might cause 832 * numerical issues 833 */ 834 if( SQRT( a ) <= d ) 835 { 836 sol = SCIPinfinity(scip); 837 838 return sol; 839 } 840 841 SCIPintervalSetBounds(&bounds, 0.0, SCIPinfinity(scip)); 842 843 /* SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar finds all x such that a x^2 + b x >= c and x in bounds. 844 * it is known that if tsol is the root we are looking for, then gamma(zlp + t * ray) <= 0 between 0 and tsol, thus 845 * tsol is the smallest t such that (A - D^2) t^2 + (B - 2 D*E) t + (C - E^2) >= 0 846 */ 847 SCIPintervalSolveUnivariateQuadExpressionPositiveAllScalar(SCIP_INTERVAL_INFINITY, &result, a - d * d, b - 2.0 * d * 848 e, -(c - e * e), bounds); 849 850 /* it can still be empty because of our infinity, I guess... */ 851 sol = SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, result) ? SCIPinfinity(scip) : SCIPintervalGetInf(result); 852 853 /* check that solution is acceptable, ideally it should be <= 0, however when it is positive, we trigger a binary 854 * search to make it negative. This binary search might return a solution point that is not at accurately 0 as the 855 * one obtained from the function above. Thus, it might fail to satisfy the condition of case 4b in some cases, e.g., 856 * ex8_3_1, bchoco05, etc 857 */ 858 if( evalPhiAtRay(scip, sol, a, b, c, d, e) <= 1e-10 ) 859 { 860 #ifdef INTERCUT_MOREDEBUG 861 printf("interval solution returned %g -> phival = %g, believe it\n", sol, evalPhiAtRay(sol, a, b, c, d, e)); 862 printf("don't do bin search\n"); 863 #endif 864 865 return sol; 866 } 867 else 868 { 869 /* perform a binary search to make it negative: this might correct a wrong infinity (e.g. crudeoil_lee1_05) */ 870 #ifdef INTERCUT_MOREDEBUG 871 printf("do bin search because phival is %g\n", evalPhiAtRay(scip, sol, a, b, c, d, e)); 872 #endif 873 doBinarySearch(scip, a, b, c, d, e, &sol); 874 } 875 876 return sol; 877 } 878 879 /** The maximal S-free set is gamma(z) <= 0; we find the intersection point of the ray `ray` starting from zlp with the 880 * boundary of the S-free set. 881 * That is, we find t >= 0 such that gamma(zlp + t * ray) = 0. 882 * 883 * In cases 1,2, and 3, gamma is of the form 884 * gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) 885 * 886 * In the case 4 gamma is of the form 887 * gamma(zlp + t * ray) = SQRT(A t^2 + B t + C) - (D t + E) if some condition holds 888 * SQRT(A' t^2 + B' t + C') - (D' t + E') otherwise 889 * 890 * It can be shown (given the special properties of gamma) that the smallest positive root of each function of the form 891 * SQRT(a t^2 + b t + c) - (d t + e) 892 * is the same as the smallest positive root of the quadratic equation: 893 * (SQRT(a t^2 + b t + c) - (d t + e)) * (SQRT(a t^2 + b t + c) + (d t + e)) = 0 894 * <==> (a - d^2) t^2 + (b - 2 d*e) t + (c - e^2) = 0 895 * 896 * So, in cases 1, 2, and 3, this function just returns the solution of the above equation. 897 * In case 4, it first solves the equation assuming we are in the first piece. 898 * If there is no solution, then the second piece can't have a solution (first piece >= second piece for all t) 899 * Then we check if the solution satisfies the condition. 900 * If it doesn't then we solve the equation for the second piece. 901 * If it has a solution, then it _has_ to be the solution. 902 */ 903 static 904 SCIP_Real computeIntersectionPoint( 905 SCIP* scip, /**< SCIP data structure */ 906 SCIP_Bool usebounds, /**< whether we are in case 4 or not */ 907 SCIP_Real* coefs, /**< values of A, B, C, D, and E of cases 1, 2, 3, or 4a */ 908 SCIP_Real* coefs4b, /**< values of A, B, C, D, and E of case 4b */ 909 SCIP_Real* coefscondition /**< values needed to evaluate condition of case 4 */ 910 ) 911 { 912 SCIP_Real sol; 913 SCIP_Real sol4b; 914 915 assert(coefs != NULL); 916 917 if( ! usebounds ) 918 return computeRoot(scip, coefs); 919 920 assert(coefs4b != NULL); 921 assert(coefscondition != NULL); 922 923 /* compute solution of first piece */ 924 sol = computeRoot(scip, coefs); 925 926 /* if there is no solution --> second piece doesn't have solution */ 927 if( SCIPisInfinity(scip, sol) ) 928 { 929 /* this assert fails on multiplants_mtg5 the problem is that sqrt(A) <= D in 4a but not in 4b, 930 * now, this is impossible since the phi4a >= phi4b, so actually sqrt(A) is 10e-15 away from 931 * D in 4b 932 */ 933 /* assert(SCIPisInfinity(scip, computeRoot(scip, coefs4b))); */ 934 return sol; 935 } 936 937 /* if solution of 4a is in 4a, then return */ 938 if( isCase4a(sol, coefs, coefscondition) ) 939 return sol; 940 941 /* not on 4a --> then the intersection point is whatever 4b says: as phi4a >= phi4b, the solution of phi4b should 942 * always be larger (but shouldn't be equal at this point given that isCase4a failed, and the condition function 943 * evaluates to 0 when phi4a == phi4b) than the solution of phi4a; However, because of numerics (or limits in the 944 * binary search) we can find a slightly smaller solution; thus, we just keep the larger one 945 */ 946 sol4b = computeRoot(scip, coefs4b); 947 948 return MAX(sol, sol4b); 949 } 950 951 /** adds cutcoef * (col - col*) to rowprep */ 952 static 953 SCIP_RETCODE addColToCut( 954 SCIP* scip, /**< SCIP data structure */ 955 SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */ 956 SCIP_Real cutcoef, /**< cut coefficient */ 957 SCIP_COL* col /**< column to add to rowprep */ 958 ) 959 { 960 assert(col != NULL); 961 962 #ifdef DEBUG_INTERCUTS_NUMERICS 963 SCIPinfoMessage(scip, NULL, "adding col %s to cut. %g <= col <= %g\n", SCIPvarGetName(SCIPcolGetVar(col)), 964 SCIPvarGetLbLocal(SCIPcolGetVar(col)), SCIPvarGetUbLocal(SCIPcolGetVar(col))); 965 SCIPinfoMessage(scip, NULL, "col is active at %s. Value %.15f\n", SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ? "lower bound" : 966 "upper bound" , SCIPcolGetPrimsol(col)); 967 #endif 968 969 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(col), cutcoef) ); 970 SCIProwprepAddConstant(rowprep, -cutcoef * SCIPcolGetPrimsol(col) ); 971 972 return SCIP_OKAY; 973 } 974 975 /** adds cutcoef * (slack - slack*) to rowprep 976 * 977 * row is lhs <= <coefs, vars> + constant <= rhs, thus slack is defined by 978 * slack + <coefs, vars> + constant = side 979 * If row (slack) is at upper, it means that <coefs,vars*> + constant = rhs, and so 980 * slack* = side - rhs --> slack - slack* = rhs - <coefs, vars> - constant. 981 * If row (slack) is at lower, then <coefs,vars*> + constant = lhs, and so 982 * slack* = side - lhs --> slack - slack* = lhs - <coefs, vars> - constant. 983 */ 984 static 985 SCIP_RETCODE addRowToCut( 986 SCIP* scip, /**< SCIP data structure */ 987 SCIP_ROWPREP* rowprep, /**< rowprep to store intersection cut */ 988 SCIP_Real cutcoef, /**< cut coefficient */ 989 SCIP_ROW* row, /**< row, whose slack we are adding to rowprep */ 990 SCIP_Bool* success /**< buffer to store whether the row is nonbasic enough */ 991 ) 992 { 993 int i; 994 SCIP_COL** rowcols; 995 SCIP_Real* rowcoefs; 996 int nnonz; 997 998 assert(row != NULL); 999 1000 rowcols = SCIProwGetCols(row); 1001 rowcoefs = SCIProwGetVals(row); 1002 nnonz = SCIProwGetNLPNonz(row); 1003 1004 #ifdef DEBUG_INTERCUTS_NUMERICS 1005 SCIPinfoMessage(scip, NULL, "adding slack var row_%d to cut. %g <= row <= %g\n", SCIProwGetLPPos(row), SCIProwGetLhs(row), SCIProwGetRhs(row)); 1006 SCIPinfoMessage(scip, NULL, "row is active at %s = %.15f Activity %.15f\n", SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? "lhs" : 1007 "rhs" , SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ? SCIProwGetLhs(row) : SCIProwGetRhs(row), 1008 SCIPgetRowActivity(scip, row)); 1009 #endif 1010 1011 if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ) 1012 { 1013 assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row))); 1014 if( ! SCIPisFeasEQ(scip, SCIProwGetLhs(row), SCIPgetRowActivity(scip, row)) ) 1015 { 1016 *success = FALSE; 1017 return SCIP_OKAY; 1018 } 1019 1020 SCIProwprepAddConstant(rowprep, SCIProwGetLhs(row) * cutcoef); 1021 } 1022 else 1023 { 1024 assert(!SCIPisInfinity(scip, SCIProwGetRhs(row))); 1025 if( ! SCIPisFeasEQ(scip, SCIProwGetRhs(row), SCIPgetRowActivity(scip, row)) ) 1026 { 1027 *success = FALSE; 1028 return SCIP_OKAY; 1029 } 1030 1031 SCIProwprepAddConstant(rowprep, SCIProwGetRhs(row) * cutcoef); 1032 } 1033 1034 for( i = 0; i < nnonz; i++ ) 1035 { 1036 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPcolGetVar(rowcols[i]), -rowcoefs[i] * cutcoef) ); 1037 } 1038 1039 SCIProwprepAddConstant(rowprep, -SCIProwGetConstant(row) * cutcoef); 1040 1041 return SCIP_OKAY; 1042 } 1043 1044 /** get the tableau rows of the variables in vars */ 1045 static 1046 SCIP_RETCODE getTableauRows( 1047 SCIP* scip, /**< SCIP data structure */ 1048 SCIP_VAR** vars, /**< variables in the minor */ 1049 int* basicvarpos2tableaurow,/**< map between basic var and its tableau row */ 1050 SCIP_HASHMAP* tableau, /**< map between var an its tableau row */ 1051 SCIP_Real** tableaurows, /**< buffer to store tableau row */ 1052 SCIP_Bool* success /**< set to TRUE if no variable had basisstat = ZERO */ 1053 ) 1054 { 1055 int v; 1056 int nrows; 1057 int ncols; 1058 1059 *success = TRUE; 1060 1061 nrows = SCIPgetNLPRows(scip); 1062 ncols = SCIPgetNLPCols(scip); 1063 1064 /* check if we have the tableau row of the variable and if not compute it */ 1065 for( v = 0; v < 4; ++v ) 1066 { 1067 if( ! SCIPhashmapExists(tableau, (void*)vars[v]) ) 1068 { 1069 SCIP_COL* col; 1070 1071 /* get column of variable */ 1072 col = SCIPvarGetCol(vars[v]); 1073 1074 /* if variable is basic, then get its tableau row and insert it in the hashmap */ 1075 if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC ) 1076 { 1077 int lppos; 1078 SCIP_Real* densetableaurow; 1079 1080 lppos = SCIPcolGetLPPos(col); 1081 SCIP_CALL( SCIPallocBufferArray(scip, &densetableaurow, ncols + nrows) ); 1082 1083 SCIP_CALL( SCIPgetLPBInvRow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], NULL, NULL) ); 1084 SCIP_CALL( SCIPgetLPBInvARow(scip, basicvarpos2tableaurow[lppos], &densetableaurow[ncols], densetableaurow, NULL, NULL) ); 1085 1086 /* insert tableau row in hashmap*/ 1087 SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)densetableaurow) ); 1088 } 1089 else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO ) 1090 { 1091 *success = FALSE; 1092 return SCIP_OKAY; /* don't even bother */ 1093 } 1094 else 1095 { 1096 SCIP_CALL( SCIPhashmapInsert(tableau, (void*)vars[v], (void *)NULL) ); 1097 } 1098 } 1099 1100 /* get tableau row of var */ 1101 tableaurows[v] = (SCIP_Real *)SCIPhashmapGetImage(tableau, (void*)vars[v]); 1102 } 1103 return SCIP_OKAY; 1104 } 1105 1106 /** computes the cut coefs of the non-basic (non-slack) variables (correspond to cols) and adds them to the 1107 * intersection cut 1108 */ 1109 static 1110 SCIP_RETCODE addCols( 1111 SCIP* scip, /**< SCIP data structure */ 1112 SCIP_VAR** vars, /**< variables */ 1113 SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */ 1114 SCIP_ROWPREP* rowprep, /**< store cut */ 1115 SCIP_Real* rays, /**< buffer to store rays */ 1116 int* nrays, /**< pointer to store number of nonzero rays */ 1117 int* rayslppos, /**< buffer to store lppos of nonzero rays */ 1118 SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */ 1119 SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */ 1120 SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */ 1121 SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */ 1122 ) 1123 { 1124 int i; 1125 int ncols; 1126 SCIP_COL** cols; 1127 1128 *success = TRUE; 1129 1130 /* loop over non-basic (non-slack) variables */ 1131 cols = SCIPgetLPCols(scip); 1132 ncols = SCIPgetNLPCols(scip); 1133 for( i = 0; i < ncols; ++i ) 1134 { 1135 SCIP_COL* col; 1136 SCIP_Real coefs[5]; 1137 SCIP_Real coefs4b[5]; 1138 SCIP_Real coefscondition[3]; 1139 SCIP_Real factor; 1140 SCIP_Bool israynonzero; 1141 SCIP_Real cutcoef; 1142 SCIP_Real interpoint; 1143 int v; 1144 1145 col = cols[i]; 1146 1147 /* set factor to store entries of ray as = [-BinvL, BinvU] */ 1148 if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ) 1149 factor = -1.0; 1150 else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ) 1151 factor = 1.0; 1152 else if( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_ZERO ) 1153 { 1154 *success = FALSE; 1155 return SCIP_OKAY; 1156 } 1157 else 1158 continue; 1159 1160 /* build the ray */ 1161 israynonzero = FALSE; 1162 for( v = 0; v < 4; ++v ) 1163 { 1164 if( tableaurows[v] != NULL ) 1165 rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][i]) ? 0.0 : tableaurows[v][i]); 1166 else 1167 { 1168 if( col == SCIPvarGetCol(vars[v]) ) 1169 rays[(*nrays) * 4 + v] = -factor; 1170 else 1171 rays[(*nrays) * 4 + v] = 0.0; 1172 } 1173 1174 israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0); 1175 } 1176 1177 /* do nothing if ray is 0 */ 1178 if( ! israynonzero ) 1179 continue; 1180 1181 /* compute the cut */ 1182 SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds, 1183 ad, success) ); 1184 1185 if( *success == FALSE ) 1186 return SCIP_OKAY; 1187 1188 /* compute intersection point */ 1189 interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition); 1190 1191 /* store intersection points */ 1192 interpoints[*nrays] = interpoint; 1193 1194 /* remember lppos */ 1195 rayslppos[*nrays] = i; 1196 1197 /* count nonzero rays */ 1198 *nrays += 1; 1199 1200 /* compute cut coef */ 1201 cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint; 1202 1203 /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */ 1204 assert(SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER); 1205 SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ? -cutcoef : 1206 cutcoef, col) ); 1207 } 1208 1209 return SCIP_OKAY; 1210 } 1211 1212 /** computes the cut coefs of the non-basic slack variables (correspond to rows) and adds them to the 1213 * intersection cut 1214 */ 1215 static 1216 SCIP_RETCODE addRows( 1217 SCIP* scip, /**< SCIP data structure */ 1218 SCIP_VAR** vars, /**< variables */ 1219 SCIP_Real** tableaurows, /**< tableau rows corresponding to the variables in vars */ 1220 SCIP_ROWPREP* rowprep, /**< store cut */ 1221 SCIP_Real* rays, /**< buffer to store rays */ 1222 int* nrays, /**< pointer to store number of nonzero rays */ 1223 int* rayslppos, /**< buffer to store lppos of nonzero rays */ 1224 SCIP_Real* interpoints, /**< buffer to store intersection points or NULL if not needed */ 1225 SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */ 1226 SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */ 1227 SCIP_Bool* success /**< pointer to store whether the generation of cutcoefs was successful */ 1228 ) 1229 { 1230 int i; 1231 int nrows; 1232 int ncols; 1233 SCIP_ROW** rows; 1234 1235 nrows = SCIPgetNLPRows(scip); 1236 ncols = SCIPgetNLPCols(scip); 1237 1238 *success = TRUE; 1239 1240 /* loop over non-basic slack variables */ 1241 rows = SCIPgetLPRows(scip); 1242 for( i = 0; i < nrows; ++i ) 1243 { 1244 SCIP_ROW* row; 1245 SCIP_Real coefs[5]; 1246 SCIP_Real coefs4b[5]; 1247 SCIP_Real coefscondition[3]; 1248 SCIP_Real factor; 1249 SCIP_Bool israynonzero; 1250 SCIP_Real cutcoef; 1251 SCIP_Real interpoint; 1252 int v; 1253 1254 row = rows[i]; 1255 1256 /* set factor to store entries of ray as = [BinvL, -BinvU] */ 1257 if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER ) 1258 factor = 1.0; 1259 else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ) 1260 factor = -1.0; 1261 else if( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_ZERO ) 1262 { 1263 *success = FALSE; 1264 return SCIP_OKAY; 1265 } 1266 else 1267 continue; 1268 1269 /* build the ray */ 1270 israynonzero = FALSE; 1271 for( v = 0; v < 4; ++v ) 1272 { 1273 int idx; 1274 1275 idx = ncols + i; 1276 1277 if( tableaurows[v] != NULL ) 1278 rays[(*nrays) * 4 + v] = factor * (SCIPisZero(scip, tableaurows[v][idx]) ? 0.0 : tableaurows[v][idx]); 1279 else 1280 { 1281 /* TODO: We assume that slack variables can never occure in the minor. This is correct, right? */ 1282 rays[(*nrays) * 4 + v] = 0.0; 1283 } 1284 1285 israynonzero = israynonzero || (rays[(*nrays) * 4 + v] != 0.0); 1286 } 1287 1288 /* do nothing if ray is 0 */ 1289 if( ! israynonzero ) 1290 continue; 1291 1292 /* compute the cut */ 1293 SCIP_CALL( computeRestrictionToRay(scip, &rays[(*nrays) * 4], vars, coefs, coefs4b, coefscondition, usebounds, 1294 ad, success) ); 1295 1296 if( *success == FALSE ) 1297 return SCIP_OKAY; 1298 1299 /* compute intersection point */ 1300 interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition); 1301 1302 /* store intersection points */ 1303 interpoints[*nrays] = interpoint; 1304 1305 /* store lppos of ray, make it negative so we can differentiate between cols and rows */ 1306 rayslppos[*nrays] = -i - 1; 1307 1308 /* count nonzero rays */ 1309 *nrays += 1; 1310 1311 /* compute cut coef */ 1312 cutcoef = SCIPisInfinity(scip, interpoint) ? 0.0 : 1.0 / interpoint; 1313 1314 /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */ 1315 assert(SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER); 1316 1317 SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ? cutcoef : 1318 -cutcoef, row, success) ); /* rows have flipper base status! */ 1319 } 1320 1321 return SCIP_OKAY; 1322 } 1323 1324 /* checks if two rays are linearly dependent */ 1325 static 1326 SCIP_Bool raysAreDependent( 1327 SCIP* scip, /**< SCIP data structure */ 1328 SCIP_Real* ray1, /**< coefficients of ray 1 */ 1329 SCIP_Real* ray2, /**< coefficients of ray 2 */ 1330 SCIP_Real* coef /**< pointer to store coef (s.t. r1 = coef * r2) in case rays are dependent */ 1331 ) 1332 { 1333 int i; 1334 1335 *coef = 0.0; 1336 1337 for( i = 0; i < 4; ++i ) 1338 { 1339 /* rays cannot be dependent if one ray has zero entry and the other one doesn't */ 1340 if( (SCIPisZero(scip, ray1[i]) && ! SCIPisZero(scip, ray2[i])) || 1341 (! SCIPisZero(scip, ray1[i]) && SCIPisZero(scip, ray2[i])) ) 1342 { 1343 return FALSE; 1344 } 1345 1346 if( *coef != 0.0 ) 1347 { 1348 /* cannot be dependent if the coefs aren't equal for all entries */ 1349 if( ! SCIPisFeasEQ(scip, *coef, ray1[i] / ray2[i]) ) 1350 return FALSE; 1351 } 1352 else 1353 *coef = ray1[i] / ray2[i]; 1354 } 1355 1356 return TRUE; 1357 } 1358 1359 /** finds the smallest negative steplength for the current ray r_idx such that the combination 1360 * of r_idx with all rays not in the recession cone is in the recession cone 1361 */ 1362 static 1363 SCIP_RETCODE findRho( 1364 SCIP* scip, /**< SCIP data structure */ 1365 SCIP_Real* rays, /**< rays */ 1366 int nrays, /**< number of nonzero rays */ 1367 int idx, /**< index of current ray we want to find rho for */ 1368 SCIP_Real* interpoints, /**< intersection points of nonzero rays */ 1369 SCIP_VAR** vars, /**< variables */ 1370 SCIP_Real* rho, /**< pointer to store the optimal rho */ 1371 SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */ 1372 SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */ 1373 SCIP_Bool* success /**< TRUE if computation of rho was successful */ 1374 ) 1375 { 1376 int i; 1377 1378 *success = TRUE; 1379 1380 /* go through all rays not in the recession cone and compute the largest negative steplength possible. The 1381 * smallest of them is then the steplength rho we use for the current ray */ 1382 *rho = 0; 1383 for( i = 0; i < nrays; ++i ) 1384 { 1385 SCIP_Real currentrho; 1386 SCIP_Real coef; 1387 1388 if( SCIPisInfinity(scip, interpoints[i]) ) 1389 continue; 1390 1391 /* if the rays are linearly independent, we don't need to search for rho */ 1392 if( raysAreDependent(scip, &rays[4 * i], &rays[4 * idx], &coef) ) 1393 currentrho = coef * interpoints[i]; 1394 else 1395 { 1396 SCIP_Real lb; 1397 SCIP_Real ub; 1398 SCIP_Real alpha; 1399 int j; 1400 1401 /* do binary search by lookig at the convex combinations of r_i and r_j */ 1402 lb = 0.0; 1403 ub = 1.0; 1404 1405 for( j = 0; j < BINSEARCH_MAXITERS; ++j ) 1406 { 1407 SCIP_Real coefs[5]; 1408 SCIP_Real coefs4b[5]; 1409 SCIP_Real coefscondition[3]; 1410 SCIP_Real newray[4]; 1411 SCIP_Real interpoint; 1412 int k; 1413 1414 alpha = (lb + ub) / 2.0; 1415 1416 /* build the ray alpha * ray_i + (1 - alpha) * ray_idx */ 1417 for( k = 0; k < 4; ++k ) 1418 newray[k] = alpha * rays[4 * i + k] - (1 - alpha) * rays[4 * idx + k]; 1419 1420 /* restrict phi to the "new" ray */ 1421 SCIP_CALL( computeRestrictionToRay(scip, newray, vars, coefs, coefs4b, coefscondition, usebounds, 1422 ad, success) ); 1423 1424 if( ! *success ) 1425 return SCIP_OKAY; 1426 1427 /* check if restriction to "new" ray is numerically nasty. If so, treat the corresponding rho as if phi is 1428 * positive 1429 */ 1430 1431 /* compute intersection point */ 1432 interpoint = computeIntersectionPoint(scip, usebounds, coefs, coefs4b, coefscondition); 1433 1434 /* no root exists */ 1435 if( SCIPisInfinity(scip, interpoint) ) 1436 { 1437 lb = alpha; 1438 if( SCIPisEQ(scip, ub, lb) ) 1439 break; 1440 } 1441 else 1442 ub = alpha; 1443 } 1444 1445 /* now we found the best convex combination which we use to derive the corresponding coef. If alpha = 0, we 1446 * cannot move the ray in the recession cone, i.e. strengthening is not possible */ 1447 if( SCIPisZero(scip, alpha) ) 1448 { 1449 *rho = -SCIPinfinity(scip); 1450 return SCIP_OKAY; 1451 } 1452 else 1453 currentrho = (alpha - 1) * interpoints[i] / alpha; 1454 } 1455 1456 if( currentrho < *rho ) 1457 *rho = currentrho; 1458 } 1459 1460 return SCIP_OKAY; 1461 } 1462 1463 /** computes negative steplengths for the rays that are in the recession cone of the S-free set, i.e., 1464 * which have an infinite intersection point. 1465 */ 1466 static 1467 SCIP_RETCODE computeNegCutcoefs( 1468 SCIP* scip, /**< SCIP data structure */ 1469 SCIP_VAR** vars, /**< variables */ 1470 SCIP_Real* rays, /**< rays */ 1471 int nrays, /**< number of nonzero rays */ 1472 int* rayslppos, /**< lppos of nonzero rays */ 1473 SCIP_Real* interpoints, /**< intersection points */ 1474 SCIP_ROWPREP* rowprep, /**< rowprep for the generated cut */ 1475 SCIP_Bool usebounds, /**< TRUE if we want to separate non-negative bound */ 1476 SCIP_Real* ad, /**< coefs a and d for the hyperplane aTx + dTy <= 0 */ 1477 SCIP_Bool* success /**< if a cut candidate could be computed */ 1478 ) 1479 { 1480 SCIP_COL** cols; 1481 SCIP_ROW** rows; 1482 int i; 1483 1484 *success = TRUE; 1485 1486 cols = SCIPgetLPCols(scip); 1487 rows = SCIPgetLPRows(scip); 1488 1489 /* go through all intersection points that are equal to infinity -> these correspond to the rays which are in the 1490 * recession cone of C, i.e. the rays for which we (possibly) can compute a negative steplength */ 1491 for( i = 0; i < nrays ; ++i ) 1492 { 1493 SCIP_Real rho; 1494 SCIP_Real cutcoef; 1495 int lppos; 1496 1497 if( !SCIPisInfinity(scip, interpoints[i]) ) 1498 continue; 1499 1500 /* compute the smallest rho */ 1501 SCIP_CALL( findRho(scip, rays, nrays, i, interpoints, vars, &rho, usebounds, ad, success) ); 1502 1503 if( ! *success ) 1504 continue; 1505 1506 /* compute cut coef */ 1507 cutcoef = SCIPisInfinity(scip, -rho) ? 0.0 : 1.0 / rho; 1508 1509 /* add var to cut: if variable is nonbasic at upper we have to flip sign of cutcoef */ 1510 lppos = rayslppos[i]; 1511 if( lppos < 0 ) 1512 { 1513 lppos = -lppos - 1; 1514 1515 assert(SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(rows[lppos]) == 1516 SCIP_BASESTAT_UPPER); 1517 1518 SCIP_CALL( addRowToCut(scip, rowprep, SCIProwGetBasisStatus(rows[lppos]) == SCIP_BASESTAT_UPPER ? cutcoef : 1519 -cutcoef, rows[lppos], success) ); /* rows have flipped base status! */ 1520 1521 if( ! *success ) 1522 return SCIP_OKAY; 1523 } 1524 else 1525 { 1526 assert(SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER || SCIPcolGetBasisStatus(cols[lppos]) == 1527 SCIP_BASESTAT_LOWER); 1528 SCIP_CALL( addColToCut(scip, rowprep, SCIPcolGetBasisStatus(cols[lppos]) == SCIP_BASESTAT_UPPER ? -cutcoef : 1529 cutcoef, cols[lppos]) ); 1530 } 1531 } 1532 1533 return SCIP_OKAY; 1534 } 1535 1536 /** separates cuts for stored principal minors */ 1537 static 1538 SCIP_RETCODE separateDeterminant( 1539 SCIP* scip, /**< SCIP data structure */ 1540 SCIP_SEPA* sepa, /**< separator */ 1541 SCIP_SEPADATA* sepadata, /**< separator data */ 1542 SCIP_VAR* xik, /**< variable X_ik = x_i * x_k */ 1543 SCIP_VAR* xil, /**< variable X_il = x_i * x_l */ 1544 SCIP_VAR* xjk, /**< variable X_jk = x_j * x_k */ 1545 SCIP_VAR* xjl, /**< variable X_jl = x_j * x_l */ 1546 SCIP_Bool* isxikdiag, /**< is X_ik diagonal? (i.e. i = k) */ 1547 SCIP_Bool* isxildiag, /**< is X_il diagonal? (i.e. i = l) */ 1548 SCIP_Bool* isxjkdiag, /**< is X_jk diagonal? (i.e. j = k) */ 1549 SCIP_Bool* isxjldiag, /**< is X_jl diagonal? (i.e. j = l) */ 1550 int* basicvarpos2tableaurow,/**< map from basic var to its tableau row */ 1551 SCIP_HASHMAP* tableau, /**< map from var to its tableau row */ 1552 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 1553 ) 1554 { 1555 SCIP_ROWPREP* rowprep; 1556 SCIP_VAR* vars[4] = {xik, xjl, xil, xjk}; 1557 SCIP_Real* tableaurows[4]; 1558 SCIP_Real* interpoints; 1559 SCIP_Real* rays; 1560 int nrays; 1561 int* rayslppos; 1562 int ncols; 1563 int nrows; 1564 SCIP_Bool success; 1565 SCIP_Real ad[4] = {0.0, 0.0, 0.0, 0.0}; 1566 SCIP_Real solxik; 1567 SCIP_Real solxil; 1568 SCIP_Real solxjk; 1569 SCIP_Real solxjl; 1570 1571 ncols = SCIPgetNLPCols(scip); 1572 nrows = SCIPgetNLPRows(scip); 1573 1574 /* allocate memory for intersection points */ 1575 SCIP_CALL( SCIPallocBufferArray(scip, &interpoints, ncols + nrows) ); 1576 1577 /* allocate memory for rays */ 1578 SCIP_CALL( SCIPallocBufferArray(scip, &rays, 4 * (ncols + nrows)) ); 1579 SCIP_CALL( SCIPallocBufferArray(scip, &rayslppos, ncols + nrows) ); 1580 1581 /* cut (in the nonbasic space) is of the form alpha^T x >= 1 */ 1582 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, SCIP_SIDETYPE_LEFT, TRUE) ); 1583 SCIProwprepAddSide(rowprep, 1.0); 1584 1585 /* check if we have the tableau row of the variable and if not compute it */ 1586 SCIP_CALL( getTableauRows(scip, vars, basicvarpos2tableaurow, tableau, tableaurows, &success) ); 1587 1588 if( ! success ) 1589 goto CLEANUP; 1590 1591 /* if we want to enforce bounds, set the right a and d to enforce aTx + dTy <= 0 */ 1592 if( sepadata->usebounds ) 1593 { 1594 solxik = SCIPvarGetLPSol(xik); 1595 solxil = SCIPvarGetLPSol(xil); 1596 solxjk = SCIPvarGetLPSol(xjk); 1597 solxjl = SCIPvarGetLPSol(xjl); 1598 1599 if( isxikdiag && SCIPisFeasNegative(scip, solxik) ) 1600 { 1601 ad[0] = -1.0; 1602 ad[2] = 1.0; 1603 } 1604 else if( isxjldiag && SCIPisFeasNegative(scip, solxjl) ) 1605 { 1606 ad[0] = -1.0; 1607 ad[2] = -1.0; 1608 } 1609 else if( isxildiag && SCIPisFeasNegative(scip, solxil) ) 1610 { 1611 ad[1] = 1.0; 1612 ad[3] = -1.0; 1613 } 1614 else if( isxjkdiag && SCIPisFeasNegative(scip, solxjk) ) 1615 { 1616 ad[1] = -1.0; 1617 ad[3] = -1.0; 1618 } 1619 } 1620 1621 nrays = 0; 1622 /* loop over each non-basic var; get the ray; compute cut coefficient */ 1623 SCIP_CALL( addCols(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) ); 1624 1625 if( ! success ) 1626 goto CLEANUP; 1627 1628 /* loop over non-basic slack variables */ 1629 SCIP_CALL( addRows(scip, vars, tableaurows, rowprep, rays, &nrays, rayslppos, interpoints, sepadata->usebounds, ad, &success) ); 1630 1631 if( ! success ) 1632 goto CLEANUP; 1633 1634 /* do strengthening */ 1635 if( sepadata->usestrengthening ) 1636 { 1637 SCIP_CALL( computeNegCutcoefs(scip, vars, rays, nrays, rayslppos, interpoints, rowprep, sepadata->usebounds, ad, &success) ); 1638 1639 if( ! success ) 1640 goto CLEANUP; 1641 } 1642 1643 /* merge coefficients that belong to same variable */ 1644 SCIPmergeRowprepTerms(scip, rowprep); 1645 1646 SCIP_CALL( SCIPcleanupRowprep(scip, rowprep, NULL, sepadata->mincutviol, NULL, &success) ); 1647 1648 /* if cleanup was successfull, create row out of rowprep and add it */ 1649 if( success ) 1650 { 1651 SCIP_ROW* row; 1652 SCIP_Bool infeasible; 1653 1654 /* create row */ 1655 SCIP_CALL( SCIPgetRowprepRowSepa(scip, &row, rowprep, sepa) ); 1656 1657 assert(SCIPgetCutEfficacy(scip, NULL, row) > 0.0); 1658 1659 /* add row */ 1660 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) ); 1661 1662 if( infeasible ) 1663 *result = SCIP_CUTOFF; 1664 else 1665 *result = SCIP_SEPARATED; 1666 1667 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 1668 } 1669 1670 CLEANUP: 1671 SCIPfreeRowprep(scip, &rowprep); 1672 SCIPfreeBuffer(scip, &rayslppos); 1673 SCIPfreeBuffer(scip, &rays); 1674 SCIPfreeBuffer(scip, &interpoints); 1675 1676 return SCIP_OKAY; 1677 } 1678 1679 1680 /** separates cuts for stored principal minors */ 1681 static 1682 SCIP_RETCODE separatePoint( 1683 SCIP* scip, /**< SCIP data structure */ 1684 SCIP_SEPA* sepa, /**< separator */ 1685 SCIP_RESULT* result /**< pointer to store the result of the separation call */ 1686 ) 1687 { 1688 SCIP_SEPADATA* sepadata; 1689 SCIP_HASHMAP* tableau = NULL; 1690 int* basicvarpos2tableaurow = NULL; /* map between basic var and its tableau row */ 1691 int i; 1692 1693 assert(sepa != NULL); 1694 assert(result != NULL); 1695 1696 *result = SCIP_DIDNOTRUN; 1697 1698 sepadata = SCIPsepaGetData(sepa); 1699 assert(sepadata != NULL); 1700 1701 /* check whether there are some minors available */ 1702 if( sepadata->nminors == 0 ) 1703 return SCIP_OKAY; 1704 1705 *result = SCIP_DIDNOTFIND; 1706 1707 /* loop over the minors and if they are violated build cut */ 1708 for( i = 0; i < sepadata->nminors && (*result != SCIP_CUTOFF); ++i ) 1709 { 1710 SCIP_VAR* auxvarxik; 1711 SCIP_VAR* auxvarxil; 1712 SCIP_VAR* auxvarxjk; 1713 SCIP_VAR* auxvarxjl; 1714 SCIP_Bool isauxvarxikdiag; 1715 SCIP_Bool isauxvarxildiag; 1716 SCIP_Bool isauxvarxjkdiag; 1717 SCIP_Bool isauxvarxjldiag; 1718 SCIP_Real solxik; 1719 SCIP_Real solxil; 1720 SCIP_Real solxjk; 1721 SCIP_Real solxjl; 1722 SCIP_Real det; 1723 1724 /* get variables of the i-th minor */ 1725 SCIP_CALL( getMinorVars(sepadata, i, &auxvarxik, &auxvarxil, &auxvarxjk, &auxvarxjl, &isauxvarxikdiag, 1726 &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag) ); 1727 1728 /* get current solution values */ 1729 solxik = SCIPvarGetLPSol(auxvarxik); 1730 solxil = SCIPvarGetLPSol(auxvarxil); 1731 solxjk = SCIPvarGetLPSol(auxvarxjk); 1732 solxjl = SCIPvarGetLPSol(auxvarxjl); 1733 1734 det = solxik * solxjl - solxil * solxjk; 1735 1736 if( SCIPisFeasZero(scip, det) ) 1737 continue; 1738 1739 if( basicvarpos2tableaurow == NULL ) 1740 { 1741 /* allocate memory */ 1742 SCIP_CALL( SCIPallocBufferArray(scip, &basicvarpos2tableaurow, SCIPgetNLPCols(scip)) ); 1743 SCIP_CALL( SCIPhashmapCreate(&tableau, SCIPblkmem(scip), SCIPgetNVars(scip)) ); 1744 1745 /* construct basicvar to tableau row map */ 1746 SCIP_CALL( constructBasicVars2TableauRowMap(scip, basicvarpos2tableaurow) ); 1747 } 1748 assert(tableau != NULL); 1749 1750 if( SCIPisFeasPositive(scip, det) ) 1751 { 1752 SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxik, auxvarxil, auxvarxjk, auxvarxjl, &isauxvarxikdiag, 1753 &isauxvarxildiag, &isauxvarxjkdiag, &isauxvarxjldiag, basicvarpos2tableaurow, tableau, result) ); 1754 } 1755 else 1756 { 1757 assert(SCIPisFeasNegative(scip, det)); 1758 SCIP_CALL( separateDeterminant(scip, sepa, sepadata, auxvarxil, auxvarxik, auxvarxjl, auxvarxjk, &isauxvarxildiag, 1759 &isauxvarxikdiag, &isauxvarxjldiag, &isauxvarxjkdiag, basicvarpos2tableaurow, tableau, result) ); 1760 } 1761 } 1762 1763 /* all minors were feasible, so no memory to free */ 1764 if( basicvarpos2tableaurow == NULL ) 1765 return SCIP_OKAY; 1766 1767 /* free memory */ 1768 for( i = 0; i < SCIPhashmapGetNEntries(tableau); ++i ) 1769 { 1770 SCIP_HASHMAPENTRY* entry; 1771 1772 entry = SCIPhashmapGetEntry(tableau, i); 1773 1774 if( entry != NULL ) 1775 { 1776 SCIP_Real* tableaurow; 1777 1778 tableaurow = (SCIP_Real *) SCIPhashmapEntryGetImage(entry); 1779 1780 SCIPfreeBufferArrayNull(scip, &tableaurow); 1781 } 1782 } 1783 SCIPhashmapFree(&tableau); 1784 SCIPfreeBufferArray(scip, &basicvarpos2tableaurow); 1785 1786 return SCIP_OKAY; 1787 } 1788 1789 /* 1790 * Callback methods of separator 1791 */ 1792 1793 /** copy method for separator plugins (called when SCIP copies plugins) */ 1794 static 1795 SCIP_DECL_SEPACOPY(sepaCopyMinor) 1796 { /*lint --e{715}*/ 1797 assert(scip != NULL); 1798 assert(sepa != NULL); 1799 assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0); 1800 1801 /* call inclusion method of constraint handler */ 1802 SCIP_CALL( SCIPincludeSepaInterminor(scip) ); 1803 1804 return SCIP_OKAY; 1805 } 1806 1807 1808 /** destructor of separator to free user data (called when SCIP is exiting) */ 1809 static 1810 SCIP_DECL_SEPAFREE(sepaFreeMinor) 1811 { /*lint --e{715}*/ 1812 SCIP_SEPADATA* sepadata; 1813 1814 sepadata = SCIPsepaGetData(sepa); 1815 assert(sepadata != NULL); 1816 assert(sepadata->minors == NULL); 1817 assert(sepadata->nminors == 0); 1818 assert(sepadata->minorssize == 0); 1819 1820 /* free separator data */ 1821 SCIPfreeBlockMemory(scip, &sepadata); 1822 SCIPsepaSetData(sepa, NULL); 1823 1824 return SCIP_OKAY; 1825 } 1826 1827 1828 /** initialization method of separator (called after problem was transformed) */ 1829 static 1830 SCIP_DECL_SEPAINIT(sepaInitMinor) 1831 { /*lint --e{715}*/ 1832 SCIP_SEPADATA* sepadata; 1833 1834 /* get separator data */ 1835 sepadata = SCIPsepaGetData(sepa); 1836 assert(sepadata != NULL); 1837 assert(sepadata->randnumgen == NULL); 1838 1839 /* create random number generator */ 1840 SCIP_CALL( SCIPcreateRandom(scip, &sepadata->randnumgen, DEFAULT_RANDSEED, TRUE) ); 1841 1842 return SCIP_OKAY; 1843 } 1844 1845 1846 /** deinitialization method of separator (called before transformed problem is freed) */ 1847 static 1848 SCIP_DECL_SEPAEXIT(sepaExitMinor) 1849 { /*lint --e{715}*/ 1850 SCIP_SEPADATA* sepadata; 1851 1852 /* get separator data */ 1853 sepadata = SCIPsepaGetData(sepa); 1854 assert(sepadata != NULL); 1855 assert(sepadata->randnumgen != NULL); 1856 1857 /* free random number generator */ 1858 SCIPfreeRandom(scip, &sepadata->randnumgen); 1859 1860 return SCIP_OKAY; 1861 } 1862 1863 1864 /** solving process initialization method of separator (called when branch and bound process is about to begin) */ 1865 static 1866 SCIP_DECL_SEPAINITSOL(sepaInitsolMinor) 1867 { /*lint --e{715}*/ 1868 return SCIP_OKAY; 1869 } 1870 1871 1872 /** solving process deinitialization method of separator (called before branch and bound process data is freed) */ 1873 static 1874 SCIP_DECL_SEPAEXITSOL(sepaExitsolMinor) 1875 { /*lint --e{715}*/ 1876 SCIP_SEPADATA* sepadata; 1877 1878 sepadata = SCIPsepaGetData(sepa); 1879 assert(sepadata != NULL); 1880 1881 /* clear separation data */ 1882 SCIP_CALL( sepadataClear(scip, sepadata) ); 1883 1884 return SCIP_OKAY; 1885 } 1886 1887 1888 /** LP solution separation method of separator */ 1889 static 1890 SCIP_DECL_SEPAEXECLP(sepaExeclpMinor) 1891 { /*lint --e{715}*/ 1892 SCIP_SEPADATA* sepadata; 1893 int ncalls; 1894 int currentdepth; 1895 1896 /* need routine to compute eigenvalues/eigenvectors */ 1897 if( !SCIPisIpoptAvailableIpopt() ) 1898 return SCIP_OKAY; 1899 1900 sepadata = SCIPsepaGetData(sepa); 1901 assert(sepadata != NULL); 1902 currentdepth = SCIPgetDepth(scip); 1903 ncalls = SCIPsepaGetNCallsAtNode(sepa); 1904 1905 /* only call the separator a given number of times at each node */ 1906 if( (currentdepth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 1907 || (currentdepth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 1908 { 1909 SCIPdebugMsg(scip, "reached round limit for node\n"); 1910 return SCIP_OKAY; 1911 } 1912 1913 /* try to detect minors */ 1914 SCIP_CALL( detectMinors(scip, sepadata) ); 1915 1916 /* call separation method */ 1917 SCIP_CALL( separatePoint(scip, sepa, result) ); 1918 1919 return SCIP_OKAY; 1920 } 1921 1922 1923 /* 1924 * separator specific interface methods 1925 */ 1926 1927 /** creates the minor separator and includes it in SCIP */ 1928 SCIP_RETCODE SCIPincludeSepaInterminor( 1929 SCIP* scip /**< SCIP data structure */ 1930 ) 1931 { 1932 SCIP_SEPADATA* sepadata = NULL; 1933 SCIP_SEPA* sepa = NULL; 1934 1935 /* create minor separator data */ 1936 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) ); 1937 BMSclearMemory(sepadata); 1938 1939 /* include separator */ 1940 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 1941 SEPA_USESSUBSCIP, SEPA_DELAY, 1942 sepaExeclpMinor, NULL, 1943 sepadata) ); 1944 1945 assert(sepa != NULL); 1946 1947 /* set non fundamental callbacks via setter functions */ 1948 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyMinor) ); 1949 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeMinor) ); 1950 SCIP_CALL( SCIPsetSepaInit(scip, sepa, sepaInitMinor) ); 1951 SCIP_CALL( SCIPsetSepaExit(scip, sepa, sepaExitMinor) ); 1952 SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolMinor) ); 1953 SCIP_CALL( SCIPsetSepaExitsol(scip, sepa, sepaExitsolMinor) ); 1954 1955 /* add minor separator parameters */ 1956 SCIP_CALL( SCIPaddBoolParam(scip, 1957 "separating/" SEPA_NAME "/usestrengthening", 1958 "whether to use strengthened intersection cuts to separate minors", 1959 &sepadata->usestrengthening, FALSE, DEFAULT_USESTRENGTHENING, NULL, NULL) ); 1960 1961 SCIP_CALL( SCIPaddBoolParam(scip, 1962 "separating/" SEPA_NAME "/usebounds", 1963 "whether to also enforce nonegativity bounds of principle minors", 1964 &sepadata->usebounds, FALSE, DEFAULT_USEBOUNDS, NULL, NULL) ); 1965 1966 SCIP_CALL( SCIPaddRealParam(scip, 1967 "separating/" SEPA_NAME "/mincutviol", 1968 "minimum required violation of a cut", 1969 &sepadata->mincutviol, FALSE, DEFAULT_MINCUTVIOL, 0.0, SCIP_REAL_MAX, NULL, NULL) ); 1970 1971 SCIP_CALL( SCIPaddIntParam(scip, 1972 "separating/" SEPA_NAME "/maxrounds", 1973 "maximal number of separation rounds per node (-1: unlimited)", 1974 &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) ); 1975 1976 SCIP_CALL( SCIPaddIntParam(scip, 1977 "separating/" SEPA_NAME "/maxroundsroot", 1978 "maximal number of separation rounds in the root node (-1: unlimited)", 1979 &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) ); 1980 1981 return SCIP_OKAY; 1982 } 1983