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_zirounding.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief zirounding primal heuristic 28 * @author Gregor Hendel 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_zirounding.h" 35 #include "scip/pub_heur.h" 36 #include "scip/pub_lp.h" 37 #include "scip/pub_message.h" 38 #include "scip/pub_var.h" 39 #include "scip/scip_branch.h" 40 #include "scip/scip_heur.h" 41 #include "scip/scip_lp.h" 42 #include "scip/scip_mem.h" 43 #include "scip/scip_message.h" 44 #include "scip/scip_numerics.h" 45 #include "scip/scip_param.h" 46 #include "scip/scip_sol.h" 47 #include "scip/scip_solvingstats.h" 48 #include <string.h> 49 50 #define HEUR_NAME "zirounding" 51 #define HEUR_DESC "LP rounding heuristic as suggested by C. Wallace taking row slacks and bounds into account" 52 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_ROUNDING 53 #define HEUR_PRIORITY -500 54 #define HEUR_FREQ 1 55 #define HEUR_FREQOFS 0 56 #define HEUR_MAXDEPTH -1 57 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPNODE 58 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */ 59 60 #define DEFAULT_MAXROUNDINGLOOPS 2 /**< delimits the number of main loops, can be set to -1 for no limit */ 61 #define DEFAULT_STOPZIROUND TRUE /**< deactivation check is enabled by default */ 62 #define DEFAULT_STOPPERCENTAGE 0.02 /**< the tolerance percentage after which zirounding will not be executed anymore */ 63 #define DEFAULT_MINSTOPNCALLS 1000 /**< number of heuristic calls before deactivation check */ 64 65 #define MINSHIFT 1e-4 /**< minimum shift value for every single step */ 66 67 /* 68 * Data structures 69 */ 70 71 /** primal heuristic data */ 72 struct SCIP_HeurData 73 { 74 SCIP_SOL* sol; /**< working solution */ 75 SCIP_Longint lastlp; /**< the number of the last LP for which ZIRounding was called */ 76 int maxroundingloops; /**< limits rounding loops in execution */ 77 SCIP_Bool stopziround; /**< sets deactivation check */ 78 SCIP_Real stoppercentage; /**< threshold for deactivation check */ 79 int minstopncalls; /**< number of heuristic calls before deactivation check */ 80 }; 81 82 enum Direction 83 { 84 DIRECTION_UP = 1, 85 DIRECTION_DOWN = -1 86 }; 87 typedef enum Direction DIRECTION; 88 89 /* 90 * Local methods 91 */ 92 93 /** returns the fractionality of a value x, which is calculated as zivalue(x) = min(x-floor(x), ceil(x)-x) */ 94 static 95 SCIP_Real getZiValue( 96 SCIP* scip, /**< pointer to current SCIP data structure */ 97 SCIP_Real val /**< the value for which the fractionality should be computed */ 98 ) 99 { 100 SCIP_Real upgap; /* the gap between val and ceil(val) */ 101 SCIP_Real downgap; /* the gap between val and floor(val) */ 102 103 assert(scip != NULL); 104 105 upgap = SCIPfeasCeil(scip, val) - val; 106 downgap = val - SCIPfeasFloor(scip, val); 107 108 return MIN(upgap, downgap); 109 } 110 111 /** determines shifting bounds for variable */ 112 static 113 void calculateBounds( 114 SCIP* scip, /**< pointer to current SCIP data structure */ 115 SCIP_VAR* var, /**< the variable for which lb and ub have to be calculated */ 116 SCIP_Real currentvalue, /**< the current value of var in the working solution */ 117 SCIP_Real* upperbound, /**< pointer to store the calculated upper bound on the variable shift */ 118 SCIP_Real* lowerbound, /**< pointer to store the calculated lower bound on the variable shift */ 119 SCIP_Real* upslacks, /**< array that contains the slacks between row activities and the right hand sides of the rows */ 120 SCIP_Real* downslacks, /**< array that contains lhs slacks */ 121 int nslacks, /**< current number of slacks */ 122 SCIP_Bool* numericalerror /**< flag to determine whether a numerical error occurred */ 123 ) 124 { 125 SCIP_COL* col; 126 SCIP_ROW** colrows; 127 SCIP_Real* colvals; 128 int ncolvals; 129 int i; 130 131 assert(scip != NULL); 132 assert(var != NULL); 133 assert(upslacks != NULL); 134 assert(downslacks != NULL); 135 assert(upperbound != NULL); 136 assert(lowerbound != NULL); 137 138 /* get the column associated to the variable, the nonzero rows and the nonzero coefficients */ 139 col = SCIPvarGetCol(var); 140 colrows = SCIPcolGetRows(col); 141 colvals = SCIPcolGetVals(col); 142 ncolvals = SCIPcolGetNLPNonz(col); 143 144 /* only proceed, when variable has nonzero coefficients */ 145 if( ncolvals == 0 ) 146 return; 147 148 assert(colvals != NULL); 149 assert(colrows != NULL); 150 151 /* initialize the bounds on the shift to be the gap of the current solution value to the bounds of the variable */ 152 if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) ) 153 *upperbound = SCIPinfinity(scip); 154 else 155 *upperbound = SCIPvarGetUbGlobal(var) - currentvalue; 156 157 if( SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) ) 158 *lowerbound = SCIPinfinity(scip); 159 else 160 *lowerbound = currentvalue - SCIPvarGetLbGlobal(var); 161 162 /* go through every nonzero row coefficient corresponding to var to determine bounds for shifting 163 * in such a way that shifting maintains feasibility in every LP row. 164 * a lower or upper bound as it is calculated in zirounding always has to be >= 0.0. 165 * if one of these values is significantly < 0.0, this will cause the abort of execution of the heuristic so that 166 * infeasible solutions are avoided 167 */ 168 for( i = 0; i < ncolvals && MAX(*lowerbound, *upperbound) >= MINSHIFT; ++i ) 169 { 170 SCIP_ROW* row; 171 int rowpos; 172 173 row = colrows[i]; 174 rowpos = SCIProwGetLPPos(row); 175 176 /* the row might currently not be in the LP, ignore it! */ 177 if( rowpos == -1 ) 178 continue; 179 180 assert(0 <= rowpos && rowpos < nslacks); 181 182 /* all bounds and slacks as they are calculated in zirounding always have to be greater equal zero. 183 * It might however be due to numerical issues, e.g. with scaling, that they are not. Better abort in this case. 184 */ 185 if( SCIPisFeasLT(scip, *lowerbound, 0.0) || SCIPisFeasLT(scip, *upperbound, 0.0) 186 || SCIPisFeasLT(scip, upslacks[rowpos], 0.0) || SCIPisFeasLT(scip, downslacks[rowpos] , 0.0) ) 187 { 188 *numericalerror = TRUE; 189 return; 190 } 191 192 SCIPdebugMsg(scip, "colval: %15.8g, downslack: %15.8g, upslack: %5.2g, lb: %5.2g, ub: %5.2g\n", colvals[i], downslacks[rowpos], upslacks[rowpos], 193 *lowerbound, *upperbound); 194 195 /* if coefficient > 0, rounding up might violate up slack and rounding down might violate down slack 196 * thus search for the minimum so that no constraint is violated; vice versa for coefficient < 0 197 */ 198 if( colvals[i] > 0 ) 199 { 200 if( !SCIPisInfinity(scip, upslacks[rowpos]) ) 201 { 202 SCIP_Real upslack; 203 upslack = MAX(upslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */ 204 *upperbound = MIN(*upperbound, upslack/colvals[i]); 205 } 206 207 if( !SCIPisInfinity(scip, downslacks[rowpos]) ) 208 { 209 SCIP_Real downslack; 210 downslack = MAX(downslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */ 211 *lowerbound = MIN(*lowerbound, downslack/colvals[i]); 212 } 213 } 214 else 215 { 216 assert(colvals[i] != 0.0); 217 218 if( !SCIPisInfinity(scip, upslacks[rowpos]) ) 219 { 220 SCIP_Real upslack; 221 upslack = MAX(upslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */ 222 *lowerbound = MIN(*lowerbound, -upslack/colvals[i]); 223 } 224 225 if( !SCIPisInfinity(scip, downslacks[rowpos]) ) 226 { 227 SCIP_Real downslack; 228 downslack = MAX(downslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */ 229 *upperbound = MIN(*upperbound, -downslack/colvals[i]); 230 } 231 } 232 } 233 } 234 235 /** when a variable is shifted, the activities and slacks of all rows it appears in have to be updated */ 236 static 237 SCIP_RETCODE updateSlacks( 238 SCIP* scip, /**< pointer to current SCIP data structure */ 239 SCIP_SOL* sol, /**< working solution */ 240 SCIP_VAR* var, /**< pointer to variable to be modified */ 241 SCIP_Real shiftvalue, /**< the value by which the variable is shifted */ 242 SCIP_Real* upslacks, /**< upslacks of all rows the variable appears in */ 243 SCIP_Real* downslacks, /**< downslacks of all rows the variable appears in */ 244 SCIP_Real* activities, /**< activities of the LP rows */ 245 SCIP_VAR** slackvars, /**< the slack variables for equality rows */ 246 SCIP_Real* slackcoeffs, /**< the slack variable coefficients */ 247 int nslacks /**< size of the arrays */ 248 ) 249 { 250 SCIP_COL* col; /* the corresponding column of variable var */ 251 SCIP_ROW** rows; /* pointer to the nonzero coefficient rows for variable var */ 252 int nrows; /* the number of nonzeros */ 253 SCIP_Real* colvals; /* array to store the nonzero coefficients */ 254 int i; 255 256 assert(scip != NULL); 257 assert(sol != NULL); 258 assert(var != NULL); 259 assert(upslacks != NULL); 260 assert(downslacks != NULL); 261 assert(activities != NULL); 262 assert(nslacks >= 0); 263 264 col = SCIPvarGetCol(var); 265 assert(col != NULL); 266 267 rows = SCIPcolGetRows(col); 268 nrows = SCIPcolGetNLPNonz(col); 269 colvals = SCIPcolGetVals(col); 270 assert(nrows == 0 || (rows != NULL && colvals != NULL)); 271 272 /* go through all rows the shifted variable appears in */ 273 for( i = 0; i < nrows; ++i ) 274 { 275 int rowpos; 276 277 rowpos = SCIProwGetLPPos(rows[i]); 278 assert(-1 <= rowpos && rowpos < nslacks); 279 280 /* if the row is in the LP, update its activity, up and down slack */ 281 if( rowpos >= 0 ) 282 { 283 SCIP_Real val; 284 SCIP_Real lhs; 285 SCIP_Real rhs; 286 SCIP_ROW* row; 287 288 val = colvals[i] * shiftvalue; 289 row = rows[i]; 290 lhs = SCIProwGetLhs(row); 291 rhs = SCIProwGetRhs(row); 292 293 /* if the row is an equation, we update its slack variable instead of its activities */ 294 if( SCIPisFeasEQ(scip, lhs, rhs) ) 295 { 296 SCIP_Real slackvarshiftval; 297 SCIP_Real slackvarsolval; 298 299 assert(slackvars[rowpos] != NULL); 300 assert(!SCIPisZero(scip, slackcoeffs[rowpos])); 301 302 slackvarsolval = SCIPgetSolVal(scip, sol, slackvars[rowpos]); 303 slackvarshiftval = -val / slackcoeffs[rowpos]; 304 305 assert(SCIPisFeasGE(scip, slackvarsolval + slackvarshiftval, SCIPvarGetLbGlobal(slackvars[rowpos]))); 306 assert(SCIPisFeasLE(scip, slackvarsolval + slackvarshiftval, SCIPvarGetUbGlobal(slackvars[rowpos]))); 307 308 SCIP_CALL( SCIPsetSolVal(scip, sol, slackvars[rowpos], slackvarsolval + slackvarshiftval) ); 309 } 310 else if( !SCIPisInfinity(scip, REALABS(activities[rowpos])) ) 311 activities[rowpos] += val; 312 313 /* the slacks of the row now can be updated independently of its type */ 314 if( !SCIPisInfinity(scip, upslacks[rowpos]) ) 315 upslacks[rowpos] -= val; 316 if( !SCIPisInfinity(scip, -downslacks[rowpos]) ) 317 downslacks[rowpos] += val; 318 319 assert(SCIPisInfinity(scip, -lhs) || SCIPisFeasGE(scip, activities[rowpos], lhs)); 320 assert(SCIPisInfinity(scip, rhs) || SCIPisFeasLE(scip, activities[rowpos], rhs)); 321 } 322 } 323 return SCIP_OKAY; 324 } 325 326 /** finds a continuous slack variable for an equation row, NULL if none exists */ 327 static 328 void rowFindSlackVar( 329 SCIP* scip, /**< pointer to current SCIP data structure */ 330 SCIP_ROW* row, /**< the row for which a slack variable is searched */ 331 SCIP_VAR** varpointer, /**< pointer to store the slack variable */ 332 SCIP_Real* coeffpointer /**< pointer to store the coefficient of the slack variable */ 333 ) 334 { 335 int v; 336 SCIP_COL** rowcols; 337 SCIP_Real* rowvals; 338 int nrowvals; 339 340 assert(row != NULL); 341 assert(varpointer != NULL); 342 assert(coeffpointer != NULL); 343 344 rowcols = SCIProwGetCols(row); 345 rowvals = SCIProwGetVals(row); 346 nrowvals = SCIProwGetNNonz(row); 347 348 assert(nrowvals == 0 || rowvals != NULL); 349 assert(nrowvals == 0 || rowcols != NULL); 350 351 /* iterate over the row variables. Stop after the first unfixed continuous variable was found. */ 352 for( v = nrowvals - 1; v >= 0; --v ) 353 { 354 SCIP_VAR* colvar; 355 356 assert(rowcols[v] != NULL); 357 if( SCIPcolGetLPPos(rowcols[v]) == -1 ) 358 continue; 359 360 colvar = SCIPcolGetVar(rowcols[v]); 361 362 if( SCIPvarGetType(colvar) == SCIP_VARTYPE_CONTINUOUS 363 && !SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(colvar), SCIPvarGetUbGlobal(colvar)) 364 && SCIPcolGetNLPNonz(rowcols[v]) == 1 ) 365 { 366 SCIPdebugMsg(scip, " slack variable for row %s found: %s\n", SCIProwGetName(row), SCIPvarGetName(colvar)); 367 368 *coeffpointer = rowvals[v]; 369 *varpointer = colvar; 370 371 return; 372 } 373 } 374 375 *varpointer = NULL; 376 *coeffpointer = 0.0; 377 378 SCIPdebugMsg(scip, "No slack variable for row %s found. \n", SCIProwGetName(row)); 379 } 380 381 /* 382 * Callback methods of primal heuristic 383 */ 384 385 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 386 static 387 SCIP_DECL_HEURCOPY(heurCopyZirounding) 388 { /*lint --e{715}*/ 389 assert(scip != NULL); 390 assert(heur != NULL); 391 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 392 393 /* call inclusion method of primal heuristic */ 394 SCIP_CALL( SCIPincludeHeurZirounding(scip) ); 395 396 return SCIP_OKAY; 397 } 398 399 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 400 static 401 SCIP_DECL_HEURFREE(heurFreeZirounding) 402 { /*lint --e{715}*/ 403 SCIP_HEURDATA* heurdata; 404 405 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 406 407 heurdata = SCIPheurGetData(heur); 408 assert(heurdata != NULL); 409 410 /* free heuristic data */ 411 SCIPfreeBlockMemory(scip, &heurdata); 412 SCIPheurSetData(heur, NULL); 413 414 return SCIP_OKAY; 415 } 416 417 /** initialization method of primal heuristic (called after problem was transformed) */ 418 static 419 SCIP_DECL_HEURINIT(heurInitZirounding) 420 { /*lint --e{715}*/ 421 SCIP_HEURDATA* heurdata; 422 423 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 424 425 heurdata = SCIPheurGetData(heur); 426 assert(heurdata != NULL); 427 428 /* create working solution */ 429 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) ); 430 431 return SCIP_OKAY; 432 } 433 434 /** deinitialization method of primal heuristic (called before transformed problem is freed) */ 435 static 436 SCIP_DECL_HEUREXIT(heurExitZirounding) /*lint --e{715}*/ 437 { /*lint --e{715}*/ 438 SCIP_HEURDATA* heurdata; 439 440 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 441 442 heurdata = SCIPheurGetData(heur); 443 assert(heurdata != NULL); 444 445 /* free working solution */ 446 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) ); 447 448 return SCIP_OKAY; 449 } 450 451 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */ 452 static 453 SCIP_DECL_HEURINITSOL(heurInitsolZirounding) 454 { /*lint --e{715}*/ 455 SCIP_HEURDATA* heurdata; 456 457 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 458 459 heurdata = SCIPheurGetData(heur); 460 assert(heurdata != NULL); 461 462 heurdata->lastlp = -1; 463 464 return SCIP_OKAY; 465 } 466 467 468 /** execution method of primal heuristic */ 469 static 470 SCIP_DECL_HEUREXEC(heurExecZirounding) 471 { /*lint --e{715}*/ 472 SCIP_HEURDATA* heurdata; 473 SCIP_SOL* sol; 474 SCIP_VAR** lpcands; 475 SCIP_VAR** zilpcands; 476 477 SCIP_VAR** slackvars; 478 SCIP_Real* upslacks; 479 SCIP_Real* downslacks; 480 SCIP_Real* activities; 481 SCIP_Real* slackvarcoeffs; 482 SCIP_Bool* rowneedsslackvar; 483 484 SCIP_ROW** rows; 485 SCIP_Real* lpcandssol; 486 SCIP_Real* solarray; 487 488 SCIP_Longint nlps; 489 int currentlpcands; 490 int nlpcands; 491 int nimplfracs; 492 int i; 493 int c; 494 int nslacks; 495 int nroundings; 496 497 SCIP_Bool improvementfound; 498 SCIP_Bool numericalerror; 499 500 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 501 assert(result != NULL); 502 assert(SCIPhasCurrentNodeLP(scip)); 503 504 *result = SCIP_DIDNOTRUN; 505 506 /* do not call heuristic of node was already detected to be infeasible */ 507 if( nodeinfeasible ) 508 return SCIP_OKAY; 509 510 /* only call heuristic if an optimal LP-solution is at hand */ 511 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) 512 return SCIP_OKAY; 513 514 /* only call heuristic, if the LP objective value is smaller than the cutoff bound */ 515 if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) ) 516 return SCIP_OKAY; 517 518 /* get heuristic data */ 519 heurdata = SCIPheurGetData(heur); 520 assert(heurdata != NULL); 521 522 /* Do not call heuristic if deactivation check is enabled and percentage of found solutions in relation 523 * to number of calls falls below heurdata->stoppercentage */ 524 if( heurdata->stopziround && SCIPheurGetNCalls(heur) >= heurdata->minstopncalls 525 && SCIPheurGetNSolsFound(heur)/(SCIP_Real)SCIPheurGetNCalls(heur) < heurdata->stoppercentage ) 526 return SCIP_OKAY; 527 528 /* assure that heuristic has not already been called after the last LP had been solved */ 529 nlps = SCIPgetNLPs(scip); 530 if( nlps == heurdata->lastlp ) 531 return SCIP_OKAY; 532 533 heurdata->lastlp = nlps; 534 535 /* get fractional variables */ 536 SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, &nimplfracs) ); 537 nlpcands = nlpcands + nimplfracs; 538 /* make sure that there is at least one fractional variable that should be integral */ 539 if( nlpcands == 0 ) 540 return SCIP_OKAY; 541 542 assert(nlpcands > 0); 543 assert(lpcands != NULL); 544 assert(lpcandssol != NULL); 545 546 /* get LP rows data */ 547 rows = SCIPgetLPRows(scip); 548 nslacks = SCIPgetNLPRows(scip); 549 550 /* cannot do anything if LP is empty */ 551 if( nslacks == 0 ) 552 return SCIP_OKAY; 553 554 assert(rows != NULL); 555 assert(nslacks > 0); 556 557 /* get the working solution from heuristic's local data */ 558 sol = heurdata->sol; 559 assert(sol != NULL); 560 561 *result = SCIP_DIDNOTFIND; 562 563 solarray = NULL; 564 zilpcands = NULL; 565 566 /* copy the current LP solution to the working solution and allocate memory for local data */ 567 SCIP_CALL( SCIPlinkLPSol(scip, sol) ); 568 SCIP_CALL( SCIPallocBufferArray(scip, &solarray, nlpcands) ); 569 SCIP_CALL( SCIPallocBufferArray(scip, &zilpcands, nlpcands) ); 570 571 /* copy necessary data to local arrays */ 572 BMScopyMemoryArray(solarray, lpcandssol, nlpcands); 573 BMScopyMemoryArray(zilpcands, lpcands, nlpcands); 574 575 /* allocate buffer data arrays */ 576 SCIP_CALL( SCIPallocBufferArray(scip, &slackvars, nslacks) ); 577 SCIP_CALL( SCIPallocBufferArray(scip, &upslacks, nslacks) ); 578 SCIP_CALL( SCIPallocBufferArray(scip, &downslacks, nslacks) ); 579 SCIP_CALL( SCIPallocBufferArray(scip, &slackvarcoeffs, nslacks) ); 580 SCIP_CALL( SCIPallocBufferArray(scip, &rowneedsslackvar, nslacks) ); 581 SCIP_CALL( SCIPallocBufferArray(scip, &activities, nslacks) ); 582 583 BMSclearMemoryArray(slackvars, nslacks); 584 BMSclearMemoryArray(slackvarcoeffs, nslacks); 585 BMSclearMemoryArray(rowneedsslackvar, nslacks); 586 587 numericalerror = FALSE; 588 nroundings = 0; 589 590 /* loop over fractional variables and involved LP rows to find all rows which require a slack variable */ 591 for( c = 0; c < nlpcands; ++c ) 592 { 593 SCIP_VAR* cand; 594 SCIP_ROW** candrows; 595 int r; 596 int ncandrows; 597 598 cand = zilpcands[c]; 599 assert(cand != NULL); 600 assert(SCIPcolGetLPPos(SCIPvarGetCol(cand)) >= 0); 601 602 candrows = SCIPcolGetRows(SCIPvarGetCol(cand)); 603 ncandrows = SCIPcolGetNLPNonz(SCIPvarGetCol(cand)); 604 605 assert(candrows == NULL || ncandrows > 0); 606 607 for( r = 0; r < ncandrows; ++r ) 608 { 609 int rowpos; 610 611 assert(candrows != NULL); /* to please flexelint */ 612 assert(candrows[r] != NULL); 613 rowpos = SCIProwGetLPPos(candrows[r]); 614 615 if( rowpos >= 0 && SCIPisFeasEQ(scip, SCIProwGetLhs(candrows[r]), SCIProwGetRhs(candrows[r])) ) 616 { 617 rowneedsslackvar[rowpos] = TRUE; 618 SCIPdebugMsg(scip, " Row %s needs slack variable for variable %s\n", SCIProwGetName(candrows[r]), SCIPvarGetName(cand)); 619 } 620 } 621 } 622 623 /* calculate row slacks for every every row that belongs to the current LP and ensure, that the current solution 624 * has no violated constraint -- if any constraint is violated, i.e. a slack is significantly smaller than zero, 625 * this will cause the termination of the heuristic because Zirounding does not provide feasibility recovering 626 */ 627 for( i = 0; i < nslacks; ++i ) 628 { 629 SCIP_ROW* row; 630 SCIP_Real lhs; 631 SCIP_Real rhs; 632 633 row = rows[i]; 634 635 assert(row != NULL); 636 637 lhs = SCIProwGetLhs(row); 638 rhs = SCIProwGetRhs(row); 639 640 /* get row activity */ 641 activities[i] = SCIPgetRowActivity(scip, row); 642 assert(SCIPisFeasLE(scip, lhs, activities[i]) && SCIPisFeasLE(scip, activities[i], rhs)); 643 644 /* in special case if LHS or RHS is (-)infinity slacks have to be initialized as infinity */ 645 if( SCIPisInfinity(scip, -lhs) ) 646 downslacks[i] = SCIPinfinity(scip); 647 else 648 downslacks[i] = activities[i] - lhs; 649 650 if( SCIPisInfinity(scip, rhs) ) 651 upslacks[i] = SCIPinfinity(scip); 652 else 653 upslacks[i] = rhs - activities[i]; 654 655 SCIPdebugMsg(scip, "lhs:%5.2f <= act:%5.2g <= rhs:%5.2g --> down: %5.2g, up:%5.2g\n", lhs, activities[i], rhs, downslacks[i], upslacks[i]); 656 657 /* row is an equation. Try to find a slack variable in the row, i.e., 658 * a continuous variable which occurs only in this row. If no such variable exists, 659 * there is no hope for an IP-feasible solution in this round 660 */ 661 if( SCIPisFeasEQ(scip, lhs, rhs) && rowneedsslackvar[i] ) 662 { 663 /* @todo: This is only necessary for rows containing fractional variables. */ 664 rowFindSlackVar(scip, row, &(slackvars[i]), &(slackvarcoeffs[i])); 665 666 if( slackvars[i] == NULL ) 667 { 668 SCIPdebugMsg(scip, "No slack variable found for equation %s, terminating ZI Round heuristic\n", SCIProwGetName(row)); 669 goto TERMINATE; 670 } 671 else 672 { 673 SCIP_Real ubslackvar; 674 SCIP_Real lbslackvar; 675 SCIP_Real solvalslackvar; 676 SCIP_Real coeffslackvar; 677 SCIP_Real ubgap; 678 SCIP_Real lbgap; 679 680 assert(SCIPvarGetType(slackvars[i]) == SCIP_VARTYPE_CONTINUOUS); 681 solvalslackvar = SCIPgetSolVal(scip, sol, slackvars[i]); 682 ubslackvar = SCIPvarGetUbGlobal(slackvars[i]); 683 lbslackvar = SCIPvarGetLbGlobal(slackvars[i]); 684 685 coeffslackvar = slackvarcoeffs[i]; 686 assert(!SCIPisZero(scip, coeffslackvar)); 687 688 ubgap = MAX(0.0, ubslackvar - solvalslackvar); 689 lbgap = MAX(0.0, solvalslackvar - lbslackvar); 690 691 if( SCIPisPositive(scip, coeffslackvar) ) 692 { 693 if( !SCIPisInfinity(scip, lbslackvar) ) 694 upslacks[i] += coeffslackvar * lbgap; 695 else 696 upslacks[i] = SCIPinfinity(scip); 697 if( !SCIPisInfinity(scip, ubslackvar) ) 698 downslacks[i] += coeffslackvar * ubgap; 699 else 700 downslacks[i] = SCIPinfinity(scip); 701 } 702 else 703 { 704 if( !SCIPisInfinity(scip, ubslackvar) ) 705 upslacks[i] -= coeffslackvar * ubgap; 706 else 707 upslacks[i] = SCIPinfinity(scip); 708 if( !SCIPisInfinity(scip, lbslackvar) ) 709 downslacks[i] -= coeffslackvar * lbgap; 710 else 711 downslacks[i] = SCIPinfinity(scip); 712 } 713 SCIPdebugMsg(scip, " Slack variable for row %s at pos %d: %g <= %s = %g <= %g; Coeff %g, upslack = %g, downslack = %g \n", 714 SCIProwGetName(row), SCIProwGetLPPos(row), lbslackvar, SCIPvarGetName(slackvars[i]), solvalslackvar, ubslackvar, coeffslackvar, 715 upslacks[i], downslacks[i]); 716 } 717 } 718 /* due to numerical inaccuracies, the rows might be feasible, even if the slacks are 719 * significantly smaller than zero -> terminate 720 */ 721 if( SCIPisFeasLT(scip, upslacks[i], 0.0) || SCIPisFeasLT(scip, downslacks[i], 0.0) ) 722 goto TERMINATE; 723 } 724 725 assert(nslacks == 0 || (upslacks != NULL && downslacks != NULL && activities != NULL)); 726 727 /* initialize number of remaining variables and flag to enter the main loop */ 728 currentlpcands = nlpcands; 729 improvementfound = TRUE; 730 731 /* iterate over variables as long as there are fractional variables left */ 732 while( currentlpcands > 0 && improvementfound && (heurdata->maxroundingloops == -1 || nroundings < heurdata->maxroundingloops) ) 733 { /*lint --e{850}*/ 734 improvementfound = FALSE; 735 nroundings++; 736 SCIPdebugMsg(scip, "zirounding enters while loop for %d time with %d candidates left. \n", nroundings, currentlpcands); 737 738 /* check for every remaining fractional variable if a shifting decreases ZI-value of the variable */ 739 for( c = 0; c < currentlpcands; ++c ) 740 { 741 SCIP_VAR* var; 742 SCIP_Real oldsolval; 743 SCIP_Real upperbound; 744 SCIP_Real lowerbound; 745 SCIP_Real up; 746 SCIP_Real down; 747 SCIP_Real ziup; 748 SCIP_Real zidown; 749 SCIP_Real zicurrent; 750 SCIP_Real shiftval; 751 752 DIRECTION direction; 753 754 /* get values from local data */ 755 oldsolval = solarray[c]; 756 var = zilpcands[c]; 757 758 assert(!SCIPisFeasIntegral(scip, oldsolval)); 759 assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN); 760 761 /* calculate bounds for variable and make sure that there are no numerical inconsistencies */ 762 upperbound = SCIPinfinity(scip); 763 lowerbound = SCIPinfinity(scip); 764 calculateBounds(scip, var, oldsolval, &upperbound, &lowerbound, upslacks, downslacks, nslacks, &numericalerror); 765 766 if( numericalerror ) 767 goto TERMINATE; 768 769 /* continue if only marginal shifts are possible */ 770 if( MAX(upperbound, lowerbound) < MINSHIFT ) 771 { 772 /* stop immediately if a variable has not been rounded during the last round */ 773 if( nroundings == heurdata->maxroundingloops ) 774 break; 775 else 776 continue; 777 } 778 779 /* calculate the possible values after shifting */ 780 up = oldsolval + upperbound; 781 down = oldsolval - lowerbound; 782 783 /* if the variable is integer or implicit binary, do not shift further than the nearest integer */ 784 if( SCIPvarGetType(var) != SCIP_VARTYPE_BINARY) 785 { 786 SCIP_Real ceilx; 787 SCIP_Real floorx; 788 789 ceilx = SCIPfeasCeil(scip, oldsolval); 790 floorx = SCIPfeasFloor(scip, oldsolval); 791 up = MIN(up, ceilx); 792 down = MAX(down, floorx); 793 } 794 795 /* calculate necessary values */ 796 ziup = getZiValue(scip, up); 797 zidown = getZiValue(scip, down); 798 zicurrent = getZiValue(scip, oldsolval); 799 800 /* calculate the shifting direction that reduces ZI-value the most, 801 * if both directions improve ZI-value equally, take the direction which improves the objective 802 */ 803 if( SCIPisFeasLT(scip, zidown, zicurrent) || SCIPisFeasLT(scip, ziup, zicurrent) ) 804 { 805 if( SCIPisFeasEQ(scip,ziup, zidown) ) 806 direction = SCIPisFeasGE(scip, SCIPvarGetObj(var), 0.0) ? DIRECTION_DOWN : DIRECTION_UP; 807 else if( SCIPisFeasLT(scip, zidown, ziup) ) 808 direction = DIRECTION_DOWN; 809 else 810 direction = DIRECTION_UP; 811 812 /* once a possible shifting direction and value have been found, variable value is updated */ 813 shiftval = (direction == DIRECTION_UP ? up - oldsolval : down - oldsolval); 814 815 /* this improves numerical stability in some cases */ 816 if( direction == DIRECTION_UP ) 817 shiftval = MIN(shiftval, upperbound); 818 else 819 shiftval = MIN(shiftval, lowerbound); 820 /* update the solution */ 821 solarray[c] = direction == DIRECTION_UP ? up : down; 822 SCIP_CALL( SCIPsetSolVal(scip, sol, var, solarray[c]) ); 823 824 /* update the rows activities and slacks */ 825 SCIP_CALL( updateSlacks(scip, sol, var, shiftval, upslacks, 826 downslacks, activities, slackvars, slackvarcoeffs, nslacks) ); 827 828 SCIPdebugMsg(scip, "zirounding update step : %d var index, oldsolval=%g, shiftval=%g\n", 829 SCIPvarGetIndex(var), oldsolval, shiftval); 830 /* since at least one improvement has been found, heuristic will enter main loop for another time because the improvement 831 * might affect many LP rows and their current slacks and thus make further rounding steps possible */ 832 improvementfound = TRUE; 833 } 834 835 /* if solution value of variable has become feasibly integral due to rounding step, 836 * variable is put at the end of remaining candidates array so as not to be considered in future loops 837 */ 838 if( SCIPisFeasIntegral(scip, solarray[c]) ) 839 { 840 zilpcands[c] = zilpcands[currentlpcands - 1]; 841 solarray[c] = solarray[currentlpcands - 1]; 842 currentlpcands--; 843 844 /* counter is decreased if end of candidates array has not been reached yet */ 845 if( c < currentlpcands ) 846 c--; 847 } 848 else if( nroundings == heurdata->maxroundingloops ) 849 goto TERMINATE; 850 } 851 } 852 853 /* in case that no candidate is left for rounding after the final main loop 854 * the found solution has to be checked for feasibility in the original problem 855 */ 856 if( currentlpcands == 0 ) 857 { 858 SCIP_Bool stored; 859 SCIP_CALL(SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, TRUE, FALSE, &stored)); 860 if( stored ) 861 { 862 #ifdef SCIP_DEBUG 863 SCIPdebugMsg(scip, "found feasible rounded solution:\n"); 864 SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) ); 865 #endif 866 SCIPstatisticMessage(" ZI Round solution value: %g \n", SCIPgetSolOrigObj(scip, sol)); 867 868 *result = SCIP_FOUNDSOL; 869 } 870 } 871 872 /* free memory for all locally allocated data */ 873 TERMINATE: 874 SCIPfreeBufferArrayNull(scip, &activities); 875 SCIPfreeBufferArrayNull(scip, &rowneedsslackvar); 876 SCIPfreeBufferArrayNull(scip, &slackvarcoeffs); 877 SCIPfreeBufferArrayNull(scip, &downslacks); 878 SCIPfreeBufferArrayNull(scip, &upslacks); 879 SCIPfreeBufferArrayNull(scip, &slackvars); 880 SCIPfreeBufferArrayNull(scip, &zilpcands); 881 SCIPfreeBufferArrayNull(scip, &solarray); 882 883 return SCIP_OKAY; 884 } 885 886 /* 887 * primal heuristic specific interface methods 888 */ 889 890 /** creates the zirounding primal heuristic and includes it in SCIP */ 891 SCIP_RETCODE SCIPincludeHeurZirounding( 892 SCIP* scip /**< SCIP data structure */ 893 ) 894 { 895 SCIP_HEURDATA* heurdata; 896 SCIP_HEUR* heur; 897 898 /* create zirounding primal heuristic data */ 899 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 900 901 /* include primal heuristic */ 902 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 903 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 904 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecZirounding, heurdata) ); 905 906 assert(heur != NULL); 907 908 /* set non-NULL pointers to callback methods */ 909 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyZirounding) ); 910 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeZirounding) ); 911 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitZirounding) ); 912 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitZirounding) ); 913 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolZirounding) ); 914 915 /* add zirounding primal heuristic parameters */ 916 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/zirounding/maxroundingloops", 917 "determines maximum number of rounding loops", 918 &heurdata->maxroundingloops, TRUE, DEFAULT_MAXROUNDINGLOOPS, -1, INT_MAX, NULL, NULL) ); 919 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/zirounding/stopziround", 920 "flag to determine if Zirounding is deactivated after a certain percentage of unsuccessful calls", 921 &heurdata->stopziround, TRUE, DEFAULT_STOPZIROUND, NULL, NULL) ); 922 SCIP_CALL( SCIPaddRealParam(scip,"heuristics/zirounding/stoppercentage", 923 "if percentage of found solutions falls below this parameter, Zirounding will be deactivated", 924 &heurdata->stoppercentage, TRUE, DEFAULT_STOPPERCENTAGE, 0.0, 1.0, NULL, NULL) ); 925 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/zirounding/minstopncalls", 926 "determines the minimum number of calls before percentage-based deactivation of Zirounding is applied", 927 &heurdata->minstopncalls, TRUE, DEFAULT_MINSTOPNCALLS, 1, INT_MAX, NULL, NULL) ); 928 929 return SCIP_OKAY; 930 } 931