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_shiftandpropagate.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief shiftandpropagate primal heuristic 28 * @author Timo Berthold 29 * @author Gregor Hendel 30 */ 31 32 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 33 34 #include "blockmemshell/memory.h" 35 #include "scip/heur_shiftandpropagate.h" 36 #include "scip/pub_event.h" 37 #include "scip/pub_heur.h" 38 #include "scip/pub_lp.h" 39 #include "scip/pub_message.h" 40 #include "scip/pub_misc.h" 41 #include "scip/pub_misc_sort.h" 42 #include "scip/pub_sol.h" 43 #include "scip/pub_var.h" 44 #include "scip/scip_event.h" 45 #include "scip/scip_general.h" 46 #include "scip/scip_heur.h" 47 #include "scip/scip_lp.h" 48 #include "scip/scip_mem.h" 49 #include "scip/scip_message.h" 50 #include "scip/scip_numerics.h" 51 #include "scip/scip_param.h" 52 #include "scip/scip_prob.h" 53 #include "scip/scip_probing.h" 54 #include "scip/scip_randnumgen.h" 55 #include "scip/scip_sol.h" 56 #include "scip/scip_solvingstats.h" 57 #include "scip/scip_tree.h" 58 #include "scip/scip_var.h" 59 #include <string.h> 60 61 #define HEUR_NAME "shiftandpropagate" 62 #define HEUR_DESC "Pre-root heuristic to expand an auxiliary branch-and-bound tree and apply propagation techniques" 63 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_PROP 64 #define HEUR_PRIORITY 1000 65 #define HEUR_FREQ 0 66 #define HEUR_FREQOFS 0 67 #define HEUR_MAXDEPTH -1 68 #define HEUR_TIMING SCIP_HEURTIMING_BEFORENODE 69 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */ 70 71 #define DEFAULT_WEIGHT_INEQUALITY 1 /**< the heuristic row weight for inequalities */ 72 #define DEFAULT_WEIGHT_EQUALITY 3 /**< the heuristic row weight for equations */ 73 #define DEFAULT_RELAX TRUE /**< Should continuous variables be relaxed from the problem? */ 74 #define DEFAULT_PROBING TRUE /**< Is propagation of solution values enabled? */ 75 #define DEFAULT_ONLYWITHOUTSOL TRUE /**< Should heuristic only be executed if no primal solution was found, yet? */ 76 #define DEFAULT_NPROPROUNDS 10 /**< The default number of propagation rounds for each propagation used */ 77 #define DEFAULT_PROPBREAKER 65000 /**< fixed maximum number of propagations */ 78 #define DEFAULT_CUTOFFBREAKER 15 /**< fixed maximum number of allowed cutoffs before the heuristic stops */ 79 #define DEFAULT_RANDSEED 29 /**< the default random seed for random number generation */ 80 #define DEFAULT_SORTKEY 'v' /**< the default key for variable sorting */ 81 #define DEFAULT_SORTVARS TRUE /**< should variables be processed in sorted order? */ 82 #define DEFAULT_COLLECTSTATS TRUE /**< should variable statistics be collected during probing? */ 83 #define DEFAULT_STOPAFTERFEASIBLE TRUE /**< Should the heuristic stop calculating optimal shift values when no more rows are violated? */ 84 #define DEFAULT_PREFERBINARIES TRUE /**< Should binary variables be shifted first? */ 85 #define DEFAULT_SELECTBEST FALSE /**< should the heuristic choose the best candidate in every round? (set to FALSE for static order)? */ 86 #define DEFAULT_MAXCUTOFFQUOT 0.0 /**< maximum percentage of allowed cutoffs before stopping the heuristic */ 87 #define SORTKEYS "nrtuv"/**< options sorting key: (n)orms down, norms (u)p, (v)iolated rows decreasing, 88 * viola(t)ed rows increasing, or (r)andom */ 89 #define DEFAULT_NOZEROFIXING FALSE /**< should variables with a zero shifting value be delayed instead of being fixed? */ 90 #define DEFAULT_FIXBINLOCKS TRUE /**< should binary variables with no locks in one direction be fixed to that direction? */ 91 #define DEFAULT_BINLOCKSFIRST FALSE /**< should binary variables with no locks be preferred in the ordering? */ 92 #define DEFAULT_NORMALIZE TRUE /**< should coefficients and left/right hand sides be normalized by max row coeff? */ 93 #define DEFAULT_UPDATEWEIGHTS FALSE /**< should row weight be increased every time the row is violated? */ 94 #define DEFAULT_IMPLISCONTINUOUS TRUE /**< should implicit integer variables be treated as continuous variables? */ 95 #define DEFAULT_MINFIXINGRATELP 0.0 /**< minimum fixing rate over all variables (including continuous) to solve LP */ 96 97 #define EVENTHDLR_NAME "eventhdlrshiftandpropagate" 98 #define EVENTHDLR_DESC "event handler to catch bound changes" 99 #define EVENTTYPE_SHIFTANDPROPAGATE (SCIP_EVENTTYPE_BOUNDCHANGED | SCIP_EVENTTYPE_GBDCHANGED) 100 101 102 /* 103 * Data structures 104 */ 105 106 /** primal heuristic data */ 107 struct SCIP_HeurData 108 { 109 SCIP_COL** lpcols; /**< stores lp columns with discrete variables before cont. variables */ 110 SCIP_RANDNUMGEN* randnumgen; /**< random number generation */ 111 int* rowweights; /**< row weight storage */ 112 SCIP_Bool relax; /**< should continuous variables be relaxed from the problem */ 113 SCIP_Bool probing; /**< should probing be executed? */ 114 SCIP_Bool onlywithoutsol; /**< Should heuristic only be executed if no primal solution was found, yet? */ 115 int nlpcols; /**< the number of lp columns */ 116 int nproprounds; /**< The default number of propagation rounds for each propagation used */ 117 int cutoffbreaker; /**< the number of cutoffs before heuristic execution is stopped, or -1 for no 118 * limit */ 119 SCIP_EVENTHDLR* eventhdlr; /**< event handler to register and process variable bound changes */ 120 121 SCIP_Real maxcutoffquot; /**< maximum percentage of allowed cutoffs before stopping the heuristic */ 122 SCIP_Real minfixingratelp; /**< minimum fixing rate over all variables (including continuous) to solve LP */ 123 char sortkey; /**< the key by which variables are sorted */ 124 SCIP_Bool sortvars; /**< should variables be processed in sorted order? */ 125 SCIP_Bool collectstats; /**< should variable statistics be collected during probing? */ 126 SCIP_Bool stopafterfeasible; /**< Should the heuristic stop calculating optimal shift values when no 127 * more rows are violated? */ 128 SCIP_Bool preferbinaries; /**< Should binary variables be shifted first? */ 129 SCIP_Bool nozerofixing; /**< should variables with a zero shifting value be delayed instead of being fixed? */ 130 SCIP_Bool fixbinlocks; /**< should binary variables with no locks in one direction be fixed to that direction? */ 131 SCIP_Bool binlocksfirst; /**< should binary variables with no locks be preferred in the ordering? */ 132 SCIP_Bool normalize; /**< should coefficients and left/right hand sides be normalized by max row coeff? */ 133 SCIP_Bool updateweights; /**< should row weight be increased every time the row is violated? */ 134 SCIP_Bool impliscontinuous; /**< should implicit integer variables be treated as continuous variables? */ 135 SCIP_Bool selectbest; /**< should the heuristic choose the best candidate in every round? (set to FALSE for static order)? */ 136 SCIPstatistic( 137 SCIP_LPSOLSTAT lpsolstat; /**< the probing status after probing */ 138 SCIP_Longint ntotaldomredsfound; /**< the total number of domain reductions during heuristic */ 139 SCIP_Longint nlpiters; /**< number of LP iterations which the heuristic needed */ 140 int nremainingviols; /**< the number of remaining violations */ 141 int nprobings; /**< how many probings has the heuristic executed? */ 142 int ncutoffs; /**< has the probing node been cutoff? */ 143 ) 144 }; 145 146 /** status of a variable in heuristic transformation */ 147 enum TransformStatus 148 { 149 TRANSFORMSTATUS_NONE = 0, /**< variable has not been transformed yet */ 150 TRANSFORMSTATUS_LB = 1, /**< variable has been shifted by using lower bound (x-lb) */ 151 TRANSFORMSTATUS_NEG = 2, /**< variable has been negated by using upper bound (ub-x) */ 152 TRANSFORMSTATUS_FREE = 3 /**< variable does not have to be shifted */ 153 }; 154 typedef enum TransformStatus TRANSFORMSTATUS; 155 156 /** information about the matrix after its heuristic transformation */ 157 struct ConstraintMatrix 158 { 159 SCIP_Real* rowmatvals; /**< matrix coefficients row by row */ 160 int* rowmatind; /**< the indices of the corresponding variables */ 161 int* rowmatbegin; /**< the starting indices of each row */ 162 SCIP_Real* colmatvals; /**< matrix coefficients column by column */ 163 int* colmatind; /**< the indices of the corresponding rows for each coefficient */ 164 int* colmatbegin; /**< the starting indices of each column */ 165 int* violrows; /**< the number of violated rows for every variable */ 166 TRANSFORMSTATUS* transformstatus; /**< information about transform status of every discrete variable */ 167 SCIP_Real* lhs; /**< left hand side vector after normalization */ 168 SCIP_Real* rhs; /**< right hand side vector after normalization */ 169 SCIP_Real* colnorms; /**< vector norms of all discrete problem variables after normalization */ 170 SCIP_Real* upperbounds; /**< the upper bounds of every non-continuous variable after transformation*/ 171 SCIP_Real* transformshiftvals; /**< values by which original discrete variable bounds were shifted */ 172 int nnonzs; /**< number of nonzero column entries */ 173 int nrows; /**< number of rows of matrix */ 174 int ncols; /**< the number of columns in matrix (including continuous vars) */ 175 int ndiscvars; /**< number of discrete problem variables */ 176 SCIP_Bool normalized; /**< indicates if the matrix data has already been normalized */ 177 }; 178 typedef struct ConstraintMatrix CONSTRAINTMATRIX; 179 180 struct SCIP_EventhdlrData 181 { 182 CONSTRAINTMATRIX* matrix; /**< the constraint matrix of the heuristic */ 183 SCIP_HEURDATA* heurdata; /**< heuristic data */ 184 int* violatedrows; /**< all currently violated LP rows */ 185 int* violatedrowpos; /**< position in violatedrows array for every row */ 186 int* nviolatedrows; /**< pointer to the total number of currently violated rows */ 187 }; 188 189 struct SCIP_EventData 190 { 191 int colpos; /**< column position of the event-related variable */ 192 }; 193 /* 194 * Local methods 195 */ 196 197 /** returns whether a given variable is counted as discrete, depending on the parameter impliscontinuous */ 198 static 199 SCIP_Bool varIsDiscrete( 200 SCIP_VAR* var, /**< variable to check for discreteness */ 201 SCIP_Bool impliscontinuous /**< should implicit integer variables be counted as continuous? */ 202 ) 203 { 204 return SCIPvarIsIntegral(var) && (SCIPvarGetType(var) != SCIP_VARTYPE_IMPLINT || !impliscontinuous); 205 } 206 207 /** returns whether a given column is counted as discrete, depending on the parameter impliscontinuous */ 208 static 209 SCIP_Bool colIsDiscrete( 210 SCIP_COL* col, /**< column to check for discreteness */ 211 SCIP_Bool impliscontinuous /**< should implicit integer variables be counted as continuous? */ 212 ) 213 { 214 return SCIPcolIsIntegral(col) && (!impliscontinuous || SCIPvarGetType(SCIPcolGetVar(col)) != SCIP_VARTYPE_IMPLINT); 215 } 216 217 /** returns nonzero values and corresponding columns of given row */ 218 static 219 void getRowData( 220 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 221 int rowindex, /**< index of the desired row */ 222 SCIP_Real** valpointer, /**< pointer to store the nonzero coefficients of the row */ 223 SCIP_Real* lhs, /**< lhs of the row */ 224 SCIP_Real* rhs, /**< rhs of the row */ 225 int** indexpointer, /**< pointer to store column indices which belong to the nonzeros */ 226 int* nrowvals /**< pointer to store number of nonzeros in the desired row (or NULL) */ 227 ) 228 { 229 int arrayposition; 230 231 assert(matrix != NULL); 232 assert(0 <= rowindex && rowindex < matrix->nrows); 233 234 arrayposition = matrix->rowmatbegin[rowindex]; 235 236 if ( nrowvals != NULL ) 237 { 238 if( rowindex == matrix->nrows - 1 ) 239 *nrowvals = matrix->nnonzs - arrayposition; 240 else 241 *nrowvals = matrix->rowmatbegin[rowindex + 1] - arrayposition; /*lint !e679*/ 242 } 243 244 if( valpointer != NULL ) 245 *valpointer = &(matrix->rowmatvals[arrayposition]); 246 if( indexpointer != NULL ) 247 *indexpointer = &(matrix->rowmatind[arrayposition]); 248 249 if( lhs != NULL ) 250 *lhs = matrix->lhs[rowindex]; 251 252 if( rhs != NULL ) 253 *rhs = matrix->rhs[rowindex]; 254 } 255 256 /** returns nonzero values and corresponding rows of given column */ 257 static 258 void getColumnData( 259 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 260 int colindex, /**< the index of the desired column */ 261 SCIP_Real** valpointer, /**< pointer to store the nonzero coefficients of the column */ 262 int** indexpointer, /**< pointer to store row indices which belong to the nonzeros */ 263 int* ncolvals /**< pointer to store number of nonzeros in the desired column */ 264 ) 265 { 266 int arrayposition; 267 268 assert(matrix != NULL); 269 assert(0 <= colindex && colindex < matrix->ncols); 270 271 arrayposition = matrix->colmatbegin[colindex]; 272 273 if( ncolvals != NULL ) 274 { 275 if( colindex == matrix->ncols - 1 ) 276 *ncolvals = matrix->nnonzs - arrayposition; 277 else 278 *ncolvals = matrix->colmatbegin[colindex + 1] - arrayposition; /*lint !e679*/ 279 } 280 if( valpointer != NULL ) 281 *valpointer = &(matrix->colmatvals[arrayposition]); 282 283 if( indexpointer != NULL ) 284 *indexpointer = &(matrix->colmatind[arrayposition]); 285 } 286 287 /** relaxes a continuous variable from all its rows, which has influence 288 * on both the left and right hand side of the constraint. 289 */ 290 static 291 void relaxVar( 292 SCIP* scip, /**< current scip instance */ 293 SCIP_VAR* var, /**< variable which is relaxed from the problem */ 294 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 295 SCIP_Bool normalize /**< should coefficients and be normalized by rows maximum norms? */ 296 ) 297 { 298 SCIP_ROW** colrows; 299 SCIP_COL* varcol; 300 SCIP_Real* colvals; 301 SCIP_Real ub; 302 SCIP_Real lb; 303 int ncolvals; 304 int r; 305 306 assert(var != NULL); 307 assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN); 308 309 varcol = SCIPvarGetCol(var); 310 assert(varcol != NULL); 311 312 /* get nonzero values and corresponding rows of variable */ 313 colvals = SCIPcolGetVals(varcol); 314 ncolvals = SCIPcolGetNLPNonz(varcol); 315 colrows = SCIPcolGetRows(varcol); 316 317 ub = SCIPvarGetUbGlobal(var); 318 lb = SCIPvarGetLbGlobal(var); 319 320 assert(colvals != NULL || ncolvals == 0); 321 322 SCIPdebugMsg(scip, "Relaxing variable <%s> with lb <%g> and ub <%g>\n", 323 SCIPvarGetName(var), lb, ub); 324 325 assert(matrix->normalized); 326 /* relax variable from all its constraints */ 327 for( r = 0; r < ncolvals; ++r ) 328 { 329 SCIP_ROW* colrow; 330 SCIP_Real lhs; 331 SCIP_Real rhs; 332 SCIP_Real lhsvarbound; 333 SCIP_Real rhsvarbound; 334 SCIP_Real rowabs; 335 SCIP_Real colval; 336 int rowindex; 337 338 colrow = colrows[r]; 339 rowindex = SCIProwGetLPPos(colrow); 340 341 if( rowindex == -1 ) 342 break; 343 344 rowabs = SCIPgetRowMaxCoef(scip, colrow); 345 assert(colvals != NULL); /* to please flexelint */ 346 colval = colvals[r]; 347 if( normalize && SCIPisFeasGT(scip, rowabs, 0.0) ) 348 colval /= rowabs; 349 350 assert(0 <= rowindex && rowindex < matrix->nrows); 351 getRowData(matrix, rowindex, NULL, &lhs, &rhs, NULL, NULL); 352 /* variables bound influence the lhs and rhs of current row depending on the sign 353 * of the variables coefficient. 354 */ 355 if( SCIPisFeasPositive(scip, colval) ) 356 { 357 lhsvarbound = ub; 358 rhsvarbound = lb; 359 } 360 else if( SCIPisFeasNegative(scip, colval) ) 361 { 362 lhsvarbound = lb; 363 rhsvarbound = ub; 364 } 365 else 366 continue; 367 368 /* relax variable from the current row */ 369 if( !SCIPisInfinity(scip, -matrix->lhs[rowindex]) && !SCIPisInfinity(scip, ABS(lhsvarbound)) ) 370 matrix->lhs[rowindex] -= colval * lhsvarbound; 371 else 372 matrix->lhs[rowindex] = -SCIPinfinity(scip); 373 374 if( !SCIPisInfinity(scip, matrix->rhs[rowindex]) && !SCIPisInfinity(scip, ABS(rhsvarbound)) ) 375 matrix->rhs[rowindex] -= colval * rhsvarbound; 376 else 377 matrix->rhs[rowindex] = SCIPinfinity(scip); 378 379 SCIPdebugMsg(scip, "Row <%s> changed:Coefficient <%g>, LHS <%g> --> <%g>, RHS <%g> --> <%g>\n", 380 SCIProwGetName(colrow), colval, lhs, matrix->lhs[rowindex], rhs, matrix->rhs[rowindex]); 381 } /*lint !e438*/ 382 } 383 384 /** transforms bounds of a given variable s.t. its lower bound equals zero afterwards. 385 * If the variable already has lower bound zero, the variable is not transformed, 386 * if not, the variable's bounds are changed w.r.t. the smaller absolute value of its 387 * bounds in order to avoid numerical inaccuracies. If both lower and upper bound 388 * of the variable differ from infinity, there are two cases. If |lb| <= |ub|, 389 * the bounds are shifted by -lb, else a new variable ub - x replaces x. 390 * The transformation is memorized by the transform status of the variable s.t. 391 * retransformation is possible. 392 */ 393 static 394 void transformVariable( 395 SCIP* scip, /**< current scip instance */ 396 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 397 SCIP_HEURDATA* heurdata, /**< heuristic data */ 398 int colpos /**< position of variable column in matrix */ 399 ) 400 { 401 SCIP_COL* col; 402 SCIP_VAR* var; 403 SCIP_Real lb; 404 SCIP_Real ub; 405 406 SCIP_Bool negatecoeffs; /* do the row coefficients need to be negated? */ 407 SCIP_Real deltashift; /* difference from previous transformation */ 408 409 assert(matrix != NULL); 410 assert(0 <= colpos && colpos < heurdata->nlpcols); 411 col = heurdata->lpcols[colpos]; 412 assert(col != NULL); 413 assert(SCIPcolIsInLP(col)); 414 415 var = SCIPcolGetVar(col); 416 assert(var != NULL); 417 assert(SCIPvarIsIntegral(var)); 418 lb = SCIPvarGetLbLocal(var); 419 ub = SCIPvarGetUbLocal(var); 420 421 negatecoeffs = FALSE; 422 /* if both lower and upper bound are -infinity and infinity, resp., this is reflected by a free transform status. 423 * If the lower bound is already zero, this is reflected by identity transform status. In both cases, none of the 424 * corresponding rows needs to be modified. 425 */ 426 if( SCIPisInfinity(scip, -lb) && SCIPisInfinity(scip, ub) ) 427 { 428 if( matrix->transformstatus[colpos] == TRANSFORMSTATUS_NEG ) 429 negatecoeffs = TRUE; 430 431 deltashift = matrix->transformshiftvals[colpos]; 432 matrix->transformshiftvals[colpos] = 0.0; 433 matrix->transformstatus[colpos] = TRANSFORMSTATUS_FREE; 434 } 435 else if( SCIPisLE(scip, REALABS(lb), REALABS(ub)) ) 436 { 437 assert(!SCIPisInfinity(scip, REALABS(lb))); 438 439 matrix->transformstatus[colpos] = TRANSFORMSTATUS_LB; 440 deltashift = lb; 441 matrix->transformshiftvals[colpos] = lb; 442 } 443 else 444 { 445 assert(!SCIPisInfinity(scip, ub)); 446 if( matrix->transformstatus[colpos] != TRANSFORMSTATUS_NEG ) 447 negatecoeffs = TRUE; 448 matrix->transformstatus[colpos] = TRANSFORMSTATUS_NEG; 449 deltashift = ub; 450 matrix->transformshiftvals[colpos] = ub; 451 } 452 453 /* determine the upper bound for this variable in heuristic transformation (lower bound is implicit; always 0) */ 454 if( !SCIPisInfinity(scip, ub) && !SCIPisInfinity(scip, lb) ) 455 matrix->upperbounds[colpos] = MIN(ub - lb, SCIPinfinity(scip)); /*lint !e666*/ 456 else 457 matrix->upperbounds[colpos] = SCIPinfinity(scip); 458 459 /* a real transformation is necessary. The variable x is either shifted by -lb or 460 * replaced by ub - x, depending on the smaller absolute of lb and ub. 461 */ 462 if( !SCIPisFeasZero(scip, deltashift) || negatecoeffs ) 463 { 464 SCIP_Real* vals; 465 int* rows; 466 int nrows; 467 int i; 468 469 assert(!SCIPisInfinity(scip, deltashift)); 470 471 /* get nonzero values and corresponding rows of column */ 472 getColumnData(matrix, colpos, &vals, &rows, &nrows); 473 assert(nrows == 0 ||(vals != NULL && rows != NULL)); 474 475 /* go through rows and modify its lhs, rhs and the variable coefficient, if necessary */ 476 for( i = 0; i < nrows; ++i ) 477 { 478 int rowpos = rows[i]; 479 assert(rowpos >= 0); 480 assert(rowpos < matrix->nrows); 481 482 if( !SCIPisInfinity(scip, -(matrix->lhs[rowpos])) ) 483 matrix->lhs[rowpos] -= (vals[i]) * deltashift; 484 485 if( !SCIPisInfinity(scip, matrix->rhs[rowpos]) ) 486 matrix->rhs[rowpos] -= (vals[i]) * deltashift; 487 488 if( negatecoeffs ) 489 (vals[i]) = -(vals[i]); 490 491 assert(SCIPisFeasLE(scip, matrix->lhs[rowpos], matrix->rhs[rowpos])); 492 } 493 } 494 SCIPdebugMsg(scip, "Variable <%s> at colpos %d transformed. Status %d LB <%g> --> <%g>, UB <%g> --> <%g>\n", 495 SCIPvarGetName(var), colpos, matrix->transformstatus[colpos], lb, 0.0, ub, matrix->upperbounds[colpos]); 496 } 497 498 /** initializes copy of the original coefficient matrix and applies heuristic specific adjustments: normalizing row 499 * vectors, transforming variable domains such that lower bound is zero, and relaxing continuous variables. 500 */ 501 static 502 SCIP_RETCODE initMatrix( 503 SCIP* scip, /**< current scip instance */ 504 CONSTRAINTMATRIX* matrix, /**< constraint matrix object to be initialized */ 505 SCIP_HEURDATA* heurdata, /**< heuristic data */ 506 int* colposs, /**< position of columns according to variable type sorting */ 507 SCIP_Bool normalize, /**< should coefficients and be normalized by rows maximum norms? */ 508 int* nmaxrows, /**< maximum number of rows a variable appears in */ 509 SCIP_Bool relax, /**< should continuous variables be relaxed from the problem? */ 510 SCIP_Bool* initialized, /**< was the initialization successful? */ 511 SCIP_Bool* infeasible /**< is the problem infeasible? */ 512 ) 513 { 514 SCIP_ROW** lprows; 515 SCIP_COL** lpcols; 516 SCIP_Bool impliscontinuous; 517 int i; 518 int j; 519 int currentpointer; 520 521 int nrows; 522 int ncols; 523 524 assert(scip != NULL); 525 assert(matrix != NULL); 526 assert(initialized!= NULL); 527 assert(infeasible != NULL); 528 assert(nmaxrows != NULL); 529 530 SCIPdebugMsg(scip, "entering Matrix Initialization method of SHIFTANDPROPAGATE heuristic!\n"); 531 532 /* get LP row data; column data is already initialized in heurdata */ 533 SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nrows) ); 534 lpcols = heurdata->lpcols; 535 ncols = heurdata->nlpcols; 536 537 matrix->nrows = nrows; 538 matrix->nnonzs = 0; 539 matrix->normalized = FALSE; 540 matrix->ndiscvars = 0; 541 *nmaxrows = 0; 542 impliscontinuous = heurdata->impliscontinuous; 543 544 /* count the number of nonzeros of the LP constraint matrix */ 545 for( j = 0; j < ncols; ++j ) 546 { 547 assert(lpcols[j] != NULL); 548 assert(SCIPcolGetLPPos(lpcols[j]) >= 0); 549 550 if( colIsDiscrete(lpcols[j], impliscontinuous) ) 551 { 552 matrix->nnonzs += SCIPcolGetNLPNonz(lpcols[j]); 553 ++matrix->ndiscvars; 554 } 555 } 556 557 matrix->ncols = matrix->ndiscvars; 558 559 if( matrix->nnonzs == 0 ) 560 { 561 SCIPdebugMsg(scip, "No matrix entries - Terminating initialization of matrix.\n"); 562 563 *initialized = FALSE; 564 565 return SCIP_OKAY; 566 } 567 568 /* allocate memory for the members of heuristic matrix */ 569 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatvals, matrix->nnonzs) ); 570 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatind, matrix->nnonzs) ); 571 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatvals, matrix->nnonzs) ); 572 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatind, matrix->nnonzs) ); 573 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rowmatbegin, matrix->nrows) ); 574 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colmatbegin, matrix->ncols) ); 575 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->lhs, matrix->nrows) ); 576 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->rhs, matrix->nrows) ); 577 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->colnorms, matrix->ncols) ); 578 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->violrows, matrix->ncols) ); 579 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->transformstatus, matrix->ndiscvars) ); 580 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->upperbounds, matrix->ndiscvars) ); 581 SCIP_CALL( SCIPallocBufferArray(scip, &matrix->transformshiftvals, matrix->ndiscvars) ); 582 583 /* set transform status of variables */ 584 for( j = 0; j < matrix->ndiscvars; ++j ) 585 matrix->transformstatus[j] = TRANSFORMSTATUS_NONE; 586 587 currentpointer = 0; 588 *infeasible = FALSE; 589 590 /* initialize the rows vector of the heuristic matrix together with its corresponding 591 * lhs, rhs. 592 */ 593 for( i = 0; i < nrows; ++i ) 594 { 595 SCIP_COL** cols; 596 SCIP_ROW* row; 597 SCIP_Real* rowvals; 598 SCIP_Real constant; 599 SCIP_Real maxval; 600 int nrowlpnonz; 601 602 /* get LP row information */ 603 row = lprows[i]; 604 rowvals = SCIProwGetVals(row); 605 nrowlpnonz = SCIProwGetNLPNonz(row); 606 maxval = SCIPgetRowMaxCoef(scip, row); 607 cols = SCIProwGetCols(row); 608 constant = SCIProwGetConstant(row); 609 610 SCIPdebugMsg(scip, " %s : maxval=%g \n", SCIProwGetName(row), maxval); 611 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 612 assert(!SCIPisInfinity(scip, constant)); 613 614 matrix->rowmatbegin[i] = currentpointer; 615 616 /* modify the lhs and rhs w.r.t to the rows constant and normalize by 1-norm, i.e divide the lhs and rhs by the 617 * maximum absolute value of the row 618 */ 619 if( !SCIPisInfinity(scip, -SCIProwGetLhs(row)) ) 620 matrix->lhs[i] = SCIProwGetLhs(row) - constant; 621 else 622 matrix->lhs[i] = -SCIPinfinity(scip); 623 624 if( !SCIPisInfinity(scip, SCIProwGetRhs(row)) ) 625 matrix->rhs[i] = SCIProwGetRhs(row) - constant; 626 else 627 matrix->rhs[i] = SCIPinfinity(scip); 628 629 /* make sure that maxval is larger than zero before normalization. 630 * Maxval may be zero if the constraint contains no variables but is modifiable, hence not redundant 631 */ 632 if( normalize && !SCIPisFeasZero(scip, maxval) ) 633 { 634 if( !SCIPisInfinity(scip, -matrix->lhs[i]) ) 635 matrix->lhs[i] /= maxval; 636 if( !SCIPisInfinity(scip, matrix->rhs[i]) ) 637 matrix->rhs[i] /= maxval; 638 } 639 640 /* in case of empty rows with a 0 < lhs <= 0.0 or 0.0 <= rhs < 0 we deduce the infeasibility of the problem */ 641 if( nrowlpnonz == 0 && (SCIPisFeasPositive(scip, matrix->lhs[i]) || SCIPisFeasNegative(scip, matrix->rhs[i])) ) 642 { 643 *infeasible = TRUE; 644 SCIPdebugMsg(scip, " Matrix initialization stopped because of row infeasibility! \n"); 645 break; 646 } 647 648 /* row coefficients are normalized and copied to heuristic matrix */ 649 for( j = 0; j < nrowlpnonz; ++j ) 650 { 651 if( !colIsDiscrete(cols[j], impliscontinuous) ) 652 continue; 653 assert(SCIPcolGetLPPos(cols[j]) >= 0); 654 assert(currentpointer < matrix->nnonzs); 655 656 matrix->rowmatvals[currentpointer] = rowvals[j]; 657 if( normalize && SCIPisFeasGT(scip, maxval, 0.0) ) 658 matrix->rowmatvals[currentpointer] /= maxval; 659 660 matrix->rowmatind[currentpointer] = colposs[SCIPcolGetLPPos(cols[j])]; 661 662 ++currentpointer; 663 } 664 } 665 666 matrix->normalized = TRUE; 667 668 if( *infeasible ) 669 return SCIP_OKAY; 670 671 assert(currentpointer == matrix->nnonzs); 672 673 currentpointer = 0; 674 675 /* copy the nonzero coefficient data column by column to heuristic matrix */ 676 for( j = 0; j < matrix->ncols; ++j ) 677 { 678 SCIP_COL* currentcol; 679 SCIP_ROW** rows; 680 SCIP_Real* colvals; 681 int ncolnonz; 682 683 assert(SCIPcolGetLPPos(lpcols[j]) >= 0); 684 685 currentcol = lpcols[j]; 686 assert(colIsDiscrete(currentcol, impliscontinuous)); 687 688 colvals = SCIPcolGetVals(currentcol); 689 rows = SCIPcolGetRows(currentcol); 690 ncolnonz = SCIPcolGetNLPNonz(currentcol); 691 matrix->colnorms[j] = ncolnonz; 692 693 *nmaxrows = MAX(*nmaxrows, ncolnonz); 694 695 /* loop over all rows with nonzero coefficients in the column, transform them and add them to the heuristic matrix */ 696 matrix->colmatbegin[j] = currentpointer; 697 698 for( i = 0; i < ncolnonz; ++i ) 699 { 700 SCIP_Real maxval; 701 702 assert(rows[i] != NULL); 703 assert(0 <= SCIProwGetLPPos(rows[i])); 704 assert(SCIProwGetLPPos(rows[i]) < nrows); 705 assert(currentpointer < matrix->nnonzs); 706 707 /* rows are normalized by maximum norm */ 708 maxval = SCIPgetRowMaxCoef(scip, rows[i]); 709 710 assert(maxval > 0); 711 712 matrix->colmatvals[currentpointer] = colvals[i]; 713 if( normalize && SCIPisFeasGT(scip, maxval, 0.0) ) 714 matrix->colmatvals[currentpointer] /= maxval; 715 716 matrix->colmatind[currentpointer] = SCIProwGetLPPos(rows[i]); 717 718 /* update the column norm */ 719 matrix->colnorms[j] += ABS(matrix->colmatvals[currentpointer]); 720 ++currentpointer; 721 } 722 } 723 assert(currentpointer == matrix->nnonzs); 724 725 /* each variable is either transformed, if it supposed to be integral, or relaxed */ 726 for( j = 0; j < (relax ? ncols : matrix->ndiscvars); ++j ) 727 { 728 SCIP_COL* col; 729 730 col = lpcols[j]; 731 if( colIsDiscrete(col, impliscontinuous) ) 732 { 733 matrix->transformshiftvals[j] = 0.0; 734 transformVariable(scip, matrix, heurdata, j); 735 } 736 else 737 { 738 SCIP_VAR* var; 739 var = SCIPcolGetVar(col); 740 assert(!varIsDiscrete(var, impliscontinuous)); 741 relaxVar(scip, var, matrix, normalize); 742 } 743 } 744 *initialized = TRUE; 745 746 SCIPdebugMsg(scip, "Matrix initialized for %d discrete variables with %d cols, %d rows and %d nonzero entries\n", 747 matrix->ndiscvars, matrix->ncols, matrix->nrows, matrix->nnonzs); 748 return SCIP_OKAY; 749 } 750 751 /** frees all members of the heuristic matrix */ 752 static 753 void freeMatrix( 754 SCIP* scip, /**< current SCIP instance */ 755 CONSTRAINTMATRIX** matrix /**< constraint matrix object */ 756 ) 757 { 758 assert(scip != NULL); 759 assert(matrix != NULL); 760 761 /* all fields are only allocated, if problem is not empty */ 762 if( (*matrix)->nnonzs > 0 ) 763 { 764 assert((*matrix) != NULL); 765 assert((*matrix)->rowmatbegin != NULL); 766 assert((*matrix)->rowmatvals != NULL); 767 assert((*matrix)->rowmatind != NULL); 768 assert((*matrix)->colmatbegin != NULL); 769 assert((*matrix)->colmatvals!= NULL); 770 assert((*matrix)->colmatind != NULL); 771 assert((*matrix)->lhs != NULL); 772 assert((*matrix)->rhs != NULL); 773 assert((*matrix)->transformstatus != NULL); 774 assert((*matrix)->transformshiftvals != NULL); 775 776 /* free all fields */ 777 SCIPfreeBufferArray(scip, &((*matrix)->transformshiftvals)); 778 SCIPfreeBufferArray(scip, &((*matrix)->upperbounds)); 779 SCIPfreeBufferArray(scip, &((*matrix)->transformstatus)); 780 SCIPfreeBufferArray(scip, &((*matrix)->violrows)); 781 SCIPfreeBufferArray(scip, &((*matrix)->colnorms)); 782 SCIPfreeBufferArray(scip, &((*matrix)->rhs)); 783 SCIPfreeBufferArray(scip, &((*matrix)->lhs)); 784 SCIPfreeBufferArray(scip, &((*matrix)->colmatbegin)); 785 SCIPfreeBufferArray(scip, &((*matrix)->rowmatbegin)); 786 SCIPfreeBufferArray(scip, &((*matrix)->colmatind)); 787 SCIPfreeBufferArray(scip, &((*matrix)->colmatvals)); 788 SCIPfreeBufferArray(scip, &((*matrix)->rowmatind)); 789 SCIPfreeBufferArray(scip, &((*matrix)->rowmatvals)); 790 791 (*matrix)->nrows = 0; 792 (*matrix)->ncols = 0; 793 } 794 795 /* free matrix */ 796 SCIPfreeBuffer(scip, matrix); 797 } 798 799 /** updates the information about a row whenever violation status changes */ 800 static 801 void checkRowViolation( 802 SCIP* scip, /**< current SCIP instance */ 803 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 804 int rowindex, /**< index of the row */ 805 int* violatedrows, /**< contains all violated rows */ 806 int* violatedrowpos, /**< positions of rows in the violatedrows array */ 807 int* nviolatedrows, /**< pointer to update total number of violated rows */ 808 int* rowweights, /**< row weight storage */ 809 SCIP_Bool updateweights /**< should row weight be increased every time the row is violated? */ 810 ) 811 { 812 int* cols; 813 int ncols; 814 int c; 815 int violadd; 816 assert(matrix != NULL); 817 assert(violatedrows != NULL); 818 assert(violatedrowpos != NULL); 819 assert(nviolatedrows != NULL); 820 821 getRowData(matrix, rowindex, NULL, NULL, NULL, &cols, &ncols); 822 violadd = 0; 823 824 /* row is now violated. Enqueue it in the set of violated rows. */ 825 if( violatedrowpos[rowindex] == -1 && (SCIPisFeasGT(scip, matrix->lhs[rowindex], 0.0) || SCIPisFeasLT(scip, matrix->rhs[rowindex], 0.0)) ) 826 { 827 assert(*nviolatedrows < matrix->nrows); 828 829 violatedrows[*nviolatedrows] = rowindex; 830 violatedrowpos[rowindex] = *nviolatedrows; 831 ++(*nviolatedrows); 832 if( updateweights ) 833 ++rowweights[rowindex]; 834 835 violadd = 1; 836 } 837 /* row is now feasible. Remove it from the set of violated rows. */ 838 else if( violatedrowpos[rowindex] >= 0 && SCIPisFeasLE(scip, matrix->lhs[rowindex], 0.0) && SCIPisFeasGE(scip, matrix->rhs[rowindex], 0.0) ) 839 { 840 /* swap the row with last violated row */ 841 if( violatedrowpos[rowindex] != *nviolatedrows - 1 ) 842 { 843 assert(*nviolatedrows - 1 >= 0); 844 violatedrows[violatedrowpos[rowindex]] = violatedrows[*nviolatedrows - 1]; 845 violatedrowpos[violatedrows[*nviolatedrows - 1]] = violatedrowpos[rowindex]; 846 } 847 848 /* unlink the row from its position in the array and decrease number of violated rows */ 849 violatedrowpos[rowindex] = -1; 850 --(*nviolatedrows); 851 violadd = -1; 852 } 853 854 /* increase or decrease the column violation counter */ 855 for( c = 0; c < ncols; ++c ) 856 { 857 matrix->violrows[cols[c]] += violadd; 858 assert(matrix->violrows[cols[c]] >= 0); 859 } 860 } 861 862 /** collects the necessary information about row violations for the zero-solution. That is, 863 * all solution values in heuristic transformation are zero. 864 */ 865 static 866 void checkViolations( 867 SCIP* scip, /**< current scip instance */ 868 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 869 int colidx, /**< column index for specific column, or -1 for all rows */ 870 int* violatedrows, /**< violated rows */ 871 int* violatedrowpos, /**< row positions of violated rows */ 872 int* nviolatedrows, /**< pointer to store the number of violated rows */ 873 int* rowweights, /**< weight array for every row */ 874 SCIP_Bool updateweights /**< should row weight be increased every time the row is violated? */ 875 ) 876 { 877 int nrows; 878 int* rowindices; 879 int i; 880 881 assert(matrix != NULL); 882 assert(violatedrows != NULL); 883 assert(violatedrowpos != NULL); 884 assert(nviolatedrows != NULL); 885 assert(-1 <= colidx && colidx < matrix->ncols); 886 887 /* check if we requested an update for a single variable, or if we want to (re)-initialize the whole violation info */ 888 if( colidx >= 0 ) 889 getColumnData(matrix, colidx, NULL, &rowindices, &nrows); 890 else 891 { 892 nrows = matrix->nrows; 893 rowindices = NULL; 894 *nviolatedrows = 0; 895 896 /* reinitialize the violated rows */ 897 for( i = 0; i < nrows; ++i ) 898 violatedrowpos[i] = -1; 899 900 /* clear the violated row counters for all variables */ 901 BMSclearMemoryArray(matrix->violrows, matrix->ndiscvars); 902 } 903 904 assert(colidx < 0 || *nviolatedrows >= 0); 905 SCIPdebugMsg(scip, "Entering violation check for %d rows! \n", nrows); 906 /* loop over rows and check if it is violated */ 907 for( i = 0; i < nrows; ++i ) 908 { 909 int rowpos; 910 if( colidx >= 0 ) 911 { 912 assert(rowindices != NULL); 913 rowpos = rowindices[i]; 914 } 915 else 916 rowpos = i; 917 /* check, if zero solution violates this row */ 918 checkRowViolation(scip, matrix, rowpos, violatedrows, violatedrowpos, nviolatedrows, rowweights, updateweights); 919 920 assert((violatedrowpos[rowpos] == -1 && SCIPisFeasGE(scip, matrix->rhs[rowpos], 0.0) && SCIPisFeasLE(scip, matrix->lhs[rowpos], 0.0)) 921 || (violatedrowpos[rowpos] >= 0 &&(SCIPisFeasLT(scip, matrix->rhs[rowpos], 0.0) || SCIPisFeasGT(scip, matrix->lhs[rowpos], 0.0)))); 922 } 923 } 924 925 /** retransforms solution values of variables according to their transformation status */ 926 static 927 SCIP_Real retransformVariable( 928 SCIP* scip, /**< current scip instance */ 929 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 930 SCIP_VAR* var, /**< variable whose solution value has to be retransformed */ 931 int varindex, /**< permutation of variable indices according to sorting */ 932 SCIP_Real solvalue /**< solution value of the variable */ 933 ) 934 { 935 TRANSFORMSTATUS status; 936 937 assert(matrix != NULL); 938 assert(var != NULL); 939 940 status = matrix->transformstatus[varindex]; 941 assert(status != TRANSFORMSTATUS_NONE); 942 943 /* check if original variable has different bounds and transform solution value correspondingly */ 944 if( status == TRANSFORMSTATUS_LB ) 945 { 946 assert(!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var))); 947 948 return solvalue + matrix->transformshiftvals[varindex]; 949 } 950 else if( status == TRANSFORMSTATUS_NEG ) 951 { 952 assert(!SCIPisInfinity(scip, SCIPvarGetUbLocal(var))); 953 return matrix->transformshiftvals[varindex] - solvalue; 954 } 955 return solvalue; 956 } 957 958 /** determines the best shifting value of a variable 959 * @todo if there is already an incumbent solution, try considering the objective cutoff as additional constraint */ 960 static 961 SCIP_RETCODE getOptimalShiftingValue( 962 SCIP* scip, /**< current scip instance */ 963 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 964 int varindex, /**< index of variable which should be shifted */ 965 int direction, /**< the direction for this variable */ 966 int* rowweights, /**< weighting of rows for best shift calculation */ 967 SCIP_Real* steps, /**< buffer array to store the individual steps for individual rows */ 968 int* violationchange, /**< buffer array to store the individual change of feasibility of row */ 969 SCIP_Real* beststep, /**< pointer to store optimal shifting step */ 970 int* rowviolations /**< pointer to store new weighted sum of row violations, i.e, v - f */ 971 ) 972 { 973 SCIP_Real* vals; 974 int* rows; 975 976 SCIP_Real slacksurplus; 977 SCIP_Real upperbound; 978 979 int nrows; 980 int sum; 981 int i; 982 983 SCIP_Bool allzero; 984 985 assert(beststep != NULL); 986 assert(rowviolations != NULL); 987 assert(rowweights != NULL); 988 assert(steps != NULL); 989 assert(violationchange != NULL); 990 assert(direction == 1 || direction == -1); 991 992 upperbound = matrix->upperbounds[varindex]; 993 994 /* get nonzero values and corresponding rows of variable */ 995 getColumnData(matrix, varindex, &vals, &rows, &nrows); 996 997 /* loop over rows and calculate, which is the minimum shift to make this row feasible 998 * or the minimum shift to violate this row 999 */ 1000 allzero = TRUE; 1001 slacksurplus = 0.0; 1002 for( i = 0; i < nrows; ++i ) 1003 { 1004 SCIP_Real lhs; 1005 SCIP_Real rhs; 1006 SCIP_Real val; 1007 int rowpos; 1008 SCIP_Bool rowisviolated; 1009 int rowweight; 1010 1011 /* get the row data */ 1012 rowpos = rows[i]; 1013 assert(rowpos >= 0); 1014 lhs = matrix->lhs[rowpos]; 1015 rhs = matrix->rhs[rowpos]; 1016 rowweight = rowweights[rowpos]; 1017 val = direction * vals[i]; 1018 1019 /* determine if current row is violated or not */ 1020 rowisviolated =(SCIPisFeasLT(scip, rhs, 0.0) || SCIPisFeasLT(scip, -lhs, 0.0)); 1021 1022 /* for a feasible row, determine the minimum integer value within the bounds of the variable by which it has to be 1023 * shifted to make row infeasible. 1024 */ 1025 if( !rowisviolated ) 1026 { 1027 SCIP_Real maxfeasshift; 1028 1029 maxfeasshift = SCIPinfinity(scip); 1030 1031 /* feasibility can only be violated if the variable has a lock in the corresponding direction, 1032 * i.e. a positive coefficient for a "<="-constraint, a negative coefficient for a ">="-constraint. 1033 */ 1034 if( SCIPisFeasGT(scip, val, 0.0) && !SCIPisInfinity(scip, rhs) ) 1035 maxfeasshift = SCIPfeasFloor(scip, rhs/val); 1036 else if( SCIPisFeasLT(scip, val, 0.0) && !SCIPisInfinity(scip, -lhs) ) 1037 maxfeasshift = SCIPfeasFloor(scip, lhs/val); 1038 1039 /* if the variable has no lock in the current row, it can still help to increase the slack of this row; 1040 * we measure slack increase for shifting by one 1041 */ 1042 if( SCIPisFeasGT(scip, val, 0.0) && SCIPisInfinity(scip, rhs) ) 1043 slacksurplus += val; 1044 if( SCIPisFeasLT(scip, val, 0.0) && SCIPisInfinity(scip, -lhs) ) 1045 slacksurplus -= val; 1046 1047 /* check if the least violating shift lies within variable bounds and set corresponding array values */ 1048 if( !SCIPisInfinity(scip, maxfeasshift) && SCIPisFeasLE(scip, maxfeasshift + 1.0, upperbound) ) 1049 { 1050 steps[i] = maxfeasshift + 1.0; 1051 violationchange[i] = rowweight; 1052 allzero = FALSE; 1053 } 1054 else 1055 { 1056 steps[i] = upperbound; 1057 violationchange[i] = 0; 1058 } 1059 } 1060 /* for a violated row, determine the minimum integral value within the bounds of the variable by which it has to be 1061 * shifted to make row feasible. 1062 */ 1063 else 1064 { 1065 SCIP_Real minfeasshift; 1066 1067 minfeasshift = SCIPinfinity(scip); 1068 1069 /* if coefficient has the right sign to make row feasible, determine the minimum integer to shift variable 1070 * to obtain feasibility 1071 */ 1072 if( SCIPisFeasLT(scip, -lhs, 0.0) && SCIPisFeasGT(scip, val, 0.0) ) 1073 minfeasshift = SCIPfeasCeil(scip, lhs/val); 1074 else if( SCIPisFeasLT(scip, rhs,0.0) && SCIPisFeasLT(scip, val, 0.0) ) 1075 minfeasshift = SCIPfeasCeil(scip, rhs/val); 1076 1077 /* check if the minimum feasibility recovery shift lies within variable bounds and set corresponding array 1078 * values 1079 */ 1080 if( !SCIPisInfinity(scip, minfeasshift) && SCIPisFeasLE(scip, minfeasshift, upperbound) ) 1081 { 1082 steps[i] = minfeasshift; 1083 violationchange[i] = -rowweight; 1084 allzero = FALSE; 1085 } 1086 else 1087 { 1088 steps[i] = upperbound; 1089 violationchange[i] = 0; 1090 } 1091 } 1092 } 1093 1094 /* in case that the variable cannot affect the feasibility of any row, in particular it cannot violate 1095 * a single row, but we can add slack to already feasible rows, we will do this 1096 */ 1097 if( allzero ) 1098 { 1099 if( ! SCIPisInfinity(scip, upperbound) && SCIPisGT(scip, slacksurplus, 0.0) ) 1100 *beststep = direction * upperbound; 1101 else 1102 *beststep = 0.0; 1103 1104 return SCIP_OKAY; 1105 } 1106 1107 /* sorts rows by increasing value of steps */ 1108 SCIPsortRealInt(steps, violationchange, nrows); 1109 1110 *beststep = 0.0; 1111 *rowviolations = 0; 1112 sum = 0; 1113 1114 /* best shifting step is calculated by summing up the violation changes for each relevant step and 1115 * taking the one which leads to the minimum sum. This sum measures the balance of feasibility recovering and 1116 * violating changes which will be obtained by shifting the variable by this step 1117 * note, the sums for smaller steps have to be taken into account for all bigger steps, i.e., the sums can be 1118 * computed iteratively 1119 */ 1120 for( i = 0; i < nrows && !SCIPisInfinity(scip, steps[i]); ++i ) 1121 { 1122 sum += violationchange[i]; 1123 1124 /* if we reached the last entry for the current step value, we have finished computing its sum and 1125 * update the step defining the minimum sum 1126 */ 1127 if( (i == nrows-1 || steps[i+1] > steps[i]) && sum < *rowviolations ) /*lint !e679*/ 1128 { 1129 *rowviolations = sum; 1130 *beststep = direction * steps[i]; 1131 } 1132 } 1133 assert(*rowviolations <= 0); 1134 assert(!SCIPisInfinity(scip, *beststep)); 1135 1136 return SCIP_OKAY; 1137 } 1138 1139 /** updates transformation of a given variable by taking into account current local bounds. if the bounds have changed 1140 * since last update, updating the heuristic specific upper bound of the variable, its current transformed solution value 1141 * and all affected rows is necessary. 1142 */ 1143 static 1144 SCIP_RETCODE updateTransformation( 1145 SCIP* scip, /**< current scip */ 1146 CONSTRAINTMATRIX* matrix, /**< constraint matrix object */ 1147 SCIP_HEURDATA* heurdata, /**< heuristic data */ 1148 int varindex, /**< index of variable in matrix */ 1149 SCIP_Real lb, /**< local lower bound of the variable */ 1150 SCIP_Real ub, /**< local upper bound of the variable */ 1151 int* violatedrows, /**< violated rows */ 1152 int* violatedrowpos, /**< violated row positions */ 1153 int* nviolatedrows /**< pointer to store number of violated rows */ 1154 ) 1155 { 1156 TRANSFORMSTATUS status; 1157 SCIP_Real deltashift; 1158 SCIP_Bool checkviolations; 1159 1160 assert(scip != NULL); 1161 assert(matrix != NULL); 1162 assert(0 <= varindex && varindex < matrix->ndiscvars); 1163 1164 /* deltashift is the difference between the old and new transformation value. */ 1165 deltashift = 0.0; 1166 status = matrix->transformstatus[varindex]; 1167 1168 SCIPdebugMsg(scip, " Variable <%d> [%g,%g], status %d(%g), ub %g \n", varindex, lb, ub, status, 1169 matrix->transformshiftvals[varindex], matrix->upperbounds[varindex]); 1170 1171 checkviolations = FALSE; 1172 /* depending on the variable status, deltashift is calculated differently. */ 1173 switch( status ) 1174 { 1175 case TRANSFORMSTATUS_LB: 1176 if( SCIPisInfinity(scip, -lb) ) 1177 { 1178 transformVariable(scip, matrix, heurdata, varindex); 1179 checkviolations = TRUE; 1180 } 1181 else 1182 { 1183 deltashift = lb - (matrix->transformshiftvals[varindex]); 1184 matrix->transformshiftvals[varindex] = lb; 1185 if( !SCIPisInfinity(scip, ub) ) 1186 matrix->upperbounds[varindex] = ub - lb; 1187 else 1188 matrix->upperbounds[varindex] = SCIPinfinity(scip); 1189 } 1190 break; 1191 case TRANSFORMSTATUS_NEG: 1192 if( SCIPisInfinity(scip, ub) ) 1193 { 1194 transformVariable(scip, matrix, heurdata, varindex); 1195 checkviolations = TRUE; 1196 } 1197 else 1198 { 1199 deltashift = (matrix->transformshiftvals[varindex]) - ub; 1200 matrix->transformshiftvals[varindex] = ub; 1201 1202 if( !SCIPisInfinity(scip, -lb) ) 1203 matrix->upperbounds[varindex] = MIN(ub - lb, SCIPinfinity(scip)); /*lint !e666*/ 1204 else 1205 matrix->upperbounds[varindex] = SCIPinfinity(scip); 1206 } 1207 break; 1208 case TRANSFORMSTATUS_FREE: 1209 /* in case of a free transform status, if one of the bounds has become finite, we want 1210 * to transform this variable to a variable with a lowerbound or a negated transform status */ 1211 if( !SCIPisInfinity(scip, -lb) || !SCIPisInfinity(scip, ub) ) 1212 { 1213 transformVariable(scip, matrix, heurdata, varindex); 1214 1215 /* violations have to be rechecked for rows in which variable appears */ 1216 checkviolations = TRUE; 1217 1218 assert(matrix->transformstatus[varindex] == TRANSFORMSTATUS_LB || matrix->transformstatus[varindex] == TRANSFORMSTATUS_NEG); 1219 assert(SCIPisLE(scip, ABS(lb), ABS(ub)) || matrix->transformstatus[varindex] == TRANSFORMSTATUS_NEG); 1220 } 1221 break; 1222 1223 case TRANSFORMSTATUS_NONE: 1224 default: 1225 SCIPerrorMessage("Error: Invalid variable status <%d> in shift and propagagate heuristic, aborting!\n", status); 1226 SCIPABORT(); 1227 return SCIP_INVALIDDATA; /*lint !e527*/ 1228 } 1229 /* if the bound, by which the variable was shifted, has changed, deltashift is different from zero, which requires 1230 * an update of all affected rows 1231 */ 1232 if( !SCIPisFeasZero(scip, deltashift) ) 1233 { 1234 int i; 1235 int* rows; 1236 SCIP_Real* vals; 1237 int nrows; 1238 1239 /* get nonzero values and corresponding rows of variable */ 1240 getColumnData(matrix, varindex, &vals, &rows, &nrows); 1241 1242 /* go through rows, update the rows w.r.t. the influence of the changed transformation of the variable */ 1243 for( i = 0; i < nrows; ++i ) 1244 { 1245 SCIPdebugMsg(scip, " update slacks of row<%d>: coefficient <%g>, %g <= 0 <= %g \n", 1246 rows[i], vals[i], matrix->lhs[rows[i]], matrix->rhs[rows[i]]); 1247 1248 if( !SCIPisInfinity(scip, -(matrix->lhs[rows[i]])) ) 1249 matrix->lhs[rows[i]] -= (vals[i]) * deltashift; 1250 1251 if( !SCIPisInfinity(scip, matrix->rhs[rows[i]]) ) 1252 matrix->rhs[rows[i]] -= (vals[i]) * deltashift; 1253 } 1254 checkviolations = TRUE; 1255 } 1256 1257 /* check and update information about violated rows, if necessary */ 1258 if( checkviolations ) 1259 checkViolations(scip, matrix, varindex, violatedrows, violatedrowpos, nviolatedrows, heurdata->rowweights, heurdata->updateweights); 1260 1261 SCIPdebugMsg(scip, " Variable <%d> [%g,%g], status %d(%g), ub %g \n", varindex, lb, ub, status, 1262 matrix->transformshiftvals[varindex], matrix->upperbounds[varindex]); 1263 1264 return SCIP_OKAY; 1265 } 1266 1267 /** comparison method for columns; binary < integer < implicit < continuous variables */ 1268 static 1269 SCIP_DECL_SORTPTRCOMP(heurSortColsShiftandpropagate) 1270 { 1271 SCIP_COL* col1; 1272 SCIP_COL* col2; 1273 SCIP_VAR* var1; 1274 SCIP_VAR* var2; 1275 SCIP_VARTYPE vartype1; 1276 SCIP_VARTYPE vartype2; 1277 int key1; 1278 int key2; 1279 1280 col1 = (SCIP_COL*)elem1; 1281 col2 = (SCIP_COL*)elem2; 1282 var1 = SCIPcolGetVar(col1); 1283 var2 = SCIPcolGetVar(col2); 1284 assert(var1 != NULL && var2 != NULL); 1285 1286 vartype1 = SCIPvarGetType(var1); 1287 vartype2 = SCIPvarGetType(var2); 1288 1289 switch (vartype1) 1290 { 1291 case SCIP_VARTYPE_BINARY: 1292 key1 = 1; 1293 break; 1294 case SCIP_VARTYPE_INTEGER: 1295 key1 = 2; 1296 break; 1297 case SCIP_VARTYPE_IMPLINT: 1298 key1 = 3; 1299 break; 1300 case SCIP_VARTYPE_CONTINUOUS: 1301 key1 = 4; 1302 break; 1303 default: 1304 key1 = -1; 1305 SCIPerrorMessage("unknown variable type\n"); 1306 SCIPABORT(); 1307 break; 1308 } 1309 switch (vartype2) 1310 { 1311 case SCIP_VARTYPE_BINARY: 1312 key2 = 1; 1313 break; 1314 case SCIP_VARTYPE_INTEGER: 1315 key2 = 2; 1316 break; 1317 case SCIP_VARTYPE_IMPLINT: 1318 key2 = 3; 1319 break; 1320 case SCIP_VARTYPE_CONTINUOUS: 1321 key2 = 4; 1322 break; 1323 default: 1324 key2 = -1; 1325 SCIPerrorMessage("unknown variable type\n"); 1326 SCIPABORT(); 1327 break; 1328 } 1329 return key1 - key2; 1330 } 1331 1332 /* 1333 * Callback methods of primal heuristic 1334 */ 1335 1336 /** deinitialization method of primal heuristic(called before transformed problem is freed) */ 1337 static 1338 SCIP_DECL_HEUREXIT(heurExitShiftandpropagate) 1339 { /*lint --e{715}*/ 1340 SCIP_HEURDATA* heurdata; 1341 1342 heurdata = SCIPheurGetData(heur); 1343 assert(heurdata != NULL); 1344 1345 /* free random number generator */ 1346 SCIPfreeRandom(scip, &heurdata->randnumgen); 1347 1348 /* if statistic mode is enabled, statistics are printed to console */ 1349 SCIPstatistic( 1350 SCIPstatisticMessage( 1351 " DETAILS : %d violations left, %d probing status\n", 1352 heurdata->nremainingviols, 1353 heurdata->lpsolstat 1354 ); 1355 SCIPstatisticMessage( 1356 " SHIFTANDPROPAGATE PROBING : %d probings, %" SCIP_LONGINT_FORMAT " domain reductions, ncutoffs: %d , LP iterations: %" SCIP_LONGINT_FORMAT " \n ", 1357 heurdata->nprobings, 1358 heurdata->ntotaldomredsfound, 1359 heurdata->ncutoffs, 1360 heurdata->nlpiters); 1361 ); 1362 1363 return SCIP_OKAY; 1364 } 1365 1366 /** initialization method of primal heuristic(called after problem was transformed). We only need this method for 1367 * statistic mode of heuristic. 1368 */ 1369 static 1370 SCIP_DECL_HEURINIT(heurInitShiftandpropagate) 1371 { /*lint --e{715}*/ 1372 SCIP_HEURDATA* heurdata; 1373 1374 heurdata = SCIPheurGetData(heur); 1375 1376 assert(heurdata != NULL); 1377 1378 /* create random number generator */ 1379 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, 1380 DEFAULT_RANDSEED, TRUE) ); 1381 1382 SCIPstatistic( 1383 heurdata->lpsolstat = SCIP_LPSOLSTAT_NOTSOLVED; 1384 heurdata->nremainingviols = 0; 1385 heurdata->nprobings = 0; 1386 heurdata->ntotaldomredsfound = 0; 1387 heurdata->ncutoffs = 0; 1388 heurdata->nlpiters = 0; 1389 ) 1390 return SCIP_OKAY; 1391 } 1392 1393 /** destructor of primal heuristic to free user data(called when SCIP is exiting) */ 1394 static 1395 SCIP_DECL_HEURFREE(heurFreeShiftandpropagate) 1396 { /*lint --e{715}*/ 1397 SCIP_HEURDATA* heurdata; 1398 SCIP_EVENTHDLR* eventhdlr; 1399 SCIP_EVENTHDLRDATA* eventhdlrdata; 1400 1401 heurdata = SCIPheurGetData(heur); 1402 assert(heurdata != NULL); 1403 eventhdlr = heurdata->eventhdlr; 1404 assert(eventhdlr != NULL); 1405 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 1406 1407 SCIPfreeBlockMemoryNull(scip, &eventhdlrdata); 1408 1409 /* free heuristic data */ 1410 SCIPfreeBlockMemory(scip, &heurdata); 1411 1412 SCIPheurSetData(heur, NULL); 1413 1414 return SCIP_OKAY; 1415 } 1416 1417 1418 /** copy method for primal heuristic plugins(called when SCIP copies plugins) */ 1419 static 1420 SCIP_DECL_HEURCOPY(heurCopyShiftandpropagate) 1421 { /*lint --e{715}*/ 1422 assert(scip != NULL); 1423 assert(heur != NULL); 1424 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 1425 1426 /* call inclusion method of primal heuristic */ 1427 SCIP_CALL( SCIPincludeHeurShiftandpropagate(scip) ); 1428 1429 return SCIP_OKAY; 1430 } 1431 1432 /** execution method of primal heuristic */ 1433 static 1434 SCIP_DECL_HEUREXEC(heurExecShiftandpropagate) 1435 { /*lint --e{715}*/ 1436 SCIP_HEURDATA* heurdata; /* heuristic data */ 1437 SCIP_EVENTHDLR* eventhdlr; /* shiftandpropagate event handler */ 1438 SCIP_EVENTHDLRDATA* eventhdlrdata; /* event handler data */ 1439 SCIP_EVENTDATA** eventdatas; /* event data for every variable */ 1440 1441 CONSTRAINTMATRIX* matrix; /* constraint matrix object */ 1442 SCIP_COL** lpcols; /* lp columns */ 1443 SCIP_SOL* sol; /* solution pointer */ 1444 SCIP_Real* colnorms; /* contains Euclidean norms of column vectors */ 1445 1446 SCIP_Real* steps; /* buffer arrays for best shift selection in main loop */ 1447 int* violationchange; 1448 1449 int* violatedrows; /* the violated rows */ 1450 int* violatedrowpos; /* the array position of a violated row, or -1 */ 1451 int* permutation; /* reflects the position of the variables after sorting */ 1452 int* violatedvarrows; /* number of violated rows for each variable */ 1453 int* colposs; /* position of columns according to variable type sorting */ 1454 int nlpcols; /* number of lp columns */ 1455 int nviolatedrows; /* number of violated rows */ 1456 int ndiscvars; /* number of non-continuous variables of the problem */ 1457 int lastindexofsusp; /* last variable which has been swapped due to a cutoff */ 1458 int nbinvars; /* number of binary variables */ 1459 #ifndef NDEBUG 1460 int nintvars = 0; /* number of integer variables */ 1461 #endif 1462 int i; 1463 int r; 1464 int v; 1465 int c; 1466 int ncutoffs; /* counts the number of cutoffs for this execution */ 1467 int nprobings; /* counts the number of probings */ 1468 int nlprows; /* the number LP rows */ 1469 int nmaxrows; /* maximum number of LP rows of a variable */ 1470 1471 SCIP_Bool initialized; /* has the matrix been initialized? */ 1472 SCIP_Bool cutoff; /* has current probing node been cutoff? */ 1473 SCIP_Bool probing; /* should probing be applied or not? */ 1474 SCIP_Bool infeasible; /* FALSE as long as currently infeasible rows have variables left */ 1475 SCIP_Bool impliscontinuous; 1476 1477 heurdata = SCIPheurGetData(heur); 1478 assert(heurdata != NULL); 1479 1480 eventhdlr = heurdata->eventhdlr; 1481 assert(eventhdlr != NULL); 1482 1483 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 1484 assert(eventhdlrdata != NULL); 1485 1486 *result = SCIP_DIDNOTRUN; 1487 SCIPdebugMsg(scip, "entering execution method of shift and propagate heuristic\n"); 1488 1489 /* heuristic is obsolete if there are only continuous variables */ 1490 if( SCIPgetNVars(scip) - SCIPgetNContVars(scip) == 0 ) 1491 return SCIP_OKAY; 1492 1493 /* stop execution method if there is already a primarily feasible solution at hand */ 1494 if( SCIPgetBestSol(scip) != NULL && heurdata->onlywithoutsol ) 1495 return SCIP_OKAY; 1496 1497 /* stop if there is no LP available */ 1498 if ( ! SCIPhasCurrentNodeLP(scip) ) 1499 return SCIP_OKAY; 1500 1501 if( !SCIPisLPConstructed(scip) ) 1502 { 1503 /* @note this call can have the side effect that variables are created */ 1504 SCIP_CALL( SCIPconstructLP(scip, &cutoff) ); 1505 1506 /* manually cut off the node if the LP construction detected infeasibility (heuristics cannot return such a result) */ 1507 if( cutoff ) 1508 { 1509 SCIP_CALL( SCIPcutoffNode(scip, SCIPgetCurrentNode(scip)) ); 1510 return SCIP_OKAY; 1511 } 1512 1513 SCIP_CALL( SCIPflushLP(scip) ); 1514 } 1515 1516 SCIPstatistic( heurdata->nlpiters = SCIPgetNLPIterations(scip) ); 1517 1518 nlprows = SCIPgetNLPRows(scip); 1519 1520 SCIP_CALL( SCIPgetLPColsData(scip, &lpcols, &nlpcols) ); 1521 assert(nlpcols == 0 || lpcols != NULL); 1522 1523 /* we need an LP */ 1524 if( nlprows == 0 || nlpcols == 0 ) 1525 return SCIP_OKAY; 1526 1527 *result = SCIP_DIDNOTFIND; 1528 initialized = FALSE; 1529 1530 /* allocate lp column array */ 1531 SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->lpcols, nlpcols) ); 1532 heurdata->nlpcols = nlpcols; 1533 1534 impliscontinuous = heurdata->impliscontinuous; 1535 1536 #ifndef NDEBUG 1537 BMSclearMemoryArray(heurdata->lpcols, nlpcols); 1538 #endif 1539 1540 /* copy and sort the columns by their variable types (binary before integer before implicit integer before continuous) */ 1541 BMScopyMemoryArray(heurdata->lpcols, lpcols, nlpcols); 1542 1543 SCIPsortPtr((void**)heurdata->lpcols, heurSortColsShiftandpropagate, nlpcols); 1544 1545 SCIP_CALL( SCIPallocBufferArray(scip, &colposs, nlpcols) ); 1546 1547 /* we have to collect the number of different variable types before we start probing since during probing variable 1548 * can be created (e.g., cons_xor.c) 1549 */ 1550 ndiscvars = 0; 1551 nbinvars = 0; 1552 for( c = 0; c < nlpcols; ++c ) 1553 { 1554 SCIP_COL* col; 1555 SCIP_VAR* colvar; 1556 1557 col = heurdata->lpcols[c]; 1558 assert(col != NULL); 1559 colvar = SCIPcolGetVar(col); 1560 assert(colvar != NULL); 1561 1562 if( varIsDiscrete(colvar, impliscontinuous) ) 1563 ++ndiscvars; 1564 if( SCIPvarGetType(colvar) == SCIP_VARTYPE_BINARY ) 1565 ++nbinvars; 1566 #ifndef NDEBUG 1567 else if( SCIPvarGetType(colvar) == SCIP_VARTYPE_INTEGER ) 1568 ++nintvars; 1569 #endif 1570 1571 /* save the position of this column in the array such that it can be accessed as the "true" column position */ 1572 assert(SCIPcolGetLPPos(col) >= 0); 1573 colposs[SCIPcolGetLPPos(col)] = c; 1574 } 1575 assert(nbinvars + nintvars <= ndiscvars); 1576 1577 /* start probing mode */ 1578 SCIP_CALL( SCIPstartProbing(scip) ); 1579 1580 /* enables collection of variable statistics during probing */ 1581 if( heurdata->collectstats ) 1582 SCIPenableVarHistory(scip); 1583 else 1584 SCIPdisableVarHistory(scip); 1585 1586 /* this should always be fulfilled because we perform shift and propagate only at the root node */ 1587 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(scip)); 1588 1589 /* @todo check if this node is necessary (I don't think so) */ 1590 SCIP_CALL( SCIPnewProbingNode(scip) ); 1591 ncutoffs = 0; 1592 nprobings = 0; 1593 nmaxrows = 0; 1594 infeasible = FALSE; 1595 1596 /* initialize heuristic matrix and working solution */ 1597 SCIP_CALL( SCIPallocBuffer(scip, &matrix) ); 1598 SCIP_CALL( initMatrix(scip, matrix, heurdata, colposs, heurdata->normalize, &nmaxrows, heurdata->relax, &initialized, &infeasible) ); 1599 1600 /* could not initialize matrix */ 1601 if( !initialized || infeasible ) 1602 { 1603 SCIPdebugMsg(scip, " MATRIX not initialized -> Execution of heuristic stopped! \n"); 1604 goto TERMINATE; 1605 } 1606 1607 /* the number of discrete LP column variables can be less than the actual number of variables, if, e.g., there 1608 * are nonlinearities in the problem. The heuristic execution can be terminated in that case. 1609 */ 1610 if( matrix->ndiscvars < ndiscvars ) 1611 { 1612 SCIPdebugMsg(scip, "Not all discrete variables are in the current LP. Shiftandpropagate execution terminated.\n"); 1613 goto TERMINATE; 1614 } 1615 1616 assert(nmaxrows > 0); 1617 1618 eventhdlrdata->matrix = matrix; 1619 eventhdlrdata->heurdata = heurdata; 1620 1621 SCIP_CALL( SCIPcreateSol(scip, &sol, heur) ); 1622 SCIPsolSetHeur(sol, heur); 1623 1624 /* allocate arrays for execution method */ 1625 SCIP_CALL( SCIPallocBufferArray(scip, &permutation, ndiscvars) ); 1626 SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->rowweights, matrix->nrows) ); 1627 1628 /* allocate necessary memory for best shift search */ 1629 SCIP_CALL( SCIPallocBufferArray(scip, &steps, nmaxrows) ); 1630 SCIP_CALL( SCIPallocBufferArray(scip, &violationchange, nmaxrows) ); 1631 1632 /* allocate arrays to store information about infeasible rows */ 1633 SCIP_CALL( SCIPallocBufferArray(scip, &violatedrows, matrix->nrows) ); 1634 SCIP_CALL( SCIPallocBufferArray(scip, &violatedrowpos, matrix->nrows) ); 1635 1636 eventhdlrdata->violatedrows = violatedrows; 1637 eventhdlrdata->violatedrowpos = violatedrowpos; 1638 eventhdlrdata->nviolatedrows = &nviolatedrows; 1639 1640 /* initialize arrays. Before sorting, permutation is the identity permutation */ 1641 for( i = 0; i < ndiscvars; ++i ) 1642 permutation[i] = i; 1643 1644 /* initialize row weights */ 1645 for( r = 0; r < matrix->nrows; ++r ) 1646 { 1647 if( !SCIPisInfinity(scip, -(matrix->lhs[r])) && !SCIPisInfinity(scip, matrix->rhs[r]) ) 1648 heurdata->rowweights[r] = DEFAULT_WEIGHT_EQUALITY; 1649 else 1650 heurdata->rowweights[r] = DEFAULT_WEIGHT_INEQUALITY; 1651 } 1652 colnorms = matrix->colnorms; 1653 1654 assert(nbinvars >= 0); 1655 assert(nintvars >= 0); 1656 1657 /* check rows for infeasibility */ 1658 checkViolations(scip, matrix, -1, violatedrows, violatedrowpos, &nviolatedrows, heurdata->rowweights, heurdata->updateweights); 1659 1660 /* allocate memory for violatedvarrows array only if variable ordering relies on it */ 1661 if( heurdata->sortvars && (heurdata->sortkey == 't' || heurdata->sortkey == 'v') ) 1662 { 1663 SCIP_CALL( SCIPallocBufferArray(scip, &violatedvarrows, ndiscvars) ); 1664 BMScopyMemoryArray(violatedvarrows, matrix->violrows, ndiscvars); 1665 } 1666 else 1667 violatedvarrows = NULL; 1668 1669 /* sort variables w.r.t. the sorting key parameter. Sorting is indirect, all matrix column data 1670 * stays in place, but permutation array gives access to the sorted order of variables 1671 */ 1672 if( heurdata->sortvars ) 1673 { 1674 switch (heurdata->sortkey) 1675 { 1676 case 'n': 1677 /* variable ordering w.r.t. column norms nonincreasing */ 1678 if( heurdata->preferbinaries ) 1679 { 1680 if( nbinvars > 0 ) 1681 SCIPsortDownRealInt(colnorms, permutation, nbinvars); 1682 if( nbinvars < ndiscvars ) 1683 SCIPsortDownRealInt(&colnorms[nbinvars], &permutation[nbinvars], ndiscvars - nbinvars); 1684 } 1685 else 1686 { 1687 SCIPsortDownRealInt(colnorms, permutation, ndiscvars); 1688 } 1689 SCIPdebugMsg(scip, "Variables sorted down w.r.t their normalized columns!\n"); 1690 break; 1691 case 'u': 1692 /* variable ordering w.r.t. column norms nondecreasing */ 1693 if( heurdata->preferbinaries ) 1694 { 1695 if( nbinvars > 0 ) 1696 SCIPsortRealInt(colnorms, permutation, nbinvars); 1697 if( nbinvars < ndiscvars ) 1698 SCIPsortRealInt(&colnorms[nbinvars], &permutation[nbinvars], ndiscvars - nbinvars); 1699 } 1700 else 1701 { 1702 SCIPsortRealInt(colnorms, permutation, ndiscvars); 1703 } 1704 SCIPdebugMsg(scip, "Variables sorted w.r.t their normalized columns!\n"); 1705 break; 1706 case 'v': 1707 /* variable ordering w.r.t. nonincreasing number of violated rows */ 1708 assert(violatedvarrows != NULL); 1709 if( heurdata->preferbinaries ) 1710 { 1711 if( nbinvars > 0 ) 1712 SCIPsortDownIntInt(violatedvarrows, permutation, nbinvars); 1713 if( nbinvars < ndiscvars ) 1714 SCIPsortDownIntInt(&violatedvarrows[nbinvars], &permutation[nbinvars], ndiscvars - nbinvars); 1715 } 1716 else 1717 { 1718 SCIPsortDownIntInt(violatedvarrows, permutation, ndiscvars); 1719 } 1720 1721 SCIPdebugMsg(scip, "Variables sorted down w.r.t their number of currently infeasible rows!\n"); 1722 break; 1723 case 't': 1724 /* variable ordering w.r.t. nondecreasing number of violated rows */ 1725 assert(violatedvarrows != NULL); 1726 if( heurdata->preferbinaries ) 1727 { 1728 if( nbinvars > 0 ) 1729 SCIPsortIntInt(violatedvarrows, permutation, nbinvars); 1730 if( nbinvars < ndiscvars ) 1731 SCIPsortIntInt(&violatedvarrows[nbinvars], &permutation[nbinvars], ndiscvars - nbinvars); 1732 } 1733 else 1734 { 1735 SCIPsortIntInt(violatedvarrows, permutation, ndiscvars); 1736 } 1737 1738 SCIPdebugMsg(scip, "Variables sorted (upwards) w.r.t their number of currently infeasible rows!\n"); 1739 break; 1740 case 'r': 1741 /* random sorting */ 1742 if( heurdata->preferbinaries ) 1743 { 1744 if( nbinvars > 0 ) 1745 SCIPrandomPermuteIntArray(heurdata->randnumgen, permutation, 0, nbinvars - 1); 1746 if( nbinvars < ndiscvars ) 1747 SCIPrandomPermuteIntArray(heurdata->randnumgen, &permutation[nbinvars], nbinvars - 1, 1748 ndiscvars - nbinvars - 1); 1749 } 1750 else 1751 { 1752 SCIPrandomPermuteIntArray(heurdata->randnumgen, permutation, 0, ndiscvars - 1); 1753 } 1754 SCIPdebugMsg(scip, "Variables permuted randomly!\n"); 1755 break; 1756 default: 1757 SCIPdebugMsg(scip, "No variable permutation applied\n"); 1758 break; 1759 } 1760 } 1761 1762 /* should binary variables without locks be treated first? */ 1763 if( heurdata->binlocksfirst ) 1764 { 1765 SCIP_VAR* var; 1766 int nbinwithoutlocks = 0; 1767 1768 /* count number of binaries without locks */ 1769 if( heurdata->preferbinaries ) 1770 { 1771 for( c = 0; c < nbinvars; ++c ) 1772 { 1773 var = SCIPcolGetVar(heurdata->lpcols[permutation[c]]); 1774 if( SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 0 1775 || SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 0 ) 1776 ++nbinwithoutlocks; 1777 } 1778 } 1779 else 1780 { 1781 for( c = 0; c < ndiscvars; ++c ) 1782 { 1783 var = SCIPcolGetVar(heurdata->lpcols[permutation[c]]); 1784 if( SCIPvarIsBinary(var) ) 1785 { 1786 if( SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 0 1787 || SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 0 ) 1788 ++nbinwithoutlocks; 1789 } 1790 } 1791 } 1792 1793 if( nbinwithoutlocks > 0 ) 1794 { 1795 SCIP_VAR* binvar; 1796 int b = 1; 1797 int tmp; 1798 c = 0; 1799 1800 /* if c reaches nbinwithoutlocks, then all binary variables without locks were sorted to the beginning of the array */ 1801 while( c < nbinwithoutlocks && b < ndiscvars ) 1802 { 1803 assert(c < b); 1804 assert(c < ndiscvars); 1805 assert(b < ndiscvars); 1806 var = SCIPcolGetVar(heurdata->lpcols[permutation[c]]); 1807 binvar = SCIPcolGetVar(heurdata->lpcols[permutation[b]]); 1808 1809 /* search for next variable which is not a binary variable without locks */ 1810 while( SCIPvarIsBinary(var) && (SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 0 1811 || SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 0) ) 1812 { 1813 ++c; 1814 if( c >= nbinwithoutlocks ) 1815 break; 1816 var = SCIPcolGetVar(heurdata->lpcols[permutation[c]]); 1817 } 1818 if( c >= nbinwithoutlocks ) 1819 break; 1820 1821 /* search for next binary variable without locks (with position > c) */ 1822 if( b <= c ) 1823 { 1824 b = c + 1; 1825 binvar = SCIPcolGetVar(heurdata->lpcols[permutation[b]]); 1826 } 1827 while( !SCIPvarIsBinary(binvar) || (SCIPvarGetNLocksUpType(binvar, SCIP_LOCKTYPE_MODEL) > 0 1828 && SCIPvarGetNLocksDownType(binvar, SCIP_LOCKTYPE_MODEL) > 0) ) 1829 { 1830 ++b; 1831 assert(b < ndiscvars); 1832 binvar = SCIPcolGetVar(heurdata->lpcols[permutation[b]]); 1833 } 1834 1835 /* swap the two variables */ 1836 tmp = permutation[b]; 1837 permutation[b] = permutation[c]; 1838 permutation[c] = tmp; 1839 1840 /* increase counters */ 1841 ++c; 1842 ++b; 1843 } 1844 } 1845 1846 #ifndef NDEBUG 1847 for( c = 0; c < ndiscvars; ++c ) 1848 { 1849 assert((c < nbinwithoutlocks) == (SCIPvarIsBinary(SCIPcolGetVar(heurdata->lpcols[permutation[c]])) 1850 && (SCIPvarGetNLocksUpType(SCIPcolGetVar(heurdata->lpcols[permutation[c]]), SCIP_LOCKTYPE_MODEL) == 0 1851 || SCIPvarGetNLocksDownType(SCIPcolGetVar(heurdata->lpcols[permutation[c]]), SCIP_LOCKTYPE_MODEL) == 0))); 1852 } 1853 #endif 1854 } 1855 1856 SCIP_CALL( SCIPallocBufferArray(scip, &eventdatas, matrix->ndiscvars) ); 1857 BMSclearMemoryArray(eventdatas, matrix->ndiscvars); 1858 1859 /* initialize variable events to catch bound changes during propagation */ 1860 for( c = 0; c < matrix->ndiscvars; ++c ) 1861 { 1862 SCIP_VAR* var; 1863 1864 var = SCIPcolGetVar(heurdata->lpcols[c]); 1865 assert(var != NULL); 1866 assert(SCIPvarIsIntegral(var)); 1867 assert(eventdatas[c] == NULL); 1868 1869 SCIP_CALL( SCIPallocBuffer(scip, &(eventdatas[c])) ); /*lint !e866*/ 1870 1871 eventdatas[c]->colpos = c; 1872 1873 SCIP_CALL( SCIPcatchVarEvent(scip, var, EVENTTYPE_SHIFTANDPROPAGATE, eventhdlr, eventdatas[c], NULL) ); 1874 } 1875 1876 cutoff = FALSE; 1877 1878 lastindexofsusp = -1; 1879 probing = heurdata->probing; 1880 infeasible = FALSE; 1881 1882 SCIPdebugMsg(scip, "SHIFT_AND_PROPAGATE heuristic starts main loop with %d violations and %d remaining variables!\n", 1883 nviolatedrows, ndiscvars); 1884 1885 assert(matrix->ndiscvars == ndiscvars); 1886 1887 /* loop over variables, shift them according to shifting criteria and try to reduce the global infeasibility */ 1888 for( c = 0; c < ndiscvars; ++c ) 1889 { 1890 SCIP_VAR* var; 1891 SCIP_Longint ndomredsfound; 1892 SCIP_Real optimalshiftvalue; 1893 SCIP_Real origsolval; 1894 SCIP_Real lb; 1895 SCIP_Real ub; 1896 int nviolations; 1897 int permutedvarindex; 1898 int j; 1899 SCIP_Bool marksuspicious; 1900 1901 if( heurdata->selectbest ) 1902 { /* search for best candidate */ 1903 j = c + 1; 1904 while( j < ndiscvars ) 1905 { 1906 /* run through remaining variables and search for best candidate */ 1907 if( matrix->violrows[permutation[c]] < matrix->violrows[permutation[j]] ) 1908 { 1909 int tmp; 1910 tmp = permutation[c]; 1911 permutation[c] = permutation[j]; 1912 permutation[j] = tmp; 1913 } 1914 ++j; 1915 } 1916 } 1917 permutedvarindex = permutation[c]; 1918 optimalshiftvalue = 0.0; 1919 nviolations = 0; 1920 var = SCIPcolGetVar(heurdata->lpcols[permutedvarindex]); 1921 lb = SCIPvarGetLbLocal(var); 1922 ub = SCIPvarGetUbLocal(var); 1923 assert(SCIPcolGetLPPos(SCIPvarGetCol(var)) >= 0); 1924 assert(SCIPvarIsIntegral(var)); 1925 1926 /* check whether we hit some limit, e.g. the time limit, in between 1927 * since the check itself consumes some time, we only do it every tenth iteration 1928 */ 1929 if( c % 10 == 0 && SCIPisStopped(scip) ) 1930 goto TERMINATE2; 1931 1932 /* if propagation is enabled, check if propagation has changed the variables bounds 1933 * and update the transformed upper bound correspondingly 1934 * @todo this should not be necessary 1935 */ 1936 if( heurdata->probing ) 1937 { 1938 SCIP_CALL( updateTransformation(scip, matrix, heurdata, permutedvarindex,lb, ub, violatedrows, violatedrowpos, 1939 &nviolatedrows) ); 1940 } 1941 1942 SCIPdebugMsg(scip, "Variable %s with local bounds [%g,%g], status <%d>, matrix bound <%g>\n", 1943 SCIPvarGetName(var), lb, ub, matrix->transformstatus[permutedvarindex], matrix->upperbounds[permutedvarindex]); 1944 1945 /* ignore variable if propagation fixed it (lb and ub will be zero) */ 1946 if( SCIPisFeasZero(scip, matrix->upperbounds[permutedvarindex]) ) 1947 { 1948 assert(!SCIPisInfinity(scip, ub)); 1949 assert(SCIPisFeasEQ(scip, lb, ub)); 1950 1951 SCIP_CALL( SCIPsetSolVal(scip, sol, var, ub) ); 1952 1953 continue; 1954 } 1955 1956 marksuspicious = FALSE; 1957 1958 /* check whether the variable is binary and has no locks in one direction, so that we want to fix it to the 1959 * respective bound (only enabled by parameter) 1960 */ 1961 if( heurdata->fixbinlocks && SCIPvarIsBinary(var) 1962 && (SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 0 1963 || SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 0) ) 1964 { 1965 if( SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) == 0 ) 1966 origsolval = SCIPvarGetUbLocal(var); 1967 else 1968 { 1969 assert(SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) == 0); 1970 origsolval = SCIPvarGetLbLocal(var); 1971 } 1972 } 1973 else 1974 { 1975 /* only apply the computationally expensive best shift selection, if there is a violated row left */ 1976 if( !heurdata->stopafterfeasible || nviolatedrows > 0 ) 1977 { 1978 /* compute optimal shift value for variable */ 1979 SCIP_CALL( getOptimalShiftingValue(scip, matrix, permutedvarindex, 1, heurdata->rowweights, steps, violationchange, 1980 &optimalshiftvalue, &nviolations) ); 1981 assert(SCIPisFeasGE(scip, optimalshiftvalue, 0.0)); 1982 1983 /* Variables with FREE transform have to be dealt with twice */ 1984 if( matrix->transformstatus[permutedvarindex] == TRANSFORMSTATUS_FREE ) 1985 { 1986 SCIP_Real downshiftvalue; 1987 int ndownviolations; 1988 1989 downshiftvalue = 0.0; 1990 ndownviolations = 0; 1991 SCIP_CALL( getOptimalShiftingValue(scip, matrix, permutedvarindex, -1, heurdata->rowweights, steps, violationchange, 1992 &downshiftvalue, &ndownviolations) ); 1993 1994 assert(SCIPisLE(scip, downshiftvalue, 0.0)); 1995 1996 /* compare to positive direction and select the direction which makes more rows feasible */ 1997 if( ndownviolations < nviolations ) 1998 { 1999 optimalshiftvalue = downshiftvalue; 2000 } 2001 } 2002 } 2003 else 2004 optimalshiftvalue = 0.0; 2005 2006 /* if zero optimal shift values are forbidden by the user parameter, delay the variable by marking it suspicious */ 2007 if( heurdata->nozerofixing && nviolations > 0 && SCIPisFeasZero(scip, optimalshiftvalue) ) 2008 marksuspicious = TRUE; 2009 2010 /* retransform the solution value from the heuristic transformation space */ 2011 assert(varIsDiscrete(var, impliscontinuous)); 2012 origsolval = retransformVariable(scip, matrix, var, permutedvarindex, optimalshiftvalue); 2013 } 2014 assert(SCIPisFeasGE(scip, origsolval, lb) && SCIPisFeasLE(scip, origsolval, ub)); 2015 2016 /* check if propagation should still be performed 2017 * @todo do we need the hard coded value? we could use SCIP_MAXTREEDEPTH 2018 */ 2019 if( nprobings > DEFAULT_PROPBREAKER ) 2020 probing = FALSE; 2021 2022 /* if propagation is enabled, fix the variable to the new solution value and propagate the fixation 2023 * (to fix other variables and to find out early whether solution is already infeasible) 2024 */ 2025 if( !marksuspicious && probing ) 2026 { 2027 /* this assert should be always fulfilled because we run this heuristic at the root node only and do not 2028 * perform probing if nprobings is less than DEFAULT_PROPBREAKER (currently: 65000) 2029 */ 2030 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(scip)); 2031 2032 SCIP_CALL( SCIPnewProbingNode(scip) ); 2033 SCIP_CALL( SCIPfixVarProbing(scip, var, origsolval) ); 2034 ndomredsfound = 0; 2035 2036 SCIPdebugMsg(scip, " Shift %g(%g originally) is optimal, propagate solution\n", optimalshiftvalue, origsolval); 2037 SCIP_CALL( SCIPpropagateProbing(scip, heurdata->nproprounds, &cutoff, &ndomredsfound) ); 2038 2039 ++nprobings; 2040 SCIPstatistic( heurdata->ntotaldomredsfound += ndomredsfound ); 2041 SCIPdebugMsg(scip, "Propagation finished! <%" SCIP_LONGINT_FORMAT "> domain reductions %s, <%d> probing depth\n", ndomredsfound, cutoff ? "CUTOFF" : "", 2042 SCIPgetProbingDepth(scip)); 2043 } 2044 assert(!cutoff || probing); 2045 2046 /* propagation led to an empty domain, hence we backtrack and postpone the variable */ 2047 if( cutoff ) 2048 { 2049 assert(probing); 2050 2051 ++ncutoffs; 2052 2053 /* only continue heuristic if number of cutoffs occured so far is reasonably small */ 2054 if( heurdata->cutoffbreaker >= 0 && ncutoffs >= ((heurdata->maxcutoffquot * SCIPgetProbingDepth(scip)) + heurdata->cutoffbreaker) ) 2055 break; 2056 2057 cutoff = FALSE; 2058 2059 /* backtrack to the parent of the current node */ 2060 assert(SCIPgetProbingDepth(scip) >= 1); 2061 SCIP_CALL( SCIPbacktrackProbing(scip, SCIPgetProbingDepth(scip) - 1) ); 2062 2063 /* this assert should be always fulfilled because we run this heuristic at the root node only and do not 2064 * perform probing if nprobings is less than DEFAULT_PROPBREAKER (currently: 65000) 2065 */ 2066 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(scip)); 2067 2068 /* if the variable upper and lower bound are equal to the solution value to which we tried to fix the variable, 2069 * we are trapped at an infeasible node and break; this can only happen due to an intermediate global bound change of the variable, 2070 * I guess 2071 */ 2072 if( SCIPisFeasEQ(scip, SCIPvarGetUbLocal(var), origsolval) && SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), origsolval) ) 2073 { 2074 cutoff = TRUE; 2075 break; 2076 } 2077 else if( SCIPisFeasEQ(scip, SCIPvarGetLbLocal(var), origsolval) && REALABS( origsolval ) < 1.0 / SCIPepsilon(scip) ) 2078 { 2079 /* if the variable was set to one of its bounds, repropagate by tightening this bound by 1.0 into the 2080 * direction of the other bound, if possible; if the bound is too large (in abs value) do not even bother 2081 */ 2082 assert(SCIPisFeasGE(scip, SCIPvarGetUbLocal(var), origsolval + 1.0)); 2083 2084 ndomredsfound = 0; 2085 SCIP_CALL( SCIPnewProbingNode(scip) ); 2086 SCIP_CALL( SCIPchgVarLbProbing(scip, var, origsolval + 1.0) ); 2087 SCIP_CALL( SCIPpropagateProbing(scip, heurdata->nproprounds, &cutoff, &ndomredsfound) ); 2088 2089 SCIPstatistic( heurdata->ntotaldomredsfound += ndomredsfound ); 2090 } 2091 else if( SCIPisFeasEQ(scip, SCIPvarGetUbLocal(var), origsolval) && REALABS( origsolval ) < 1.0 / SCIPepsilon(scip) ) 2092 { 2093 /* if the variable was set to one of its bounds, repropagate by tightening this bound by 1.0 into the 2094 * direction of the other bound, if possible; if the bound is too large (in abs value) do not even bother 2095 */ 2096 assert(SCIPisFeasLE(scip, SCIPvarGetLbLocal(var), origsolval - 1.0)); 2097 2098 ndomredsfound = 0; 2099 2100 SCIP_CALL( SCIPnewProbingNode(scip) ); 2101 SCIP_CALL( SCIPchgVarUbProbing(scip, var, origsolval - 1.0) ); 2102 SCIP_CALL( SCIPpropagateProbing(scip, heurdata->nproprounds, &cutoff, &ndomredsfound) ); 2103 2104 SCIPstatistic( heurdata->ntotaldomredsfound += ndomredsfound ); 2105 } 2106 2107 /* if the tightened bound again leads to a cutoff, both subproblems are proven infeasible and the heuristic 2108 * can be stopped */ 2109 if( cutoff ) 2110 { 2111 break; 2112 } 2113 else 2114 { 2115 /* since repropagation was successful, we indicate that this variable led to a cutoff in one direction */ 2116 marksuspicious = TRUE; 2117 } 2118 } 2119 2120 if( marksuspicious ) 2121 { 2122 /* mark the variable as suspicious */ 2123 assert(permutedvarindex == permutation[c]); 2124 2125 ++lastindexofsusp; 2126 assert(lastindexofsusp >= 0 && lastindexofsusp <= c); 2127 2128 permutation[c] = permutation[lastindexofsusp]; 2129 permutation[lastindexofsusp] = permutedvarindex; 2130 2131 SCIPdebugMsg(scip, " Suspicious variable! Postponed from pos <%d> to position <%d>\n", c, lastindexofsusp); 2132 } 2133 else 2134 { 2135 SCIPdebugMsg(scip, "Variable <%d><%s> successfully shifted by value <%g>!\n", permutedvarindex, 2136 SCIPvarGetName(var), optimalshiftvalue); 2137 2138 /* update solution */ 2139 SCIP_CALL( SCIPsetSolVal(scip, sol, var, origsolval) ); 2140 2141 /* only to ensure that some assertions can be made later on */ 2142 if( !probing ) 2143 { 2144 SCIP_CALL( SCIPfixVarProbing(scip, var, origsolval) ); 2145 } 2146 } 2147 } 2148 SCIPdebugMsg(scip, "Heuristic finished with %d remaining violations and %d remaining variables!\n", 2149 nviolatedrows, lastindexofsusp + 1); 2150 2151 /* if constructed solution might be feasible, go through the queue of suspicious variables and set the solution 2152 * values 2153 */ 2154 if( nviolatedrows == 0 && !cutoff ) 2155 { 2156 SCIP_Bool stored; 2157 SCIP_Bool trysol; 2158 SCIP_Bool solvelp; 2159 2160 for( v = 0; v <= lastindexofsusp; ++v ) 2161 { 2162 SCIP_VAR* var; 2163 SCIP_Real origsolval; 2164 int permutedvarindex; 2165 2166 /* get the column position of the variable */ 2167 permutedvarindex = permutation[v]; 2168 var = SCIPcolGetVar(heurdata->lpcols[permutedvarindex]); 2169 assert(varIsDiscrete(var, impliscontinuous)); 2170 2171 /* update the transformation of the variable, since the bound might have changed after the last update. */ 2172 if( heurdata->probing ) 2173 SCIP_CALL( updateTransformation(scip, matrix, heurdata, permutedvarindex, SCIPvarGetLbLocal(var), 2174 SCIPvarGetUbLocal(var), violatedrows, violatedrowpos, &nviolatedrows) ); 2175 2176 /* retransform the solution value from the heuristic transformed space, set the solution value accordingly */ 2177 assert(varIsDiscrete(var, impliscontinuous)); 2178 origsolval = retransformVariable(scip, matrix, var, permutedvarindex, 0.0); 2179 assert(SCIPisFeasGE(scip, origsolval, SCIPvarGetLbLocal(var)) 2180 && SCIPisFeasLE(scip, origsolval, SCIPvarGetUbLocal(var))); 2181 SCIP_CALL( SCIPsetSolVal(scip, sol, var, origsolval) ); 2182 SCIP_CALL( SCIPfixVarProbing(scip, var, origsolval) ); /* only to ensure that some assertions can be made later */ 2183 2184 SCIPdebugMsg(scip, " Remaining variable <%s> set to <%g>; %d Violations\n", SCIPvarGetName(var), origsolval, 2185 nviolatedrows); 2186 } 2187 2188 /* Fixing of remaining variables led to infeasibility */ 2189 if( nviolatedrows > 0 ) 2190 goto TERMINATE2; 2191 2192 trysol = TRUE; 2193 2194 /* check if enough variables have been fixed (including continuous) to solve the remaining LP */ 2195 if( nlpcols != matrix->ndiscvars ) 2196 { 2197 SCIP_VAR** vars; 2198 int nvars = SCIPgetNVars(scip); 2199 int nminfixings = (int)(SCIPceil(scip, heurdata->minfixingratelp * nvars)); 2200 int nfixedvars = ndiscvars; 2201 2202 vars = SCIPgetVars(scip); 2203 2204 /* count fixed variables */ 2205 for( v = ndiscvars; v < nvars && nfixedvars < nminfixings; ++v ) 2206 { 2207 if( SCIPisEQ(scip, SCIPvarGetLbLocal(vars[v]), SCIPvarGetUbLocal(vars[v])) ) 2208 ++nfixedvars; 2209 } 2210 2211 solvelp = (nfixedvars >= nminfixings); 2212 trysol = solvelp; 2213 SCIPdebugMsg(scip, "Fixed %d of %d (%.1f %%) variables after probing -> %s\n", 2214 nfixedvars, nvars, (100.0 * nfixedvars / (SCIP_Real)nvars), 2215 solvelp ? "continue and solve LP for remaining variables" : "terminate without LP"); 2216 } 2217 else /* no need to solve an LP */ 2218 solvelp = FALSE; 2219 2220 /* if the constructed solution might still be extendable to a feasible solution, try this by 2221 * solving the remaining LP 2222 */ 2223 if( solvelp ) 2224 { 2225 char strbuf[SCIP_MAXSTRLEN]; 2226 /* case that remaining LP has to be solved */ 2227 SCIP_Bool lperror; 2228 int ncols; 2229 2230 #ifndef NDEBUG 2231 { 2232 SCIP_VAR** vars; 2233 2234 vars = SCIPgetVars(scip); 2235 assert(vars != NULL); 2236 /* ensure that all discrete variables in the remaining LP are fixed */ 2237 for( v = 0; v < ndiscvars; ++v ) 2238 { 2239 if( SCIPvarIsInLP(vars[v]) ) 2240 { 2241 assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(vars[v]), SCIPvarGetUbLocal(vars[v]))); 2242 } 2243 } 2244 } 2245 #endif 2246 2247 /* print message if relatively large LP is solved from scratch, since this could lead to a longer period during 2248 * which the user sees no output; more detailed probing stats only in debug mode */ 2249 ncols = SCIPgetNLPCols(scip); 2250 if( !SCIPisLPSolBasic(scip) && ncols > 1000 ) 2251 { 2252 int nunfixedcols = SCIPgetNUnfixedLPCols(scip); 2253 2254 if( nunfixedcols > 0.5 * ncols ) 2255 { 2256 SCIPverbMessage(scip, SCIP_VERBLEVEL_FULL, NULL, 2257 "Heuristic " HEUR_NAME " solving LP from scratch with %.1f %% unfixed columns (%d of %d) ...\n", 2258 100.0 * (nunfixedcols / (SCIP_Real)ncols), nunfixedcols, ncols); 2259 } 2260 } 2261 SCIPdebugMsg(scip, "Heuristic " HEUR_NAME " probing LP: %s\n", 2262 SCIPsnprintfProbingStats(scip, strbuf, SCIP_MAXSTRLEN)); 2263 SCIPdebugMsg(scip, " -> old LP iterations: %" SCIP_LONGINT_FORMAT "\n", SCIPgetNLPIterations(scip)); 2264 2265 #ifdef SCIP_DEBUG 2266 SCIP_CALL( SCIPwriteLP(scip, "shiftandpropagatelp.mps") ); 2267 #endif 2268 /* solve LP; 2269 * errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic. 2270 * hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. 2271 */ 2272 #ifdef NDEBUG 2273 { 2274 SCIP_RETCODE retstat; 2275 retstat = SCIPsolveProbingLP(scip, -1, &lperror, NULL); 2276 if( retstat != SCIP_OKAY ) 2277 { 2278 SCIPwarningMessage(scip, "Error while solving LP in SHIFTANDPROPAGATE heuristic; LP solve terminated with code <%d>\n", 2279 retstat); 2280 } 2281 } 2282 #else 2283 SCIP_CALL( SCIPsolveProbingLP(scip, -1, &lperror, NULL) ); 2284 #endif 2285 2286 SCIPdebugMsg(scip, " -> new LP iterations: %" SCIP_LONGINT_FORMAT "\n", SCIPgetNLPIterations(scip)); 2287 SCIPdebugMsg(scip, " -> error=%u, status=%d\n", lperror, SCIPgetLPSolstat(scip)); 2288 2289 /* check if this is a feasible solution */ 2290 if( !lperror && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ) 2291 { 2292 /* copy the current LP solution to the working solution */ 2293 SCIP_CALL( SCIPlinkLPSol(scip, sol) ); 2294 } 2295 else 2296 trysol = FALSE; 2297 2298 SCIPstatistic( heurdata->lpsolstat = SCIPgetLPSolstat(scip) ); 2299 } 2300 2301 /* check solution for feasibility, and add it to solution store if possible. 2302 * None of integrality, feasibility of LP rows, variable bounds have to be checked, because they 2303 * are guaranteed by the heuristic at this stage. 2304 */ 2305 if( trysol ) 2306 { 2307 SCIP_Bool printreason; 2308 SCIP_Bool completely; 2309 #ifdef SCIP_DEBUG 2310 printreason = TRUE; 2311 #else 2312 printreason = FALSE; 2313 #endif 2314 #ifndef NDEBUG 2315 completely = TRUE; /*lint !e838*/ 2316 #else 2317 completely = FALSE; 2318 #endif 2319 2320 /* we once also checked the variable bounds which should not be necessary */ 2321 SCIP_CALL( SCIPtrySol(scip, sol, printreason, completely, FALSE, FALSE, FALSE, &stored) ); 2322 2323 if( stored ) 2324 { 2325 SCIPdebugMsg(scip, "found feasible shifted solution:\n"); 2326 SCIPdebug( SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) ) ); 2327 *result = SCIP_FOUNDSOL; 2328 2329 SCIPstatisticMessage(" Shiftandpropagate solution value: %16.9g \n", SCIPgetSolOrigObj(scip, sol)); 2330 } 2331 } 2332 } 2333 else 2334 { 2335 SCIPdebugMsg(scip, "Solution constructed by heuristic is already known to be infeasible\n"); 2336 } 2337 2338 SCIPstatistic( heurdata->nremainingviols = nviolatedrows; ); 2339 2340 TERMINATE2: 2341 /* free allocated memory in reverse order of allocation */ 2342 for( c = matrix->ndiscvars - 1; c >= 0; --c ) 2343 { 2344 SCIP_VAR* var; 2345 2346 var = SCIPcolGetVar(heurdata->lpcols[c]); 2347 assert(var != NULL); 2348 assert(eventdatas[c] != NULL); 2349 2350 SCIP_CALL( SCIPdropVarEvent(scip, var, EVENTTYPE_SHIFTANDPROPAGATE, eventhdlr, eventdatas[c], -1) ); 2351 SCIPfreeBuffer(scip, &(eventdatas[c])); 2352 } 2353 SCIPfreeBufferArray(scip, &eventdatas); 2354 2355 if( violatedvarrows != NULL ) 2356 { 2357 assert(heurdata->sortkey == 'v' || heurdata->sortkey == 't'); 2358 SCIPfreeBufferArray(scip, &violatedvarrows); 2359 } 2360 /* free all allocated memory */ 2361 SCIPfreeBufferArray(scip, &violatedrowpos); 2362 SCIPfreeBufferArray(scip, &violatedrows); 2363 SCIPfreeBufferArray(scip, &violationchange); 2364 SCIPfreeBufferArray(scip, &steps); 2365 SCIPfreeBufferArray(scip, &heurdata->rowweights); 2366 SCIPfreeBufferArray(scip, &permutation); 2367 SCIP_CALL( SCIPfreeSol(scip, &sol) ); 2368 2369 eventhdlrdata->nviolatedrows = NULL; 2370 eventhdlrdata->violatedrowpos = NULL; 2371 eventhdlrdata->violatedrows = NULL; 2372 2373 TERMINATE: 2374 /* terminate probing mode and free the remaining memory */ 2375 SCIPstatistic( 2376 heurdata->ncutoffs += ncutoffs; 2377 heurdata->nprobings += nprobings; 2378 heurdata->nlpiters = SCIPgetNLPIterations(scip) - heurdata->nlpiters; 2379 ); 2380 2381 SCIP_CALL( SCIPendProbing(scip) ); 2382 freeMatrix(scip, &matrix); 2383 SCIPfreeBufferArray(scip, &colposs); 2384 SCIPfreeBufferArray(scip, &heurdata->lpcols); 2385 eventhdlrdata->matrix = NULL; 2386 2387 return SCIP_OKAY; 2388 } 2389 2390 /** event handler execution method for the heuristic which catches all 2391 * events in which a lower or upper bound were tightened */ 2392 static 2393 SCIP_DECL_EVENTEXEC(eventExecShiftandpropagate) 2394 { /*lint --e{715}*/ 2395 SCIP_EVENTHDLRDATA* eventhdlrdata; 2396 SCIP_VAR* var; 2397 SCIP_COL* col; 2398 SCIP_Real lb; 2399 SCIP_Real ub; 2400 int colpos; 2401 CONSTRAINTMATRIX* matrix; 2402 SCIP_HEURDATA* heurdata; 2403 2404 assert(scip != NULL); 2405 assert(eventhdlr != NULL); 2406 assert(strcmp(EVENTHDLR_NAME, SCIPeventhdlrGetName(eventhdlr)) == 0); 2407 2408 eventhdlrdata = SCIPeventhdlrGetData(eventhdlr); 2409 assert(eventhdlrdata != NULL); 2410 2411 matrix = eventhdlrdata->matrix; 2412 2413 heurdata = eventhdlrdata->heurdata; 2414 assert(heurdata != NULL && heurdata->lpcols != NULL); 2415 2416 colpos = eventdata->colpos; 2417 2418 assert(0 <= colpos && colpos < matrix->ndiscvars); 2419 2420 col = heurdata->lpcols[colpos]; 2421 var = SCIPcolGetVar(col); 2422 2423 lb = SCIPvarGetLbLocal(var); 2424 ub = SCIPvarGetUbLocal(var); 2425 2426 SCIP_CALL( updateTransformation(scip, matrix, eventhdlrdata->heurdata, colpos, lb, ub, eventhdlrdata->violatedrows, 2427 eventhdlrdata->violatedrowpos, eventhdlrdata->nviolatedrows) ); 2428 2429 return SCIP_OKAY; 2430 } 2431 2432 /* 2433 * primal heuristic specific interface methods 2434 */ 2435 2436 /** creates the shiftandpropagate primal heuristic and includes it in SCIP */ 2437 SCIP_RETCODE SCIPincludeHeurShiftandpropagate( 2438 SCIP* scip /**< SCIP data structure */ 2439 ) 2440 { 2441 SCIP_HEURDATA* heurdata; 2442 SCIP_HEUR* heur; 2443 SCIP_EVENTHDLRDATA* eventhandlerdata; 2444 SCIP_EVENTHDLR* eventhdlr; 2445 2446 SCIP_CALL( SCIPallocBlockMemory(scip, &eventhandlerdata) ); 2447 eventhandlerdata->matrix = NULL; 2448 2449 eventhdlr = NULL; 2450 SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC, 2451 eventExecShiftandpropagate, eventhandlerdata) ); 2452 assert(eventhdlr != NULL); 2453 2454 /* create Shiftandpropagate primal heuristic data */ 2455 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 2456 heurdata->rowweights = NULL; 2457 heurdata->nlpcols = 0; 2458 heurdata->eventhdlr = eventhdlr; 2459 2460 /* include primal heuristic */ 2461 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 2462 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 2463 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecShiftandpropagate, heurdata) ); 2464 2465 assert(heur != NULL); 2466 2467 /* set non-NULL pointers to callback methods */ 2468 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyShiftandpropagate) ); 2469 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeShiftandpropagate) ); 2470 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitShiftandpropagate) ); 2471 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitShiftandpropagate) ); 2472 2473 /* add shiftandpropagate primal heuristic parameters */ 2474 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/nproprounds", 2475 "The number of propagation rounds used for each propagation", 2476 &heurdata->nproprounds, TRUE, DEFAULT_NPROPROUNDS, -1, 1000, NULL, NULL) ); 2477 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/relax", "Should continuous variables be relaxed?", 2478 &heurdata->relax, TRUE, DEFAULT_RELAX, NULL, NULL) ); 2479 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/probing", "Should domains be reduced by probing?", 2480 &heurdata->probing, TRUE, DEFAULT_PROBING, NULL, NULL) ); 2481 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/onlywithoutsol", 2482 "Should heuristic only be executed if no primal solution was found, yet?", 2483 &heurdata->onlywithoutsol, TRUE, DEFAULT_ONLYWITHOUTSOL, NULL, NULL) ); 2484 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/cutoffbreaker", "The number of cutoffs before heuristic stops", 2485 &heurdata->cutoffbreaker, TRUE, DEFAULT_CUTOFFBREAKER, -1, 1000000, NULL, NULL) ); 2486 SCIP_CALL( SCIPaddCharParam(scip, "heuristics/" HEUR_NAME "/sortkey", 2487 "the key for variable sorting: (n)orms down, norms (u)p, (v)iolations down, viola(t)ions up, or (r)andom", 2488 &heurdata->sortkey, TRUE, DEFAULT_SORTKEY, SORTKEYS, NULL, NULL) ); 2489 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/sortvars", "Should variables be sorted for the heuristic?", 2490 &heurdata->sortvars, TRUE, DEFAULT_SORTVARS, NULL, NULL)); 2491 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/collectstats", "should variable statistics be collected during probing?", 2492 &heurdata->collectstats, TRUE, DEFAULT_COLLECTSTATS, NULL, NULL) ); 2493 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/stopafterfeasible", 2494 "Should the heuristic stop calculating optimal shift values when no more rows are violated?", 2495 &heurdata->stopafterfeasible, TRUE, DEFAULT_STOPAFTERFEASIBLE, NULL, NULL) ); 2496 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/preferbinaries", 2497 "Should binary variables be shifted first?", 2498 &heurdata->preferbinaries, TRUE, DEFAULT_PREFERBINARIES, NULL, NULL) ); 2499 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/nozerofixing", 2500 "should variables with a zero shifting value be delayed instead of being fixed?", 2501 &heurdata->nozerofixing, TRUE, DEFAULT_NOZEROFIXING, NULL, NULL) ); 2502 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/fixbinlocks", 2503 "should binary variables with no locks in one direction be fixed to that direction?", 2504 &heurdata->fixbinlocks, TRUE, DEFAULT_FIXBINLOCKS, NULL, NULL) ); 2505 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/binlocksfirst", 2506 "should binary variables with no locks be preferred in the ordering?", 2507 &heurdata->binlocksfirst, TRUE, DEFAULT_BINLOCKSFIRST, NULL, NULL) ); 2508 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/normalize", 2509 "should coefficients and left/right hand sides be normalized by max row coeff?", 2510 &heurdata->normalize, TRUE, DEFAULT_NORMALIZE, NULL, NULL) ); 2511 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/updateweights", 2512 "should row weight be increased every time the row is violated?", 2513 &heurdata->updateweights, TRUE, DEFAULT_UPDATEWEIGHTS, NULL, NULL) ); 2514 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/impliscontinuous", 2515 "should implicit integer variables be treated as continuous variables?", 2516 &heurdata->impliscontinuous, TRUE, DEFAULT_IMPLISCONTINUOUS, NULL, NULL) ); 2517 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/shiftandpropagate/selectbest", 2518 "should the heuristic choose the best candidate in every round? (set to FALSE for static order)?", 2519 &heurdata->selectbest, TRUE, DEFAULT_SELECTBEST, NULL, NULL) ); 2520 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/maxcutoffquot", 2521 "maximum percentage of allowed cutoffs before stopping the heuristic", 2522 &heurdata->maxcutoffquot, TRUE, DEFAULT_MAXCUTOFFQUOT, 0.0, 2.0, NULL, NULL) ); 2523 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minfixingratelp", 2524 "minimum fixing rate over all variables (including continuous) to solve LP", 2525 &heurdata->minfixingratelp, TRUE, DEFAULT_MINFIXINGRATELP, 0.0, 1.0, NULL, NULL) ); 2526 2527 return SCIP_OKAY; 2528 } 2529