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_disjunctive.c 26 * @ingroup DEFPLUGINS_SEPA 27 * @brief disjunctive cut separator 28 * @author Tobias Fischer 29 * @author Marc Pfetsch 30 * 31 * We separate disjunctive cuts for two term disjunctions of the form \f$x_1 = 0 \vee x_2 = 0\f$. They can be generated 32 * directly from the simplex tableau. For further information, we refer to@n 33 * "A complementarity-based partitioning and disjunctive cut algorithm for mathematical programming problems with 34 * equilibrium constraints"@n 35 * Júdice, J.J., Sherali, H.D., Ribeiro, I.M., Faustino, A.M., Journal of Global Optimization 36(1), 89–114 (2006) 36 * 37 * Cut coefficients belonging to integer variables can be strengthened by the 'monoidal cut strengthening' procedure, see@n 38 * "Strengthening cuts for mixed integer programs"@n 39 * Egon Balas, Robert G. Jeroslow, European Journal of Operational Research, Volume 4, Issue 4, 1980, Pages 224-234 40 */ 41 42 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 43 44 #include "blockmemshell/memory.h" 45 #include "scip/cons_sos1.h" 46 #include "scip/pub_cons.h" 47 #include "scip/pub_lp.h" 48 #include "scip/pub_message.h" 49 #include "scip/pub_misc.h" 50 #include "scip/pub_misc_sort.h" 51 #include "scip/pub_sepa.h" 52 #include "scip/pub_var.h" 53 #include "scip/scip_cons.h" 54 #include "scip/scip_cut.h" 55 #include "scip/scip_general.h" 56 #include "scip/scip_lp.h" 57 #include "scip/scip_mem.h" 58 #include "scip/scip_message.h" 59 #include "scip/scip_numerics.h" 60 #include "scip/scip_param.h" 61 #include "scip/scip_sepa.h" 62 #include "scip/scip_sol.h" 63 #include "scip/scip_solvingstats.h" 64 #include "scip/scip_tree.h" 65 #include "scip/sepa_disjunctive.h" 66 #include <string.h> 67 68 69 #define SEPA_NAME "disjunctive" 70 #define SEPA_DESC "disjunctive cut separator" 71 #define SEPA_PRIORITY 10 /**< priority for separation */ 72 #define SEPA_FREQ 0 /**< frequency for separating cuts; zero means to separate only in the root node */ 73 #define SEPA_MAXBOUNDDIST 0.0 /**< maximal relative distance from the current node's dual bound to primal bound 74 * compared to best node's dual bound for applying separation.*/ 75 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */ 76 #define SEPA_DELAY TRUE /**< should separation method be delayed, if other separators found cuts? */ 77 78 #define DEFAULT_MAXRANK 20 /**< maximal rank of a cut that could not be scaled to integral coefficients (-1: unlimited) */ 79 #define DEFAULT_MAXRANKINTEGRAL -1 /**< maximal rank of a cut that could be scaled to integral coefficients (-1: unlimited) */ 80 #define DEFAULT_MAXWEIGHTRANGE 1e+03 /**< maximal valid range max(|weights|)/min(|weights|) of row weights */ 81 #define DEFAULT_STRENGTHEN TRUE /**< strengthen cut if integer variables are present */ 82 83 #define DEFAULT_MAXDEPTH -1 /**< node depth of separating cuts (-1: no limit) */ 84 #define DEFAULT_MAXROUNDS 25 /**< maximal number of separation rounds in a branching node (-1: no limit) */ 85 #define DEFAULT_MAXROUNDSROOT 100 /**< maximal number of separation rounds in the root node (-1: no limit) */ 86 #define DEFAULT_MAXINVCUTS 50 /**< maximal number of cuts investigated per iteration in a branching node */ 87 #define DEFAULT_MAXINVCUTSROOT 250 /**< maximal number of cuts investigated per iteration in the root node */ 88 #define DEFAULT_MAXCONFSDELAY 100000 /**< delay separation if number of conflict graph edges is larger than predefined value (-1: no limit) */ 89 90 #define MAKECONTINTEGRAL FALSE /**< convert continuous variable to integral variables in SCIPmakeRowIntegral() */ 91 92 93 /** separator data */ 94 struct SCIP_SepaData 95 { 96 SCIP_Bool strengthen; /**< strengthen cut if integer variables are present */ 97 SCIP_CONSHDLR* conshdlr; /**< SOS1 constraint handler */ 98 SCIP_Real maxweightrange; /**< maximal valid range max(|weights|)/min(|weights|) of row weights */ 99 int maxrank; /**< maximal rank of a cut that could not be scaled to integral coefficients (-1: unlimited) */ 100 int maxrankintegral; /**< maximal rank of a cut that could be scaled to integral coefficients (-1: unlimited) */ 101 int maxdepth; /**< node depth of separating cuts (-1: no limit) */ 102 int maxrounds; /**< maximal number of separation rounds in a branching node (-1: no limit) */ 103 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: no limit) */ 104 int maxinvcuts; /**< maximal number of cuts separated per iteration in a branching node */ 105 int maxinvcutsroot; /**< maximal number of cuts separated per iteration in the root node */ 106 int maxconfsdelay; /**< delay separation if number of conflict graph edges is larger than predefined value (-1: no limit) */ 107 int lastncutsfound; /**< total number of cuts found after last call of separator */ 108 }; 109 110 111 /** gets rank of variable corresponding to row of \f$B^{-1}\f$ */ 112 static 113 int getVarRank( 114 SCIP* scip, /**< SCIP pointer */ 115 SCIP_Real* binvrow, /**< row of \f$B^{-1}\f$ */ 116 SCIP_Real* rowsmaxval, /**< maximal row multiplicator from nonbasic matrix A_N */ 117 SCIP_Real maxweightrange, /**< maximal valid range max(|weights|)/min(|weights|) of row weights */ 118 SCIP_ROW** rows, /**< rows of LP relaxation of scip */ 119 int nrows /**< number of rows */ 120 ) 121 { 122 SCIP_Real maxweight = 0.0; 123 int maxrank = 0; 124 int r; 125 126 assert( scip != NULL ); 127 assert( binvrow != NULL || nrows == 0 ); 128 assert( rowsmaxval != NULL || nrows == 0 ); 129 assert( rows != NULL || nrows == 0 ); 130 131 /* compute maximum row weights resulting from multiplication */ 132 for (r = 0; r < nrows; ++r) 133 { 134 SCIP_Real val; 135 136 val = REALABS(binvrow[r] * rowsmaxval[r]);/*lint !e613*/ 137 if ( SCIPisGT(scip, val, maxweight) ) 138 maxweight = val; 139 } 140 141 /* compute rank */ 142 for (r = 0; r < nrows; ++r) 143 { 144 SCIP_Real val; 145 int rank; 146 147 val = REALABS(binvrow[r] * rowsmaxval[r]);/*lint !e613*/ 148 rank = SCIProwGetRank(rows[r]);/*lint !e613*/ 149 if ( rank > maxrank && SCIPisGT(scip, val * maxweightrange, maxweight) ) 150 maxrank = rank; 151 } 152 153 return maxrank; 154 } 155 156 157 /** gets the nonbasic coefficients of a simplex row */ 158 static 159 SCIP_RETCODE getSimplexCoefficients( 160 SCIP* scip, /**< SCIP pointer */ 161 SCIP_ROW** rows, /**< LP rows */ 162 int nrows, /**< number LP rows */ 163 SCIP_COL** cols, /**< LP columns */ 164 int ncols, /**< number of LP columns */ 165 SCIP_Real* coef, /**< row of \f$B^{-1} \cdot A\f$ */ 166 SCIP_Real* binvrow, /**< row of \f$B^{-1}\f$ */ 167 SCIP_Real* simplexcoefs, /**< pointer to store the nonbasic simplex-coefficients */ 168 int* nonbasicnumber /**< pointer to store the number of nonbasic simplex-coefficients */ 169 ) 170 { 171 int r; 172 int c; 173 174 assert( scip != NULL ); 175 assert( rows != NULL ); 176 assert( nonbasicnumber != NULL ); 177 assert( simplexcoefs != NULL ); 178 assert( cols != NULL ); 179 180 *nonbasicnumber = 0; 181 182 /* note: the non-slack variables have to be added first (see the function generateDisjCutSOS1()) */ 183 184 /* get simplex-coefficients of the non-basic non-slack variables */ 185 for (c = 0; c < ncols; ++c) 186 { 187 SCIP_COL* col; 188 189 col = cols[c]; 190 assert( col != NULL ); 191 if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER || SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ) 192 simplexcoefs[(*nonbasicnumber)++] = coef[c]; 193 } 194 195 /* get simplex-coefficients of the non-basic slack variables */ 196 for (r = 0; r < nrows; ++r) 197 { 198 SCIP_ROW* row; 199 row = rows[r]; 200 assert( row != NULL ); 201 202 if ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER || SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER ) 203 { 204 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, row) - SCIProwGetRhs(row)) || SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, row) - SCIProwGetLhs(row)) ); 205 206 simplexcoefs[(*nonbasicnumber)++] = binvrow[r]; 207 } 208 } 209 210 return SCIP_OKAY; 211 } 212 213 214 /** computes a disjunctive cut inequality based on two simplex tableau rows */ 215 static 216 SCIP_RETCODE generateDisjCutSOS1( 217 SCIP* scip, /**< SCIP pointer */ 218 SCIP_SEPA* sepa, /**< separator */ 219 int depth, /**< current depth */ 220 SCIP_ROW** rows, /**< LP rows */ 221 int nrows, /**< number of LP rows */ 222 SCIP_COL** cols, /**< LP columns */ 223 int ncols, /**< number of LP columns */ 224 int ndisjcuts, /**< number of disjunctive cuts found so far */ 225 SCIP_Bool scale, /**< should cut be scaled */ 226 SCIP_Bool strengthen, /**< should cut be strengthened if integer variables are present */ 227 SCIP_Real cutlhs1, /**< left hand side of the first simplex row */ 228 SCIP_Real cutlhs2, /**< left hand side of the second simplex row */ 229 SCIP_Real bound1, /**< bound of first simplex row */ 230 SCIP_Real bound2, /**< bound of second simplex row */ 231 SCIP_Real* simplexcoefs1, /**< simplex coefficients of first row */ 232 SCIP_Real* simplexcoefs2, /**< simplex coefficients of second row */ 233 SCIP_Real* cutcoefs, /**< pointer to store cut coefficients (length: nscipvars) */ 234 SCIP_ROW** row, /**< pointer to store disjunctive cut inequality */ 235 SCIP_Bool* madeintegral /**< pointer to store whether cut has been scaled to integral values */ 236 ) 237 { 238 char cutname[SCIP_MAXSTRLEN]; 239 SCIP_COL** rowcols; 240 SCIP_COL* col; 241 SCIP_Real* rowvals; 242 SCIP_Real lhsrow; 243 SCIP_Real rhsrow; 244 SCIP_Real cutlhs; 245 SCIP_Real sgn; 246 SCIP_Real lb; 247 SCIP_Real ub; 248 int nonbasicnumber = 0; 249 int rownnonz; 250 int ind; 251 int r; 252 int c; 253 254 assert( scip != NULL ); 255 assert( row != NULL ); 256 assert( rows != NULL ); 257 assert( cols != NULL ); 258 assert( simplexcoefs1 != NULL ); 259 assert( simplexcoefs2 != NULL ); 260 assert( cutcoefs != NULL ); 261 assert( sepa != NULL ); 262 assert( madeintegral != NULL ); 263 264 *madeintegral = FALSE; 265 266 /* check signs */ 267 if ( SCIPisFeasPositive(scip, cutlhs1) == SCIPisFeasPositive(scip, cutlhs2) ) 268 sgn = 1.0; 269 else 270 sgn = -1.0; 271 272 /* check bounds */ 273 if ( SCIPisInfinity(scip, REALABS(bound1)) || SCIPisInfinity(scip, REALABS(bound2)) ) 274 strengthen = FALSE; 275 276 /* compute left hand side of row (a later update is possible, see below) */ 277 cutlhs = sgn * cutlhs1 * cutlhs2; 278 279 /* add cut-coefficients of the non-basic non-slack variables */ 280 for (c = 0; c < ncols; ++c) 281 { 282 col = cols[c]; 283 assert( col != NULL ); 284 ind = SCIPcolGetLPPos(col); 285 assert( ind >= 0 ); 286 287 if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_LOWER ) 288 { 289 lb = SCIPcolGetLb(col); 290 291 /* for integer variables we may obtain stronger coefficients */ 292 if ( strengthen && SCIPcolIsIntegral(col) ) 293 { 294 SCIP_Real mval; 295 SCIP_Real mvalfloor; 296 SCIP_Real mvalceil; 297 298 mval = (cutlhs2 * simplexcoefs1[nonbasicnumber] - cutlhs1 * simplexcoefs2[nonbasicnumber]) / (cutlhs2 * bound1 + cutlhs1 * bound2); 299 mvalfloor = SCIPfloor(scip, mval); 300 mvalceil = SCIPceil(scip, mval); 301 302 cutcoefs[ind] = MIN(sgn * cutlhs2 * (simplexcoefs1[nonbasicnumber] - mvalfloor * bound1), sgn * cutlhs1 * (simplexcoefs2[nonbasicnumber] + mvalceil * bound2)); 303 assert( SCIPisFeasLE(scip, cutcoefs[ind], MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber])) ); 304 } 305 else 306 cutcoefs[ind] = MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]); 307 308 cutlhs += cutcoefs[ind] * lb; 309 ++nonbasicnumber; 310 } 311 else if ( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_UPPER ) 312 { 313 ub = SCIPcolGetUb(col); 314 315 /* for integer variables we may obtain stronger coefficients */ 316 if ( strengthen && SCIPcolIsIntegral(col) ) 317 { 318 SCIP_Real mval; 319 SCIP_Real mvalfloor; 320 SCIP_Real mvalceil; 321 322 mval = (cutlhs2 * simplexcoefs1[nonbasicnumber] - cutlhs1 * simplexcoefs2[nonbasicnumber]) / (cutlhs2 * bound1 + cutlhs1 * bound2); 323 mvalfloor = SCIPfloor(scip, -mval); 324 mvalceil = SCIPceil(scip, -mval); 325 326 cutcoefs[ind] = MAX(sgn * cutlhs2 * (simplexcoefs1[nonbasicnumber] + mvalfloor * bound1), sgn * cutlhs1 * (simplexcoefs2[nonbasicnumber] - mvalceil * bound2)); 327 assert( SCIPisFeasLE(scip, -cutcoefs[ind], -MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber])) ); 328 } 329 else 330 cutcoefs[ind] = MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]); 331 332 cutlhs += cutcoefs[ind] * ub; 333 ++nonbasicnumber; 334 } 335 else 336 { 337 assert( SCIPcolGetBasisStatus(col) != SCIP_BASESTAT_ZERO ); 338 cutcoefs[ind] = 0.0; 339 } 340 } 341 342 /* add cut-coefficients of the non-basic slack variables */ 343 for (r = 0; r < nrows; ++r) 344 { 345 rhsrow = SCIProwGetRhs(rows[r]) - SCIProwGetConstant(rows[r]); 346 lhsrow = SCIProwGetLhs(rows[r]) - SCIProwGetConstant(rows[r]); 347 348 assert( SCIProwGetBasisStatus(rows[r]) != SCIP_BASESTAT_ZERO ); 349 assert( SCIPisFeasZero(scip, lhsrow - rhsrow) || SCIPisNegative(scip, lhsrow - rhsrow) ); 350 assert( SCIProwIsInLP(rows[r]) ); 351 352 if ( SCIProwGetBasisStatus(rows[r]) != SCIP_BASESTAT_BASIC ) 353 { 354 SCIP_Real cutcoef; 355 356 if ( SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_UPPER ) 357 { 358 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, rows[r]) - SCIProwGetRhs(rows[r])) ); 359 360 cutcoef = MAX(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]); 361 cutlhs -= cutcoef * rhsrow; 362 ++nonbasicnumber; 363 } 364 else /* SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_LOWER */ 365 { 366 assert( SCIProwGetBasisStatus(rows[r]) == SCIP_BASESTAT_LOWER ); 367 assert( SCIPisFeasZero(scip, SCIPgetRowLPActivity(scip, rows[r]) - SCIProwGetLhs(rows[r])) ); 368 369 cutcoef = MIN(sgn * cutlhs2 * simplexcoefs1[nonbasicnumber], sgn * cutlhs1 * simplexcoefs2[nonbasicnumber]); 370 cutlhs -= cutcoef * lhsrow; 371 ++nonbasicnumber; 372 } 373 374 rownnonz = SCIProwGetNNonz(rows[r]); 375 rowvals = SCIProwGetVals(rows[r]); 376 rowcols = SCIProwGetCols(rows[r]); 377 378 for (c = 0; c < rownnonz; ++c) 379 { 380 ind = SCIPcolGetLPPos(rowcols[c]); 381 382 /* if column is not in LP, then return without generating cut */ 383 if ( ind < 0 ) 384 { 385 *row = NULL; 386 return SCIP_OKAY; 387 } 388 389 cutcoefs[ind] -= cutcoef * rowvals[c]; 390 } 391 } 392 } 393 394 /* create cut */ 395 (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "%s_%" SCIP_LONGINT_FORMAT "_%d", SCIPsepaGetName(sepa), SCIPgetNLPs(scip), ndisjcuts); 396 397 /* we create the cut as locally valid, SCIP will make it globally valid if we are at the root node */ 398 SCIP_CALL( SCIPcreateEmptyRowSepa(scip, row, sepa, cutname, cutlhs, SCIPinfinity(scip), TRUE, FALSE, TRUE) ); 399 400 SCIP_CALL( SCIPcacheRowExtensions(scip, *row) ); 401 for (c = 0; c < ncols; ++c) 402 { 403 ind = SCIPcolGetLPPos(cols[c]); 404 assert( ind >= 0 ); 405 if ( ! SCIPisFeasZero(scip, cutcoefs[ind]) ) 406 { 407 SCIP_CALL( SCIPaddVarToRow(scip, *row, SCIPcolGetVar(cols[c]), cutcoefs[ind] ) ); 408 } 409 } 410 SCIP_CALL( SCIPflushRowExtensions(scip, *row) ); 411 412 /* try to scale the cut to integral values 413 * @todo find better but still stable disjunctive cut settings 414 */ 415 if ( scale ) 416 { 417 SCIP_Longint maxdnom; 418 SCIP_Real maxscale; 419 420 assert( depth >= 0 ); 421 if( depth == 0 ) 422 { 423 maxdnom = 100; 424 maxscale = 100.0; 425 } 426 else 427 { 428 maxdnom = 10; 429 maxscale = 10.0; 430 } 431 432 SCIP_CALL( SCIPmakeRowIntegral(scip, *row, -SCIPepsilon(scip), SCIPsumepsilon(scip), maxdnom, maxscale, TRUE, madeintegral) ); 433 } 434 435 return SCIP_OKAY; 436 } 437 438 439 /* 440 * Callback methods 441 */ 442 443 /** copy method for separator plugins (called when SCIP copies plugins) */ 444 static 445 SCIP_DECL_SEPACOPY(sepaCopyDisjunctive) 446 { 447 assert( scip != NULL ); 448 assert( sepa != NULL ); 449 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 ); 450 451 /* call inclusion method of constraint handler */ 452 SCIP_CALL( SCIPincludeSepaDisjunctive(scip) ); 453 454 return SCIP_OKAY; 455 } 456 457 458 /** destructor of separator to free user data (called when SCIP is exiting) */ 459 static 460 SCIP_DECL_SEPAFREE(sepaFreeDisjunctive)/*lint --e{715}*/ 461 { 462 SCIP_SEPADATA* sepadata; 463 464 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 ); 465 466 /* free separator data */ 467 sepadata = SCIPsepaGetData(sepa); 468 assert( sepadata != NULL ); 469 470 SCIPfreeBlockMemory(scip, &sepadata); 471 472 SCIPsepaSetData(sepa, NULL); 473 474 return SCIP_OKAY; 475 } 476 477 478 /** solving process initialization method of separator */ 479 static 480 SCIP_DECL_SEPAEXITSOL(sepaInitsolDisjunctive) 481 { /*lint --e{715}*/ 482 SCIP_SEPADATA* sepadata; 483 484 sepadata = SCIPsepaGetData(sepa); 485 assert(sepadata != NULL); 486 487 sepadata->conshdlr = SCIPfindConshdlr(scip, "SOS1"); 488 489 return SCIP_OKAY; 490 } 491 492 493 /** LP solution separation method for disjunctive cuts */ 494 static 495 SCIP_DECL_SEPAEXECLP(sepaExeclpDisjunctive) 496 { 497 SCIP_SEPADATA* sepadata; 498 SCIP_CONSHDLR* conshdlr; 499 SCIP_DIGRAPH* conflictgraph; 500 SCIP_ROW** rows; 501 SCIP_COL** cols; 502 SCIP_Real* cutcoefs = NULL; 503 SCIP_Real* simplexcoefs1 = NULL; 504 SCIP_Real* simplexcoefs2 = NULL; 505 SCIP_Real* coef = NULL; 506 SCIP_Real* binvrow = NULL; 507 SCIP_Real* rowsmaxval = NULL; 508 SCIP_Real* violationarray = NULL; 509 int* fixings1 = NULL; 510 int* fixings2 = NULL; 511 int* basisind = NULL; 512 int* basisrow = NULL; 513 int* varrank = NULL; 514 int* edgearray = NULL; 515 int nedges; 516 int ndisjcuts; 517 int nrelevantedges; 518 int nsos1vars; 519 int nconss; 520 int maxcuts; 521 int ncalls; 522 int ncols; 523 int nrows; 524 int ind; 525 int j; 526 int i; 527 528 assert( sepa != NULL ); 529 assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 ); 530 assert( scip != NULL ); 531 assert( result != NULL ); 532 533 *result = SCIP_DIDNOTRUN; 534 535 if( !allowlocal ) 536 return SCIP_OKAY; 537 538 /* only generate disjunctive cuts if we are not close to terminating */ 539 if ( SCIPisStopped(scip) ) 540 return SCIP_OKAY; 541 542 /* only generate disjunctive cuts if an optimal LP solution is at hand */ 543 if ( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 544 return SCIP_OKAY; 545 546 /* only generate disjunctive cuts if the LP solution is basic */ 547 if ( ! SCIPisLPSolBasic(scip) ) 548 return SCIP_OKAY; 549 550 /* get LP data */ 551 SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) ); 552 SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) ); 553 554 /* return if LP has no columns or no rows */ 555 if ( ncols == 0 || nrows == 0 ) 556 return SCIP_OKAY; 557 558 assert( cols != NULL ); 559 assert( rows != NULL ); 560 561 /* get sepa data */ 562 sepadata = SCIPsepaGetData(sepa); 563 assert( sepadata != NULL ); 564 565 /* get constraint handler */ 566 conshdlr = sepadata->conshdlr; 567 if ( conshdlr == NULL ) 568 return SCIP_OKAY; 569 570 /* get number of constraints */ 571 nconss = SCIPconshdlrGetNConss(conshdlr); 572 if ( nconss == 0 ) 573 return SCIP_OKAY; 574 575 /* check for maxdepth < depth, maxinvcutsroot = 0 and maxinvcuts = 0 */ 576 if ( ( sepadata->maxdepth >= 0 && sepadata->maxdepth < depth ) 577 || ( depth == 0 && sepadata->maxinvcutsroot == 0 ) 578 || ( depth > 0 && sepadata->maxinvcuts == 0 ) ) 579 return SCIP_OKAY; 580 581 /* only call the cut separator a given number of times at each node */ 582 ncalls = SCIPsepaGetNCallsAtNode(sepa); 583 if ( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot) 584 || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) ) 585 return SCIP_OKAY; 586 587 /* get conflict graph and number of conflict graph edges (note that the digraph arcs were added in both directions) */ 588 conflictgraph = SCIPgetConflictgraphSOS1(conshdlr); 589 if( conflictgraph == NULL ) 590 return SCIP_OKAY; 591 592 nedges = (int)SCIPceil(scip, (SCIP_Real)SCIPdigraphGetNArcs(conflictgraph)/2); 593 594 /* if too many conflict graph edges, the separator can be slow: delay it until no other cuts have been found */ 595 if ( sepadata->maxconfsdelay >= 0 && nedges >= sepadata->maxconfsdelay ) 596 { 597 int ncutsfound; 598 599 ncutsfound = SCIPgetNCutsFound(scip); 600 if ( ncutsfound > sepadata->lastncutsfound || ! SCIPsepaWasLPDelayed(sepa) ) 601 { 602 sepadata->lastncutsfound = ncutsfound; 603 *result = SCIP_DELAYED; 604 return SCIP_OKAY; 605 } 606 } 607 608 /* check basis status */ 609 for (j = 0; j < ncols; ++j) 610 { 611 if ( SCIPcolGetBasisStatus(cols[j]) == SCIP_BASESTAT_ZERO ) 612 return SCIP_OKAY; 613 } 614 615 /* check accuracy of rows */ 616 for (j = 0; j < nrows; ++j) 617 { 618 SCIP_ROW* row; 619 620 row = rows[j]; 621 622 if ( ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_UPPER && ! SCIPisEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetRhs(row)) ) 623 || ( SCIProwGetBasisStatus(row) == SCIP_BASESTAT_LOWER && ! SCIPisEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row)) ) ) 624 return SCIP_OKAY; 625 } 626 627 /* get number of SOS1 variables */ 628 nsos1vars = SCIPgetNSOS1Vars(conshdlr); 629 630 /* allocate buffer arrays */ 631 SCIP_CALL( SCIPallocBufferArray(scip, &edgearray, nedges) ); 632 SCIP_CALL( SCIPallocBufferArray(scip, &fixings1, nedges) ); 633 SCIP_CALL( SCIPallocBufferArray(scip, &fixings2, nedges) ); 634 SCIP_CALL( SCIPallocBufferArray(scip, &violationarray, nedges) ); 635 636 /* get all violated conflicts {i, j} in the conflict graph and sort them based on the degree of a violation value */ 637 nrelevantedges = 0; 638 for (j = 0; j < nsos1vars; ++j) 639 { 640 SCIP_VAR* var; 641 642 var = SCIPnodeGetVarSOS1(conflictgraph, j); 643 644 if ( SCIPvarIsActive(var) && ! SCIPisFeasZero(scip, SCIPcolGetPrimsol(SCIPvarGetCol(var))) && SCIPcolGetBasisStatus(SCIPvarGetCol(var)) == SCIP_BASESTAT_BASIC ) 645 { 646 int* succ; 647 int nsucc; 648 649 /* get successors and number of successors */ 650 nsucc = SCIPdigraphGetNSuccessors(conflictgraph, j); 651 succ = SCIPdigraphGetSuccessors(conflictgraph, j); 652 653 for (i = 0; i < nsucc; ++i) 654 { 655 SCIP_VAR* varsucc; 656 int succind; 657 658 succind = succ[i]; 659 varsucc = SCIPnodeGetVarSOS1(conflictgraph, succind); 660 if ( SCIPvarIsActive(varsucc) && succind < j && ! SCIPisFeasZero(scip, SCIPgetSolVal(scip, NULL, varsucc) ) && 661 SCIPcolGetBasisStatus(SCIPvarGetCol(varsucc)) == SCIP_BASESTAT_BASIC ) 662 { 663 fixings1[nrelevantedges] = j; 664 fixings2[nrelevantedges] = succind; 665 edgearray[nrelevantedges] = nrelevantedges; 666 violationarray[nrelevantedges++] = SCIPgetSolVal(scip, NULL, var) * SCIPgetSolVal(scip, NULL, varsucc); 667 } 668 } 669 } 670 } 671 672 /* sort violation score values */ 673 if ( nrelevantedges > 0) 674 SCIPsortDownRealInt(violationarray, edgearray, nrelevantedges); 675 else 676 { 677 SCIPfreeBufferArrayNull(scip, &violationarray); 678 SCIPfreeBufferArrayNull(scip, &fixings2); 679 SCIPfreeBufferArrayNull(scip, &fixings1); 680 SCIPfreeBufferArrayNull(scip, &edgearray); 681 682 return SCIP_OKAY; 683 } 684 SCIPfreeBufferArrayNull(scip, &violationarray); 685 686 /* compute maximal number of cuts */ 687 if ( depth == 0 ) 688 maxcuts = MIN(sepadata->maxinvcutsroot, nrelevantedges); 689 else 690 maxcuts = MIN(sepadata->maxinvcuts, nrelevantedges); 691 assert( maxcuts > 0 ); 692 693 /* allocate buffer arrays */ 694 SCIP_CALL( SCIPallocBufferArray(scip, &varrank, ncols) ); 695 SCIP_CALL( SCIPallocBufferArray(scip, &rowsmaxval, nrows) ); 696 SCIP_CALL( SCIPallocBufferArray(scip, &basisrow, ncols) ); 697 SCIP_CALL( SCIPallocBufferArray(scip, &binvrow, nrows) ); 698 SCIP_CALL( SCIPallocBufferArray(scip, &coef, ncols) ); 699 SCIP_CALL( SCIPallocBufferArray(scip, &simplexcoefs1, ncols) ); 700 SCIP_CALL( SCIPallocBufferArray(scip, &simplexcoefs2, ncols) ); 701 SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, ncols) ); 702 SCIP_CALL( SCIPallocBufferArray(scip, &basisind, nrows) ); 703 704 /* get basis indices */ 705 SCIP_CALL( SCIPgetLPBasisInd(scip, basisind) ); 706 707 /* create vector "basisrow" with basisrow[column of non-slack basis variable] = corresponding row of B^-1; 708 * compute maximum absolute value of nonbasic row coefficients */ 709 for (j = 0; j < nrows; ++j) 710 { 711 SCIP_COL** rowcols; 712 SCIP_Real* rowvals; 713 SCIP_ROW* row; 714 SCIP_Real val; 715 SCIP_Real max = 0.0; 716 int nnonz; 717 718 /* fill basisrow vector */ 719 ind = basisind[j]; 720 if ( ind >= 0 ) 721 basisrow[ind] = j; 722 723 /* compute maximum absolute value of nonbasic row coefficients */ 724 row = rows[j]; 725 assert( row != NULL ); 726 rowvals = SCIProwGetVals(row); 727 nnonz = SCIProwGetNNonz(row); 728 rowcols = SCIProwGetCols(row); 729 730 for (i = 0; i < nnonz; ++i) 731 { 732 if ( SCIPcolGetBasisStatus(rowcols[i]) == SCIP_BASESTAT_LOWER || SCIPcolGetBasisStatus(rowcols[i]) == SCIP_BASESTAT_UPPER ) 733 { 734 val = REALABS(rowvals[i]); 735 if ( SCIPisFeasGT(scip, val, max) ) 736 max = REALABS(val); 737 } 738 } 739 740 /* handle slack variable coefficient and save maximum value */ 741 rowsmaxval[j] = MAX(max, 1.0); 742 } 743 744 /* initialize variable ranks with -1 */ 745 for (j = 0; j < ncols; ++j) 746 varrank[j] = -1; 747 748 /* free buffer array */ 749 SCIPfreeBufferArrayNull(scip, &basisind); 750 751 /* for the most promising disjunctions: try to generate disjunctive cuts */ 752 ndisjcuts = 0; 753 for (i = 0; i < maxcuts; ++i) 754 { 755 SCIP_Bool madeintegral; 756 SCIP_Real cutlhs1; 757 SCIP_Real cutlhs2; 758 SCIP_Real bound1; 759 SCIP_Real bound2; 760 SCIP_ROW* row = NULL; 761 SCIP_VAR* var; 762 SCIP_COL* col; 763 764 int nonbasicnumber; 765 int cutrank = 0; 766 int edgenumber; 767 int rownnonz; 768 769 edgenumber = edgearray[i]; 770 771 /* determine first simplex row */ 772 var = SCIPnodeGetVarSOS1(conflictgraph, fixings1[edgenumber]); 773 col = SCIPvarGetCol(var); 774 ind = SCIPcolGetLPPos(col); 775 assert( ind >= 0 ); 776 assert( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC ); 777 778 /* get the 'ind'th row of B^-1 and B^-1 \cdot A */ 779 SCIP_CALL( SCIPgetLPBInvRow(scip, basisrow[ind], binvrow, NULL, NULL) ); 780 SCIP_CALL( SCIPgetLPBInvARow(scip, basisrow[ind], binvrow, coef, NULL, NULL) ); 781 782 /* get the simplex-coefficients of the non-basic variables */ 783 SCIP_CALL( getSimplexCoefficients(scip, rows, nrows, cols, ncols, coef, binvrow, simplexcoefs1, &nonbasicnumber) ); 784 785 /* get rank of variable if not known already */ 786 if ( varrank[ind] < 0 ) 787 varrank[ind] = getVarRank(scip, binvrow, rowsmaxval, sepadata->maxweightrange, rows, nrows); 788 cutrank = MAX(cutrank, varrank[ind]); 789 790 /* get right hand side and bound of simplex talbeau row */ 791 cutlhs1 = SCIPcolGetPrimsol(col); 792 if ( SCIPisFeasPositive(scip, cutlhs1) ) 793 bound1 = SCIPcolGetUb(col); 794 else 795 bound1 = SCIPcolGetLb(col); 796 797 /* determine second simplex row */ 798 var = SCIPnodeGetVarSOS1(conflictgraph, fixings2[edgenumber]); 799 col = SCIPvarGetCol(var); 800 ind = SCIPcolGetLPPos(col); 801 assert( ind >= 0 ); 802 assert( SCIPcolGetBasisStatus(col) == SCIP_BASESTAT_BASIC ); 803 804 /* get the 'ind'th row of B^-1 and B^-1 \cdot A */ 805 SCIP_CALL( SCIPgetLPBInvRow(scip, basisrow[ind], binvrow, NULL, NULL) ); 806 SCIP_CALL( SCIPgetLPBInvARow(scip, basisrow[ind], binvrow, coef, NULL, NULL) ); 807 808 /* get the simplex-coefficients of the non-basic variables */ 809 SCIP_CALL( getSimplexCoefficients(scip, rows, nrows, cols, ncols, coef, binvrow, simplexcoefs2, &nonbasicnumber) ); 810 811 /* get rank of variable if not known already */ 812 if ( varrank[ind] < 0 ) 813 varrank[ind] = getVarRank(scip, binvrow, rowsmaxval, sepadata->maxweightrange, rows, nrows); 814 cutrank = MAX(cutrank, varrank[ind]); 815 816 /* get right hand side and bound of simplex talbeau row */ 817 cutlhs2 = SCIPcolGetPrimsol(col); 818 if ( SCIPisFeasPositive(scip, cutlhs2) ) 819 bound2 = SCIPcolGetUb(col); 820 else 821 bound2 = SCIPcolGetLb(col); 822 823 /* add coefficients to cut */ 824 SCIP_CALL( generateDisjCutSOS1(scip, sepa, depth, rows, nrows, cols, ncols, ndisjcuts, MAKECONTINTEGRAL, sepadata->strengthen, cutlhs1, cutlhs2, bound1, bound2, simplexcoefs1, simplexcoefs2, cutcoefs, &row, &madeintegral) ); 825 if ( row == NULL ) 826 continue; 827 828 /* raise cutrank for present cut */ 829 ++cutrank; 830 831 /* check if there are numerical evidences */ 832 if ( ( madeintegral && ( sepadata->maxrankintegral == -1 || cutrank <= sepadata->maxrankintegral ) ) 833 || ( ! madeintegral && ( sepadata->maxrank == -1 || cutrank <= sepadata->maxrank ) ) ) 834 { 835 /* possibly add cut to LP if it is useful; in case the lhs of the cut is minus infinity (due to scaling) the cut is useless */ 836 rownnonz = SCIProwGetNNonz(row); 837 if ( rownnonz > 0 && ! SCIPisInfinity(scip, -SCIProwGetLhs(row)) && ! SCIProwIsInLP(row) && SCIPisCutEfficacious(scip, NULL, row) ) 838 { 839 SCIP_Bool infeasible; 840 841 /* set cut rank */ 842 SCIProwChgRank(row, cutrank); 843 844 /* add cut */ 845 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) ); 846 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 847 if ( infeasible ) 848 { 849 *result = SCIP_CUTOFF; 850 break; 851 } 852 ++ndisjcuts; 853 } 854 } 855 856 /* release row */ 857 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 858 } 859 860 /* save total number of cuts found so far */ 861 sepadata->lastncutsfound = SCIPgetNCutsFound(scip); 862 863 /* evaluate the result of the separation */ 864 if ( *result != SCIP_CUTOFF ) 865 { 866 if ( ndisjcuts > 0 ) 867 *result = SCIP_SEPARATED; 868 else 869 *result = SCIP_DIDNOTFIND; 870 } 871 872 SCIPdebugMsg(scip, "Number of found disjunctive cuts: %d.\n", ndisjcuts); 873 874 /* free buffer arrays */ 875 SCIPfreeBufferArrayNull(scip, &cutcoefs); 876 SCIPfreeBufferArrayNull(scip, &simplexcoefs2); 877 SCIPfreeBufferArrayNull(scip, &simplexcoefs1); 878 SCIPfreeBufferArrayNull(scip, &coef); 879 SCIPfreeBufferArrayNull(scip, &binvrow); 880 SCIPfreeBufferArrayNull(scip, &basisrow); 881 SCIPfreeBufferArrayNull(scip, &rowsmaxval); 882 SCIPfreeBufferArrayNull(scip, &varrank); 883 SCIPfreeBufferArrayNull(scip, &fixings2); 884 SCIPfreeBufferArrayNull(scip, &fixings1); 885 SCIPfreeBufferArrayNull(scip, &edgearray); 886 887 return SCIP_OKAY; 888 } 889 890 891 /** creates the disjunctive cut separator and includes it in SCIP */ 892 SCIP_RETCODE SCIPincludeSepaDisjunctive( 893 SCIP* scip /**< SCIP data structure */ 894 ) 895 { 896 SCIP_SEPADATA* sepadata = NULL; 897 SCIP_SEPA* sepa = NULL; 898 899 /* create separator data */ 900 SCIP_CALL( SCIPallocBlockMemory(scip, &sepadata) ); 901 sepadata->conshdlr = NULL; 902 sepadata->lastncutsfound = 0; 903 904 /* include separator */ 905 SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST, 906 SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpDisjunctive, NULL, sepadata) ); 907 908 assert( sepa != NULL ); 909 910 /* set non fundamental callbacks via setter functions */ 911 SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyDisjunctive) ); 912 SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeDisjunctive) ); 913 SCIP_CALL( SCIPsetSepaInitsol(scip, sepa, sepaInitsolDisjunctive) ); 914 915 /* add separator parameters */ 916 SCIP_CALL( SCIPaddBoolParam(scip, "separating/" SEPA_NAME "/strengthen", 917 "strengthen cut if integer variables are present.", 918 &sepadata->strengthen, TRUE, DEFAULT_STRENGTHEN, NULL, NULL) ); 919 920 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxdepth", 921 "node depth of separating bipartite disjunctive cuts (-1: no limit)", 922 &sepadata->maxdepth, TRUE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) ); 923 924 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrounds", 925 "maximal number of separation rounds per iteration in a branching node (-1: no limit)", 926 &sepadata->maxrounds, TRUE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) ); 927 928 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxroundsroot", 929 "maximal number of separation rounds in the root node (-1: no limit)", 930 &sepadata->maxroundsroot, TRUE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) ); 931 932 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxinvcuts", 933 "maximal number of cuts investigated per iteration in a branching node", 934 &sepadata->maxinvcuts, TRUE, DEFAULT_MAXINVCUTS, 0, INT_MAX, NULL, NULL) ); 935 936 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxinvcutsroot", 937 "maximal number of cuts investigated per iteration in the root node", 938 &sepadata->maxinvcutsroot, TRUE, DEFAULT_MAXINVCUTSROOT, 0, INT_MAX, NULL, NULL) ); 939 940 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxconfsdelay", 941 "delay separation if number of conflict graph edges is larger than predefined value (-1: no limit)", 942 &sepadata->maxconfsdelay, TRUE, DEFAULT_MAXCONFSDELAY, -1, INT_MAX, NULL, NULL) ); 943 944 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrank", 945 "maximal rank of a disj. cut that could not be scaled to integral coefficients (-1: unlimited)", 946 &sepadata->maxrank, FALSE, DEFAULT_MAXRANK, -1, INT_MAX, NULL, NULL) ); 947 948 SCIP_CALL( SCIPaddIntParam(scip, "separating/" SEPA_NAME "/maxrankintegral", 949 "maximal rank of a disj. cut that could be scaled to integral coefficients (-1: unlimited)", 950 &sepadata->maxrankintegral, FALSE, DEFAULT_MAXRANKINTEGRAL, -1, INT_MAX, NULL, NULL) ); 951 952 SCIP_CALL( SCIPaddRealParam(scip, "separating/" SEPA_NAME "/maxweightrange", 953 "maximal valid range max(|weights|)/min(|weights|) of row weights", 954 &sepadata->maxweightrange, TRUE, DEFAULT_MAXWEIGHTRANGE, 1.0, SCIP_REAL_MAX, NULL, NULL) ); 955 956 return SCIP_OKAY; 957 } 958