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_actconsdiving.c
26   	 * @ingroup DEFPLUGINS_HEUR
27   	 * @brief  LP diving heuristic that chooses fixings w.r.t. the active constraints the variable appear in
28   	 * @author Tobias Achterberg
29   	 */
30   	
31   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32   	
33   	#include "scip/heur_actconsdiving.h"
34   	#include "scip/heuristics.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_numerics.h"
45   	#include "scip/scip_prob.h"
46   	#include "scip/scip_sol.h"
47   	#include <string.h>
48   	
49   	#define HEUR_NAME             "actconsdiving"
50   	#define HEUR_DESC             "LP diving heuristic that chooses fixings w.r.t. the active constraints"
51   	#define HEUR_DISPCHAR         SCIP_HEURDISPCHAR_DIVING
52   	#define HEUR_PRIORITY         -1003700
53   	#define HEUR_FREQ             -1
54   	#define HEUR_FREQOFS          5
55   	#define HEUR_MAXDEPTH         -1
56   	#define HEUR_TIMING           SCIP_HEURTIMING_AFTERLPPLUNGE
57   	#define HEUR_USESSUBSCIP      FALSE  /**< does the heuristic use a secondary SCIP instance? */
58   	#define DIVESET_DIVETYPES     SCIP_DIVETYPE_INTEGRALITY /**< bit mask that represents all supported dive types */
59   	#define DIVESET_ISPUBLIC      TRUE   /**< is this dive set publicly available (ie., can be used by other primal heuristics?) */
60   	
61   	/*
62   	 * Default parameter settings
63   	 */
64   	
65   	#define DEFAULT_MINRELDEPTH         0.0 /**< minimal relative depth to start diving */
66   	#define DEFAULT_MAXRELDEPTH         1.0 /**< maximal relative depth to start diving */
67   	#define DEFAULT_MAXLPITERQUOT      0.05 /**< maximal fraction of diving LP iterations compared to node LP iterations */
68   	#define DEFAULT_MAXLPITEROFS       1000 /**< additional number of allowed LP iterations */
69   	#define DEFAULT_MAXDIVEUBQUOT       0.8 /**< maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound)
70   	                                         *   where diving is performed (0.0: no limit) */
71   	#define DEFAULT_MAXDIVEAVGQUOT      0.0 /**< maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound)
72   	                                         *   where diving is performed (0.0: no limit) */
73   	#define DEFAULT_MAXDIVEUBQUOTNOSOL  1.0 /**< maximal UBQUOT when no solution was found yet (0.0: no limit) */
74   	#define DEFAULT_MAXDIVEAVGQUOTNOSOL 1.0 /**< maximal AVGQUOT when no solution was found yet (0.0: no limit) */
75   	#define DEFAULT_BACKTRACK          TRUE /**< use one level of backtracking if infeasibility is encountered? */
76   	#define DEFAULT_LPRESOLVEDOMCHGQUOT 0.15 /**< percentage of immediate domain changes during probing to trigger LP resolve */
77   	#define DEFAULT_LPSOLVEFREQ           0 /**< LP solve frequency for diving heuristics */
78   	#define DEFAULT_ONLYLPBRANCHCANDS  TRUE /**< should only LP branching candidates be considered instead of the slower but
79   	                                         *   more general constraint handler diving variable selection? */
80   	#define DEFAULT_RANDSEED            149 /**< default random seed */
81   	
82   	/* locally defined heuristic data */
83   	struct SCIP_HeurData
84   	{
85   	   SCIP_SOL*             sol;                /**< working solution */
86   	};
87   	
88   	
89   	/*
90   	 * local methods
91   	 */
92   	
93   	/** returns a score value for the given variable based on the active constraints that the variable appears in */
94   	static
95   	SCIP_Real getNActiveConsScore(
96   	   SCIP*                 scip,               /**< SCIP data structure */
97   	   SCIP_SOL*             sol,                /**< working solution */
98   	   SCIP_VAR*             var,                /**< variable to get the score value for */
99   	   SCIP_Real*            downscore,          /**< pointer to store the score for branching downwards */
100  	   SCIP_Real*            upscore             /**< pointer to store the score for branching upwards */
101  	   )
102  	{
103  	   SCIP_COL* col;
104  	   SCIP_ROW** rows;
105  	   SCIP_Real* vals;
106  	   int nrows;
107  	   int r;
108  	   int nactrows;
109  	   SCIP_Real nlprows;
110  	   SCIP_Real downcoefsum;
111  	   SCIP_Real upcoefsum;
112  	   SCIP_Real score;
113  	
114  	   assert(downscore != NULL);
115  	   assert(upscore != NULL);
116  	
117  	   *downscore = 0.0;
118  	   *upscore = 0.0;
119  	   if( SCIPvarGetStatus(var) != SCIP_VARSTATUS_COLUMN )
120  	      return 0.0;
121  	
122  	   col = SCIPvarGetCol(var);
123  	   assert(col != NULL);
124  	
125  	   rows = SCIPcolGetRows(col);
126  	   vals = SCIPcolGetVals(col);
127  	   nrows = SCIPcolGetNLPNonz(col);
128  	   nactrows = 0;
129  	   downcoefsum = 0.0;
130  	   upcoefsum = 0.0;
131  	   for( r = 0; r < nrows; ++r )
132  	   {
133  	      SCIP_ROW* row;
134  	      SCIP_Real activity;
135  	      SCIP_Real lhs;
136  	      SCIP_Real rhs;
137  	      SCIP_Real dualsol;
138  	
139  	      row = rows[r];
140  	      /* calculate number of active constraint sides, i.e., count equations as two */
141  	      lhs = SCIProwGetLhs(row);
142  	      rhs = SCIProwGetRhs(row);
143  	
144  	      /* @todo this is suboptimal because activity is calculated by looping over all nonzeros of this row, need to
145  	       * store LP activities instead (which cannot be retrieved if no LP was solved at this node)
146  	       */
147  	      activity = SCIPgetRowSolActivity(scip, row, sol);
148  	
149  	      dualsol = SCIProwGetDualsol(row);
150  	      if( SCIPisFeasEQ(scip, activity, lhs) )
151  	      {
152  	         SCIP_Real coef;
153  	
154  	         nactrows++;
155  	         coef = vals[r] / SCIProwGetNorm(row);
156  	         if( SCIPisFeasPositive(scip, dualsol) )
157  	         {
158  	            if( coef > 0.0 )
159  	               downcoefsum += coef;
160  	            else
161  	               upcoefsum -= coef;
162  	         }
163  	      }
164  	      else if( SCIPisFeasEQ(scip, activity, rhs) )
165  	      {
166  	         SCIP_Real coef;
167  	
168  	         nactrows++;
169  	         coef = vals[r] / SCIProwGetNorm(row);
170  	         if( SCIPisFeasNegative(scip, dualsol) )
171  	         {
172  	            if( coef > 0.0 )
173  	               upcoefsum += coef;
174  	            else
175  	               downcoefsum -= coef;
176  	         }
177  	      }
178  	   }
179  	
180  	   /* use the number of LP rows for normalization */
181  	   nlprows = (SCIP_Real)SCIPgetNLPRows(scip);
182  	   upcoefsum /= nlprows;
183  	   downcoefsum /= nlprows;
184  	
185  	   /* calculate the score using SCIP's branch score. Pass NULL as variable to not have the var branch factor influence
186  	    * the result
187  	    */
188  	   score = nactrows / nlprows + SCIPgetBranchScore(scip, NULL, downcoefsum, upcoefsum);
189  	
190  	   assert(score <= 3.0);
191  	   assert(score >= 0.0);
192  	
193  	   *downscore = downcoefsum;
194  	   *upscore = upcoefsum;
195  	
196  	   return score;
197  	}
198  	
199  	
200  	/*
201  	 * Callback methods
202  	 */
203  	
204  	/** copy method for primal heuristic plugins (called when SCIP copies plugins) */
205  	static
206  	SCIP_DECL_HEURCOPY(heurCopyActconsdiving)
207  	{  /*lint --e{715}*/
208  	   assert(scip != NULL);
209  	   assert(heur != NULL);
210  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
211  	
212  	   /* call inclusion method of primal heuristic */
213  	   SCIP_CALL( SCIPincludeHeurActconsdiving(scip) );
214  	
215  	   return SCIP_OKAY;
216  	}
217  	
218  	/** destructor of primal heuristic to free user data (called when SCIP is exiting) */
219  	static
220  	SCIP_DECL_HEURFREE(heurFreeActconsdiving) /*lint --e{715}*/
221  	{  /*lint --e{715}*/
222  	   SCIP_HEURDATA* heurdata;
223  	
224  	   assert(heur != NULL);
225  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
226  	   assert(scip != NULL);
227  	
228  	   /* free heuristic data */
229  	   heurdata = SCIPheurGetData(heur);
230  	   assert(heurdata != NULL);
231  	
232  	   SCIPfreeBlockMemory(scip, &heurdata);
233  	   SCIPheurSetData(heur, NULL);
234  	
235  	   return SCIP_OKAY;
236  	}
237  	
238  	
239  	/** initialization method of primal heuristic (called after problem was transformed) */
240  	static
241  	SCIP_DECL_HEURINIT(heurInitActconsdiving) /*lint --e{715}*/
242  	{  /*lint --e{715}*/
243  	   SCIP_HEURDATA* heurdata;
244  	
245  	   assert(heur != NULL);
246  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
247  	
248  	   /* get heuristic data */
249  	   heurdata = SCIPheurGetData(heur);
250  	   assert(heurdata != NULL);
251  	
252  	   /* create working solution */
253  	   SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
254  	
255  	   return SCIP_OKAY;
256  	}
257  	
258  	
259  	/** deinitialization method of primal heuristic (called before transformed problem is freed) */
260  	static
261  	SCIP_DECL_HEUREXIT(heurExitActconsdiving) /*lint --e{715}*/
262  	{  /*lint --e{715}*/
263  	   SCIP_HEURDATA* heurdata;
264  	
265  	   assert(heur != NULL);
266  	   assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
267  	
268  	   /* get heuristic data */
269  	   heurdata = SCIPheurGetData(heur);
270  	   assert(heurdata != NULL);
271  	
272  	   /* free working solution */
273  	   SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
274  	
275  	   return SCIP_OKAY;
276  	}
277  	
278  	
279  	/** execution method of primal heuristic */
280  	static
281  	SCIP_DECL_HEUREXEC(heurExecActconsdiving) /*lint --e{715}*/
282  	{  /*lint --e{715}*/
283  	   SCIP_HEURDATA* heurdata;
284  	   SCIP_DIVESET* diveset;
285  	
286  	   heurdata = SCIPheurGetData(heur);
287  	
288  	   assert(SCIPheurGetNDivesets(heur) > 0);
289  	   assert(SCIPheurGetDivesets(heur) != NULL);
290  	   diveset = SCIPheurGetDivesets(heur)[0];
291  	   assert(diveset != NULL);
292  	
293  	   *result = SCIP_DIDNOTRUN;
294  	
295  	   /* if there are no integer variables, stop execution */
296  	   if( SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip) == 0 )
297  	      return SCIP_OKAY;
298  	
299  	   SCIP_CALL( SCIPperformGenericDivingAlgorithm(scip, diveset, heurdata->sol, heur, result, nodeinfeasible, -1L, -1, -1.0, SCIP_DIVECONTEXT_SINGLE) );
300  	
301  	   return SCIP_OKAY;
302  	}
303  	
304  	/** calculate score and preferred rounding direction for the candidate variable; the best candidate maximizes the
305  	 *  score
306  	 */
307  	static
308  	SCIP_DECL_DIVESETGETSCORE(divesetGetScoreActconsdiving)
309  	{
310  	   SCIP_Bool mayrounddown;
311  	   SCIP_Bool mayroundup;
312  	   SCIP_Real downscore;
313  	   SCIP_Real upscore;
314  	
315  	   mayrounddown = SCIPvarMayRoundDown(cand);
316  	   mayroundup = SCIPvarMayRoundUp(cand);
317  	
318  	   /* first, calculate the variable score */
319  	   assert(SCIPdivesetGetWorkSolution(diveset) != NULL);
320  	   *score = getNActiveConsScore(scip, SCIPdivesetGetWorkSolution(diveset), cand, &downscore, &upscore);
321  	
322  	   /* get the rounding direction: prefer an unroundable direction */
323  	   if( mayrounddown && mayroundup )
324  	   {
325  	      /* try to avoid variability; decide randomly if the LP solution can contain some noise */
326  	      if( SCIPisEQ(scip, candsfrac, 0.5) )
327  	         *roundup = (SCIPrandomGetInt(SCIPdivesetGetRandnumgen(diveset), 0, 1) == 0);
328  	      else
329  	         *roundup = (candsfrac > 0.5);
330  	   }
331  	   else if( mayrounddown || mayroundup )
332  	      *roundup = mayrounddown;
333  	   else
334  	      *roundup = (downscore > upscore);
335  	
336  	   if( *roundup )
337  	      candsfrac = 1.0 - candsfrac;
338  	
339  	   /* penalize too small fractions */
340  	   if( SCIPisEQ(scip, candsfrac, 0.01) )
341  	   {
342  	      /* try to avoid variability; decide randomly if the LP solution can contain some noise.
343  	       * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for scaling the score
344  	       */
345  	      if( SCIPrandomGetInt(SCIPdivesetGetRandnumgen(diveset), 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
346  	         (*score) *= 0.01;
347  	   }
348  	   else if( candsfrac < 0.01 )
349  	      (*score) *= 0.01;
350  	
351  	   /* prefer decisions on binary variables */
352  	   if( !SCIPvarIsBinary(cand) )
353  	      (*score) *= 0.01;
354  	
355  	   /* penalize variable if it may be rounded */
356  	   if( mayrounddown || mayroundup )
357  	      *score -= 3.0;
358  	
359  	   assert(!(mayrounddown || mayroundup) || *score <= 0.0);
360  	
361  	   return SCIP_OKAY;
362  	}
363  	
364  	#define divesetAvailableActconsdiving NULL
365  	
366  	/*
367  	 * heuristic specific interface methods
368  	 */
369  	
370  	/** creates the actconsdiving heuristic and includes it in SCIP */
371  	SCIP_RETCODE SCIPincludeHeurActconsdiving(
372  	   SCIP*                 scip                /**< SCIP data structure */
373  	   )
374  	{
375  	   SCIP_HEURDATA* heurdata;
376  	   SCIP_HEUR* heur;
377  	
378  	   /* create actconsdiving primal heuristic data */
379  	   SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
380  	
381  	   /* include primal heuristic */
382  	   SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
383  	         HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
384  	         HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecActconsdiving, heurdata) );
385  	
386  	   assert(heur != NULL);
387  	
388  	   /* set non-NULL pointers to callback methods */
389  	   SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyActconsdiving) );
390  	   SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeActconsdiving) );
391  	   SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitActconsdiving) );
392  	   SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitActconsdiving) );
393  	
394  	   /* create a diveset (this will automatically install some additional parameters for the heuristic)*/
395  	   SCIP_CALL( SCIPcreateDiveset(scip, NULL, heur, HEUR_NAME, DEFAULT_MINRELDEPTH, DEFAULT_MAXRELDEPTH, DEFAULT_MAXLPITERQUOT,
396  	         DEFAULT_MAXDIVEUBQUOT, DEFAULT_MAXDIVEAVGQUOT, DEFAULT_MAXDIVEUBQUOTNOSOL, DEFAULT_MAXDIVEAVGQUOTNOSOL, DEFAULT_LPRESOLVEDOMCHGQUOT,
397  	         DEFAULT_LPSOLVEFREQ, DEFAULT_MAXLPITEROFS, DEFAULT_RANDSEED,
398  	         DEFAULT_BACKTRACK, DEFAULT_ONLYLPBRANCHCANDS,
399  	         DIVESET_ISPUBLIC, DIVESET_DIVETYPES, divesetGetScoreActconsdiving, divesetAvailableActconsdiving) );
400  	
401  	   return SCIP_OKAY;
402  	}
403  	
404