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_twoopt.c 26 * @ingroup DEFPLUGINS_HEUR 27 * @brief primal heuristic to improve incumbent solution by flipping pairs of variables 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_twoopt.h" 36 #include "scip/pub_heur.h" 37 #include "scip/pub_lp.h" 38 #include "scip/pub_message.h" 39 #include "scip/pub_misc.h" 40 #include "scip/pub_misc_sort.h" 41 #include "scip/pub_sol.h" 42 #include "scip/pub_var.h" 43 #include "scip/scip_heur.h" 44 #include "scip/scip_lp.h" 45 #include "scip/scip_mem.h" 46 #include "scip/scip_message.h" 47 #include "scip/scip_numerics.h" 48 #include "scip/scip_param.h" 49 #include "scip/scip_prob.h" 50 #include "scip/scip_randnumgen.h" 51 #include "scip/scip_sol.h" 52 #include "scip/scip_solvingstats.h" 53 #include <string.h> 54 55 #define HEUR_NAME "twoopt" 56 #define HEUR_DESC "primal heuristic to improve incumbent solution by flipping pairs of variables" 57 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_ITERATIVE 58 #define HEUR_PRIORITY -20100 59 #define HEUR_FREQ -1 60 #define HEUR_FREQOFS 0 61 #define HEUR_MAXDEPTH -1 62 63 #define HEUR_TIMING SCIP_HEURTIMING_AFTERNODE 64 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */ 65 66 /* default parameter values */ 67 #define DEFAULT_INTOPT FALSE /**< optional integer optimization is applied by default */ 68 #define DEFAULT_WAITINGNODES 0 /**< default number of nodes to wait after current best solution before calling heuristic */ 69 #define DEFAULT_MATCHINGRATE 0.5 /**< default percentage by which two variables have to match in their LP-row set to be 70 * associated as pair by heuristic */ 71 #define DEFAULT_MAXNSLAVES 199 /**< default number of slave candidates for a master variable */ 72 #define DEFAULT_ARRAYSIZE 10 /**< the default array size for temporary arrays */ 73 #define DEFAULT_RANDSEED 37 /**< initial random seed */ 74 75 /* 76 * Data structures 77 */ 78 79 /** primal heuristic data */ 80 struct SCIP_HeurData 81 { 82 int lastsolindex; /**< index of last solution for which heuristic was performed */ 83 SCIP_Real matchingrate; /**< percentage by which two variables have have to match in their LP-row 84 * set to be associated as pair by heuristic */ 85 SCIP_VAR** binvars; /**< Array of binary variables which are sorted with respect to their occurrence 86 * in the LP-rows */ 87 int nbinvars; /**< number of binary variables stored in heuristic array */ 88 int waitingnodes; /**< user parameter to determine number of nodes to wait after last best solution 89 * before calling heuristic */ 90 SCIP_Bool presolved; /**< flag to indicate whether presolving has already been executed */ 91 int* binblockstart; /**< array to store the start indices of each binary block */ 92 int* binblockend; /**< array to store the end indices of each binary block */ 93 int nbinblocks; /**< number of blocks */ 94 95 /* integer variable twoopt data */ 96 SCIP_Bool intopt; /**< parameter to determine if integer 2-opt should be applied */ 97 SCIP_VAR** intvars; /**< array to store the integer variables in non-decreasing order 98 * with respect to their objective coefficient */ 99 int nintvars; /**< the number of integer variables stored in array intvars */ 100 int* intblockstart; /**< array to store the start indices of each binary block */ 101 int* intblockend; /**< array to store the end indices of each binary block */ 102 int nintblocks; /**< number of blocks */ 103 104 SCIP_Bool execute; /**< has presolveTwoOpt detected necessary structure for execution of heuristic? */ 105 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */ 106 int maxnslaves; /**< delimits the maximum number of slave candidates for a master variable */ 107 108 #ifdef SCIP_STATISTIC 109 /* statistics */ 110 int ntotalbinvars; /**< total number of binary variables over all runs */ 111 int ntotalintvars; /**< total number of Integer variables over all runs */ 112 int nruns; /**< counts the number of runs, i.e. the number of initialized 113 * branch and bound processes */ 114 int maxbinblocksize; /**< maximum size of a binary block */ 115 int maxintblocksize; /**< maximum size of an integer block */ 116 int binnblockvars; /**< number of binary variables that appear in blocks */ 117 int binnblocks; /**< number of blocks with at least two variables */ 118 int intnblockvars; /**< number of Integer variables that appear in blocks */ 119 int intnblocks; /**< number of blocks with at least two variables */ 120 int binnexchanges; /**< number of executed changes of binary solution values leading to 121 * improvement in objective function */ 122 int intnexchanges; /**< number of executed changes of Integer solution values leading to improvement in 123 * objective function */ 124 #endif 125 }; 126 127 /** indicator for optimizing for binaries or integer variables */ 128 enum Opttype 129 { 130 OPTTYPE_BINARY = 1, 131 OPTTYPE_INTEGER = 2 132 }; 133 typedef enum Opttype OPTTYPE; 134 135 /** indicator for direction of shifting variables */ 136 enum Direction 137 { 138 DIRECTION_UP = 1, 139 DIRECTION_DOWN = -1, 140 DIRECTION_NONE = 0 141 }; 142 typedef enum Direction DIRECTION; 143 144 /* 145 * Local methods 146 */ 147 148 /** Tries to switch the values of two binary or integer variables and checks feasibility with respect to the LP. 149 * 150 * @todo Adapt method not to copy entire activities array, but only the relevant region. 151 */ 152 static 153 SCIP_RETCODE shiftValues( 154 SCIP* scip, /**< scip instance */ 155 SCIP_VAR* master, /**< first variable of variable pair */ 156 SCIP_VAR* slave, /**< second variable of pair */ 157 SCIP_Real mastersolval, /**< current value of variable1 in solution */ 158 DIRECTION masterdir, /**< the direction into which the master variable has to be shifted */ 159 SCIP_Real slavesolval, /**< current value of variable2 in solution */ 160 DIRECTION slavedir, /**< the direction into which the slave variable has to be shifted */ 161 SCIP_Real shiftval, /**< the value that variables should be shifted by */ 162 SCIP_Real* activities, /**< the LP-row activities */ 163 int nrows, /**< size of activities array */ 164 SCIP_Bool* feasible /**< set to true if method has successfully switched the variable values */ 165 ) 166 { /*lint --e{715}*/ 167 SCIP_COL* col; 168 SCIP_ROW** masterrows; 169 SCIP_ROW** slaverows; 170 SCIP_Real* mastercolvals; 171 SCIP_Real* slavecolvals; 172 int ncolmasterrows; 173 int ncolslaverows; 174 175 assert(scip != NULL); 176 assert(master != NULL); 177 assert(slave != NULL); 178 assert(activities != NULL); 179 assert(SCIPisFeasGT(scip, shiftval, 0.0)); 180 181 assert(SCIPisFeasGE(scip, mastersolval + (int)masterdir * shiftval, SCIPvarGetLbGlobal(master))); 182 assert(SCIPisFeasLE(scip, mastersolval + (int)masterdir * shiftval, SCIPvarGetUbGlobal(master))); 183 184 assert(SCIPisFeasGE(scip, slavesolval + (int)slavedir * shiftval, SCIPvarGetLbGlobal(slave))); 185 assert(SCIPisFeasLE(scip, slavesolval + (int)slavedir * shiftval, SCIPvarGetUbGlobal(slave))); 186 187 /* get variable specific rows and coefficients for both master and slave. */ 188 col = SCIPvarGetCol(master); 189 masterrows = SCIPcolGetRows(col); 190 mastercolvals = SCIPcolGetVals(col); 191 ncolmasterrows = SCIPcolGetNNonz(col); 192 assert(ncolmasterrows == 0 || masterrows != NULL); 193 194 col = SCIPvarGetCol(slave); 195 slaverows = SCIPcolGetRows(col); 196 slavecolvals = SCIPcolGetVals(col); 197 ncolslaverows = SCIPcolGetNNonz(col); 198 assert(ncolslaverows == 0 || slaverows != NULL); 199 200 /* update the activities of the LP rows of the master variable */ 201 for( int i = 0; i < ncolmasterrows && SCIProwGetLPPos(masterrows[i]) >= 0; ++i ) 202 { 203 int rowpos; 204 205 rowpos = SCIProwGetLPPos(masterrows[i]); 206 assert(rowpos < nrows); 207 208 /* skip local rows */ 209 if( rowpos >= 0 && ! SCIProwIsLocal(masterrows[i]) ) 210 activities[rowpos] += mastercolvals[i] * (int)masterdir * shiftval; 211 } 212 213 /* update the activities of the LP rows of the slave variable */ 214 for( int j = 0; j < ncolslaverows && SCIProwGetLPPos(slaverows[j]) >= 0; ++j ) 215 { 216 int rowpos; 217 218 rowpos = SCIProwGetLPPos(slaverows[j]); 219 assert(rowpos < nrows); 220 221 /* skip local rows */ 222 if( rowpos >= 0 && ! SCIProwIsLocal(slaverows[j]) ) 223 { 224 activities[rowpos] += slavecolvals[j] * (int)slavedir * shiftval; 225 assert(SCIPisFeasGE(scip, activities[rowpos], SCIProwGetLhs(slaverows[j]))); 226 assert(SCIPisFeasLE(scip, activities[rowpos], SCIProwGetRhs(slaverows[j]))); 227 } 228 } 229 230 /* in debug mode, the master rows are checked for feasibility which should be granted by the 231 * decision for a shift value */ 232 #ifndef NDEBUG 233 for( int i = 0; i < ncolmasterrows && SCIProwGetLPPos(masterrows[i]) >= 0; ++i ) 234 { 235 /* local rows can be skipped */ 236 if( SCIProwIsLocal(masterrows[i]) ) 237 continue; 238 239 assert(SCIPisFeasGE(scip, activities[SCIProwGetLPPos(masterrows[i])], SCIProwGetLhs(masterrows[i]))); 240 assert(SCIPisFeasLE(scip, activities[SCIProwGetLPPos(masterrows[i])], SCIProwGetRhs(masterrows[i]))); 241 } 242 #endif 243 244 *feasible = TRUE; 245 246 return SCIP_OKAY; 247 } 248 249 /** Compare two variables with respect to their columns. 250 * 251 * Columns are treated as {0,1} vector, where every nonzero entry is treated as '1', and compared to each other 252 * lexicographically. I.e. var1 is < var2 if the corresponding column of var2 has the smaller single nonzero index of 253 * the two columns. This comparison costs O(constraints) in the worst case 254 */ 255 static 256 int varColCompare( 257 SCIP_VAR* var1, /**< left argument of comparison */ 258 SCIP_VAR* var2 /**< right argument of comparison */ 259 ) 260 { 261 SCIP_COL* col1; 262 SCIP_COL* col2; 263 SCIP_ROW** rows1; 264 SCIP_ROW** rows2; 265 int nnonzeros1; 266 int nnonzeros2; 267 268 assert(var1 != NULL); 269 assert(var2 != NULL); 270 271 /* get the necessary row and column data */ 272 col1 = SCIPvarGetCol(var1); 273 col2 = SCIPvarGetCol(var2); 274 rows1 = SCIPcolGetRows(col1); 275 rows2 = SCIPcolGetRows(col2); 276 nnonzeros1 = SCIPcolGetNNonz(col1); 277 nnonzeros2 = SCIPcolGetNNonz(col2); 278 279 assert(nnonzeros1 == 0 || rows1 != NULL); 280 assert(nnonzeros2 == 0 || rows2 != NULL); 281 282 /* loop over the rows, stopped as soon as they differ in one index, 283 * or if counter reaches the end of a variables row set */ 284 for( int i = 0; i < nnonzeros1 && i < nnonzeros2; ++i ) 285 { 286 if( SCIProwGetIndex(rows1[i]) != SCIProwGetIndex(rows2[i]) ) 287 return SCIProwGetIndex(rows1[i]) - SCIProwGetIndex(rows2[i]); 288 } 289 290 /* loop is finished, without differing in one of common row indices, due to loop invariant 291 * variable i reached either nnonzeros1 or nnonzeros2 or both. 292 * one can easily check that the difference of these two numbers always has the desired sign for comparison. */ 293 return nnonzeros2 - nnonzeros1 ; 294 } 295 296 /** implements a comparator to compare two variables with respect to their column entries */ 297 static 298 SCIP_DECL_SORTPTRCOMP(SCIPvarcolComp) 299 { 300 return varColCompare((SCIP_VAR*) elem1, (SCIP_VAR*) elem2); 301 } 302 303 /** checks if two given variables are contained in common LP rows, 304 * returns true if variables share the necessary percentage (matchingrate) of rows. 305 */ 306 static 307 SCIP_Bool checkConstraintMatching( 308 SCIP* scip, /**< current SCIP instance */ 309 SCIP_VAR* var1, /**< first variable */ 310 SCIP_VAR* var2, /**< second variable */ 311 SCIP_Real matchingrate /**< determines the ratio of shared LP rows compared to the total number of 312 * LP-rows each variable appears in */ 313 ) 314 { 315 SCIP_COL* col1; 316 SCIP_COL* col2; 317 SCIP_ROW** rows1; 318 SCIP_ROW** rows2; 319 int nnonzeros1; 320 int nnonzeros2; 321 int i; 322 int j; 323 int nrows1not2; /* the number of LP-rows of variable 1 which variable 2 doesn't appear in */ 324 int nrows2not1; /* vice versa */ 325 int nrowmaximum; 326 int nrowabs; 327 328 assert(var1 != NULL); 329 assert(var2 != NULL); 330 331 /* get the necessary row and column data */ 332 col1 = SCIPvarGetCol(var1); 333 col2 = SCIPvarGetCol(var2); 334 rows1 = SCIPcolGetRows(col1); 335 rows2 = SCIPcolGetRows(col2); 336 nnonzeros1 = SCIPcolGetNNonz(col1); 337 nnonzeros2 = SCIPcolGetNNonz(col2); 338 339 assert(nnonzeros1 == 0 || rows1 != NULL); 340 assert(nnonzeros2 == 0 || rows2 != NULL); 341 342 if( nnonzeros1 == 0 && nnonzeros2 == 0 ) 343 return TRUE; 344 345 /* if matching rate is 0.0, we don't need to check anything */ 346 if( matchingrate == 0.0 ) 347 return TRUE; 348 349 /* initialize the counters for the number of rows not shared. */ 350 nrowmaximum = MAX(nnonzeros1, nnonzeros2); 351 352 nrowabs = ABS(nnonzeros1 - nnonzeros2); 353 nrows1not2 = nrowmaximum - nnonzeros2; 354 nrows2not1 = nrowmaximum - nnonzeros1; 355 356 /* if the numbers of nonzero rows differs too much, w.r.t.matching ratio, the more expensive check over the rows 357 * doesn't have to be applied anymore because the counters for not shared rows can only increase. 358 */ 359 assert(nrowmaximum > 0); 360 361 if( (nrowmaximum - nrowabs) / (SCIP_Real) nrowmaximum < matchingrate ) 362 return FALSE; 363 364 i = 0; 365 j = 0; 366 367 /* loop over all rows and determine number of non-shared rows */ 368 while( i < nnonzeros1 && j < nnonzeros2 ) 369 { 370 /* variables share a common row */ 371 if( SCIProwGetIndex(rows1[i]) == SCIProwGetIndex(rows2[j]) ) 372 { 373 ++i; 374 ++j; 375 } 376 /* variable 1 appears in rows1[i], variable 2 doesn't */ 377 else if( SCIProwGetIndex(rows1[i]) < SCIProwGetIndex(rows2[j]) ) 378 { 379 ++i; 380 ++nrows1not2; 381 } 382 /* variable 2 appears in rows2[j], variable 1 doesn't */ 383 else 384 { 385 ++j; 386 ++nrows2not1; 387 } 388 } 389 390 /* now apply the ratio based comparison, that is if the ratio of shared rows is greater or equal the matching rate 391 * for each variable */ 392 /* nnonzeros1 = 0 or nnonzeros2 = 0 iff matching rate is 0, but in this case, we return TRUE at the beginning */ 393 /* coverity[divide_by_zero] */ 394 return ( SCIPisFeasLE(scip, matchingrate, (nnonzeros1 - nrows1not2) / (SCIP_Real)(nnonzeros1)) || 395 SCIPisFeasLE(scip, matchingrate, (nnonzeros2 - nrows2not1) / (SCIP_Real)(nnonzeros2)) ); /*lint !e795 */ 396 } 397 398 /** Determines a bound by which the absolute solution value of two integer variables can be shifted at most. 399 * 400 * The criterion is the maintenance of feasibility of any global LP row. 401 * The first implementation only considers shifting proportion 1:1, i.e. if master value is shifted by a certain 402 * integer value k downwards, the value of slave is simultaneously shifted by k upwards. 403 */ 404 static 405 SCIP_Real determineBound( 406 SCIP* scip, /**< current scip instance */ 407 SCIP_SOL* sol, /**< current incumbent */ 408 SCIP_VAR* master, /**< current master variable */ 409 DIRECTION masterdirection, /**< the shifting direction of the master variable */ 410 SCIP_VAR* slave, /**< slave variable with same LP_row set as master variable */ 411 DIRECTION slavedirection, /**< the shifting direction of the slave variable */ 412 SCIP_Real* activities, /**< array of LP row activities */ 413 int nrows /**< the number of rows in LP and the size of the activities array */ 414 ) 415 { /*lint --e{715}*/ 416 SCIP_Real masterbound; 417 SCIP_Real slavebound; 418 SCIP_Real bound; 419 SCIP_Real mastersolval; 420 SCIP_Real slavesolval; 421 422 SCIP_COL* col; 423 SCIP_ROW** slaverows; 424 SCIP_ROW** masterrows; 425 SCIP_Real* mastercolvals; 426 SCIP_Real* slavecolvals; 427 int nslaverows; 428 int nmasterrows; 429 int i; 430 int j; 431 432 assert(scip != NULL); 433 assert(sol != NULL); 434 assert(master != NULL); 435 assert(slave != NULL); 436 assert(SCIPvarIsIntegral(master) && SCIPvarIsIntegral(slave)); 437 assert(masterdirection == DIRECTION_UP || masterdirection == DIRECTION_DOWN); 438 assert(slavedirection == DIRECTION_UP || slavedirection == DIRECTION_DOWN); 439 440 mastersolval = SCIPgetSolVal(scip, sol, master); 441 slavesolval = SCIPgetSolVal(scip, sol, slave); 442 443 /* determine the trivial variable bounds for shift */ 444 if( masterdirection == DIRECTION_UP ) 445 { 446 bound = SCIPvarGetUbGlobal(master); 447 masterbound = bound - mastersolval; 448 masterbound = SCIPisFeasLE(scip, mastersolval + ceil(masterbound), bound) ? ceil(masterbound) : floor(masterbound); 449 } 450 else 451 { 452 bound = SCIPvarGetLbGlobal(master); 453 masterbound = mastersolval - bound; 454 masterbound = SCIPisFeasGE(scip, mastersolval - ceil(masterbound), bound) ? ceil(masterbound) : floor(masterbound); 455 } 456 457 if( slavedirection == DIRECTION_UP ) 458 { 459 bound = SCIPvarGetUbGlobal(slave); 460 slavebound = bound - slavesolval; 461 slavebound = SCIPisFeasLE(scip, slavesolval + ceil(slavebound), bound) ? ceil(slavebound) : floor(slavebound); 462 } 463 else 464 { 465 bound = SCIPvarGetLbGlobal(slave); 466 slavebound = slavesolval - bound; 467 slavebound = SCIPisFeasGE(scip, slavesolval - ceil(slavebound), bound) ? ceil(slavebound) : floor(slavebound); 468 } 469 470 bound = MIN(masterbound, slavebound); 471 472 /* due to numerical reasons, bound can be negative -> Return value zero */ 473 if( bound <= 0.0 ) 474 return 0.0; 475 476 /* get the necessary row and and column data for each variable */ 477 col = SCIPvarGetCol(slave); 478 slaverows = SCIPcolGetRows(col); 479 slavecolvals = SCIPcolGetVals(col); 480 nslaverows = SCIPcolGetNNonz(col); 481 482 col = SCIPvarGetCol(master); 483 mastercolvals = SCIPcolGetVals(col); 484 masterrows = SCIPcolGetRows(col); 485 nmasterrows = SCIPcolGetNNonz(col); 486 487 assert(nslaverows == 0 || slavecolvals != NULL); 488 assert(nmasterrows == 0 || mastercolvals != NULL); 489 490 SCIPdebugMsg(scip, " Master: %s with direction %d and %d rows, Slave: %s with direction %d and %d rows \n", SCIPvarGetName(master), 491 (int)masterdirection, nmasterrows, SCIPvarGetName(slave), (int)slavedirection, nslaverows); 492 493 /* loop over all LP rows and determine the maximum integer bound by which both variables 494 * can be shifted without loss of feasibility 495 */ 496 i = 0; 497 j = 0; 498 while( i < nslaverows || j < nmasterrows ) 499 { 500 SCIP_ROW* row; 501 int rowpos; 502 int masterindex; 503 int slaveindex; 504 SCIP_Bool slaveincrement; 505 SCIP_Bool masterincrement; 506 507 /* check if one pointer already reached the end of the respective array */ 508 if( i < nslaverows && SCIProwGetLPPos(slaverows[i]) == -1 ) 509 { 510 SCIPdebugMsg(scip, " Slaverow %s is not in LP (i=%d, j=%d)\n", SCIProwGetName(slaverows[i]), i, j); 511 i = nslaverows; 512 continue; 513 } 514 if( j < nmasterrows && SCIProwGetLPPos(masterrows[j]) == -1 ) 515 { 516 SCIPdebugMsg(scip, " Masterrow %s is not in LP (i=%d, j=%d) \n", SCIProwGetName(masterrows[j]), i, j); 517 j = nmasterrows; 518 continue; 519 } 520 521 slaveincrement = FALSE; 522 /* If one counter has already reached its limit, assign a huge number to the corresponding 523 * row index to simulate an always greater row position. */ 524 if( i < nslaverows ) 525 slaveindex = SCIProwGetIndex(slaverows[i]); 526 else 527 slaveindex = INT_MAX; 528 529 if( j < nmasterrows ) 530 masterindex = SCIProwGetIndex(masterrows[j]); 531 else 532 masterindex = INT_MAX; 533 534 assert(0 <= slaveindex && 0 <= masterindex); 535 536 assert(slaveindex < INT_MAX || masterindex < INT_MAX); 537 538 /* the current row is the one with the smaller index */ 539 if( slaveindex <= masterindex ) 540 { 541 rowpos = SCIProwGetLPPos(slaverows[i]); 542 row = slaverows[i]; 543 slaveincrement = TRUE; 544 masterincrement = (slaveindex == masterindex); 545 } 546 else 547 { 548 assert(j < nmasterrows); 549 550 rowpos = SCIProwGetLPPos(masterrows[j]); 551 row = masterrows[j]; 552 masterincrement = TRUE; 553 } 554 assert(row != NULL); 555 556 /* only global rows need to be valid */ 557 if( rowpos >= 0 && !SCIProwIsLocal(row) ) 558 { 559 SCIP_Real effect; 560 SCIP_Real side; 561 SCIP_Bool left; 562 563 /* effect is the effect on the row activity by shifting the variables by 1 in the respective directions */ 564 effect = 0.0; 565 if( slaveindex <= masterindex ) 566 effect += (slavecolvals[i] * (int)slavedirection); 567 if( masterindex <= slaveindex ) 568 effect += (mastercolvals[j] * (int)masterdirection); 569 left = effect < 0.0; 570 side = left ? SCIProwGetLhs(row) : SCIProwGetRhs(row); 571 572 /* only non-zero effects and finite bounds need to be considered */ 573 if( !SCIPisZero(scip, effect) && !SCIPisInfinity(scip, left ? -side : side) ) 574 { 575 SCIP_Real newval; 576 577 /* effect does not equal zero, the bound is determined as maximum positive integer such that 578 * feasibility of this constraint is maintained 579 */ 580 assert( rowpos < nrows ); 581 assert( SCIPisFeasGE(scip, activities[rowpos], SCIProwGetLhs(row)) && SCIPisFeasLE(scip, activities[rowpos], SCIProwGetRhs(row)) ); 582 assert( effect ); 583 584 SCIPdebugMsg(scip, " %g <= %g <= %g, bound = %g, effect = %g (%g * %d + %g * %d) (i=%d,j=%d)\n", 585 SCIProwGetLhs(row), activities[rowpos], SCIProwGetRhs(row), bound, effect, 586 slaveindex <= masterindex ? slavecolvals[i] : 0.0, (int)slavedirection, 587 masterindex <= slaveindex ? mastercolvals[j] : 0.0, (int)masterdirection, i, j); 588 589 newval = (side - activities[rowpos]) / effect; 590 591 /* update shifting distance */ 592 if( newval < bound ) 593 { 594 SCIP_Real activity; 595 596 activity = activities[rowpos] + effect * ceil(newval); 597 598 /* ensure that shifting preserves feasibility */ 599 if( ( left && SCIPisFeasGE(scip, activity, side) ) || ( !left && SCIPisFeasLE(scip, activity, side) ) ) 600 bound = ceil(newval); 601 else 602 bound = floor(newval); 603 604 /* due to numerical reasons, bound can be negative. A variable shift by a negative bound is not desired by 605 * the heuristic -> Return value zero */ 606 if( bound <= 0.0 ) 607 return 0.0; 608 } 609 610 assert( SCIPisFeasGE(scip, activities[rowpos] + effect * bound, SCIProwGetLhs(row)) && SCIPisFeasLE(scip, activities[rowpos] + effect * bound, SCIProwGetRhs(row)) ); 611 assert( SCIPisFeasIntegral(scip, bound) ); 612 } 613 else 614 { 615 SCIPdebugMsg(scip, " No influence of row %s, effect %g, master coeff: %g slave coeff: %g (i=%d, j=%d)\n", 616 SCIProwGetName(row), effect, mastercolvals[j], slavecolvals[i], i, j); 617 } 618 } 619 else 620 { 621 SCIPdebugMsg(scip, " Row %s is local.\n", SCIProwGetName(row)); 622 } 623 624 /* increase the counters which belong to the corresponding row. Both counters are increased by 625 * 1 iff rowpos1 <= rowpos2 <= rowpos1 */ 626 if( slaveincrement ) 627 ++i; 628 if( masterincrement ) 629 ++j; 630 } 631 632 /* we must not shift variables to infinity */ 633 return SCIPisInfinity(scip, bound + MAX((int)masterdirection * mastersolval, (int)slavedirection * slavesolval)) ? 0.0 : bound; 634 } 635 636 637 /** Disposes variable with no heuristic relevancy, e.g., due to a fixed solution value, from its neighborhood block. 638 * 639 * The affected neighborhood block is reduced by 1. 640 */ 641 static 642 void disposeVariable( 643 SCIP_VAR** vars, /**< problem variables */ 644 int* blockend, /**< contains end index of block */ 645 int varindex /**< variable index */ 646 ) 647 { 648 assert(blockend != NULL); 649 assert(varindex <= *blockend); 650 651 vars[varindex] = vars[*blockend]; 652 --(*blockend); 653 } 654 655 /** realizes the presolve independently from type of variables it's applied to */ 656 static 657 SCIP_RETCODE innerPresolve( 658 SCIP* scip, /**< current scip */ 659 SCIP_VAR** vars, /**< problem vars */ 660 SCIP_VAR*** varspointer, /**< pointer to heuristic specific variable memory */ 661 int nvars, /**< the number of variables */ 662 int* nblocks, /**< pointer to store the number of detected blocks */ 663 int* maxblocksize, /**< maximum size of a block */ 664 int* nblockvars, /**< pointer to store the number of block variables */ 665 int** blockstart, /**< pointer to store the array of block start indices */ 666 int** blockend, /**< pointer to store the array of block end indices */ 667 SCIP_HEUR* heur, /**< the heuristic */ 668 SCIP_HEURDATA* heurdata /**< the heuristic data */ 669 ) 670 { 671 int startindex; 672 673 assert(scip != NULL); 674 assert(vars != NULL); 675 assert(nvars >= 2); 676 assert(nblocks != NULL); 677 assert(nblockvars != NULL); 678 assert(blockstart != NULL); 679 assert(blockend != NULL); 680 assert(heur != NULL); 681 assert(heurdata != NULL); 682 683 /* allocate the heuristic specific variables */ 684 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, varspointer, vars, nvars)); 685 686 /* sort the variables with respect to their columns */ 687 SCIPsortPtr((void**)(*varspointer), SCIPvarcolComp, nvars); 688 689 /* start determining blocks, i.e. a set of at least two variables which share most of their row set. 690 * If there is none, heuristic does not need to be executed. 691 */ 692 startindex = 0; 693 *nblocks = 0; 694 *nblockvars = 0; 695 696 SCIP_CALL( SCIPallocBlockMemoryArray(scip, blockstart, nvars/2) ); 697 SCIP_CALL( SCIPallocBlockMemoryArray(scip, blockend, nvars/2) ); 698 699 /* loop over variables and compare neighbors */ 700 for( int v = 1; v < nvars; ++v ) 701 { 702 if( !checkConstraintMatching(scip, (*varspointer)[startindex], (*varspointer)[v], heurdata->matchingrate) ) 703 { 704 /* current block has its last variable at position v-1. If v differs from startindex by at least 2, 705 * a block is detected. Update the data correspondingly */ 706 if( v - startindex >= 2 ) 707 { 708 assert(*nblocks < nvars/2); 709 (*nblockvars) += v - startindex; 710 (*maxblocksize) = MAX((*maxblocksize), v - startindex); 711 (*blockstart)[*nblocks] = startindex; 712 (*blockend)[*nblocks] = v - 1; 713 (*nblocks)++; 714 } 715 startindex = v; 716 } 717 else if( v == nvars - 1 && v - startindex >= 1 ) 718 { 719 assert(*nblocks < nvars/2); 720 (*nblockvars) += v - startindex + 1; 721 (*maxblocksize) = MAX((*maxblocksize), v - startindex + 1); 722 (*blockstart)[*nblocks] = startindex; 723 (*blockend)[*nblocks] = v; 724 (*nblocks)++; 725 } 726 } 727 728 /* reallocate memory with respect to the number of found blocks; if there were none, free the memory */ 729 if( *nblocks > 0 ) 730 { 731 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, blockstart, nvars/2, *nblocks) ); 732 SCIP_CALL( SCIPreallocBlockMemoryArray(scip, blockend, nvars/2, *nblocks) ); 733 } 734 else 735 { 736 SCIPfreeBlockMemoryArray(scip, blockstart, nvars/2); 737 SCIPfreeBlockMemoryArray(scip, blockend, nvars/2); 738 739 *blockstart = NULL; 740 *blockend = NULL; 741 } 742 743 return SCIP_OKAY; 744 } 745 746 /** initializes the required structures for execution of heuristic. 747 * 748 * If objective coefficient functions are not all equal, each Binary and Integer variables are sorted 749 * into heuristic-specific arrays with respect to their lexicographical column order, 750 * where every zero in a column is interpreted as zero and every nonzero as '1'. 751 * After the sorting, the variables are compared with respect to user parameter matchingrate and 752 * the heuristic specific blocks are determined. 753 */ 754 static 755 SCIP_RETCODE presolveTwoOpt( 756 SCIP* scip, /**< current scip instance */ 757 SCIP_HEUR* heur, /**< heuristic */ 758 SCIP_HEURDATA* heurdata /**< the heuristic data */ 759 ) 760 { 761 int nbinvars; 762 int nintvars; 763 int nvars; 764 SCIP_VAR** vars; 765 int nbinblockvars = 0; 766 int nintblockvars; 767 int maxbinblocksize = 0; 768 int maxintblocksize; 769 770 assert(scip != NULL); 771 assert(heurdata != NULL); 772 773 /* ensure that method is not executed if presolving was already applied once in current branch and bound process */ 774 if( heurdata->presolved ) 775 return SCIP_OKAY; 776 777 /* get necessary variable information, i.e. number of binary and integer variables */ 778 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) ); 779 780 /* if number of binary problem variables exceeds 2, they are subject to 2-optimization algorithm, hence heuristic 781 * calls innerPresolve method to detect necessary structures. */ 782 if( nbinvars >= 2 ) 783 { 784 SCIP_CALL( innerPresolve(scip, vars, &(heurdata->binvars), nbinvars, &(heurdata->nbinblocks), &maxbinblocksize, 785 &nbinblockvars, &(heurdata->binblockstart), &(heurdata->binblockend), heur, heurdata) ); 786 } 787 788 heurdata->nbinvars = nbinvars; 789 heurdata->execute = nbinvars > 1 && heurdata->nbinblocks > 0; 790 791 #ifdef SCIP_STATISTIC 792 /* update statistics */ 793 heurdata->binnblocks += (heurdata->nbinblocks); 794 heurdata->binnblockvars += nbinblockvars; 795 heurdata->ntotalbinvars += nbinvars; 796 heurdata->maxbinblocksize = MAX(maxbinblocksize, heurdata->maxbinblocksize); 797 798 SCIPstatisticMessage(" Twoopt BINARY presolving finished with <%d> blocks, <%d> block variables \n", 799 heurdata->nbinblocks, nbinblockvars); 800 #endif 801 802 if( heurdata->intopt && nintvars > 1 ) 803 { 804 SCIP_CALL( innerPresolve(scip, &(vars[nbinvars]), &(heurdata->intvars), nintvars, &(heurdata->nintblocks), &maxintblocksize, 805 &nintblockvars, &(heurdata->intblockstart), &(heurdata->intblockend), 806 heur, heurdata) ); 807 808 heurdata->execute = heurdata->execute || heurdata->nintblocks > 0; 809 810 #ifdef SCIP_STATISTIC 811 /* update statistics */ 812 heurdata->intnblocks += heurdata->nintblocks; 813 heurdata->intnblockvars += nintblockvars; 814 heurdata->ntotalintvars += nintvars; 815 heurdata->maxintblocksize = MAX(maxintblocksize, heurdata->maxintblocksize); 816 SCIPstatisticMessage(" Twoopt Integer presolving finished with <%d> blocks, <%d> block variables \n", 817 heurdata->nintblocks, nintblockvars); 818 SCIPstatisticMessage(" INTEGER coefficients are all equal \n"); 819 #endif 820 } 821 heurdata->nintvars = nintvars; 822 823 /* presolving is finished, heuristic data is updated*/ 824 heurdata->presolved = TRUE; 825 SCIPheurSetData(heur, heurdata); 826 827 return SCIP_OKAY; 828 } 829 830 /* 831 * Callback methods of primal heuristic 832 */ 833 834 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */ 835 static 836 SCIP_DECL_HEURCOPY(heurCopyTwoopt) 837 { /*lint --e{715}*/ 838 assert(scip != NULL); 839 assert(heur != NULL); 840 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 841 842 /* call inclusion method of primal heuristic */ 843 SCIP_CALL( SCIPincludeHeurTwoopt(scip) ); 844 845 return SCIP_OKAY; 846 } 847 848 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */ 849 static 850 SCIP_DECL_HEURFREE(heurFreeTwoopt) 851 { /*lint --e{715}*/ 852 SCIP_HEURDATA* heurdata; 853 854 assert(heur != NULL); 855 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 856 assert(scip != NULL); 857 858 /* free heuristic data */ 859 heurdata = SCIPheurGetData(heur); 860 assert(heurdata != NULL); 861 862 SCIPfreeBlockMemory(scip, &heurdata); 863 SCIPheurSetData(heur, NULL); 864 865 return SCIP_OKAY; 866 } 867 868 /** initialization method of primal heuristic (called after problem was transformed) */ 869 static 870 SCIP_DECL_HEURINIT(heurInitTwoopt) 871 { 872 SCIP_HEURDATA* heurdata; 873 assert(heur != NULL); 874 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 875 assert(scip != NULL); 876 877 heurdata = SCIPheurGetData(heur); 878 assert(heurdata != NULL); 879 880 /* heuristic has not run yet, all heuristic specific data is set to initial values */ 881 heurdata->nbinvars = 0; 882 heurdata->nintvars = 0; 883 heurdata->lastsolindex = -1; 884 heurdata->presolved = FALSE; 885 heurdata->nbinblocks = 0; 886 heurdata->nintblocks = 0; 887 888 /* create random number generator */ 889 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, 890 DEFAULT_RANDSEED, TRUE) ); 891 892 #ifdef SCIP_STATISTIC 893 /* initialize statistics */ 894 heurdata->binnexchanges = 0; 895 heurdata->intnexchanges = 0; 896 heurdata->binnblockvars = 0; 897 heurdata->intnblockvars = 0; 898 heurdata->binnblocks = 0; 899 heurdata->intnblocks = 0; 900 901 heurdata->maxbinblocksize = 0; 902 heurdata->maxintblocksize = 0; 903 904 heurdata->ntotalbinvars = 0; 905 heurdata->ntotalintvars = 0; 906 heurdata->nruns = 0; 907 #endif 908 909 /* all pointers are initially set to NULL. Since presolving 910 * of the heuristic requires a lot of calculation time on some instances, 911 * but might not be needed e.g. if problem is infeasible, presolving is applied 912 * when heuristic is executed for the first time 913 */ 914 heurdata->binvars = NULL; 915 heurdata->intvars = NULL; 916 heurdata->binblockstart = NULL; 917 heurdata->binblockend = NULL; 918 heurdata->intblockstart = NULL; 919 heurdata->intblockend = NULL; 920 921 SCIPheurSetData(heur, heurdata); 922 923 return SCIP_OKAY; 924 } 925 926 /* Realizes the 2-optimization algorithm, which tries to improve incumbent solution 927 * by shifting pairs of variables which share a common row set. 928 */ 929 static 930 SCIP_RETCODE optimize( 931 SCIP* scip, /**< current SCIP instance */ 932 SCIP_SOL* worksol, /**< working solution */ 933 SCIP_VAR** vars, /**< binary or integer variables */ 934 int* blockstart, /**< contains start indices of blocks */ 935 int* blockend, /**< contains end indices of blocks */ 936 int nblocks, /**< the number of blocks */ 937 OPTTYPE opttype, /**< are binaries or integers optimized */ 938 SCIP_Real* activities, /**< the LP-row activities */ 939 int nrows, /**< the number of LP rows */ 940 SCIP_Bool* improvement, /**< was there a successful shift? */ 941 SCIP_Bool* varboundserr, /**< has the current incumbent already been cut off */ 942 SCIP_HEURDATA* heurdata /**< the heuristic data */ 943 ) 944 { /*lint --e{715}*/ 945 SCIP_Real* objchanges; 946 SCIP_VAR** bestmasters; 947 SCIP_VAR** bestslaves; 948 int* bestdirections; 949 int arraysize; 950 int npairs = 0; 951 952 assert(scip != NULL); 953 assert(nblocks > 0); 954 assert(blockstart != NULL && blockend != NULL); 955 assert(varboundserr != NULL); 956 assert(activities != NULL); 957 assert(worksol != NULL); 958 assert(improvement != NULL); 959 960 *varboundserr = FALSE; 961 962 SCIP_CALL( SCIPallocBufferArray(scip, &bestmasters, DEFAULT_ARRAYSIZE) ); 963 SCIP_CALL( SCIPallocBufferArray(scip, &bestslaves, DEFAULT_ARRAYSIZE) ); 964 SCIP_CALL( SCIPallocBufferArray(scip, &objchanges, DEFAULT_ARRAYSIZE) ); 965 SCIP_CALL( SCIPallocBufferArray(scip, &bestdirections, DEFAULT_ARRAYSIZE) ); 966 arraysize = DEFAULT_ARRAYSIZE; 967 968 /* iterate over blocks */ 969 for( int b = 0; b < nblocks; ++b ) 970 { 971 int blocklen; 972 973 blocklen = blockend[b] - blockstart[b] + 1; 974 975 /* iterate over variables in current block */ 976 for( int m = 0; m < blocklen; ++m ) 977 { 978 /* determine the new master variable for heuristic's optimization method */ 979 SCIP_VAR* master; 980 SCIP_Real masterobj; 981 SCIP_Real mastersolval; 982 SCIP_Real bestimprovement; 983 SCIP_Real bestbound; 984 int bestslavepos; 985 int firstslave; 986 int nslaves; 987 int bestdirection; 988 DIRECTION bestmasterdir; 989 DIRECTION bestslavedir; 990 991 master = vars[blockstart[b] + m]; /*lint !e679*/ 992 masterobj = SCIPvarGetObj(master); 993 mastersolval = SCIPgetSolVal(scip, worksol, master); 994 995 /* due to cuts or fixings of solution values, worksol might not be feasible w.r.t. its bounds. 996 * Exit method in that case. */ 997 if( SCIPisFeasGT(scip, mastersolval, SCIPvarGetUbGlobal(master)) || SCIPisFeasLT(scip, mastersolval, SCIPvarGetLbGlobal(master)) ) 998 { 999 *varboundserr = TRUE; 1000 SCIPdebugMsg(scip, "Solution has violated variable bounds for var %s: %g <= %g <= %g \n", 1001 SCIPvarGetName(master), SCIPvarGetLbGlobal(master), mastersolval, SCIPvarGetUbGlobal(master)); 1002 goto TERMINATE; 1003 } 1004 1005 /* if variable has fixed solution value, it is deleted from heuristic array */ 1006 if( SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(master), SCIPvarGetLbGlobal(master)) ) 1007 { 1008 disposeVariable(vars, &(blockend[b]), blockstart[b] + m); 1009 --blocklen; 1010 continue; 1011 } 1012 else if( SCIPvarGetStatus(master) != SCIP_VARSTATUS_COLUMN ) 1013 continue; 1014 1015 assert(SCIPisFeasIntegral(scip, mastersolval)); 1016 1017 assert(opttype == OPTTYPE_INTEGER || (SCIPisFeasLE(scip, mastersolval, 1.0) || SCIPisFeasGE(scip, mastersolval, 0.0))); 1018 1019 /* initialize the data of the best available shift */ 1020 bestimprovement = 0.0; 1021 bestslavepos = -1; 1022 bestbound = 0.0; 1023 bestmasterdir = DIRECTION_NONE; 1024 bestslavedir = DIRECTION_NONE; 1025 bestdirection = -1; 1026 1027 /* in blocks with more than heurdata->maxnslaves variables, a slave candidate region is chosen */ 1028 if( heurdata->maxnslaves >= 0 && blocklen > heurdata->maxnslaves ) 1029 firstslave = SCIPrandomGetInt(heurdata->randnumgen, blockstart[b] + m, blockend[b]); 1030 else 1031 firstslave = blockstart[b] + m + 1; 1032 1033 nslaves = MIN((heurdata->maxnslaves == -1 ? INT_MAX : heurdata->maxnslaves), blocklen); 1034 1035 /* Loop over block and determine a slave shift candidate for master variable. 1036 * If more than one candidate is available, choose the shift which improves objective function 1037 * the most. */ 1038 for( int s = 0; s < nslaves; ++s ) 1039 { 1040 SCIP_VAR* slave; 1041 SCIP_Real slaveobj; 1042 SCIP_Real slavesolval; 1043 SCIP_Real changedobj; 1044 SCIP_Real diffdirbound; 1045 SCIP_Real equaldirbound; 1046 int directions; 1047 int slaveindex; 1048 1049 slaveindex = (firstslave + s - blockstart[b]) % blocklen; 1050 slaveindex += blockstart[b]; 1051 1052 /* in case of a small block, we do not want to try possible pairings twice */ 1053 if( (blocklen <= heurdata->maxnslaves || heurdata->maxnslaves == -1) && slaveindex < blockstart[b] + m ) 1054 break; 1055 /* master and slave should not be the same variable */ 1056 if( slaveindex == blockstart[b] + m ) 1057 continue; 1058 1059 /* get the next slave variable */ 1060 slave = vars[slaveindex]; 1061 slaveobj = SCIPvarGetObj(slave); 1062 slavesolval = SCIPgetSolVal(scip, worksol, slave); 1063 changedobj = 0.0; 1064 1065 assert(SCIPvarGetType(master) == SCIPvarGetType(slave)); 1066 assert(SCIPisFeasIntegral(scip, slavesolval)); 1067 assert(opttype == OPTTYPE_INTEGER || (SCIPisFeasLE(scip, mastersolval, 1.0) || SCIPisFeasGE(scip, mastersolval, 0.0))); 1068 1069 /* solution is not feasible w.r.t. the variable bounds, stop optimization in this case */ 1070 if( SCIPisFeasGT(scip, slavesolval, SCIPvarGetUbGlobal(slave)) || SCIPisFeasLT(scip, slavesolval, SCIPvarGetLbGlobal(slave)) ) 1071 { 1072 *varboundserr = TRUE; 1073 SCIPdebugMsg(scip, "Solution has violated variable bounds for var %s: %g <= %g <= %g \n", 1074 SCIPvarGetName(slave), SCIPvarGetLbGlobal(slave), slavesolval, SCIPvarGetUbGlobal(slave)); 1075 goto TERMINATE; 1076 } 1077 1078 /* if solution value of the variable is fixed, delete it from the remaining candidates in the block */ 1079 if( SCIPisFeasEQ(scip, SCIPvarGetUbGlobal(slave), SCIPvarGetLbGlobal(slave) ) ) 1080 { 1081 disposeVariable(vars, &(blockend[b]), slaveindex); 1082 --blocklen; 1083 continue; 1084 } 1085 else if( SCIPvarGetStatus(master) != SCIP_VARSTATUS_COLUMN ) 1086 continue; 1087 1088 /* determine the shifting direction to improve the objective function */ 1089 /* The heuristic chooses the shifting direction and the corresponding maximum nonnegative 1090 * integer shift value for the two variables which preserves feasibility and improves 1091 * the objective function value. */ 1092 directions = -1; 1093 diffdirbound = 0.0; 1094 equaldirbound = 0.0; 1095 1096 if( SCIPisPositive(scip, slaveobj - masterobj) ) 1097 { 1098 diffdirbound = determineBound(scip, worksol, master, DIRECTION_UP, slave, DIRECTION_DOWN, activities, nrows); 1099 directions = 2; 1100 /* the improvement of objective function is calculated */ 1101 changedobj = (masterobj - slaveobj) * diffdirbound; 1102 } 1103 else if( SCIPisPositive(scip, masterobj - slaveobj) ) 1104 { 1105 diffdirbound = determineBound(scip, worksol, master, DIRECTION_DOWN, slave, DIRECTION_UP, activities, nrows); 1106 directions = 1; 1107 changedobj = (slaveobj - masterobj) * diffdirbound; 1108 } 1109 1110 if( SCIPisPositive(scip, -(masterobj + slaveobj)) ) 1111 { 1112 equaldirbound = determineBound(scip, worksol, master, DIRECTION_UP, slave, DIRECTION_UP, activities, nrows); 1113 if( (masterobj + slaveobj) * equaldirbound < changedobj ) 1114 { 1115 changedobj = (masterobj + slaveobj) * equaldirbound; 1116 directions = 3; 1117 } 1118 } 1119 else if( SCIPisPositive(scip, masterobj + slaveobj) ) 1120 { 1121 equaldirbound = determineBound(scip, worksol, master, DIRECTION_DOWN, slave, DIRECTION_DOWN, activities, nrows); 1122 if( -(masterobj + slaveobj) * equaldirbound < changedobj ) 1123 { 1124 changedobj = -(masterobj + slaveobj) * equaldirbound; 1125 directions = 0; 1126 } 1127 } 1128 assert(SCIPisFeasIntegral(scip, equaldirbound)); 1129 assert(SCIPisFeasIntegral(scip, diffdirbound)); 1130 assert(SCIPisFeasGE(scip, equaldirbound, 0.0)); 1131 assert(SCIPisFeasGE(scip, diffdirbound, 0.0)); 1132 1133 /* choose the candidate which improves the objective function the most */ 1134 if( (SCIPisFeasGT(scip, equaldirbound, 0.0) || SCIPisFeasGT(scip, diffdirbound, 0.0)) 1135 && changedobj < bestimprovement ) 1136 { 1137 bestimprovement = changedobj; 1138 bestslavepos = slaveindex; 1139 bestdirection = directions; 1140 1141 /* the most promising shift, i.e., the one which can improve the objective 1142 * the most, is recorded by the integer 'directions'. It is recovered via the use 1143 * of a binary representation of the four different combinations for the shifting directions 1144 * of two variables */ 1145 if( directions / 2 == 1 ) 1146 bestmasterdir = DIRECTION_UP; 1147 else 1148 bestmasterdir = DIRECTION_DOWN; 1149 1150 if( directions % 2 == 1 ) 1151 bestslavedir = DIRECTION_UP; 1152 else 1153 bestslavedir = DIRECTION_DOWN; 1154 1155 if( bestmasterdir == bestslavedir ) 1156 bestbound = equaldirbound; 1157 else 1158 bestbound = diffdirbound; 1159 } 1160 } 1161 1162 /* choose the most promising candidate, if one exists */ 1163 if( bestslavepos >= 0 ) 1164 { 1165 if( npairs == arraysize ) 1166 { 1167 SCIP_CALL( SCIPreallocBufferArray(scip, &bestmasters, 2 * arraysize) ); 1168 SCIP_CALL( SCIPreallocBufferArray(scip, &bestslaves, 2 * arraysize) ); 1169 SCIP_CALL( SCIPreallocBufferArray(scip, &objchanges, 2 * arraysize) ); 1170 SCIP_CALL( SCIPreallocBufferArray(scip, &bestdirections, 2 * arraysize) ); 1171 arraysize = 2 * arraysize; 1172 } 1173 assert( npairs < arraysize ); 1174 1175 bestmasters[npairs] = master; 1176 bestslaves[npairs] = vars[bestslavepos]; 1177 objchanges[npairs] = ((int)bestslavedir * SCIPvarGetObj(bestslaves[npairs]) + (int)bestmasterdir * masterobj) * bestbound; 1178 bestdirections[npairs] = bestdirection; 1179 1180 assert(objchanges[npairs] < 0); 1181 1182 SCIPdebugMsg(scip, " Saved candidate pair {%s=%g, %s=%g} with objectives <%g>, <%g> to be set to {%g, %g} %d\n", 1183 SCIPvarGetName(master), mastersolval, SCIPvarGetName(bestslaves[npairs]), SCIPgetSolVal(scip, worksol, bestslaves[npairs]) , 1184 masterobj, SCIPvarGetObj(bestslaves[npairs]), mastersolval + (int)bestmasterdir * bestbound, SCIPgetSolVal(scip, worksol, bestslaves[npairs]) 1185 + (int)bestslavedir * bestbound, bestdirections[npairs]); 1186 1187 ++npairs; 1188 } 1189 } 1190 } 1191 1192 if( npairs == 0 ) 1193 goto TERMINATE; 1194 1195 SCIPsortRealPtrPtrInt(objchanges, (void**)bestmasters, (void**)bestslaves, bestdirections, npairs); 1196 1197 for( int b = 0; b < npairs; ++b ) 1198 { 1199 SCIP_VAR* master; 1200 SCIP_VAR* slave; 1201 SCIP_Real mastersolval; 1202 SCIP_Real slavesolval; 1203 SCIP_Real masterobj; 1204 SCIP_Real slaveobj; 1205 SCIP_Real bound; 1206 DIRECTION masterdir; 1207 DIRECTION slavedir; 1208 1209 master = bestmasters[b]; 1210 slave = bestslaves[b]; 1211 mastersolval = SCIPgetSolVal(scip, worksol, master); 1212 slavesolval = SCIPgetSolVal(scip, worksol, slave); 1213 masterobj =SCIPvarGetObj(master); 1214 slaveobj = SCIPvarGetObj(slave); 1215 1216 assert(0 <= bestdirections[b] && bestdirections[b] < 4); 1217 1218 if( bestdirections[b] / 2 == 1 ) 1219 masterdir = DIRECTION_UP; 1220 else 1221 masterdir = DIRECTION_DOWN; 1222 1223 if( bestdirections[b] % 2 == 1 ) 1224 slavedir = DIRECTION_UP; 1225 else 1226 slavedir = DIRECTION_DOWN; 1227 1228 bound = determineBound(scip, worksol, master, masterdir, slave, slavedir, activities, nrows); 1229 1230 if( !SCIPisZero(scip, bound) ) 1231 { 1232 SCIP_Bool feasible; 1233 #ifndef NDEBUG 1234 SCIP_Real changedobj; 1235 #endif 1236 1237 SCIPdebugMsg(scip, " Promising candidates {%s=%g, %s=%g} with objectives <%g>, <%g> to be set to {%g, %g}\n", 1238 SCIPvarGetName(master), mastersolval, SCIPvarGetName(slave), slavesolval, 1239 masterobj, slaveobj, mastersolval + (int)masterdir * bound, slavesolval + (int)slavedir * bound); 1240 1241 #ifndef NDEBUG 1242 /* the improvement of objective function is calculated */ 1243 changedobj = ((int)slavedir * slaveobj + (int)masterdir * masterobj) * bound; 1244 assert( SCIPisPositive(scip, -changedobj) ); 1245 #endif 1246 1247 assert(SCIPvarGetStatus(master) == SCIP_VARSTATUS_COLUMN && SCIPvarGetStatus(slave) == SCIP_VARSTATUS_COLUMN); 1248 /* try to change the solution values of the variables */ 1249 feasible = FALSE; 1250 SCIP_CALL( shiftValues(scip, master, slave, mastersolval, masterdir, slavesolval, slavedir, bound, 1251 activities, nrows, &feasible) ); 1252 1253 if( feasible ) 1254 { 1255 /* The variables' solution values were successfully shifted and can hence be updated. */ 1256 assert(SCIPisFeasIntegral(scip, mastersolval + ((int)masterdir * bound))); 1257 assert(SCIPisFeasIntegral(scip, slavesolval + ((int)slavedir * bound))); 1258 1259 SCIP_CALL( SCIPsetSolVal(scip, worksol, master, mastersolval + (int)masterdir * bound) ); 1260 SCIP_CALL( SCIPsetSolVal(scip, worksol, slave, slavesolval + (int)slavedir * bound) ); 1261 SCIPdebugMsg(scip, " Feasible shift: <%s>[%g, %g] (obj: %f) <%f> --> <%f>\n", 1262 SCIPvarGetName(master), SCIPvarGetLbGlobal(master), SCIPvarGetUbGlobal(master), masterobj, mastersolval, SCIPgetSolVal(scip, worksol, master)); 1263 SCIPdebugMsg(scip, " <%s>[%g, %g] (obj: %f) <%f> --> <%f>\n", 1264 SCIPvarGetName(slave), SCIPvarGetLbGlobal(slave), SCIPvarGetUbGlobal(slave), slaveobj, slavesolval, SCIPgetSolVal(scip, worksol, slave)); 1265 1266 #ifdef SCIP_STATISTIC 1267 /* update statistics */ 1268 if( opttype == OPTTYPE_BINARY ) 1269 ++(heurdata->binnexchanges); 1270 else 1271 ++(heurdata->intnexchanges); 1272 #endif 1273 1274 *improvement = TRUE; 1275 } 1276 } 1277 } 1278 TERMINATE: 1279 SCIPfreeBufferArray(scip, &bestdirections); 1280 SCIPfreeBufferArray(scip, &objchanges); 1281 SCIPfreeBufferArray(scip, &bestslaves); 1282 SCIPfreeBufferArray(scip, &bestmasters); 1283 1284 return SCIP_OKAY; 1285 } 1286 1287 /** deinitialization method of primal heuristic (called before transformed problem is freed) */ 1288 static 1289 SCIP_DECL_HEUREXIT(heurExitTwoopt) 1290 { 1291 SCIP_HEURDATA* heurdata; 1292 1293 heurdata = SCIPheurGetData(heur); 1294 assert(heurdata != NULL); 1295 1296 /*ensure that initialization was successful */ 1297 assert(heurdata->nbinvars <= 1 || heurdata->binvars != NULL); 1298 1299 #ifdef SCIP_STATISTIC 1300 /* print relevant statistics to console */ 1301 SCIPstatisticMessage( 1302 " Twoopt Binary Statistics : " 1303 "%6.2g %6.2g %4.2g %4.0g %6d (blocks/run, variables/run, varpercentage, avg. block size, max block size) \n", 1304 heurdata->nruns == 0 ? 0.0 : (SCIP_Real)heurdata->binnblocks/(heurdata->nruns), 1305 heurdata->nruns == 0 ? 0.0 : (SCIP_Real)heurdata->binnblockvars/(heurdata->nruns), 1306 heurdata->ntotalbinvars == 0 ? 0.0 : (SCIP_Real)heurdata->binnblockvars/(heurdata->ntotalbinvars) * 100.0, 1307 heurdata->binnblocks == 0 ? 0.0 : heurdata->binnblockvars/(SCIP_Real)(heurdata->binnblocks), 1308 heurdata->maxbinblocksize); 1309 1310 SCIPstatisticMessage( 1311 " Twoopt Integer statistics : " 1312 "%6.2g %6.2g %4.2g %4.0g %6d (blocks/run, variables/run, varpercentage, avg block size, max block size) \n", 1313 heurdata->nruns == 0 ? 0.0 : (SCIP_Real)heurdata->intnblocks/(heurdata->nruns), 1314 heurdata->nruns == 0 ? 0.0 : (SCIP_Real)heurdata->intnblockvars/(heurdata->nruns), 1315 heurdata->ntotalintvars == 0 ? 0.0 : (SCIP_Real)heurdata->intnblockvars/(heurdata->ntotalintvars) * 100.0, 1316 heurdata->intnblocks == 0 ? 0.0 : heurdata->intnblockvars/(SCIP_Real)(heurdata->intnblocks), 1317 heurdata->maxintblocksize); 1318 1319 SCIPstatisticMessage( 1320 " Twoopt results : " 1321 "%6d %6d %4d %4.2g (runs, binary exchanges, Integer shiftings, matching rate)\n", 1322 heurdata->nruns, 1323 heurdata->binnexchanges, 1324 heurdata->intnexchanges, 1325 heurdata->matchingrate); 1326 1327 /* set statistics to initial values*/ 1328 heurdata->binnblockvars = 0; 1329 heurdata->binnblocks = 0; 1330 heurdata->intnblocks = 0; 1331 heurdata->intnblockvars = 0; 1332 heurdata->binnexchanges = 0; 1333 heurdata->intnexchanges = 0; 1334 #endif 1335 1336 /* free the allocated memory for the binary variables */ 1337 if( heurdata->binvars != NULL ) 1338 { 1339 SCIPfreeBlockMemoryArray(scip, &heurdata->binvars, heurdata->nbinvars); 1340 } 1341 1342 if( heurdata->nbinblocks > 0 ) 1343 { 1344 assert(heurdata->binblockstart != NULL); 1345 assert(heurdata->binblockend != NULL); 1346 1347 SCIPfreeBlockMemoryArray(scip, &heurdata->binblockstart, heurdata->nbinblocks); 1348 SCIPfreeBlockMemoryArray(scip, &heurdata->binblockend, heurdata->nbinblocks); 1349 } 1350 heurdata->nbinvars = 0; 1351 heurdata->nbinblocks = 0; 1352 1353 if( heurdata->nintblocks > 0 ) 1354 { 1355 assert(heurdata->intblockstart != NULL); 1356 assert(heurdata->intblockend != NULL); 1357 1358 SCIPfreeBlockMemoryArray(scip, &heurdata->intblockstart, heurdata->nintblocks); 1359 SCIPfreeBlockMemoryArray(scip, &heurdata->intblockend, heurdata->nintblocks); 1360 } 1361 1362 /* free the allocated memory for the integers */ 1363 if( heurdata->intvars != NULL ) 1364 { 1365 SCIPfreeBlockMemoryArray(scip, &heurdata->intvars, heurdata->nintvars); 1366 } 1367 1368 heurdata->nbinblocks = 0; 1369 heurdata->nintblocks = 0; 1370 heurdata->nbinvars = 0; 1371 heurdata->nintvars = 0; 1372 1373 assert(heurdata->binvars == NULL); 1374 assert(heurdata->intvars == NULL); 1375 1376 /* free random number generator */ 1377 SCIPfreeRandom(scip, &heurdata->randnumgen); 1378 1379 SCIPheurSetData(heur, heurdata); 1380 1381 return SCIP_OKAY; 1382 } 1383 1384 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */ 1385 static 1386 SCIP_DECL_HEURINITSOL(heurInitsolTwoopt) 1387 { 1388 SCIP_HEURDATA* heurdata; 1389 assert(heur != NULL); 1390 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 1391 assert(scip != NULL); 1392 1393 /* get heuristic data */ 1394 heurdata = SCIPheurGetData(heur); 1395 1396 assert(heurdata != NULL); 1397 assert(heurdata->binvars == NULL && heurdata->intvars == NULL); 1398 assert(heurdata->binblockstart == NULL && heurdata->binblockend == NULL); 1399 assert(heurdata->intblockstart == NULL && heurdata->intblockend == NULL); 1400 1401 /* set heuristic data to initial values, but increase the total number of runs */ 1402 heurdata->nbinvars = 0; 1403 heurdata->nintvars = 0; 1404 heurdata->lastsolindex = -1; 1405 heurdata->presolved = FALSE; 1406 1407 #ifdef SCIP_STATISTIC 1408 ++(heurdata->nruns); 1409 #endif 1410 1411 SCIPheurSetData(heur, heurdata); 1412 1413 return SCIP_OKAY; 1414 } 1415 1416 1417 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */ 1418 static 1419 SCIP_DECL_HEUREXITSOL(heurExitsolTwoopt) 1420 { 1421 SCIP_HEURDATA* heurdata; 1422 int nbinvars; 1423 int nintvars; 1424 1425 assert(heur != NULL); 1426 assert(scip != NULL); 1427 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0); 1428 assert(scip != NULL); 1429 1430 /* get heuristic data */ 1431 heurdata = SCIPheurGetData(heur); 1432 1433 assert(heurdata != NULL); 1434 1435 nbinvars = heurdata->nbinvars; 1436 nintvars = heurdata->nintvars; 1437 1438 /* free the allocated memory for the binary variables */ 1439 if( heurdata->binvars != NULL ) 1440 { 1441 SCIPfreeBlockMemoryArray(scip, &heurdata->binvars, nbinvars); 1442 } 1443 if( heurdata->binblockstart != NULL ) 1444 { 1445 assert(heurdata->binblockend != NULL); 1446 1447 SCIPfreeBlockMemoryArray(scip, &heurdata->binblockstart, heurdata->nbinblocks); 1448 SCIPfreeBlockMemoryArray(scip, &heurdata->binblockend, heurdata->nbinblocks); 1449 } 1450 heurdata->nbinvars = 0; 1451 heurdata->nbinblocks = 0; 1452 1453 if( heurdata->intblockstart != NULL ) 1454 { 1455 assert(heurdata->intblockend != NULL); 1456 1457 SCIPfreeBlockMemoryArray(scip, &heurdata->intblockstart, heurdata->nintblocks); 1458 SCIPfreeBlockMemoryArray(scip, &heurdata->intblockend, heurdata->nintblocks); 1459 } 1460 heurdata->nintblocks = 0; 1461 1462 /* free the allocated memory for the integers */ 1463 if( heurdata->intvars != NULL ) 1464 { 1465 SCIPfreeBlockMemoryArray(scip, &heurdata->intvars, nintvars); 1466 } 1467 1468 heurdata->nintvars = 0; 1469 1470 assert(heurdata->binvars == NULL && heurdata->intvars == NULL); 1471 assert(heurdata->binblockstart == NULL && heurdata->binblockend == NULL); 1472 assert(heurdata->intblockstart == NULL && heurdata->intblockend == NULL); 1473 1474 /* set heuristic data */ 1475 SCIPheurSetData(heur, heurdata); 1476 1477 return SCIP_OKAY; 1478 } 1479 1480 /** execution method of primal heuristic */ 1481 static 1482 SCIP_DECL_HEUREXEC(heurExecTwoopt) 1483 { /*lint --e{715}*/ 1484 SCIP_HEURDATA* heurdata; 1485 SCIP_SOL* bestsol; 1486 SCIP_SOL* worksol; 1487 SCIP_ROW** lprows; 1488 SCIP_Real* activities; 1489 SCIP_COL** cols; 1490 int ncols; 1491 int nbinvars; 1492 int nintvars; 1493 int ndiscvars; 1494 int nlprows; 1495 int ncolsforsorting; 1496 SCIP_Bool improvement; 1497 SCIP_Bool presolthiscall; 1498 SCIP_Bool varboundserr; 1499 1500 assert(heur != NULL); 1501 assert(scip != NULL); 1502 assert(result != NULL); 1503 1504 /* get heuristic data */ 1505 heurdata = SCIPheurGetData(heur); 1506 assert(heurdata != NULL); 1507 1508 *result = SCIP_DIDNOTRUN; 1509 1510 /* we need an LP */ 1511 if( SCIPgetNLPRows(scip) == 0 ) 1512 return SCIP_OKAY; 1513 1514 bestsol = SCIPgetBestSol(scip); 1515 1516 /* ensure that heuristic has not already been processed on current incumbent */ 1517 if( bestsol == NULL || heurdata->lastsolindex == SCIPsolGetIndex(bestsol) ) 1518 return SCIP_OKAY; 1519 1520 heurdata->lastsolindex = SCIPsolGetIndex(bestsol); 1521 1522 /* we can only work on solutions valid in the transformed space */ 1523 if( SCIPsolIsOriginal(bestsol) ) 1524 return SCIP_OKAY; 1525 1526 #ifdef SCIP_DEBUG 1527 SCIP_CALL( SCIPprintSol(scip, bestsol, NULL, TRUE) ); 1528 #endif 1529 1530 /* ensure that the user defined number of nodes after last best solution has been reached, return otherwise */ 1531 if( (SCIPgetNNodes(scip) - SCIPsolGetNodenum(bestsol)) < heurdata->waitingnodes ) 1532 return SCIP_OKAY; 1533 1534 presolthiscall = FALSE; 1535 SCIP_CALL( SCIPgetLPColsData(scip,&cols, &ncols) ); 1536 ndiscvars = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip); 1537 ncolsforsorting = MIN(ncols, ndiscvars); 1538 1539 /* ensure that heuristic specific presolve is applied when heuristic is executed first */ 1540 if( !heurdata->presolved ) 1541 { 1542 SCIP_CALL( SCIPgetLPColsData(scip,&cols, &ncols) ); 1543 1544 for( int i = 0; i < ncolsforsorting; ++i ) 1545 SCIPcolSort(cols[i]); 1546 1547 SCIP_CALL( presolveTwoOpt(scip, heur, heurdata) ); 1548 presolthiscall = TRUE; 1549 } 1550 1551 assert(heurdata->presolved); 1552 1553 SCIPdebugMsg(scip, " Twoopt heuristic is %sexecuting.\n", heurdata->execute ? "" : "not "); 1554 /* ensure that presolve has detected structures in the problem to which the 2-optimization can be applied. 1555 * That is if variables exist which share a common set of LP-rows. */ 1556 if( !heurdata->execute ) 1557 return SCIP_OKAY; 1558 1559 nbinvars = heurdata->nbinvars; 1560 nintvars = heurdata->nintvars; 1561 ndiscvars = nbinvars + nintvars; 1562 1563 /* we need to be able to start diving from current node in order to resolve the LP 1564 * with continuous or implicit integer variables 1565 */ 1566 if( SCIPgetNVars(scip) > ndiscvars && ( !SCIPhasCurrentNodeLP(scip) || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL ) ) 1567 return SCIP_OKAY; 1568 1569 /* problem satisfies all necessary conditions for 2-optimization heuristic, execute heuristic! */ 1570 *result = SCIP_DIDNOTFIND; 1571 1572 /* initialize a working solution as a copy of the current incumbent to be able to store 1573 * possible improvements obtained by 2-optimization */ 1574 SCIP_CALL( SCIPcreateSolCopy(scip, &worksol, bestsol) ); 1575 SCIPsolSetHeur(worksol, heur); 1576 1577 /* get the LP row activities from current incumbent bestsol */ 1578 SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nlprows) ); 1579 SCIP_CALL( SCIPallocBufferArray(scip, &activities, nlprows) ); 1580 1581 for( int i = 0; i < nlprows; ++i ) 1582 { 1583 SCIP_ROW* row; 1584 1585 row = lprows[i]; 1586 assert(row != NULL); 1587 assert(SCIProwGetLPPos(row) == i); 1588 SCIPdebugMsg(scip, " Row <%d> is %sin LP: \n", i, SCIProwGetLPPos(row) >= 0 ? "" : "not "); 1589 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) ); 1590 activities[i] = SCIPgetRowSolActivity(scip, row, bestsol); 1591 1592 /* Heuristic does not provide infeasibility recovery, thus if any constraint is violated, 1593 * execution has to be terminated. 1594 */ 1595 if( !SCIProwIsLocal(row) && (SCIPisFeasLT(scip, activities[i], SCIProwGetLhs(row)) 1596 || SCIPisFeasGT(scip, activities[i], SCIProwGetRhs(row))) ) 1597 goto TERMINATE; 1598 } 1599 1600 if( !presolthiscall ) 1601 { 1602 for( int i = 0; i < ncolsforsorting; ++i ) 1603 SCIPcolSort(cols[i]); 1604 } 1605 SCIPdebugMsg(scip, " Twoopt heuristic has initialized activities and sorted rows! \n"); 1606 1607 /* start with binary optimization */ 1608 improvement = FALSE; 1609 varboundserr = FALSE; 1610 1611 if( heurdata->nbinblocks > 0 ) 1612 { 1613 SCIP_CALL( optimize(scip, worksol, heurdata->binvars, heurdata->binblockstart, heurdata->binblockend, heurdata->nbinblocks, 1614 OPTTYPE_BINARY, activities, nlprows, &improvement, &varboundserr, heurdata) ); 1615 1616 SCIPdebugMsg(scip, " Binary Optimization finished!\n"); 1617 } 1618 1619 if( varboundserr ) 1620 goto TERMINATE; 1621 1622 /* ensure that their are at least two integer variables which do not have the same coefficient 1623 * in the objective function. In one of these cases, the heuristic will automatically skip the 1624 * integer variable optimization */ 1625 if( heurdata->nintblocks > 0 ) 1626 { 1627 assert(heurdata->intopt); 1628 SCIP_CALL( optimize(scip, worksol, heurdata->intvars, heurdata->intblockstart, heurdata->intblockend, heurdata->nintblocks, 1629 OPTTYPE_INTEGER, activities, nlprows, &improvement, &varboundserr, heurdata) ); 1630 1631 SCIPdebugMsg(scip, " Integer Optimization finished!\n"); 1632 } 1633 1634 if( ! improvement || varboundserr ) 1635 goto TERMINATE; 1636 1637 if( SCIPgetNVars(scip) == ndiscvars ) 1638 { 1639 /* the problem is a pure IP, hence, no continuous or implicit variables are left for diving. 1640 * try if new working solution is feasible in original problem */ 1641 SCIP_Bool success; 1642 #ifndef NDEBUG 1643 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) ); 1644 #else 1645 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, TRUE, &success) ); 1646 #endif 1647 1648 if( success ) 1649 { 1650 SCIPdebugMsg(scip, "found feasible shifted solution:\n"); 1651 SCIPdebug( SCIP_CALL( SCIPprintSol(scip, worksol, NULL, FALSE) ) ); 1652 heurdata->lastsolindex = SCIPsolGetIndex(bestsol); 1653 *result = SCIP_FOUNDSOL; 1654 1655 #ifdef SCIP_STATISTIC 1656 SCIPstatisticMessage("***Twoopt improved solution found by %10s . \n", 1657 SCIPsolGetHeur(bestsol) != NULL ? SCIPheurGetName(SCIPsolGetHeur(bestsol)) :"Tree"); 1658 #endif 1659 } 1660 } 1661 /* fix the integer variables and start diving to optimize continuous variables with respect to reduced domain */ 1662 else 1663 { 1664 SCIP_VAR** allvars; 1665 SCIP_Bool lperror; 1666 #ifdef NDEBUG 1667 SCIP_RETCODE retstat; 1668 #endif 1669 1670 SCIPdebugMsg(scip, "shifted solution should be feasible -> solve LP to fix continuous variables to best values\n"); 1671 1672 allvars = SCIPgetVars(scip); 1673 1674 #ifdef SCIP_DEBUG 1675 for( int i = ndiscvars; i < SCIPgetNVars(scip); ++i ) 1676 { 1677 SCIPdebugMsg(scip, " Cont. variable <%s>, status %d with bounds [%g <= %g <= x <= %g <= %g]\n", 1678 SCIPvarGetName(allvars[i]), SCIPvarGetStatus(allvars[i]), SCIPvarGetLbGlobal(allvars[i]), SCIPvarGetLbLocal(allvars[i]), SCIPvarGetUbLocal(allvars[i]), 1679 SCIPvarGetUbGlobal(allvars[i])); 1680 } 1681 #endif 1682 1683 /* start diving to calculate the LP relaxation */ 1684 SCIP_CALL( SCIPstartDive(scip) ); 1685 1686 /* set the bounds of the variables: fixed for integers, global bounds for continuous */ 1687 for( int i = 0; i < SCIPgetNVars(scip); ++i ) 1688 { 1689 if( SCIPvarGetStatus(allvars[i]) == SCIP_VARSTATUS_COLUMN ) 1690 { 1691 SCIP_CALL( SCIPchgVarLbDive(scip, allvars[i], SCIPvarGetLbGlobal(allvars[i])) ); 1692 SCIP_CALL( SCIPchgVarUbDive(scip, allvars[i], SCIPvarGetUbGlobal(allvars[i])) ); 1693 } 1694 } 1695 1696 /* apply this after global bounds to not cause an error with intermediate empty domains */ 1697 for( int i = 0; i < ndiscvars; ++i ) 1698 { 1699 if( SCIPvarGetStatus(allvars[i]) == SCIP_VARSTATUS_COLUMN ) 1700 { 1701 SCIP_Real solval; 1702 1703 solval = SCIPgetSolVal(scip, worksol, allvars[i]); 1704 1705 assert(SCIPvarGetType(allvars[i]) != SCIP_VARTYPE_CONTINUOUS && SCIPisFeasIntegral(scip, solval)); 1706 1707 SCIP_CALL( SCIPchgVarLbDive(scip, allvars[i], solval) ); 1708 SCIP_CALL( SCIPchgVarUbDive(scip, allvars[i], solval) ); 1709 } 1710 } 1711 for( int i = 0; i < ndiscvars; ++i ) 1712 { 1713 assert( SCIPisFeasEQ(scip, SCIPgetVarLbDive(scip, allvars[i]), SCIPgetVarUbDive(scip, allvars[i])) ); 1714 } 1715 /* solve LP */ 1716 SCIPdebugMsg(scip, " -> old LP iterations: %" SCIP_LONGINT_FORMAT "\n", SCIPgetNLPIterations(scip)); 1717 1718 /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic. 1719 * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. */ 1720 #ifdef NDEBUG 1721 retstat = SCIPsolveDiveLP(scip, -1, &lperror, NULL); 1722 if( retstat != SCIP_OKAY ) 1723 { 1724 SCIPwarningMessage(scip, "Error while solving LP in Twoopt heuristic; LP solve terminated with code <%d>\n",retstat); 1725 } 1726 #else 1727 SCIP_CALL( SCIPsolveDiveLP(scip, -1, &lperror, NULL) ); 1728 #endif 1729 1730 SCIPdebugMsg(scip, " -> new LP iterations: %" SCIP_LONGINT_FORMAT "\n", SCIPgetNLPIterations(scip)); 1731 SCIPdebugMsg(scip, " -> error=%u, status=%d\n", lperror, SCIPgetLPSolstat(scip)); 1732 1733 /* check if this is a feasible solution */ 1734 if( !lperror && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL ) 1735 { 1736 SCIP_Bool success; 1737 1738 /* copy the current LP solution to the working solution */ 1739 SCIP_CALL( SCIPlinkLPSol(scip, worksol) ); 1740 1741 /* check solution for feasibility */ 1742 #ifndef NDEBUG 1743 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) ); 1744 #else 1745 SCIP_CALL( SCIPtrySol(scip, worksol, FALSE, FALSE, FALSE, FALSE, TRUE, &success) ); 1746 #endif 1747 1748 if( success ) 1749 { 1750 SCIPdebugMsg(scip, "found feasible shifted solution:\n"); 1751 SCIPdebug( SCIP_CALL( SCIPprintSol(scip, worksol, NULL, FALSE) ) ); 1752 heurdata->lastsolindex = SCIPsolGetIndex(bestsol); 1753 *result = SCIP_FOUNDSOL; 1754 1755 #ifdef SCIP_STATISTIC 1756 SCIPstatisticMessage("*** Twoopt improved solution found by %10s . \n", 1757 SCIPsolGetHeur(bestsol) != NULL ? SCIPheurGetName(SCIPsolGetHeur(bestsol)) :"Tree"); 1758 #endif 1759 } 1760 } 1761 1762 /* terminate the diving */ 1763 SCIP_CALL( SCIPendDive(scip) ); 1764 } 1765 1766 TERMINATE: 1767 SCIPdebugMsg(scip, "Termination of Twoopt heuristic\n"); 1768 SCIPfreeBufferArray(scip, &activities); 1769 SCIP_CALL( SCIPfreeSol(scip, &worksol) ); 1770 1771 return SCIP_OKAY; 1772 } 1773 1774 /* 1775 * primal heuristic specific interface methods 1776 */ 1777 1778 /** creates the twoopt primal heuristic and includes it in SCIP */ 1779 SCIP_RETCODE SCIPincludeHeurTwoopt( 1780 SCIP* scip /**< SCIP data structure */ 1781 ) 1782 { 1783 SCIP_HEURDATA* heurdata; 1784 SCIP_HEUR* heur; 1785 1786 /* create Twoopt primal heuristic data */ 1787 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) ); 1788 1789 /* include primal heuristic */ 1790 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, 1791 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS, 1792 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecTwoopt, heurdata) ); 1793 1794 assert(heur != NULL); 1795 1796 /* set non-NULL pointers to callback methods */ 1797 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyTwoopt) ); 1798 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeTwoopt) ); 1799 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitTwoopt) ); 1800 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitTwoopt) ); 1801 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolTwoopt) ); 1802 SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolTwoopt) ); 1803 1804 /* include boolean flag intopt */ 1805 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/twoopt/intopt", " Should Integer-2-Optimization be applied or not?", 1806 &heurdata->intopt, TRUE, DEFAULT_INTOPT, NULL, NULL) ); 1807 1808 /* include parameter waitingnodes */ 1809 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/twoopt/waitingnodes", "user parameter to determine number of " 1810 "nodes to wait after last best solution before calling heuristic", 1811 &heurdata->waitingnodes, TRUE, DEFAULT_WAITINGNODES, 0, 10000, NULL, NULL)); 1812 1813 /* include parameter maxnslaves */ 1814 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/twoopt/maxnslaves", "maximum number of slaves for one master variable", 1815 &heurdata->maxnslaves, TRUE, DEFAULT_MAXNSLAVES, -1, 1000000, NULL, NULL) ); 1816 1817 /* include parameter matchingrate */ 1818 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/twoopt/matchingrate", 1819 "parameter to determine the percentage of rows two variables have to share before they are considered equal", 1820 &heurdata->matchingrate, TRUE, DEFAULT_MATCHINGRATE, 0.0, 1.0, NULL, NULL) ); 1821 1822 return SCIP_OKAY; 1823 } 1824