1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 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 visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file heur_nlpdiving.c
17 * @ingroup DEFPLUGINS_HEUR
18 * @brief NLP diving heuristic that chooses fixings w.r.t. the fractionalities
19 * @author Timo Berthold
20 * @author Stefan Vigerske
21 */
22
23 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
24
25 #include "blockmemshell/memory.h"
26 #include "nlpi/nlpi.h"
27 #include "scip/heur_nlpdiving.h"
28 #include "scip/heur_subnlp.h"
29 #include "scip/heur_undercover.h"
30 #include "scip/pub_event.h"
31 #include "scip/pub_heur.h"
32 #include "scip/pub_message.h"
33 #include "scip/pub_misc.h"
34 #include "scip/pub_sol.h"
35 #include "scip/pub_var.h"
36 #include "scip/scip_branch.h"
37 #include "scip/scip_copy.h"
38 #include "scip/scip_event.h"
39 #include "scip/scip_general.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_nlp.h"
45 #include "scip/scip_nodesel.h"
46 #include "scip/scip_numerics.h"
47 #include "scip/scip_param.h"
48 #include "scip/scip_prob.h"
49 #include "scip/scip_probing.h"
50 #include "scip/scip_randnumgen.h"
51 #include "scip/scip_sol.h"
52 #include "scip/scip_solve.h"
53 #include "scip/scip_solvingstats.h"
54 #include "scip/scip_timing.h"
55 #include "scip/scip_tree.h"
56 #include "scip/scip_var.h"
57 #include <string.h>
58
59 #define HEUR_NAME "nlpdiving"
60 #define HEUR_DESC "NLP diving heuristic that chooses fixings w.r.t. the fractionalities"
61 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_DIVING
62 #define HEUR_PRIORITY -1003000
63 #define HEUR_FREQ 10
64 #define HEUR_FREQOFS 3
65 #define HEUR_MAXDEPTH -1
66 #define HEUR_TIMING SCIP_HEURTIMING_AFTERLPPLUNGE
67 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
68
69 /* event handler properties */
70 #define EVENTHDLR_NAME "Nlpdiving"
71 #define EVENTHDLR_DESC "bound change event handler for " HEUR_NAME " heuristic"
72
73
74 /*
75 * Default parameter settings
76 */
77
78 #define DEFAULT_MINRELDEPTH 0.0 /**< minimal relative depth to start diving */
79 #define DEFAULT_MAXRELDEPTH 1.0 /**< maximal relative depth to start diving */
80 #define DEFAULT_MAXNLPITERABS 200 /**< minimial absolute number of allowed NLP iterations */
81 #define DEFAULT_MAXNLPITERREL 10 /**< additional allowed number of NLP iterations relative to successfully found solutions */
82 #define DEFAULT_MAXDIVEUBQUOT 0.8 /**< maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound)
83 * where diving is performed (0.0: no limit) */
84 #define DEFAULT_MAXDIVEAVGQUOT 0.0 /**< maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound)
85 * where diving is performed (0.0: no limit) */
86 #define DEFAULT_MAXDIVEUBQUOTNOSOL 0.1 /**< maximal UBQUOT when no solution was found yet (0.0: no limit) */
87 #define DEFAULT_MAXDIVEAVGQUOTNOSOL 0.0 /**< maximal AVGQUOT when no solution was found yet (0.0: no limit) */
88 #define DEFAULT_MINSUCCQUOT 0.1 /**< heuristic will not run if less then this percentage of calls succeeded (0.0: no limit) */
89 #define DEFAULT_MAXFEASNLPS 10 /**< maximal number of NLPs with feasible solution to solve during one dive */
90 #define DEFAULT_FIXQUOT 0.2 /**< percentage of fractional variables that should be fixed before the next NLP solve */
91 #define DEFAULT_BACKTRACK TRUE /**< use one level of backtracking if infeasibility is encountered? */
92 #define DEFAULT_LP FALSE /**< should the LP relaxation be solved before the NLP relaxation? */
93 #define DEFAULT_PREFERLPFRACS FALSE /**< prefer variables that are also fractional in LP solution? */
94 #define DEFAULT_PREFERCOVER TRUE /**< should variables in a minimal cover be preferred? */
95 #define DEFAULT_SOLVESUBMIP FALSE /**< should a sub-MIP be solved if all cover variables are fixed? */
96 #define DEFAULT_NLPSTART 's' /**< which point should be used as starting point for the NLP solver? */
97 #define DEFAULT_VARSELRULE 'd' /**< which variable selection should be used? ('f'ractionality, 'c'oefficient,
98 * 'p'seudocost, 'g'uided, 'd'ouble)
99 */
100 #define DEFAULT_NLPFASTFAIL TRUE /**< should the NLP solver stop early if it converges slow? */
101 #define DEFAULT_RANDSEED 97 /**< initial random seed */
102
103 #define MINNLPITER 10 /**< minimal number of NLP iterations allowed in each NLP solving call */
104
105 /* locally defined heuristic data */
106 struct SCIP_HeurData
107 {
108 SCIP_SOL* sol; /**< working solution */
109 SCIP_Real minreldepth; /**< minimal relative depth to start diving */
110 SCIP_Real maxreldepth; /**< maximal relative depth to start diving */
111 int maxnlpiterabs; /**< minimial absolute number of allowed NLP iterations */
112 int maxnlpiterrel; /**< additional allowed number of NLP iterations relative to successfully found solutions */
113 SCIP_Real maxdiveubquot; /**< maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound)
114 * where diving is performed (0.0: no limit) */
115 SCIP_Real maxdiveavgquot; /**< maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound)
116 * where diving is performed (0.0: no limit) */
117 SCIP_Real maxdiveubquotnosol; /**< maximal UBQUOT when no solution was found yet (0.0: no limit) */
118 SCIP_Real maxdiveavgquotnosol;/**< maximal AVGQUOT when no solution was found yet (0.0: no limit) */
119 int maxfeasnlps; /**< maximal number of NLPs with feasible solution to solve during one dive */
120 SCIP_Real minsuccquot; /**< heuristic will not run if less then this percentage of calls succeeded (0.0: no limit) */
121 SCIP_Real fixquot; /**< percentage of fractional variables that should be fixed before the next NLP solve */
122 SCIP_Bool backtrack; /**< use one level of backtracking if infeasibility is encountered? */
123 SCIP_Bool lp; /**< should the LP relaxation be solved before the NLP relaxation? */
124 SCIP_Bool preferlpfracs; /**< prefer variables that are also fractional in LP solution? */
125 SCIP_Bool prefercover; /**< should variables in a minimal cover be preferred? */
126 SCIP_Bool solvesubmip; /**< should a sub-MIP be solved if all cover variables are fixed? */
127 SCIP_Bool nlpfastfail; /**< should the NLP solver stop early if it converges slow? */
128 char nlpstart; /**< which point should be used as starting point for the NLP solver? */
129 char varselrule; /**< which variable selection should be used? ('f'ractionality, 'c'oefficient,
130 * 'p'seudocost, 'g'uided, 'd'ouble)
131 */
132
133 int nnlpiterations; /**< NLP iterations used in this heuristic */
134 int nsuccess; /**< number of runs that produced at least one feasible solution */
135 int nfixedcovervars; /**< number of variables in the cover that are already fixed */
136 #ifdef SCIP_STATISTIC
137 int nnlpsolves; /**< number of NLP solves */
138 int nfailcutoff; /**< number of fails due to cutoff */
139 int nfaildepth; /**< number of fails due to too deep */
140 int nfailnlperror; /**< number of fails due to NLP error */
141 #endif
142 SCIP_EVENTHDLR* eventhdlr; /**< event handler for bound change events */
143 SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
144 };
145
146
147 /*
148 * local methods
149 */
150
151 /** gets fractional variables of last NLP solution along with solution values and fractionalities
152 *
153 * @return \ref SCIP_OKAY is returned if everything worked. Otherwise a suitable error code is passed. See \ref
154 * SCIP_Retcode "SCIP_RETCODE" for a complete list of error codes.
155 *
156 * @pre This method can be called if SCIP is in one of the following stages:
157 * - \ref SCIP_STAGE_INITSOLVE
158 * - \ref SCIP_STAGE_SOLVING
159 */
160 static
getNLPFracVars(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR *** nlpcands,SCIP_Real ** nlpcandssol,SCIP_Real ** nlpcandsfrac,int * nnlpcands)161 SCIP_RETCODE getNLPFracVars(
162 SCIP* scip, /**< SCIP data structure */
163 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
164 SCIP_VAR*** nlpcands, /**< pointer to store the array of NLP fractional variables, or NULL */
165 SCIP_Real** nlpcandssol, /**< pointer to store the array of NLP fractional variables solution values, or NULL */
166 SCIP_Real** nlpcandsfrac, /**< pointer to store the array of NLP fractional variables fractionalities, or NULL */
167 int* nnlpcands /**< pointer to store the number of NLP fractional variables , or NULL */
168 )
169 {
170 int c;
171
172 assert(scip != NULL);
173 assert(heurdata != NULL);
174 assert(nlpcands != NULL);
175 assert(nlpcandssol != NULL);
176 assert(nlpcandsfrac != NULL);
177 assert(nnlpcands != NULL);
178
179 /* get fractional variables that should be integral */
180 SCIP_CALL( SCIPgetNLPFracVars(scip, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, NULL) );
181
182 /* values may be outside the domain in exact arithmetic, but inside the domain within relative tolerance, and still
183 * slightly fractional, because SCIPisFeasIntegral() uses absolute tolerance; project value onto domain to avoid this
184 * (example: primsol=29.99999853455704, lower bound = 30)
185 */
186 for( c = 0; c < *nnlpcands; ++c )
187 {
188 assert(!SCIPisFeasIntegral(scip, (*nlpcandssol)[c]));
189
190 if( (*nlpcandssol)[c] < SCIPvarGetLbLocal((*nlpcands)[c]) || (*nlpcandssol)[c] > SCIPvarGetUbLocal((*nlpcands)[c]) )
191 {
192 SCIP_Real newval;
193
194 newval = ((*nlpcandssol)[c] < SCIPvarGetLbLocal((*nlpcands)[c]))
195 ? SCIPvarGetLbLocal((*nlpcands)[c]) - 0.5*SCIPfeastol(scip)
196 : SCIPvarGetUbLocal((*nlpcands)[c]) + 0.5*SCIPfeastol(scip);
197
198 assert(SCIPisFeasIntegral(scip, newval));
199
200 SCIP_CALL( SCIPsetSolVal(scip, heurdata->sol, (*nlpcands)[c], newval) );
201
202 (*nnlpcands)--;
203
204 if( c < *nnlpcands )
205 {
206 (*nlpcands)[c] = (*nlpcands)[*nnlpcands];
207 (*nlpcandssol)[c] = (*nlpcandssol)[*nnlpcands];
208 (*nlpcandsfrac)[c] = (*nlpcandsfrac)[*nnlpcands];
209 }
210 }
211 }
212
213 /* prefer decisions on variables which are also fractional in LP solution */
214 if( heurdata->preferlpfracs && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
215 {
216 for( c = 0; c < *nnlpcands; ++c )
217 {
218 if( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, NULL, (*nlpcands)[c])) )
219 (*nlpcandsfrac)[c] *= 100.0;
220 }
221 }
222
223 return SCIP_OKAY;
224 }
225
226 /** finds best candidate variable w.r.t. fractionality:
227 * - prefer variables that may not be rounded without destroying NLP feasibility:
228 * - of these variables, round least fractional variable in corresponding direction
229 * - if all remaining fractional variables may be rounded without destroying NLP feasibility:
230 * - round variable with least increasing objective value
231 * - binary variables are prefered
232 * - variables in a minimal cover or variables that are also fractional in an optimal LP solution might
233 * also be prefered if a correpsonding parameter is set
234 */
235 static
chooseFracVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** nlpcands,SCIP_Real * nlpcandssol,SCIP_Real * nlpcandsfrac,int nnlpcands,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)236 SCIP_RETCODE chooseFracVar(
237 SCIP* scip, /**< original SCIP data structure */
238 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
239 SCIP_VAR** nlpcands, /**< array of NLP fractional variables */
240 SCIP_Real* nlpcandssol, /**< array of NLP fractional variables solution values */
241 SCIP_Real* nlpcandsfrac, /**< array of NLP fractional variables fractionalities */
242 int nnlpcands, /**< number of NLP fractional variables */
243 SCIP_HASHMAP* varincover, /**< hash map for variables */
244 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
245 int* bestcand, /**< pointer to store the index of the best candidate variable */
246 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
247 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
248 )
249 {
250 SCIP_Real bestobjgain;
251 SCIP_Real bestfrac;
252 SCIP_Bool bestcandmayrounddown;
253 SCIP_Bool bestcandmayroundup;
254 int c;
255
256 /* check preconditions */
257 assert(scip != NULL);
258 assert(heurdata != NULL);
259 assert(nlpcands != NULL);
260 assert(nlpcandssol != NULL);
261 assert(nlpcandsfrac != NULL);
262 assert(covercomputed == (varincover != NULL));
263 assert(bestcand != NULL);
264 assert(bestcandmayround != NULL);
265 assert(bestcandroundup != NULL);
266
267 bestcandmayrounddown = TRUE;
268 bestcandmayroundup = TRUE;
269 bestobjgain = SCIPinfinity(scip);
270 bestfrac = SCIP_INVALID;
271
272 for( c = 0; c < nnlpcands; ++c )
273 {
274 SCIP_VAR* var;
275 SCIP_Bool mayrounddown;
276 SCIP_Bool mayroundup;
277 SCIP_Bool roundup;
278 SCIP_Real frac;
279 SCIP_Real obj;
280 SCIP_Real objgain;
281
282 var = nlpcands[c];
283
284 mayrounddown = SCIPvarMayRoundDown(var);
285 mayroundup = SCIPvarMayRoundUp(var);
286 frac = nlpcandsfrac[c];
287 obj = SCIPvarGetObj(var);
288
289 /* since we are not solving the NLP after each fixing, the old NLP solution might be outside the propagated bounds */
290 if( SCIPisLT(scip, nlpcandssol[c], SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpcandssol[c], SCIPvarGetUbLocal(var)) )
291 continue;
292
293 if( mayrounddown || mayroundup )
294 {
295 /* the candidate may be rounded: choose this candidate only, if the best candidate may also be rounded */
296 if( bestcandmayrounddown || bestcandmayroundup )
297 {
298 /* choose rounding direction:
299 * - if variable may be rounded in both directions, round corresponding to the fractionality
300 * - otherwise, round in the infeasible direction, because feasible direction is tried by rounding
301 * the current fractional solution
302 */
303 if( mayrounddown && mayroundup )
304 {
305 if( SCIPisEQ(scip, frac, 0.5) )
306 roundup = (SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0);
307 else
308 roundup = (frac > 0.5);
309 }
310 else
311 roundup = mayrounddown;
312
313 if( roundup )
314 {
315 frac = 1.0 - frac;
316 objgain = frac*obj;
317 }
318 else
319 objgain = -frac*obj;
320
321 /* penalize too small fractions */
322 if( SCIPisEQ(scip, frac, 0.01) )
323 {
324 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
325 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
326 */
327 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
328 objgain *= 1000.0;
329 }
330 else if( frac < 0.01 )
331 objgain *= 1000.0;
332
333 /* prefer decisions on binary variables */
334 if( !SCIPvarIsBinary(var) )
335 objgain *= 1000.0;
336
337 /* prefer decisions on cover variables */
338 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
339 objgain *= 1000.0;
340
341 /* check, if candidate is new best candidate */
342 if( SCIPisLT(scip, objgain, bestobjgain) || (SCIPisEQ(scip, objgain, bestobjgain) && frac < bestfrac) )
343 {
344 *bestcand = c;
345 bestobjgain = objgain;
346 bestfrac = frac;
347 bestcandmayrounddown = mayrounddown;
348 bestcandmayroundup = mayroundup;
349 *bestcandroundup = roundup;
350 }
351 }
352 }
353 else
354 {
355 /* the candidate may not be rounded */
356 if( SCIPisEQ(scip, frac, 0.5) )
357 roundup = (SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0);
358 else if( frac < 0.5 )
359 roundup = FALSE;
360 else
361 roundup = TRUE;
362
363 /* adjust fractional part */
364 if( roundup )
365 frac = 1.0 - frac;
366
367 /* penalize too small fractions */
368 if( SCIPisEQ(scip, frac, 0.01) )
369 {
370 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
371 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
372 */
373 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
374 frac += 10.0;
375 }
376 else if( frac < 0.01 )
377 frac += 10.0;
378
379 /* prefer decisions on binary variables */
380 if( !SCIPvarIsBinary(var) )
381 frac *= 1000.0;
382
383 /* prefer decisions on cover variables */
384 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
385 frac *= 1000.0;
386
387 /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
388 if( bestcandmayrounddown || bestcandmayroundup || frac < bestfrac )
389 {
390 *bestcand = c;
391 bestfrac = frac;
392 bestcandmayrounddown = FALSE;
393 bestcandmayroundup = FALSE;
394 *bestcandroundup = roundup;
395 }
396 assert(bestfrac < SCIP_INVALID);
397 }
398 }
399
400 *bestcandmayround = bestcandmayroundup || bestcandmayrounddown;
401
402 return SCIP_OKAY;
403 }
404
405 /** finds best candidate variable w.r.t. vector length:
406 * - round variable with a small ratio between the increase in the objective and the locking numbers
407 * - binary variables are prefered
408 * - variables in a minimal cover or variables that are also fractional in an optimal LP solution might
409 * also be prefered if a corresponding parameter is set
410 */
411 static
chooseVeclenVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** nlpcands,SCIP_Real * nlpcandssol,SCIP_Real * nlpcandsfrac,int nnlpcands,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)412 SCIP_RETCODE chooseVeclenVar(
413 SCIP* scip, /**< original SCIP data structure */
414 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
415 SCIP_VAR** nlpcands, /**< array of NLP fractional variables */
416 SCIP_Real* nlpcandssol, /**< array of NLP fractional variables solution values */
417 SCIP_Real* nlpcandsfrac, /**< array of NLP fractional variables fractionalities */
418 int nnlpcands, /**< number of NLP fractional variables */
419 SCIP_HASHMAP* varincover, /**< hash map for variables */
420 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
421 int* bestcand, /**< pointer to store the index of the best candidate variable */
422 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
423 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
424 )
425 {
426 SCIP_Real bestscore;
427 int c;
428
429 /* check preconditions */
430 assert(scip != NULL);
431 assert(heurdata != NULL);
432 assert(nlpcands != NULL);
433 assert(nlpcandsfrac != NULL);
434 assert(nlpcandssol != NULL);
435 assert(bestcand != NULL);
436 assert(bestcandmayround != NULL);
437 assert(bestcandroundup != NULL);
438
439 *bestcandmayround = TRUE;
440 bestscore = SCIP_REAL_MAX;
441
442 /* get best candidate */
443 for( c = 0; c < nnlpcands; ++c )
444 {
445 SCIP_VAR* var;
446
447 SCIP_Real obj;
448 SCIP_Real objdelta;
449 SCIP_Real frac;
450 SCIP_Real score;
451
452 int nlocks;
453
454 SCIP_Bool roundup;
455
456 var = nlpcands[c];
457
458 /* since we are not solving the NLP after each fixing, the old NLP solution might be outside the propagated bounds */
459 if( SCIPisLT(scip, nlpcandssol[c], SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpcandssol[c], SCIPvarGetUbLocal(var)) )
460 continue;
461
462 frac = nlpcandsfrac[c];
463 obj = SCIPvarGetObj(var);
464 roundup = (obj >= 0.0);
465 objdelta = (roundup ? (1.0-frac)*obj : -frac * obj);
466 assert(objdelta >= 0.0);
467
468 /* check whether the variable is roundable */
469 *bestcandmayround = *bestcandmayround && (SCIPvarMayRoundDown(var) || SCIPvarMayRoundUp(var));
470 nlocks = SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL) + SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL);
471
472 /* smaller score is better */
473 score = (objdelta + SCIPsumepsilon(scip))/((SCIP_Real)nlocks+1.0);
474
475 /* prefer decisions on binary variables */
476 if( SCIPvarGetType(var) != SCIP_VARTYPE_BINARY )
477 score *= 1000.0;
478
479 /* prefer decisions on cover variables */
480 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
481 score *= 1000;
482
483 /* check, if candidate is new best candidate */
484 if( score < bestscore )
485 {
486 *bestcand = c;
487 bestscore = score;
488 *bestcandroundup = roundup;
489 }
490 }
491
492 return SCIP_OKAY;
493 }
494
495
496 /** finds best candidate variable w.r.t. locking numbers:
497 * - prefer variables that may not be rounded without destroying LP feasibility:
498 * - of these variables, round variable with least number of locks in corresponding direction
499 * - if all remaining fractional variables may be rounded without destroying LP feasibility:
500 * - round variable with least number of locks in opposite of its feasible rounding direction
501 * - binary variables are prefered
502 * - variables in a minimal cover or variables that are also fractional in an optimal LP solution might
503 * also be prefered if a correpsonding parameter is set
504 */
505 static
chooseCoefVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** nlpcands,SCIP_Real * nlpcandssol,SCIP_Real * nlpcandsfrac,int nnlpcands,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)506 SCIP_RETCODE chooseCoefVar(
507 SCIP* scip, /**< original SCIP data structure */
508 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
509 SCIP_VAR** nlpcands, /**< array of NLP fractional variables */
510 SCIP_Real* nlpcandssol, /**< array of NLP fractional variables solution values */
511 SCIP_Real* nlpcandsfrac, /**< array of NLP fractional variables fractionalities */
512 int nnlpcands, /**< number of NLP fractional variables */
513 SCIP_HASHMAP* varincover, /**< hash map for variables */
514 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
515 int* bestcand, /**< pointer to store the index of the best candidate variable */
516 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
517 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
518 )
519 {
520 SCIP_Bool bestcandmayrounddown;
521 SCIP_Bool bestcandmayroundup;
522 int bestnviolrows; /* number of violated rows for best candidate */
523 SCIP_Real bestcandfrac; /* fractionality of best candidate */
524 int c;
525
526 /* check preconditions */
527 assert(scip != NULL);
528 assert(heurdata != NULL);
529 assert(nlpcands != NULL);
530 assert(nlpcandsfrac != NULL);
531 assert(nlpcandssol != NULL);
532 assert(bestcand != NULL);
533 assert(bestcandmayround != NULL);
534 assert(bestcandroundup != NULL);
535
536 bestcandmayrounddown = TRUE;
537 bestcandmayroundup = TRUE;
538 bestnviolrows = INT_MAX;
539 bestcandfrac = SCIP_INVALID;
540
541 /* get best candidate */
542 for( c = 0; c < nnlpcands; ++c )
543 {
544 SCIP_VAR* var;
545
546 int nlocksdown;
547 int nlocksup;
548 int nviolrows;
549
550 SCIP_Bool mayrounddown;
551 SCIP_Bool mayroundup;
552 SCIP_Bool roundup;
553 SCIP_Real frac;
554
555 var = nlpcands[c];
556 mayrounddown = SCIPvarMayRoundDown(var);
557 mayroundup = SCIPvarMayRoundUp(var);
558 frac = nlpcandsfrac[c];
559
560 /* since we are not solving the NLP after each fixing, the old NLP solution might be outside the propagated bounds */
561 if( SCIPisLT(scip, nlpcandssol[c], SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpcandssol[c], SCIPvarGetUbLocal(var)) )
562 continue;
563
564 if( mayrounddown || mayroundup )
565 {
566 /* the candidate may be rounded: choose this candidate only, if the best candidate may also be rounded */
567 if( bestcandmayrounddown || bestcandmayroundup )
568 {
569 /* choose rounding direction:
570 * - if variable may be rounded in both directions, round corresponding to the fractionality
571 * - otherwise, round in the infeasible direction, because feasible direction is tried by rounding
572 * the current fractional solution
573 */
574 if( mayrounddown && mayroundup )
575 {
576 if( SCIPisEQ(scip, frac, 0.5) )
577 roundup = (SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0);
578 else
579 roundup = (frac > 0.5);
580 }
581 else
582 roundup = mayrounddown;
583
584 if( roundup )
585 {
586 frac = 1.0 - frac;
587 nviolrows = SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL);
588 }
589 else
590 nviolrows = SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL);
591
592 /* penalize too small fractions */
593 if( SCIPisEQ(scip, frac, 0.01) )
594 {
595 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
596 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
597 */
598 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
599 nviolrows *= 100;
600 }
601 else if( frac < 0.01 )
602 nviolrows *= 100;
603
604 /* prefer decisions on binary variables */
605 if( !SCIPvarIsBinary(var) )
606 nviolrows *= 1000;
607
608 /* prefer decisions on cover variables */
609 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
610 nviolrows *= 1000;
611
612 /* check, if candidate is new best candidate */
613 assert( (0.0 < frac && frac < 1.0) || SCIPvarIsBinary(var) );
614 if( nviolrows + frac < bestnviolrows + bestcandfrac )
615 {
616 *bestcand = c;
617 bestnviolrows = nviolrows;
618 bestcandfrac = frac;
619 bestcandmayrounddown = mayrounddown;
620 bestcandmayroundup = mayroundup;
621 *bestcandroundup = roundup;
622 }
623 }
624 }
625 else
626 {
627 /* the candidate may not be rounded */
628 nlocksdown = SCIPvarGetNLocksDownType(var, SCIP_LOCKTYPE_MODEL);
629 nlocksup = SCIPvarGetNLocksUpType(var, SCIP_LOCKTYPE_MODEL);
630
631 roundup = (nlocksdown > nlocksup);
632 if( !roundup )
633 {
634 roundup = (nlocksdown == nlocksup);
635 if( SCIPisEQ(scip, frac, 0.5) )
636 roundup = (roundup && (SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0));
637 else
638 roundup = (roundup && frac > 0.5);
639 }
640
641 if( roundup )
642 {
643 nviolrows = nlocksup;
644 frac = 1.0 - frac;
645 }
646 else
647 nviolrows = nlocksdown;
648
649 /* penalize too small fractions */
650 if( SCIPisEQ(scip, frac, 0.01) )
651 {
652 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
653 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
654 */
655 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
656 nviolrows *= 100;
657 }
658 else if( frac < 0.01 )
659 nviolrows *= 100;
660
661 /* prefer decisions on binary variables */
662 if( !SCIPvarIsBinary(var) )
663 nviolrows *= 100;
664
665 /* prefer decisions on cover variables */
666 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
667 nviolrows *= 1000;
668
669 /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
670 assert((0.0 < frac && frac < 1.0) || SCIPvarIsBinary(var));
671 if( bestcandmayrounddown || bestcandmayroundup || nviolrows + frac < bestnviolrows + bestcandfrac )
672 {
673 *bestcand = c;
674 bestnviolrows = nviolrows;
675 bestcandfrac = frac;
676 bestcandmayrounddown = FALSE;
677 bestcandmayroundup = FALSE;
678 *bestcandroundup = roundup;
679 }
680 assert(bestcandfrac < SCIP_INVALID);
681 }
682 }
683
684 *bestcandmayround = bestcandmayroundup || bestcandmayrounddown;
685
686 return SCIP_OKAY;
687 }
688
689 /** calculates the pseudocost score for a given variable w.r.t. a given solution value and a given rounding direction */
690 static
calcPscostQuot(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR * var,SCIP_Real primsol,SCIP_Real frac,int rounddir,SCIP_Real * pscostquot,SCIP_Bool * roundup,SCIP_Bool prefvar)691 void calcPscostQuot(
692 SCIP* scip, /**< SCIP data structure */
693 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
694 SCIP_VAR* var, /**< problem variable */
695 SCIP_Real primsol, /**< primal solution of variable */
696 SCIP_Real frac, /**< fractionality of variable */
697 int rounddir, /**< -1: round down, +1: round up, 0: select due to pseudo cost values */
698 SCIP_Real* pscostquot, /**< pointer to store pseudo cost quotient */
699 SCIP_Bool* roundup, /**< pointer to store whether the variable should be rounded up */
700 SCIP_Bool prefvar /**< should this variable be preferred because it is in a minimal cover? */
701 )
702 {
703 SCIP_Real pscostdown;
704 SCIP_Real pscostup;
705
706 assert(heurdata != NULL);
707 assert(pscostquot != NULL);
708 assert(roundup != NULL);
709 assert(SCIPisEQ(scip, frac, primsol - SCIPfeasFloor(scip, primsol)));
710
711 /* bound fractions to not prefer variables that are nearly integral */
712 frac = MAX(frac, 0.1);
713 frac = MIN(frac, 0.9);
714
715 /* get pseudo cost quotient */
716 pscostdown = SCIPgetVarPseudocostVal(scip, var, 0.0-frac);
717 pscostup = SCIPgetVarPseudocostVal(scip, var, 1.0-frac);
718 assert(pscostdown >= 0.0 && pscostup >= 0.0);
719
720 /* choose rounding direction
721 *
722 * to avoid performance variability caused by numerics we use random numbers to decide whether we want to roundup or
723 * round down if the values to compare are equal within tolerances.
724 */
725 if( rounddir == -1 )
726 *roundup = FALSE;
727 else if( rounddir == +1 )
728 *roundup = TRUE;
729 else if( SCIPisLT(scip, primsol, SCIPvarGetRootSol(var) - 0.4)
730 || (SCIPisEQ(scip, primsol, SCIPvarGetRootSol(var) - 0.4) && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0) )
731 *roundup = FALSE;
732 else if( SCIPisGT(scip, primsol, SCIPvarGetRootSol(var) + 0.4)
733 || (SCIPisEQ(scip, primsol, SCIPvarGetRootSol(var) + 0.4) && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0) )
734 *roundup = TRUE;
735 else if( SCIPisLT(scip, frac, 0.3) || (SCIPisEQ(scip, frac, 0.3) && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0) )
736 *roundup = FALSE;
737 else if( SCIPisGT(scip, frac, 0.7) || (SCIPisEQ(scip, frac, 0.7) && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0) )
738 *roundup = TRUE;
739 else if( SCIPisLT(scip, pscostdown, pscostup)
740 || (SCIPisEQ(scip, pscostdown, pscostup) && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0))
741 *roundup = FALSE;
742 else
743 *roundup = TRUE;
744
745 /* calculate pseudo cost quotient */
746 if( *roundup )
747 *pscostquot = sqrt(frac) * (1.0+pscostdown) / (1.0+pscostup);
748 else
749 *pscostquot = sqrt(1.0-frac) * (1.0+pscostup) / (1.0+pscostdown);
750
751 /* prefer decisions on binary variables */
752 if( SCIPvarIsBinary(var) )
753 (*pscostquot) *= 1000.0;
754
755 /* prefer decisions on cover variables */
756 if( prefvar )
757 (*pscostquot) *= 1000.0;
758 }
759
760 /** finds best candidate variable w.r.t. pseudo costs:
761 * - prefer variables that may not be rounded without destroying LP feasibility:
762 * - of these variables, round variable with largest rel. difference of pseudo cost values in corresponding
763 * direction
764 * - if all remaining fractional variables may be rounded without destroying LP feasibility:
765 * - round variable in the objective value direction
766 * - binary variables are prefered
767 * - variables in a minimal cover or variables that are also fractional in an optimal LP solution might
768 * also be prefered if a correpsonding parameter is set
769 */
770 static
choosePscostVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** nlpcands,SCIP_Real * nlpcandssol,SCIP_Real * nlpcandsfrac,int nnlpcands,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)771 SCIP_RETCODE choosePscostVar(
772 SCIP* scip, /**< original SCIP data structure */
773 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
774 SCIP_VAR** nlpcands, /**< array of NLP fractional variables */
775 SCIP_Real* nlpcandssol, /**< array of NLP fractional variables solution values */
776 SCIP_Real* nlpcandsfrac, /**< array of NLP fractional variables fractionalities */
777 int nnlpcands, /**< number of NLP fractional variables */
778 SCIP_HASHMAP* varincover, /**< hash map for variables */
779 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
780 int* bestcand, /**< pointer to store the index of the best candidate variable */
781 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
782 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
783 )
784 {
785 SCIP_Bool bestcandmayrounddown;
786 SCIP_Bool bestcandmayroundup;
787 SCIP_Real bestpscostquot;
788 int c;
789
790 /* check preconditions */
791 assert(scip != NULL);
792 assert(heurdata != NULL);
793 assert(nlpcands != NULL);
794 assert(nlpcandsfrac != NULL);
795 assert(nlpcandssol != NULL);
796 assert(bestcand != NULL);
797 assert(bestcandmayround != NULL);
798 assert(bestcandroundup != NULL);
799
800 bestcandmayrounddown = TRUE;
801 bestcandmayroundup = TRUE;
802 bestpscostquot = -1.0;
803
804 for( c = 0; c < nnlpcands; ++c )
805 {
806 SCIP_VAR* var;
807 SCIP_Real primsol;
808
809 SCIP_Bool mayrounddown;
810 SCIP_Bool mayroundup;
811 SCIP_Bool roundup;
812 SCIP_Bool prefvar;
813 SCIP_Real frac;
814 SCIP_Real pscostquot;
815
816 var = nlpcands[c];
817 mayrounddown = SCIPvarMayRoundDown(var);
818 mayroundup = SCIPvarMayRoundUp(var);
819 primsol = nlpcandssol[c];
820 frac = nlpcandsfrac[c];
821 prefvar = covercomputed && heurdata->prefercover && SCIPhashmapExists(varincover, var);
822 pscostquot = SCIP_INVALID;
823
824 if( SCIPisLT(scip, nlpcandssol[c], SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpcandssol[c], SCIPvarGetUbLocal(var)) )
825 continue;
826
827 if( mayrounddown || mayroundup )
828 {
829 /* the candidate may be rounded: choose this candidate only, if the best candidate may also be rounded */
830 if( bestcandmayrounddown || bestcandmayroundup )
831 {
832 /* choose rounding direction:
833 * - if variable may be rounded in both directions, round corresponding to the pseudo cost values
834 * - otherwise, round in the infeasible direction, because feasible direction is tried by rounding
835 * the current fractional solution
836 */
837 roundup = FALSE;
838 if( mayrounddown && mayroundup )
839 calcPscostQuot(scip, heurdata, var, primsol, frac, 0, &pscostquot, &roundup, prefvar);
840 else if( mayrounddown )
841 calcPscostQuot(scip, heurdata, var, primsol, frac, +1, &pscostquot, &roundup, prefvar);
842 else
843 calcPscostQuot(scip, heurdata, var, primsol, frac, -1, &pscostquot, &roundup, prefvar);
844
845 assert(!SCIPisInfinity(scip,ABS(pscostquot)));
846
847 /* check, if candidate is new best candidate */
848 if( pscostquot > bestpscostquot )
849 {
850 *bestcand = c;
851 bestpscostquot = pscostquot;
852 bestcandmayrounddown = mayrounddown;
853 bestcandmayroundup = mayroundup;
854 *bestcandroundup = roundup;
855 }
856 }
857 }
858 else
859 {
860 /* the candidate may not be rounded: calculate pseudo cost quotient and preferred direction */
861 calcPscostQuot(scip, heurdata, var, primsol, frac, 0, &pscostquot, &roundup, prefvar);
862 assert(!SCIPisInfinity(scip,ABS(pscostquot)));
863
864 /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
865 if( bestcandmayrounddown || bestcandmayroundup || pscostquot > bestpscostquot )
866 {
867 *bestcand = c;
868 bestpscostquot = pscostquot;
869 bestcandmayrounddown = FALSE;
870 bestcandmayroundup = FALSE;
871 *bestcandroundup = roundup;
872 }
873 }
874 }
875
876 *bestcandmayround = bestcandmayroundup || bestcandmayrounddown;
877
878 return SCIP_OKAY;
879 }
880
881 /** finds best candidate variable w.r.t. the incumbent solution:
882 * - prefer variables that may not be rounded without destroying LP feasibility:
883 * - of these variables, round a variable to its value in direction of incumbent solution, and choose the
884 * variable that is closest to its rounded value
885 * - if all remaining fractional variables may be rounded without destroying LP feasibility:
886 * - round variable in direction that destroys LP feasibility (other direction is checked by SCIProundSol())
887 * - round variable with least increasing objective value
888 * - binary variables are prefered
889 * - variables in a minimal cover or variables that are also fractional in an optimal LP solution might
890 * also be prefered if a correpsonding parameter is set
891 */
892 static
chooseGuidedVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** nlpcands,SCIP_Real * nlpcandssol,SCIP_Real * nlpcandsfrac,int nnlpcands,SCIP_SOL * bestsol,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)893 SCIP_RETCODE chooseGuidedVar(
894 SCIP* scip, /**< original SCIP data structure */
895 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
896 SCIP_VAR** nlpcands, /**< array of NLP fractional variables */
897 SCIP_Real* nlpcandssol, /**< array of NLP fractional variables solution values */
898 SCIP_Real* nlpcandsfrac, /**< array of NLP fractional variables fractionalities */
899 int nnlpcands, /**< number of NLP fractional variables */
900 SCIP_SOL* bestsol, /**< incumbent solution */
901 SCIP_HASHMAP* varincover, /**< hash map for variables */
902 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
903 int* bestcand, /**< pointer to store the index of the best candidate variable */
904 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
905 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
906 )
907 {
908 SCIP_Real bestobjgain;
909 SCIP_Real bestfrac;
910 SCIP_Bool bestcandmayrounddown;
911 SCIP_Bool bestcandmayroundup;
912 int c;
913
914 /* check preconditions */
915 assert(scip != NULL);
916 assert(heurdata != NULL);
917 assert(nlpcands != NULL);
918 assert(nlpcandsfrac != NULL);
919 assert(nlpcandssol != NULL);
920 assert(bestcand != NULL);
921 assert(bestcandmayround != NULL);
922 assert(bestcandroundup != NULL);
923
924 bestcandmayrounddown = TRUE;
925 bestcandmayroundup = TRUE;
926 bestobjgain = SCIPinfinity(scip);
927 bestfrac = SCIP_INVALID;
928
929 for( c = 0; c < nnlpcands; ++c )
930 {
931 SCIP_VAR* var;
932 SCIP_Real bestsolval;
933 SCIP_Real solval;
934 SCIP_Real obj;
935 SCIP_Real frac;
936 SCIP_Real objgain;
937
938 SCIP_Bool mayrounddown;
939 SCIP_Bool mayroundup;
940 SCIP_Bool roundup;
941
942 var = nlpcands[c];
943 mayrounddown = SCIPvarMayRoundDown(var);
944 mayroundup = SCIPvarMayRoundUp(var);
945 solval = nlpcandssol[c];
946 frac = nlpcandsfrac[c];
947 obj = SCIPvarGetObj(var);
948 bestsolval = SCIPgetSolVal(scip, bestsol, var);
949
950 /* since we are not solving the NLP after each fixing, the old NLP solution might be outside the propagated bounds */
951 if( SCIPisLT(scip, solval, SCIPvarGetLbLocal(var)) || SCIPisGT(scip, solval, SCIPvarGetUbLocal(var)) )
952 continue;
953
954 /* select default rounding direction
955 * try to avoid variability; decide randomly if the LP solution can contain some noise
956 */
957 if( SCIPisEQ(scip, solval, bestsolval) )
958 roundup = (SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 0);
959 else
960 roundup = (solval < bestsolval);
961
962 if( mayrounddown || mayroundup )
963 {
964 /* the candidate may be rounded: choose this candidate only, if the best candidate may also be rounded */
965 if( bestcandmayrounddown || bestcandmayroundup )
966 {
967 /* choose rounding direction:
968 * - if variable may be rounded in both directions, round corresponding to its value in incumbent solution
969 * - otherwise, round in the infeasible direction, because feasible direction is tried by rounding
970 * the current fractional solution with SCIProundSol()
971 */
972 if( !mayrounddown || !mayroundup )
973 roundup = mayrounddown;
974
975 if( roundup )
976 {
977 frac = 1.0 - frac;
978 objgain = frac*obj;
979 }
980 else
981 objgain = -frac*obj;
982
983 /* penalize too small fractions */
984 if( SCIPisEQ(scip, frac, 0.01) )
985 {
986 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
987 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
988 */
989 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
990 objgain *= 1000.0;
991 }
992 else if( frac < 0.01 )
993 objgain *= 1000.0;
994
995 /* prefer decisions on binary variables */
996 if( !SCIPvarIsBinary(var) )
997 objgain *= 1000.0;
998
999 /* prefer decisions on cover variables */
1000 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
1001 objgain *= 1000.0;
1002
1003 /* check, if candidate is new best candidate */
1004 if( SCIPisLT(scip, objgain, bestobjgain) || (SCIPisEQ(scip, objgain, bestobjgain) && frac < bestfrac) )
1005 {
1006 *bestcand = c;
1007 bestobjgain = objgain;
1008 bestfrac = frac;
1009 bestcandmayrounddown = mayrounddown;
1010 bestcandmayroundup = mayroundup;
1011 *bestcandroundup = roundup;
1012 }
1013 }
1014 }
1015 else
1016 {
1017 /* the candidate may not be rounded */
1018 if( roundup )
1019 frac = 1.0 - frac;
1020
1021 /* penalize too small fractions */
1022 if( SCIPisEQ(scip, frac, 0.01) )
1023 {
1024 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
1025 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
1026 */
1027 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
1028 frac += 10.0;
1029 }
1030 else if( frac < 0.01 )
1031 frac += 10.0;
1032
1033 /* prefer decisions on binary variables */
1034 if( !SCIPvarIsBinary(var) )
1035 frac *= 1000.0;
1036
1037 /* prefer decisions on cover variables */
1038 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
1039 frac *= 1000.0;
1040
1041 /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
1042 if( bestcandmayrounddown || bestcandmayroundup || frac < bestfrac )
1043 {
1044 *bestcand = c;
1045 bestfrac = frac;
1046 bestcandmayrounddown = FALSE;
1047 bestcandmayroundup = FALSE;
1048 *bestcandroundup = roundup;
1049 }
1050 }
1051 }
1052
1053 *bestcandmayround = bestcandmayroundup || bestcandmayrounddown;
1054
1055 return SCIP_OKAY;
1056 }
1057
1058 /** finds best candidate variable w.r.t. both, the LP and the NLP solution:
1059 * - choose a variable for which the sum of the distances from the relaxations' solutions to a common
1060 * integer value is minimal
1061 * - binary variables are prefered
1062 * - variables in a minimal cover might be prefered if a corresponding parameter is set
1063 */
1064 static
chooseDoubleVar(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** pseudocands,SCIP_Real * pseudocandsnlpsol,SCIP_Real * pseudocandslpsol,int npseudocands,SCIP_HASHMAP * varincover,SCIP_Bool covercomputed,int * bestcand,SCIP_Real * bestboundval,SCIP_Bool * bestcandmayround,SCIP_Bool * bestcandroundup)1065 SCIP_RETCODE chooseDoubleVar(
1066 SCIP* scip, /**< original SCIP data structure */
1067 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
1068 SCIP_VAR** pseudocands, /**< array of non-fixed variables */
1069 SCIP_Real* pseudocandsnlpsol, /**< array of NLP solution values */
1070 SCIP_Real* pseudocandslpsol, /**< array of LP solution values */
1071 int npseudocands, /**< number of NLP fractional variables */
1072 SCIP_HASHMAP* varincover, /**< hash map for variables */
1073 SCIP_Bool covercomputed, /**< has a minimal cover been computed? */
1074 int* bestcand, /**< pointer to store the index of the best candidate variable */
1075 SCIP_Real* bestboundval, /**< pointer to store the bound, the best candidate should be rounded to */
1076 SCIP_Bool* bestcandmayround, /**< pointer to store whether best candidate is trivially roundable */
1077 SCIP_Bool* bestcandroundup /**< pointer to store whether best candidate should be rounded up */
1078 )
1079 {
1080 SCIP_Real bestfrac;
1081 int c;
1082
1083 /* check preconditions */
1084 assert(scip != NULL);
1085 assert(heurdata != NULL);
1086 assert(pseudocands != NULL);
1087 assert(pseudocandsnlpsol != NULL);
1088 assert(pseudocandslpsol != NULL);
1089 assert(covercomputed == (varincover != NULL));
1090 assert(bestcand != NULL);
1091 assert(bestcandmayround != NULL);
1092 assert(bestcandroundup != NULL);
1093
1094 bestfrac = SCIP_INVALID;
1095
1096 for( c = 0; c < npseudocands; ++c )
1097 {
1098 SCIP_VAR* var;
1099 SCIP_Bool mayround;
1100 SCIP_Bool roundup;
1101
1102 SCIP_Real frac;
1103 SCIP_Real lpsol;
1104 SCIP_Real nlpsol;
1105 SCIP_Real lpsolfloor;
1106 SCIP_Real nlpsolfloor;
1107 SCIP_Real lpsolceil;
1108 SCIP_Real nlpsolceil;
1109 SCIP_Real boundval;
1110 SCIP_Real floorval;
1111 SCIP_Real ceilval;
1112
1113 var = pseudocands[c];
1114 lpsol = pseudocandslpsol[c];
1115 nlpsol = pseudocandsnlpsol[c];
1116
1117 assert(SCIPvarGetUbLocal(var)-SCIPvarGetLbLocal(var) > 0.5);
1118 assert(SCIPisLE(scip, SCIPvarGetLbLocal(var), lpsol) && SCIPisLE(scip, lpsol, SCIPvarGetUbLocal(var)));
1119
1120 /* since we are not solving the NLP after each fixing, the old NLP solution might be outside the propagated bounds */
1121 if( SCIPisLT(scip, nlpsol, SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpsol, SCIPvarGetUbLocal(var)) )
1122 continue;
1123
1124 mayround = SCIPvarMayRoundDown(var) || SCIPvarMayRoundUp(var);
1125
1126 /* if this candidate is trivially roundable, and we already know a candidate that is not, continue */
1127 if( mayround && !(*bestcandmayround) )
1128 continue;
1129
1130 if( SCIPisFeasEQ(scip, lpsol, nlpsol) && SCIPisFeasIntegral(scip, lpsol))
1131 continue;
1132
1133 lpsolfloor = SCIPfeasFloor(scip, lpsol);
1134 nlpsolfloor = SCIPfeasFloor(scip, nlpsol);
1135 lpsolceil = SCIPfeasCeil(scip, lpsol);
1136 nlpsolceil = SCIPfeasCeil(scip, nlpsol);
1137 floorval = MIN(lpsolfloor,nlpsolfloor);
1138 ceilval = MAX(lpsolceil,nlpsolceil);
1139
1140 /* if both values are in the same interval, find out which integer is (in sum) the closer one, this will be the
1141 * new bound. The minima and maxima are necessary since one or both values with be integer
1142 */
1143 if( SCIPvarIsBinary(var) || ceilval-floorval < 1.5 )
1144 {
1145 frac = 0.33*(lpsol-floorval) + 0.67*(nlpsol-floorval);
1146 if( frac < 0.5 )
1147 {
1148 roundup = FALSE;
1149 boundval = MIN(lpsolfloor,nlpsolfloor);
1150 }
1151 else
1152 {
1153 roundup = TRUE;
1154 frac = 1.0-frac;
1155 boundval = MAX(nlpsolceil,lpsolceil);
1156 }
1157 }
1158 else
1159 {
1160 /* determine new bound in the middle of both relaxations, such that the NLP stays feasible */
1161 SCIP_Real midval;
1162 midval = (nlpsol+lpsol)/2.0;
1163 roundup = nlpsol > lpsol;
1164 frac = ABS(nlpsol-lpsol);
1165
1166 if( roundup )
1167 boundval = SCIPfeasCeil(scip, midval);
1168 else
1169 boundval = SCIPfeasFloor(scip, midval);
1170
1171 assert(roundup == SCIPisGT(scip, nlpsol, boundval));
1172 }
1173
1174 /* penalize too small fractions */
1175 if( SCIPisEQ(scip, frac, 0.01) )
1176 {
1177 /* try to avoid variability; decide randomly if the LP solution can contain some noise.
1178 * use a 1:SCIP_PROBINGSCORE_PENALTYRATIO chance for increasing the fractionality, i.e., the score.
1179 */
1180 if( SCIPrandomGetInt(heurdata->randnumgen, 0, SCIP_PROBINGSCORE_PENALTYRATIO) == 0 )
1181 frac += 10.0;
1182 }
1183 else if( frac < 0.01 )
1184 frac += 10.0;
1185
1186 /* prefer decisions on binary variables */
1187 if( !SCIPvarIsBinary(var) )
1188 frac *= 1000.0;
1189
1190 /* prefer decisions on cover variables */
1191 if( covercomputed && heurdata->prefercover && !SCIPhashmapExists(varincover, var) )
1192 frac *= 1000.0;
1193
1194 /* check, if candidate is new best candidate: prefer unroundable candidates in any case */
1195 if( frac < bestfrac || (*bestcandmayround && !mayround) )
1196 {
1197 *bestcand = c;
1198 bestfrac = frac;
1199 *bestcandmayround = FALSE;
1200 *bestcandroundup = roundup;
1201 *bestboundval = boundval;
1202 }
1203 assert(bestfrac < SCIP_INVALID);
1204 }
1205
1206 if( *bestcandroundup )
1207 *bestboundval -= 0.5;
1208 else
1209 *bestboundval += 0.5;
1210
1211 return SCIP_OKAY;
1212 }
1213
1214 /** creates a new solution for the original problem by copying the solution of the subproblem */
1215 static
createNewSol(SCIP * scip,SCIP * subscip,SCIP_HEUR * heur,SCIP_HASHMAP * varmap,SCIP_SOL * subsol,SCIP_Bool * success)1216 SCIP_RETCODE createNewSol(
1217 SCIP* scip, /**< original SCIP data structure */
1218 SCIP* subscip, /**< SCIP structure of the subproblem */
1219 SCIP_HEUR* heur, /**< heuristic structure */
1220 SCIP_HASHMAP* varmap, /**< hash map for variables */
1221 SCIP_SOL* subsol, /**< solution of the subproblem */
1222 SCIP_Bool* success /**< used to store whether new solution was found or not */
1223 )
1224 {
1225 SCIP_VAR** vars; /* the original problem's variables */
1226 SCIP_Real* subsolvals; /* solution values of the subproblem */
1227 SCIP_SOL* newsol; /* solution to be created for the original problem */
1228 int nvars; /* the original problem's number of variables */
1229 SCIP_VAR* subvar;
1230 int i;
1231
1232 assert(scip != NULL);
1233 assert(subscip != NULL);
1234 assert(subsol != NULL);
1235
1236 /* get variables' data */
1237 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
1238
1239 SCIP_CALL( SCIPallocBufferArray(scip, &subsolvals, nvars) );
1240
1241 /* copy the solution */
1242 for( i = 0; i < nvars; ++i )
1243 {
1244 subvar = (SCIP_VAR*) SCIPhashmapGetImage(varmap, vars[i]);
1245 if( subvar == NULL )
1246 subsolvals[i] = MIN(MAX(0.0, SCIPvarGetLbLocal(vars[i])), SCIPvarGetUbLocal(vars[i])); /*lint !e666*/
1247 else
1248 subsolvals[i] = SCIPgetSolVal(subscip, subsol, subvar);
1249 }
1250
1251 /* create new solution for the original problem */
1252 SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) );
1253 SCIP_CALL( SCIPsetSolVals(scip, newsol, nvars, vars, subsolvals) );
1254
1255 /* try to add new solution to scip and free it immediately */
1256 SCIP_CALL( SCIPtrySolFree(scip, &newsol, FALSE, FALSE, TRUE, TRUE, TRUE, success) );
1257
1258 SCIPfreeBufferArray(scip, &subsolvals);
1259
1260 return SCIP_OKAY;
1261 }
1262
1263 /** todo setup and solve the subMIP */
1264 static
doSolveSubMIP(SCIP * scip,SCIP * subscip,SCIP_HEUR * heur,SCIP_VAR ** covervars,int ncovervars,SCIP_Bool * success)1265 SCIP_RETCODE doSolveSubMIP(
1266 SCIP* scip, /**< SCIP data structure */
1267 SCIP* subscip, /**< NLP diving subscip */
1268 SCIP_HEUR* heur, /**< heuristic data structure */
1269 SCIP_VAR** covervars, /**< variables in the cover, should be fixed locally */
1270 int ncovervars, /**< number of variables in the cover */
1271 SCIP_Bool* success /**< pointer to store whether a solution was found */
1272 )
1273 {
1274 SCIP_HASHMAP* varmap;
1275 SCIP_SOL** subsols;
1276 int c;
1277 int nsubsols;
1278
1279 assert(subscip != NULL);
1280 assert(scip != NULL);
1281 assert(heur != NULL);
1282
1283 /* create the variable mapping hash map */
1284 SCIP_CALL( SCIPhashmapCreate(&varmap, SCIPblkmem(subscip), SCIPgetNVars(scip)) );
1285
1286 *success = FALSE;
1287
1288 /* copy original problem to subproblem; do not copy pricers */
1289 SCIP_CALL( SCIPcopyConsCompression(scip, subscip, varmap, NULL, "undercoversub", NULL, NULL, 0, FALSE, FALSE, FALSE,
1290 TRUE, NULL) );
1291
1292 /* assert that cover variables are fixed in source and target SCIP */
1293 for( c = 0; c < ncovervars; c++)
1294 {
1295 assert(SCIPhashmapGetImage(varmap, covervars[c]) != NULL); /* cover variable cannot be relaxation-only, thus must have been copied */
1296 assert(SCIPisFeasEQ(scip, SCIPvarGetLbLocal(covervars[c]), SCIPvarGetUbLocal(covervars[c])));
1297 assert(SCIPisFeasEQ(scip, SCIPvarGetLbGlobal((SCIP_VAR*) SCIPhashmapGetImage(varmap, covervars[c])),
1298 SCIPvarGetUbGlobal((SCIP_VAR*) SCIPhashmapGetImage(varmap, covervars[c]))));
1299 }
1300
1301 /* set parameters for sub-SCIP */
1302
1303 /* do not abort subproblem on CTRL-C */
1304 SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", FALSE) );
1305
1306 #ifdef SCIP_DEBUG
1307 /* for debugging, enable full output */
1308 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
1309 SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 100000000) );
1310 #else
1311 /* disable statistic timing inside sub SCIP and output to console */
1312 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
1313 SCIP_CALL( SCIPsetBoolParam(subscip, "timing/statistictiming", FALSE) );
1314 #endif
1315
1316 /* set limits for the subproblem */
1317 SCIP_CALL( SCIPcopyLimits(scip, subscip) );
1318 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", (SCIP_Longint)100) );
1319 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", (SCIP_Longint)500) );
1320
1321 /* forbid recursive call of heuristics and separators solving sub-SCIPs */
1322 SCIP_CALL( SCIPsetSubscipsOff(subscip, TRUE) );
1323
1324 /* disable cutting plane separation */
1325 SCIP_CALL( SCIPsetSeparating(subscip, SCIP_PARAMSETTING_OFF, TRUE) );
1326
1327 /* disable expensive presolving */
1328 SCIP_CALL( SCIPsetPresolving(subscip, SCIP_PARAMSETTING_FAST, TRUE) );
1329
1330 /* use best estimate node selection */
1331 if( SCIPfindNodesel(subscip, "estimate") != NULL && !SCIPisParamFixed(subscip, "nodeselection/estimate/stdpriority") )
1332 {
1333 SCIP_CALL( SCIPsetIntParam(subscip, "nodeselection/estimate/stdpriority", INT_MAX/4) );
1334 }
1335
1336 /* use inference branching */
1337 if( SCIPfindBranchrule(subscip, "inference") != NULL && !SCIPisParamFixed(subscip, "branching/inference/priority") )
1338 {
1339 SCIP_CALL( SCIPsetIntParam(subscip, "branching/inference/priority", INT_MAX/4) );
1340 }
1341
1342 /* enable conflict analysis, disable analysis of boundexceeding LPs, and restrict conflict pool */
1343 if( !SCIPisParamFixed(subscip, "conflict/enable") )
1344 {
1345 SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/enable", TRUE) );
1346 }
1347 if( !SCIPisParamFixed(subscip, "conflict/useboundlp") )
1348 {
1349 SCIP_CALL( SCIPsetCharParam(subscip, "conflict/useboundlp", 'o') );
1350 }
1351 if( !SCIPisParamFixed(subscip, "conflict/maxstoresize") )
1352 {
1353 SCIP_CALL( SCIPsetIntParam(subscip, "conflict/maxstoresize", 100) );
1354 }
1355
1356 if( SCIPgetNSols(scip) > 0 )
1357 {
1358 SCIP_Real upperbound;
1359 SCIP_Real cutoffbound;
1360 SCIP_Real minimprove;
1361
1362 assert( !SCIPisInfinity(scip,SCIPgetUpperbound(scip)) );
1363
1364 upperbound = SCIPgetUpperbound(scip) - SCIPsumepsilon(scip);
1365 minimprove = 0.01;
1366
1367 if( !SCIPisInfinity(scip,-1.0*SCIPgetLowerbound(scip)) )
1368 {
1369 cutoffbound = (1-minimprove)*SCIPgetUpperbound(scip) + minimprove*SCIPgetLowerbound(scip);
1370 }
1371 else
1372 {
1373 if( SCIPgetUpperbound(scip) >= 0 )
1374 cutoffbound = (1 - minimprove)*SCIPgetUpperbound(scip);
1375 else
1376 cutoffbound = (1 + minimprove)*SCIPgetUpperbound(scip);
1377 }
1378 cutoffbound = MIN(upperbound, cutoffbound);
1379 SCIP_CALL( SCIPsetObjlimit(subscip, cutoffbound) );
1380 }
1381
1382 SCIP_CALL_ABORT( SCIPsolve(subscip) );
1383
1384 /* check, whether a solution was found;
1385 * due to numerics, it might happen that not all solutions are feasible -> try all solutions until one was accepted
1386 */
1387 nsubsols = SCIPgetNSols(subscip);
1388 subsols = SCIPgetSols(subscip);
1389 for( c = 0; c < nsubsols && !(*success); ++c )
1390 {
1391 SCIP_CALL( createNewSol(scip, subscip, heur, varmap, subsols[c], success) );
1392 }
1393
1394 SCIPhashmapFree(&varmap);
1395
1396 return SCIP_OKAY;
1397 }
1398
1399
1400 /** solves subproblem and passes best feasible solution to original SCIP instance */
1401 static
solveSubMIP(SCIP * scip,SCIP_HEUR * heur,SCIP_VAR ** covervars,int ncovervars,SCIP_Bool * success)1402 SCIP_RETCODE solveSubMIP(
1403 SCIP* scip, /**< SCIP data structure of the original problem */
1404 SCIP_HEUR* heur, /**< heuristic data structure */
1405 SCIP_VAR** covervars, /**< variables in the cover, should be fixed locally */
1406 int ncovervars, /**< number of variables in the cover */
1407 SCIP_Bool* success /**< pointer to store whether a solution was found */
1408 )
1409 {
1410 SCIP* subscip;
1411 SCIP_RETCODE retcode;
1412
1413 /* check whether there is enough time and memory left */
1414 SCIP_CALL( SCIPcheckCopyLimits(scip, success) );
1415
1416 if( !(*success) )
1417 return SCIP_OKAY;
1418
1419 /* create subproblem */
1420 SCIP_CALL( SCIPcreate(&subscip) );
1421
1422 retcode = doSolveSubMIP(scip, subscip, heur, covervars, ncovervars, success);
1423
1424 /* free sub-SCIP even if an error occurred during the subscip solve */
1425 SCIP_CALL( SCIPfree(&subscip) );
1426
1427 SCIP_CALL( retcode );
1428
1429 return SCIP_OKAY;
1430 }
1431
1432 /* ---------------- Callback methods of event handler ---------------- */
1433
1434 /* exec the event handler
1435 *
1436 * We update the number of variables fixed in the cover
1437 */
1438 static
SCIP_DECL_EVENTEXEC(eventExecNlpdiving)1439 SCIP_DECL_EVENTEXEC(eventExecNlpdiving)
1440 {
1441 SCIP_EVENTTYPE eventtype;
1442 SCIP_HEURDATA* heurdata;
1443 SCIP_VAR* var;
1444
1445 SCIP_Real oldbound;
1446 SCIP_Real newbound;
1447 SCIP_Real otherbound;
1448
1449 assert(eventhdlr != NULL);
1450 assert(eventdata != NULL);
1451 assert(strcmp(SCIPeventhdlrGetName(eventhdlr), EVENTHDLR_NAME) == 0);
1452 assert(event != NULL);
1453
1454 heurdata = (SCIP_HEURDATA*)eventdata;
1455 assert(heurdata != NULL);
1456 assert(0 <= heurdata->nfixedcovervars && heurdata->nfixedcovervars <= SCIPgetNVars(scip));
1457
1458 oldbound = SCIPeventGetOldbound(event);
1459 newbound = SCIPeventGetNewbound(event);
1460 var = SCIPeventGetVar(event);
1461
1462 eventtype = SCIPeventGetType(event);
1463 otherbound = (eventtype & SCIP_EVENTTYPE_LBCHANGED) ? SCIPvarGetUbLocal(var) : SCIPvarGetLbLocal(var);
1464
1465 switch( eventtype )
1466 {
1467 case SCIP_EVENTTYPE_LBTIGHTENED:
1468 case SCIP_EVENTTYPE_UBTIGHTENED:
1469 /* if cover variable is now fixed */
1470 if( SCIPisFeasEQ(scip, newbound, otherbound) && !SCIPisFeasEQ(scip, oldbound, otherbound) )
1471 {
1472 assert(!SCIPisEQ(scip, oldbound, otherbound));
1473 ++(heurdata->nfixedcovervars);
1474 }
1475 break;
1476 case SCIP_EVENTTYPE_LBRELAXED:
1477 case SCIP_EVENTTYPE_UBRELAXED:
1478 /* if cover variable is now unfixed */
1479 if( SCIPisFeasEQ(scip, oldbound, otherbound) && !SCIPisFeasEQ(scip, newbound, otherbound) )
1480 {
1481 assert(!SCIPisEQ(scip, newbound, otherbound));
1482 --(heurdata->nfixedcovervars);
1483 }
1484 break;
1485 default:
1486 SCIPerrorMessage("invalid event type.\n");
1487 return SCIP_INVALIDDATA;
1488 }
1489 assert(0 <= heurdata->nfixedcovervars && heurdata->nfixedcovervars <= SCIPgetNVars(scip));
1490
1491 /* SCIPdebugMsg(scip, "changed bound of cover variable <%s> from %f to %f (nfixedcovervars: %d).\n", SCIPvarGetName(var),
1492 oldbound, newbound, heurdata->nfixedcovervars); */
1493
1494 return SCIP_OKAY;
1495 }
1496
1497
1498 /*
1499 * Callback methods
1500 */
1501
1502 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
1503 static
SCIP_DECL_HEURCOPY(heurCopyNlpdiving)1504 SCIP_DECL_HEURCOPY(heurCopyNlpdiving)
1505 { /*lint --e{715}*/
1506 assert(scip != NULL);
1507 assert(heur != NULL);
1508 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1509
1510 /* call inclusion method of primal heuristic */
1511 SCIP_CALL( SCIPincludeHeurNlpdiving(scip) );
1512
1513 return SCIP_OKAY;
1514 }
1515
1516 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
1517 static
SCIP_DECL_HEURFREE(heurFreeNlpdiving)1518 SCIP_DECL_HEURFREE(heurFreeNlpdiving) /*lint --e{715}*/
1519 { /*lint --e{715}*/
1520 SCIP_HEURDATA* heurdata;
1521
1522 assert(heur != NULL);
1523 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1524 assert(scip != NULL);
1525
1526 /* free heuristic data */
1527 heurdata = SCIPheurGetData(heur);
1528 assert(heurdata != NULL);
1529 SCIPfreeBlockMemory(scip, &heurdata);
1530 SCIPheurSetData(heur, NULL);
1531
1532 return SCIP_OKAY;
1533 }
1534
1535
1536 /** initialization method of primal heuristic (called after problem was transformed) */
1537 static
SCIP_DECL_HEURINIT(heurInitNlpdiving)1538 SCIP_DECL_HEURINIT(heurInitNlpdiving) /*lint --e{715}*/
1539 { /*lint --e{715}*/
1540 SCIP_HEURDATA* heurdata;
1541
1542 assert(heur != NULL);
1543 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1544
1545 /* get heuristic data */
1546 heurdata = SCIPheurGetData(heur);
1547 assert(heurdata != NULL);
1548
1549 /* create working solution */
1550 SCIP_CALL( SCIPcreateSol(scip, &heurdata->sol, heur) );
1551
1552 /* create random number generator */
1553 SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, DEFAULT_RANDSEED, TRUE) );
1554
1555 /* initialize data */
1556 heurdata->nnlpiterations = 0;
1557 heurdata->nsuccess = 0;
1558 heurdata->nfixedcovervars = 0;
1559 SCIPstatistic(
1560 heurdata->nnlpsolves = 0;
1561 heurdata->nfailcutoff = 0;
1562 heurdata->nfaildepth = 0;
1563 heurdata->nfailnlperror = 0;
1564 );
1565
1566 return SCIP_OKAY;
1567 }
1568
1569
1570 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
1571 static
SCIP_DECL_HEUREXIT(heurExitNlpdiving)1572 SCIP_DECL_HEUREXIT(heurExitNlpdiving) /*lint --e{715}*/
1573 { /*lint --e{715}*/
1574 SCIP_HEURDATA* heurdata;
1575
1576 assert(heur != NULL);
1577 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1578
1579 /* get heuristic data */
1580 heurdata = SCIPheurGetData(heur);
1581 assert(heurdata != NULL);
1582
1583 /* free random number generator */
1584 SCIPfreeRandom(scip, &heurdata->randnumgen);
1585
1586 /* free working solution */
1587 SCIP_CALL( SCIPfreeSol(scip, &heurdata->sol) );
1588
1589 SCIPstatistic(
1590 if( strstr(SCIPgetProbName(scip), "_covering") == NULL && SCIPheurGetNCalls(heur) > 0 )
1591 {
1592 SCIPstatisticMessage("%-30s %5" SCIP_LONGINT_FORMAT " sols in %5" SCIP_LONGINT_FORMAT " runs, %6.1fs, %7d NLP iters in %5d NLP solves, %5.1f avg., %3d%% success %3d%% cutoff %3d%% depth %3d%% nlperror\n",
1593 SCIPgetProbName(scip), SCIPheurGetNSolsFound(heur), SCIPheurGetNCalls(heur), SCIPheurGetTime(heur),
1594 heurdata->nnlpiterations, heurdata->nnlpsolves, heurdata->nnlpiterations/MAX(1.0,(SCIP_Real)heurdata->nnlpsolves),
1595 (100*heurdata->nsuccess) / (int)SCIPheurGetNCalls(heur), (100*heurdata->nfailcutoff) / (int)SCIPheurGetNCalls(heur), (100*heurdata->nfaildepth) / (int)SCIPheurGetNCalls(heur), (100*heurdata->nfailnlperror) / (int)SCIPheurGetNCalls(heur)
1596 );
1597 }
1598 );
1599
1600 return SCIP_OKAY;
1601 }
1602
1603
1604 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
1605 static
SCIP_DECL_HEURINITSOL(heurInitsolNlpdiving)1606 SCIP_DECL_HEURINITSOL(heurInitsolNlpdiving)
1607 { /*lint --e{715}*/
1608 SCIP_HEUR* nlpheur;
1609
1610 if( !SCIPisNLPConstructed(scip) )
1611 return SCIP_OKAY;
1612
1613 /* find NLP local search heuristic */
1614 nlpheur = SCIPfindHeur(scip, "subnlp");
1615
1616 /* add global linear constraints to NLP relaxation */
1617 if( nlpheur != NULL )
1618 {
1619 SCIP_CALL( SCIPaddLinearConsToNlpHeurSubNlp(scip, nlpheur, TRUE, TRUE) );
1620 }
1621
1622 return SCIP_OKAY;
1623 }
1624
1625
1626 /** execution method of primal heuristic */
1627 static
SCIP_DECL_HEUREXEC(heurExecNlpdiving)1628 SCIP_DECL_HEUREXEC(heurExecNlpdiving)
1629 { /*lint --e{715}*/
1630 SCIP_HEURDATA* heurdata;
1631 SCIP_NLPSOLSTAT nlpsolstat;
1632 SCIP_LPSOLSTAT lpsolstat;
1633 SCIP_SOL* nlpstartsol;
1634 SCIP_SOL* bestsol;
1635 SCIP_VAR** nlpcands;
1636 SCIP_VAR** covervars;
1637 SCIP_Real* nlpcandssol;
1638 SCIP_Real* nlpcandsfrac;
1639 SCIP_Real* pseudocandslpsol;
1640 SCIP_Real* pseudocandsnlpsol;
1641 SCIP_HASHMAP* varincover;
1642 SCIP_Real searchubbound;
1643 SCIP_Real searchavgbound;
1644 SCIP_Real searchbound;
1645 SCIP_Real objval;
1646 SCIP_Real oldobjval;
1647 SCIP_Real fixquot;
1648 SCIP_Real bestboundval;
1649 SCIP_Real timelim;
1650 SCIP_Bool bestcandmayround;
1651 SCIP_Bool bestcandroundup;
1652 SCIP_Bool nlperror;
1653 SCIP_Bool lperror;
1654 SCIP_Bool cutoff;
1655 SCIP_Bool backtracked;
1656 SCIP_Bool solvenlp;
1657 SCIP_Bool covercomputed;
1658 SCIP_Bool solvesubmip;
1659 SCIP_Longint ncalls;
1660 SCIP_Longint nsolsfound;
1661 int avgnnlpiterations;
1662 int maxnnlpiterations;
1663 int npseudocands;
1664 int nlpbranchcands;
1665 int ncovervars;
1666 int nnlpcands;
1667 int startnnlpcands;
1668 int depth;
1669 int maxdepth;
1670 int maxdivedepth;
1671 int divedepth;
1672 int lastnlpsolvedepth;
1673 int nfeasnlps;
1674 int bestcand;
1675 int origiterlim;
1676 int origfastfail;
1677 int c;
1678 int backtrackdepth; /* depth where to go when backtracking */
1679 SCIP_VAR* backtrackvar; /* (first) variable to fix differently in backtracking */
1680 SCIP_Real backtrackvarval; /* (fractional) value of backtrack variable */
1681 SCIP_Bool backtrackroundup; /* whether variable should be rounded up in backtracking */
1682
1683 backtrackdepth = -1;
1684 backtrackvar = NULL;
1685 backtrackvarval = 0.0;
1686 backtrackroundup = FALSE;
1687 bestsol = NULL;
1688 pseudocandsnlpsol = NULL;
1689 pseudocandslpsol = NULL;
1690 covervars = NULL;
1691
1692 assert(heur != NULL);
1693 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1694 assert(scip != NULL);
1695 assert(result != NULL);
1696 /* assert(SCIPhasCurrentNodeLP(scip)); */
1697
1698 *result = SCIP_DIDNOTRUN;
1699
1700 /* do not call heuristic of node was already detected to be infeasible */
1701 if( nodeinfeasible )
1702 return SCIP_OKAY;
1703
1704 /* only call heuristic, if an NLP relaxation has been constructed */
1705 if( !SCIPisNLPConstructed(scip) || SCIPgetNNlpis(scip) == 0 )
1706 return SCIP_OKAY;
1707
1708 /* only call heuristic, if the current node will not be cutoff, e.g., due to a (integer and NLP-)feasible LP solution */
1709 if( SCIPisFeasGE(scip, SCIPgetLocalLowerbound(scip), SCIPgetUpperbound(scip)) )
1710 return SCIP_OKAY;
1711
1712 /* get heuristic's data */
1713 heurdata = SCIPheurGetData(heur);
1714 assert(heurdata != NULL);
1715
1716 /* do not call heuristic, if it barely succeded */
1717 if( (SCIPheurGetNSolsFound(heur) + 1.0) / (SCIP_Real)(SCIPheurGetNCalls(heur) + 1.0) < heurdata->minsuccquot )
1718 return SCIP_OKAY;
1719
1720 *result = SCIP_DELAYED;
1721
1722 /* don't dive two times at the same node */
1723 if( SCIPgetLastDivenode(scip) == SCIPgetNNodes(scip) && SCIPgetDepth(scip) > 0 )
1724 return SCIP_OKAY;
1725
1726 *result = SCIP_DIDNOTRUN;
1727
1728 /* only try to dive, if we are in the correct part of the tree, given by minreldepth and maxreldepth */
1729 depth = SCIPgetDepth(scip);
1730 maxdepth = SCIPgetMaxDepth(scip);
1731 maxdepth = MAX(maxdepth, 30);
1732 if( depth < heurdata->minreldepth*maxdepth || depth > heurdata->maxreldepth*maxdepth )
1733 return SCIP_OKAY;
1734
1735 /* calculate the maximal number of NLP iterations until heuristic is aborted
1736 * maximal number is maxnlpiterabs plus a success-depending multiplier of maxnlpiterrel
1737 */
1738 ncalls = SCIPheurGetNCalls(heur);
1739 nsolsfound = 10*SCIPheurGetNBestSolsFound(heur) + heurdata->nsuccess;
1740 maxnnlpiterations = heurdata->maxnlpiterabs;
1741 maxnnlpiterations += (int)((1.0 + 10.0*(nsolsfound+1.0)/(ncalls+1.0)) * heurdata->maxnlpiterrel);
1742
1743 /* don't try to dive, if we took too many NLP iterations during diving */
1744 if( heurdata->nnlpiterations >= maxnnlpiterations )
1745 return SCIP_OKAY;
1746
1747 /* allow at least a bit more than the so far average number of NLP iterations per dive */
1748 avgnnlpiterations = (int)(heurdata->nnlpiterations / MAX(ncalls, 1.0));
1749 maxnnlpiterations = (int)MAX((SCIP_Real) maxnnlpiterations, (SCIP_Real) heurdata->nnlpiterations + 1.2*avgnnlpiterations);
1750
1751 /* don't try to dive, if there are no unfixed discrete variables */
1752 SCIP_CALL( SCIPgetPseudoBranchCands(scip, NULL, &npseudocands, NULL) );
1753 if( npseudocands == 0 )
1754 return SCIP_OKAY;
1755
1756 /* set time limit for NLP solver */
1757 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelim) );
1758 if( !SCIPisInfinity(scip, timelim) )
1759 timelim -= SCIPgetSolvingTime(scip);
1760 /* possibly exit if time is up (need to check here, since the paramter has to be >= 0) */
1761 if ( timelim <= 0.0 )
1762 return SCIP_OKAY;
1763 SCIP_CALL( SCIPsetNLPRealPar(scip, SCIP_NLPPAR_TILIM, timelim) );
1764
1765 *result = SCIP_DIDNOTFIND;
1766
1767 #ifdef SCIP_DEBUG
1768 /* SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_VERBLEVEL, 1) ); */
1769 #endif
1770
1771 /* set iteration limit */
1772 SCIP_CALL( SCIPgetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, &origiterlim) );
1773 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, maxnnlpiterations - heurdata->nnlpiterations) );
1774
1775 /* set whether NLP solver should fail fast */
1776 SCIP_CALL( SCIPgetNLPIntPar(scip, SCIP_NLPPAR_FASTFAIL, &origfastfail) );
1777 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_FASTFAIL, (int)heurdata->nlpfastfail) );
1778
1779 /* set starting point to lp solution */
1780 SCIP_CALL( SCIPsetNLPInitialGuessSol(scip, NULL) );
1781
1782 /* solve NLP relaxation, if not solved already */
1783 nlpsolstat = SCIPgetNLPSolstat(scip);
1784 if( nlpsolstat > SCIP_NLPSOLSTAT_FEASIBLE )
1785 {
1786 SCIP_NLPSTATISTICS* nlpstatistics;
1787
1788 SCIP_CALL( SCIPsolveNLP(scip) );
1789 SCIPstatistic( ++heurdata->nnlpsolves );
1790
1791 /* update iteration count */
1792 if( SCIPgetNLPTermstat(scip) < SCIP_NLPTERMSTAT_NUMERR )
1793 {
1794 SCIP_CALL( SCIPnlpStatisticsCreate(SCIPblkmem(scip), &nlpstatistics) );
1795 SCIP_CALL( SCIPgetNLPStatistics(scip, nlpstatistics) );
1796 heurdata->nnlpiterations += SCIPnlpStatisticsGetNIterations(nlpstatistics);
1797 SCIPnlpStatisticsFree(SCIPblkmem(scip), &nlpstatistics);
1798 }
1799
1800 nlpsolstat = SCIPgetNLPSolstat(scip);
1801
1802 /* give up, if no feasible solution found */
1803 if( nlpsolstat >= SCIP_NLPSOLSTAT_LOCINFEASIBLE )
1804 {
1805 SCIPdebugMsg(scip, "initial NLP infeasible or not solvable --> stop\n");
1806
1807 SCIPstatistic(
1808 if( SCIPgetNLPTermstat(scip) < SCIP_NLPTERMSTAT_NUMERR )
1809 heurdata->nfailcutoff++;
1810 else
1811 heurdata->nfailnlperror++;
1812 )
1813
1814 /* reset changed NLP parameters */
1815 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, origiterlim) );
1816 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_FASTFAIL, origfastfail) );
1817
1818 return SCIP_OKAY;
1819 }
1820 }
1821
1822 /* get NLP solution */
1823 SCIP_CALL( SCIPlinkNLPSol(scip, heurdata->sol) );
1824
1825 /* get fractional variables that should be integral */
1826 SCIP_CALL( getNLPFracVars(scip, heurdata, &nlpcands, &nlpcandssol, &nlpcandsfrac, &nnlpcands) );
1827 assert(nnlpcands <= npseudocands);
1828
1829 /* get LP candidates if LP solution is optimal */
1830 lpsolstat = SCIPgetLPSolstat(scip);
1831 if( lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
1832 nlpbranchcands = SCIPgetNLPBranchCands(scip);
1833 else
1834 nlpbranchcands = 0;
1835
1836 /* don't try to dive, if there are no fractional variables */
1837 if( nnlpcands == 0 )
1838 {
1839 SCIP_Bool success;
1840
1841 /* check, if solution was feasible and good enough
1842 *
1843 * Note that even if the NLP solver found a feasible solution it does not mean that is satisfy the integrality
1844 * conditions for fixed variables. This happens because the NLP solver uses relative tolerances for the bound
1845 * constraints but SCIP uses absolute tolerances for checking the integrality conditions.
1846 */
1847 #ifdef SCIP_DEBUG
1848 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, TRUE, TRUE, FALSE, TRUE, TRUE, &success) );
1849 #else
1850 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, TRUE, TRUE, &success) );
1851 #endif
1852 if( success )
1853 {
1854 SCIPdebugMsg(scip, " -> solution of first NLP was integral, feasible, and good enough\n");
1855 *result = SCIP_FOUNDSOL;
1856 }
1857
1858 /* reset changed NLP parameters */
1859 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, origiterlim) );
1860 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_FASTFAIL, origfastfail) );
1861
1862 return SCIP_OKAY;
1863 }
1864
1865 /* for guided diving: don't dive, if no feasible solutions exist */
1866 if( heurdata->varselrule == 'g' && SCIPgetNSols(scip) == 0 )
1867 return SCIP_OKAY;
1868
1869 /* for guided diving: get best solution that should guide the search; if this solution lives in the original variable space,
1870 * we cannot use it since it might violate the global bounds of the current problem
1871 */
1872 if( heurdata->varselrule == 'g' && SCIPsolIsOriginal(SCIPgetBestSol(scip)) )
1873 return SCIP_OKAY;
1874
1875 nlpstartsol = NULL;
1876 assert(nlpcandsfrac != NULL);
1877 assert(nlpcands != NULL);
1878 assert(nlpcandssol != NULL);
1879
1880 /* save solution of first NLP, if we may use it later */
1881 if( heurdata->nlpstart != 'n' )
1882 {
1883 SCIP_CALL( SCIPcreateNLPSol(scip, &nlpstartsol, heur) );
1884 SCIP_CALL( SCIPunlinkSol(scip, nlpstartsol) );
1885 }
1886
1887 /* calculate the objective search bound */
1888 if( SCIPgetNSolsFound(scip) == 0 )
1889 {
1890 if( heurdata->maxdiveubquotnosol > 0.0 )
1891 searchubbound = SCIPgetLowerbound(scip)
1892 + heurdata->maxdiveubquotnosol * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip));
1893 else
1894 searchubbound = SCIPinfinity(scip);
1895 if( heurdata->maxdiveavgquotnosol > 0.0 )
1896 searchavgbound = SCIPgetLowerbound(scip)
1897 + heurdata->maxdiveavgquotnosol * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip));
1898 else
1899 searchavgbound = SCIPinfinity(scip);
1900 }
1901 else
1902 {
1903 if( heurdata->maxdiveubquot > 0.0 )
1904 searchubbound = SCIPgetLowerbound(scip)
1905 + heurdata->maxdiveubquot * (SCIPgetCutoffbound(scip) - SCIPgetLowerbound(scip));
1906 else
1907 searchubbound = SCIPinfinity(scip);
1908 if( heurdata->maxdiveavgquot > 0.0 )
1909 searchavgbound = SCIPgetLowerbound(scip)
1910 + heurdata->maxdiveavgquot * (SCIPgetAvgLowerbound(scip) - SCIPgetLowerbound(scip));
1911 else
1912 searchavgbound = SCIPinfinity(scip);
1913 }
1914 searchbound = MIN(searchubbound, searchavgbound);
1915 if( SCIPisObjIntegral(scip) )
1916 searchbound = SCIPceil(scip, searchbound);
1917
1918 /* calculate the maximal diving depth: 10 * min{number of integer variables, max depth} */
1919 maxdivedepth = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
1920 maxdivedepth = MIN(maxdivedepth, maxdepth);
1921 maxdivedepth *= 10;
1922
1923 covercomputed = FALSE;
1924 varincover = NULL;
1925
1926 /* compute cover, if required */
1927 if( heurdata->prefercover || heurdata->solvesubmip )
1928 {
1929 SCIP_Real timelimit;
1930 SCIP_Real memorylimit;
1931
1932 /* get limits */
1933 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
1934 SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &memorylimit) );
1935 if( !SCIPisInfinity(scip, timelimit) )
1936 timelimit -= SCIPgetSolvingTime(scip);
1937
1938 /* substract the memory already used by the main SCIP and the estimated memory usage of external software */
1939 if( !SCIPisInfinity(scip, memorylimit) )
1940 {
1941 memorylimit -= SCIPgetMemUsed(scip)/1048576.0;
1942 memorylimit -= SCIPgetMemExternEstim(scip)/1048576.0;
1943 }
1944
1945 /* compute cover */
1946 ncovervars = -1;
1947 SCIP_CALL( SCIPallocBufferArray(scip, &covervars, SCIPgetNVars(scip)) );
1948 if( memorylimit > 2.0*SCIPgetMemExternEstim(scip)/1048576.0 && timelimit > 0.05 )
1949 {
1950 SCIP_CALL( SCIPcomputeCoverUndercover(scip, &ncovervars, covervars, timelimit, memorylimit, SCIPinfinity(scip), FALSE, FALSE, FALSE, 'u', &covercomputed) );
1951 }
1952
1953 if( covercomputed )
1954 {
1955 /* a cover can be empty, if the cover computation reveals that all nonlinear constraints are linear w.r.t. current variable fixations */
1956 assert(ncovervars >= 0);
1957
1958 /* create hash map */
1959 SCIP_CALL( SCIPhashmapCreate(&varincover, SCIPblkmem(scip), ncovervars) );
1960
1961 /* process variables in the cover */
1962 for( c = 0; c < ncovervars; c++ )
1963 {
1964 /* insert variable into hash map */
1965 if( SCIPvarGetType(covervars[c]) < SCIP_VARTYPE_IMPLINT )
1966 {
1967 assert(!SCIPhashmapExists(varincover, covervars[c]));
1968 SCIP_CALL( SCIPhashmapInsertInt(varincover, covervars[c], c+1) );
1969 }
1970
1971 /* catch bound change events of cover variables */
1972 assert(heurdata->eventhdlr != NULL);
1973 SCIP_CALL( SCIPcatchVarEvent(scip, covervars[c], SCIP_EVENTTYPE_BOUNDCHANGED, heurdata->eventhdlr,
1974 (SCIP_EVENTDATA*) heurdata, NULL) );
1975 assert(!SCIPisFeasEQ(scip, SCIPvarGetLbLocal(covervars[c]), SCIPvarGetUbLocal(covervars[c])));
1976 }
1977 }
1978 }
1979 else
1980 {
1981 covervars = NULL;
1982 ncovervars = 0;
1983 }
1984
1985 /* start diving */
1986 SCIP_CALL( SCIPstartProbing(scip) );
1987
1988 /* enables collection of variable statistics during probing */
1989 SCIPenableVarHistory(scip);
1990
1991 /* get NLP objective value*/
1992 objval = SCIPgetNLPObjval(scip);
1993
1994 SCIPdebugMsg(scip, "(node %" SCIP_LONGINT_FORMAT ") executing nlpdiving heuristic: depth=%d, %d fractionals, dualbound=%g, searchbound=%g\n",
1995 SCIPgetNNodes(scip), SCIPgetDepth(scip), nnlpcands, SCIPgetDualbound(scip), SCIPretransformObj(scip, searchbound));
1996
1997 /* store a copy of the best solution, if guided diving should be used */
1998 if( heurdata->varselrule == 'g' )
1999 {
2000 assert(SCIPgetNSols(scip) > 0);
2001 assert(!SCIPsolIsOriginal(SCIPgetBestSol(scip)));
2002
2003 SCIP_CALL( SCIPcreateSolCopy(scip, &bestsol, SCIPgetBestSol(scip)) );
2004 }
2005
2006 /* if double diving should be used, create arrays to hold to entire LP and NLP solution */
2007 if( heurdata->varselrule == 'd' )
2008 {
2009 SCIP_CALL( SCIPallocBufferArray(scip, &pseudocandslpsol, npseudocands) );
2010 SCIP_CALL( SCIPallocBufferArray(scip, &pseudocandsnlpsol, npseudocands) );
2011 }
2012
2013 /* dive as long we are in the given objective, depth and iteration limits and fractional variables exist, but
2014 * - if possible, we dive at least with the depth 10
2015 * - if the number of fractional variables decreased at least with 1 variable per 2 dive depths, we continue diving
2016 */
2017 nlperror = FALSE;
2018 lperror = FALSE;
2019 cutoff = FALSE;
2020 divedepth = 0;
2021 lastnlpsolvedepth = 0;
2022 backtracked = FALSE; /* whether we are in backtracking */
2023 fixquot = heurdata->fixquot;
2024 nfeasnlps = 1;
2025 startnnlpcands = nnlpcands;
2026 solvesubmip = heurdata->solvesubmip;
2027
2028 while( !nlperror && !cutoff && (nlpsolstat <= SCIP_NLPSOLSTAT_FEASIBLE || nlpsolstat == SCIP_NLPSOLSTAT_UNKNOWN) && nnlpcands > 0
2029 && (nfeasnlps < heurdata->maxfeasnlps
2030 || nnlpcands <= startnnlpcands - divedepth/2
2031 || (nfeasnlps < maxdivedepth && heurdata->nnlpiterations < maxnnlpiterations && objval < searchbound))
2032 && !SCIPisStopped(scip) )
2033 {
2034 SCIP_VAR* var;
2035 SCIP_Bool updatepscost;
2036
2037 /* open a new probing node if this will not exceed the maximal tree depth, otherwise stop here */
2038 if( SCIPgetDepth(scip) < SCIP_MAXTREEDEPTH )
2039 {
2040 SCIP_CALL( SCIPnewProbingNode(scip) );
2041 divedepth++;
2042 }
2043 else
2044 break;
2045
2046 bestcand = -1;
2047 bestcandmayround = TRUE;
2048 bestcandroundup = FALSE;
2049 bestboundval = SCIP_INVALID;
2050 updatepscost = TRUE;
2051 var = NULL;
2052
2053 /* find best candidate variable */
2054 switch( heurdata->varselrule )
2055 {
2056 case 'c':
2057 SCIP_CALL( chooseCoefVar(scip, heurdata, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, varincover, covercomputed,
2058 &bestcand, &bestcandmayround, &bestcandroundup) );
2059 if( bestcand >= 0 )
2060 {
2061 var = nlpcands[bestcand];
2062 bestboundval = nlpcandssol[bestcand];
2063 }
2064 break;
2065 case 'v':
2066 SCIP_CALL( chooseVeclenVar(scip, heurdata, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, varincover, covercomputed,
2067 &bestcand, &bestcandmayround, &bestcandroundup) );
2068 if( bestcand >= 0 )
2069 {
2070 var = nlpcands[bestcand];
2071 bestboundval = nlpcandssol[bestcand];
2072 }
2073 break;
2074 case 'p':
2075 SCIP_CALL( choosePscostVar(scip, heurdata, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, varincover, covercomputed,
2076 &bestcand, &bestcandmayround, &bestcandroundup) );
2077 if( bestcand >= 0 )
2078 {
2079 var = nlpcands[bestcand];
2080 bestboundval = nlpcandssol[bestcand];
2081 }
2082 break;
2083 case 'g':
2084 SCIP_CALL( chooseGuidedVar(scip, heurdata, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, bestsol, varincover, covercomputed,
2085 &bestcand, &bestcandmayround, &bestcandroundup) );
2086 if( bestcand >= 0 )
2087 {
2088 var = nlpcands[bestcand];
2089 bestboundval = nlpcandssol[bestcand];
2090 }
2091 break;
2092 case 'd':
2093 /* double diving only works if we have both relaxations at hand, otherwise we fall back to fractional diving */
2094 if( lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
2095 {
2096 SCIP_VAR** pseudocands;
2097
2098 SCIP_CALL( SCIPgetPseudoBranchCands(scip, &pseudocands, &npseudocands, NULL) );
2099 assert(backtrackdepth > 0 || nnlpcands <= npseudocands);
2100 assert(SCIPgetNLPBranchCands(scip) <= npseudocands);
2101 SCIP_CALL( SCIPgetSolVals(scip, NULL, npseudocands, pseudocands, pseudocandslpsol) );
2102 SCIP_CALL( SCIPgetSolVals(scip, heurdata->sol, npseudocands, pseudocands, pseudocandsnlpsol) );
2103 SCIP_CALL( chooseDoubleVar(scip, heurdata, pseudocands, pseudocandsnlpsol, pseudocandslpsol, npseudocands,
2104 varincover, covercomputed, &bestcand, &bestboundval, &bestcandmayround, &bestcandroundup) );
2105 if( bestcand >= 0 )
2106 var = pseudocands[bestcand];
2107 break;
2108 }
2109 else
2110 updatepscost = FALSE;
2111 /*lint -fallthrough*/
2112 case 'f':
2113 SCIP_CALL( chooseFracVar(scip, heurdata, nlpcands, nlpcandssol, nlpcandsfrac, nnlpcands, varincover, covercomputed,
2114 &bestcand, &bestcandmayround, &bestcandroundup) );
2115 if( bestcand >= 0 )
2116 {
2117 var = nlpcands[bestcand];
2118 bestboundval = nlpcandssol[bestcand];
2119 }
2120 break;
2121 default:
2122 SCIPerrorMessage("invalid variable selection rule\n");
2123 return SCIP_INVALIDDATA;
2124 }
2125
2126 /* if all candidates are roundable, try to round the solution
2127 * if var == NULL (i.e., bestcand == -1), then all solution candidates are outside bounds
2128 * this should only happen if they are slightly outside bounds (i.e., still within feastol, relative tolerance),
2129 * but far enough out to be considered as fractional (within feastol, but using absolute tolerance)
2130 * in this case, we also try our luck with rounding
2131 */
2132 if( (var == NULL || bestcandmayround) && backtrackdepth == -1 )
2133 {
2134 SCIP_Bool success;
2135
2136 /* create solution from diving NLP and try to round it */
2137 SCIP_CALL( SCIProundSol(scip, heurdata->sol, &success) );
2138
2139 if( success )
2140 {
2141 SCIPdebugMsg(scip, "nlpdiving found roundable primal solution: obj=%g\n", SCIPgetSolOrigObj(scip, heurdata->sol));
2142
2143 /* try to add solution to SCIP */
2144 #ifdef SCIP_DEBUG
2145 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, TRUE, TRUE, FALSE, FALSE, TRUE, &success) );
2146 #else
2147 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, FALSE, TRUE, &success) );
2148 #endif
2149
2150 /* check, if solution was feasible and good enough */
2151 if( success )
2152 {
2153 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n");
2154 *result = SCIP_FOUNDSOL;
2155 }
2156 }
2157 }
2158
2159 /* if all variables have been found to be essentially integral (even though there is some numerical doubt, see comment above), then stop */
2160 if( var == NULL )
2161 break;
2162
2163 do
2164 {
2165 SCIP_Real frac;
2166 frac = SCIP_INVALID;
2167
2168 if( backtracked && backtrackdepth > 0 )
2169 {
2170 assert(backtrackvar != NULL);
2171
2172 /* if the variable is already fixed or if the solution value is outside the domain, numerical troubles may have
2173 * occured or variable was fixed by propagation while backtracking => Abort diving!
2174 */
2175 if( SCIPvarGetLbLocal(backtrackvar) >= SCIPvarGetUbLocal(backtrackvar) - 0.5 )
2176 {
2177 SCIPdebugMsg(scip, "Selected variable <%s> already fixed to [%g,%g] (solval: %.9f), diving aborted \n",
2178 SCIPvarGetName(backtrackvar), SCIPvarGetLbLocal(backtrackvar), SCIPvarGetUbLocal(backtrackvar), backtrackvarval);
2179 cutoff = TRUE;
2180 break;
2181 }
2182 if( SCIPisFeasLT(scip, backtrackvarval, SCIPvarGetLbLocal(backtrackvar)) || SCIPisFeasGT(scip, backtrackvarval, SCIPvarGetUbLocal(backtrackvar)) )
2183 {
2184 SCIPdebugMsg(scip, "selected variable's <%s> solution value is outside the domain [%g,%g] (solval: %.9f), diving aborted\n",
2185 SCIPvarGetName(backtrackvar), SCIPvarGetLbLocal(backtrackvar), SCIPvarGetUbLocal(backtrackvar), backtrackvarval);
2186 assert(backtracked);
2187 break;
2188 }
2189
2190 /* round backtrack variable up or down */
2191 if( backtrackroundup )
2192 {
2193 SCIPdebugMsg(scip, " dive %d/%d, NLP iter %d/%d: var <%s>, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
2194 divedepth, maxdivedepth, heurdata->nnlpiterations, maxnnlpiterations,
2195 SCIPvarGetName(backtrackvar), backtrackvarval, SCIPvarGetLbLocal(backtrackvar), SCIPvarGetUbLocal(backtrackvar),
2196 SCIPfeasCeil(scip, backtrackvarval), SCIPvarGetUbLocal(backtrackvar));
2197 SCIP_CALL( SCIPchgVarLbProbing(scip, backtrackvar, SCIPfeasCeil(scip, backtrackvarval)) );
2198 }
2199 else
2200 {
2201 SCIPdebugMsg(scip, " dive %d/%d, NLP iter %d/%d: var <%s>, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
2202 divedepth, maxdivedepth, heurdata->nnlpiterations, maxnnlpiterations,
2203 SCIPvarGetName(backtrackvar), backtrackvarval, SCIPvarGetLbLocal(backtrackvar), SCIPvarGetUbLocal(backtrackvar),
2204 SCIPvarGetLbLocal(backtrackvar), SCIPfeasFloor(scip, backtrackvarval));
2205 SCIP_CALL( SCIPchgVarUbProbing(scip, backtrackvar, SCIPfeasFloor(scip, backtrackvarval)) );
2206 }
2207
2208 /* forget about backtrack variable */
2209 backtrackdepth = -1;
2210
2211 /* for pseudo cost computation */
2212 bestcandroundup = backtrackroundup;
2213 frac = SCIPfrac(scip, backtrackvarval);
2214 var = backtrackvar;
2215 }
2216 else
2217 {
2218 assert(var != NULL);
2219
2220 /* if the variable is already fixed or if the solution value is outside the domain, numerical troubles may have
2221 * occured or variable was fixed by propagation while backtracking => Abort diving!
2222 */
2223 if( SCIPvarGetLbLocal(var) >= SCIPvarGetUbLocal(var) - 0.5 )
2224 {
2225 SCIPdebugMsg(scip, "Selected variable <%s> already fixed to [%g,%g] (solval: %.9f), diving aborted \n",
2226 SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), bestboundval);
2227 cutoff = TRUE;
2228 break;
2229 }
2230 if( SCIPisFeasLT(scip, bestboundval, SCIPvarGetLbLocal(var)) || SCIPisFeasGT(scip, bestboundval, SCIPvarGetUbLocal(var)) )
2231 {
2232 SCIPdebugMsg(scip, "selected variable's <%s> solution value is outside the domain [%g,%g] (solval: %.9f), diving aborted\n",
2233 SCIPvarGetName(var), SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var), bestboundval);
2234 assert(backtracked);
2235 break;
2236 }
2237
2238 /* apply rounding of best candidate */
2239 if( bestcandroundup == !backtracked )
2240 {
2241 /* round variable up */
2242 SCIPdebugMsg(scip, " dive %d/%d, NLP iter %d/%d: var <%s>, round=%u, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
2243 divedepth, maxdivedepth, heurdata->nnlpiterations, maxnnlpiterations,
2244 SCIPvarGetName(var), bestcandmayround,
2245 bestboundval, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var),
2246 SCIPfeasCeil(scip, bestboundval), SCIPvarGetUbLocal(var));
2247 SCIP_CALL( SCIPchgVarLbProbing(scip, var, SCIPfeasCeil(scip, bestboundval)) );
2248
2249 /* remember variable for backtracking, if we have none yet (e.g., we are just after NLP solve) or we are half way to the next NLP solve */
2250 if( backtrackdepth == -1 || (divedepth - lastnlpsolvedepth == (int)(MIN(fixquot * nnlpcands, nlpbranchcands)/2.0)) )
2251 {
2252 backtrackdepth = divedepth;
2253 backtrackvar = var;
2254 backtrackvarval = bestboundval;
2255 backtrackroundup = FALSE;
2256 }
2257 }
2258 else
2259 {
2260 /* round variable down */
2261 SCIPdebugMsg(scip, " dive %d/%d, NLP iter %d/%d: var <%s>, round=%u, sol=%g, oldbounds=[%g,%g], newbounds=[%g,%g]\n",
2262 divedepth, maxdivedepth, heurdata->nnlpiterations, maxnnlpiterations,
2263 SCIPvarGetName(var), bestcandmayround,
2264 bestboundval, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var),
2265 SCIPvarGetLbLocal(var), SCIPfeasFloor(scip, bestboundval));
2266 SCIP_CALL( SCIPchgVarUbProbing(scip, var, SCIPfeasFloor(scip, bestboundval)) );
2267
2268 /* remember variable for backtracking, if we have none yet (e.g., we are just after NLP solve) or we are half way to the next NLP solve */
2269 if( backtrackdepth == -1 || (divedepth - lastnlpsolvedepth == (int)(MIN(fixquot * nnlpcands, nlpbranchcands)/2.0)) )
2270 {
2271 backtrackdepth = divedepth;
2272 backtrackvar = var;
2273 backtrackvarval = bestboundval;
2274 backtrackroundup = TRUE;
2275 }
2276 }
2277
2278 /* for pseudo-cost computation */
2279 if( updatepscost && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL )
2280 {
2281 if( heurdata->varselrule == 'd' )
2282 {
2283 assert(pseudocandsnlpsol != NULL);
2284 assert(0 <= bestcand && bestcand < npseudocands);
2285 frac = SCIPfrac(scip, pseudocandsnlpsol[bestcand]);
2286 }
2287 else
2288 frac = nlpcandsfrac[bestcand];
2289 }
2290 }
2291
2292 /* apply domain propagation */
2293 SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, NULL) );
2294 if( cutoff )
2295 {
2296 SCIPdebugMsg(scip, " *** cutoff detected in propagation at level %d\n", SCIPgetProbingDepth(scip));
2297 }
2298
2299 /* if all variables in the cover are fixed or there is no fractional variable in the cover,
2300 * then solve a sub-MIP
2301 */
2302 if( !cutoff && solvesubmip && covercomputed &&
2303 (heurdata->nfixedcovervars == ncovervars ||
2304 (heurdata->nfixedcovervars >= (ncovervars+1)/2 && !SCIPhashmapExists(varincover, var))) )
2305 {
2306 int probingdepth;
2307
2308 solvesubmip = FALSE;
2309 probingdepth = SCIPgetProbingDepth(scip);
2310 assert(probingdepth >= 1);
2311 assert(covervars != NULL);
2312
2313 if( heurdata->nfixedcovervars != ncovervars )
2314 {
2315 /* fix all remaining cover variables */
2316 for( c = 0; c < ncovervars && !cutoff ; c++ )
2317 {
2318 SCIP_Real lb;
2319 SCIP_Real ub;
2320 lb = SCIPvarGetLbLocal(covervars[c]);
2321 ub = SCIPvarGetUbLocal(covervars[c]);
2322 if( !SCIPisFeasEQ(scip, lb, ub) )
2323 {
2324 SCIP_Real nlpsolval;
2325
2326 /* adopt lpsolval w.r.t. intermediate bound changes by propagation */
2327 nlpsolval = SCIPvarGetNLPSol(covervars[c]);
2328 nlpsolval = MIN(nlpsolval,ub);
2329 nlpsolval = MAX(nlpsolval,lb);
2330 assert(SCIPvarGetType(covervars[c]) >= SCIP_VARTYPE_IMPLINT || SCIPisFeasIntegral(scip, nlpsolval));
2331
2332 /* open a new probing node if this will not exceed the maximal tree depth,
2333 * otherwise fix all the remaining variables at the same probing node
2334 * @todo do we need a new probing node for each fixing? if one of these fixings leads to a cutoff
2335 * we backtrack to the last probing node before we started to fix the covervars (and we do
2336 * not solve the probing LP). thus, it would be less work load in SCIPendProbing
2337 * and SCIPbacktrackProbing.
2338 */
2339 if( SCIP_MAXTREEDEPTH > SCIPgetDepth(scip) )
2340 {
2341 SCIP_CALL( SCIPnewProbingNode(scip) );
2342 }
2343
2344 /* fix and propagate */
2345 assert(SCIPisLbBetter(scip, nlpsolval, lb, ub) || SCIPisUbBetter(scip, nlpsolval, lb, ub));
2346
2347 if( SCIPisLbBetter(scip, nlpsolval, lb, ub) )
2348 {
2349 SCIP_CALL( SCIPchgVarLbProbing(scip, covervars[c], nlpsolval) );
2350 /* if covervar was implicit integer and fractional, then nlpsolval may be below lower bound now, so adjust to new bound */
2351 nlpsolval = MAX(nlpsolval, SCIPvarGetLbLocal(covervars[c])); /*lint !e666*/
2352 }
2353 if( SCIPisUbBetter(scip, nlpsolval, lb, ub) )
2354 {
2355 SCIP_CALL( SCIPchgVarUbProbing(scip, covervars[c], nlpsolval) );
2356 }
2357
2358 SCIP_CALL( SCIPpropagateProbing(scip, 0, &cutoff, NULL) );
2359 }
2360 }
2361 }
2362
2363 /* solve sub-MIP or return to standard diving */
2364 if( cutoff )
2365 {
2366 SCIP_CALL( SCIPbacktrackProbing(scip, probingdepth) );
2367 }
2368 else
2369 {
2370 SCIP_Bool success;
2371 success = FALSE;
2372
2373 SCIP_CALL( solveSubMIP(scip, heur, covervars, ncovervars, &success));
2374 if( success )
2375 *result = SCIP_FOUNDSOL;
2376 backtracked = TRUE; /* to avoid backtracking */
2377 nnlpcands = 0; /* to force termination */
2378 cutoff = TRUE;
2379 }
2380 }
2381
2382 /* resolve the diving LP */
2383 if( !cutoff && !lperror && (heurdata->lp || heurdata->varselrule == 'd')
2384 && SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL && SCIPisLPSolBasic(scip) )
2385 {
2386 SCIP_CALL( SCIPsolveProbingLP(scip, 100, &lperror, &cutoff) );
2387
2388 /* get LP solution status, objective value, and fractional variables, that should be integral */
2389 lpsolstat = SCIPgetLPSolstat(scip);
2390 assert(cutoff || (lpsolstat != SCIP_LPSOLSTAT_OBJLIMIT && lpsolstat != SCIP_LPSOLSTAT_INFEASIBLE &&
2391 (lpsolstat != SCIP_LPSOLSTAT_OPTIMAL || SCIPisLT(scip, SCIPgetLPObjval(scip), SCIPgetCutoffbound(scip)))));
2392
2393 if( lpsolstat == SCIP_LPSOLSTAT_OPTIMAL )
2394 {
2395 nlpbranchcands = SCIPgetNLPBranchCands(scip);
2396
2397 /* get new objective value */
2398 oldobjval = objval;
2399 objval = SCIPgetLPObjval(scip);
2400
2401 /* update pseudo cost values */
2402 if( updatepscost && SCIPisGT(scip, objval, oldobjval) )
2403 {
2404 assert(frac != SCIP_INVALID); /*lint !e777*/
2405 if( bestcandroundup )
2406 {
2407 SCIP_CALL( SCIPupdateVarPseudocost(scip, var, 1.0-frac, objval - oldobjval, 1.0) );
2408 }
2409 else
2410 {
2411 SCIP_CALL( SCIPupdateVarPseudocost(scip, var, 0.0-frac, objval - oldobjval, 1.0) );
2412 }
2413 }
2414 }
2415 else
2416 {
2417 nlpbranchcands = 0;
2418 }
2419
2420 if( cutoff )
2421 {
2422 SCIPdebugMsg(scip, " *** cutoff detected in LP solving at level %d, lpsolstat = %d\n", SCIPgetProbingDepth(scip), lpsolstat);
2423 }
2424 }
2425 else
2426 lpsolstat = SCIP_LPSOLSTAT_NOTSOLVED;
2427
2428 /* check whether we want to solve the NLP, which is the case if
2429 * - we are in backtracking, or
2430 * - we have (actively) fixed/rounded fixquot*nnlpcands variables
2431 * - all fractional variables were rounded/fixed (due to fixing and domain propagation)
2432 */
2433 solvenlp = backtracked;
2434 if( !solvenlp && !cutoff )
2435 {
2436 solvenlp = (lastnlpsolvedepth < divedepth - fixquot * nnlpcands);
2437 if( !solvenlp )
2438 {
2439 /* check if fractional NLP variables are left (some may have been fixed by propagation) */
2440 for( c = 0; c < nnlpcands; ++c )
2441 {
2442 var = nlpcands[c];
2443 if( SCIPisLT(scip, nlpcandssol[c], SCIPvarGetLbLocal(var)) || SCIPisGT(scip, nlpcandssol[c], SCIPvarGetUbLocal(var)) )
2444 continue;
2445 else
2446 break;
2447 }
2448 if( c == nnlpcands )
2449 solvenlp = TRUE;
2450 }
2451 }
2452
2453 nlpsolstat = SCIP_NLPSOLSTAT_UNKNOWN;
2454
2455 /* resolve the diving NLP */
2456 if( !cutoff && solvenlp )
2457 {
2458 SCIP_NLPTERMSTAT termstat;
2459 SCIP_NLPSTATISTICS* nlpstatistics;
2460
2461 /* set iteration limit, allow at least MINNLPITER many iterations */
2462 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, MAX(maxnnlpiterations - heurdata->nnlpiterations, MINNLPITER)) );
2463
2464 /* set time limit for NLP solver */
2465 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelim) );
2466 if( !SCIPisInfinity(scip, timelim) )
2467 timelim = MAX(0.0, timelim-SCIPgetSolvingTime(scip));/*lint !e666*/
2468 SCIP_CALL( SCIPsetNLPRealPar(scip, SCIP_NLPPAR_TILIM, timelim) );
2469
2470 /* set start solution, if we are in backtracking (previous NLP solve was infeasible) */
2471 if( heurdata->nlpstart != 'n' && backtracked )
2472 {
2473 assert(nlpstartsol != NULL);
2474
2475 SCIPdebugMsg(scip, "setting NLP initial guess\n");
2476
2477 SCIP_CALL( SCIPsetNLPInitialGuessSol(scip, nlpstartsol) );
2478 }
2479
2480 SCIP_CALL( SCIPsolveNLP(scip) );
2481 SCIPstatistic( ++heurdata->nnlpsolves );
2482
2483 termstat = SCIPgetNLPTermstat(scip);
2484 if( termstat >= SCIP_NLPTERMSTAT_NUMERR )
2485 {
2486 if( termstat >= SCIP_NLPTERMSTAT_LICERR )
2487 {
2488 SCIPverbMessage(scip, SCIP_VERBLEVEL_MINIMAL, NULL,
2489 "Error while solving NLP in nlpdiving heuristic; NLP solve terminated with code <%d>\n", termstat);
2490 }
2491 nlperror = TRUE;
2492 break;
2493 }
2494
2495 /* update iteration count */
2496 SCIP_CALL( SCIPnlpStatisticsCreate(SCIPblkmem(scip), &nlpstatistics) );
2497 SCIP_CALL( SCIPgetNLPStatistics(scip, nlpstatistics) );
2498 heurdata->nnlpiterations += SCIPnlpStatisticsGetNIterations(nlpstatistics);
2499 SCIPnlpStatisticsFree(SCIPblkmem(scip), &nlpstatistics);
2500
2501 /* get NLP solution status, objective value, and fractional variables, that should be integral */
2502 nlpsolstat = SCIPgetNLPSolstat(scip);
2503 cutoff = (nlpsolstat > SCIP_NLPSOLSTAT_FEASIBLE);
2504
2505 if( cutoff )
2506 {
2507 SCIPdebugMsg(scip, " *** cutoff detected in NLP solving at level %d, nlpsolstat: %d\n", SCIPgetProbingDepth(scip), nlpsolstat);
2508 }
2509 else
2510 {
2511 SCIP_CALL( SCIPlinkNLPSol(scip, heurdata->sol) );
2512
2513 /* remember that we have solve NLP on this depth successfully */
2514 lastnlpsolvedepth = divedepth;
2515 /* forget previous backtrack variable, we will never go back to a depth before the current one */
2516 backtrackdepth = -1;
2517 /* store NLP solution for warmstarting, if nlpstart is 'f' */
2518 if( heurdata->nlpstart == 'f' )
2519 {
2520 assert(nlpstartsol != NULL);
2521
2522 /* copy NLP solution values into nlpstartsol, is there a better way to do this???? */
2523 SCIP_CALL( SCIPlinkNLPSol(scip, nlpstartsol) );
2524 SCIP_CALL( SCIPunlinkSol(scip, nlpstartsol) );
2525 }
2526 /* increase counter on number of NLP solves with feasible solution */
2527 ++nfeasnlps;
2528 }
2529 }
2530
2531 /* perform backtracking if a cutoff was detected */
2532 if( cutoff && !backtracked && heurdata->backtrack )
2533 {
2534 if( backtrackdepth == -1 )
2535 {
2536 /* backtrack one step */
2537 SCIPdebugMsg(scip, " *** cutoff detected at level %d - backtracking one step\n", SCIPgetProbingDepth(scip));
2538 SCIP_CALL( SCIPbacktrackProbing(scip, SCIPgetProbingDepth(scip)-1) );
2539
2540 /* after backtracking there has to be at least one open node without exceeding the maximal tree depth */
2541 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(scip));
2542
2543 SCIP_CALL( SCIPnewProbingNode(scip) );
2544 }
2545 else
2546 {
2547 /* if we have a stored a depth for backtracking, go there */
2548 SCIPdebugMsg(scip, " *** cutoff detected at level %d - backtracking to depth %d\n", SCIPgetProbingDepth(scip), backtrackdepth);
2549 SCIP_CALL( SCIPbacktrackProbing(scip, backtrackdepth-1) );
2550
2551 /* after backtracking there has to be at least one open node without exceeding the maximal tree depth */
2552 assert(SCIP_MAXTREEDEPTH > SCIPgetDepth(scip));
2553
2554 SCIP_CALL( SCIPnewProbingNode(scip) );
2555 divedepth = backtrackdepth;
2556
2557 /* do not update pseudocosts if backtracking by more than one level */
2558 updatepscost = FALSE;
2559
2560 /* in case, we are feasible after backtracking, fix less variables at once in continuing diving
2561 * @todo should we remember the fixquot in heurdata for the next run?
2562 */
2563 fixquot *= 0.5;
2564 }
2565 /* remember that we are backtracking now */
2566 backtracked = TRUE;
2567 }
2568 else
2569 backtracked = FALSE;
2570 }
2571 while( backtracked );
2572
2573 if( !nlperror && !cutoff && nlpsolstat <= SCIP_NLPSOLSTAT_FEASIBLE )
2574 {
2575 /* get new fractional variables */
2576 SCIP_CALL( getNLPFracVars(scip, heurdata, &nlpcands, &nlpcandssol, &nlpcandsfrac, &nnlpcands) );
2577 }
2578 SCIPdebugMsg(scip, " -> nlpsolstat=%d, objval=%g/%g, nfrac nlp=%d lp=%d\n", nlpsolstat, objval, searchbound, nnlpcands, nlpbranchcands);
2579 }
2580
2581 /*lint --e{774}*/
2582 SCIPdebugMsg(scip, "NLP nlpdiving ABORT due to ");
2583 if( nlperror || (nlpsolstat > SCIP_NLPSOLSTAT_LOCINFEASIBLE && nlpsolstat != SCIP_NLPSOLSTAT_UNKNOWN) )
2584 {
2585 SCIPdebugMsgPrint(scip, "NLP bad status - nlperror: %ud nlpsolstat: %d \n", nlperror, nlpsolstat);
2586 SCIPstatistic( heurdata->nfailnlperror++ );
2587 }
2588 else if( SCIPisStopped(scip) || cutoff )
2589 {
2590 SCIPdebugMsgPrint(scip, "LIMIT hit - stop: %ud cutoff: %ud \n", SCIPisStopped(scip), cutoff);
2591 SCIPstatistic( heurdata->nfailcutoff++ );
2592 }
2593 else if(! (divedepth < 10
2594 || nnlpcands <= startnnlpcands - divedepth/2
2595 || (divedepth < maxdivedepth && heurdata->nnlpiterations < maxnnlpiterations && objval < searchbound) ) )
2596 {
2597 SCIPdebugMsgPrint(scip, "TOO DEEP - divedepth: %4d cands halfed: %d ltmaxdepth: %d ltmaxiter: %d bound: %d\n", divedepth,
2598 (nnlpcands > startnnlpcands - divedepth/2), (divedepth >= maxdivedepth), (heurdata->nnlpiterations >= maxnnlpiterations),
2599 (objval >= searchbound));
2600 SCIPstatistic( heurdata->nfaildepth++ );
2601 }
2602 else if( nnlpcands == 0 && !nlperror && !cutoff && nlpsolstat <= SCIP_NLPSOLSTAT_FEASIBLE )
2603 {
2604 SCIPdebugMsgPrint(scip, "SUCCESS\n");
2605 }
2606 else
2607 {
2608 SCIPdebugMsgPrint(scip, "UNKNOWN, very mysterical reason\n"); /* see also special case var == NULL (bestcand == -1) after choose*Var above */
2609 }
2610
2611 /* check if a solution has been found */
2612 if( nnlpcands == 0 && !nlperror && !cutoff && nlpsolstat <= SCIP_NLPSOLSTAT_FEASIBLE )
2613 {
2614 SCIP_Bool success;
2615
2616 /* create solution from diving NLP */
2617 SCIPdebugMsg(scip, "nlpdiving found primal solution: obj=%g\n", SCIPgetSolOrigObj(scip, heurdata->sol));
2618
2619 /* try to add solution to SCIP
2620 *
2621 * Note that even if the NLP solver found a feasible solution it does not mean that is satisfy the integrality
2622 * conditions for fixed variables. This happens because the NLP solver uses relative tolerances for the bound
2623 * constraints but SCIP uses absolute tolerances for checking the integrality conditions.
2624 */
2625 #ifdef SCIP_DEBUG
2626 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, TRUE, TRUE, FALSE, TRUE, TRUE, &success) );
2627 #else
2628 SCIP_CALL( SCIPtrySol(scip, heurdata->sol, FALSE, FALSE, FALSE, TRUE, TRUE, &success) );
2629 #endif
2630
2631 /* check, if solution was feasible and good enough */
2632 if( success )
2633 {
2634 SCIPdebugMsg(scip, " -> solution was feasible and good enough\n");
2635 *result = SCIP_FOUNDSOL;
2636 }
2637 else
2638 {
2639 SCIPdebugMsg(scip, " -> solution was not accepted\n");
2640 }
2641 }
2642
2643 /* end diving */
2644 SCIP_CALL( SCIPendProbing(scip) );
2645
2646 /* free hash map and drop variable bound change events */
2647 if( covercomputed )
2648 {
2649 assert(heurdata->eventhdlr != NULL);
2650 assert(heurdata->nfixedcovervars >= 0); /* variables might have been globally fixed in propagation */
2651 assert(varincover != NULL);
2652 assert(covervars != NULL);
2653
2654 SCIPhashmapFree(&varincover);
2655
2656 /* drop bound change events of cover variables */
2657 for( c = 0; c < ncovervars; c++ )
2658 {
2659 SCIP_CALL( SCIPdropVarEvent(scip, covervars[c], SCIP_EVENTTYPE_BOUNDCHANGED, heurdata->eventhdlr, (SCIP_EVENTDATA*)heurdata, -1) );
2660 }
2661 }
2662 else
2663 assert(varincover == NULL);
2664
2665 /* free NLP start solution */
2666 if( nlpstartsol != NULL )
2667 {
2668 SCIP_CALL( SCIPfreeSol(scip, &nlpstartsol) );
2669 }
2670
2671 /* reset changed NLP parameters */
2672 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_ITLIM, origiterlim) );
2673 SCIP_CALL( SCIPsetNLPIntPar(scip, SCIP_NLPPAR_FASTFAIL, origfastfail) );
2674
2675 /* free copied best solution */
2676 if( heurdata->varselrule == 'g' )
2677 {
2678 assert(bestsol != NULL);
2679 SCIP_CALL( SCIPfreeSol(scip, &bestsol) );
2680 }
2681 else
2682 assert(bestsol == NULL);
2683
2684 /* free arrays of LP and NLP solution */
2685 if( heurdata->varselrule == 'd' )
2686 {
2687 assert(pseudocandsnlpsol != NULL);
2688 assert(pseudocandsnlpsol != NULL);
2689 SCIPfreeBufferArray(scip, &pseudocandsnlpsol);
2690 SCIPfreeBufferArray(scip, &pseudocandslpsol);
2691 }
2692 else
2693 {
2694 assert(pseudocandsnlpsol == NULL);
2695 assert(pseudocandsnlpsol == NULL);
2696 }
2697
2698 /* free array of cover variables */
2699 if( heurdata->prefercover || heurdata->solvesubmip )
2700 {
2701 assert(covervars != NULL || !covercomputed);
2702 if( covervars != NULL )
2703 SCIPfreeBufferArray(scip, &covervars);
2704 }
2705 else
2706 assert(covervars == NULL);
2707
2708 if( *result == SCIP_FOUNDSOL )
2709 heurdata->nsuccess++;
2710
2711 SCIPdebugMsg(scip, "nlpdiving heuristic finished\n");
2712
2713 return SCIP_OKAY; /*lint !e438*/
2714 }
2715
2716
2717 /*
2718 * heuristic specific interface methods
2719 */
2720
2721 /** creates the nlpdiving heuristic and includes it in SCIP */
SCIPincludeHeurNlpdiving(SCIP * scip)2722 SCIP_RETCODE SCIPincludeHeurNlpdiving(
2723 SCIP* scip /**< SCIP data structure */
2724 )
2725 {
2726 SCIP_HEURDATA* heurdata;
2727 SCIP_HEUR* heur = NULL;
2728
2729 /* create heuristic data */
2730 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
2731
2732 /* include heuristic */
2733 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur, HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
2734 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecNlpdiving, heurdata) );
2735
2736 assert(heur != NULL);
2737 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyNlpdiving) );
2738 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeNlpdiving) );
2739 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitNlpdiving) );
2740 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitNlpdiving) );
2741 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolNlpdiving) );
2742
2743 /* get event handler for bound change events */
2744 heurdata->eventhdlr = NULL;
2745 /* create event handler for bound change events */
2746 SCIP_CALL( SCIPincludeEventhdlrBasic(scip, &heurdata->eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC,
2747 eventExecNlpdiving, NULL) );
2748 if ( heurdata->eventhdlr == NULL )
2749 {
2750 SCIPerrorMessage("event handler for " HEUR_NAME " heuristic not found.\n");
2751 return SCIP_PLUGINNOTFOUND;
2752 }
2753
2754 /* nlpdiving heuristic parameters */
2755 SCIP_CALL( SCIPaddRealParam(scip,
2756 "heuristics/" HEUR_NAME "/minreldepth",
2757 "minimal relative depth to start diving",
2758 &heurdata->minreldepth, TRUE, DEFAULT_MINRELDEPTH, 0.0, 1.0, NULL, NULL) );
2759 SCIP_CALL( SCIPaddRealParam(scip,
2760 "heuristics/" HEUR_NAME "/maxreldepth",
2761 "maximal relative depth to start diving",
2762 &heurdata->maxreldepth, TRUE, DEFAULT_MAXRELDEPTH, 0.0, 1.0, NULL, NULL) );
2763 SCIP_CALL( SCIPaddIntParam(scip,
2764 "heuristics/" HEUR_NAME "/maxnlpiterabs",
2765 "minimial absolute number of allowed NLP iterations",
2766 &heurdata->maxnlpiterabs, FALSE, DEFAULT_MAXNLPITERABS, 0, INT_MAX, NULL, NULL) );
2767 SCIP_CALL( SCIPaddIntParam(scip,
2768 "heuristics/" HEUR_NAME "/maxnlpiterrel",
2769 "additional allowed number of NLP iterations relative to successfully found solutions",
2770 &heurdata->maxnlpiterrel, FALSE, DEFAULT_MAXNLPITERREL, 0, INT_MAX, NULL, NULL) );
2771 SCIP_CALL( SCIPaddRealParam(scip,
2772 "heuristics/" HEUR_NAME "/maxdiveubquot",
2773 "maximal quotient (curlowerbound - lowerbound)/(cutoffbound - lowerbound) where diving is performed (0.0: no limit)",
2774 &heurdata->maxdiveubquot, TRUE, DEFAULT_MAXDIVEUBQUOT, 0.0, 1.0, NULL, NULL) );
2775 SCIP_CALL( SCIPaddRealParam(scip,
2776 "heuristics/" HEUR_NAME "/maxdiveavgquot",
2777 "maximal quotient (curlowerbound - lowerbound)/(avglowerbound - lowerbound) where diving is performed (0.0: no limit)",
2778 &heurdata->maxdiveavgquot, TRUE, DEFAULT_MAXDIVEAVGQUOT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2779 SCIP_CALL( SCIPaddRealParam(scip,
2780 "heuristics/" HEUR_NAME "/maxdiveubquotnosol",
2781 "maximal UBQUOT when no solution was found yet (0.0: no limit)",
2782 &heurdata->maxdiveubquotnosol, TRUE, DEFAULT_MAXDIVEUBQUOTNOSOL, 0.0, 1.0, NULL, NULL) );
2783 SCIP_CALL( SCIPaddRealParam(scip,
2784 "heuristics/" HEUR_NAME "/maxdiveavgquotnosol",
2785 "maximal AVGQUOT when no solution was found yet (0.0: no limit)",
2786 &heurdata->maxdiveavgquotnosol, TRUE, DEFAULT_MAXDIVEAVGQUOTNOSOL, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2787 SCIP_CALL( SCIPaddIntParam(scip,
2788 "heuristics/" HEUR_NAME "/maxfeasnlps",
2789 "maximal number of NLPs with feasible solution to solve during one dive",
2790 &heurdata->maxfeasnlps, FALSE, DEFAULT_MAXFEASNLPS, 1, INT_MAX, NULL, NULL) );
2791 SCIP_CALL( SCIPaddBoolParam(scip,
2792 "heuristics/" HEUR_NAME "/backtrack",
2793 "use one level of backtracking if infeasibility is encountered?",
2794 &heurdata->backtrack, FALSE, DEFAULT_BACKTRACK, NULL, NULL) );
2795 SCIP_CALL( SCIPaddBoolParam(scip,
2796 "heuristics/" HEUR_NAME "/lp",
2797 "should the LP relaxation be solved before the NLP relaxation?",
2798 &heurdata->lp, TRUE, DEFAULT_LP, NULL, NULL) );
2799 SCIP_CALL( SCIPaddBoolParam(scip,
2800 "heuristics/" HEUR_NAME "/preferlpfracs",
2801 "prefer variables that are also fractional in LP solution?",
2802 &heurdata->preferlpfracs, TRUE, DEFAULT_PREFERLPFRACS, NULL, NULL) );
2803 SCIP_CALL( SCIPaddRealParam(scip,
2804 "heuristics/" HEUR_NAME "/minsuccquot",
2805 "heuristic will not run if less then this percentage of calls succeeded (0.0: no limit)",
2806 &heurdata->minsuccquot, FALSE, DEFAULT_MINSUCCQUOT, 0.0, 1.0, NULL, NULL) );
2807 SCIP_CALL( SCIPaddRealParam(scip,
2808 "heuristics/" HEUR_NAME "/fixquot",
2809 "percentage of fractional variables that should be fixed before the next NLP solve",
2810 &heurdata->fixquot, FALSE, DEFAULT_FIXQUOT, 0.0, 1.0, NULL, NULL) );
2811 SCIP_CALL( SCIPaddBoolParam(scip,
2812 "heuristics/" HEUR_NAME "/prefercover",
2813 "should variables in a minimal cover be preferred?",
2814 &heurdata->prefercover, FALSE, DEFAULT_PREFERCOVER, NULL, NULL) );
2815 SCIP_CALL( SCIPaddBoolParam(scip,
2816 "heuristics/" HEUR_NAME "/solvesubmip",
2817 "should a sub-MIP be solved if all cover variables are fixed?",
2818 &heurdata->solvesubmip, FALSE, DEFAULT_SOLVESUBMIP, NULL, NULL) );
2819 SCIP_CALL( SCIPaddBoolParam(scip,
2820 "heuristics/" HEUR_NAME "/nlpfastfail",
2821 "should the NLP solver stop early if it converges slow?",
2822 &heurdata->nlpfastfail, FALSE, DEFAULT_NLPFASTFAIL, NULL, NULL) );
2823 SCIP_CALL( SCIPaddCharParam(scip,
2824 "heuristics/" HEUR_NAME "/nlpstart",
2825 "which point should be used as starting point for the NLP solver? ('n'one, last 'f'easible, from dive's'tart)",
2826 &heurdata->nlpstart, TRUE, DEFAULT_NLPSTART, "fns", NULL, NULL) );
2827 SCIP_CALL( SCIPaddCharParam(scip,
2828 "heuristics/" HEUR_NAME "/varselrule",
2829 "which variable selection should be used? ('f'ractionality, 'c'oefficient, 'p'seudocost, 'g'uided, 'd'ouble, 'v'eclen)",
2830 &heurdata->varselrule, FALSE, DEFAULT_VARSELRULE, "fcpgdv", NULL, NULL) );
2831
2832 return SCIP_OKAY;
2833 }
2834