Scippy

SCIP

Solving Constraint Integer Programs

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