Scippy

SCIP

Solving Constraint Integer Programs

heur_pscostdiving.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file heur_pscostdiving.c
17  * @brief LP diving heuristic that chooses fixings w.r.t. the pseudo cost values
18  * @author Tobias Achterberg
19  */
20 
21 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
22 
23 #include <assert.h>
24 #include <string.h>
25 
26 #include "scip/heur_pscostdiving.h"
27 
28 
29 #define HEUR_NAME "pscostdiving"
30 #define HEUR_DESC "LP diving heuristic that chooses fixings w.r.t. the pseudo cost values"
31 #define HEUR_DISPCHAR 'p'
32 #define HEUR_PRIORITY -1002000
33 #define HEUR_FREQ 10
34 #define HEUR_FREQOFS 2
35 #define HEUR_MAXDEPTH -1
36 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPPLUNGE
37 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
38 
39 
40 /*
41  * Default parameter settings
42  */
43 
44 #define DEFAULT_MINRELDEPTH 0.0 /**< minimal relative depth to start diving */
45 #define DEFAULT_MAXRELDEPTH 1.0 /**< maximal relative depth to start diving */
46 #define DEFAULT_MAXLPITERQUOT 0.05 /**< maximal fraction of diving LP iterations compared to node LP iterations */
47 #define DEFAULT_MAXLPITEROFS 1000 /**< additional number of allowed LP iterations */
48 #define DEFAULT_MAXDIVEUBQUOT 0.8 /**< maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound)
49  * where diving is performed (0.0: no limit) */
50 #define DEFAULT_MAXDIVEAVGQUOT 0.0 /**< maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound)
51  * where diving is performed (0.0: no limit) */
52 #define DEFAULT_MAXDIVEUBQUOTNOSOL 0.1 /**< maximal UBQUOT when no solution was found yet (0.0: no limit) */
53 #define DEFAULT_MAXDIVEAVGQUOTNOSOL 0.0 /**< maximal AVGQUOT when no solution was found yet (0.0: no limit) */
54 #define DEFAULT_BACKTRACK TRUE /**< use one level of backtracking if infeasibility is encountered? */
55 
56 #define MINLPITER 10000 /**< minimal number of LP iterations allowed in each LP solving call */
57 
58 
59 /* locally defined heuristic data */
60 struct SCIP_HeurData
61 {
62  SCIP_SOL* sol; /**< working solution */
63  SCIP_Real minreldepth; /**< minimal relative depth to start diving */
64  SCIP_Real maxreldepth; /**< maximal relative depth to start diving */
65  SCIP_Real maxlpiterquot; /**< maximal fraction of diving LP iterations compared to node LP iterations */
66  int maxlpiterofs; /**< additional number of allowed LP iterations */
67  SCIP_Real maxdiveubquot; /**< maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound)
68  * where diving is performed (0.0: no limit) */
69  SCIP_Real maxdiveavgquot; /**< maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound)
70  * where diving is performed (0.0: no limit) */
71  SCIP_Real maxdiveubquotnosol; /**< maximal UBQUOT when no solution was found yet (0.0: no limit) */
72  SCIP_Real maxdiveavgquotnosol;/**< maximal AVGQUOT when no solution was found yet (0.0: no limit) */
73  SCIP_Bool backtrack; /**< use one level of backtracking if infeasibility is encountered? */
74  SCIP_Longint nlpiterations; /**< LP iterations used in this heuristic */
75  int nsuccess; /**< number of runs that produced at least one feasible solution */
76 };
77 
78 
79 /*
80  * local methods
81  */
82 
83 /** calculates the pseudocost score for a given variable w.r.t. a given solution value and a given rounding direction */
84 static
85 void calcPscostQuot(
86  SCIP* scip, /**< SCIP data structure */
87  SCIP_VAR* var, /**< problem variable */
88  SCIP_Real primsol, /**< primal solution of variable */
89  SCIP_Real frac, /**< fractionality of variable */
90  int rounddir, /**< -1: round down, +1: round up, 0: select due to pseudo cost values */
91  SCIP_Real* pscostquot, /**< pointer to store pseudo cost quotient */
92  SCIP_Bool* roundup /**< pointer to store whether the variable should be rounded up */
93  )
94 {
95  SCIP_Real pscostdown;
96  SCIP_Real pscostup;
97 
98  assert(pscostquot != NULL);
99  assert(roundup != NULL);
100  assert(SCIPisEQ(scip, frac, primsol - SCIPfeasFloor(scip, primsol)));
101 
102  /* bound fractions to not prefer variables that are nearly integral */
103  frac = MAX(frac, 0.1);
104  frac = MIN(frac, 0.9);
105 
106  /* get pseudo cost quotient */
107  pscostdown = SCIPgetVarPseudocostVal(scip, var, 0.0-frac);
108  pscostup = SCIPgetVarPseudocostVal(scip, var, 1.0-frac);
109  assert(pscostdown >= 0.0 && pscostup >= 0.0);
110 
111  /* choose rounding direction */
112  if( rounddir == -1 )
113  *roundup = FALSE;
114  else if( rounddir == +1 )
115  *roundup = TRUE;
116  else if( primsol < SCIPvarGetRootSol(var) - 0.4 )
117  *roundup = FALSE;
118  else if( primsol > SCIPvarGetRootSol(var) + 0.4 )
119  *roundup = TRUE;
120  else if( frac < 0.3 )
121  *roundup = FALSE;
122  else if( frac > 0.7 )
123  *roundup = TRUE;
124  else if( pscostdown < pscostup )
125  *roundup = FALSE;
126  else
127  *roundup = TRUE;
128 
129  /* calculate pseudo cost quotient */
130  if( *roundup )
131  *pscostquot = sqrt(frac) * (1.0+pscostdown) / (1.0+pscostup);
132  else
133  *pscostquot = sqrt(1.0-frac) * (1.0+pscostup) / (1.0+pscostdown);
134 
135  /* prefer decisions on binary variables */
136  if( SCIPvarIsBinary(var) )
137  (*pscostquot) *= 1000.0;
138 }
139 
140 
141 /*
142  * Callback methods
143  */
144 
145 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
146 static
147 SCIP_DECL_HEURCOPY(heurCopyPscostdiving)
148 { /*lint --e{715}*/
149  assert(scip != NULL);
150  assert(heur != NULL);
151  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
152 
153  /* call inclusion method of primal heuristic */
154  SCIP_CALL( SCIPincludeHeurPscostdiving(scip) );
155 
156  return SCIP_OKAY;
157 }
158 
159 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
160 static
161 SCIP_DECL_HEURFREE(heurFreePscostdiving) /*lint --e{715}*/
162 { /*lint --e{715}*/
163  SCIP_HEURDATA* heurdata;
164 
165  assert(heur != NULL);
166  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
167  assert(scip != NULL);
168 
169  /* free heuristic data */
170  heurdata = SCIPheurGetData(heur);
171  assert(heurdata != NULL);
172  SCIPfreeMemory(scip, &heurdata);
173  SCIPheurSetData(heur, NULL);
174 
175  return SCIP_OKAY;
176 }
177 
178 
179 /** initialization method of primal heuristic (called after problem was transformed) */
180 static
181 SCIP_DECL_HEURINIT(heurInitPscostdiving) /*lint --e{715}*/
182 { /*lint --e{715}*/
183  SCIP_HEURDATA* heurdata;
184 
185  assert(heur != NULL);
186  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
187 
188  /* get heuristic data */
189  heurdata = SCIPheurGetData(heur);
190  assert(heurdata != NULL);
191 
192  /* create working solution */
193  SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
194 
195  /* initialize data */
196  heurdata->nlpiterations = 0;
197  heurdata->nsuccess = 0;
198 
199  return SCIP_OKAY;
200 }
201 
202 
203 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
204 static
205 SCIP_DECL_HEUREXIT(heurExitPscostdiving) /*lint --e{715}*/
206 { /*lint --e{715}*/
207  SCIP_HEURDATA* heurdata;
208 
209  assert(heur != NULL);
210  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
211 
212  /* get heuristic data */
213  heurdata = SCIPheurGetData(heur);
214  assert(heurdata != NULL);
215 
216  /* free working solution */
217  SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
218 
219  return SCIP_OKAY;
220 }
221 
222 
223 /** execution method of primal heuristic */
224 static
225 SCIP_DECL_HEUREXEC(heurExecPscostdiving) /*lint --e{715}*/
226 { /*lint --e{715}*/
227  SCIP_HEURDATA* heurdata;
228  SCIP_LPSOLSTAT lpsolstat;
229  SCIP_VAR* var;
230  SCIP_VAR** lpcands;
231  SCIP_Real* lpcandssol;
232  SCIP_Real* lpcandsfrac;
233  SCIP_Real searchubbound;
234  SCIP_Real searchavgbound;
235  SCIP_Real searchbound;
236  SCIP_Real objval;
237  SCIP_Real oldobjval;
238  SCIP_Real primsol;
239  SCIP_Real frac;
240  SCIP_Real pscostquot;
241  SCIP_Real bestpscostquot;
242  SCIP_Bool bestcandmayrounddown;
243  SCIP_Bool bestcandmayroundup;
244  SCIP_Bool bestcandroundup;
245  SCIP_Bool mayrounddown;
246  SCIP_Bool mayroundup;
247  SCIP_Bool roundup;
248  SCIP_Bool lperror;
249  SCIP_Bool cutoff;
250  SCIP_Bool backtracked;
251  SCIP_Longint ncalls;
252  SCIP_Longint nsolsfound;
253  SCIP_Longint nlpiterations;
254  SCIP_Longint maxnlpiterations;
255  int nlpcands;
256  int startnlpcands;
257  int depth;
258  int maxdepth;
259  int maxdivedepth;
260  int divedepth;
261  int bestcand;
262  int c;
263 
264  assert(heur != NULL);
265  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
266  assert(scip != NULL);
267  assert(result != NULL);
268  assert(SCIPhasCurrentNodeLP(scip));
269 
270  *result = SCIP_DELAYED;
271 
272  /* do not call heuristic of node was already detected to be infeasible */
273  if( nodeinfeasible )
274  return SCIP_OKAY;
275 
276  /* only call heuristic, if an optimal LP solution is at hand */
277  if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
278  return SCIP_OKAY;
279 
280  /* only call heuristic, if the LP objective value is smaller than the cutoff bound */
281  if( SCIPisGE(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)) )
282  return SCIP_OKAY;
283 
284  /* only call heuristic, if the LP solution is basic (which allows fast resolve in diving) */
285  if( !SCIPisLPSolBasic(scip) )
286  return SCIP_OKAY;
287 
288  /* don't dive two times at the same node */
289  if( SCIPgetLastDivenode(scip) == SCIPgetNNodes(scip) && SCIPgetDepth(scip) > 0 )
290  return SCIP_OKAY;
291 
292  *result = SCIP_DIDNOTRUN;
293 
294  /* get heuristic's data */
295  heurdata = SCIPheurGetData(heur);
296  assert(heurdata != NULL);
297 
298  /* only try to dive, if we are in the correct part of the tree, given by minreldepth and maxreldepth */
299  depth = SCIPgetDepth(scip);
300  maxdepth = SCIPgetMaxDepth(scip);
301  maxdepth = MAX(maxdepth, 30);
302  if( depth < heurdata->minreldepth*maxdepth || depth > heurdata->maxreldepth*maxdepth )
303  return SCIP_OKAY;
304 
305  /* calculate the maximal number of LP iterations until heuristic is aborted */
306  nlpiterations = SCIPgetNNodeLPIterations(scip);
307  ncalls = SCIPheurGetNCalls(heur);
308  nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + heurdata->nsuccess;
309  maxnlpiterations = (SCIP_Longint)((1.0 + 10.0*(nsolsfound+1.0)/(ncalls+1.0)) * heurdata->maxlpiterquot * nlpiterations);
310  maxnlpiterations += heurdata->maxlpiterofs;
311 
312  /* don't try to dive, if we took too many LP iterations during diving */
313  if( heurdata->nlpiterations >= maxnlpiterations )
314  return SCIP_OKAY;
315 
316  /* allow at least a certain number of LP iterations in this dive */
317  maxnlpiterations = MAX(maxnlpiterations, heurdata->nlpiterations + MINLPITER);
318 
319  /* get fractional variables that should be integral */
320  SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, &lpcandsfrac, &nlpcands, NULL, NULL) );
321 
322  /* don't try to dive, if there are no fractional variables */
323  if( nlpcands == 0 )
324  return SCIP_OKAY;
325 
326  /* calculate the objective search bound */
327  if( SCIPgetNSolsFound(scip) == 0 )
328  {
329  if( heurdata->maxdiveubquotnosol > 0.0 )
330  searchubbound = SCIPgetLowerbound(scip)
331  + heurdata->maxdiveubquotnosol * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip));
332  else
333  searchubbound = SCIPinfinity(scip);
334  if( heurdata->maxdiveavgquotnosol > 0.0 )
335  searchavgbound = SCIPgetLowerbound(scip)
336  + heurdata->maxdiveavgquotnosol * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip));
337  else
338  searchavgbound = SCIPinfinity(scip);
339  }
340  else
341  {
342  if( heurdata->maxdiveubquot > 0.0 )
343  searchubbound = SCIPgetLowerbound(scip)
344  + heurdata->maxdiveubquot * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip));
345  else
346  searchubbound = SCIPinfinity(scip);
347  if( heurdata->maxdiveavgquot > 0.0 )
348  searchavgbound = SCIPgetLowerbound(scip)
349  + heurdata->maxdiveavgquot * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip));
350  else
351  searchavgbound = SCIPinfinity(scip);
352  }
353  searchbound = MIN(searchubbound, searchavgbound);
354  if( SCIPisObjIntegral(scip) )
355  searchbound = SCIPceil(scip, searchbound);
356 
357  /* calculate the maximal diving depth: 10 * min{number of integer variables, max depth} */
358  maxdivedepth = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
359  maxdivedepth = MIN(maxdivedepth, maxdepth);
360  maxdivedepth *= 10;
361 
362 
363  *result = SCIP_DIDNOTFIND;
364 
365  /* start diving */
366  SCIP_CALL( SCIPstartProbing(scip) );
367 
368  /* enables collection of variable statistics during probing */
369  SCIPenableVarHistory(scip);
370 
371  /* get LP objective value */
372  lpsolstat = SCIP_LPSOLSTAT_OPTIMAL;
373  objval = SCIPgetLPObjval(scip);
374 
375  SCIPdebugMessage("(node %"SCIP_LONGINT_FORMAT") executing pscostdiving heuristic: depth=%d, %d fractionals, dualbound=%g, searchbound=%g\n",
376  SCIPgetNNodes(scip), SCIPgetDepth(scip), nlpcands, SCIPgetDualbound(scip), SCIPretransformObj(scip, searchbound));
377 
378  /* dive as long we are in the given objective, depth and iteration limits and fractional variables exist, but
379  * - if possible, we dive at least with the depth 10
380  * - if the number of fractional variables decreased at least with 1 variable per 2 dive depths, we continue diving
381  */
382  lperror = FALSE;
383  cutoff = FALSE;
384  divedepth = 0;
385  bestcandmayrounddown = FALSE;
386  bestcandmayroundup = FALSE;
387  startnlpcands = nlpcands;
388  while( !lperror && !cutoff && lpsolstat == SCIP_LPSOLSTAT_OPTIMAL && nlpcands > 0
389  && (divedepth < 10
390  || nlpcands <= startnlpcands - divedepth/2
391  || (divedepth < maxdivedepth && heurdata->nlpiterations < maxnlpiterations && objval < searchbound))
392  && !SCIPisStopped(scip) )
393  {
394  SCIP_CALL( SCIPnewProbingNode(scip) );
395  divedepth++;
396 
397  /* choose variable fixing:
398  * - prefer variables that may not be rounded without destroying LP feasibility:
399  * - of these variables, round variable with largest rel. difference of pseudo cost values in corresponding
400  * direction
401  * - if all remaining fractional variables may be rounded without destroying LP feasibility:
402  * - round variable in the objective value direction
403  */
404  bestcand = -1;
405  bestpscostquot = -1.0;
406  bestcandmayrounddown = TRUE;
407  bestcandmayroundup = TRUE;
408  bestcandroundup = FALSE;
409  for( c = 0; c < nlpcands; ++c )
410  {
411  var = lpcands[c];
412  mayrounddown = SCIPvarMayRoundDown(var);
413  mayroundup = SCIPvarMayRoundUp(var);
414  primsol = lpcandssol[c];
415  frac = lpcandsfrac[c];
416  if( mayrounddown || mayroundup )
417  {
418  /* the candidate may be rounded: choose this candidate only, if the best candidate may also be rounded */
419  if( bestcandmayrounddown || bestcandmayroundup )
420  {
421  /* choose rounding direction:
422  * - if variable may be rounded in both directions, round corresponding to the pseudo cost values
423  * - otherwise, round in the infeasible direction, because feasible direction is tried by rounding
424  * the current fractional solution
425  */
426  roundup = FALSE;
427  if( mayrounddown && mayroundup )
428  calcPscostQuot(scip, var, primsol, frac, 0, &pscostquot, &roundup);
429  else if( mayrounddown )
430  calcPscostQuot(scip, var, primsol, frac, +1, &pscostquot, &roundup);
431  else
432  calcPscostQuot(scip, var, primsol, frac, -1, &pscostquot, &roundup);
433 
434  /* check, if candidate is new best candidate */
435  if( pscostquot > bestpscostquot )
436  {
437  bestcand = c;
438  bestpscostquot = pscostquot;
439  bestcandmayrounddown = mayrounddown;
440  bestcandmayroundup = mayroundup;
441  bestcandroundup = roundup;
442  }
443  }
444  }
445  else
446  {
447  /* the candidate may not be rounded: calculate pseudo cost quotient and preferred direction */
448  calcPscostQuot(scip, var, primsol, frac, 0, &pscostquot, &roundup);
449 
450  /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
451  if( bestcandmayrounddown || bestcandmayroundup || pscostquot > bestpscostquot )
452  {
453  bestcand = c;
454  bestpscostquot = pscostquot;
455  bestcandmayrounddown = FALSE;
456  bestcandmayroundup = FALSE;
457  bestcandroundup = roundup;
458  }
459  }
460  }
461  assert(bestcand != -1);
462 
463  /* if all candidates are roundable, try to round the solution */
464  if( bestcandmayrounddown || bestcandmayroundup )
465  {
466  SCIP_Bool success;
467 
468  /* create solution from diving LP and try to round it */
469  SCIP_CALL( SCIPlinkLPSol(scip, heurdata->sol) );
470  SCIP_CALL( SCIProundSol(scip, heurdata->sol, &success) );
471 
472  if( success )
473  {
474  SCIPdebugMessage("pscostdiving found roundable primal solution: obj=%g\n", SCIPgetSolOrigObj(scip, heurdata->sol));
475 
476  /* try to add solution to SCIP */
477  SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, &success) );
478 
479  /* check, if solution was feasible and good enough */
480  if( success )
481  {
482  SCIPdebugMessage(" -> solution was feasible and good enough\n");
483  *result = SCIP_FOUNDSOL;
484  }
485  }
486  }
487 
488  var = lpcands[bestcand];
489 
490  backtracked = FALSE;
491  do
492  {
493  /* if the variable is already fixed or if the solution value is outside the domain, numerical troubles may have
494  * occured or variable was fixed by propagation while backtracking => Abort diving!
495  */
496  if( SCIPvarGetLbLocal(var) >= SCIPvarGetUbLocal(var) - 0.5 )
497  {
498  SCIPdebugMessage("Selected variable <%s> already fixed to [%g,%g] (solval: %.9f), diving aborted \n",
499  SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), lpcandssol[bestcand]);
500  cutoff = TRUE;
501  break;
502  }
503  if( SCIPisFeasLT(scip, lpcandssol[bestcand], SCIPvarGetLbLocal(var)) || SCIPisFeasGT(scip, lpcandssol[bestcand], SCIPvarGetUbLocal(var)) )
504  {
505  SCIPdebugMessage("selected variable's <%s> solution value is outside the domain [%g,%g] (solval: %.9f), diving aborted\n",
506  SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), lpcandssol[bestcand]);
507  assert(backtracked);
508  break;
509  }
510 
511  /* apply rounding of best candidate */
512  if( bestcandroundup == !backtracked )
513  {
514  /* round variable up */
515  SCIPdebugMessage(" dive %d/%d, LP iter %"SCIP_LONGINT_FORMAT"/%"SCIP_LONGINT_FORMAT": var <%s>, round=%u/%u, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
516  divedepth, maxdivedepth, heurdata->nlpiterations, maxnlpiterations,
517  SCIPvarGetName(var), bestcandmayrounddown, bestcandmayroundup,
518  lpcandssol[bestcand], SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var),
519  SCIPfeasCeil(scip, lpcandssol[bestcand]), SCIPvarGetUbLocal(var));
520  SCIP_CALL( SCIPchgVarLbProbing(scip, var, SCIPfeasCeil(scip, lpcandssol[bestcand])) );
521  }
522  else
523  {
524  /* round variable down */
525  SCIPdebugMessage(" dive %d/%d, LP iter %"SCIP_LONGINT_FORMAT"/%"SCIP_LONGINT_FORMAT": var <%s>, round=%u/%u, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
526  divedepth, maxdivedepth, heurdata->nlpiterations, maxnlpiterations,
527  SCIPvarGetName(var), bestcandmayrounddown, bestcandmayroundup,
528  lpcandssol[bestcand], SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var),
529  SCIPvarGetLbLocal(var), SCIPfeasFloor(scip, lpcandssol[bestcand]));
530  SCIP_CALL( SCIPchgVarUbProbing(scip, lpcands[bestcand], SCIPfeasFloor(scip, lpcandssol[bestcand])) );
531  }
532 
533  /* apply domain propagation */
534  SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, NULL) );
535  if( !cutoff )
536  {
537  /* resolve the diving LP */
538  /* Errors in the LP solver should not kill the overall solving process, if the LP is just needed for a heuristic.
539  * Hence in optimized mode, the return code is caught and a warning is printed, only in debug mode, SCIP will stop.
540  */
541 #ifdef NDEBUG
542  SCIP_RETCODE retstat;
543  nlpiterations = SCIPgetNLPIterations(scip);
544  retstat = SCIPsolveProbingLP(scip, MAX((int)(maxnlpiterations - heurdata->nlpiterations), MINLPITER), &lperror, &cutoff);
545  if( retstat != SCIP_OKAY )
546  {
547  SCIPwarningMessage(scip, "Error while solving LP in Pscostdiving heuristic; LP solve terminated with code <%d>\n",retstat);
548  }
549 #else
550  nlpiterations = SCIPgetNLPIterations(scip);
551  SCIP_CALL( SCIPsolveProbingLP(scip, MAX((int)(maxnlpiterations - heurdata->nlpiterations), MINLPITER), &lperror, &cutoff) );
552 #endif
553 
554  if( lperror )
555  break;
556 
557  /* update iteration count */
558  heurdata->nlpiterations += SCIPgetNLPIterations(scip) - nlpiterations;
559 
560  /* get LP solution status, objective value, and fractional variables, that should be integral */
561  lpsolstat = SCIPgetLPSolstat(scip);
562  assert(cutoff || (lpsolstat != SCIP_LPSOLSTAT_OBJLIMIT && lpsolstat != SCIP_LPSOLSTAT_INFEASIBLE &&
563  (lpsolstat != SCIP_LPSOLSTAT_OPTIMAL || SCIPisLT(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)))));
564  }
565 
566  /* perform backtracking if a cutoff was detected */
567  if( cutoff && !backtracked && heurdata->backtrack )
568  {
569  SCIPdebugMessage(" *** cutoff detected at level %d - backtracking\n", SCIPgetProbingDepth(scip));
570  SCIP_CALL( SCIPbacktrackProbing(scip, SCIPgetProbingDepth(scip)-1) );
571  SCIP_CALL( SCIPnewProbingNode(scip) );
572  backtracked = TRUE;
573  }
574  else
575  backtracked = FALSE;
576  }
577  while( backtracked );
578 
579  if( !lperror && !cutoff && lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
580  {
581  /* get new objective value */
582  oldobjval = objval;
583  objval = SCIPgetLPObjval(scip);
584 
585  /* update pseudo cost values */
586  if( SCIPisGT(scip, objval, oldobjval) )
587  {
588  if( bestcandroundup )
589  {
590  SCIP_CALL( SCIPupdateVarPseudocost(scip, lpcands[bestcand], 1.0-lpcandsfrac[bestcand],
591  objval - oldobjval, 1.0) );
592  }
593  else
594  {
595  SCIP_CALL( SCIPupdateVarPseudocost(scip, lpcands[bestcand], 0.0-lpcandsfrac[bestcand],
596  objval - oldobjval, 1.0) );
597  }
598  }
599 
600  /* get new fractional variables */
601  SCIP_CALL( SCIPgetLPBranchCands(scip, &lpcands, &lpcandssol, &lpcandsfrac, &nlpcands, NULL, NULL) );
602  }
603  SCIPdebugMessage(" -> lpsolstat=%d, objval=%g/%g, nfrac=%d, lpiterations=%"SCIP_LONGINT_FORMAT"/%"SCIP_LONGINT_FORMAT"\n",
604  lpsolstat, objval, searchbound, nlpcands, heurdata->nlpiterations, maxnlpiterations);
605  }
606 
607  /* check if a solution has been found */
608  if( nlpcands == 0 && !lperror && !cutoff && lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
609  {
610  SCIP_Bool success;
611 
612  /* create solution from diving LP */
613  SCIP_CALL( SCIPlinkLPSol(scip, heurdata->sol) );
614  SCIPdebugMessage("pscostdiving found primal solution: obj=%g\n", SCIPgetSolOrigObj(scip, heurdata->sol));
615 
616  /* try to add solution to SCIP */
617  SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, &success) );
618 
619  /* check, if solution was feasible and good enough */
620  if( success )
621  {
622  SCIPdebugMessage(" -> solution was feasible and good enough\n");
623  *result = SCIP_FOUNDSOL;
624  }
625  }
626 
627  /* end diving */
628  SCIP_CALL( SCIPendProbing(scip) );
629 
630  if( *result == SCIP_FOUNDSOL )
631  heurdata->nsuccess++;
632 
633  SCIPdebugMessage("pscostdiving heuristic finished\n");
634 
635  return SCIP_OKAY;
636 }
637 
638 
639 /*
640  * heuristic specific interface methods
641  */
642 
643 /** creates the pscostdiving heuristic and includes it in SCIP */
644 SCIP_RETCODE SCIPincludeHeurPscostdiving(
645  SCIP* scip /**< SCIP data structure */
646  )
647 {
648  SCIP_HEURDATA* heurdata;
649  SCIP_HEUR* heur;
650 
651  /* create Pscostdiving primal heuristic data */
652  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
653 
654  /* include primal heuristic */
655  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
656  HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
657  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecPscostdiving, heurdata) );
658 
659  assert(heur != NULL);
660 
661  /* set non-NULL pointers to callback methods */
662  SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyPscostdiving) );
663  SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreePscostdiving) );
664  SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitPscostdiving) );
665  SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitPscostdiving) );
666 
667  /* pscostdiving heuristic parameters */
668  SCIP_CALL( SCIPaddRealParam(scip,
669  "heuristics/pscostdiving/minreldepth",
670  "minimal relative depth to start diving",
671  &heurdata->minreldepth, TRUE, DEFAULT_MINRELDEPTH, 0.0, 1.0, NULL, NULL) );
672  SCIP_CALL( SCIPaddRealParam(scip,
673  "heuristics/pscostdiving/maxreldepth",
674  "maximal relative depth to start diving",
675  &heurdata->maxreldepth, TRUE, DEFAULT_MAXRELDEPTH, 0.0, 1.0, NULL, NULL) );
676  SCIP_CALL( SCIPaddRealParam(scip,
677  "heuristics/pscostdiving/maxlpiterquot",
678  "maximal fraction of diving LP iterations compared to node LP iterations",
679  &heurdata->maxlpiterquot, FALSE, DEFAULT_MAXLPITERQUOT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
680  SCIP_CALL( SCIPaddIntParam(scip,
681  "heuristics/pscostdiving/maxlpiterofs",
682  "additional number of allowed LP iterations",
683  &heurdata->maxlpiterofs, FALSE, DEFAULT_MAXLPITEROFS, 0, INT_MAX, NULL, NULL) );
684  SCIP_CALL( SCIPaddRealParam(scip,
685  "heuristics/pscostdiving/maxdiveubquot",
686  "maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound) where diving is performed (0.0: no limit)",
687  &heurdata->maxdiveubquot, TRUE, DEFAULT_MAXDIVEUBQUOT, 0.0, 1.0, NULL, NULL) );
688  SCIP_CALL( SCIPaddRealParam(scip,
689  "heuristics/pscostdiving/maxdiveavgquot",
690  "maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound) where diving is performed (0.0: no limit)",
691  &heurdata->maxdiveavgquot, TRUE, DEFAULT_MAXDIVEAVGQUOT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
692  SCIP_CALL( SCIPaddRealParam(scip,
693  "heuristics/pscostdiving/maxdiveubquotnosol",
694  "maximal UBQUOT when no solution was found yet (0.0: no limit)",
695  &heurdata->maxdiveubquotnosol, TRUE, DEFAULT_MAXDIVEUBQUOTNOSOL, 0.0, 1.0, NULL, NULL) );
696  SCIP_CALL( SCIPaddRealParam(scip,
697  "heuristics/pscostdiving/maxdiveavgquotnosol",
698  "maximal AVGQUOT when no solution was found yet (0.0: no limit)",
699  &heurdata->maxdiveavgquotnosol, TRUE, DEFAULT_MAXDIVEAVGQUOTNOSOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
700  SCIP_CALL( SCIPaddBoolParam(scip,
701  "heuristics/pscostdiving/backtrack",
702  "use one level of backtracking if infeasibility is encountered?",
703  &heurdata->backtrack, FALSE, DEFAULT_BACKTRACK, NULL, NULL) );
704 
705  return SCIP_OKAY;
706 }
707 
708