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_zirounding.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  zirounding primal heuristic
28   	 * @author Gregor Hendel
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_zirounding.h"
35   	#include "scip/pub_heur.h"
36   	#include "scip/pub_lp.h"
37   	#include "scip/pub_message.h"
38   	#include "scip/pub_var.h"
39   	#include "scip/scip_branch.h"
40   	#include "scip/scip_heur.h"
41   	#include "scip/scip_lp.h"
42   	#include "scip/scip_mem.h"
43   	#include "scip/scip_message.h"
44   	#include "scip/scip_numerics.h"
45   	#include "scip/scip_param.h"
46   	#include "scip/scip_sol.h"
47   	#include "scip/scip_solvingstats.h"
48   	#include <string.h>
49   	
50   	#define HEUR_NAME             "zirounding"
51   	#define HEUR_DESC             "LP rounding heuristic as suggested by C. Wallace taking row slacks and bounds into account"
52   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_ROUNDING
53   	#define HEUR_PRIORITY         -500
54   	#define HEUR_FREQ             1
55   	#define HEUR_FREQOFS          0
56   	#define HEUR_MAXDEPTH         -1
57   	#define HEUR_TIMING           SCIP_HEURTIMING_AFTERLPNODE
58   	#define HEUR_USESSUBSCIP      FALSE  /**< does the heuristic use a secondary SCIP instance? */
59   	
60   	#define DEFAULT_MAXROUNDINGLOOPS   2      /**< delimits the number of main loops, can be set to -1 for no limit */
61   	#define DEFAULT_STOPZIROUND        TRUE   /**< deactivation check is enabled by default */
62   	#define DEFAULT_STOPPERCENTAGE     0.02   /**< the tolerance percentage after which zirounding will not be executed anymore */
63   	#define DEFAULT_MINSTOPNCALLS      1000   /**< number of heuristic calls before deactivation check */
64   	
65   	#define MINSHIFT                   1e-4   /**< minimum shift value for every single step */
66   	
67   	/*
68   	 * Data structures
69   	 */
70   	
71   	/** primal heuristic data */
72   	struct SCIP_HeurData
73   	{
74   	   SCIP_SOL*             sol;                /**< working solution */
75   	   SCIP_Longint          lastlp;             /**< the number of the last LP for which ZIRounding was called */
76   	   int                   maxroundingloops;   /**< limits rounding loops in execution */
77   	   SCIP_Bool             stopziround;        /**< sets deactivation check */
78   	   SCIP_Real             stoppercentage;     /**< threshold for deactivation check */
79   	   int                   minstopncalls;      /**< number of heuristic calls before deactivation check */
80   	};
81   	
82   	enum Direction
83   	{
84   	   DIRECTION_UP          =  1,
85   	   DIRECTION_DOWN        = -1
86   	};
87   	typedef enum Direction DIRECTION;
88   	
89   	/*
90   	 * Local methods
91   	 */
92   	
93   	/** returns the fractionality of a value x, which is calculated as zivalue(x) = min(x-floor(x), ceil(x)-x) */
94   	static
95   	SCIP_Real getZiValue(
96   	   SCIP*                 scip,               /**< pointer to current SCIP data structure */
97   	   SCIP_Real             val                 /**< the value for which the fractionality should be computed */
98   	   )
99   	{
100  	   SCIP_Real upgap;     /* the gap between val and ceil(val) */
101  	   SCIP_Real downgap;   /* the gap between val and floor(val) */
102  	
103  	   assert(scip != NULL);
104  	
105  	   upgap   = SCIPfeasCeil(scip, val) - val;
106  	   downgap = val - SCIPfeasFloor(scip, val);
107  	
108  	   return MIN(upgap, downgap);
109  	}
110  	
111  	/** determines shifting bounds for variable */
112  	static
113  	void calculateBounds(
114  	   SCIP*                 scip,               /**< pointer to current SCIP data structure */
115  	   SCIP_VAR*             var,                /**< the variable for which lb and ub have to be calculated */
116  	   SCIP_Real             currentvalue,       /**< the current value of var in the working solution */
117  	   SCIP_Real*            upperbound,         /**< pointer to store the calculated upper bound on the variable shift */
118  	   SCIP_Real*            lowerbound,         /**< pointer to store the calculated lower bound on the variable shift */
119  	   SCIP_Real*            upslacks,           /**< array that contains the slacks between row activities and the right hand sides of the rows */
120  	   SCIP_Real*            downslacks,         /**< array that contains lhs slacks */
121  	   int                   nslacks,            /**< current number of slacks */
122  	   SCIP_Bool*            numericalerror      /**< flag to determine whether a numerical error occurred */
123  	   )
124  	{
125  	   SCIP_COL*      col;
126  	   SCIP_ROW**     colrows;
127  	   SCIP_Real*     colvals;
128  	   int            ncolvals;
129  	   int i;
130  	
131  	   assert(scip != NULL);
132  	   assert(var != NULL);
133  	   assert(upslacks != NULL);
134  	   assert(downslacks != NULL);
135  	   assert(upperbound != NULL);
136  	   assert(lowerbound != NULL);
137  	
138  	   /* get the column associated to the variable, the nonzero rows and the nonzero coefficients */
139  	   col       = SCIPvarGetCol(var);
140  	   colrows   = SCIPcolGetRows(col);
141  	   colvals   = SCIPcolGetVals(col);
142  	   ncolvals  = SCIPcolGetNLPNonz(col);
143  	
144  	   /* only proceed, when variable has nonzero coefficients */
145  	   if( ncolvals == 0 )
146  	      return;
147  	
148  	   assert(colvals != NULL);
149  	   assert(colrows != NULL);
150  	
151  	   /* initialize the bounds on the shift to be the gap of the current solution value to the bounds of the variable */
152  	   if( SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) )
153  	      *upperbound = SCIPinfinity(scip);
154  	   else
155  	      *upperbound = SCIPvarGetUbGlobal(var) - currentvalue;
156  	
157  	   if( SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) )
158  	      *lowerbound = SCIPinfinity(scip);
159  	   else
160  	      *lowerbound = currentvalue - SCIPvarGetLbGlobal(var);
161  	
162  	   /* go through every nonzero row coefficient corresponding to var to determine bounds for shifting
163  	    * in such a way that shifting maintains feasibility in every LP row.
164  	    * a lower or upper bound as it is calculated in zirounding always has to be >= 0.0.
165  	    * if one of these values is significantly < 0.0, this will cause the abort of execution of the heuristic so that
166  	    * infeasible solutions are avoided
167  	    */
168  	   for( i = 0; i < ncolvals && MAX(*lowerbound, *upperbound) >= MINSHIFT; ++i )
169  	   {
170  	      SCIP_ROW* row;
171  	      int       rowpos;
172  	
173  	      row = colrows[i];
174  	      rowpos = SCIProwGetLPPos(row);
175  	
176  	      /* the row might currently not be in the LP, ignore it! */
177  	      if( rowpos == -1 )
178  	         continue;
179  	
180  	      assert(0 <= rowpos && rowpos < nslacks);
181  	
182  	      /* all bounds and slacks as they are calculated in zirounding always have to be greater equal zero.
183  	       * It might however be due to numerical issues, e.g. with scaling, that they are not. Better abort in this case.
184  	       */
185  	      if( SCIPisFeasLT(scip, *lowerbound, 0.0) || SCIPisFeasLT(scip, *upperbound, 0.0)
186  	         || SCIPisFeasLT(scip, upslacks[rowpos], 0.0) || SCIPisFeasLT(scip, downslacks[rowpos] , 0.0) )
187  	      {
188  	         *numericalerror = TRUE;
189  	         return;
190  	      }
191  	
192  	      SCIPdebugMsg(scip, "colval: %15.8g, downslack: %15.8g, upslack: %5.2g, lb: %5.2g, ub: %5.2g\n", colvals[i], downslacks[rowpos], upslacks[rowpos],
193  	         *lowerbound, *upperbound);
194  	
195  	      /* if coefficient > 0, rounding up might violate up slack and rounding down might violate down slack
196  	       * thus search for the minimum so that no constraint is violated; vice versa for coefficient < 0
197  	       */
198  	      if( colvals[i] > 0 )
199  	      {
200  	         if( !SCIPisInfinity(scip, upslacks[rowpos]) )
201  	         {
202  	            SCIP_Real upslack;
203  	            upslack = MAX(upslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */
204  	            *upperbound = MIN(*upperbound, upslack/colvals[i]);
205  	         }
206  	
207  	         if( !SCIPisInfinity(scip, downslacks[rowpos]) )
208  	         {
209  	            SCIP_Real downslack;
210  	            downslack = MAX(downslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */
211  	            *lowerbound = MIN(*lowerbound, downslack/colvals[i]);
212  	         }
213  	      }
214  	      else
215  	      {
216  	         assert(colvals[i] != 0.0);
217  	
218  	         if( !SCIPisInfinity(scip, upslacks[rowpos]) )
219  	         {
220  	            SCIP_Real upslack;
221  	            upslack = MAX(upslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */
222  	            *lowerbound = MIN(*lowerbound, -upslack/colvals[i]);
223  	         }
224  	
225  	         if( !SCIPisInfinity(scip, downslacks[rowpos]) )
226  	         {
227  	            SCIP_Real downslack;
228  	            downslack = MAX(downslacks[rowpos], 0.0); /* avoid errors due to numerically slightly infeasible rows */
229  	            *upperbound = MIN(*upperbound, -downslack/colvals[i]);
230  	         }
231  	      }
232  	   }
233  	}
234  	
235  	/**  when a variable is shifted, the activities and slacks of all rows it appears in have to be updated */
236  	static
237  	SCIP_RETCODE updateSlacks(
238  	   SCIP*                 scip,               /**< pointer to current SCIP data structure */
239  	   SCIP_SOL*             sol,                /**< working solution */
240  	   SCIP_VAR*             var,                /**< pointer to variable to be modified */
241  	   SCIP_Real             shiftvalue,         /**< the value by which the variable is shifted */
242  	   SCIP_Real*            upslacks,           /**< upslacks of all rows the variable appears in */
243  	   SCIP_Real*            downslacks,         /**< downslacks of all rows the variable appears in */
244  	   SCIP_Real*            activities,         /**< activities of the LP rows */
245  	   SCIP_VAR**            slackvars,          /**< the slack variables for equality rows */
246  	   SCIP_Real*            slackcoeffs,        /**< the slack variable coefficients */
247  	   int                   nslacks             /**< size of the arrays */
248  	   )
249  	{
250  	   SCIP_COL*    col;        /* the corresponding column of variable var */
251  	   SCIP_ROW**   rows;       /* pointer to the nonzero coefficient rows for variable var */
252  	   int          nrows;      /* the number of nonzeros */
253  	   SCIP_Real*   colvals;    /* array to store the nonzero coefficients */
254  	   int i;
255  	
256  	   assert(scip != NULL);
257  	   assert(sol != NULL);
258  	   assert(var != NULL);
259  	   assert(upslacks != NULL);
260  	   assert(downslacks != NULL);
261  	   assert(activities != NULL);
262  	   assert(nslacks >= 0);
263  	
264  	   col = SCIPvarGetCol(var);
265  	   assert(col != NULL);
266  	
267  	   rows     = SCIPcolGetRows(col);
268  	   nrows    = SCIPcolGetNLPNonz(col);
269  	   colvals  = SCIPcolGetVals(col);
270  	   assert(nrows == 0 || (rows != NULL && colvals != NULL));
271  	
272  	   /* go through all rows the shifted variable appears in */
273  	   for( i = 0; i < nrows; ++i )
274  	   {
275  	      int rowpos;
276  	
277  	      rowpos = SCIProwGetLPPos(rows[i]);
278  	      assert(-1 <= rowpos && rowpos < nslacks);
279  	
280  	      /* if the row is in the LP, update its activity, up and down slack */
281  	      if( rowpos >= 0 )
282  	      {
283  	         SCIP_Real val;
284  	         SCIP_Real lhs;
285  	         SCIP_Real rhs;
286  	         SCIP_ROW* row;
287  	
288  	         val = colvals[i] * shiftvalue;
289  	         row = rows[i];
290  	         lhs = SCIProwGetLhs(row);
291  	         rhs = SCIProwGetRhs(row);
292  	
293  	         /* if the row is an equation, we update its slack variable instead of its activities */
294  	         if( SCIPisFeasEQ(scip, lhs, rhs) )
295  	         {
296  	            SCIP_Real slackvarshiftval;
297  	            SCIP_Real slackvarsolval;
298  	
299  	            assert(slackvars[rowpos] != NULL);
300  	            assert(!SCIPisZero(scip, slackcoeffs[rowpos]));
301  	
302  	            slackvarsolval = SCIPgetSolVal(scip, sol, slackvars[rowpos]);
303  	            slackvarshiftval = -val / slackcoeffs[rowpos];
304  	
305  	            assert(SCIPisFeasGE(scip, slackvarsolval + slackvarshiftval, SCIPvarGetLbGlobal(slackvars[rowpos])));
306  	            assert(SCIPisFeasLE(scip, slackvarsolval + slackvarshiftval, SCIPvarGetUbGlobal(slackvars[rowpos])));
307  	
308  	            SCIP_CALL( SCIPsetSolVal(scip, sol, slackvars[rowpos], slackvarsolval + slackvarshiftval) );
309  	         }
310  	         else if( !SCIPisInfinity(scip, REALABS(activities[rowpos])) )
311  	            activities[rowpos] += val;
312  	
313  	         /* the slacks of the row now can be updated independently of its type */
314  	         if( !SCIPisInfinity(scip, upslacks[rowpos]) )
315  	            upslacks[rowpos] -= val;
316  	         if( !SCIPisInfinity(scip, -downslacks[rowpos]) )
317  	            downslacks[rowpos] += val;
318  	
319  	         assert(SCIPisInfinity(scip, -lhs) || SCIPisFeasGE(scip, activities[rowpos], lhs));
320  	         assert(SCIPisInfinity(scip, rhs) || SCIPisFeasLE(scip, activities[rowpos], rhs));
321  	      }
322  	   }
323  	   return SCIP_OKAY;
324  	}
325  	
326  	/** finds a continuous slack variable for an equation row, NULL if none exists */
327  	static
328  	void rowFindSlackVar(
329  	   SCIP*                 scip,               /**< pointer to current SCIP data structure */
330  	   SCIP_ROW*             row,                /**< the row for which a slack variable is searched */
331  	   SCIP_VAR**            varpointer,         /**< pointer to store the slack variable */
332  	   SCIP_Real*            coeffpointer        /**< pointer to store the coefficient of the slack variable */
333  	   )
334  	{
335  	   int v;
336  	   SCIP_COL** rowcols;
337  	   SCIP_Real* rowvals;
338  	   int nrowvals;
339  	
340  	   assert(row != NULL);
341  	   assert(varpointer != NULL);
342  	   assert(coeffpointer != NULL);
343  	
344  	   rowcols = SCIProwGetCols(row);
345  	   rowvals = SCIProwGetVals(row);
346  	   nrowvals = SCIProwGetNNonz(row);
347  	
348  	   assert(nrowvals == 0 || rowvals != NULL);
349  	   assert(nrowvals == 0 || rowcols != NULL);
350  	
351  	   /* iterate over the row variables. Stop after the first unfixed continuous variable was found. */
352  	   for( v = nrowvals - 1; v >= 0; --v )
353  	   {
354  	      SCIP_VAR* colvar;
355  	
356  	      assert(rowcols[v] != NULL);
357  	      if( SCIPcolGetLPPos(rowcols[v]) == -1 )
358  	         continue;
359  	
360  	      colvar = SCIPcolGetVar(rowcols[v]);
361  	
362  	      if( SCIPvarGetType(colvar) == SCIP_VARTYPE_CONTINUOUS
363  	         && !SCIPisFeasEQ(scip, SCIPvarGetLbGlobal(colvar), SCIPvarGetUbGlobal(colvar))
364  	         && SCIPcolGetNLPNonz(rowcols[v]) == 1 )
365  	      {
366  	         SCIPdebugMsg(scip, "  slack variable for row %s found: %s\n", SCIProwGetName(row), SCIPvarGetName(colvar));
367  	
368  	         *coeffpointer = rowvals[v];
369  	         *varpointer = colvar;
370  	
371  	         return;
372  	      }
373  	   }
374  	
375  	   *varpointer = NULL;
376  	   *coeffpointer = 0.0;
377  	
378  	   SCIPdebugMsg(scip, "No slack variable for row %s found. \n", SCIProwGetName(row));
379  	}
380  	
381  	/*
382  	 * Callback methods of primal heuristic
383  	 */
384  	
385  	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
386  	static
387  	SCIP_DECL_HEURCOPY(heurCopyZirounding)
388  	{  /*lint --e{715}*/
389  	   assert(scip != NULL);
390  	   assert(heur != NULL);
391  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
392  	
393  	   /* call inclusion method of primal heuristic */
394  	   SCIP_CALL( SCIPincludeHeurZirounding(scip) );
395  	
396  	   return SCIP_OKAY;
397  	}
398  	
399  	/** destructor of primal heuristic to free user data (called when SCIP is exiting) */
400  	static
401  	SCIP_DECL_HEURFREE(heurFreeZirounding)
402  	{  /*lint --e{715}*/
403  	   SCIP_HEURDATA* heurdata;
404  	
405  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
406  	
407  	   heurdata = SCIPheurGetData(heur);
408  	   assert(heurdata != NULL);
409  	
410  	   /* free heuristic data */
411  	   SCIPfreeBlockMemory(scip, &heurdata);
412  	   SCIPheurSetData(heur, NULL);
413  	
414  	   return SCIP_OKAY;
415  	}
416  	
417  	/** initialization method of primal heuristic (called after problem was transformed) */
418  	static
419  	SCIP_DECL_HEURINIT(heurInitZirounding)
420  	{  /*lint --e{715}*/
421  	   SCIP_HEURDATA* heurdata;
422  	
423  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
424  	
425  	   heurdata = SCIPheurGetData(heur);
426  	   assert(heurdata != NULL);
427  	
428  	   /* create working solution */
429  	   SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
430  	
431  	   return SCIP_OKAY;
432  	}
433  	
434  	/** deinitialization method of primal heuristic (called before transformed problem is freed) */
435  	static
436  	SCIP_DECL_HEUREXIT(heurExitZirounding)  /*lint --e{715}*/
437  	{  /*lint --e{715}*/
438  	   SCIP_HEURDATA* heurdata;
439  	
440  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
441  	
442  	   heurdata = SCIPheurGetData(heur);
443  	   assert(heurdata != NULL);
444  	
445  	   /* free working solution */
446  	   SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
447  	
448  	   return SCIP_OKAY;
449  	}
450  	
451  	/** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
452  	static
453  	SCIP_DECL_HEURINITSOL(heurInitsolZirounding)
454  	{  /*lint --e{715}*/
455  	   SCIP_HEURDATA* heurdata;
456  	
457  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
458  	
459  	   heurdata = SCIPheurGetData(heur);
460  	   assert(heurdata != NULL);
461  	
462  	   heurdata->lastlp = -1;
463  	
464  	   return SCIP_OKAY;
465  	}
466  	
467  	
468  	/** execution method of primal heuristic */
469  	static
470  	SCIP_DECL_HEUREXEC(heurExecZirounding)
471  	{  /*lint --e{715}*/
472  	   SCIP_HEURDATA*     heurdata;
473  	   SCIP_SOL*          sol;
474  	   SCIP_VAR**         lpcands;
475  	   SCIP_VAR**         zilpcands;
476  	
477  	   SCIP_VAR**         slackvars;
478  	   SCIP_Real*         upslacks;
479  	   SCIP_Real*         downslacks;
480  	   SCIP_Real*         activities;
481  	   SCIP_Real*         slackvarcoeffs;
482  	   SCIP_Bool*         rowneedsslackvar;
483  	
484  	   SCIP_ROW**         rows;
485  	   SCIP_Real*         lpcandssol;
486  	   SCIP_Real*         solarray;
487  	
488  	   SCIP_Longint       nlps;
489  	   int                currentlpcands;
490  	   int                nlpcands;
491  	   int                nimplfracs;
492  	   int                i;
493  	   int                c;
494  	   int                nslacks;
495  	   int                nroundings;
496  	
497  	   SCIP_Bool          improvementfound;
498  	   SCIP_Bool          numericalerror;
499  	
500  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
501  	   assert(result != NULL);
502  	   assert(SCIPhasCurrentNodeLP(scip));
503  	
504  	   *result = SCIP_DIDNOTRUN;
505  	
506  	   /* do not call heuristic of node was already detected to be infeasible */
507  	   if( nodeinfeasible )
508  	      return SCIP_OKAY;
509  	
510  	   /* only call heuristic if an optimal LP-solution is at hand */
511  	   if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
512  	      return SCIP_OKAY;
513  	
514  	   /* only call heuristic, if the LP objective value is smaller than the cutoff bound */
515  	   if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) )
516  	      return SCIP_OKAY;
517  	
518  	   /* get heuristic data */
519  	   heurdata = SCIPheurGetData(heur);
520  	   assert(heurdata != NULL);
521  	
522  	   /* Do not call heuristic if deactivation check is enabled and percentage of found solutions in relation
523  	    * to number of calls falls below heurdata->stoppercentage */
524  	   if( heurdata->stopziround && SCIPheurGetNCalls(heur) >= heurdata->minstopncalls
525  	      && SCIPheurGetNSolsFound(heur)/(SCIP_Real)SCIPheurGetNCalls(heur) < heurdata->stoppercentage )
526  	      return SCIP_OKAY;
527  	
528  	   /* assure that heuristic has not already been called after the last LP had been solved */
529  	   nlps = SCIPgetNLPs(scip);
530  	   if( nlps == heurdata->lastlp )
531  	      return SCIP_OKAY;
532  	
533  	   heurdata->lastlp = nlps;
534  	
535  	   /* get fractional variables */
536  	   SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, NULL, &nlpcands, NULL, &nimplfracs) );
537  	   nlpcands = nlpcands + nimplfracs;
538  	   /* make sure that there is at least one fractional variable that should be integral */
539  	   if( nlpcands == 0 )
540  	      return SCIP_OKAY;
541  	
542  	   assert(nlpcands > 0);
543  	   assert(lpcands != NULL);
544  	   assert(lpcandssol != NULL);
545  	
546  	   /* get LP rows data */
547  	   rows    = SCIPgetLPRows(scip);
548  	   nslacks = SCIPgetNLPRows(scip);
549  	
550  	   /* cannot do anything if LP is empty */
551  	   if( nslacks == 0 )
552  	      return SCIP_OKAY;
553  	
554  	   assert(rows != NULL);
555  	   assert(nslacks > 0);
556  	
557  	   /* get the working solution from heuristic's local data */
558  	   sol = heurdata->sol;
559  	   assert(sol != NULL);
560  	
561  	   *result = SCIP_DIDNOTFIND;
562  	
563  	   solarray = NULL;
564  	   zilpcands = NULL;
565  	
566  	   /* copy the current LP solution to the working solution and allocate memory for local data */
567  	   SCIP_CALL( SCIPlinkLPSol(scip, sol) );
568  	   SCIP_CALL( SCIPallocBufferArray(scip, &solarray, nlpcands) );
569  	   SCIP_CALL( SCIPallocBufferArray(scip, &zilpcands, nlpcands) );
570  	
571  	   /* copy necessary data to local arrays */
572  	   BMScopyMemoryArray(solarray, lpcandssol, nlpcands);
573  	   BMScopyMemoryArray(zilpcands, lpcands, nlpcands);
574  	
575  	   /* allocate buffer data arrays */
576  	   SCIP_CALL( SCIPallocBufferArray(scip, &slackvars, nslacks) );
577  	   SCIP_CALL( SCIPallocBufferArray(scip, &upslacks, nslacks) );
578  	   SCIP_CALL( SCIPallocBufferArray(scip, &downslacks, nslacks) );
579  	   SCIP_CALL( SCIPallocBufferArray(scip, &slackvarcoeffs, nslacks) );
580  	   SCIP_CALL( SCIPallocBufferArray(scip, &rowneedsslackvar, nslacks) );
581  	   SCIP_CALL( SCIPallocBufferArray(scip, &activities, nslacks) );
582  	
583  	   BMSclearMemoryArray(slackvars, nslacks);
584  	   BMSclearMemoryArray(slackvarcoeffs, nslacks);
585  	   BMSclearMemoryArray(rowneedsslackvar, nslacks);
586  	
587  	   numericalerror = FALSE;
588  	   nroundings = 0;
589  	
590  	   /* loop over fractional variables and involved LP rows to find all rows which require a slack variable */
591  	   for( c = 0; c < nlpcands; ++c )
592  	   {
593  	      SCIP_VAR* cand;
594  	      SCIP_ROW** candrows;
595  	      int r;
596  	      int ncandrows;
597  	
598  	      cand = zilpcands[c];
599  	      assert(cand != NULL);
600  	      assert(SCIPcolGetLPPos(SCIPvarGetCol(cand)) >= 0);
601  	
602  	      candrows = SCIPcolGetRows(SCIPvarGetCol(cand));
603  	      ncandrows = SCIPcolGetNLPNonz(SCIPvarGetCol(cand));
604  	
605  	      assert(candrows == NULL || ncandrows > 0);
606  	
607  	      for( r = 0; r < ncandrows; ++r )
608  	      {
609  	         int rowpos;
610  	
611  	         assert(candrows != NULL); /* to please flexelint */
612  	         assert(candrows[r] != NULL);
613  	         rowpos = SCIProwGetLPPos(candrows[r]);
614  	
615  	         if( rowpos >= 0 && SCIPisFeasEQ(scip, SCIProwGetLhs(candrows[r]), SCIProwGetRhs(candrows[r])) )
616  	         {
617  	            rowneedsslackvar[rowpos] = TRUE;
618  	            SCIPdebugMsg(scip, "  Row %s needs slack variable for variable %s\n", SCIProwGetName(candrows[r]), SCIPvarGetName(cand));
619  	         }
620  	      }
621  	   }
622  	
623  	   /* calculate row slacks for every every row that belongs to the current LP and ensure, that the current solution
624  	    * has no violated constraint -- if any constraint is violated, i.e. a slack is significantly smaller than zero,
625  	    * this will cause the termination of the heuristic because Zirounding does not provide feasibility recovering
626  	    */
627  	   for( i = 0; i < nslacks; ++i )
628  	   {
629  	      SCIP_ROW*          row;
630  	      SCIP_Real          lhs;
631  	      SCIP_Real          rhs;
632  	
633  	      row = rows[i];
634  	
635  	      assert(row != NULL);
636  	
637  	      lhs = SCIProwGetLhs(row);
638  	      rhs = SCIProwGetRhs(row);
639  	
640  	      /* get row activity */
641  	      activities[i] = SCIPgetRowActivity(scip, row);
642  	      assert(SCIPisFeasLE(scip, lhs, activities[i]) && SCIPisFeasLE(scip, activities[i], rhs));
643  	
644  	      /* in special case if LHS or RHS is (-)infinity slacks have to be initialized as infinity */
645  	      if( SCIPisInfinity(scip, -lhs) )
646  	         downslacks[i] = SCIPinfinity(scip);
647  	      else
648  	         downslacks[i] = activities[i] - lhs;
649  	
650  	      if( SCIPisInfinity(scip, rhs) )
651  	         upslacks[i] = SCIPinfinity(scip);
652  	      else
653  	         upslacks[i] = rhs - activities[i];
654  	
655  	      SCIPdebugMsg(scip, "lhs:%5.2f <= act:%5.2g <= rhs:%5.2g --> down: %5.2g, up:%5.2g\n", lhs, activities[i], rhs, downslacks[i], upslacks[i]);
656  	
657  	      /* row is an equation. Try to find a slack variable in the row, i.e.,
658  	       * a continuous variable which occurs only in this row. If no such variable exists,
659  	       * there is no hope for an IP-feasible solution in this round
660  	       */
661  	      if( SCIPisFeasEQ(scip, lhs, rhs) && rowneedsslackvar[i] )
662  	      {
663  	         /* @todo: This is only necessary for rows containing fractional variables. */
664  	         rowFindSlackVar(scip, row, &(slackvars[i]), &(slackvarcoeffs[i]));
665  	
666  	         if( slackvars[i] == NULL )
667  	         {
668  	            SCIPdebugMsg(scip, "No slack variable found for equation %s, terminating ZI Round heuristic\n", SCIProwGetName(row));
669  	            goto TERMINATE;
670  	         }
671  	         else
672  	         {
673  	            SCIP_Real ubslackvar;
674  	            SCIP_Real lbslackvar;
675  	            SCIP_Real solvalslackvar;
676  	            SCIP_Real coeffslackvar;
677  	            SCIP_Real ubgap;
678  	            SCIP_Real lbgap;
679  	
680  	            assert(SCIPvarGetType(slackvars[i]) == SCIP_VARTYPE_CONTINUOUS);
681  	            solvalslackvar = SCIPgetSolVal(scip, sol, slackvars[i]);
682  	            ubslackvar = SCIPvarGetUbGlobal(slackvars[i]);
683  	            lbslackvar = SCIPvarGetLbGlobal(slackvars[i]);
684  	
685  	            coeffslackvar = slackvarcoeffs[i];
686  	            assert(!SCIPisZero(scip, coeffslackvar));
687  	
688  	            ubgap = MAX(0.0, ubslackvar - solvalslackvar);
689  	            lbgap = MAX(0.0, solvalslackvar - lbslackvar);
690  	
691  	            if( SCIPisPositive(scip, coeffslackvar) )
692  	            {
693  	              if( !SCIPisInfinity(scip, lbslackvar) )
694  	                upslacks[i] += coeffslackvar * lbgap;
695  	              else
696  	                upslacks[i] = SCIPinfinity(scip);
697  	              if( !SCIPisInfinity(scip, ubslackvar) )
698  	                downslacks[i] += coeffslackvar * ubgap;
699  	              else
700  	                downslacks[i] = SCIPinfinity(scip);
701  	            }
702  	            else
703  	            {
704  	               if( !SCIPisInfinity(scip, ubslackvar) )
705  	                  upslacks[i] -= coeffslackvar * ubgap;
706  	               else
707  	                  upslacks[i] = SCIPinfinity(scip);
708  	               if( !SCIPisInfinity(scip, lbslackvar) )
709  	                  downslacks[i] -= coeffslackvar * lbgap;
710  	               else
711  	                  downslacks[i] = SCIPinfinity(scip);
712  	            }
713  	            SCIPdebugMsg(scip, "  Slack variable for row %s at pos %d: %g <= %s = %g <= %g; Coeff %g, upslack = %g, downslack = %g  \n",
714  	               SCIProwGetName(row), SCIProwGetLPPos(row), lbslackvar, SCIPvarGetName(slackvars[i]), solvalslackvar, ubslackvar, coeffslackvar,
715  	               upslacks[i], downslacks[i]);
716  	         }
717  	      }
718  	      /* due to numerical inaccuracies, the rows might be feasible, even if the slacks are
719  	       * significantly smaller than zero -> terminate
720  	       */
721  	      if( SCIPisFeasLT(scip, upslacks[i], 0.0) || SCIPisFeasLT(scip, downslacks[i], 0.0) )
722  	         goto TERMINATE;
723  	   }
724  	
725  	   assert(nslacks == 0 || (upslacks != NULL && downslacks != NULL && activities != NULL));
726  	
727  	   /* initialize number of remaining variables and flag to enter the main loop */
728  	   currentlpcands = nlpcands;
729  	   improvementfound = TRUE;
730  	
731  	   /* iterate over variables as long as there are fractional variables left */
732  	   while( currentlpcands > 0 && improvementfound && (heurdata->maxroundingloops == -1 || nroundings < heurdata->maxroundingloops) )
733  	   {  /*lint --e{850}*/
734  	      improvementfound = FALSE;
735  	      nroundings++;
736  	      SCIPdebugMsg(scip, "zirounding enters while loop for %d time with %d candidates left. \n", nroundings, currentlpcands);
737  	
738  	      /* check for every remaining fractional variable if a shifting decreases ZI-value of the variable */
739  	      for( c = 0; c < currentlpcands; ++c )
740  	      {
741  	         SCIP_VAR* var;
742  	         SCIP_Real oldsolval;
743  	         SCIP_Real upperbound;
744  	         SCIP_Real lowerbound;
745  	         SCIP_Real up;
746  	         SCIP_Real down;
747  	         SCIP_Real ziup;
748  	         SCIP_Real zidown;
749  	         SCIP_Real zicurrent;
750  	         SCIP_Real shiftval;
751  	
752  	         DIRECTION direction;
753  	
754  	         /* get values from local data */
755  	         oldsolval = solarray[c];
756  	         var = zilpcands[c];
757  	
758  	         assert(!SCIPisFeasIntegral(scip, oldsolval));
759  	         assert(SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN);
760  	
761  	         /* calculate bounds for variable and make sure that there are no numerical inconsistencies */
762  	         upperbound = SCIPinfinity(scip);
763  	         lowerbound = SCIPinfinity(scip);
764  	         calculateBounds(scip, var, oldsolval, &upperbound, &lowerbound, upslacks, downslacks, nslacks, &numericalerror);
765  	
766  	         if( numericalerror )
767  	            goto TERMINATE;
768  	
769  	         /* continue if only marginal shifts are possible */
770  	         if( MAX(upperbound, lowerbound) < MINSHIFT )
771  	         {
772  	            /* stop immediately if a variable has not been rounded during the last round */
773  	            if( nroundings == heurdata->maxroundingloops )
774  	               break;
775  	            else
776  	               continue;
777  	         }
778  	
779  	         /* calculate the possible values after shifting */
780  	         up   = oldsolval + upperbound;
781  	         down = oldsolval - lowerbound;
782  	
783  	         /* if the variable is integer or implicit binary, do not shift further than the nearest integer */
784  	         if( SCIPvarGetType(var) != SCIP_VARTYPE_BINARY)
785  	         {
786  	            SCIP_Real ceilx;
787  	            SCIP_Real floorx;
788  	
789  	            ceilx = SCIPfeasCeil(scip, oldsolval);
790  	            floorx = SCIPfeasFloor(scip, oldsolval);
791  	            up   = MIN(up, ceilx);
792  	            down = MAX(down, floorx);
793  	         }
794  	
795  	         /* calculate necessary values */
796  	         ziup      = getZiValue(scip, up);
797  	         zidown    = getZiValue(scip, down);
798  	         zicurrent = getZiValue(scip, oldsolval);
799  	
800  	         /* calculate the shifting direction that reduces ZI-value the most,
801  	          * if both directions improve ZI-value equally, take the direction which improves the objective
802  	          */
803  	         if( SCIPisFeasLT(scip, zidown, zicurrent) || SCIPisFeasLT(scip, ziup, zicurrent) )
804  	         {
805  	            if( SCIPisFeasEQ(scip,ziup, zidown) )
806  	               direction  = SCIPisFeasGE(scip, SCIPvarGetObj(var), 0.0) ? DIRECTION_DOWN : DIRECTION_UP;
807  	            else if( SCIPisFeasLT(scip, zidown, ziup) )
808  	               direction = DIRECTION_DOWN;
809  	            else
810  	               direction = DIRECTION_UP;
811  	
812  	            /* once a possible shifting direction and value have been found, variable value is updated */
813  	            shiftval = (direction == DIRECTION_UP ? up - oldsolval : down - oldsolval);
814  	
815  	            /* this improves numerical stability in some cases */
816  	            if( direction == DIRECTION_UP )
817  	               shiftval = MIN(shiftval, upperbound);
818  	            else
819  	               shiftval = MIN(shiftval, lowerbound);
820  	            /* update the solution */
821  	            solarray[c] = direction == DIRECTION_UP ? up : down;
822  	            SCIP_CALL( SCIPsetSolVal(scip, sol, var, solarray[c]) );
823  	
824  	            /* update the rows activities and slacks */
825  	            SCIP_CALL( updateSlacks(scip, sol, var, shiftval, upslacks,
826  	                  downslacks, activities, slackvars, slackvarcoeffs, nslacks) );
827  	
828  	            SCIPdebugMsg(scip, "zirounding update step : %d var index, oldsolval=%g, shiftval=%g\n",
829  	               SCIPvarGetIndex(var), oldsolval, shiftval);
830  	            /* since at least one improvement has been found, heuristic will enter main loop for another time because the improvement
831  	             * might affect many LP rows and their current slacks and thus make further rounding steps possible */
832  	            improvementfound = TRUE;
833  	         }
834  	
835  	         /* if solution value of variable has become feasibly integral due to rounding step,
836  	          * variable is put at the end of remaining candidates array so as not to be considered in future loops
837  	          */
838  	         if( SCIPisFeasIntegral(scip, solarray[c]) )
839  	         {
840  	            zilpcands[c] = zilpcands[currentlpcands - 1];
841  	            solarray[c] = solarray[currentlpcands - 1];
842  	            currentlpcands--;
843  	
844  	            /* counter is decreased if end of candidates array has not been reached yet */
845  	            if( c < currentlpcands )
846  	               c--;
847  	         }
848  	         else if( nroundings == heurdata->maxroundingloops )
849  	            goto TERMINATE;
850  	      }
851  	   }
852  	
853  	   /* in case that no candidate is left for rounding after the final main loop
854  	    * the found solution has to be checked for feasibility in the original problem
855  	    */
856  	   if( currentlpcands == 0 )
857  	   {
858  	      SCIP_Bool stored;
859  	      SCIP_CALL(SCIPtrySol(scip, sol, FALSE, FALSE, FALSE, TRUE, FALSE, &stored));
860  	      if( stored )
861  	      {
862  	#ifdef SCIP_DEBUG
863  	         SCIPdebugMsg(scip, "found feasible rounded solution:\n");
864  	         SCIP_CALL( SCIPprintSol(scip, sol, NULL, FALSE) );
865  	#endif
866  	         SCIPstatisticMessage("  ZI Round solution value: %g \n", SCIPgetSolOrigObj(scip, sol));
867  	
868  	         *result = SCIP_FOUNDSOL;
869  	      }
870  	   }
871  	
872  	   /* free memory for all locally allocated data */
873  	 TERMINATE:
874  	   SCIPfreeBufferArrayNull(scip, &activities);
875  	   SCIPfreeBufferArrayNull(scip, &rowneedsslackvar);
876  	   SCIPfreeBufferArrayNull(scip, &slackvarcoeffs);
877  	   SCIPfreeBufferArrayNull(scip, &downslacks);
878  	   SCIPfreeBufferArrayNull(scip, &upslacks);
879  	   SCIPfreeBufferArrayNull(scip, &slackvars);
880  	   SCIPfreeBufferArrayNull(scip, &zilpcands);
881  	   SCIPfreeBufferArrayNull(scip, &solarray);
882  	
883  	   return SCIP_OKAY;
884  	}
885  	
886  	/*
887  	 * primal heuristic specific interface methods
888  	 */
889  	
890  	/** creates the zirounding primal heuristic and includes it in SCIP */
891  	SCIP_RETCODE SCIPincludeHeurZirounding(
892  	   SCIP*                 scip                /**< SCIP data structure */
893  	   )
894  	{
895  	   SCIP_HEURDATA* heurdata;
896  	   SCIP_HEUR* heur;
897  	
898  	   /* create zirounding primal heuristic data */
899  	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
900  	
901  	   /* include primal heuristic */
902  	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
903  	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
904  	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecZirounding, heurdata) );
905  	
906  	   assert(heur != NULL);
907  	
908  	   /* set non-NULL pointers to callback methods */
909  	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyZirounding) );
910  	   SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeZirounding) );
911  	   SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitZirounding) );
912  	   SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitZirounding) );
913  	   SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolZirounding) );
914  	
915  	   /* add zirounding primal heuristic parameters */
916  	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/zirounding/maxroundingloops",
917  	         "determines maximum number of rounding loops",
918  	         &heurdata->maxroundingloops, TRUE, DEFAULT_MAXROUNDINGLOOPS, -1, INT_MAX, NULL, NULL) );
919  	   SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/zirounding/stopziround",
920  	         "flag to determine if Zirounding is deactivated after a certain percentage of unsuccessful calls",
921  	         &heurdata->stopziround, TRUE, DEFAULT_STOPZIROUND, NULL, NULL) );
922  	   SCIP_CALL( SCIPaddRealParam(scip,"heuristics/zirounding/stoppercentage",
923  	         "if percentage of found solutions falls below this parameter, Zirounding will be deactivated",
924  	         &heurdata->stoppercentage, TRUE, DEFAULT_STOPPERCENTAGE, 0.0, 1.0, NULL, NULL) );
925  	   SCIP_CALL( SCIPaddIntParam(scip, "heuristics/zirounding/minstopncalls",
926  	         "determines the minimum number of calls before percentage-based deactivation of Zirounding is applied",
927  	         &heurdata->minstopncalls, TRUE, DEFAULT_MINSTOPNCALLS, 1, INT_MAX, NULL, NULL) );
928  	
929  	   return SCIP_OKAY;
930  	}
931