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 heur_shifting.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief LP rounding heuristic that tries to recover from intermediate infeasibilities and shifts continuous variables 28 * @author Tobias Achterberg 29 */ 30 31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 32 33 #include "blockmemshell/memory.h" 34 #include "scip/heur_shifting.h" 35 #include "scip/pub_heur.h" 36 #include "scip/pub_lp.h" 37 #include "scip/pub_message.h" 38 #include "scip/pub_misc.h" 39 #include "scip/pub_var.h" 40 #include "scip/scip_branch.h" 41 #include "scip/scip_heur.h" 42 #include "scip/scip_lp.h" 43 #include "scip/scip_mem.h" 44 #include "scip/scip_message.h" 45 #include "scip/scip_numerics.h" 46 #include "scip/scip_prob.h" 47 #include "scip/scip_randnumgen.h" 48 #include "scip/scip_sol.h" 49 #include "scip/scip_solvingstats.h" 50 #include <string.h> 51 52 #define HEUR_NAME "shifting" 53 #define HEUR_DESC "LP rounding heuristic with infeasibility recovering also using continuous variables" 54 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_ROUNDING 55 #define HEUR_PRIORITY -5000 56 #define HEUR_FREQ 10 57 #define HEUR_FREQOFS 0 58 #define HEUR_MAXDEPTH -1 59 #define HEUR_TIMING SCIP_HEURTIMING_DURINGLPLOOP 60 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */ 61 62 #define MAXSHIFTINGS 50 /**< maximal number of non improving shiftings */ 63 #define WEIGHTFACTOR 1.1 64 #define DEFAULT_RANDSEED 31 /**< initial random seed */ 65 66 67 /* locally defined heuristic data */ 68 struct SCIP_HeurData 69 { 70 SCIP_SOL* sol; /**< working solution */ 71 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */ 72 SCIP_Longint lastlp; /**< last LP number where the heuristic was applied */ 73 }; 74 75 76 /* 77 * local methods 78 */ 79 80 /** update row violation arrays after a row's activity value changed */ 81 static 82 void updateViolations( 83 SCIP* scip, /**< SCIP data structure */ 84 SCIP_ROW* row, /**< LP row */ 85 SCIP_ROW** violrows, /**< array with currently violated rows */ 86 int* violrowpos, /**< position of LP rows in violrows array */ 87 int* nviolrows, /**< pointer to the number of currently violated rows */ 88 int* nviolfracrows, /**< pointer to the number of violated rows with fractional candidates */ 89 int* nfracsinrow, /**< array with number of fractional variables for every row */ 90 SCIP_Real oldactivity, /**< old activity value of LP row */ 91 SCIP_Real newactivity /**< new activity value of LP row */ 92 ) 93 { 94 SCIP_Real lhs; 95 SCIP_Real rhs; 96 SCIP_Bool oldviol; 97 SCIP_Bool newviol; 98 99 assert(violrows != NULL); 100 assert(violrowpos != NULL); 101 assert(nviolrows != NULL); 102 103 lhs = SCIProwGetLhs(row); 104 rhs = SCIProwGetRhs(row); 105 106 /* SCIPisFeasLT cannot handle comparing different infinities. To prevent this, we make a case distinction. */ 107 if( !(SCIPisInfinity(scip, oldactivity) || SCIPisInfinity(scip, -oldactivity)) ) 108 { 109 oldviol = (SCIPisFeasLT(scip, oldactivity, lhs) || SCIPisFeasGT(scip, oldactivity, rhs)); 110 } 111 else 112 { 113 oldviol = (SCIPisInfinity(scip, -oldactivity) && !SCIPisInfinity(scip, -lhs)) || 114 (SCIPisInfinity(scip, oldactivity) && !SCIPisInfinity(scip, rhs)); 115 } 116 117 /* SCIPisFeasLT cannot handle comparing different infinities. To prevent this, we make a case distinction. */ 118 if( !(SCIPisInfinity(scip, newactivity) || SCIPisInfinity(scip, -newactivity)) ) 119 { 120 newviol = (SCIPisFeasLT(scip, newactivity, lhs) || SCIPisFeasGT(scip, newactivity, rhs)); 121 } 122 else 123 { 124 newviol = (SCIPisInfinity(scip, -newactivity) && !SCIPisInfinity(scip, -lhs)) || 125 (SCIPisInfinity(scip, newactivity) && !SCIPisInfinity(scip, rhs)); 126 } 127 128 if( oldviol != newviol ) 129 { 130 int rowpos; 131 int violpos; 132 133 rowpos = SCIProwGetLPPos(row); 134 assert(rowpos >= 0); 135 136 if( oldviol ) 137 { 138 /* the row violation was repaired: remove row from violrows array, decrease violation count */ 139 violpos = violrowpos[rowpos]; 140 assert(0 <= violpos && violpos < *nviolrows); 141 assert(violrows[violpos] == row); 142 violrowpos[rowpos] = -1; 143 144 /* first, move the row to the end of the subset of violated rows with fractional variables */ 145 if( nfracsinrow[rowpos] > 0 ) 146 { 147 assert(violpos < *nviolfracrows); 148 149 /* replace with last violated row containing fractional variables */ 150 if( violpos != *nviolfracrows - 1 ) 151 { 152 violrows[violpos] = violrows[*nviolfracrows - 1]; 153 violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos; 154 violpos = *nviolfracrows - 1; 155 } 156 (*nviolfracrows)--; 157 } 158 159 assert(violpos >= *nviolfracrows); 160 161 /* swap row at the end of the violated array to the position of this row and decrease the counter */ 162 if( violpos != *nviolrows - 1 ) 163 { 164 violrows[violpos] = violrows[*nviolrows - 1]; 165 violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos; 166 } 167 (*nviolrows)--; 168 } 169 else 170 { 171 /* the row is now violated: add row to violrows array, increase violation count */ 172 assert(violrowpos[rowpos] == -1); 173 violrows[*nviolrows] = row; 174 violrowpos[rowpos] = *nviolrows; 175 (*nviolrows)++; 176 177 /* if the row contains fractional variables, swap with the last violated row that has no fractional variables 178 * at position *nviolfracrows 179 */ 180 if( nfracsinrow[rowpos] > 0 ) 181 { 182 if( *nviolfracrows < *nviolrows - 1 ) 183 { 184 assert(nfracsinrow[SCIProwGetLPPos(violrows[*nviolfracrows])] == 0); 185 186 violrows[*nviolrows - 1] = violrows[*nviolfracrows]; 187 violrowpos[SCIProwGetLPPos(violrows[*nviolrows - 1])] = *nviolrows - 1; 188 189 violrows[*nviolfracrows] = row; 190 violrowpos[rowpos] = *nviolfracrows; 191 } 192 (*nviolfracrows)++; 193 } 194 } 195 } 196 } 197 198 /** update row activities after a variable's solution value changed */ 199 static 200 SCIP_RETCODE updateActivities( 201 SCIP* scip, /**< SCIP data structure */ 202 SCIP_Real* activities, /**< LP row activities */ 203 SCIP_ROW** violrows, /**< array with currently violated rows */ 204 int* violrowpos, /**< position of LP rows in violrows array */ 205 int* nviolrows, /**< pointer to the number of currently violated rows */ 206 int* nviolfracrows, /**< pointer to the number of violated rows with fractional candidates */ 207 int* nfracsinrow, /**< array with number of fractional variables for every row */ 208 int nlprows, /**< number of rows in current LP */ 209 SCIP_VAR* var, /**< variable that has been changed */ 210 SCIP_Real oldsolval, /**< old solution value of variable */ 211 SCIP_Real newsolval /**< new solution value of variable */ 212 ) 213 { 214 SCIP_COL* col; 215 SCIP_ROW** colrows; 216 SCIP_Real* colvals; 217 SCIP_Real delta; 218 int ncolrows; 219 int r; 220 221 assert(activities != NULL); 222 assert(nviolrows != NULL); 223 assert(0 <= *nviolrows && *nviolrows <= nlprows); 224 225 delta = newsolval - oldsolval; 226 col = SCIPvarGetCol(var); 227 colrows = SCIPcolGetRows(col); 228 colvals = SCIPcolGetVals(col); 229 ncolrows = SCIPcolGetNLPNonz(col); 230 assert(ncolrows == 0 || (colrows != NULL && colvals != NULL)); 231 232 for( r = 0; r < ncolrows; ++r ) 233 { 234 SCIP_ROW* row; 235 int rowpos; 236 237 row = colrows[r]; 238 rowpos = SCIProwGetLPPos(row); 239 assert(-1 <= rowpos && rowpos < nlprows); 240 241 if( rowpos >= 0 && !SCIProwIsLocal(row) ) 242 { 243 SCIP_Real oldactivity; 244 SCIP_Real newactivity; 245 246 assert(SCIProwIsInLP(row)); 247 248 /* update row activity */ 249 oldactivity = activities[rowpos]; 250 if( !SCIPisInfinity(scip, -oldactivity) && !SCIPisInfinity(scip, oldactivity) ) 251 { 252 newactivity = oldactivity + delta * colvals[r]; 253 if( SCIPisInfinity(scip, newactivity) ) 254 newactivity = SCIPinfinity(scip); 255 else if( SCIPisInfinity(scip, -newactivity) ) 256 newactivity = -SCIPinfinity(scip); 257 activities[rowpos] = newactivity; 258 259 /* update row violation arrays */ 260 updateViolations(scip, row, violrows, violrowpos, nviolrows, nviolfracrows, nfracsinrow, oldactivity, newactivity); 261 } 262 } 263 } 264 265 return SCIP_OKAY; 266 } 267 268 /** returns a variable, that pushes activity of the row in the given direction with minimal negative impact on other rows; 269 * if variables have equal impact, chooses the one with best objective value improvement in corresponding direction; 270 * prefer fractional integers over other variables in order to become integral during the process; 271 * shifting in a direction is forbidden, if this forces the objective value over the upper bound, or if the variable 272 * was already shifted in the opposite direction 273 */ 274 static 275 SCIP_RETCODE selectShifting( 276 SCIP* scip, /**< SCIP data structure */ 277 SCIP_SOL* sol, /**< primal solution */ 278 SCIP_ROW* row, /**< LP row */ 279 SCIP_Real rowactivity, /**< activity of LP row */ 280 int direction, /**< should the activity be increased (+1) or decreased (-1)? */ 281 SCIP_Real* nincreases, /**< array with weighted number of increasings per variables */ 282 SCIP_Real* ndecreases, /**< array with weighted number of decreasings per variables */ 283 SCIP_Real increaseweight, /**< current weight of increase/decrease updates */ 284 SCIP_VAR** shiftvar, /**< pointer to store the shifting variable, returns NULL if impossible */ 285 SCIP_Real* oldsolval, /**< pointer to store old solution value of shifting variable */ 286 SCIP_Real* newsolval /**< pointer to store new (shifted) solution value of shifting variable */ 287 ) 288 { 289 SCIP_COL** rowcols; 290 SCIP_Real* rowvals; 291 int nrowcols; 292 SCIP_Real activitydelta; 293 SCIP_Real bestshiftscore; 294 SCIP_Real bestdeltaobj; 295 int c; 296 297 assert(direction == +1 || direction == -1); 298 assert(nincreases != NULL); 299 assert(ndecreases != NULL); 300 assert(shiftvar != NULL); 301 assert(oldsolval != NULL); 302 assert(newsolval != NULL); 303 304 /* get row entries */ 305 rowcols = SCIProwGetCols(row); 306 rowvals = SCIProwGetVals(row); 307 nrowcols = SCIProwGetNLPNonz(row); 308 309 /* calculate how much the activity must be shifted in order to become feasible */ 310 activitydelta = (direction == +1 ? SCIProwGetLhs(row) - rowactivity : SCIProwGetRhs(row) - rowactivity); 311 assert((direction == +1 && SCIPisPositive(scip, activitydelta)) 312 || (direction == -1 && SCIPisNegative(scip, activitydelta))); 313 314 /* select shifting variable */ 315 bestshiftscore = SCIP_REAL_MAX; 316 bestdeltaobj = SCIPinfinity(scip); 317 *shiftvar = NULL; 318 *newsolval = 0.0; 319 *oldsolval = 0.0; 320 for( c = 0; c < nrowcols; ++c ) 321 { 322 SCIP_COL* col; 323 SCIP_VAR* var; 324 SCIP_Real val; 325 SCIP_Real solval; 326 SCIP_Real shiftval; 327 SCIP_Real shiftscore; 328 SCIP_Bool isinteger; 329 SCIP_Bool isfrac; 330 SCIP_Bool increase; 331 332 col = rowcols[c]; 333 var = SCIPcolGetVar(col); 334 val = rowvals[c]; 335 assert(!SCIPisZero(scip, val)); 336 solval = SCIPgetSolVal(scip, sol, var); 337 338 isinteger = (SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER); 339 isfrac = (isinteger && !SCIPisFeasIntegral(scip, solval)); 340 increase = (direction * val > 0.0); 341 342 /* calculate the score of the shifting (prefer smaller values) */ 343 if( isfrac ) 344 shiftscore = increase ? -1.0 / (SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) + 1.0) : 345 -1.0 / (SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) + 1.0); 346 else 347 { 348 int probindex; 349 probindex = SCIPvarGetProbindex(var); 350 351 if( increase ) 352 shiftscore = ndecreases[probindex]/increaseweight; 353 else 354 shiftscore = nincreases[probindex]/increaseweight; 355 if( isinteger ) 356 shiftscore += 1.0; 357 } 358 359 if( shiftscore <= bestshiftscore ) 360 { 361 SCIP_Real deltaobj; 362 363 if( !increase ) 364 { 365 /* shifting down */ 366 assert(direction * val < 0.0); 367 if( isfrac ) 368 shiftval = SCIPfeasFloor(scip, solval); 369 else 370 { 371 SCIP_Real lb; 372 373 assert(activitydelta/val < 0.0); 374 shiftval = solval + activitydelta/val; 375 assert(shiftval <= solval); /* may be equal due to numerical digit erasement in the subtraction */ 376 if( SCIPvarIsIntegral(var) ) 377 shiftval = SCIPfeasFloor(scip, shiftval); 378 lb = SCIPvarGetLbGlobal(var); 379 shiftval = MAX(shiftval, lb); 380 } 381 } 382 else 383 { 384 /* shifting up */ 385 assert(direction * val > 0.0); 386 if( isfrac ) 387 shiftval = SCIPfeasCeil(scip, solval); 388 else 389 { 390 SCIP_Real ub; 391 392 assert(activitydelta/val > 0.0); 393 shiftval = solval + activitydelta/val; 394 assert(shiftval >= solval); /* may be equal due to numerical digit erasement in the subtraction */ 395 if( SCIPvarIsIntegral(var) ) 396 shiftval = SCIPfeasCeil(scip, shiftval); 397 ub = SCIPvarGetUbGlobal(var); 398 shiftval = MIN(shiftval, ub); 399 } 400 } 401 402 if( (SCIPisInfinity(scip, shiftval) && SCIPisInfinity(scip, SCIPvarGetUbGlobal(var))) 403 || (SCIPisInfinity(scip, -shiftval) && SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var))) 404 || SCIPisEQ(scip, shiftval, solval) ) 405 continue; 406 407 deltaobj = SCIPvarGetObj(var) * (shiftval - solval); 408 if( shiftscore < bestshiftscore || deltaobj < bestdeltaobj ) 409 { 410 bestshiftscore = shiftscore; 411 bestdeltaobj = deltaobj; 412 *shiftvar = var; 413 *oldsolval = solval; 414 *newsolval = shiftval; 415 } 416 } 417 } 418 419 return SCIP_OKAY; 420 } 421 422 /** returns a fractional variable, that has most impact on rows in opposite direction, i.e. that is most crucial to 423 * fix in the other direction; 424 * if variables have equal impact, chooses the one with best objective value improvement in corresponding direction; 425 * shifting in a direction is forbidden, if this forces the objective value over the upper bound 426 */ 427 static 428 SCIP_RETCODE selectEssentialRounding( 429 SCIP* scip, /**< SCIP data structure */ 430 SCIP_SOL* sol, /**< primal solution */ 431 SCIP_Real minobj, /**< minimal objective value possible after shifting remaining fractional vars */ 432 SCIP_VAR** lpcands, /**< fractional variables in LP */ 433 int nlpcands, /**< number of fractional variables in LP */ 434 SCIP_VAR** shiftvar, /**< pointer to store the shifting variable, returns NULL if impossible */ 435 SCIP_Real* oldsolval, /**< old (fractional) solution value of shifting variable */ 436 SCIP_Real* newsolval /**< new (shifted) solution value of shifting variable */ 437 ) 438 { 439 SCIP_VAR* var; 440 SCIP_Real solval; 441 SCIP_Real shiftval; 442 SCIP_Real obj; 443 SCIP_Real deltaobj; 444 SCIP_Real bestdeltaobj; 445 int maxnlocks; 446 int nlocks; 447 int v; 448 449 assert(shiftvar != NULL); 450 assert(oldsolval != NULL); 451 assert(newsolval != NULL); 452 453 /* select shifting variable */ 454 maxnlocks = -1; 455 bestdeltaobj = SCIPinfinity(scip); 456 *shiftvar = NULL; 457 for( v = 0; v < nlpcands; ++v ) 458 { 459 var = lpcands[v]; 460 assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER); 461 462 solval = SCIPgetSolVal(scip, sol, var); 463 if( !SCIPisFeasIntegral(scip, solval) ) 464 { 465 obj = SCIPvarGetObj(var); 466 467 /* shifting down */ 468 nlocks = SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL); 469 if( nlocks >= maxnlocks ) 470 { 471 shiftval = SCIPfeasFloor(scip, solval); 472 deltaobj = obj * (shiftval - solval); 473 if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj - obj < SCIPgetCutoffbound(scip) ) 474 { 475 maxnlocks = nlocks; 476 bestdeltaobj = deltaobj; 477 *shiftvar = var; 478 *oldsolval = solval; 479 *newsolval = shiftval; 480 } 481 } 482 483 /* shifting up */ 484 nlocks = SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL); 485 if( nlocks >= maxnlocks ) 486 { 487 shiftval = SCIPfeasCeil(scip, solval); 488 deltaobj = obj * (shiftval - solval); 489 if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj + obj < SCIPgetCutoffbound(scip) ) 490 { 491 maxnlocks = nlocks; 492 bestdeltaobj = deltaobj; 493 *shiftvar = var; 494 *oldsolval = solval; 495 *newsolval = shiftval; 496 } 497 } 498 } 499 } 500 501 return SCIP_OKAY; 502 } 503 504 /** adds a given value to the fractionality counters of the rows in which the given variable appears */ 505 static 506 void addFracCounter( 507 int* nfracsinrow, /**< array to store number of fractional variables per row */ 508 SCIP_ROW** violrows, /**< array with currently violated rows */ 509 int* violrowpos, /**< position of LP rows in violrows array */ 510 int* nviolfracrows, /**< pointer to store the number of violated rows with fractional variables */ 511 int nviolrows, /**< the number of currently violated rows (stays unchanged in this method) */ 512 int nlprows, /**< number of rows in LP */ 513 SCIP_VAR* var, /**< variable for which the counting should be updated */ 514 int incval /**< value that should be added to the corresponding array entries */ 515 ) 516 { 517 SCIP_COL* col; 518 SCIP_ROW** rows; 519 int nrows; 520 int r; 521 522 assert(incval != 0); 523 assert(nviolrows >= *nviolfracrows); 524 col = SCIPvarGetCol(var); 525 rows = SCIPcolGetRows(col); 526 nrows = SCIPcolGetNLPNonz(col); 527 for( r = 0; r < nrows; ++r ) 528 { 529 int rowlppos; 530 int theviolrowpos; 531 532 rowlppos = SCIProwGetLPPos(rows[r]); 533 assert(0 <= rowlppos && rowlppos < nlprows); 534 assert(!SCIProwIsLocal(rows[r]) || violrowpos[rowlppos] == -1); 535 536 /* skip local rows */ 537 if( SCIProwIsLocal(rows[r]) ) 538 continue; 539 540 nfracsinrow[rowlppos] += incval; 541 assert(nfracsinrow[rowlppos] >= 0); 542 theviolrowpos = violrowpos[rowlppos]; 543 544 /* swap positions in violrows array if fractionality has changed to 0 */ 545 if( theviolrowpos >= 0 ) 546 { 547 /* first case: the number of fractional variables has become zero: swap row in violrows array to the 548 * second part, containing only violated rows without fractional variables 549 */ 550 if( nfracsinrow[rowlppos] == 0 ) 551 { 552 assert(theviolrowpos <= *nviolfracrows - 1); 553 554 /* nothing to do if row is already at the end of the first part, otherwise, swap it to the last position 555 * and decrease the counter */ 556 if( theviolrowpos < *nviolfracrows - 1 ) 557 { 558 violrows[theviolrowpos] = violrows[*nviolfracrows - 1]; 559 violrows[*nviolfracrows - 1] = rows[r]; 560 561 violrowpos[SCIProwGetLPPos(violrows[theviolrowpos])] = theviolrowpos; 562 violrowpos[rowlppos] = *nviolfracrows - 1; 563 } 564 (*nviolfracrows)--; 565 } 566 /* second interesting case: the number of fractional variables was zero before, swap it with the first 567 * violated row without fractional variables 568 */ 569 else if( nfracsinrow[rowlppos] == incval ) 570 { 571 assert(theviolrowpos >= *nviolfracrows); 572 573 /* nothing to do if the row is exactly located at index *nviolfracrows */ 574 if( theviolrowpos > *nviolfracrows ) 575 { 576 violrows[theviolrowpos] = violrows[*nviolfracrows]; 577 violrows[*nviolfracrows] = rows[r]; 578 579 violrowpos[SCIProwGetLPPos(violrows[theviolrowpos])] = theviolrowpos; 580 violrowpos[rowlppos] = *nviolfracrows; 581 } 582 (*nviolfracrows)++; 583 } 584 } 585 } 586 } 587 588 589 /* 590 * Callback methods 591 */ 592 593 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 594 static 595 SCIP_DECL_HEURCOPY(heurCopyShifting) 596 { /*lint --e{715}*/ 597 assert(scip != NULL); 598 assert(heur != NULL); 599 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 600 601 /* call inclusion method of primal heuristic */ 602 SCIP_CALL( SCIPincludeHeurShifting(scip) ); 603 604 return SCIP_OKAY; 605 } 606 607 608 /** initialization method of primal heuristic (called after problem was transformed) */ 609 static 610 SCIP_DECL_HEURINIT(heurInitShifting) /*lint --e{715}*/ 611 { /*lint --e{715}*/ 612 SCIP_HEURDATA* heurdata; 613 614 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 615 assert(SCIPheurGetData(heur) == NULL); 616 617 /* create heuristic data */ 618 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 619 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) ); 620 heurdata->lastlp = -1; 621 622 /* create random number generator */ 623 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, 624 DEFAULT_RANDSEED, TRUE) ); 625 626 SCIPheurSetData(heur, heurdata); 627 628 return SCIP_OKAY; 629 } 630 631 /** deinitialization method of primal heuristic (called before transformed problem is freed) */ 632 static 633 SCIP_DECL_HEUREXIT(heurExitShifting) /*lint --e{715}*/ 634 { /*lint --e{715}*/ 635 SCIP_HEURDATA* heurdata; 636 637 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 638 639 /* free heuristic data */ 640 heurdata = SCIPheurGetData(heur); 641 assert(heurdata != NULL); 642 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) ); 643 644 /* free random number generator */ 645 SCIPfreeRandom(scip, &heurdata->randnumgen); 646 647 SCIPfreeBlockMemory(scip, &heurdata); 648 SCIPheurSetData(heur, NULL); 649 650 return SCIP_OKAY; 651 } 652 653 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */ 654 static 655 SCIP_DECL_HEURINITSOL(heurInitsolShifting) 656 { 657 SCIP_HEURDATA* heurdata; 658 659 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 660 661 heurdata = SCIPheurGetData(heur); 662 assert(heurdata != NULL); 663 heurdata->lastlp = -1; 664 665 return SCIP_OKAY; 666 } 667 668 669 /** execution method of primal heuristic */ 670 static 671 SCIP_DECL_HEUREXEC(heurExecShifting) /*lint --e{715}*/ 672 { /*lint --e{715}*/ 673 SCIP_HEURDATA* heurdata; 674 SCIP_SOL* sol; 675 SCIP_VAR** lpcands; 676 SCIP_Real* lpcandssol; 677 SCIP_ROW** lprows; 678 SCIP_Real* activities; 679 SCIP_ROW** violrows; 680 SCIP_Real* nincreases; 681 SCIP_Real* ndecreases; 682 int* violrowpos; 683 int* nfracsinrow; 684 SCIP_Real increaseweight; 685 SCIP_Real obj; 686 SCIP_Real bestshiftval; 687 SCIP_Real minobj; 688 int nlpcands; 689 int nlprows; 690 int nvars; 691 int nfrac; 692 int nviolrows; 693 int nviolfracrows; 694 int nprevviolrows; 695 int minnviolrows; 696 int nnonimprovingshifts; 697 int c; 698 int r; 699 SCIP_Longint nlps; 700 SCIP_Longint ncalls; 701 SCIP_Longint nsolsfound; 702 SCIP_Longint nnodes; 703 704 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 705 assert(scip != NULL); 706 assert(result != NULL); 707 assert(SCIPhasCurrentNodeLP(scip)); 708 709 *result = SCIP_DIDNOTRUN; 710 711 /* only call heuristic, if an optimal LP solution is at hand */ 712 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 713 return SCIP_OKAY; 714 715 /* only call heuristic, if the LP objective value is smaller than the cutoff bound */ 716 if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) ) 717 return SCIP_OKAY; 718 719 /* get heuristic data */ 720 heurdata = SCIPheurGetData(heur); 721 assert(heurdata != NULL); 722 723 /* don't call heuristic, if we have already processed the current LP solution */ 724 nlps = SCIPgetNLPs(scip); 725 if( nlps == heurdata->lastlp ) 726 return SCIP_OKAY; 727 heurdata->lastlp = nlps; 728 729 /* don't call heuristic, if it was not successful enough in the past */ 730 ncalls = SCIPheurGetNCalls(heur); 731 nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + SCIPheurGetNSolsFound(heur); 732 nnodes = SCIPgetNNodes(scip); 733 if( nnodes % ((ncalls/100)/(nsolsfound+1)+1) != 0 ) 734 return SCIP_OKAY; 735 736 /* get fractional variables, that should be integral */ 737 /* todo check if heuristic should include implicit integer variables for its calculations */ 738 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL) ); 739 nfrac = nlpcands; 740 741 /* only call heuristic, if LP solution is fractional */ 742 if( nfrac == 0 ) 743 return SCIP_OKAY; 744 745 *result = SCIP_DIDNOTFIND; 746 747 /* get LP rows */ 748 SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nlprows) ); 749 750 SCIPdebugMsg(scip, "executing shifting heuristic: %d LP rows, %d fractionals\n", nlprows, nfrac); 751 752 /* get memory for activities, violated rows, and row violation positions */ 753 nvars = SCIPgetNVars(scip); 754 SCIP_CALL( SCIPallocBufferArray(scip, &activities, nlprows) ); 755 SCIP_CALL( SCIPallocBufferArray(scip, &violrows, nlprows) ); 756 SCIP_CALL( SCIPallocBufferArray(scip, &violrowpos, nlprows) ); 757 SCIP_CALL( SCIPallocBufferArray(scip, &nfracsinrow, nlprows) ); 758 SCIP_CALL( SCIPallocBufferArray(scip, &nincreases, nvars) ); 759 SCIP_CALL( SCIPallocBufferArray(scip, &ndecreases, nvars) ); 760 BMSclearMemoryArray(nfracsinrow, nlprows); 761 BMSclearMemoryArray(nincreases, nvars); 762 BMSclearMemoryArray(ndecreases, nvars); 763 764 /* get the activities for all globally valid rows; 765 * the rows should be feasible, but due to numerical inaccuracies in the LP solver, they can be violated 766 */ 767 nviolrows = 0; 768 for( r = 0; r < nlprows; ++r ) 769 { 770 SCIP_ROW* row; 771 772 row = lprows[r]; 773 assert(SCIProwGetLPPos(row) == r); 774 775 if( !SCIProwIsLocal(row) ) 776 { 777 activities[r] = SCIPgetRowActivity(scip, row); 778 if( SCIPisFeasLT(scip, activities[r], SCIProwGetLhs(row)) 779 || SCIPisFeasGT(scip, activities[r], SCIProwGetRhs(row)) ) 780 { 781 violrows[nviolrows] = row; 782 violrowpos[r] = nviolrows; 783 nviolrows++; 784 } 785 else 786 violrowpos[r] = -1; 787 } 788 else 789 violrowpos[r] = -1; 790 } 791 792 nviolfracrows = 0; 793 /* calc the current number of fractional variables in rows */ 794 for( c = 0; c < nlpcands; ++c ) 795 addFracCounter(nfracsinrow, violrows, violrowpos, &nviolfracrows, nviolrows, nlprows, lpcands[c], +1); 796 797 /* get the working solution from heuristic's local data */ 798 sol = heurdata->sol; 799 assert(sol != NULL); 800 801 /* copy the current LP solution to the working solution */ 802 SCIP_CALL( SCIPlinkLPSol(scip, sol) ); 803 804 /* calculate the minimal objective value possible after rounding fractional variables */ 805 minobj = SCIPgetSolTransObj(scip, sol); 806 assert(minobj < SCIPgetCutoffbound(scip)); 807 for( c = 0; c < nlpcands; ++c ) 808 { 809 obj = SCIPvarGetObj(lpcands[c]); 810 bestshiftval = obj > 0.0 ? SCIPfeasFloor(scip, lpcandssol[c]) : SCIPfeasCeil(scip, lpcandssol[c]); 811 minobj += obj * (bestshiftval - lpcandssol[c]); 812 } 813 814 /* try to shift remaining variables in order to become/stay feasible */ 815 nnonimprovingshifts = 0; 816 minnviolrows = INT_MAX; 817 increaseweight = 1.0; 818 while( (nfrac > 0 || nviolrows > 0) && nnonimprovingshifts < MAXSHIFTINGS ) 819 { 820 SCIP_VAR* shiftvar; 821 SCIP_Real oldsolval; 822 SCIP_Real newsolval; 823 SCIP_Bool oldsolvalisfrac; 824 int probindex; 825 826 SCIPdebugMsg(scip, "shifting heuristic: nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g), cutoff=%g\n", 827 nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj), 828 SCIPretransformObj(scip, SCIPgetCutoffbound(scip))); 829 830 nprevviolrows = nviolrows; 831 832 /* choose next variable to process: 833 * - if a violated row exists, shift a variable decreasing the violation, that has least impact on other rows 834 * - otherwise, shift a variable, that has strongest devastating impact on rows in opposite direction 835 */ 836 shiftvar = NULL; 837 oldsolval = 0.0; 838 newsolval = 0.0; 839 if( nviolrows > 0 && (nfrac == 0 || nnonimprovingshifts < MAXSHIFTINGS-1) ) 840 { 841 SCIP_ROW* row; 842 int rowidx; 843 int rowpos; 844 int direction; 845 846 assert(nviolfracrows == 0 || nfrac > 0); 847 /* violated rows containing fractional variables are preferred; if such a row exists, choose the last one from the list 848 * (at position nviolfracrows - 1) because removing this row will cause one swapping operation less than other rows 849 */ 850 if( nviolfracrows > 0 ) 851 rowidx = nviolfracrows - 1; 852 else 853 /* there is no violated row containing a fractional variable, select a violated row uniformly at random */ 854 rowidx = SCIPrandomGetInt(heurdata->randnumgen, 0, nviolrows-1); 855 856 assert(0 <= rowidx && rowidx < nviolrows); 857 row = violrows[rowidx]; 858 rowpos = SCIProwGetLPPos(row); 859 assert(0 <= rowpos && rowpos < nlprows); 860 assert(violrowpos[rowpos] == rowidx); 861 assert(nfracsinrow[rowpos] == 0 || rowidx == nviolfracrows - 1); 862 863 SCIPdebugMsg(scip, "shifting heuristic: try to fix violated row <%s>: %g <= %g <= %g\n", 864 SCIProwGetName(row), SCIProwGetLhs(row), activities[rowpos], SCIProwGetRhs(row)); 865 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 866 867 /* get direction in which activity must be shifted */ 868 assert(SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row)) 869 || SCIPisFeasGT(scip, activities[rowpos], SCIProwGetRhs(row))); 870 direction = SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row)) ? +1 : -1; 871 872 /* search a variable that can shift the activity in the necessary direction */ 873 SCIP_CALL( selectShifting(scip, sol, row, activities[rowpos], direction, 874 nincreases, ndecreases, increaseweight, &shiftvar, &oldsolval, &newsolval) ); 875 } 876 877 if( shiftvar == NULL && nfrac > 0 ) 878 { 879 SCIPdebugMsg(scip, "shifting heuristic: search rounding variable and try to stay feasible\n"); 880 SCIP_CALL( selectEssentialRounding(scip, sol, minobj, lpcands, nlpcands, &shiftvar, &oldsolval, &newsolval) ); 881 } 882 883 /* check, whether shifting was possible */ 884 if( shiftvar == NULL || SCIPisEQ(scip, oldsolval, newsolval) ) 885 { 886 SCIPdebugMsg(scip, "shifting heuristic: -> didn't find a shifting variable\n"); 887 break; 888 } 889 890 SCIPdebugMsg(scip, "shifting heuristic: -> shift var <%s>[%g,%g], type=%d, oldval=%g, newval=%g, obj=%g\n", 891 SCIPvarGetName(shiftvar), SCIPvarGetLbGlobal(shiftvar), SCIPvarGetUbGlobal(shiftvar), SCIPvarGetType(shiftvar), 892 oldsolval, newsolval, SCIPvarGetObj(shiftvar)); 893 894 /* update row activities of globally valid rows */ 895 SCIP_CALL( updateActivities(scip, activities, violrows, violrowpos, &nviolrows, &nviolfracrows, nfracsinrow, nlprows, 896 shiftvar, oldsolval, newsolval) ); 897 if( nviolrows >= nprevviolrows ) 898 nnonimprovingshifts++; 899 else if( nviolrows < minnviolrows ) 900 { 901 minnviolrows = nviolrows; 902 nnonimprovingshifts = 0; 903 } 904 905 /* store new solution value and decrease fractionality counter */ 906 SCIP_CALL( SCIPsetSolVal(scip, sol, shiftvar, newsolval) ); 907 908 /* update fractionality counter and minimal objective value possible after shifting remaining variables */ 909 oldsolvalisfrac = !SCIPisFeasIntegral(scip, oldsolval) 910 && (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER); 911 obj = SCIPvarGetObj(shiftvar); 912 if( (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER) 913 && oldsolvalisfrac ) 914 { 915 assert(SCIPisFeasIntegral(scip, newsolval)); 916 nfrac--; 917 nnonimprovingshifts = 0; 918 minnviolrows = INT_MAX; 919 addFracCounter(nfracsinrow, violrows, violrowpos, &nviolfracrows, nviolrows, nlprows, shiftvar, -1); 920 921 /* the rounding was already calculated into the minobj -> update only if rounding in "wrong" direction */ 922 if( obj > 0.0 && newsolval > oldsolval ) 923 minobj += obj; 924 else if( obj < 0.0 && newsolval < oldsolval ) 925 minobj -= obj; 926 } 927 else 928 { 929 /* update minimal possible objective value */ 930 minobj += obj * (newsolval - oldsolval); 931 } 932 933 /* update increase/decrease arrays */ 934 if( !oldsolvalisfrac ) 935 { 936 probindex = SCIPvarGetProbindex(shiftvar); 937 assert(0 <= probindex && probindex < nvars); 938 increaseweight *= WEIGHTFACTOR; 939 if( newsolval < oldsolval ) 940 ndecreases[probindex] += increaseweight; 941 else 942 nincreases[probindex] += increaseweight; 943 if( increaseweight >= 1e+09 ) 944 { 945 int i; 946 947 for( i = 0; i < nvars; ++i ) 948 { 949 nincreases[i] /= increaseweight; 950 ndecreases[i] /= increaseweight; 951 } 952 increaseweight = 1.0; 953 } 954 } 955 956 SCIPdebugMsg(scip, "shifting heuristic: -> nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g)\n", 957 nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj)); 958 } 959 960 /* check, if the new solution is feasible */ 961 if( nfrac == 0 && nviolrows == 0 ) 962 { 963 SCIP_Bool stored; 964 965 /* check solution for feasibility, and add it to solution store if possible 966 * neither integrality nor feasibility of LP rows has to be checked, because this is already 967 * done in the shifting heuristic itself; however, we better check feasibility of LP rows, 968 * because of numerical problems with activity updating 969 */ 970 SCIP_CALL( SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, FALSE, TRUE, &stored) ); 971 972 if( stored ) 973 { 974 SCIPdebugMsg(scip, "found feasible shifted solution:\n"); 975 SCIPdebug( SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) ) ); 976 *result = SCIP_FOUNDSOL; 977 } 978 } 979 980 /* free memory buffers */ 981 SCIPfreeBufferArray(scip, &ndecreases); 982 SCIPfreeBufferArray(scip, &nincreases); 983 SCIPfreeBufferArray(scip, &nfracsinrow); 984 SCIPfreeBufferArray(scip, &violrowpos); 985 SCIPfreeBufferArray(scip, &violrows); 986 SCIPfreeBufferArray(scip, &activities); 987 988 return SCIP_OKAY; 989 } 990 991 992 /* 993 * heuristic specific interface methods 994 */ 995 996 /** creates the shifting heuristic with infeasibility recovering and includes it in SCIP */ 997 SCIP_RETCODE SCIPincludeHeurShifting( 998 SCIP* scip /**< SCIP data structure */ 999 ) 1000 { 1001 SCIP_HEUR* heur; 1002 1003 /* include primal heuristic */ 1004 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 1005 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 1006 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecShifting, NULL) ); 1007 1008 assert(heur != NULL); 1009 1010 /* set non-NULL pointers to callback methods */ 1011 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyShifting) ); 1012 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitShifting) ); 1013 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitShifting) ); 1014 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolShifting) ); 1015 1016 return SCIP_OKAY; 1017 } 1018 1019