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