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