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_shifting.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  LP rounding heuristic that tries to recover from intermediate infeasibilities and shifts continuous variables
28   	 * @author Tobias Achterberg
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include "blockmemshell/memory.h"
34   	#include "scip/heur_shifting.h"
35   	#include "scip/pub_heur.h"
36   	#include "scip/pub_lp.h"
37   	#include "scip/pub_message.h"
38   	#include "scip/pub_misc.h"
39   	#include "scip/pub_var.h"
40   	#include "scip/scip_branch.h"
41   	#include "scip/scip_heur.h"
42   	#include "scip/scip_lp.h"
43   	#include "scip/scip_mem.h"
44   	#include "scip/scip_message.h"
45   	#include "scip/scip_numerics.h"
46   	#include "scip/scip_prob.h"
47   	#include "scip/scip_randnumgen.h"
48   	#include "scip/scip_sol.h"
49   	#include "scip/scip_solvingstats.h"
50   	#include <string.h>
51   	
52   	#define HEUR_NAME             "shifting"
53   	#define HEUR_DESC             "LP rounding heuristic with infeasibility recovering also using continuous variables"
54   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_ROUNDING
55   	#define HEUR_PRIORITY         -5000
56   	#define HEUR_FREQ             10
57   	#define HEUR_FREQOFS          0
58   	#define HEUR_MAXDEPTH         -1
59   	#define HEUR_TIMING           SCIP_HEURTIMING_DURINGLPLOOP
60   	#define HEUR_USESSUBSCIP      FALSE  /**< does the heuristic use a secondary SCIP instance? */
61   	
62   	#define MAXSHIFTINGS          50     /**< maximal number of non improving shiftings */
63   	#define WEIGHTFACTOR          1.1
64   	#define DEFAULT_RANDSEED      31     /**< initial random seed */
65   	
66   	
67   	/* locally defined heuristic data */
68   	struct SCIP_HeurData
69   	{
70   	   SCIP_SOL*             sol;                /**< working solution */
71   	   SCIP_RANDNUMGEN*      randnumgen;         /**< random number generator */
72   	   SCIP_Longint          lastlp;             /**< last LP number where the heuristic was applied */
73   	};
74   	
75   	
76   	/*
77   	 * local methods
78   	 */
79   	
80   	/** update row violation arrays after a row's activity value changed */
81   	static
82   	void updateViolations(
83   	   SCIP*                 scip,               /**< SCIP data structure */
84   	   SCIP_ROW*             row,                /**< LP row */
85   	   SCIP_ROW**            violrows,           /**< array with currently violated rows */
86   	   int*                  violrowpos,         /**< position of LP rows in violrows array */
87   	   int*                  nviolrows,          /**< pointer to the number of currently violated rows */
88   	   int*                  nviolfracrows,      /**< pointer to the number of violated rows with fractional candidates */
89   	   int*                  nfracsinrow,        /**< array with number of fractional variables for every row */
90   	   SCIP_Real             oldactivity,        /**< old activity value of LP row */
91   	   SCIP_Real             newactivity         /**< new activity value of LP row */
92   	   )
93   	{
94   	   SCIP_Real lhs;
95   	   SCIP_Real rhs;
96   	   SCIP_Bool oldviol;
97   	   SCIP_Bool newviol;
98   	
99   	   assert(violrows != NULL);
100  	   assert(violrowpos != NULL);
101  	   assert(nviolrows != NULL);
102  	
103  	   lhs = SCIProwGetLhs(row);
104  	   rhs = SCIProwGetRhs(row);
105  	
106  	   /* SCIPisFeasLT cannot handle comparing different infinities. To prevent this, we make a case distinction. */
107  	   if( !(SCIPisInfinity(scip, oldactivity) || SCIPisInfinity(scip, -oldactivity)) )
108  	   {
109  	      oldviol = (SCIPisFeasLT(scip, oldactivity, lhs) || SCIPisFeasGT(scip, oldactivity, rhs));
110  	   }
111  	   else
112  	   {
113  	      oldviol = (SCIPisInfinity(scip, -oldactivity) && !SCIPisInfinity(scip, -lhs)) ||
114  	         (SCIPisInfinity(scip, oldactivity) && !SCIPisInfinity(scip, rhs));
115  	   }
116  	
117  	   /* SCIPisFeasLT cannot handle comparing different infinities. To prevent this, we make a case distinction. */
118  	   if( !(SCIPisInfinity(scip, newactivity) || SCIPisInfinity(scip, -newactivity)) )
119  	   {
120  	      newviol = (SCIPisFeasLT(scip, newactivity, lhs) || SCIPisFeasGT(scip, newactivity, rhs));
121  	   }
122  	   else
123  	   {
124  	      newviol = (SCIPisInfinity(scip, -newactivity) && !SCIPisInfinity(scip, -lhs)) ||
125  	         (SCIPisInfinity(scip, newactivity) && !SCIPisInfinity(scip, rhs));
126  	   }
127  	
128  	   if( oldviol != newviol )
129  	   {
130  	      int rowpos;
131  	      int violpos;
132  	
133  	      rowpos = SCIProwGetLPPos(row);
134  	      assert(rowpos >= 0);
135  	
136  	      if( oldviol )
137  	      {
138  	         /* the row violation was repaired: remove row from violrows array, decrease violation count */
139  	         violpos = violrowpos[rowpos];
140  	         assert(0 <= violpos && violpos < *nviolrows);
141  	         assert(violrows[violpos] == row);
142  	         violrowpos[rowpos] = -1;
143  	
144  	         /* first, move the row to the end of the subset of violated rows with fractional variables */
145  	         if( nfracsinrow[rowpos] > 0 )
146  	         {
147  	            assert(violpos < *nviolfracrows);
148  	
149  	            /* replace with last violated row containing fractional variables */
150  	            if( violpos != *nviolfracrows - 1 )
151  	            {
152  	               violrows[violpos] = violrows[*nviolfracrows - 1];
153  	               violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos;
154  	               violpos = *nviolfracrows - 1;
155  	            }
156  	            (*nviolfracrows)--;
157  	         }
158  	
159  	         assert(violpos >= *nviolfracrows);
160  	
161  	         /* swap row at the end of the violated array to the position of this row and decrease the counter */
162  	         if( violpos != *nviolrows - 1 )
163  	         {
164  	            violrows[violpos] = violrows[*nviolrows - 1];
165  	            violrowpos[SCIProwGetLPPos(violrows[violpos])] = violpos;
166  	         }
167  	         (*nviolrows)--;
168  	      }
169  	      else
170  	      {
171  	         /* the row is now violated: add row to violrows array, increase violation count */
172  	         assert(violrowpos[rowpos] == -1);
173  	         violrows[*nviolrows] = row;
174  	         violrowpos[rowpos] = *nviolrows;
175  	         (*nviolrows)++;
176  	
177  	         /* if the row contains fractional variables, swap with the last violated row that has no fractional variables
178  	          * at position *nviolfracrows
179  	          */
180  	         if( nfracsinrow[rowpos] > 0 )
181  	         {
182  	            if( *nviolfracrows < *nviolrows - 1 )
183  	            {
184  	               assert(nfracsinrow[SCIProwGetLPPos(violrows[*nviolfracrows])] == 0);
185  	
186  	               violrows[*nviolrows - 1] = violrows[*nviolfracrows];
187  	               violrowpos[SCIProwGetLPPos(violrows[*nviolrows - 1])] = *nviolrows - 1;
188  	
189  	               violrows[*nviolfracrows] = row;
190  	               violrowpos[rowpos] = *nviolfracrows;
191  	            }
192  	            (*nviolfracrows)++;
193  	         }
194  	      }
195  	   }
196  	}
197  	
198  	/** update row activities after a variable's solution value changed */
199  	static
200  	SCIP_RETCODE updateActivities(
201  	   SCIP*                 scip,               /**< SCIP data structure */
202  	   SCIP_Real*            activities,         /**< LP row activities */
203  	   SCIP_ROW**            violrows,           /**< array with currently violated rows */
204  	   int*                  violrowpos,         /**< position of LP rows in violrows array */
205  	   int*                  nviolrows,          /**< pointer to the number of currently violated rows */
206  	   int*                  nviolfracrows,      /**< pointer to the number of violated rows with fractional candidates */
207  	   int*                  nfracsinrow,        /**< array with number of fractional variables for every row */
208  	   int                   nlprows,            /**< number of rows in current LP */
209  	   SCIP_VAR*             var,                /**< variable that has been changed */
210  	   SCIP_Real             oldsolval,          /**< old solution value of variable */
211  	   SCIP_Real             newsolval           /**< new solution value of variable */
212  	   )
213  	{
214  	   SCIP_COL* col;
215  	   SCIP_ROW** colrows;
216  	   SCIP_Real* colvals;
217  	   SCIP_Real delta;
218  	   int ncolrows;
219  	   int r;
220  	
221  	   assert(activities != NULL);
222  	   assert(nviolrows != NULL);
223  	   assert(0 <= *nviolrows && *nviolrows <= nlprows);
224  	
225  	   delta = newsolval - oldsolval;
226  	   col = SCIPvarGetCol(var);
227  	   colrows = SCIPcolGetRows(col);
228  	   colvals = SCIPcolGetVals(col);
229  	   ncolrows = SCIPcolGetNLPNonz(col);
230  	   assert(ncolrows == 0 || (colrows != NULL && colvals != NULL));
231  	
232  	   for( r = 0; r < ncolrows; ++r )
233  	   {
234  	      SCIP_ROW* row;
235  	      int rowpos;
236  	
237  	      row = colrows[r];
238  	      rowpos = SCIProwGetLPPos(row);
239  	      assert(-1 <= rowpos && rowpos < nlprows);
240  	
241  	      if( rowpos >= 0 && !SCIProwIsLocal(row) )
242  	      {
243  	         SCIP_Real oldactivity;
244  	         SCIP_Real newactivity;
245  	
246  	         assert(SCIProwIsInLP(row));
247  	
248  	         /* update row activity */
249  	         oldactivity = activities[rowpos];
250  	         if( !SCIPisInfinity(scip, -oldactivity) && !SCIPisInfinity(scip, oldactivity) )
251  	         {
252  	            newactivity = oldactivity + delta * colvals[r];
253  	            if( SCIPisInfinity(scip, newactivity) )
254  	               newactivity = SCIPinfinity(scip);
255  	            else if( SCIPisInfinity(scip, -newactivity) )
256  	               newactivity = -SCIPinfinity(scip);
257  	            activities[rowpos] = newactivity;
258  	
259  	            /* update row violation arrays */
260  	            updateViolations(scip, row, violrows, violrowpos, nviolrows, nviolfracrows, nfracsinrow, oldactivity, newactivity);
261  	         }
262  	      }
263  	   }
264  	
265  	   return SCIP_OKAY;
266  	}
267  	
268  	/** returns a variable, that pushes activity of the row in the given direction with minimal negative impact on other rows;
269  	 *  if variables have equal impact, chooses the one with best objective value improvement in corresponding direction;
270  	 *  prefer fractional integers over other variables in order to become integral during the process;
271  	 *  shifting in a direction is forbidden, if this forces the objective value over the upper bound, or if the variable
272  	 *  was already shifted in the opposite direction
273  	 */
274  	static
275  	SCIP_RETCODE selectShifting(
276  	   SCIP*                 scip,               /**< SCIP data structure */
277  	   SCIP_SOL*             sol,                /**< primal solution */
278  	   SCIP_ROW*             row,                /**< LP row */
279  	   SCIP_Real             rowactivity,        /**< activity of LP row */
280  	   int                   direction,          /**< should the activity be increased (+1) or decreased (-1)? */
281  	   SCIP_Real*            nincreases,         /**< array with weighted number of increasings per variables */
282  	   SCIP_Real*            ndecreases,         /**< array with weighted number of decreasings per variables */
283  	   SCIP_Real             increaseweight,     /**< current weight of increase/decrease updates */
284  	   SCIP_VAR**            shiftvar,           /**< pointer to store the shifting variable, returns NULL if impossible */
285  	   SCIP_Real*            oldsolval,          /**< pointer to store old solution value of shifting variable */
286  	   SCIP_Real*            newsolval           /**< pointer to store new (shifted) solution value of shifting variable */
287  	   )
288  	{
289  	   SCIP_COL** rowcols;
290  	   SCIP_Real* rowvals;
291  	   int nrowcols;
292  	   SCIP_Real activitydelta;
293  	   SCIP_Real bestshiftscore;
294  	   SCIP_Real bestdeltaobj;
295  	   int c;
296  	
297  	   assert(direction == +1 || direction == -1);
298  	   assert(nincreases != NULL);
299  	   assert(ndecreases != NULL);
300  	   assert(shiftvar != NULL);
301  	   assert(oldsolval != NULL);
302  	   assert(newsolval != NULL);
303  	
304  	   /* get row entries */
305  	   rowcols = SCIProwGetCols(row);
306  	   rowvals = SCIProwGetVals(row);
307  	   nrowcols = SCIProwGetNLPNonz(row);
308  	
309  	   /* calculate how much the activity must be shifted in order to become feasible */
310  	   activitydelta = (direction == +1 ? SCIProwGetLhs(row) - rowactivity : SCIProwGetRhs(row) - rowactivity);
311  	   assert((direction == +1 && SCIPisPositive(scip, activitydelta))
312  	      || (direction == -1 && SCIPisNegative(scip, activitydelta)));
313  	
314  	   /* select shifting variable */
315  	   bestshiftscore = SCIP_REAL_MAX;
316  	   bestdeltaobj = SCIPinfinity(scip);
317  	   *shiftvar = NULL;
318  	   *newsolval = 0.0;
319  	   *oldsolval = 0.0;
320  	   for( c = 0; c < nrowcols; ++c )
321  	   {
322  	      SCIP_COL* col;
323  	      SCIP_VAR* var;
324  	      SCIP_Real val;
325  	      SCIP_Real solval;
326  	      SCIP_Real shiftval;
327  	      SCIP_Real shiftscore;
328  	      SCIP_Bool isinteger;
329  	      SCIP_Bool isfrac;
330  	      SCIP_Bool increase;
331  	
332  	      col = rowcols[c];
333  	      var = SCIPcolGetVar(col);
334  	      val = rowvals[c];
335  	      assert(!SCIPisZero(scip, val));
336  	      solval = SCIPgetSolVal(scip, sol, var);
337  	
338  	      isinteger = (SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
339  	      isfrac = (isinteger && !SCIPisFeasIntegral(scip, solval));
340  	      increase = (direction * val > 0.0);
341  	
342  	      /* calculate the score of the shifting (prefer smaller values) */
343  	      if( isfrac )
344  	         shiftscore = increase ? -1.0 / (SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL) + 1.0) :
345  	            -1.0 / (SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) + 1.0);
346  	      else
347  	      {
348  	         int probindex;
349  	         probindex = SCIPvarGetProbindex(var);
350  	
351  	         if( increase )
352  	            shiftscore = ndecreases[probindex]/increaseweight;
353  	         else
354  	            shiftscore = nincreases[probindex]/increaseweight;
355  	         if( isinteger )
356  	            shiftscore += 1.0;
357  	      }
358  	
359  	      if( shiftscore <= bestshiftscore )
360  	      {
361  	         SCIP_Real deltaobj;
362  	
363  	         if( !increase )
364  	         {
365  	            /* shifting down */
366  	            assert(direction * val < 0.0);
367  	            if( isfrac )
368  	               shiftval = SCIPfeasFloor(scip, solval);
369  	            else
370  	            {
371  	               SCIP_Real lb;
372  	
373  	               assert(activitydelta/val < 0.0);
374  	               shiftval = solval + activitydelta/val;
375  	               assert(shiftval <= solval); /* may be equal due to numerical digit erasement in the subtraction */
376  	               if( SCIPvarIsIntegral(var) )
377  	                  shiftval = SCIPfeasFloor(scip, shiftval);
378  	               lb = SCIPvarGetLbGlobal(var);
379  	               shiftval = MAX(shiftval, lb);
380  	            }
381  	         }
382  	         else
383  	         {
384  	            /* shifting up */
385  	            assert(direction * val > 0.0);
386  	            if( isfrac )
387  	               shiftval = SCIPfeasCeil(scip, solval);
388  	            else
389  	            {
390  	               SCIP_Real ub;
391  	
392  	               assert(activitydelta/val > 0.0);
393  	               shiftval = solval + activitydelta/val;
394  	               assert(shiftval >= solval); /* may be equal due to numerical digit erasement in the subtraction */
395  	               if( SCIPvarIsIntegral(var) )
396  	                  shiftval = SCIPfeasCeil(scip, shiftval);
397  	               ub = SCIPvarGetUbGlobal(var);
398  	               shiftval = MIN(shiftval, ub);
399  	            }
400  	         }
401  	
402  	         if( (SCIPisInfinity(scip, shiftval) && SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)))
403  	            || (SCIPisInfinity(scip, -shiftval) && SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)))
404  	            || SCIPisEQ(scip, shiftval, solval) )
405  	            continue;
406  	
407  	         deltaobj = SCIPvarGetObj(var) * (shiftval - solval);
408  	         if( shiftscore < bestshiftscore || deltaobj < bestdeltaobj )
409  	         {
410  	            bestshiftscore = shiftscore;
411  	            bestdeltaobj = deltaobj;
412  	            *shiftvar = var;
413  	            *oldsolval = solval;
414  	            *newsolval = shiftval;
415  	         }
416  	      }
417  	   }
418  	
419  	   return SCIP_OKAY;
420  	}
421  	
422  	/** returns a fractional variable, that has most impact on rows in opposite direction, i.e. that is most crucial to
423  	 *  fix in the other direction;
424  	 *  if variables have equal impact, chooses the one with best objective value improvement in corresponding direction;
425  	 *  shifting in a direction is forbidden, if this forces the objective value over the upper bound
426  	 */
427  	static
428  	SCIP_RETCODE selectEssentialRounding(
429  	   SCIP*                 scip,               /**< SCIP data structure */
430  	   SCIP_SOL*             sol,                /**< primal solution */
431  	   SCIP_Real             minobj,             /**< minimal objective value possible after shifting remaining fractional vars */
432  	   SCIP_VAR**            lpcands,            /**< fractional variables in LP */
433  	   int                   nlpcands,           /**< number of fractional variables in LP */
434  	   SCIP_VAR**            shiftvar,           /**< pointer to store the shifting variable, returns NULL if impossible */
435  	   SCIP_Real*            oldsolval,          /**< old (fractional) solution value of shifting variable */
436  	   SCIP_Real*            newsolval           /**< new (shifted) solution value of shifting variable */
437  	   )
438  	{
439  	   SCIP_VAR* var;
440  	   SCIP_Real solval;
441  	   SCIP_Real shiftval;
442  	   SCIP_Real obj;
443  	   SCIP_Real deltaobj;
444  	   SCIP_Real bestdeltaobj;
445  	   int maxnlocks;
446  	   int nlocks;
447  	   int v;
448  	
449  	   assert(shiftvar != NULL);
450  	   assert(oldsolval != NULL);
451  	   assert(newsolval != NULL);
452  	
453  	   /* select shifting variable */
454  	   maxnlocks = -1;
455  	   bestdeltaobj = SCIPinfinity(scip);
456  	   *shiftvar = NULL;
457  	   for( v = 0; v < nlpcands; ++v )
458  	   {
459  	      var = lpcands[v];
460  	      assert(SCIPvarGetType(var) == SCIP_VARTYPE_BINARY || SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
461  	
462  	      solval = SCIPgetSolVal(scip, sol, var);
463  	      if( !SCIPisFeasIntegral(scip, solval) )
464  	      {
465  	         obj = SCIPvarGetObj(var);
466  	
467  	         /* shifting down */
468  	         nlocks = SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL);
469  	         if( nlocks >= maxnlocks )
470  	         {
471  	            shiftval = SCIPfeasFloor(scip, solval);
472  	            deltaobj = obj * (shiftval - solval);
473  	            if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj - obj < SCIPgetCutoffbound(scip) )
474  	            {
475  	               maxnlocks = nlocks;
476  	               bestdeltaobj = deltaobj;
477  	               *shiftvar = var;
478  	               *oldsolval = solval;
479  	               *newsolval = shiftval;
480  	            }
481  	         }
482  	
483  	         /* shifting up */
484  	         nlocks = SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL);
485  	         if( nlocks >= maxnlocks )
486  	         {
487  	            shiftval = SCIPfeasCeil(scip, solval);
488  	            deltaobj = obj * (shiftval - solval);
489  	            if( (nlocks > maxnlocks || deltaobj < bestdeltaobj) && minobj + obj < SCIPgetCutoffbound(scip) )
490  	            {
491  	               maxnlocks = nlocks;
492  	               bestdeltaobj = deltaobj;
493  	               *shiftvar = var;
494  	               *oldsolval = solval;
495  	               *newsolval = shiftval;
496  	            }
497  	         }
498  	      }
499  	   }
500  	
501  	   return SCIP_OKAY;
502  	}
503  	
504  	/** adds a given value to the fractionality counters of the rows in which the given variable appears */
505  	static
506  	void addFracCounter(
507  	   int*                  nfracsinrow,        /**< array to store number of fractional variables per row */
508  	   SCIP_ROW**            violrows,           /**< array with currently violated rows */
509  	   int*                  violrowpos,         /**< position of LP rows in violrows array */
510  	   int*                  nviolfracrows,      /**< pointer to store the number of violated rows with fractional variables */
511  	   int                   nviolrows,          /**< the number of currently violated rows (stays unchanged in this method) */
512  	   int                   nlprows,            /**< number of rows in LP */
513  	   SCIP_VAR*             var,                /**< variable for which the counting should be updated */
514  	   int                   incval              /**< value that should be added to the corresponding array entries */
515  	   )
516  	{
517  	   SCIP_COL* col;
518  	   SCIP_ROW** rows;
519  	   int nrows;
520  	   int r;
521  	
522  	   assert(incval != 0);
523  	   assert(nviolrows >= *nviolfracrows);
524  	   col = SCIPvarGetCol(var);
525  	   rows = SCIPcolGetRows(col);
526  	   nrows = SCIPcolGetNLPNonz(col);
527  	   for( r = 0; r < nrows; ++r )
528  	   {
529  	      int rowlppos;
530  	      int theviolrowpos;
531  	
532  	      rowlppos = SCIProwGetLPPos(rows[r]);
533  	      assert(0 <= rowlppos && rowlppos < nlprows);
534  	      assert(!SCIProwIsLocal(rows[r]) || violrowpos[rowlppos] == -1);
535  	
536  	      /* skip local rows */
537  	      if( SCIProwIsLocal(rows[r]) )
538  	         continue;
539  	
540  	      nfracsinrow[rowlppos] += incval;
541  	      assert(nfracsinrow[rowlppos] >= 0);
542  	      theviolrowpos = violrowpos[rowlppos];
543  	
544  	      /* swap positions in violrows array if fractionality has changed to 0 */
545  	      if( theviolrowpos >= 0 )
546  	      {
547  	         /* first case: the number of fractional variables has become zero: swap row in violrows array to the
548  	          * second part, containing only violated rows without fractional variables
549  	          */
550  	         if( nfracsinrow[rowlppos] == 0 )
551  	         {
552  	            assert(theviolrowpos <= *nviolfracrows - 1);
553  	
554  	            /* nothing to do if row is already at the end of the first part, otherwise, swap it to the last position
555  	             * and decrease the counter */
556  	            if( theviolrowpos < *nviolfracrows - 1 )
557  	            {
558  	               violrows[theviolrowpos] = violrows[*nviolfracrows - 1];
559  	               violrows[*nviolfracrows - 1] = rows[r];
560  	
561  	               violrowpos[SCIProwGetLPPos(violrows[theviolrowpos])] = theviolrowpos;
562  	               violrowpos[rowlppos] = *nviolfracrows - 1;
563  	            }
564  	            (*nviolfracrows)--;
565  	         }
566  	         /* second interesting case: the number of fractional variables was zero before, swap it with the first
567  	          * violated row without fractional variables
568  	          */
569  	         else if( nfracsinrow[rowlppos] == incval )
570  	         {
571  	            assert(theviolrowpos >= *nviolfracrows);
572  	
573  	            /* nothing to do if the row is exactly located at index *nviolfracrows */
574  	            if( theviolrowpos > *nviolfracrows )
575  	            {
576  	               violrows[theviolrowpos] = violrows[*nviolfracrows];
577  	               violrows[*nviolfracrows] = rows[r];
578  	
579  	               violrowpos[SCIProwGetLPPos(violrows[theviolrowpos])] = theviolrowpos;
580  	               violrowpos[rowlppos] = *nviolfracrows;
581  	            }
582  	            (*nviolfracrows)++;
583  	         }
584  	      }
585  	   }
586  	}
587  	
588  	
589  	/*
590  	 * Callback methods
591  	 */
592  	
593  	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
594  	static
595  	SCIP_DECL_HEURCOPY(heurCopyShifting)
596  	{  /*lint --e{715}*/
597  	   assert(scip != NULL);
598  	   assert(heur != NULL);
599  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
600  	
601  	   /* call inclusion method of primal heuristic */
602  	   SCIP_CALL( SCIPincludeHeurShifting(scip) );
603  	
604  	   return SCIP_OKAY;
605  	}
606  	
607  	
608  	/** initialization method of primal heuristic (called after problem was transformed) */
609  	static
610  	SCIP_DECL_HEURINIT(heurInitShifting) /*lint --e{715}*/
611  	{  /*lint --e{715}*/
612  	   SCIP_HEURDATA* heurdata;
613  	
614  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
615  	   assert(SCIPheurGetData(heur) == NULL);
616  	
617  	   /* create heuristic data */
618  	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
619  	   SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
620  	   heurdata->lastlp = -1;
621  	
622  	   /* create random number generator */
623  	   SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen,
624  	         DEFAULT_RANDSEED, TRUE) );
625  	
626  	   SCIPheurSetData(heur, heurdata);
627  	
628  	   return SCIP_OKAY;
629  	}
630  	
631  	/** deinitialization method of primal heuristic (called before transformed problem is freed) */
632  	static
633  	SCIP_DECL_HEUREXIT(heurExitShifting) /*lint --e{715}*/
634  	{  /*lint --e{715}*/
635  	   SCIP_HEURDATA* heurdata;
636  	
637  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
638  	
639  	   /* free heuristic data */
640  	   heurdata = SCIPheurGetData(heur);
641  	   assert(heurdata != NULL);
642  	   SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
643  	
644  	   /* free random number generator */
645  	   SCIPfreeRandom(scip, &heurdata->randnumgen);
646  	
647  	   SCIPfreeBlockMemory(scip, &heurdata);
648  	   SCIPheurSetData(heur, NULL);
649  	
650  	   return SCIP_OKAY;
651  	}
652  	
653  	/** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
654  	static
655  	SCIP_DECL_HEURINITSOL(heurInitsolShifting)
656  	{
657  	   SCIP_HEURDATA* heurdata;
658  	
659  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
660  	
661  	   heurdata = SCIPheurGetData(heur);
662  	   assert(heurdata != NULL);
663  	   heurdata->lastlp = -1;
664  	
665  	   return SCIP_OKAY;
666  	}
667  	
668  	
669  	/** execution method of primal heuristic */
670  	static
671  	SCIP_DECL_HEUREXEC(heurExecShifting) /*lint --e{715}*/
672  	{  /*lint --e{715}*/
673  	   SCIP_HEURDATA* heurdata;
674  	   SCIP_SOL* sol;
675  	   SCIP_VAR** lpcands;
676  	   SCIP_Real* lpcandssol;
677  	   SCIP_ROW** lprows;
678  	   SCIP_Real* activities;
679  	   SCIP_ROW** violrows;
680  	   SCIP_Real* nincreases;
681  	   SCIP_Real* ndecreases;
682  	   int* violrowpos;
683  	   int* nfracsinrow;
684  	   SCIP_Real increaseweight;
685  	   SCIP_Real obj;
686  	   SCIP_Real bestshiftval;
687  	   SCIP_Real minobj;
688  	   int nlpcands;
689  	   int nlprows;
690  	   int nvars;
691  	   int nfrac;
692  	   int nviolrows;
693  	   int nviolfracrows;
694  	   int nprevviolrows;
695  	   int minnviolrows;
696  	   int nnonimprovingshifts;
697  	   int c;
698  	   int r;
699  	   SCIP_Longint nlps;
700  	   SCIP_Longint ncalls;
701  	   SCIP_Longint nsolsfound;
702  	   SCIP_Longint nnodes;
703  	
704  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
705  	   assert(scip != NULL);
706  	   assert(result != NULL);
707  	   assert(SCIPhasCurrentNodeLP(scip));
708  	
709  	   *result = SCIP_DIDNOTRUN;
710  	
711  	   /* only call heuristic, if an optimal LP solution is at hand */
712  	   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
713  	      return SCIP_OKAY;
714  	
715  	   /* only call heuristic, if the LP objective value is smaller than the cutoff bound */
716  	   if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) )
717  	      return SCIP_OKAY;
718  	
719  	   /* get heuristic data */
720  	   heurdata = SCIPheurGetData(heur);
721  	   assert(heurdata != NULL);
722  	
723  	   /* don't call heuristic, if we have already processed the current LP solution */
724  	   nlps = SCIPgetNLPs(scip);
725  	   if( nlps == heurdata->lastlp )
726  	      return SCIP_OKAY;
727  	   heurdata->lastlp = nlps;
728  	
729  	   /* don't call heuristic, if it was not successful enough in the past */
730  	   ncalls = SCIPheurGetNCalls(heur);
731  	   nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + SCIPheurGetNSolsFound(heur);
732  	   nnodes = SCIPgetNNodes(scip);
733  	   if( nnodes % ((ncalls/100)/(nsolsfound+1)+1) != 0 )
734  	      return SCIP_OKAY;
735  	
736  	   /* get fractional variables, that should be integral */
737  	   /* todo check if heuristic should include implicit integer variables for its calculations */
738  	   SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, NULL) );
739  	   nfrac = nlpcands;
740  	
741  	   /* only call heuristic, if LP solution is fractional */
742  	   if( nfrac == 0 )
743  	      return SCIP_OKAY;
744  	
745  	   *result = SCIP_DIDNOTFIND;
746  	
747  	   /* get LP rows */
748  	   SCIP_CALL( SCIPgetLPRowsData(scip, &lprows, &nlprows) );
749  	
750  	   SCIPdebugMsg(scip, "executing shifting heuristic: %d LP rows, %d fractionals\n", nlprows, nfrac);
751  	
752  	   /* get memory for activities, violated rows, and row violation positions */
753  	   nvars = SCIPgetNVars(scip);
754  	   SCIP_CALL( SCIPallocBufferArray(scip, &activities, nlprows) );
755  	   SCIP_CALL( SCIPallocBufferArray(scip, &violrows, nlprows) );
756  	   SCIP_CALL( SCIPallocBufferArray(scip, &violrowpos, nlprows) );
757  	   SCIP_CALL( SCIPallocBufferArray(scip, &nfracsinrow, nlprows) );
758  	   SCIP_CALL( SCIPallocBufferArray(scip, &nincreases, nvars) );
759  	   SCIP_CALL( SCIPallocBufferArray(scip, &ndecreases, nvars) );
760  	   BMSclearMemoryArray(nfracsinrow, nlprows);
761  	   BMSclearMemoryArray(nincreases, nvars);
762  	   BMSclearMemoryArray(ndecreases, nvars);
763  	
764  	   /* get the activities for all globally valid rows;
765  	    * the rows should be feasible, but due to numerical inaccuracies in the LP solver, they can be violated
766  	    */
767  	   nviolrows = 0;
768  	   for( r = 0; r < nlprows; ++r )
769  	   {
770  	      SCIP_ROW* row;
771  	
772  	      row = lprows[r];
773  	      assert(SCIProwGetLPPos(row) == r);
774  	
775  	      if( !SCIProwIsLocal(row) )
776  	      {
777  	         activities[r] = SCIPgetRowActivity(scip, row);
778  	         if( SCIPisFeasLT(scip, activities[r], SCIProwGetLhs(row))
779  	            || SCIPisFeasGT(scip, activities[r], SCIProwGetRhs(row)) )
780  	         {
781  	            violrows[nviolrows] = row;
782  	            violrowpos[r] = nviolrows;
783  	            nviolrows++;
784  	         }
785  	         else
786  	            violrowpos[r] = -1;
787  	      }
788  	      else
789  	         violrowpos[r] = -1;
790  	   }
791  	
792  	   nviolfracrows = 0;
793  	   /* calc the current number of fractional variables in rows */
794  	   for( c = 0; c < nlpcands; ++c )
795  	      addFracCounter(nfracsinrow, violrows, violrowpos, &nviolfracrows, nviolrows, nlprows, lpcands[c], +1);
796  	
797  	   /* get the working solution from heuristic's local data */
798  	   sol = heurdata->sol;
799  	   assert(sol != NULL);
800  	
801  	   /* copy the current LP solution to the working solution */
802  	   SCIP_CALL( SCIPlinkLPSol(scip, sol) );
803  	
804  	   /* calculate the minimal objective value possible after rounding fractional variables */
805  	   minobj = SCIPgetSolTransObj(scip, sol);
806  	   assert(minobj < SCIPgetCutoffbound(scip));
807  	   for( c = 0; c < nlpcands; ++c )
808  	   {
809  	      obj = SCIPvarGetObj(lpcands[c]);
810  	      bestshiftval = obj > 0.0 ? SCIPfeasFloor(scip, lpcandssol[c]) : SCIPfeasCeil(scip, lpcandssol[c]);
811  	      minobj += obj * (bestshiftval - lpcandssol[c]);
812  	   }
813  	
814  	   /* try to shift remaining variables in order to become/stay feasible */
815  	   nnonimprovingshifts = 0;
816  	   minnviolrows = INT_MAX;
817  	   increaseweight = 1.0;
818  	   while( (nfrac > 0 || nviolrows > 0) && nnonimprovingshifts < MAXSHIFTINGS )
819  	   {
820  	      SCIP_VAR* shiftvar;
821  	      SCIP_Real oldsolval;
822  	      SCIP_Real newsolval;
823  	      SCIP_Bool oldsolvalisfrac;
824  	      int probindex;
825  	
826  	      SCIPdebugMsg(scip, "shifting heuristic: nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g), cutoff=%g\n",
827  	         nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj),
828  	         SCIPretransformObj(scip, SCIPgetCutoffbound(scip)));
829  	
830  	      nprevviolrows = nviolrows;
831  	
832  	      /* choose next variable to process:
833  	       *  - if a violated row exists, shift a variable decreasing the violation, that has least impact on other rows
834  	       *  - otherwise, shift a variable, that has strongest devastating impact on rows in opposite direction
835  	       */
836  	      shiftvar = NULL;
837  	      oldsolval = 0.0;
838  	      newsolval = 0.0;
839  	      if( nviolrows > 0 && (nfrac == 0 || nnonimprovingshifts < MAXSHIFTINGS-1) )
840  	      {
841  	         SCIP_ROW* row;
842  	         int rowidx;
843  	         int rowpos;
844  	         int direction;
845  	
846  	         assert(nviolfracrows == 0 || nfrac > 0);
847  	         /* violated rows containing fractional variables are preferred; if such a row exists, choose the last one from the list
848  	          * (at position nviolfracrows - 1) because removing this row will cause one swapping operation less than other rows
849  	          */
850  	         if( nviolfracrows > 0 )
851  	            rowidx =  nviolfracrows - 1;
852  	         else
853  	            /* there is no violated row containing a fractional variable, select a violated row uniformly at random */
854  	            rowidx = SCIPrandomGetInt(heurdata->randnumgen, 0, nviolrows-1);
855  	
856  	         assert(0 <= rowidx && rowidx < nviolrows);
857  	         row = violrows[rowidx];
858  	         rowpos = SCIProwGetLPPos(row);
859  	         assert(0 <= rowpos && rowpos < nlprows);
860  	         assert(violrowpos[rowpos] == rowidx);
861  	         assert(nfracsinrow[rowpos] == 0 || rowidx == nviolfracrows - 1);
862  	
863  	         SCIPdebugMsg(scip, "shifting heuristic: try to fix violated row <%s>: %g <= %g <= %g\n",
864  	            SCIProwGetName(row), SCIProwGetLhs(row), activities[rowpos], SCIProwGetRhs(row));
865  	         SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
866  	
867  	         /* get direction in which activity must be shifted */
868  	         assert(SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row))
869  	            || SCIPisFeasGT(scip, activities[rowpos], SCIProwGetRhs(row)));
870  	         direction = SCIPisFeasLT(scip, activities[rowpos], SCIProwGetLhs(row)) ? +1 : -1;
871  	
872  	         /* search a variable that can shift the activity in the necessary direction */
873  	         SCIP_CALL( selectShifting(scip, sol, row, activities[rowpos], direction,
874  	               nincreases, ndecreases, increaseweight, &shiftvar, &oldsolval, &newsolval) );
875  	      }
876  	
877  	      if( shiftvar == NULL && nfrac > 0 )
878  	      {
879  	         SCIPdebugMsg(scip, "shifting heuristic: search rounding variable and try to stay feasible\n");
880  	         SCIP_CALL( selectEssentialRounding(scip, sol, minobj, lpcands, nlpcands, &shiftvar, &oldsolval, &newsolval) );
881  	      }
882  	
883  	      /* check, whether shifting was possible */
884  	      if( shiftvar == NULL || SCIPisEQ(scip, oldsolval, newsolval) )
885  	      {
886  	         SCIPdebugMsg(scip, "shifting heuristic:  -> didn't find a shifting variable\n");
887  	         break;
888  	      }
889  	
890  	      SCIPdebugMsg(scip, "shifting heuristic:  -> shift var <%s>[%g,%g], type=%d, oldval=%g, newval=%g, obj=%g\n",
891  	         SCIPvarGetName(shiftvar), SCIPvarGetLbGlobal(shiftvar), SCIPvarGetUbGlobal(shiftvar), SCIPvarGetType(shiftvar),
892  	         oldsolval, newsolval, SCIPvarGetObj(shiftvar));
893  	
894  	      /* update row activities of globally valid rows */
895  	      SCIP_CALL( updateActivities(scip, activities, violrows, violrowpos, &nviolrows, &nviolfracrows, nfracsinrow, nlprows,
896  	            shiftvar, oldsolval, newsolval) );
897  	      if( nviolrows >= nprevviolrows )
898  	         nnonimprovingshifts++;
899  	      else if( nviolrows < minnviolrows )
900  	      {
901  	         minnviolrows = nviolrows;
902  	         nnonimprovingshifts = 0;
903  	      }
904  	
905  	      /* store new solution value and decrease fractionality counter */
906  	      SCIP_CALL( SCIPsetSolVal(scip, sol, shiftvar, newsolval) );
907  	
908  	      /* update fractionality counter and minimal objective value possible after shifting remaining variables */
909  	      oldsolvalisfrac = !SCIPisFeasIntegral(scip, oldsolval)
910  	         && (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER);
911  	      obj = SCIPvarGetObj(shiftvar);
912  	      if( (SCIPvarGetType(shiftvar) == SCIP_VARTYPE_BINARY || SCIPvarGetType(shiftvar) == SCIP_VARTYPE_INTEGER)
913  	         && oldsolvalisfrac )
914  	      {
915  	         assert(SCIPisFeasIntegral(scip, newsolval));
916  	         nfrac--;
917  	         nnonimprovingshifts = 0;
918  	         minnviolrows = INT_MAX;
919  	         addFracCounter(nfracsinrow, violrows, violrowpos, &nviolfracrows, nviolrows, nlprows, shiftvar, -1);
920  	
921  	         /* the rounding was already calculated into the minobj -> update only if rounding in "wrong" direction */
922  	         if( obj > 0.0 && newsolval > oldsolval )
923  	            minobj += obj;
924  	         else if( obj < 0.0 && newsolval < oldsolval )
925  	            minobj -= obj;
926  	      }
927  	      else
928  	      {
929  	         /* update minimal possible objective value */
930  	         minobj += obj * (newsolval - oldsolval);
931  	      }
932  	
933  	      /* update increase/decrease arrays */
934  	      if( !oldsolvalisfrac )
935  	      {
936  	         probindex = SCIPvarGetProbindex(shiftvar);
937  	         assert(0 <= probindex && probindex < nvars);
938  	         increaseweight *= WEIGHTFACTOR;
939  	         if( newsolval < oldsolval )
940  	            ndecreases[probindex] += increaseweight;
941  	         else
942  	            nincreases[probindex] += increaseweight;
943  	         if( increaseweight >= 1e+09 )
944  	         {
945  	            int i;
946  	
947  	            for( i = 0; i < nvars; ++i )
948  	            {
949  	               nincreases[i] /= increaseweight;
950  	               ndecreases[i] /= increaseweight;
951  	            }
952  	            increaseweight = 1.0;
953  	         }
954  	      }
955  	
956  	      SCIPdebugMsg(scip, "shifting heuristic:  -> nfrac=%d, nviolrows=%d, obj=%g (best possible obj: %g)\n",
957  	         nfrac, nviolrows, SCIPgetSolOrigObj(scip, sol), SCIPretransformObj(scip, minobj));
958  	   }
959  	
960  	   /* check, if the new solution is feasible */
961  	   if( nfrac == 0 && nviolrows == 0 )
962  	   {
963  	      SCIP_Bool stored;
964  	
965  	      /* check solution for feasibility, and add it to solution store if possible
966  	       * neither integrality nor feasibility of LP rows has to be checked, because this is already
967  	       * done in the shifting heuristic itself; however, we better check feasibility of LP rows,
968  	       * because of numerical problems with activity updating
969  	       */
970  	      SCIP_CALL( SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, FALSE, TRUE, &stored) );
971  	
972  	      if( stored )
973  	      {
974  	         SCIPdebugMsg(scip, "found feasible shifted solution:\n");
975  	         SCIPdebug( SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) ) );
976  	         *result = SCIP_FOUNDSOL;
977  	      }
978  	   }
979  	
980  	   /* free memory buffers */
981  	   SCIPfreeBufferArray(scip, &ndecreases);
982  	   SCIPfreeBufferArray(scip, &nincreases);
983  	   SCIPfreeBufferArray(scip, &nfracsinrow);
984  	   SCIPfreeBufferArray(scip, &violrowpos);
985  	   SCIPfreeBufferArray(scip, &violrows);
986  	   SCIPfreeBufferArray(scip, &activities);
987  	
988  	   return SCIP_OKAY;
989  	}
990  	
991  	
992  	/*
993  	 * heuristic specific interface methods
994  	 */
995  	
996  	/** creates the shifting heuristic with infeasibility recovering and includes it in SCIP */
997  	SCIP_RETCODE SCIPincludeHeurShifting(
998  	   SCIP*                 scip                /**< SCIP data structure */
999  	   )
1000 	{
1001 	   SCIP_HEUR* heur;
1002 	
1003 	   /* include primal heuristic */
1004 	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
1005 	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
1006 	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecShifting, NULL) );
1007 	
1008 	   assert(heur != NULL);
1009 	
1010 	   /* set non-NULL pointers to callback methods */
1011 	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyShifting) );
1012 	   SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitShifting) );
1013 	   SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitShifting) );
1014 	   SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolShifting) );
1015 	
1016 	   return SCIP_OKAY;
1017 	}
1018 	
1019