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_alns.c
17 * @ingroup DEFPLUGINS_HEUR
18 * @brief Adaptive large neighborhood search heuristic that orchestrates popular LNS heuristics
19 * @author Gregor Hendel
20 */
21
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23
24 #include "blockmemshell/memory.h"
25 #include "scip/cons_linear.h"
26 #include "scip/heur_alns.h"
27 #include "scip/heuristics.h"
28 #include "scip/pub_bandit_epsgreedy.h"
29 #include "scip/pub_bandit_exp3.h"
30 #include "scip/pub_bandit.h"
31 #include "scip/pub_bandit_ucb.h"
32 #include "scip/pub_event.h"
33 #include "scip/pub_heur.h"
34 #include "scip/pub_message.h"
35 #include "scip/pub_misc.h"
36 #include "scip/pub_misc_select.h"
37 #include "scip/pub_sol.h"
38 #include "scip/pub_var.h"
39 #include "scip/scip_bandit.h"
40 #include "scip/scip_branch.h"
41 #include "scip/scip_cons.h"
42 #include "scip/scip_copy.h"
43 #include "scip/scip_event.h"
44 #include "scip/scip_general.h"
45 #include "scip/scip_heur.h"
46 #include "scip/scip_lp.h"
47 #include "scip/scip_mem.h"
48 #include "scip/scip_message.h"
49 #include "scip/scip_nodesel.h"
50 #include "scip/scip_numerics.h"
51 #include "scip/scip_param.h"
52 #include "scip/scip_prob.h"
53 #include "scip/scip_randnumgen.h"
54 #include "scip/scip_sol.h"
55 #include "scip/scip_solve.h"
56 #include "scip/scip_solvingstats.h"
57 #include "scip/scip_table.h"
58 #include "scip/scip_timing.h"
59 #include "scip/scip_tree.h"
60 #include "scip/scip_var.h"
61 #include <string.h>
62
63 #if !defined(_WIN32) && !defined(_WIN64)
64 #include <strings.h> /*lint --e{766}*/ /* needed for strncasecmp() */
65 #endif
66
67 #define HEUR_NAME "alns"
68 #define HEUR_DESC "Large neighborhood search heuristic that orchestrates the popular neighborhoods Local Branching, RINS, RENS, DINS etc."
69 #define HEUR_DISPCHAR SCIP_HEURDISPCHAR_LNS
70 #define HEUR_PRIORITY -1100500
71 #define HEUR_FREQ 20
72 #define HEUR_FREQOFS 0
73 #define HEUR_MAXDEPTH -1
74 #define HEUR_TIMING SCIP_HEURTIMING_AFTERNODE
75 #define HEUR_USESSUBSCIP TRUE /**< does the heuristic use a secondary SCIP instance? */
76
77 #define NNEIGHBORHOODS 9
78
79 /*
80 * limit parameters for sub-SCIPs
81 */
82 #define DEFAULT_NODESQUOT 0.1
83 #define DEFAULT_NODESOFFSET 500LL
84 #define DEFAULT_NSOLSLIM 3
85 #define DEFAULT_MINNODES 50LL
86 #define DEFAULT_MAXNODES 5000LL
87 #define DEFAULT_WAITINGNODES 25LL /**< number of nodes since last incumbent solution that the heuristic should wait */
88 #define DEFAULT_TARGETNODEFACTOR 1.05
89 #define LRATEMIN 0.01 /**< lower bound for learning rate for target nodes and minimum improvement */
90 #define LPLIMFAC 4.0
91
92 /*
93 * parameters for the minimum improvement
94 */
95 #define DEFAULT_MINIMPROVELOW 0.01
96 #define DEFAULT_MINIMPROVEHIGH 0.01
97 #define MINIMPROVEFAC 1.5
98 #define DEFAULT_STARTMINIMPROVE 0.01
99 #define DEFAULT_ADJUSTMINIMPROVE FALSE
100 #define DEFAULT_ADJUSTTARGETNODES TRUE /**< should the target nodes be dynamically adjusted? */
101
102 /*
103 * bandit algorithm parameters
104 */
105 #define DEFAULT_BESTSOLWEIGHT 1
106 #define DEFAULT_BANDITALGO 'u' /**< the default bandit algorithm: (u)pper confidence bounds, (e)xp.3, epsilon (g)reedy */
107 #define DEFAULT_REWARDCONTROL 0.8 /**< reward control to increase the weight of the simple solution indicator and decrease the weight of the closed gap reward */
108 #define DEFAULT_SCALEBYEFFORT TRUE /**< should the reward be scaled by the effort? */
109 #define DEFAULT_RESETWEIGHTS TRUE /**< should the bandit algorithms be reset when a new problem is read? */
110 #define DEFAULT_SUBSCIPRANDSEEDS FALSE /**< should random seeds of sub-SCIPs be altered to increase diversification? */
111 #define DEFAULT_REWARDBASELINE 0.5 /**< the reward baseline to separate successful and failed calls */
112 #define DEFAULT_FIXTOL 0.1 /**< tolerance by which the fixing rate may be missed without generic fixing */
113 #define DEFAULT_UNFIXTOL 0.1 /**< tolerance by which the fixing rate may be exceeded without generic unfixing */
114 #define DEFAULT_USELOCALREDCOST FALSE /**< should local reduced costs be used for generic (un)fixing? */
115 #define DEFAULT_BETA 0.0 /**< default reward offset between 0 and 1 at every observation for exp3 */
116
117 /*
118 * the following 3 parameters have been tuned by a simulation experiment
119 * as described in the paper.
120 */
121 #define DEFAULT_EPS 0.4685844 /**< increase exploration in epsilon-greedy bandit algorithm */
122 #define DEFAULT_ALPHA 0.0016 /**< parameter to increase the confidence width in UCB */
123 #define DEFAULT_GAMMA 0.07041455 /**< default weight between uniform (gamma ~ 1) and weight driven (gamma ~ 0) probability distribution for exp3 */
124 /*
125 * parameters to control variable fixing
126 */
127 #define DEFAULT_USEREDCOST TRUE /**< should reduced cost scores be used for variable priorization? */
128 #define DEFAULT_USEPSCOST TRUE /**< should pseudo cost scores be used for variable priorization? */
129 #define DEFAULT_USEDISTANCES TRUE /**< should distances from fixed variables be used for variable priorization */
130 #define DEFAULT_DOMOREFIXINGS TRUE /**< should the ALNS heuristic do more fixings by itself based on variable prioritization
131 * until the target fixing rate is reached? */
132 #define DEFAULT_ADJUSTFIXINGRATE TRUE /**< should the heuristic adjust the target fixing rate based on the success? */
133 #define FIXINGRATE_DECAY 0.75 /**< geometric decay for fixing rate adjustments */
134 #define FIXINGRATE_STARTINC 0.2 /**< initial increment value for fixing rate */
135 #define DEFAULT_USESUBSCIPHEURS FALSE /**< should the heuristic activate other sub-SCIP heuristics during its search? */
136 #define DEFAULT_COPYCUTS FALSE /**< should cutting planes be copied to the sub-SCIP? */
137 #define DEFAULT_REWARDFILENAME "-" /**< file name to store all rewards and the selection of the bandit */
138
139 /* individual random seeds */
140 #define DEFAULT_SEED 113
141 #define MUTATIONSEED 121
142 #define CROSSOVERSEED 321
143
144 /* individual neighborhood parameters */
145 #define DEFAULT_MINFIXINGRATE_RENS 0.3
146 #define DEFAULT_MAXFIXINGRATE_RENS 0.9
147 #define DEFAULT_ACTIVE_RENS TRUE
148 #define DEFAULT_PRIORITY_RENS 1.0
149
150 #define DEFAULT_MINFIXINGRATE_RINS 0.3
151 #define DEFAULT_MAXFIXINGRATE_RINS 0.9
152 #define DEFAULT_ACTIVE_RINS TRUE
153 #define DEFAULT_PRIORITY_RINS 1.0
154
155 #define DEFAULT_MINFIXINGRATE_MUTATION 0.3
156 #define DEFAULT_MAXFIXINGRATE_MUTATION 0.9
157 #define DEFAULT_ACTIVE_MUTATION TRUE
158 #define DEFAULT_PRIORITY_MUTATION 1.0
159
160 #define DEFAULT_MINFIXINGRATE_LOCALBRANCHING 0.3
161 #define DEFAULT_MAXFIXINGRATE_LOCALBRANCHING 0.9
162 #define DEFAULT_ACTIVE_LOCALBRANCHING TRUE
163 #define DEFAULT_PRIORITY_LOCALBRANCHING 1.0
164
165 #define DEFAULT_MINFIXINGRATE_PROXIMITY 0.3
166 #define DEFAULT_MAXFIXINGRATE_PROXIMITY 0.9
167 #define DEFAULT_ACTIVE_PROXIMITY TRUE
168 #define DEFAULT_PRIORITY_PROXIMITY 1.0
169
170 #define DEFAULT_MINFIXINGRATE_CROSSOVER 0.3
171 #define DEFAULT_MAXFIXINGRATE_CROSSOVER 0.9
172 #define DEFAULT_ACTIVE_CROSSOVER TRUE
173 #define DEFAULT_PRIORITY_CROSSOVER 1.0
174
175 #define DEFAULT_MINFIXINGRATE_ZEROOBJECTIVE 0.3
176 #define DEFAULT_MAXFIXINGRATE_ZEROOBJECTIVE 0.9
177 #define DEFAULT_ACTIVE_ZEROOBJECTIVE TRUE
178 #define DEFAULT_PRIORITY_ZEROOBJECTIVE 1.0
179
180 #define DEFAULT_MINFIXINGRATE_DINS 0.3
181 #define DEFAULT_MAXFIXINGRATE_DINS 0.9
182 #define DEFAULT_ACTIVE_DINS TRUE
183 #define DEFAULT_PRIORITY_DINS 1.0
184
185 #define DEFAULT_MINFIXINGRATE_TRUSTREGION 0.3
186 #define DEFAULT_MAXFIXINGRATE_TRUSTREGION 0.9
187 #define DEFAULT_ACTIVE_TRUSTREGION FALSE
188 #define DEFAULT_PRIORITY_TRUSTREGION 1.0
189
190
191 #define DEFAULT_NSOLS_CROSSOVER 2 /**< parameter for the number of solutions that crossover should combine */
192 #define DEFAULT_NPOOLSOLS_DINS 5 /**< number of pool solutions where binary solution values must agree */
193 #define DEFAULT_VIOLPENALTY_TRUSTREGION 100.0 /**< the penalty for violating the trust region */
194
195 /* event handler properties */
196 #define EVENTHDLR_NAME "Alns"
197 #define EVENTHDLR_DESC "LP event handler for " HEUR_NAME " heuristic"
198 #define SCIP_EVENTTYPE_ALNS (SCIP_EVENTTYPE_LPSOLVED | SCIP_EVENTTYPE_SOLFOUND | SCIP_EVENTTYPE_BESTSOLFOUND)
199
200 /* properties of the ALNS neighborhood statistics table */
201 #define TABLE_NAME_NEIGHBORHOOD "neighborhood"
202 #define TABLE_DESC_NEIGHBORHOOD "ALNS neighborhood statistics"
203 #define TABLE_POSITION_NEIGHBORHOOD 12500 /**< the position of the statistics table */
204 #define TABLE_EARLIEST_STAGE_NEIGHBORHOOD SCIP_STAGE_TRANSFORMED /**< output of the statistics table is only printed from this stage onwards */
205
206 /*
207 * Data structures
208 */
209
210 /*
211 * additional neighborhood data structures
212 */
213
214
215 typedef struct data_crossover DATA_CROSSOVER; /**< crossover neighborhood data structure */
216
217 typedef struct data_mutation DATA_MUTATION; /**< mutation neighborhood data structure */
218
219 typedef struct data_dins DATA_DINS; /**< dins neighborhood data structure */
220
221 typedef struct data_trustregion DATA_TRUSTREGION; /**< trustregion neighborhood data structure */
222
223 typedef struct NH_FixingRate NH_FIXINGRATE; /** fixing rate data structure */
224
225 typedef struct NH_Stats NH_STATS; /**< neighborhood statistics data structure */
226
227 typedef struct Nh NH; /**< neighborhood data structure */
228
229
230 /*
231 * variable priorization data structure for sorting
232 */
233 typedef struct VarPrio VARPRIO;
234
235 /** callback to collect variable fixings of neighborhood */
236 #define DECL_VARFIXINGS(x) SCIP_RETCODE x ( \
237 SCIP* scip, /**< SCIP data structure */ \
238 NH* neighborhood, /**< ALNS neighborhood data structure */ \
239 SCIP_VAR** varbuf, /**< buffer array to collect variables to fix */\
240 SCIP_Real* valbuf, /**< buffer array to collect fixing values */ \
241 int* nfixings, /**< pointer to store the number of fixings */ \
242 SCIP_RESULT* result /**< result pointer */ \
243 )
244
245 /** callback for subproblem changes other than variable fixings
246 *
247 * this callback can be used to further modify the subproblem by changes other than variable fixings.
248 * Typical modifications include restrictions of variable domains, the formulation of additional constraints,
249 * or changed objective coefficients.
250 *
251 * The callback should set the \p success pointer to indicate whether it was successful with its modifications or not.
252 */
253 #define DECL_CHANGESUBSCIP(x) SCIP_RETCODE x ( \
254 SCIP* sourcescip, /**< source SCIP data structure */\
255 SCIP* targetscip, /**< target SCIP data structure */\
256 NH* neighborhood, /**< ALNS neighborhood data structure */\
257 SCIP_VAR** subvars, /**< array of targetscip variables in the same order as the source SCIP variables */\
258 int* ndomchgs, /**< pointer to store the number of performed domain changes */\
259 int* nchgobjs, /**< pointer to store the number of changed objective coefficients */ \
260 int* naddedconss, /**< pointer to store the number of additional constraints */\
261 SCIP_Bool* success /**< pointer to store if the sub-MIP was successfully adjusted */\
262 )
263
264 /** optional initialization callback for neighborhoods when a new problem is read */
265 #define DECL_NHINIT(x) SCIP_RETCODE x ( \
266 SCIP* scip, /**< SCIP data structure */ \
267 NH* neighborhood /**< neighborhood data structure */ \
268 )
269
270 /** deinitialization callback for neighborhoods when exiting a problem */
271 #define DECL_NHEXIT(x) SCIP_RETCODE x ( \
272 SCIP* scip, /**< SCIP data structure */ \
273 NH* neighborhood /**< neighborhood data structure */ \
274 )
275
276 /** deinitialization callback for neighborhoods before SCIP is freed */
277 #define DECL_NHFREE(x) SCIP_RETCODE x ( \
278 SCIP* scip, /**< SCIP data structure */ \
279 NH* neighborhood /**< neighborhood data structure */ \
280 )
281
282 /** callback function to return a feasible reference solution for further fixings
283 *
284 * The reference solution should be stored in the \p solptr.
285 * The \p result pointer can be used to indicate either
286 *
287 * - SCIP_SUCCESS or
288 * - SCIP_DIDNOTFIND
289 */
290 #define DECL_NHREFSOL(x) SCIP_RETCODE x ( \
291 SCIP* scip, /**< SCIP data structure */ \
292 NH* neighborhood, /**< neighborhood data structure */ \
293 SCIP_SOL** solptr, /**< pointer to store the reference solution */ \
294 SCIP_RESULT* result /**< pointer to indicate the callback success whether a reference solution is available */ \
295 )
296
297 /** callback function to deactivate neighborhoods on problems where they are irrelevant */
298 #define DECL_NHDEACTIVATE(x) SCIP_RETCODE x (\
299 SCIP* scip, /**< SCIP data structure */ \
300 SCIP_Bool* deactivate /**< pointer to store whether the neighborhood should be deactivated (TRUE) for an instance */ \
301 )
302
303 /** sub-SCIP status code enumerator */
304 enum HistIndex
305 {
306 HIDX_OPT = 0, /**< sub-SCIP was solved to optimality */
307 HIDX_USR = 1, /**< sub-SCIP was user interrupted */
308 HIDX_NODELIM = 2, /**< sub-SCIP reached the node limit */
309 HIDX_STALLNODE = 3, /**< sub-SCIP reached the stall node limit */
310 HIDX_INFEAS = 4, /**< sub-SCIP was infeasible */
311 HIDX_SOLLIM = 5, /**< sub-SCIP reached the solution limit */
312 HIDX_OTHER = 6 /**< sub-SCIP reached none of the above codes */
313 };
314 typedef enum HistIndex HISTINDEX;
315 #define NHISTENTRIES 7
316
317
318 /** statistics for a neighborhood */
319 struct NH_Stats
320 {
321 SCIP_CLOCK* setupclock; /**< clock for sub-SCIP setup time */
322 SCIP_CLOCK* submipclock; /**< clock for the sub-SCIP solve */
323 SCIP_Longint usednodes; /**< total number of used nodes */
324 SCIP_Real oldupperbound; /**< upper bound before the sub-SCIP started */
325 SCIP_Real newupperbound; /**< new upper bound for allrewards mode to work correctly */
326 int nruns; /**< number of runs of a neighborhood */
327 int nrunsbestsol; /**< number of runs that produced a new incumbent */
328 SCIP_Longint nsolsfound; /**< the total number of solutions found */
329 SCIP_Longint nbestsolsfound; /**< the total number of improving solutions found */
330 int nfixings; /**< the number of fixings in one run */
331 int statushist[NHISTENTRIES]; /**< array to count sub-SCIP statuses */
332 };
333
334
335 /** fixing rate data structure to control the amount of target fixings of a neighborhood */
336 struct NH_FixingRate
337 {
338 SCIP_Real minfixingrate; /**< the minimum fixing rate */
339 SCIP_Real targetfixingrate; /**< the current target fixing rate */
340 SCIP_Real increment; /**< the current increment by which the target fixing rate is in-/decreased */
341 SCIP_Real maxfixingrate; /**< the maximum fixing rate */
342 };
343
344 /** neighborhood data structure with callbacks, statistics, fixing rate */
345 struct Nh
346 {
347 char* name; /**< the name of this neighborhood */
348 NH_FIXINGRATE fixingrate; /**< fixing rate for this neighborhood */
349 NH_STATS stats; /**< statistics for this neighborhood */
350 DECL_VARFIXINGS ((*varfixings)); /**< variable fixings callback for this neighborhood */
351 DECL_CHANGESUBSCIP ((*changesubscip)); /**< callback for subproblem changes other than variable fixings */
352 DECL_NHINIT ((*nhinit)); /**< initialization callback when a new problem is read */
353 DECL_NHEXIT ((*nhexit)); /**< deinitialization callback when exiting a problem */
354 DECL_NHFREE ((*nhfree)); /**< deinitialization callback before SCIP is freed */
355 DECL_NHREFSOL ((*nhrefsol)); /**< callback function to return a reference solution for further fixings, or NULL */
356 DECL_NHDEACTIVATE ((*nhdeactivate)); /**< callback function to deactivate neighborhoods on problems where they are irrelevant, or NULL if it is always active */
357 SCIP_Bool active; /**< is this neighborhood active or not? */
358 SCIP_Real priority; /**< positive call priority to initialize bandit algorithms */
359 union
360 {
361 DATA_MUTATION* mutation; /**< mutation data */
362 DATA_CROSSOVER* crossover; /**< crossover data */
363 DATA_DINS* dins; /**< dins data */
364 DATA_TRUSTREGION* trustregion; /**< trustregion data */
365 } data; /**< data object for neighborhood specific data */
366 };
367
368 /** mutation neighborhood data structure */
369 struct data_mutation
370 {
371 SCIP_RANDNUMGEN* rng; /**< random number generator */
372 };
373
374 /** crossover neighborhood data structure */
375 struct data_crossover
376 {
377 int nsols; /**< the number of solutions that crossover should combine */
378 SCIP_RANDNUMGEN* rng; /**< random number generator to draw from the solution pool */
379 SCIP_SOL* selsol; /**< best selected solution by crossover as reference point */
380 };
381
382 /** dins neighborhood data structure */
383 struct data_dins
384 {
385 int npoolsols; /**< number of pool solutions where binary solution values must agree */
386 };
387
388 struct data_trustregion
389 {
390 SCIP_Real violpenalty; /**< the penalty for violating the trust region */
391 };
392
393 /** primal heuristic data */
394 struct SCIP_HeurData
395 {
396 NH** neighborhoods; /**< array of neighborhoods */
397 SCIP_BANDIT* bandit; /**< bandit algorithm */
398 char* rewardfilename; /**< file name to store all rewards and the selection of the bandit */
399 FILE* rewardfile; /**< reward file pointer, or NULL */
400 SCIP_Longint nodesoffset; /**< offset added to the nodes budget */
401 SCIP_Longint maxnodes; /**< maximum number of nodes in a single sub-SCIP */
402 SCIP_Longint targetnodes; /**< targeted number of nodes to start a sub-SCIP */
403 SCIP_Longint minnodes; /**< minimum number of nodes required to start a sub-SCIP */
404 SCIP_Longint usednodes; /**< total number of nodes already spent in sub-SCIPs */
405 SCIP_Longint waitingnodes; /**< number of nodes since last incumbent solution that the heuristic should wait */
406 SCIP_Real nodesquot; /**< fraction of nodes compared to the main SCIP for budget computation */
407 SCIP_Real startminimprove; /**< initial factor by which ALNS should at least improve the incumbent */
408 SCIP_Real minimprovelow; /**< lower threshold for the minimal improvement over the incumbent */
409 SCIP_Real minimprovehigh; /**< upper bound for the minimal improvement over the incumbent */
410 SCIP_Real minimprove; /**< factor by which ALNS should at least improve the incumbent */
411 SCIP_Real lplimfac; /**< limit fraction of LPs per node to interrupt sub-SCIP */
412 SCIP_Real exp3_gamma; /**< weight between uniform (gamma ~ 1) and weight driven (gamma ~ 0) probability distribution for exp3 */
413 SCIP_Real exp3_beta; /**< reward offset between 0 and 1 at every observation for exp3 */
414 SCIP_Real epsgreedy_eps; /**< increase exploration in epsilon-greedy bandit algorithm */
415 SCIP_Real ucb_alpha; /**< parameter to increase the confidence width in UCB */
416 SCIP_Real rewardcontrol; /**< reward control to increase the weight of the simple solution indicator
417 * and decrease the weight of the closed gap reward */
418 SCIP_Real targetnodefactor; /**< factor by which target node number is eventually increased */
419 SCIP_Real rewardbaseline; /**< the reward baseline to separate successful and failed calls */
420 SCIP_Real fixtol; /**< tolerance by which the fixing rate may be missed without generic fixing */
421 SCIP_Real unfixtol; /**< tolerance by which the fixing rate may be exceeded without generic unfixing */
422 int nneighborhoods; /**< number of neighborhoods */
423 int nactiveneighborhoods;/**< number of active neighborhoods */
424 int ninitneighborhoods; /**< neighborhoods that were used at least one time */
425 int nsolslim; /**< limit on the number of improving solutions in a sub-SCIP call */
426 int seed; /**< initial random seed for bandit algorithms and random decisions by neighborhoods */
427 int currneighborhood; /**< index of currently selected neighborhood */
428 int ndelayedcalls; /**< the number of delayed calls */
429 char banditalgo; /**< the bandit algorithm: (u)pper confidence bounds, (e)xp.3, epsilon (g)reedy */
430 SCIP_Bool useredcost; /**< should reduced cost scores be used for variable prioritization? */
431 SCIP_Bool usedistances; /**< should distances from fixed variables be used for variable prioritization */
432 SCIP_Bool usepscost; /**< should pseudo cost scores be used for variable prioritization? */
433 SCIP_Bool domorefixings; /**< should the ALNS heuristic do more fixings by itself based on variable prioritization
434 * until the target fixing rate is reached? */
435 SCIP_Bool adjustfixingrate; /**< should the heuristic adjust the target fixing rate based on the success? */
436 SCIP_Bool usesubscipheurs; /**< should the heuristic activate other sub-SCIP heuristics during its search? */
437 SCIP_Bool adjustminimprove; /**< should the factor by which the minimum improvement is bound be dynamically updated? */
438 SCIP_Bool adjusttargetnodes; /**< should the target nodes be dynamically adjusted? */
439 SCIP_Bool resetweights; /**< should the bandit algorithms be reset when a new problem is read? */
440 SCIP_Bool subsciprandseeds; /**< should random seeds of sub-SCIPs be altered to increase diversification? */
441 SCIP_Bool scalebyeffort; /**< should the reward be scaled by the effort? */
442 SCIP_Bool copycuts; /**< should cutting planes be copied to the sub-SCIP? */
443 SCIP_Bool uselocalredcost; /**< should local reduced costs be used for generic (un)fixing? */
444 };
445
446 /** event handler data */
447 struct SCIP_EventData
448 {
449 SCIP_VAR** subvars; /**< the variables of the subproblem */
450 SCIP* sourcescip; /**< original SCIP data structure */
451 SCIP_HEUR* heur; /**< alns heuristic structure */
452 SCIP_Longint nodelimit; /**< node limit of the run */
453 SCIP_Real lplimfac; /**< limit fraction of LPs per node to interrupt sub-SCIP */
454 NH_STATS* runstats; /**< run statistics for the current neighborhood */
455 SCIP_Bool allrewardsmode; /**< true if solutions should only be checked for reward comparisons */
456 };
457
458 /** represents limits for the sub-SCIP solving process */
459 struct SolveLimits
460 {
461 SCIP_Longint nodelimit; /**< maximum number of solving nodes for the sub-SCIP */
462 SCIP_Real memorylimit; /**< memory limit for the sub-SCIP */
463 SCIP_Real timelimit; /**< time limit for the sub-SCIP */
464 SCIP_Longint stallnodes; /**< maximum number of nodes without (primal) stalling */
465 };
466
467 typedef struct SolveLimits SOLVELIMITS;
468
469 /** data structure that can be used for variable prioritization for additional fixings */
470 struct VarPrio
471 {
472 SCIP* scip; /**< SCIP data structure */
473 SCIP_Real* randscores; /**< random scores for prioritization */
474 int* distances; /**< breadth-first distances from already fixed variables */
475 SCIP_Real* redcostscores; /**< reduced cost scores for fixing a variable to a reference value */
476 SCIP_Real* pscostscores; /**< pseudocost scores for fixing a variable to a reference value */
477 unsigned int useredcost:1; /**< should reduced cost scores be used for variable prioritization? */
478 unsigned int usedistances:1; /**< should distances from fixed variables be used for variable prioritization */
479 unsigned int usepscost:1; /**< should pseudo cost scores be used for variable prioritization? */
480 };
481
482 /*
483 * Local methods
484 */
485
486 /** Reset target fixing rate */
487 static
resetFixingRate(SCIP * scip,NH_FIXINGRATE * fixingrate)488 SCIP_RETCODE resetFixingRate(
489 SCIP* scip, /**< SCIP data structure */
490 NH_FIXINGRATE* fixingrate /**< heuristic fixing rate */
491 )
492 {
493 assert(scip != NULL);
494 assert(fixingrate != NULL);
495 fixingrate->increment = FIXINGRATE_STARTINC;
496
497 /* always start with the most conservative value */
498 fixingrate->targetfixingrate = fixingrate->maxfixingrate;
499
500 return SCIP_OKAY;
501 }
502
503 /** reset the currently active neighborhood */
504 static
resetCurrentNeighborhood(SCIP_HEURDATA * heurdata)505 void resetCurrentNeighborhood(
506 SCIP_HEURDATA* heurdata
507 )
508 {
509 assert(heurdata != NULL);
510 heurdata->currneighborhood = -1;
511 heurdata->ndelayedcalls = 0;
512 }
513
514 /** update increment for fixing rate */
515 static
updateFixingRateIncrement(NH_FIXINGRATE * fx)516 void updateFixingRateIncrement(
517 NH_FIXINGRATE* fx /**< fixing rate */
518 )
519 {
520 fx->increment *= FIXINGRATE_DECAY;
521 fx->increment = MAX(fx->increment, LRATEMIN);
522 }
523
524
525 /** increase fixing rate
526 *
527 * decrease also the rate by which the target fixing rate is adjusted
528 */
529 static
increaseFixingRate(NH_FIXINGRATE * fx)530 void increaseFixingRate(
531 NH_FIXINGRATE* fx /**< fixing rate */
532 )
533 {
534 fx->targetfixingrate += fx->increment;
535 fx->targetfixingrate = MIN(fx->targetfixingrate, fx->maxfixingrate);
536 }
537
538 /** decrease fixing rate
539 *
540 * decrease also the rate by which the target fixing rate is adjusted
541 */
542 static
decreaseFixingRate(NH_FIXINGRATE * fx)543 void decreaseFixingRate(
544 NH_FIXINGRATE* fx /**< fixing rate */
545 )
546 {
547 fx->targetfixingrate -= fx->increment;
548 fx->targetfixingrate = MAX(fx->targetfixingrate, fx->minfixingrate);
549 }
550
551 /** update fixing rate based on the results of the current run */
552 static
updateFixingRate(NH * neighborhood,SCIP_STATUS subscipstatus,NH_STATS * runstats)553 void updateFixingRate(
554 NH* neighborhood, /**< neighborhood */
555 SCIP_STATUS subscipstatus, /**< status of the sub-SCIP run */
556 NH_STATS* runstats /**< run statistics for this run */
557 )
558 {
559 NH_FIXINGRATE* fx;
560
561 fx = &neighborhood->fixingrate;
562
563 switch (subscipstatus)
564 {
565 case SCIP_STATUS_OPTIMAL:
566 case SCIP_STATUS_INFEASIBLE:
567 case SCIP_STATUS_INFORUNBD:
568 case SCIP_STATUS_SOLLIMIT:
569 case SCIP_STATUS_BESTSOLLIMIT:
570 /* decrease the fixing rate (make subproblem harder) */
571 decreaseFixingRate(fx);
572 break;
573 case SCIP_STATUS_STALLNODELIMIT:
574 case SCIP_STATUS_USERINTERRUPT:
575 case SCIP_STATUS_TERMINATE:
576 case SCIP_STATUS_NODELIMIT:
577 /* increase the fixing rate (make the subproblem easier) only if no solution was found */
578 if( runstats->nbestsolsfound <= 0 )
579 increaseFixingRate(fx);
580 break;
581 /* fall through cases to please lint */
582 case SCIP_STATUS_UNKNOWN:
583 case SCIP_STATUS_TOTALNODELIMIT:
584 case SCIP_STATUS_TIMELIMIT:
585 case SCIP_STATUS_MEMLIMIT:
586 case SCIP_STATUS_GAPLIMIT:
587 case SCIP_STATUS_RESTARTLIMIT:
588 case SCIP_STATUS_UNBOUNDED:
589 default:
590 break;
591 }
592
593 updateFixingRateIncrement(fx);
594 }
595
596 /** increase target node limit */
597 static
increaseTargetNodeLimit(SCIP_HEURDATA * heurdata)598 void increaseTargetNodeLimit(
599 SCIP_HEURDATA* heurdata /**< heuristic data */
600 )
601 {
602 heurdata->targetnodes = (SCIP_Longint)(heurdata->targetnodes * heurdata->targetnodefactor) + 1;
603
604 /* respect upper and lower parametrized bounds on targetnodes */
605 if( heurdata->targetnodes > heurdata->maxnodes )
606 heurdata->targetnodes = heurdata->maxnodes;
607 }
608
609 /** reset target node limit */
610 static
resetTargetNodeLimit(SCIP_HEURDATA * heurdata)611 void resetTargetNodeLimit(
612 SCIP_HEURDATA* heurdata /**< heuristic data */
613 )
614 {
615 heurdata->targetnodes = heurdata->minnodes;
616 }
617
618 /** update target node limit based on the current run results */
619 static
updateTargetNodeLimit(SCIP_HEURDATA * heurdata,NH_STATS * runstats,SCIP_STATUS subscipstatus)620 void updateTargetNodeLimit(
621 SCIP_HEURDATA* heurdata, /**< heuristic data */
622 NH_STATS* runstats, /**< statistics of the run */
623 SCIP_STATUS subscipstatus /**< status of the sub-SCIP run */
624 )
625 {
626 switch (subscipstatus)
627 {
628 case SCIP_STATUS_STALLNODELIMIT:
629 case SCIP_STATUS_NODELIMIT:
630 /* the subproblem could be explored more */
631 if( runstats->nbestsolsfound == 0 )
632 increaseTargetNodeLimit(heurdata);
633 break;
634 case SCIP_STATUS_OPTIMAL:
635 case SCIP_STATUS_INFEASIBLE:
636 case SCIP_STATUS_INFORUNBD:
637 case SCIP_STATUS_SOLLIMIT:
638 case SCIP_STATUS_BESTSOLLIMIT:
639 break;
640 case SCIP_STATUS_USERINTERRUPT:
641 case SCIP_STATUS_TERMINATE:
642 case SCIP_STATUS_UNKNOWN:
643 case SCIP_STATUS_TOTALNODELIMIT:
644 case SCIP_STATUS_TIMELIMIT:
645 case SCIP_STATUS_MEMLIMIT:
646 case SCIP_STATUS_GAPLIMIT:
647 case SCIP_STATUS_RESTARTLIMIT:
648 case SCIP_STATUS_UNBOUNDED:
649 break;
650 default:
651 break;
652 }
653 }
654
655 /** reset the minimum improvement for the sub-SCIPs */
656 static
resetMinimumImprovement(SCIP_HEURDATA * heurdata)657 void resetMinimumImprovement(
658 SCIP_HEURDATA* heurdata /**< heuristic data */
659 )
660 {
661 assert(heurdata != NULL);
662 heurdata->minimprove = heurdata->startminimprove;
663 }
664
665 /** increase minimum improvement for the sub-SCIPs */
666 static
increaseMinimumImprovement(SCIP_HEURDATA * heurdata)667 void increaseMinimumImprovement(
668 SCIP_HEURDATA* heurdata /**< heuristic data */
669 )
670 {
671 assert(heurdata != NULL);
672
673 heurdata->minimprove *= MINIMPROVEFAC;
674 heurdata->minimprove = MIN(heurdata->minimprove, heurdata->minimprovehigh);
675 }
676
677 /** decrease the minimum improvement for the sub-SCIPs */
678 static
decreaseMinimumImprovement(SCIP_HEURDATA * heurdata)679 void decreaseMinimumImprovement(
680 SCIP_HEURDATA* heurdata /**< heuristic data */
681 )
682 {
683 assert(heurdata != NULL);
684
685 heurdata->minimprove /= MINIMPROVEFAC;
686 SCIPdebugMessage("%.4f", heurdata->minimprovelow);
687 heurdata->minimprove = MAX(heurdata->minimprove, heurdata->minimprovelow);
688 }
689
690 /** update the minimum improvement based on the status of the sub-SCIP */
691 static
updateMinimumImprovement(SCIP_HEURDATA * heurdata,SCIP_STATUS subscipstatus,NH_STATS * runstats)692 void updateMinimumImprovement(
693 SCIP_HEURDATA* heurdata, /**< heuristic data */
694 SCIP_STATUS subscipstatus, /**< status of the sub-SCIP run */
695 NH_STATS* runstats /**< run statistics for this run */
696 )
697 {
698 assert(heurdata != NULL);
699
700 /* if the sub-SCIP status was infeasible, we rather want to make the sub-SCIP easier
701 * with a smaller minimum improvement.
702 *
703 * If a solution limit was reached, we may, set it higher.
704 */
705 switch (subscipstatus)
706 {
707 case SCIP_STATUS_INFEASIBLE:
708 case SCIP_STATUS_INFORUNBD:
709 /* subproblem was infeasible, probably due to the minimum improvement -> decrease minimum improvement */
710 decreaseMinimumImprovement(heurdata);
711
712 break;
713 case SCIP_STATUS_SOLLIMIT:
714 case SCIP_STATUS_BESTSOLLIMIT:
715 case SCIP_STATUS_OPTIMAL:
716 /* subproblem could be optimally solved -> try higher minimum improvement */
717 increaseMinimumImprovement(heurdata);
718 break;
719 case SCIP_STATUS_NODELIMIT:
720 case SCIP_STATUS_STALLNODELIMIT:
721 case SCIP_STATUS_USERINTERRUPT:
722 /* subproblem was too hard, decrease minimum improvement */
723 if( runstats->nbestsolsfound <= 0 )
724 decreaseMinimumImprovement(heurdata);
725 break;
726 case SCIP_STATUS_UNKNOWN:
727 case SCIP_STATUS_TOTALNODELIMIT:
728 case SCIP_STATUS_TIMELIMIT:
729 case SCIP_STATUS_MEMLIMIT:
730 case SCIP_STATUS_GAPLIMIT:
731 case SCIP_STATUS_RESTARTLIMIT:
732 case SCIP_STATUS_UNBOUNDED:
733 case SCIP_STATUS_TERMINATE:
734 default:
735 break;
736 }
737 }
738
739 /** Reset neighborhood statistics */
740 static
neighborhoodStatsReset(SCIP * scip,NH_STATS * stats)741 SCIP_RETCODE neighborhoodStatsReset(
742 SCIP* scip, /**< SCIP data structure */
743 NH_STATS* stats /**< neighborhood statistics */
744 )
745 {
746 assert(scip != NULL);
747 assert(stats != NULL);
748
749 stats->nbestsolsfound = 0;
750 stats->nruns = 0;
751 stats->nrunsbestsol = 0;
752 stats->nsolsfound = 0;
753 stats->usednodes = 0L;
754 stats->nfixings = 0L;
755
756 BMSclearMemoryArray(stats->statushist, NHISTENTRIES);
757
758 SCIP_CALL( SCIPresetClock(scip, stats->setupclock) );
759 SCIP_CALL( SCIPresetClock(scip, stats->submipclock) );
760
761 return SCIP_OKAY;
762 }
763
764 /** create a neighborhood of the specified name and include it into the ALNS heuristic */
765 static
alnsIncludeNeighborhood(SCIP * scip,SCIP_HEURDATA * heurdata,NH ** neighborhood,const char * name,SCIP_Real minfixingrate,SCIP_Real maxfixingrate,SCIP_Bool active,SCIP_Real priority,DECL_VARFIXINGS ((* varfixings)),DECL_CHANGESUBSCIP ((* changesubscip)),DECL_NHINIT ((* nhinit)),DECL_NHEXIT ((* nhexit)),DECL_NHFREE ((* nhfree)),DECL_NHREFSOL ((* nhrefsol)),DECL_NHDEACTIVATE ((* nhdeactivate)))766 SCIP_RETCODE alnsIncludeNeighborhood(
767 SCIP* scip, /**< SCIP data structure */
768 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS heuristic */
769 NH** neighborhood, /**< pointer to store the neighborhood */
770 const char* name, /**< name for this neighborhood */
771 SCIP_Real minfixingrate, /**< default value for minfixingrate parameter of this neighborhood */
772 SCIP_Real maxfixingrate, /**< default value for maxfixingrate parameter of this neighborhood */
773 SCIP_Bool active, /**< default value for active parameter of this neighborhood */
774 SCIP_Real priority, /**< positive call priority to initialize bandit algorithms */
775 DECL_VARFIXINGS ((*varfixings)), /**< variable fixing callback for this neighborhood, or NULL */
776 DECL_CHANGESUBSCIP ((*changesubscip)), /**< subscip changes callback for this neighborhood, or NULL */
777 DECL_NHINIT ((*nhinit)), /**< initialization callback for neighborhood, or NULL */
778 DECL_NHEXIT ((*nhexit)), /**< deinitialization callback for neighborhood, or NULL */
779 DECL_NHFREE ((*nhfree)), /**< deinitialization callback before SCIP is freed, or NULL */
780 DECL_NHREFSOL ((*nhrefsol)), /**< callback function to return a reference solution for further fixings, or NULL */
781 DECL_NHDEACTIVATE ((*nhdeactivate)) /**< callback function to deactivate neighborhoods on problems where they are irrelevant, or NULL if neighborhood is always active */
782 )
783 {
784 char paramname[SCIP_MAXSTRLEN];
785
786 assert(scip != NULL);
787 assert(heurdata != NULL);
788 assert(neighborhood != NULL);
789 assert(name != NULL);
790
791 SCIP_CALL( SCIPallocBlockMemory(scip, neighborhood) );
792 assert(*neighborhood != NULL);
793
794 SCIP_ALLOC( BMSduplicateMemoryArray(&(*neighborhood)->name, name, strlen(name)+1) );
795
796 SCIP_CALL( SCIPcreateClock(scip, &(*neighborhood)->stats.setupclock) );
797 SCIP_CALL( SCIPcreateClock(scip, &(*neighborhood)->stats.submipclock) );
798
799 (*neighborhood)->changesubscip = changesubscip;
800 (*neighborhood)->varfixings = varfixings;
801 (*neighborhood)->nhinit = nhinit;
802 (*neighborhood)->nhexit = nhexit;
803 (*neighborhood)->nhfree = nhfree;
804 (*neighborhood)->nhrefsol = nhrefsol;
805 (*neighborhood)->nhdeactivate = nhdeactivate;
806
807 /* add parameters for this neighborhood */
808 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "heuristics/alns/%s/minfixingrate", name);
809 SCIP_CALL( SCIPaddRealParam(scip, paramname, "minimum fixing rate for this neighborhood",
810 &(*neighborhood)->fixingrate.minfixingrate, TRUE, minfixingrate, 0.0, 1.0, NULL, NULL) );
811 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "heuristics/alns/%s/maxfixingrate", name);
812 SCIP_CALL( SCIPaddRealParam(scip, paramname, "maximum fixing rate for this neighborhood",
813 &(*neighborhood)->fixingrate.maxfixingrate, TRUE, maxfixingrate, 0.0, 1.0, NULL, NULL) );
814 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "heuristics/alns/%s/active", name);
815 SCIP_CALL( SCIPaddBoolParam(scip, paramname, "is this neighborhood active?",
816 &(*neighborhood)->active, TRUE, active, NULL, NULL) );
817 (void) SCIPsnprintf(paramname, SCIP_MAXSTRLEN, "heuristics/alns/%s/priority", name);
818 SCIP_CALL( SCIPaddRealParam(scip, paramname, "positive call priority to initialize bandit algorithms",
819 &(*neighborhood)->priority, TRUE, priority, 1e-2, 1.0, NULL, NULL) );
820
821 /* add the neighborhood to the ALNS heuristic */
822 heurdata->neighborhoods[heurdata->nneighborhoods++] = (*neighborhood);
823
824 return SCIP_OKAY;
825 }
826
827 /** release all data and free neighborhood */
828 static
alnsFreeNeighborhood(SCIP * scip,NH ** neighborhood)829 SCIP_RETCODE alnsFreeNeighborhood(
830 SCIP* scip, /**< SCIP data structure */
831 NH** neighborhood /**< pointer to neighborhood that should be freed */
832 )
833 {
834 NH* nhptr;
835 assert(scip != NULL);
836 assert(neighborhood != NULL);
837
838 nhptr = *neighborhood;
839 assert(nhptr != NULL);
840
841 BMSfreeMemoryArray(&nhptr->name);
842
843 /* release further, neighborhood specific data structures */
844 if( nhptr->nhfree != NULL )
845 {
846 SCIP_CALL( nhptr->nhfree(scip, nhptr) );
847 }
848
849 SCIP_CALL( SCIPfreeClock(scip, &nhptr->stats.setupclock) );
850 SCIP_CALL( SCIPfreeClock(scip, &nhptr->stats.submipclock) );
851
852 SCIPfreeBlockMemory(scip, neighborhood);
853 *neighborhood = NULL;
854
855 return SCIP_OKAY;
856 }
857
858 /** initialize neighborhood specific data */
859 static
neighborhoodInit(SCIP * scip,NH * neighborhood)860 SCIP_RETCODE neighborhoodInit(
861 SCIP* scip, /**< SCIP data structure */
862 NH* neighborhood /**< neighborhood to initialize */
863 )
864 {
865 assert(scip != NULL);
866 assert(neighborhood != NULL);
867
868 /* call the init callback of the neighborhood */
869 if( neighborhood->nhinit != NULL )
870 {
871 SCIP_CALL( neighborhood->nhinit(scip, neighborhood) );
872 }
873
874 return SCIP_OKAY;
875 }
876
877 /** deinitialize neighborhood specific data */
878 static
neighborhoodExit(SCIP * scip,NH * neighborhood)879 SCIP_RETCODE neighborhoodExit(
880 SCIP* scip, /**< SCIP data structure */
881 NH* neighborhood /**< neighborhood to initialize */
882 )
883 {
884 assert(scip != NULL);
885 assert(neighborhood != NULL);
886
887 if( neighborhood->nhexit != NULL )
888 {
889 SCIP_CALL( neighborhood->nhexit(scip, neighborhood) );
890 }
891
892 return SCIP_OKAY;
893 }
894
895 /** creates a new solution for the original problem by copying the solution of the subproblem */
896 static
transferSolution(SCIP * subscip,SCIP_EVENTDATA * eventdata)897 SCIP_RETCODE transferSolution(
898 SCIP* subscip, /**< SCIP data structure of the subproblem */
899 SCIP_EVENTDATA* eventdata /**< event handler data */
900 )
901 {
902 SCIP* sourcescip; /* original SCIP data structure */
903 SCIP_VAR** subvars; /* the variables of the subproblem */
904 SCIP_HEUR* heur; /* alns heuristic structure */
905 SCIP_SOL* subsol; /* solution of the subproblem */
906 SCIP_SOL* newsol; /* solution to be created for the original problem */
907 SCIP_Bool success;
908 NH_STATS* runstats;
909 SCIP_SOL* oldbestsol;
910
911 assert(subscip != NULL);
912
913 subsol = SCIPgetBestSol(subscip);
914 assert(subsol != NULL);
915
916 sourcescip = eventdata->sourcescip;
917 subvars = eventdata->subvars;
918 heur = eventdata->heur;
919 runstats = eventdata->runstats;
920 assert(sourcescip != NULL);
921 assert(sourcescip != subscip);
922 assert(heur != NULL);
923 assert(subvars != NULL);
924 assert(runstats != NULL);
925
926 SCIP_CALL( SCIPtranslateSubSol(sourcescip, subscip, subsol, heur, subvars, &newsol) );
927
928 oldbestsol = SCIPgetBestSol(sourcescip);
929
930 /* in the special, experimental all rewards mode, the solution is only checked for feasibility
931 * but not stored
932 */
933 if( eventdata->allrewardsmode )
934 {
935 SCIP_CALL( SCIPcheckSol(sourcescip, newsol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) );
936
937 if( success )
938 {
939 runstats->nsolsfound++;
940 if( SCIPgetSolTransObj(sourcescip, newsol) < SCIPgetCutoffbound(sourcescip) )
941 runstats->nbestsolsfound++;
942 }
943
944 SCIP_CALL( SCIPfreeSol(sourcescip, &newsol) );
945 }
946 else
947 {
948 /* try to add new solution to scip and free it immediately */
949 SCIP_CALL( SCIPtrySolFree(sourcescip, &newsol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) );
950
951 if( success )
952 {
953 runstats->nsolsfound++;
954 if( SCIPgetBestSol(sourcescip) != oldbestsol )
955 runstats->nbestsolsfound++;
956 }
957 }
958
959 /* update new upper bound for reward later */
960 runstats->newupperbound = SCIPgetUpperbound(sourcescip);
961
962 return SCIP_OKAY;
963 }
964
965
966 /* ---------------- Callback methods of event handler ---------------- */
967
968 /** execution callback of the event handler
969 *
970 * transfer new solutions or interrupt the solving process manually
971 */
972 static
SCIP_DECL_EVENTEXEC(eventExecAlns)973 SCIP_DECL_EVENTEXEC(eventExecAlns)
974 {
975 assert(eventhdlr != NULL);
976 assert(eventdata != NULL);
977 assert(strcmp(SCIPeventhdlrGetName(eventhdlr), EVENTHDLR_NAME) == 0);
978 assert(event != NULL);
979 assert(SCIPeventGetType(event) & SCIP_EVENTTYPE_ALNS);
980 assert(eventdata != NULL);
981
982 /* treat the different atomic events */
983 switch( SCIPeventGetType(event) )
984 {
985 case SCIP_EVENTTYPE_SOLFOUND:
986 case SCIP_EVENTTYPE_BESTSOLFOUND:
987 /* try to transfer the solution to the original SCIP */
988 SCIP_CALL( transferSolution(scip, eventdata) );
989 break;
990 case SCIP_EVENTTYPE_LPSOLVED:
991 /* interrupt solution process of sub-SCIP */
992 if( SCIPgetNLPs(scip) > eventdata->lplimfac * eventdata->nodelimit )
993 {
994 SCIPdebugMsg(scip, "interrupt after %" SCIP_LONGINT_FORMAT " LPs\n", SCIPgetNLPs(scip));
995 SCIP_CALL( SCIPinterruptSolve(scip) );
996 }
997 break;
998 default:
999 break;
1000 }
1001
1002 return SCIP_OKAY;
1003 }
1004
1005 /** initialize neighborhood statistics before the next run */
1006 static
initRunStats(SCIP * scip,NH_STATS * stats)1007 void initRunStats(
1008 SCIP* scip, /**< SCIP data structure */
1009 NH_STATS* stats /**< run statistics */
1010 )
1011 {
1012 stats->nbestsolsfound = 0;
1013 stats->nsolsfound = 0;
1014 stats->usednodes = 0L;
1015 stats->nfixings = 0;
1016 stats->oldupperbound = SCIPgetUpperbound(scip);
1017 stats->newupperbound = SCIPgetUpperbound(scip);
1018 }
1019
1020 /** update run stats after the sub SCIP was solved */
1021 static
updateRunStats(NH_STATS * stats,SCIP * subscip)1022 void updateRunStats(
1023 NH_STATS* stats, /**< run statistics */
1024 SCIP* subscip /**< sub-SCIP instance, or NULL */
1025 )
1026 {
1027 /* treat an untransformed subscip as if none was created */
1028 if( subscip != NULL && ! SCIPisTransformed(subscip) )
1029 subscip = NULL;
1030
1031 stats->usednodes = subscip != NULL ? SCIPgetNNodes(subscip) : 0L;
1032 }
1033
1034 /** get the histogram index for this status */
1035 static
getHistIndex(SCIP_STATUS subscipstatus)1036 int getHistIndex(
1037 SCIP_STATUS subscipstatus /**< sub-SCIP status */
1038 )
1039 {
1040 switch (subscipstatus)
1041 {
1042 case SCIP_STATUS_OPTIMAL:
1043 return (int)HIDX_OPT;
1044 case SCIP_STATUS_INFEASIBLE:
1045 return (int)HIDX_INFEAS;
1046 case SCIP_STATUS_NODELIMIT:
1047 return (int)HIDX_NODELIM;
1048 case SCIP_STATUS_STALLNODELIMIT:
1049 return (int)HIDX_STALLNODE;
1050 case SCIP_STATUS_SOLLIMIT:
1051 case SCIP_STATUS_BESTSOLLIMIT:
1052 return (int)HIDX_SOLLIM;
1053 case SCIP_STATUS_USERINTERRUPT:
1054 return (int)HIDX_USR;
1055 default:
1056 return (int)HIDX_OTHER;
1057 } /*lint !e788*/
1058 }
1059
1060 /** print neighborhood statistics */
1061 static
printNeighborhoodStatistics(SCIP * scip,SCIP_HEURDATA * heurdata,FILE * file)1062 void printNeighborhoodStatistics(
1063 SCIP* scip, /**< SCIP data structure */
1064 SCIP_HEURDATA* heurdata, /**< heuristic data */
1065 FILE* file /**< file handle, or NULL for standard out */
1066 )
1067 {
1068 int i;
1069 int j;
1070 HISTINDEX statusses[] = {HIDX_OPT, HIDX_INFEAS, HIDX_NODELIM, HIDX_STALLNODE, HIDX_SOLLIM, HIDX_USR, HIDX_OTHER};
1071
1072 SCIPinfoMessage(scip, file, "Neighborhoods : %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %4s %4s %4s %4s %4s %4s %4s %4s\n",
1073 "Calls", "SetupTime", "SolveTime", "SolveNodes", "Sols", "Best", "Exp3", "EpsGreedy", "UCB", "TgtFixRate",
1074 "Opt", "Inf", "Node", "Stal", "Sol", "Usr", "Othr", "Actv");
1075
1076 /* loop over neighborhoods and fill in statistics */
1077 for( i = 0; i < heurdata->nneighborhoods; ++i )
1078 {
1079 NH* neighborhood;
1080 SCIP_Real proba;
1081 SCIP_Real ucb;
1082 SCIP_Real epsgreedyweight;
1083
1084 neighborhood = heurdata->neighborhoods[i];
1085 SCIPinfoMessage(scip, file, " %-17s:", neighborhood->name);
1086 SCIPinfoMessage(scip, file, " %10d", neighborhood->stats.nruns);
1087 SCIPinfoMessage(scip, file, " %10.2f", SCIPgetClockTime(scip, neighborhood->stats.setupclock) );
1088 SCIPinfoMessage(scip, file, " %10.2f", SCIPgetClockTime(scip, neighborhood->stats.submipclock) );
1089 SCIPinfoMessage(scip, file, " %10" SCIP_LONGINT_FORMAT, neighborhood->stats.usednodes );
1090 SCIPinfoMessage(scip, file, " %10" SCIP_LONGINT_FORMAT, neighborhood->stats.nsolsfound);
1091 SCIPinfoMessage(scip, file, " %10" SCIP_LONGINT_FORMAT, neighborhood->stats.nbestsolsfound);
1092
1093 proba = 0.0;
1094 ucb = 1.0;
1095 epsgreedyweight = -1.0;
1096
1097 if( heurdata->bandit != NULL && i < heurdata->nactiveneighborhoods )
1098 {
1099 switch (heurdata->banditalgo)
1100 {
1101 case 'u':
1102 ucb = SCIPgetConfidenceBoundUcb(heurdata->bandit, i);
1103 break;
1104 case 'g':
1105 epsgreedyweight = SCIPgetWeightsEpsgreedy(heurdata->bandit)[i];
1106 break;
1107 case 'e':
1108 proba = SCIPgetProbabilityExp3(heurdata->bandit, i);
1109 break;
1110 default:
1111 break;
1112 }
1113 }
1114
1115 SCIPinfoMessage(scip, file, " %10.5f", proba);
1116 SCIPinfoMessage(scip, file, " %10.5f", epsgreedyweight);
1117 SCIPinfoMessage(scip, file, " %10.5f", ucb);
1118 SCIPinfoMessage(scip, file, " %10.3f", neighborhood->fixingrate.targetfixingrate);
1119
1120 /* loop over status histogram */
1121 for( j = 0; j < NHISTENTRIES; ++j )
1122 SCIPinfoMessage(scip, file, " %4d", neighborhood->stats.statushist[statusses[j]]);
1123
1124 SCIPinfoMessage(scip, file, " %4d", i < heurdata->nactiveneighborhoods ? 1 : 0);
1125 SCIPinfoMessage(scip, file, "\n");
1126 }
1127 }
1128
1129 /** update the statistics of the neighborhood based on the sub-SCIP run */
1130 static
updateNeighborhoodStats(NH_STATS * runstats,NH * neighborhood,SCIP_STATUS subscipstatus)1131 void updateNeighborhoodStats(
1132 NH_STATS* runstats, /**< run statistics */
1133 NH* neighborhood, /**< the selected neighborhood */
1134 SCIP_STATUS subscipstatus /**< status of the sub-SCIP solve */
1135 )
1136 { /*lint --e{715}*/
1137 NH_STATS* stats;
1138 stats = &neighborhood->stats;
1139
1140 /* copy run statistics into neighborhood statistics */
1141 stats->nbestsolsfound += runstats->nbestsolsfound;
1142 stats->nsolsfound += runstats->nsolsfound;
1143 stats->usednodes += runstats->usednodes;
1144 stats->nruns += 1;
1145
1146 if( runstats->nbestsolsfound > 0 )
1147 stats->nrunsbestsol += DEFAULT_BESTSOLWEIGHT;
1148 else if( runstats->nsolsfound > 0 )
1149 stats->nrunsbestsol++;
1150
1151 /* update the counter for the subscip status */
1152 ++stats->statushist[getHistIndex(subscipstatus)];
1153 }
1154
1155 /** sort callback for variable pointers using the ALNS variable prioritization
1156 *
1157 * the variable prioritization works hierarchically as follows. A variable
1158 * a has the higher priority over b iff
1159 *
1160 * - variable distances should be used and a has a smaller distance than b
1161 * - variable reduced costs should be used and a has a smaller score than b
1162 * - variable pseudo costs should be used and a has a smaller score than b
1163 * - based on previously assigned random scores
1164 *
1165 * @note: distances are context-based. For fixing more variables,
1166 * distances are initialized from the already fixed variables.
1167 * For unfixing variables, distances are initialized starting
1168 * from the unfixed variables
1169 */
1170 static
SCIP_DECL_SORTINDCOMP(sortIndCompAlns)1171 SCIP_DECL_SORTINDCOMP(sortIndCompAlns)
1172 { /*lint --e{715}*/
1173 VARPRIO* varprio;
1174
1175 varprio = (VARPRIO*)dataptr;
1176 assert(varprio != NULL);
1177 assert(varprio->randscores != NULL);
1178
1179 if( ind1 == ind2 )
1180 return 0;
1181
1182 /* priority is on distances, if enabled. The variable which is closer in a breadth-first search sense to
1183 * the already fixed variables has precedence */
1184 if( varprio->usedistances )
1185 {
1186 int dist1;
1187 int dist2;
1188
1189 dist1 = varprio->distances[ind1];
1190 dist2 = varprio->distances[ind2];
1191
1192 if( dist1 < 0 )
1193 dist1 = INT_MAX;
1194
1195 if( dist2 < 0 )
1196 dist2 = INT_MAX;
1197
1198 assert(varprio->distances != NULL);
1199 if( dist1 < dist2 )
1200 return -1;
1201 else if( dist1 > dist2 )
1202 return 1;
1203 }
1204
1205 assert(! varprio->usedistances || varprio->distances[ind1] == varprio->distances[ind2]);
1206
1207 /* if the indices tie considering distances or distances are disabled -> use reduced cost information instead */
1208 if( varprio->useredcost )
1209 {
1210 assert(varprio->redcostscores != NULL);
1211
1212 if( varprio->redcostscores[ind1] < varprio->redcostscores[ind2] )
1213 return -1;
1214 else if( varprio->redcostscores[ind1] > varprio->redcostscores[ind2] )
1215 return 1;
1216 }
1217
1218 /* use pseudo cost scores if reduced costs are disabled or a tie was found */
1219 if( varprio->usepscost )
1220 {
1221 assert(varprio->pscostscores != NULL);
1222
1223 /* prefer the variable with smaller pseudocost score */
1224 if( varprio->pscostscores[ind1] < varprio->pscostscores[ind2] )
1225 return -1;
1226 else if( varprio->pscostscores[ind1] > varprio->pscostscores[ind2] )
1227 return 1;
1228 }
1229
1230 if( varprio->randscores[ind1] < varprio->randscores[ind2] )
1231 return -1;
1232 else if( varprio->randscores[ind1] > varprio->randscores[ind2] )
1233 return 1;
1234
1235 return ind1 - ind2;
1236 }
1237
1238 /** Compute the reduced cost score for this variable in the reference solution */
1239 static
getVariableRedcostScore(SCIP * scip,SCIP_VAR * var,SCIP_Real refsolval,SCIP_Bool uselocalredcost)1240 SCIP_Real getVariableRedcostScore(
1241 SCIP* scip, /**< SCIP data structure */
1242 SCIP_VAR* var, /**< the variable for which the score should be computed */
1243 SCIP_Real refsolval, /**< solution value in reference solution */
1244 SCIP_Bool uselocalredcost /**< should local reduced costs be used for generic (un)fixing? */
1245 )
1246 {
1247 SCIP_Real bestbound;
1248 SCIP_Real redcost;
1249 SCIP_Real score;
1250 assert(scip != NULL);
1251 assert(var != NULL);
1252
1253 /* prefer column variables */
1254 if( SCIPvarGetStatus(var) != SCIP_VARSTATUS_COLUMN )
1255 return SCIPinfinity(scip);
1256
1257 if( ! uselocalredcost )
1258 {
1259 redcost = SCIPvarGetBestRootRedcost(var);
1260
1261 bestbound = SCIPvarGetBestRootSol(var);
1262
1263 /* using global reduced costs, the two factors yield a nonnegative score within tolerances */
1264 assert(SCIPisDualfeasZero(scip, redcost)
1265 || (SCIPisDualfeasNegative(scip, redcost) && ! SCIPisFeasPositive(scip, refsolval - bestbound))
1266 || (SCIPisDualfeasPositive(scip, redcost) && ! SCIPisFeasNegative(scip, refsolval - bestbound)));
1267 }
1268 else
1269 {
1270 /* this can be safely asserted here, since the heuristic would not reach this point, otherwise */
1271 assert(SCIPhasCurrentNodeLP(scip));
1272 assert(SCIPgetLPSolstat(scip) == SCIP_LPSOLSTAT_OPTIMAL);
1273
1274 redcost = SCIPgetVarRedcost(scip, var);
1275
1276 bestbound = SCIPvarGetLPSol(var);
1277 }
1278
1279 assert(! SCIPisInfinity(scip, REALABS(bestbound)));
1280 assert(SCIPisDualfeasZero(scip, redcost) || SCIPisFeasIntegral(scip, bestbound));
1281
1282 score = redcost * (refsolval - bestbound);
1283
1284 /* max out numerical inaccuracies from global scores */
1285 if( ! uselocalredcost )
1286 score = MAX(score, 0.0);
1287
1288 return score;
1289 }
1290
1291 /** get the pseudo cost score of this variable with respect to the reference solution */
1292 static
getVariablePscostScore(SCIP * scip,SCIP_VAR * var,SCIP_Real refsolval)1293 SCIP_Real getVariablePscostScore(
1294 SCIP* scip, /**< SCIP data structure */
1295 SCIP_VAR* var, /**< the variable for which the score should be computed */
1296 SCIP_Real refsolval /**< solution value in reference solution */
1297 )
1298 {
1299 SCIP_Real rootsolval;
1300
1301 assert(scip != NULL);
1302 assert(var != NULL);
1303
1304 /* variables that aren't LP columns have no pseudocost score */
1305 if( SCIPvarGetStatus(var) != SCIP_VARSTATUS_COLUMN )
1306 return 0.0;
1307
1308 rootsolval = SCIPvarGetRootSol(var);
1309
1310 /* the score is 0.0 if the values are equal */
1311 if( SCIPisEQ(scip, rootsolval, refsolval) )
1312 return 0.0;
1313 else
1314 return SCIPgetVarPseudocostVal(scip, var, refsolval - rootsolval);
1315 }
1316
1317 /** add variable and solution value to buffer data structure for variable fixings. The method checks if
1318 * the value still lies within the variable bounds. The value stays unfixed otherwise.
1319 */
1320 static
tryAdd2variableBuffer(SCIP * scip,SCIP_VAR * var,SCIP_Real val,SCIP_VAR ** varbuf,SCIP_Real * valbuf,int * nfixings,SCIP_Bool integer)1321 void tryAdd2variableBuffer(
1322 SCIP* scip, /**< SCIP data structure */
1323 SCIP_VAR* var, /**< (source) SCIP variable that should be added to the buffer */
1324 SCIP_Real val, /**< fixing value for this variable */
1325 SCIP_VAR** varbuf, /**< variable buffer to store variables that should be fixed */
1326 SCIP_Real* valbuf, /**< value buffer to store fixing values */
1327 int* nfixings, /**< pointer to number of fixed buffer variables, will be increased by 1 */
1328 SCIP_Bool integer /**< is this an integer variable? */
1329 )
1330 {
1331 assert(SCIPisFeasIntegral(scip, val) || ! SCIPvarIsIntegral(var));
1332 assert(*nfixings < SCIPgetNVars(scip));
1333
1334 /* round the value to its nearest integer */
1335 if( integer )
1336 val = SCIPfloor(scip, val + 0.5);
1337
1338 /* only add fixing if it is still valid within the global variable bounds. Invalidity
1339 * of this solution value may come from a dual reduction that was performed after the solution from which
1340 * this value originated was found
1341 */
1342 if( SCIPvarGetLbGlobal(var) <= val && val <= SCIPvarGetUbGlobal(var) )
1343 {
1344 varbuf[*nfixings] = var;
1345 valbuf[*nfixings] = val;
1346 ++(*nfixings);
1347 }
1348 }
1349
1350 /** query neighborhood for a reference solution for further fixings */
1351 static
neighborhoodGetRefsol(SCIP * scip,NH * neighborhood,SCIP_SOL ** solptr)1352 SCIP_RETCODE neighborhoodGetRefsol(
1353 SCIP* scip, /**< SCIP data structure */
1354 NH* neighborhood, /**< ALNS neighborhood data structure */
1355 SCIP_SOL** solptr /**< solution pointer */
1356 )
1357 {
1358 assert(solptr != NULL);
1359 assert(scip != NULL);
1360 assert(neighborhood != NULL);
1361
1362 *solptr = NULL;
1363 if( neighborhood->nhrefsol != NULL )
1364 {
1365 SCIP_RESULT result;
1366 SCIP_CALL( neighborhood->nhrefsol(scip, neighborhood, solptr, &result) );
1367
1368 if( result == SCIP_DIDNOTFIND )
1369 *solptr = NULL;
1370 else
1371 assert(*solptr != NULL);
1372 }
1373
1374 return SCIP_OKAY;
1375 }
1376
1377 /** fix additional variables found in feasible reference solution if the ones that the neighborhood found were not enough
1378 *
1379 * use not always the best solution for the values, but a reference solution provided by the neighborhood itself
1380 *
1381 * @note it may happen that the target fixing rate is not completely reached. This is the case if intermediate,
1382 * dual reductions render the solution values of the reference solution infeasible for
1383 * the current, global variable bounds.
1384 */
1385 static
alnsFixMoreVariables(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_SOL * refsol,SCIP_VAR ** varbuf,SCIP_Real * valbuf,int * nfixings,int ntargetfixings,SCIP_Bool * success)1386 SCIP_RETCODE alnsFixMoreVariables(
1387 SCIP* scip, /**< SCIP data structure */
1388 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
1389 SCIP_SOL* refsol, /**< feasible reference solution for more variable fixings */
1390 SCIP_VAR** varbuf, /**< buffer array to store variables to fix */
1391 SCIP_Real* valbuf, /**< buffer array to store fixing values */
1392 int* nfixings, /**< pointer to store the number of fixings */
1393 int ntargetfixings, /**< number of required target fixings */
1394 SCIP_Bool* success /**< pointer to store whether the target fixings have been successfully reached */
1395 )
1396 {
1397 VARPRIO varprio;
1398 SCIP_VAR** vars;
1399 SCIP_Real* redcostscores;
1400 SCIP_Real* pscostscores;
1401 SCIP_Real* solvals;
1402 SCIP_RANDNUMGEN* rng;
1403 SCIP_VAR** unfixedvars;
1404 SCIP_Bool* isfixed;
1405 int* distances;
1406 int* perm;
1407 SCIP_Real* randscores;
1408 int nbinvars;
1409 int nintvars;
1410 int nbinintvars;
1411 int nvars;
1412 int b;
1413 int nvarstoadd;
1414 int nunfixedvars;
1415
1416 assert(scip != NULL);
1417 assert(varbuf != NULL);
1418 assert(nfixings != NULL);
1419 assert(success != NULL);
1420 assert(heurdata != NULL);
1421 assert(refsol != NULL);
1422
1423 *success = FALSE;
1424
1425 /* if the user parameter forbids more fixings, return immediately */
1426 if( ! heurdata->domorefixings )
1427 return SCIP_OKAY;
1428
1429 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) );
1430
1431 nbinintvars = nbinvars + nintvars;
1432
1433 if( ntargetfixings >= nbinintvars )
1434 return SCIP_OKAY;
1435
1436 /* determine the number of required additional fixings */
1437 nvarstoadd = ntargetfixings - *nfixings;
1438 if( nvarstoadd == 0 )
1439 return SCIP_OKAY;
1440
1441 varprio.usedistances = heurdata->usedistances && (*nfixings >= 1);
1442 varprio.useredcost = heurdata->useredcost;
1443 varprio.usepscost = heurdata->usepscost;
1444 varprio.scip = scip;
1445 rng = SCIPbanditGetRandnumgen(heurdata->bandit);
1446 assert(rng != NULL);
1447
1448 SCIP_CALL( SCIPallocBufferArray(scip, &randscores, nbinintvars) );
1449 SCIP_CALL( SCIPallocBufferArray(scip, &perm, nbinintvars) );
1450 SCIP_CALL( SCIPallocBufferArray(scip, &distances, nvars) );
1451 SCIP_CALL( SCIPallocBufferArray(scip, &redcostscores, nbinintvars) );
1452 SCIP_CALL( SCIPallocBufferArray(scip, &solvals, nbinintvars) );
1453 SCIP_CALL( SCIPallocBufferArray(scip, &isfixed, nbinintvars) );
1454 SCIP_CALL( SCIPallocBufferArray(scip, &unfixedvars, nbinintvars) );
1455 SCIP_CALL( SCIPallocBufferArray(scip, &pscostscores, nbinintvars) );
1456
1457 /* initialize variable graph distances from already fixed variables */
1458 if( varprio.usedistances )
1459 {
1460 SCIP_CALL( SCIPvariablegraphBreadthFirst(scip, NULL, varbuf, *nfixings, distances, INT_MAX, INT_MAX, ntargetfixings) );
1461 }
1462 else
1463 {
1464 /* initialize all equal distances to make them irrelevant */
1465 BMSclearMemoryArray(distances, nbinintvars);
1466 }
1467
1468 BMSclearMemoryArray(isfixed, nbinintvars);
1469
1470 /* mark binary and integer variables if they are fixed */
1471 for( b = 0; b < *nfixings; ++b )
1472 {
1473 int probindex;
1474
1475 assert(varbuf[b] != NULL);
1476 probindex = SCIPvarGetProbindex(varbuf[b]);
1477 assert(probindex >= 0);
1478
1479 if( probindex < nbinintvars )
1480 isfixed[probindex] = TRUE;
1481 }
1482
1483 SCIP_CALL( SCIPgetSolVals(scip, refsol, nbinintvars, vars, solvals) );
1484
1485 /* assign scores to unfixed every discrete variable of the problem */
1486 nunfixedvars = 0;
1487 for( b = 0; b < nbinintvars; ++b )
1488 {
1489 SCIP_VAR* var = vars[b];
1490
1491 /* filter fixed variables */
1492 if( isfixed[b] )
1493 continue;
1494
1495 /* filter variables with a solution value outside its global bounds */
1496 if( solvals[b] < SCIPvarGetLbGlobal(var) - 0.5 || solvals[b] > SCIPvarGetUbGlobal(var) + 0.5 )
1497 continue;
1498
1499 redcostscores[nunfixedvars] = getVariableRedcostScore(scip, var, solvals[b], heurdata->uselocalredcost);
1500 pscostscores[nunfixedvars] = getVariablePscostScore(scip, var, solvals[b]);
1501
1502 unfixedvars[nunfixedvars] = var;
1503 perm[nunfixedvars] = nunfixedvars;
1504 randscores[nunfixedvars] = SCIPrandomGetReal(rng, 0.0, 1.0);
1505
1506 /* these assignments are based on the fact that nunfixedvars <= b */
1507 solvals[nunfixedvars] = solvals[b];
1508 distances[nunfixedvars] = distances[b];
1509
1510 SCIPdebugMsg(scip, "Var <%s> scores: dist %3d, red cost %15.9g, pscost %15.9g rand %6.4f\n",
1511 SCIPvarGetName(var), distances[nunfixedvars], redcostscores[nunfixedvars],
1512 pscostscores[nunfixedvars], randscores[nunfixedvars]);
1513
1514 nunfixedvars++;
1515 }
1516
1517 /* use selection algorithm (order of the variables does not matter) for quickly completing the fixing */
1518 varprio.randscores = randscores;
1519 varprio.distances = distances;
1520 varprio.redcostscores = redcostscores;
1521 varprio.pscostscores = pscostscores;
1522
1523 nvarstoadd = MIN(nunfixedvars, nvarstoadd);
1524
1525 /* select the first nvarstoadd many variables according to the score */
1526 if( nvarstoadd < nunfixedvars )
1527 SCIPselectInd(perm, sortIndCompAlns, &varprio, nvarstoadd, nunfixedvars);
1528
1529 /* loop over the first elements of the selection defined in permutation. They represent the best variables */
1530 for( b = 0; b < nvarstoadd; ++b )
1531 {
1532 int permindex = perm[b];
1533 assert(permindex >= 0);
1534 assert(permindex < nunfixedvars);
1535
1536 tryAdd2variableBuffer(scip, unfixedvars[permindex], solvals[permindex], varbuf, valbuf, nfixings, TRUE);
1537 }
1538
1539 *success = TRUE;
1540
1541 /* free buffer arrays */
1542 SCIPfreeBufferArray(scip, &pscostscores);
1543 SCIPfreeBufferArray(scip, &unfixedvars);
1544 SCIPfreeBufferArray(scip, &isfixed);
1545 SCIPfreeBufferArray(scip, &solvals);
1546 SCIPfreeBufferArray(scip, &redcostscores);
1547 SCIPfreeBufferArray(scip, &distances);
1548 SCIPfreeBufferArray(scip, &perm);
1549 SCIPfreeBufferArray(scip, &randscores);
1550
1551 return SCIP_OKAY;
1552 }
1553
1554 /** create the bandit algorithm for the heuristic depending on the user parameter */
1555 static
createBandit(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_Real * priorities,unsigned int initseed)1556 SCIP_RETCODE createBandit(
1557 SCIP* scip, /**< SCIP data structure */
1558 SCIP_HEURDATA* heurdata, /**< heuristic data structure */
1559 SCIP_Real* priorities, /**< call priorities for active neighborhoods */
1560 unsigned int initseed /**< initial random seed */
1561 )
1562 {
1563 switch (heurdata->banditalgo)
1564 {
1565 case 'u':
1566 SCIP_CALL( SCIPcreateBanditUcb(scip, &heurdata->bandit, priorities,
1567 heurdata->ucb_alpha, heurdata->nactiveneighborhoods, initseed) );
1568 break;
1569
1570 case 'e':
1571 SCIP_CALL( SCIPcreateBanditExp3(scip, &heurdata->bandit, priorities,
1572 heurdata->exp3_gamma, heurdata->exp3_beta, heurdata->nactiveneighborhoods, initseed) );
1573 break;
1574
1575 case 'g':
1576 SCIP_CALL( SCIPcreateBanditEpsgreedy(scip, &heurdata->bandit, priorities,
1577 heurdata->epsgreedy_eps, FALSE, 0.9, 0, heurdata->nactiveneighborhoods, initseed) );
1578 break;
1579
1580 default:
1581 SCIPerrorMessage("Unknown bandit parameter %c\n", heurdata->banditalgo);
1582 return SCIP_INVALIDDATA;
1583 }
1584
1585 return SCIP_OKAY;
1586 }
1587
1588 /*
1589 * Callback methods of primal heuristic
1590 */
1591
1592 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
1593 static
SCIP_DECL_HEURCOPY(heurCopyAlns)1594 SCIP_DECL_HEURCOPY(heurCopyAlns)
1595 { /*lint --e{715}*/
1596 assert(scip != NULL);
1597 assert(heur != NULL);
1598 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1599
1600 /* call inclusion method of primal heuristic */
1601 SCIP_CALL( SCIPincludeHeurAlns(scip) );
1602
1603 return SCIP_OKAY;
1604 }
1605
1606 /** unfix some of the variables because there are too many fixed
1607 *
1608 * a variable is ideally unfixed if it is close to other unfixed variables
1609 * and fixing it has a high reduced cost impact
1610 */
1611 static
alnsUnfixVariables(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_VAR ** varbuf,SCIP_Real * valbuf,int * nfixings,int ntargetfixings,SCIP_Bool * success)1612 SCIP_RETCODE alnsUnfixVariables(
1613 SCIP* scip, /**< SCIP data structure */
1614 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
1615 SCIP_VAR** varbuf, /**< buffer array to store variables to fix */
1616 SCIP_Real* valbuf, /**< buffer array to store fixing values */
1617 int* nfixings, /**< pointer to store the number of fixings */
1618 int ntargetfixings, /**< number of required target fixings */
1619 SCIP_Bool* success /**< pointer to store whether the target fixings have been successfully reached */
1620 )
1621 {
1622 VARPRIO varprio;
1623 SCIP_Real* redcostscores;
1624 SCIP_Real* pscostscores;
1625 SCIP_Real* randscores;
1626 SCIP_VAR** unfixedvars;
1627 SCIP_VAR** varbufcpy;
1628 SCIP_Real* valbufcpy;
1629 SCIP_Bool* isfixedvar;
1630 SCIP_VAR** vars;
1631 SCIP_RANDNUMGEN* rng;
1632 int* distances;
1633 int* fixeddistances;
1634 int* perm;
1635 int nvars;
1636 int i;
1637 int nbinintvars;
1638 int nunfixed;
1639
1640 *success = FALSE;
1641
1642 nbinintvars = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
1643 if( nbinintvars == 0 )
1644 return SCIP_OKAY;
1645
1646 assert(*nfixings > 0);
1647
1648 nvars = SCIPgetNVars(scip);
1649 SCIP_CALL( SCIPallocBufferArray(scip, &isfixedvar, nvars) );
1650 SCIP_CALL( SCIPallocBufferArray(scip, &unfixedvars, nbinintvars) );
1651 SCIP_CALL( SCIPallocBufferArray(scip, &distances, nvars) );
1652 SCIP_CALL( SCIPallocBufferArray(scip, &fixeddistances, *nfixings) );
1653 SCIP_CALL( SCIPallocBufferArray(scip, &redcostscores, *nfixings) );
1654 SCIP_CALL( SCIPallocBufferArray(scip, &randscores, *nfixings) );
1655 SCIP_CALL( SCIPallocBufferArray(scip, &perm, *nfixings) );
1656 SCIP_CALL( SCIPallocBufferArray(scip, &pscostscores, *nfixings) );
1657
1658 SCIP_CALL( SCIPduplicateBufferArray(scip, &varbufcpy, varbuf, *nfixings) );
1659 SCIP_CALL( SCIPduplicateBufferArray(scip, &valbufcpy, valbuf, *nfixings) );
1660
1661 /*
1662 * collect the unfixed binary and integer variables
1663 */
1664 BMSclearMemoryArray(isfixedvar, nvars);
1665 /* loop over fixed variables and mark their respective positions as fixed */
1666 for( i = 0; i < *nfixings; ++i )
1667 {
1668 int probindex = SCIPvarGetProbindex(varbuf[i]);
1669
1670 assert(probindex >= 0);
1671
1672 isfixedvar[probindex] = TRUE;
1673 }
1674
1675 nunfixed = 0;
1676 vars = SCIPgetVars(scip);
1677 /* collect unfixed binary and integer variables */
1678 for( i = 0; i < nbinintvars; ++i )
1679 {
1680 if( ! isfixedvar[i] )
1681 unfixedvars[nunfixed++] = vars[i];
1682 }
1683
1684 varprio.usedistances = heurdata->usedistances && nunfixed > 0;
1685
1686 /* collect distances of all fixed variables from those that are not fixed */
1687 if( varprio.usedistances )
1688 {
1689 SCIP_CALL( SCIPvariablegraphBreadthFirst(scip, NULL, unfixedvars, nunfixed, distances, INT_MAX, INT_MAX, INT_MAX) );
1690
1691 for( i = 0; i < *nfixings; ++i )
1692 {
1693 int probindex = SCIPvarGetProbindex(varbuf[i]);
1694 if( probindex >= 0 )
1695 fixeddistances[i] = distances[probindex];
1696 }
1697 }
1698 else
1699 {
1700 BMSclearMemoryArray(fixeddistances, *nfixings);
1701 }
1702
1703 /* collect reduced cost scores of the fixings and assign random scores */
1704 rng = SCIPbanditGetRandnumgen(heurdata->bandit);
1705 for( i = 0; i < *nfixings; ++i )
1706 {
1707 SCIP_VAR* fixedvar = varbuf[i];
1708 SCIP_Real fixval = valbuf[i];
1709
1710 /* use negative reduced cost and pseudo cost scores to prefer variable fixings with small score */
1711 redcostscores[i] = - getVariableRedcostScore(scip, fixedvar, fixval, heurdata->uselocalredcost);
1712 pscostscores[i] = - getVariablePscostScore(scip, fixedvar, fixval);
1713 randscores[i] = SCIPrandomGetReal(rng, 0.0, 1.0);
1714 perm[i] = i;
1715
1716 SCIPdebugMsg(scip, "Var <%s> scores: dist %3d, red cost %15.9g, pscost %15.9g rand %6.4f\n",
1717 SCIPvarGetName(fixedvar), fixeddistances[i], redcostscores[i], pscostscores[i], randscores[i]);
1718 }
1719
1720 varprio.distances = fixeddistances;
1721 varprio.randscores = randscores;
1722 varprio.redcostscores = redcostscores;
1723 varprio.pscostscores = pscostscores;
1724 varprio.useredcost = heurdata->useredcost;
1725 varprio.usepscost = heurdata->usepscost;
1726 varprio.scip = scip;
1727
1728 /* scores are assigned in such a way that variables with a smaller score should be fixed last */
1729 SCIPselectDownInd(perm, sortIndCompAlns, &varprio, ntargetfixings, *nfixings);
1730
1731 /* bring the desired variables to the front of the array */
1732 for( i = 0; i < ntargetfixings; ++i )
1733 {
1734 valbuf[i] = valbufcpy[perm[i]];
1735 varbuf[i] = varbufcpy[perm[i]];
1736 }
1737
1738 *nfixings = ntargetfixings;
1739
1740 /* free the buffer arrays in reverse order of allocation */
1741 SCIPfreeBufferArray(scip, &valbufcpy);
1742 SCIPfreeBufferArray(scip, &varbufcpy);
1743 SCIPfreeBufferArray(scip, &pscostscores);
1744 SCIPfreeBufferArray(scip, &perm);
1745 SCIPfreeBufferArray(scip, &randscores);
1746 SCIPfreeBufferArray(scip, &redcostscores);
1747 SCIPfreeBufferArray(scip, &fixeddistances);
1748 SCIPfreeBufferArray(scip, &distances);
1749 SCIPfreeBufferArray(scip, &unfixedvars);
1750 SCIPfreeBufferArray(scip, &isfixedvar);
1751
1752 *success = TRUE;
1753
1754 return SCIP_OKAY;
1755 }
1756
1757 /** call variable fixing callback for this neighborhood and orchestrate additional variable fixings, if necessary */
1758 static
neighborhoodFixVariables(SCIP * scip,SCIP_HEURDATA * heurdata,NH * neighborhood,SCIP_VAR ** varbuf,SCIP_Real * valbuf,int * nfixings,SCIP_RESULT * result)1759 SCIP_RETCODE neighborhoodFixVariables(
1760 SCIP* scip, /**< SCIP data structure */
1761 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
1762 NH* neighborhood, /**< neighborhood data structure */
1763 SCIP_VAR** varbuf, /**< buffer array to keep variables that should be fixed */
1764 SCIP_Real* valbuf, /**< buffer array to keep fixing values */
1765 int* nfixings, /**< pointer to store the number of variable fixings */
1766 SCIP_RESULT* result /**< pointer to store the result of the fixing operation */
1767 )
1768 {
1769 int ntargetfixings;
1770 int nmaxfixings;
1771 int nminfixings;
1772 int nbinintvars;
1773
1774 assert(scip != NULL);
1775 assert(neighborhood != NULL);
1776 assert(varbuf != NULL);
1777 assert(valbuf != NULL);
1778 assert(nfixings != NULL);
1779 assert(result != NULL);
1780
1781 *nfixings = 0;
1782
1783 *result = SCIP_DIDNOTRUN;
1784 ntargetfixings = (int)(neighborhood->fixingrate.targetfixingrate * (SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip)));
1785
1786 if( neighborhood->varfixings != NULL )
1787 {
1788 SCIP_CALL( neighborhood->varfixings(scip, neighborhood, varbuf, valbuf, nfixings, result) );
1789
1790 if( *result != SCIP_SUCCESS )
1791 return SCIP_OKAY;
1792 }
1793 else if( ntargetfixings == 0 )
1794 {
1795 *result = SCIP_SUCCESS;
1796
1797 return SCIP_OKAY;
1798 }
1799
1800 /* compute upper and lower target fixing limits using tolerance parameters */
1801 assert(neighborhood->varfixings == NULL || *result != SCIP_DIDNOTRUN);
1802 nbinintvars = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
1803 ntargetfixings = (int)(neighborhood->fixingrate.targetfixingrate * nbinintvars);
1804 nminfixings = (int)((neighborhood->fixingrate.targetfixingrate - heurdata->fixtol) * nbinintvars);
1805 nminfixings = MAX(nminfixings, 0);
1806 nmaxfixings = (int)((neighborhood->fixingrate.targetfixingrate + heurdata->unfixtol) * nbinintvars);
1807 nmaxfixings = MIN(nmaxfixings, nbinintvars);
1808
1809 SCIPdebugMsg(scip, "Neighborhood Fixings/Target: %d / %d <= %d <= %d\n",*nfixings, nminfixings, ntargetfixings, nmaxfixings);
1810
1811 /* if too few fixings, use a strategy to select more variable fixings: randomized, LP graph, ReducedCost based, mix */
1812 if( (*result == SCIP_SUCCESS || *result == SCIP_DIDNOTRUN) && (*nfixings < nminfixings) )
1813 {
1814 SCIP_Bool success;
1815 SCIP_SOL* refsol;
1816
1817 /* get reference solution from neighborhood */
1818 SCIP_CALL( neighborhoodGetRefsol(scip, neighborhood, &refsol) );
1819
1820 /* try to fix more variables based on the reference solution */
1821 if( refsol != NULL )
1822 {
1823 SCIP_CALL( alnsFixMoreVariables(scip, heurdata, refsol, varbuf, valbuf, nfixings, ntargetfixings, &success) );
1824 }
1825 else
1826 success = FALSE;
1827
1828 if( success )
1829 *result = SCIP_SUCCESS;
1830 else if( *result == SCIP_SUCCESS )
1831 *result = SCIP_DIDNOTFIND;
1832 else
1833 *result = SCIP_DIDNOTRUN;
1834
1835 SCIPdebugMsg(scip, "After additional fixings: %d / %d\n",*nfixings, ntargetfixings);
1836 }
1837 else if( (SCIP_Real)(*nfixings) > nmaxfixings )
1838 {
1839 SCIP_Bool success;
1840
1841 SCIP_CALL( alnsUnfixVariables(scip, heurdata, varbuf, valbuf, nfixings, ntargetfixings, &success) );
1842
1843 assert(success);
1844 *result = SCIP_SUCCESS;
1845 SCIPdebugMsg(scip, "Unfixed variables, fixed variables remaining: %d\n", ntargetfixings);
1846 }
1847 else
1848 {
1849 SCIPdebugMsg(scip, "No additional fixings performed\n");
1850 }
1851
1852 return SCIP_OKAY;
1853 }
1854
1855 /** change the sub-SCIP by restricting variable domains, changing objective coefficients, or adding constraints */
1856 static
neighborhoodChangeSubscip(SCIP * sourcescip,SCIP * targetscip,NH * neighborhood,SCIP_VAR ** targetvars,int * ndomchgs,int * nchgobjs,int * naddedconss,SCIP_Bool * success)1857 SCIP_RETCODE neighborhoodChangeSubscip(
1858 SCIP* sourcescip, /**< source SCIP data structure */
1859 SCIP* targetscip, /**< target SCIP data structure */
1860 NH* neighborhood, /**< neighborhood */
1861 SCIP_VAR** targetvars, /**< array of target SCIP variables aligned with source SCIP variables */
1862 int* ndomchgs, /**< pointer to store the number of variable domain changes */
1863 int* nchgobjs, /**< pointer to store the number of changed objective coefficients */
1864 int* naddedconss, /**< pointer to store the number of added constraints */
1865 SCIP_Bool* success /**< pointer to store whether the sub-SCIP has been successfully modified */
1866 )
1867 {
1868 assert(sourcescip != NULL);
1869 assert(targetscip != NULL);
1870 assert(neighborhood != NULL);
1871 assert(targetvars != NULL);
1872 assert(ndomchgs != NULL);
1873 assert(nchgobjs != NULL);
1874 assert(naddedconss != NULL);
1875 assert(success != NULL);
1876
1877 *success = FALSE;
1878 *ndomchgs = 0;
1879 *nchgobjs = 0;
1880 *naddedconss = 0;
1881
1882 /* call the change sub-SCIP callback of the neighborhood */
1883 if( neighborhood->changesubscip != NULL )
1884 {
1885 SCIP_CALL( neighborhood->changesubscip(sourcescip, targetscip, neighborhood, targetvars, ndomchgs, nchgobjs, naddedconss, success) );
1886 }
1887 else
1888 {
1889 *success = TRUE;
1890 }
1891
1892 return SCIP_OKAY;
1893 }
1894
1895 /** set sub-SCIP solving limits */
1896 static
setLimits(SCIP * subscip,SOLVELIMITS * solvelimits)1897 SCIP_RETCODE setLimits(
1898 SCIP* subscip, /**< SCIP data structure */
1899 SOLVELIMITS* solvelimits /**< pointer to solving limits data structure */
1900 )
1901 {
1902 assert(subscip != NULL);
1903 assert(solvelimits != NULL);
1904
1905 assert(solvelimits->nodelimit >= solvelimits->stallnodes);
1906
1907 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", solvelimits->nodelimit) );
1908 SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", solvelimits->stallnodes));
1909 SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", solvelimits->timelimit) );
1910 SCIP_CALL( SCIPsetRealParam(subscip, "limits/memory", solvelimits->memorylimit) );
1911
1912 return SCIP_OKAY;
1913 }
1914
1915 /** determine limits for a sub-SCIP */
1916 static
determineLimits(SCIP * scip,SCIP_HEUR * heur,SOLVELIMITS * solvelimits,SCIP_Bool * runagain)1917 SCIP_RETCODE determineLimits(
1918 SCIP* scip, /**< SCIP data structure */
1919 SCIP_HEUR* heur, /**< this heuristic */
1920 SOLVELIMITS* solvelimits, /**< pointer to solving limits data structure */
1921 SCIP_Bool* runagain /**< can we solve another sub-SCIP with these limits */
1922 )
1923 {
1924 SCIP_HEURDATA* heurdata;
1925 SCIP_Real initfactor;
1926 SCIP_Real nodesquot;
1927
1928 assert(scip != NULL);
1929 assert(heur != NULL);
1930 assert(solvelimits != NULL);
1931 assert(runagain != NULL);
1932
1933 heurdata = SCIPheurGetData(heur);
1934
1935 /* check whether there is enough time and memory left */
1936 SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &solvelimits->timelimit) );
1937 if( ! SCIPisInfinity(scip, solvelimits->timelimit) )
1938 solvelimits->timelimit -= SCIPgetSolvingTime(scip);
1939 SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &solvelimits->memorylimit) );
1940
1941 /* substract the memory already used by the main SCIP and the estimated memory usage of external software */
1942 if( ! SCIPisInfinity(scip, solvelimits->memorylimit) )
1943 {
1944 solvelimits->memorylimit -= SCIPgetMemUsed(scip)/1048576.0;
1945 solvelimits->memorylimit -= SCIPgetMemExternEstim(scip)/1048576.0;
1946 }
1947
1948 /* abort if no time is left or not enough memory to create a copy of SCIP, including external memory usage */
1949 if( solvelimits->timelimit <= 0.0 || solvelimits->memorylimit <= 2.0*SCIPgetMemExternEstim(scip)/1048576.0 )
1950 *runagain = FALSE;
1951
1952 nodesquot = heurdata->nodesquot;
1953
1954 /* if the heuristic is used to measure all rewards, it will always be penalized here */
1955 if( heurdata->rewardfile == NULL )
1956 nodesquot *= (SCIPheurGetNBestSolsFound(heur) + 1.0)/(SCIPheurGetNCalls(heur) + 1.0);
1957
1958 /* calculate the search node limit of the heuristic */
1959 solvelimits->stallnodes = (SCIP_Longint)(nodesquot * SCIPgetNNodes(scip));
1960 solvelimits->stallnodes += heurdata->nodesoffset;
1961 solvelimits->stallnodes -= heurdata->usednodes;
1962 solvelimits->stallnodes -= 100 * SCIPheurGetNCalls(heur);
1963 solvelimits->stallnodes = MIN(heurdata->maxnodes, solvelimits->stallnodes);
1964
1965 /* use a smaller budget if not all neighborhoods have been initialized yet */
1966 assert(heurdata->ninitneighborhoods >= 0);
1967 initfactor = (heurdata->nactiveneighborhoods - heurdata->ninitneighborhoods + 1.0) / (heurdata->nactiveneighborhoods + 1.0);
1968 solvelimits->stallnodes = (SCIP_Longint)(solvelimits->stallnodes * initfactor);
1969 solvelimits->nodelimit = (SCIP_Longint)(heurdata->maxnodes);
1970
1971 /* check whether we have enough nodes left to call subproblem solving */
1972 if( solvelimits->stallnodes < heurdata->targetnodes )
1973 *runagain = FALSE;
1974
1975 return SCIP_OKAY;
1976 }
1977
1978 /** return the bandit algorithm that should be used */
1979 static
getBandit(SCIP_HEURDATA * heurdata)1980 SCIP_BANDIT* getBandit(
1981 SCIP_HEURDATA* heurdata /**< heuristic data of the ALNS neighborhood */
1982 )
1983 {
1984 assert(heurdata != NULL);
1985 return heurdata->bandit;
1986 }
1987
1988 /** select a neighborhood depending on the selected bandit algorithm */
1989 static
selectNeighborhood(SCIP * scip,SCIP_HEURDATA * heurdata,int * neighborhoodidx)1990 SCIP_RETCODE selectNeighborhood(
1991 SCIP* scip, /**< SCIP data structure */
1992 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
1993 int* neighborhoodidx /**< pointer to store the selected neighborhood index */
1994 )
1995 {
1996 SCIP_BANDIT* bandit;
1997 assert(scip != NULL);
1998 assert(heurdata != NULL);
1999 assert(neighborhoodidx != NULL);
2000
2001 *neighborhoodidx = -1;
2002
2003 bandit = getBandit(heurdata);
2004
2005 SCIP_CALL( SCIPbanditSelect(bandit, neighborhoodidx) );
2006 assert(*neighborhoodidx >= 0);
2007
2008 return SCIP_OKAY;
2009 }
2010
2011 /** Calculate reward based on the selected reward measure */
2012 static
getReward(SCIP * scip,SCIP_HEURDATA * heurdata,NH_STATS * runstats,SCIP_Real * rewardptr)2013 SCIP_RETCODE getReward(
2014 SCIP* scip, /**< SCIP data structure */
2015 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
2016 NH_STATS* runstats, /**< run statistics */
2017 SCIP_Real* rewardptr /**< pointer to store the computed reward */
2018 )
2019 {
2020 SCIP_Real reward = 0.0;
2021 SCIP_Real effort;
2022 int ndiscretevars;
2023
2024 assert(rewardptr != NULL);
2025 assert(runstats->usednodes >= 0);
2026 assert(runstats->nfixings >= 0);
2027
2028 effort = runstats->usednodes / 100.0;
2029
2030 ndiscretevars = SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip);
2031 /* assume that every fixed variable linearly reduces the subproblem complexity */
2032 if( ndiscretevars > 0 )
2033 {
2034 effort = (1.0 - (runstats->nfixings / (SCIP_Real)ndiscretevars)) * effort;
2035 }
2036 assert(rewardptr != NULL);
2037
2038 /* a positive reward is only assigned if a new incumbent solution was found */
2039 if( runstats->nbestsolsfound > 0 )
2040 {
2041 SCIP_Real bestsolreward;
2042 SCIP_Real closedgapreward;
2043 SCIP_Real rewardcontrol = heurdata->rewardcontrol;
2044
2045 SCIP_Real lb;
2046 SCIP_Real ub;
2047
2048 /* the indicator function is simply 1.0 */
2049 bestsolreward = 1.0;
2050
2051 ub = runstats->newupperbound;
2052 lb = SCIPgetLowerbound(scip);
2053
2054 /* compute the closed gap reward */
2055 if( SCIPisEQ(scip, ub, lb) || SCIPisInfinity(scip, runstats->oldupperbound) )
2056 closedgapreward = 1.0;
2057 else
2058 {
2059 closedgapreward = (runstats->oldupperbound - ub) / (runstats->oldupperbound - lb);
2060 }
2061
2062 /* the reward is a convex combination of the best solution reward and the reward for the closed gap */
2063 reward = rewardcontrol * bestsolreward + (1.0 - rewardcontrol) * closedgapreward;
2064
2065 /* optionally, scale the reward by the involved effort */
2066 if( heurdata->scalebyeffort )
2067 reward /= (effort + 1.0);
2068
2069 /* add the baseline and rescale the reward into the interval [baseline, 1.0] */
2070 reward = heurdata->rewardbaseline + (1.0 - heurdata->rewardbaseline) * reward;
2071 }
2072 else
2073 {
2074 /* linearly decrease the reward based on the number of nodes spent */
2075 SCIP_Real maxeffort = heurdata->targetnodes;
2076 SCIP_Real usednodes = runstats->usednodes;
2077
2078 if( ndiscretevars > 0 )
2079 usednodes *= (1.0 - (runstats->nfixings / (SCIP_Real)ndiscretevars));
2080
2081 reward = heurdata->rewardbaseline - (usednodes) * heurdata->rewardbaseline / maxeffort;
2082 reward = MAX(0.0, reward);
2083 }
2084
2085 *rewardptr = reward;
2086 return SCIP_OKAY;
2087 }
2088
2089 /** update internal bandit algorithm statistics for future draws */
2090 static
updateBanditAlgorithm(SCIP * scip,SCIP_HEURDATA * heurdata,SCIP_Real reward,int neighborhoodidx)2091 SCIP_RETCODE updateBanditAlgorithm(
2092 SCIP* scip, /**< SCIP data structure */
2093 SCIP_HEURDATA* heurdata, /**< heuristic data of the ALNS neighborhood */
2094 SCIP_Real reward, /**< measured reward */
2095 int neighborhoodidx /**< the neighborhood that was chosen */
2096 )
2097 {
2098 SCIP_BANDIT* bandit;
2099 assert(scip != NULL);
2100 assert(heurdata != NULL);
2101 assert(neighborhoodidx >= 0);
2102 assert(neighborhoodidx < heurdata->nactiveneighborhoods);
2103
2104 bandit = getBandit(heurdata);
2105
2106 SCIPdebugMsg(scip, "Rewarding bandit algorithm action %d with reward %.2f\n", neighborhoodidx, reward);
2107 SCIP_CALL( SCIPbanditUpdate(bandit, neighborhoodidx, reward) );
2108
2109 return SCIP_OKAY;
2110 }
2111
2112 /** set up the sub-SCIP parameters, objective cutoff, and solution limits */
2113 static
setupSubScip(SCIP * scip,SCIP * subscip,SCIP_VAR ** subvars,SOLVELIMITS * solvelimits,SCIP_HEUR * heur,SCIP_Bool objchgd)2114 SCIP_RETCODE setupSubScip(
2115 SCIP* scip, /**< SCIP data structure */
2116 SCIP* subscip, /**< sub-SCIP data structure */
2117 SCIP_VAR** subvars, /**< array of sub-SCIP variables in the order of the main SCIP */
2118 SOLVELIMITS* solvelimits, /**< pointer to solving limits data structure */
2119 SCIP_HEUR* heur, /**< this heuristic */
2120 SCIP_Bool objchgd /**< did the objective change between the source and the target SCIP? */
2121 )
2122 {
2123 SCIP_HEURDATA* heurdata;
2124 SCIP_Real cutoff;
2125 SCIP_Real upperbound;
2126
2127 heurdata = SCIPheurGetData(heur);
2128
2129 /* do not abort subproblem on CTRL-C */
2130 SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", FALSE) );
2131
2132 /* disable output to console unless we are in debug mode */
2133 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
2134
2135 /* disable statistic timing inside sub SCIP */
2136 SCIP_CALL( SCIPsetBoolParam(subscip, "timing/statistictiming", FALSE) );
2137
2138 #ifdef ALNS_SUBSCIPOUTPUT
2139 SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
2140 SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 1) );
2141 /* enable statistic timing inside sub SCIP */
2142 SCIP_CALL( SCIPsetBoolParam(subscip, "timing/statistictiming", TRUE) );
2143 #endif
2144
2145 SCIP_CALL( SCIPsetIntParam(subscip, "limits/bestsol", heurdata->nsolslim) );
2146
2147 /* forbid recursive call of heuristics and separators solving subMIPs */
2148 if( ! heurdata->usesubscipheurs )
2149 {
2150 SCIP_CALL( SCIPsetSubscipsOff(subscip, TRUE) );
2151 }
2152
2153 /* disable cutting plane separation */
2154 SCIP_CALL( SCIPsetSeparating(subscip, SCIP_PARAMSETTING_OFF, TRUE) );
2155
2156 /* disable expensive presolving */
2157 SCIP_CALL( SCIPsetPresolving(subscip, SCIP_PARAMSETTING_FAST, TRUE) );
2158
2159 /* use best estimate node selection */
2160 if( SCIPfindNodesel(subscip, "estimate") != NULL && ! SCIPisParamFixed(subscip, "nodeselection/estimate/stdpriority") )
2161 {
2162 SCIP_CALL( SCIPsetIntParam(subscip, "nodeselection/estimate/stdpriority", INT_MAX/4) );
2163 }
2164
2165 /* use inference branching */
2166 if( SCIPfindBranchrule(subscip, "inference") != NULL && ! SCIPisParamFixed(subscip, "branching/inference/priority") )
2167 {
2168 SCIP_CALL( SCIPsetIntParam(subscip, "branching/inference/priority", INT_MAX/4) );
2169 }
2170
2171 /* enable conflict analysis and restrict conflict pool */
2172 if( ! SCIPisParamFixed(subscip, "conflict/enable") )
2173 {
2174 SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/enable", TRUE) );
2175 }
2176
2177 if( !SCIPisParamFixed(subscip, "conflict/useboundlp") )
2178 {
2179 SCIP_CALL( SCIPsetCharParam(subscip, "conflict/useboundlp", 'o') );
2180 }
2181
2182 if( ! SCIPisParamFixed(subscip, "conflict/maxstoresize") )
2183 {
2184 SCIP_CALL( SCIPsetIntParam(subscip, "conflict/maxstoresize", 100) );
2185 }
2186
2187 /* speed up sub-SCIP by not checking dual LP feasibility */
2188 SCIP_CALL( SCIPsetBoolParam(subscip, "lp/checkdualfeas", FALSE) );
2189
2190 /* employ a limit on the number of enforcement rounds in the quadratic constraint handlers; this fixes the issue that
2191 * sometimes the quadratic constraint handler needs hundreds or thousands of enforcement rounds to determine the
2192 * feasibility status of a single node without fractional branching candidates by separation (namely for uflquad
2193 * instances); however, the solution status of the sub-SCIP might get corrupted by this; hence no decutions shall be
2194 * made for the original SCIP
2195 */
2196 if( SCIPfindConshdlr(subscip, "quadratic") != NULL && ! SCIPisParamFixed(subscip, "constraints/quadratic/enfolplimit") )
2197 {
2198 SCIP_CALL( SCIPsetIntParam(subscip, "constraints/quadratic/enfolplimit", 10) );
2199 }
2200
2201 /* add an objective cutoff */
2202 if( ! SCIPisInfinity(scip, SCIPgetUpperbound(scip)) )
2203 {
2204 upperbound = SCIPgetUpperbound(scip) - SCIPsumepsilon(scip);
2205 if( ! SCIPisInfinity(scip, -1.0 * SCIPgetLowerbound(scip)) )
2206 {
2207 cutoff = (1 - heurdata->minimprove) * SCIPgetUpperbound(scip)
2208 + heurdata->minimprove * SCIPgetLowerbound(scip);
2209 }
2210 else
2211 {
2212 if( SCIPgetUpperbound(scip) >= 0 )
2213 cutoff = (1 - heurdata->minimprove) * SCIPgetUpperbound(scip);
2214 else
2215 cutoff = (1 + heurdata->minimprove) * SCIPgetUpperbound(scip);
2216 }
2217 cutoff = MIN(upperbound, cutoff);
2218
2219 if( SCIPisObjIntegral(scip) )
2220 cutoff = SCIPfloor(scip, cutoff);
2221
2222 SCIPdebugMsg(scip, "Sub-SCIP cutoff: %15.9" SCIP_REAL_FORMAT " (%15.9" SCIP_REAL_FORMAT " in original space)\n",
2223 cutoff, SCIPretransformObj(scip, cutoff));
2224
2225 /* if the objective changed between the source and the target SCIP, encode the cutoff as a constraint */
2226 if( ! objchgd )
2227 {
2228 SCIP_CALL(SCIPsetObjlimit(subscip, cutoff));
2229
2230 SCIPdebugMsg(scip, "Cutoff added as Objective Limit\n");
2231 }
2232 else
2233 {
2234 SCIP_CONS* objcons;
2235 int nvars;
2236 SCIP_VAR** vars;
2237 int i;
2238
2239 vars = SCIPgetVars(scip);
2240 nvars = SCIPgetNVars(scip);
2241
2242 SCIP_CALL( SCIPcreateConsLinear(subscip, &objcons, "objbound_of_origscip", 0, NULL, NULL, -SCIPinfinity(subscip), cutoff,
2243 TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
2244 for( i = 0; i < nvars; ++i)
2245 {
2246 if( ! SCIPisFeasZero(subscip, SCIPvarGetObj(vars[i])) )
2247 {
2248 assert(subvars[i] != NULL);
2249 SCIP_CALL( SCIPaddCoefLinear(subscip, objcons, subvars[i], SCIPvarGetObj(vars[i])) );
2250 }
2251 }
2252 SCIP_CALL( SCIPaddCons(subscip, objcons) );
2253 SCIP_CALL( SCIPreleaseCons(subscip, &objcons) );
2254
2255 SCIPdebugMsg(scip, "Cutoff added as constraint\n");
2256 }
2257 }
2258
2259 /* set solve limits for sub-SCIP */
2260 SCIP_CALL( setLimits(subscip, solvelimits) );
2261
2262 /* change random seed of sub-SCIP */
2263 if( heurdata->subsciprandseeds )
2264 {
2265 SCIP_CALL( SCIPsetIntParam(subscip, "randomization/randomseedshift", (int)SCIPheurGetNCalls(heur)) );
2266 }
2267
2268 SCIPdebugMsg(scip, "Solve Limits: %lld (%lld) nodes (stall nodes), %.1f sec., %d sols\n",
2269 solvelimits->nodelimit, solvelimits->stallnodes, solvelimits->timelimit, heurdata->nsolslim);
2270
2271 return SCIP_OKAY;
2272 }
2273
2274 /** execution method of primal heuristic */
2275 static
SCIP_DECL_HEUREXEC(heurExecAlns)2276 SCIP_DECL_HEUREXEC(heurExecAlns)
2277 { /*lint --e{715}*/
2278 SCIP_HEURDATA* heurdata;
2279 SCIP_VAR** varbuf;
2280 SCIP_Real* valbuf;
2281 SCIP_VAR** vars;
2282 SCIP_VAR** subvars;
2283 NH_STATS runstats[NNEIGHBORHOODS];
2284 SCIP_STATUS subscipstatus[NNEIGHBORHOODS];
2285 SCIP* subscip = NULL;
2286
2287 int nfixings;
2288 int nvars;
2289 int neighborhoodidx;
2290 int ntries;
2291 SCIP_Bool tryagain;
2292 NH* neighborhood;
2293 SOLVELIMITS solvelimits;
2294 SCIP_Bool success;
2295 SCIP_Bool run;
2296 SCIP_Bool allrewardsmode;
2297 SCIP_Real rewards[NNEIGHBORHOODS];
2298 int banditidx;
2299 int i;
2300
2301 heurdata = SCIPheurGetData(heur);
2302 assert(heurdata != NULL);
2303
2304 *result = SCIP_DIDNOTRUN;
2305
2306 if( heurdata->nactiveneighborhoods == 0 )
2307 return SCIP_OKAY;
2308
2309 /* wait for a sufficient number of nodes since last incumbent solution */
2310 if( SCIPgetDepth(scip) > 0 && SCIPgetBestSol(scip) != NULL
2311 && (SCIPgetNNodes(scip) - SCIPsolGetNodenum(SCIPgetBestSol(scip))) < heurdata->waitingnodes )
2312 {
2313 SCIPdebugMsg(scip, "Waiting nodes not satisfied\n");
2314 return SCIP_OKAY;
2315 }
2316
2317 run = TRUE;
2318 /* check if budget allows a run of the next selected neighborhood */
2319 SCIP_CALL( determineLimits(scip, heur, &solvelimits, &run) );
2320 SCIPdebugMsg(scip, "Budget check: %" SCIP_LONGINT_FORMAT " (%" SCIP_LONGINT_FORMAT ") %s\n", solvelimits.nodelimit, heurdata->targetnodes, run ? "passed" : "must wait");
2321
2322 if( ! run )
2323 return SCIP_OKAY;
2324
2325 /* delay the heuristic if local reduced costs should be used for generic variable unfixing */
2326 if( heurdata->uselocalredcost && (nodeinfeasible || ! SCIPhasCurrentNodeLP(scip) || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL) )
2327 {
2328 *result = SCIP_DELAYED;
2329
2330 return SCIP_OKAY;
2331 }
2332
2333 allrewardsmode = heurdata->rewardfile != NULL;
2334
2335 /* apply some other rules for a fair all rewards mode; in normal execution mode, neighborhoods are iterated through */
2336 if( allrewardsmode )
2337 {
2338 /* most neighborhoods require an incumbent solution */
2339 if( SCIPgetNSols(scip) < 2 )
2340 {
2341 SCIPdebugMsg(scip, "Not enough solutions for all rewards mode\n");
2342 return SCIP_OKAY;
2343 }
2344
2345 /* if the node is infeasible, or has no LP solution, which is required by some neighborhoods
2346 * if we are not in all rewards mode, the neighborhoods delay themselves individually
2347 */
2348 if( nodeinfeasible || ! SCIPhasCurrentNodeLP(scip) || SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
2349 {
2350 SCIPdebugMsg(scip, "Delay ALNS heuristic until a feasible node with optimally solved LP relaxation\n");
2351 *result = SCIP_DELAYED;
2352 return SCIP_OKAY;
2353 }
2354 }
2355
2356 /* use the neighborhood that requested a delay or select the next neighborhood to run based on the selected bandit algorithm */
2357 if( heurdata->currneighborhood >= 0 )
2358 {
2359 assert(! allrewardsmode);
2360 banditidx = heurdata->currneighborhood;
2361 SCIPdebugMsg(scip, "Select delayed neighborhood %d (was delayed %d times)\n", banditidx, heurdata->ndelayedcalls);
2362 }
2363 else
2364 {
2365 SCIP_CALL( selectNeighborhood(scip, heurdata, &banditidx) );
2366 SCIPdebugMsg(scip, "Selected neighborhood %d with bandit algorithm\n", banditidx);
2367 }
2368
2369 /* in all rewards mode, we simply loop over all heuristics */
2370 if( ! allrewardsmode )
2371 neighborhoodidx = banditidx;
2372 else
2373 neighborhoodidx = 0;
2374
2375 assert(0 <= neighborhoodidx && neighborhoodidx < NNEIGHBORHOODS);
2376 assert(heurdata->nactiveneighborhoods > neighborhoodidx);
2377
2378 /* allocate memory for variable fixings buffer */
2379 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2380 SCIP_CALL( SCIPallocBufferArray(scip, &varbuf, nvars) );
2381 SCIP_CALL( SCIPallocBufferArray(scip, &valbuf, nvars) );
2382 SCIP_CALL( SCIPallocBufferArray(scip, &subvars, nvars) );
2383
2384 /* initialize neighborhood statistics for a run */
2385 ntries = 1;
2386 do
2387 {
2388 SCIP_HASHMAP* varmapf;
2389 SCIP_EVENTHDLR* eventhdlr;
2390 SCIP_EVENTDATA eventdata;
2391 char probnamesuffix[SCIP_MAXSTRLEN];
2392 SCIP_Real allfixingrate;
2393 int ndomchgs;
2394 int nchgobjs;
2395 int naddedconss;
2396 int v;
2397 SCIP_RETCODE retcode;
2398 SCIP_RESULT fixresult;
2399
2400 tryagain = FALSE;
2401 neighborhood = heurdata->neighborhoods[neighborhoodidx];
2402 SCIPdebugMsg(scip, "Running '%s' neighborhood %d\n", neighborhood->name, neighborhoodidx);
2403
2404 initRunStats(scip, &runstats[neighborhoodidx]);
2405 rewards[neighborhoodidx] = 0.0;
2406
2407 subscipstatus[neighborhoodidx] = SCIP_STATUS_UNKNOWN;
2408 SCIP_CALL( SCIPstartClock(scip, neighborhood->stats.setupclock) );
2409
2410 /* determine variable fixings and objective coefficients of this neighborhood */
2411 SCIP_CALL( neighborhoodFixVariables(scip, heurdata, neighborhood, varbuf, valbuf, &nfixings, &fixresult) );
2412
2413 SCIPdebugMsg(scip, "Fix %d/%d variables, result code %d\n", nfixings, nvars,fixresult);
2414
2415 /* Fixing was not successful, either because the fixing rate was not reached (and no additional variable
2416 * prioritization was used), or the neighborhood requested a delay, e.g., because no LP relaxation solution exists
2417 * at the current node
2418 *
2419 * The ALNS heuristic keeps a delayed neighborhood active and delays itself.
2420 */
2421 if( fixresult != SCIP_SUCCESS )
2422 {
2423 SCIP_CALL( SCIPstopClock(scip, neighborhood->stats.setupclock) );
2424
2425 /* to determine all rewards, we cannot delay neighborhoods */
2426 if( allrewardsmode )
2427 {
2428 if( ntries == heurdata->nactiveneighborhoods )
2429 break;
2430
2431 neighborhoodidx = (neighborhoodidx + 1) % heurdata->nactiveneighborhoods;
2432 ntries++;
2433 tryagain = TRUE;
2434
2435 continue;
2436 }
2437
2438 /* delay the heuristic along with the selected neighborhood
2439 *
2440 * if the neighborhood has been delayed for too many consecutive calls, the delay is treated as a failure */
2441 if( fixresult == SCIP_DELAYED )
2442 {
2443 if( heurdata->ndelayedcalls > (SCIPheurGetFreq(heur) / 4 + 1) )
2444 {
2445 resetCurrentNeighborhood(heurdata);
2446
2447 /* use SCIP_DIDNOTFIND to penalize the neighborhood with a bad reward */
2448 fixresult = SCIP_DIDNOTFIND;
2449 }
2450 else if( heurdata->currneighborhood == -1 )
2451 {
2452 heurdata->currneighborhood = neighborhoodidx;
2453 heurdata->ndelayedcalls = 1;
2454 }
2455 else
2456 {
2457 heurdata->ndelayedcalls++;
2458 }
2459 }
2460
2461 if( fixresult == SCIP_DIDNOTRUN )
2462 {
2463 if( ntries < heurdata->nactiveneighborhoods )
2464 {
2465 SCIP_CALL( updateBanditAlgorithm(scip, heurdata, 0.0, neighborhoodidx) );
2466 SCIP_CALL( selectNeighborhood(scip, heurdata, &neighborhoodidx) );
2467 ntries++;
2468 tryagain = TRUE;
2469
2470 SCIPdebugMsg(scip, "Neighborhood cannot run -> try next neighborhood %d\n", neighborhoodidx);
2471 continue;
2472 }
2473 else
2474 break;
2475 }
2476
2477 assert(fixresult == SCIP_DIDNOTFIND || fixresult == SCIP_DELAYED);
2478 *result = fixresult;
2479 break;
2480 }
2481
2482 *result = SCIP_DIDNOTFIND;
2483
2484 neighborhood->stats.nfixings += nfixings;
2485 runstats[neighborhoodidx].nfixings = nfixings;
2486
2487 SCIP_CALL( SCIPcreate(&subscip) );
2488 SCIP_CALL( SCIPhashmapCreate(&varmapf, SCIPblkmem(scip), nvars) );
2489 (void) SCIPsnprintf(probnamesuffix, SCIP_MAXSTRLEN, "alns_%s", neighborhood->name);
2490
2491 /* todo later: run global propagation for this set of fixings */
2492 SCIP_CALL( SCIPcopyLargeNeighborhoodSearch(scip, subscip, varmapf, probnamesuffix, varbuf, valbuf, nfixings, FALSE, heurdata->copycuts, &success, NULL) );
2493
2494 /* store sub-SCIP variables in array for faster access */
2495 for( v = 0; v < nvars; ++v )
2496 {
2497 subvars[v] = (SCIP_VAR*)SCIPhashmapGetImage(varmapf, (void *)vars[v]);
2498 }
2499
2500 SCIPhashmapFree(&varmapf);
2501
2502 /* let the neighborhood add additional constraints, or restrict domains */
2503 SCIP_CALL( neighborhoodChangeSubscip(scip, subscip, neighborhood, subvars, &ndomchgs, &nchgobjs, &naddedconss, &success) );
2504
2505 if( ! success )
2506 {
2507 SCIP_CALL( SCIPstopClock(scip, neighborhood->stats.setupclock) );
2508
2509 if( ! allrewardsmode || ntries == heurdata->nactiveneighborhoods )
2510 break;
2511
2512 neighborhoodidx = (neighborhoodidx + 1) % heurdata->nactiveneighborhoods;
2513 ntries++;
2514 tryagain = TRUE;
2515
2516 SCIP_CALL( SCIPfree(&subscip) );
2517
2518 continue;
2519 }
2520
2521 /* set up sub-SCIP parameters */
2522 SCIP_CALL( setupSubScip(scip, subscip, subvars, &solvelimits, heur, nchgobjs > 0) );
2523
2524 /* copy the necessary data into the event data to create new solutions */
2525 eventdata.nodelimit = solvelimits.nodelimit; /*lint !e644*/
2526 eventdata.lplimfac = heurdata->lplimfac;
2527 eventdata.heur = heur;
2528 eventdata.sourcescip = scip;
2529 eventdata.subvars = subvars;
2530 eventdata.runstats = &runstats[neighborhoodidx];
2531 eventdata.allrewardsmode = allrewardsmode;
2532
2533 /* include an event handler to transfer solutions into the main SCIP */
2534 SCIP_CALL( SCIPincludeEventhdlrBasic(subscip, &eventhdlr, EVENTHDLR_NAME, EVENTHDLR_DESC, eventExecAlns, NULL) );
2535
2536 /* transform the problem before catching the events */
2537 SCIP_CALL( SCIPtransformProb(subscip) );
2538 SCIP_CALL( SCIPcatchEvent(subscip, SCIP_EVENTTYPE_ALNS, eventhdlr, &eventdata, NULL) );
2539
2540 SCIP_CALL( SCIPstopClock(scip, neighborhood->stats.setupclock) );
2541
2542 SCIP_CALL( SCIPstartClock(scip, neighborhood->stats.submipclock) );
2543
2544 /* set up sub-SCIP and run presolving */
2545 retcode = SCIPpresolve(subscip);
2546 if( retcode != SCIP_OKAY )
2547 {
2548 SCIPwarningMessage(scip, "Error while presolving subproblem in ALNS heuristic; sub-SCIP terminated with code <%d>\n", retcode);
2549 SCIP_CALL( SCIPstopClock(scip, neighborhood->stats.submipclock) );
2550
2551 SCIPABORT(); /*lint --e{527}*/
2552 break;
2553 }
2554
2555 /* was presolving successful enough regarding fixings? otherwise, terminate */
2556 allfixingrate = (SCIPgetNOrigVars(subscip) - SCIPgetNVars(subscip)) / (SCIP_Real)SCIPgetNOrigVars(subscip);
2557
2558 /* additional variables added in presolving may lead to the subSCIP having more variables than the original */
2559 allfixingrate = MAX(allfixingrate, 0.0);
2560
2561 if( allfixingrate >= neighborhood->fixingrate.targetfixingrate / 2.0 )
2562 {
2563 /* run sub-SCIP for the given budget, and collect statistics */
2564 SCIP_CALL_ABORT( SCIPsolve(subscip) );
2565 }
2566 else
2567 {
2568 SCIPdebugMsg(scip, "Fixed only %.3f of all variables after presolving -> do not solve sub-SCIP\n", allfixingrate);
2569 }
2570
2571 #ifdef ALNS_SUBSCIPOUTPUT
2572 SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2573 #endif
2574
2575 SCIP_CALL( SCIPstopClock(scip, neighborhood->stats.submipclock) );
2576
2577 /* update statistics based on the sub-SCIP run results */
2578 updateRunStats(&runstats[neighborhoodidx], subscip);
2579 subscipstatus[neighborhoodidx] = SCIPgetStatus(subscip);
2580 SCIPdebugMsg(scip, "Status of sub-SCIP run: %d\n", subscipstatus[neighborhoodidx]);
2581
2582 SCIP_CALL( getReward(scip, heurdata, &runstats[neighborhoodidx], &rewards[neighborhoodidx]) );
2583
2584 /* in all rewards mode, continue with the next neighborhood */
2585 if( allrewardsmode && ntries < heurdata->nactiveneighborhoods )
2586 {
2587 neighborhoodidx = (neighborhoodidx + 1) % heurdata->nactiveneighborhoods;
2588 ntries++;
2589 tryagain = TRUE;
2590
2591 SCIP_CALL( SCIPfree(&subscip) );
2592 }
2593 }
2594 while( tryagain && ! SCIPisStopped(scip) );
2595
2596 if( subscip != NULL )
2597 {
2598 SCIP_CALL( SCIPfree(&subscip) );
2599 }
2600
2601 SCIPfreeBufferArray(scip, &subvars);
2602 SCIPfreeBufferArray(scip, &valbuf);
2603 SCIPfreeBufferArray(scip, &varbuf);
2604
2605 /* update bandit index that may have changed unless we are in all rewards mode */
2606 if( ! allrewardsmode )
2607 banditidx = neighborhoodidx;
2608
2609 if( *result != SCIP_DELAYED )
2610 {
2611 /* decrease the number of neighborhoods that have not been initialized */
2612 if( neighborhood->stats.nruns == 0 )
2613 --heurdata->ninitneighborhoods;
2614
2615 heurdata->usednodes += runstats[banditidx].usednodes;
2616
2617 /* determine the success of this neighborhood, and update the target fixing rate for the next time */
2618 updateNeighborhoodStats(&runstats[banditidx], heurdata->neighborhoods[banditidx], subscipstatus[banditidx]);
2619
2620 /* adjust the fixing rate for this neighborhood
2621 * make no adjustments in all rewards mode, because this only affects 1 of 8 heuristics
2622 */
2623 if( heurdata->adjustfixingrate && ! allrewardsmode )
2624 {
2625 SCIPdebugMsg(scip, "Update fixing rate: %.2f\n", heurdata->neighborhoods[banditidx]->fixingrate.targetfixingrate);
2626 updateFixingRate(heurdata->neighborhoods[banditidx], subscipstatus[banditidx], &runstats[banditidx]);
2627 SCIPdebugMsg(scip, "New fixing rate: %.2f\n", heurdata->neighborhoods[banditidx]->fixingrate.targetfixingrate);
2628 }
2629 /* similarly, update the minimum improvement for the ALNS heuristic */
2630 if( heurdata->adjustminimprove )
2631 {
2632 SCIPdebugMsg(scip, "Update Minimum Improvement: %.4f\n", heurdata->minimprove);
2633 updateMinimumImprovement(heurdata, subscipstatus[banditidx], &runstats[banditidx]);
2634 SCIPdebugMsg(scip, "--> %.4f\n", heurdata->minimprove);
2635 }
2636
2637 /* update the target node limit based on the status of the selected algorithm */
2638 if( heurdata->adjusttargetnodes && SCIPheurGetNCalls(heur) >= heurdata->nactiveneighborhoods )
2639 {
2640 updateTargetNodeLimit(heurdata, &runstats[banditidx], subscipstatus[banditidx]);
2641 }
2642
2643 /* update the bandit algorithm by the measured reward */
2644 SCIP_CALL( updateBanditAlgorithm(scip, heurdata, rewards[banditidx], banditidx) );
2645
2646 resetCurrentNeighborhood(heurdata);
2647 }
2648
2649 /* write single, measured rewards and the bandit index to the reward file */
2650 if( allrewardsmode )
2651 {
2652 for( i = 0; i < heurdata->nactiveneighborhoods; ++i )
2653 {
2654 fprintf(heurdata->rewardfile, "%.4f,", rewards[i]);
2655 }
2656 fprintf(heurdata->rewardfile, "%d\n", banditidx);
2657 }
2658
2659 return SCIP_OKAY;
2660 }
2661
2662 /** callback to collect variable fixings of RENS */
2663 static
DECL_VARFIXINGS(varFixingsRens)2664 DECL_VARFIXINGS(varFixingsRens)
2665 { /*lint --e{715}*/
2666 int nbinvars;
2667 int nintvars;
2668 SCIP_VAR** vars;
2669 int i;
2670 assert(scip != NULL);
2671 assert(varbuf != NULL);
2672 assert(nfixings != NULL);
2673 assert(valbuf != NULL);
2674
2675 *result = SCIP_DELAYED;
2676
2677 if( ! SCIPhasCurrentNodeLP(scip) )
2678 return SCIP_OKAY;
2679 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
2680 return SCIP_OKAY;
2681
2682 *result = SCIP_DIDNOTRUN;
2683
2684 /* get variable information */
2685 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
2686
2687 /* return if no binary or integer variables are present */
2688 if( nbinvars + nintvars == 0 )
2689 return SCIP_OKAY;
2690
2691 /* loop over binary and integer variables; determine those that should be fixed in the sub-SCIP */
2692 for( i = 0; i < nbinvars + nintvars; ++i )
2693 {
2694 SCIP_VAR* var = vars[i];
2695 SCIP_Real lpsolval = SCIPvarGetLPSol(var);
2696 assert((i < nbinvars && SCIPvarIsBinary(var)) || (i >= nbinvars && SCIPvarIsIntegral(var)));
2697
2698 /* fix all binary and integer variables with integer LP solution value */
2699 if( SCIPisFeasIntegral(scip, lpsolval) )
2700 tryAdd2variableBuffer(scip, var, lpsolval, varbuf, valbuf, nfixings, TRUE);
2701 }
2702
2703 *result = SCIP_SUCCESS;
2704
2705 return SCIP_OKAY;
2706 }
2707
2708 /** callback for RENS subproblem changes */
2709 static
DECL_CHANGESUBSCIP(changeSubscipRens)2710 DECL_CHANGESUBSCIP(changeSubscipRens)
2711 { /*lint --e{715}*/
2712 SCIP_VAR** vars;
2713 int nintvars;
2714 int nbinvars;
2715 int i;
2716
2717 assert(SCIPhasCurrentNodeLP(sourcescip));
2718 assert(SCIPgetLPSolstat(sourcescip) == SCIP_LPSOLSTAT_OPTIMAL);
2719
2720 /* get variable information */
2721 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
2722
2723 /* restrict bounds of integer variables with fractional solution value */
2724 for( i = nbinvars; i < nbinvars + nintvars; ++i )
2725 {
2726 SCIP_VAR* var = vars[i];
2727 SCIP_Real lpsolval = SCIPgetSolVal(sourcescip, NULL, var);
2728
2729 if( subvars[i] == NULL )
2730 continue;
2731
2732 if( ! SCIPisFeasIntegral(sourcescip, lpsolval) )
2733 {
2734 SCIP_Real newlb = SCIPfloor(sourcescip, lpsolval);
2735 SCIP_Real newub = newlb + 1.0;
2736
2737 /* only count this as a domain change if the new lower and upper bound are a further restriction */
2738 if( newlb > SCIPvarGetLbGlobal(subvars[i]) + 0.5 || newub < SCIPvarGetUbGlobal(subvars[i]) - 0.5 )
2739 {
2740 SCIP_CALL( SCIPchgVarLbGlobal(targetscip, subvars[i], newlb) );
2741 SCIP_CALL( SCIPchgVarUbGlobal(targetscip, subvars[i], newub) );
2742 (*ndomchgs)++;
2743 }
2744 }
2745 }
2746
2747 *success = TRUE;
2748
2749 return SCIP_OKAY;
2750 }
2751
2752 /** collect fixings by matching solution values in a collection of solutions for all binary and integer variables,
2753 * or for a custom set of variables
2754 */
2755 static
fixMatchingSolutionValues(SCIP * scip,SCIP_SOL ** sols,int nsols,SCIP_VAR ** vars,int nvars,SCIP_VAR ** varbuf,SCIP_Real * valbuf,int * nfixings)2756 SCIP_RETCODE fixMatchingSolutionValues(
2757 SCIP* scip, /**< SCIP data structure */
2758 SCIP_SOL** sols, /**< array of 2 or more solutions. It is okay for the array to contain one element
2759 * equal to NULL to represent the current LP solution */
2760 int nsols, /**< number of solutions in the array */
2761 SCIP_VAR** vars, /**< variable array for which solution values must agree */
2762 int nvars, /**< number of variables, or -1 for all binary and integer variables */
2763 SCIP_VAR** varbuf, /**< buffer storage for variable fixings */
2764 SCIP_Real* valbuf, /**< buffer storage for fixing values */
2765 int* nfixings /**< pointer to store the number of fixings */
2766 )
2767 {
2768 int v;
2769 int nbinintvars;
2770 SCIP_SOL* firstsol;
2771
2772 assert(scip != NULL);
2773 assert(sols != NULL);
2774 assert(nsols >= 2);
2775 assert(varbuf != NULL);
2776 assert(valbuf != NULL);
2777 assert(nfixings != NULL);
2778 assert(*nfixings == 0);
2779
2780 if( nvars == -1 || vars == NULL )
2781 {
2782 int nbinvars;
2783 int nintvars;
2784 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
2785 nbinintvars = nbinvars + nintvars;
2786 nvars = nbinintvars;
2787 }
2788 firstsol = sols[0];
2789 assert(nvars > 0);
2790
2791 /* loop over integer and binary variables and check if their solution values match in all solutions */
2792 for( v = 0; v < nvars; ++v )
2793 {
2794 SCIP_Real solval;
2795 SCIP_VAR* var;
2796 int s;
2797
2798 var = vars[v];
2799 assert((v < SCIPgetNBinVars(scip) && SCIPvarIsBinary(var)) || (v >= SCIPgetNBinVars(scip) && SCIPvarIsIntegral(var)));
2800 solval = SCIPgetSolVal(scip, firstsol, var);
2801
2802 /* determine if solution values match in all given solutions */
2803 for( s = 1; s < nsols; ++s )
2804 {
2805 SCIP_Real solval2 = SCIPgetSolVal(scip, sols[s], var);
2806 if( ! SCIPisEQ(scip, solval, solval2) )
2807 break;
2808 }
2809
2810 /* if we did not break early, all solutions agree on the solution value of this variable */
2811 if( s == nsols )
2812 {
2813 tryAdd2variableBuffer(scip, var, solval, varbuf, valbuf, nfixings, TRUE);
2814 }
2815 }
2816
2817 return SCIP_OKAY;
2818 }
2819
2820 /** callback to collect variable fixings of RINS */
2821 static
DECL_VARFIXINGS(varFixingsRins)2822 DECL_VARFIXINGS(varFixingsRins)
2823 {
2824 /*lint --e{715}*/
2825 int nbinvars;
2826 int nintvars;
2827 SCIP_VAR** vars;
2828 SCIP_SOL* incumbent;
2829 SCIP_SOL* sols[2];
2830 assert(scip != NULL);
2831 assert(varbuf != NULL);
2832 assert(nfixings != NULL);
2833 assert(valbuf != NULL);
2834
2835 *result = SCIP_DELAYED;
2836
2837 if( ! SCIPhasCurrentNodeLP(scip) )
2838 return SCIP_OKAY;
2839 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
2840 return SCIP_OKAY;
2841
2842 *result = SCIP_DIDNOTRUN;
2843
2844 incumbent = SCIPgetBestSol(scip);
2845 if( incumbent == NULL )
2846 return SCIP_OKAY;
2847
2848 if( SCIPsolGetOrigin(incumbent) == SCIP_SOLORIGIN_ORIGINAL )
2849 return SCIP_OKAY;
2850
2851 /* get variable information */
2852 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
2853
2854 /* return if no binary or integer variables are present */
2855 if( nbinvars + nintvars == 0 )
2856 return SCIP_OKAY;
2857
2858 sols[0] = NULL;
2859 sols[1] = incumbent;
2860
2861 SCIP_CALL( fixMatchingSolutionValues(scip, sols, 2, vars, nbinvars + nintvars, varbuf, valbuf, nfixings) );
2862
2863 *result = SCIP_SUCCESS;
2864
2865 return SCIP_OKAY;
2866 }
2867
2868 /** initialization callback for crossover when a new problem is read */
2869 static
DECL_NHINIT(nhInitCrossover)2870 DECL_NHINIT(nhInitCrossover)
2871 { /*lint --e{715}*/
2872 DATA_CROSSOVER* data;
2873
2874 data = neighborhood->data.crossover;
2875 assert(data != NULL);
2876
2877 if( data->rng != NULL )
2878 SCIPfreeRandom(scip, &data->rng);
2879
2880 data->selsol = NULL;
2881
2882 SCIP_CALL( SCIPcreateRandom(scip, &data->rng, CROSSOVERSEED + (unsigned int)SCIPgetNVars(scip), TRUE) );
2883
2884 return SCIP_OKAY;
2885 }
2886
2887 /** deinitialization callback for crossover when exiting a problem */
2888 static
DECL_NHEXIT(nhExitCrossover)2889 DECL_NHEXIT(nhExitCrossover)
2890 { /*lint --e{715}*/
2891 DATA_CROSSOVER* data;
2892 data = neighborhood->data.crossover;
2893
2894 assert(neighborhood != NULL);
2895 assert(data->rng != NULL);
2896
2897 SCIPfreeRandom(scip, &data->rng);
2898
2899 return SCIP_OKAY;
2900 }
2901
2902 /** deinitialization callback for crossover before SCIP is freed */
2903 static
DECL_NHFREE(nhFreeCrossover)2904 DECL_NHFREE(nhFreeCrossover)
2905 { /*lint --e{715}*/
2906 assert(neighborhood->data.crossover != NULL);
2907 SCIPfreeBlockMemory(scip, &neighborhood->data.crossover);
2908
2909 return SCIP_OKAY;
2910 }
2911
2912 /** callback to collect variable fixings of crossover */
2913 static
DECL_VARFIXINGS(varFixingsCrossover)2914 DECL_VARFIXINGS(varFixingsCrossover)
2915 { /*lint --e{715}*/
2916 DATA_CROSSOVER* data;
2917 SCIP_RANDNUMGEN* rng;
2918 SCIP_SOL** sols;
2919 SCIP_SOL** scipsols;
2920 int nsols;
2921 int lastdraw;
2922 assert(scip != NULL);
2923 assert(varbuf != NULL);
2924 assert(nfixings != NULL);
2925 assert(valbuf != NULL);
2926
2927 data = neighborhood->data.crossover;
2928
2929 assert(data != NULL);
2930 nsols = data->nsols;
2931 data->selsol = NULL;
2932
2933 *result = SCIP_DIDNOTRUN;
2934
2935 /* return if the pool has not enough solutions */
2936 if( nsols > SCIPgetNSols(scip) )
2937 return SCIP_OKAY;
2938
2939 /* return if no binary or integer variables are present */
2940 if( SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip) == 0 )
2941 return SCIP_OKAY;
2942
2943 rng = data->rng;
2944 lastdraw = SCIPgetNSols(scip);
2945 SCIP_CALL( SCIPallocBufferArray(scip, &sols, nsols) );
2946 scipsols = SCIPgetSols(scip);
2947
2948 /* draw as many solutions from the pool as required by crossover, biased towards
2949 * better solutions; therefore, the sorting of the solutions by objective is implicitly used
2950 */
2951 while( nsols > 0 )
2952 {
2953 /* no need for randomization anymore, exactly nsols many solutions remain for the selection */
2954 if( lastdraw == nsols )
2955 {
2956 int s;
2957
2958 /* fill the remaining slots 0,...,nsols - 1 by the solutions at the same places */
2959 for( s = 0; s < nsols; ++s )
2960 sols[s] = scipsols[s];
2961
2962 nsols = 0;
2963 }
2964 else
2965 {
2966 int nextdraw;
2967
2968 assert(nsols < lastdraw);
2969
2970 /* draw from the lastdraw - nsols many solutions nsols - 1, ... lastdraw - 1 such that nsols many solution */
2971 nextdraw = SCIPrandomGetInt(rng, nsols - 1, lastdraw - 1);
2972 assert(nextdraw >= 0);
2973
2974 sols[nsols - 1] = scipsols[nextdraw];
2975 nsols--;
2976 lastdraw = nextdraw;
2977 }
2978 }
2979
2980 SCIP_CALL( fixMatchingSolutionValues(scip, sols, data->nsols, NULL, -1, varbuf, valbuf, nfixings) );
2981
2982 /* store best selected solution as reference solution */
2983 data->selsol = sols[0];
2984 assert(data->selsol != NULL);
2985
2986 *result = SCIP_SUCCESS;
2987
2988 SCIPfreeBufferArray(scip, &sols);
2989
2990 return SCIP_OKAY;
2991 }
2992
2993 /** callback for crossover reference solution */
2994 static
DECL_NHREFSOL(nhRefsolCrossover)2995 DECL_NHREFSOL(nhRefsolCrossover)
2996 { /*lint --e{715}*/
2997 DATA_CROSSOVER* data;
2998
2999 data = neighborhood->data.crossover;
3000
3001 if( data->selsol != NULL )
3002 {
3003 *solptr = data->selsol;
3004 *result = SCIP_SUCCESS;
3005 }
3006 else
3007 {
3008 *result = SCIP_DIDNOTFIND;
3009 }
3010
3011 return SCIP_OKAY;
3012 }
3013
3014 /** initialization callback for mutation when a new problem is read */
3015 static
DECL_NHINIT(nhInitMutation)3016 DECL_NHINIT(nhInitMutation)
3017 { /*lint --e{715}*/
3018 DATA_MUTATION* data;
3019 assert(scip != NULL);
3020 assert(neighborhood != NULL);
3021
3022 SCIP_CALL( SCIPallocBlockMemory(scip, &neighborhood->data.mutation) );
3023
3024 data = neighborhood->data.mutation;
3025 assert(data != NULL);
3026
3027 SCIP_CALL( SCIPcreateRandom(scip, &data->rng, MUTATIONSEED + (unsigned int)SCIPgetNVars(scip), TRUE) );
3028
3029 return SCIP_OKAY;
3030 }
3031
3032 /** deinitialization callback for mutation when exiting a problem */
3033 static
DECL_NHEXIT(nhExitMutation)3034 DECL_NHEXIT(nhExitMutation)
3035 { /*lint --e{715}*/
3036 DATA_MUTATION* data;
3037 assert(scip != NULL);
3038 assert(neighborhood != NULL);
3039 data = neighborhood->data.mutation;
3040 assert(data != NULL);
3041
3042 SCIPfreeRandom(scip, &data->rng);
3043
3044 SCIPfreeBlockMemory(scip, &neighborhood->data.mutation);
3045
3046 return SCIP_OKAY;
3047 }
3048
3049 /** callback to collect variable fixings of mutation */
3050 static
DECL_VARFIXINGS(varFixingsMutation)3051 DECL_VARFIXINGS(varFixingsMutation)
3052 { /*lint --e{715}*/
3053 SCIP_RANDNUMGEN* rng;
3054
3055 SCIP_VAR** vars;
3056 SCIP_VAR** varscpy;
3057 int i;
3058 int nvars;
3059 int nbinvars;
3060 int nintvars;
3061 int nbinintvars;
3062 int ntargetfixings;
3063 SCIP_SOL* incumbentsol;
3064 SCIP_Real targetfixingrate;
3065
3066 assert(scip != NULL);
3067 assert(neighborhood != NULL);
3068 assert(neighborhood->data.mutation != NULL);
3069 assert(neighborhood->data.mutation->rng != NULL);
3070 rng = neighborhood->data.mutation->rng;
3071
3072 *result = SCIP_DIDNOTRUN;
3073
3074 /* get the problem variables */
3075 SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) );
3076
3077 nbinintvars = nbinvars + nintvars;
3078 if( nbinintvars == 0 )
3079 return SCIP_OKAY;
3080
3081 incumbentsol = SCIPgetBestSol(scip);
3082 if( incumbentsol == NULL )
3083 return SCIP_OKAY;
3084
3085 targetfixingrate = neighborhood->fixingrate.targetfixingrate;
3086 ntargetfixings = (int)(targetfixingrate * nbinintvars) + 1;
3087
3088 /* don't continue if number of discrete variables is too small to reach target fixing rate */
3089 if( nbinintvars <= ntargetfixings )
3090 return SCIP_OKAY;
3091
3092 *result = SCIP_DIDNOTFIND;
3093
3094 /* copy variables into a buffer array */
3095 SCIP_CALL( SCIPduplicateBufferArray(scip, &varscpy, vars, nbinintvars) );
3096
3097 /* partially perturb the array until the number of target fixings is reached */
3098 for( i = 0; *nfixings < ntargetfixings && i < nbinintvars; ++i )
3099 {
3100 int randint = SCIPrandomGetInt(rng, i, nbinintvars - 1);
3101 assert(randint < nbinintvars);
3102
3103 if( randint > i )
3104 {
3105 SCIPswapPointers((void**)&varscpy[i], (void**)&varscpy[randint]);
3106 }
3107 /* copy the selected variables and their solution values into the buffer */
3108 tryAdd2variableBuffer(scip, varscpy[i], SCIPgetSolVal(scip, incumbentsol, varscpy[i]), varbuf, valbuf, nfixings, TRUE);
3109 }
3110
3111 assert(i == nbinintvars || *nfixings == ntargetfixings);
3112
3113 /* Not reaching the number of target fixings means that there is a significant fraction (at least 1 - targetfixingrate)
3114 * of variables for which the incumbent solution value does not lie within the global bounds anymore. This is a nonsuccess
3115 * for the neighborhood (additional fixings are not possible), which is okay because the incumbent solution is
3116 * significantly outdated
3117 */
3118 if( *nfixings == ntargetfixings )
3119 *result = SCIP_SUCCESS;
3120
3121 /* free the buffer array */
3122 SCIPfreeBufferArray(scip, &varscpy);
3123
3124 return SCIP_OKAY;
3125 }
3126
3127 /** add local branching constraint */
3128 static
addLocalBranchingConstraint(SCIP * sourcescip,SCIP * targetscip,SCIP_VAR ** subvars,int distance,SCIP_Bool * success,int * naddedconss)3129 SCIP_RETCODE addLocalBranchingConstraint(
3130 SCIP* sourcescip, /**< source SCIP data structure */
3131 SCIP* targetscip, /**< target SCIP data structure */
3132 SCIP_VAR** subvars, /**< array of sub SCIP variables in same order as source SCIP variables */
3133 int distance, /**< right hand side of the local branching constraint */
3134 SCIP_Bool* success, /**< pointer to store of a local branching constraint has been successfully added */
3135 int* naddedconss /**< pointer to increase the number of added constraints */
3136 )
3137 {
3138 int nbinvars;
3139 int i;
3140 SCIP_SOL* referencesol;
3141 SCIP_CONS* localbranchcons;
3142 SCIP_VAR** vars;
3143 SCIP_Real* consvals;
3144 SCIP_Real rhs;
3145
3146 assert(sourcescip != NULL);
3147 assert(*success == FALSE);
3148
3149 nbinvars = SCIPgetNBinVars(sourcescip);
3150 vars = SCIPgetVars(sourcescip);
3151
3152 if( nbinvars <= 3 )
3153 return SCIP_OKAY;
3154
3155 referencesol = SCIPgetBestSol(sourcescip);
3156 if( referencesol == NULL )
3157 return SCIP_OKAY;
3158
3159 rhs = (SCIP_Real)distance;
3160 rhs = MAX(rhs, 2.0);
3161
3162 SCIP_CALL( SCIPallocBufferArray(sourcescip, &consvals, nbinvars) );
3163
3164 /* loop over binary variables and fill the local branching constraint */
3165 for( i = 0; i < nbinvars; ++i )
3166 {
3167 /* skip variables that are not present in sub-SCIP */
3168 if( subvars[i] == NULL )
3169 continue;
3170
3171 if( SCIPisEQ(sourcescip, SCIPgetSolVal(sourcescip, referencesol, vars[i]), 0.0) )
3172 consvals[i] = 1.0;
3173 else
3174 {
3175 consvals[i] = -1.0;
3176 rhs -= 1.0;
3177 }
3178 }
3179
3180 /* create the local branching constraint in the target scip */
3181 SCIP_CALL( SCIPcreateConsBasicLinear(targetscip, &localbranchcons, "localbranch", nbinvars, subvars, consvals, -SCIPinfinity(sourcescip), rhs) );
3182 SCIP_CALL( SCIPaddCons(targetscip, localbranchcons) );
3183 SCIP_CALL( SCIPreleaseCons(targetscip, &localbranchcons) );
3184
3185 *naddedconss = 1;
3186 *success = TRUE;
3187
3188 SCIPfreeBufferArray(sourcescip, &consvals);
3189
3190 return SCIP_OKAY;
3191 }
3192
3193 /** callback for local branching subproblem changes */
3194 static
DECL_CHANGESUBSCIP(changeSubscipLocalbranching)3195 DECL_CHANGESUBSCIP(changeSubscipLocalbranching)
3196 { /*lint --e{715}*/
3197
3198 SCIP_CALL( addLocalBranchingConstraint(sourcescip, targetscip, subvars, (int)(0.2 * SCIPgetNBinVars(sourcescip)), success, naddedconss) );
3199
3200 return SCIP_OKAY;
3201 }
3202
3203 /** callback for proximity subproblem changes */
3204 static
DECL_CHANGESUBSCIP(changeSubscipProximity)3205 DECL_CHANGESUBSCIP(changeSubscipProximity)
3206 { /*lint --e{715}*/
3207 SCIP_SOL* referencesol;
3208 SCIP_VAR** vars;
3209 int nbinvars;
3210 int nintvars;
3211 int nvars;
3212 int i;
3213
3214 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, &nvars, &nbinvars, &nintvars, NULL, NULL) );
3215
3216 if( nbinvars == 0 )
3217 return SCIP_OKAY;
3218
3219 referencesol = SCIPgetBestSol(sourcescip);
3220 if( referencesol == NULL )
3221 return SCIP_OKAY;
3222
3223 /* loop over binary variables, set objective coefficients based on reference solution in a local branching fashion */
3224 for( i = 0; i < nbinvars; ++i )
3225 {
3226 SCIP_Real newobj;
3227
3228 /* skip variables not present in sub-SCIP */
3229 if( subvars[i] == NULL )
3230 continue;
3231
3232 if( SCIPgetSolVal(sourcescip, referencesol, vars[i]) < 0.5 )
3233 newobj = -1.0;
3234 else
3235 newobj = 1.0;
3236 SCIP_CALL( SCIPchgVarObj(targetscip, subvars[i], newobj) );
3237 }
3238
3239 /* loop over the remaining variables and change their objective coefficients to 0 */
3240 for( ; i < nvars; ++i )
3241 {
3242 /* skip variables not present in sub-SCIP */
3243 if( subvars[i] == NULL )
3244 continue;
3245
3246 SCIP_CALL( SCIPchgVarObj(targetscip, subvars[i], 0.0) );
3247 }
3248
3249 *nchgobjs = nvars;
3250 *success = TRUE;
3251
3252 return SCIP_OKAY;
3253 }
3254
3255 /** callback for zeroobjective subproblem changes */
3256 static
DECL_CHANGESUBSCIP(changeSubscipZeroobjective)3257 DECL_CHANGESUBSCIP(changeSubscipZeroobjective)
3258 { /*lint --e{715}*/
3259 SCIP_VAR** vars;
3260 int nvars;
3261 int i;
3262
3263 assert(*success == FALSE);
3264
3265 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, &nvars, NULL, NULL, NULL, NULL) );
3266
3267 /* do not run if no objective variables are present */
3268 if( SCIPgetNObjVars(sourcescip) == 0 )
3269 return SCIP_OKAY;
3270
3271 /* loop over the variables and change their objective coefficients to 0 */
3272 for( i = 0; i < nvars; ++i )
3273 {
3274 /* skip variables not present in sub-SCIP */
3275 if( subvars[i] == NULL )
3276 continue;
3277
3278 SCIP_CALL( SCIPchgVarObj(targetscip, subvars[i], 0.0) );
3279 }
3280
3281 *nchgobjs = nvars;
3282 *success = TRUE;
3283
3284 return SCIP_OKAY;
3285 }
3286
3287 /** compute tightened bounds for integer variables depending on how much the LP and the incumbent solution values differ */
3288 static
computeIntegerVariableBoundsDins(SCIP * scip,SCIP_VAR * var,SCIP_Real * lbptr,SCIP_Real * ubptr)3289 void computeIntegerVariableBoundsDins(
3290 SCIP* scip, /**< SCIP data structure of the original problem */
3291 SCIP_VAR* var, /**< the variable for which bounds should be computed */
3292 SCIP_Real* lbptr, /**< pointer to store the lower bound in the DINS sub-SCIP */
3293 SCIP_Real* ubptr /**< pointer to store the upper bound in the DINS sub-SCIP */
3294 )
3295 {
3296 SCIP_Real mipsol;
3297 SCIP_Real lpsol;
3298
3299 SCIP_Real lbglobal;
3300 SCIP_Real ubglobal;
3301 SCIP_SOL* bestsol;
3302
3303 /* get the bounds for each variable */
3304 lbglobal = SCIPvarGetLbGlobal(var);
3305 ubglobal = SCIPvarGetUbGlobal(var);
3306
3307 assert(SCIPvarGetType(var) == SCIP_VARTYPE_INTEGER);
3308 /* get the current LP solution for each variable */
3309 lpsol = SCIPvarGetLPSol(var);
3310
3311 /* get the current MIP solution for each variable */
3312 bestsol = SCIPgetBestSol(scip);
3313 mipsol = SCIPgetSolVal(scip, bestsol, var);
3314
3315 /* if the solution values differ by 0.5 or more, the variable is rebounded, otherwise it is just copied */
3316 if( REALABS(lpsol - mipsol) >= 0.5 )
3317 {
3318 SCIP_Real range;
3319
3320 *lbptr = lbglobal;
3321 *ubptr = ubglobal;
3322
3323 /* create an equally sized range around lpsol for general integers: bounds are lpsol +- (mipsol-lpsol) */
3324 range = 2 * lpsol - mipsol;
3325
3326 if( mipsol >= lpsol )
3327 {
3328 range = SCIPfeasCeil(scip, range);
3329 *lbptr = MAX(*lbptr, range);
3330
3331 /* when the bound new upper bound is equal to the current MIP solution, we set both bounds to the integral bound (without eps) */
3332 if( SCIPisFeasEQ(scip, mipsol, *lbptr) )
3333 *ubptr = *lbptr;
3334 else
3335 *ubptr = mipsol;
3336 }
3337 else
3338 {
3339 range = SCIPfeasFloor(scip, range);
3340 *ubptr = MIN(*ubptr, range);
3341
3342 /* when the bound new upper bound is equal to the current MIP solution, we set both bounds to the integral bound (without eps) */
3343 if( SCIPisFeasEQ(scip, mipsol, *ubptr) )
3344 *lbptr = *ubptr;
3345 else
3346 *lbptr = mipsol;
3347 }
3348
3349 /* the global domain of variables might have been reduced since incumbent was found: adjust lb and ub accordingly */
3350 *lbptr = MAX(*lbptr, lbglobal);
3351 *ubptr = MIN(*ubptr, ubglobal);
3352 }
3353 else
3354 {
3355 /* the global domain of variables might have been reduced since incumbent was found: adjust it accordingly */
3356 *lbptr = MAX(mipsol, lbglobal);
3357 *ubptr = MIN(mipsol, ubglobal);
3358 }
3359 }
3360
3361 /** callback to collect variable fixings of DINS */
3362 static
DECL_VARFIXINGS(varFixingsDins)3363 DECL_VARFIXINGS(varFixingsDins)
3364 {
3365 DATA_DINS* data;
3366 SCIP_SOL* rootlpsol;
3367 SCIP_SOL** sols;
3368 int nsols;
3369 int nmipsols;
3370 int nbinvars;
3371 int nintvars;
3372 SCIP_VAR** vars;
3373 int v;
3374
3375 data = neighborhood->data.dins;
3376 assert(data != NULL);
3377 nmipsols = SCIPgetNSols(scip);
3378 nmipsols = MIN(nmipsols, data->npoolsols);
3379
3380 *result = SCIP_DELAYED;
3381
3382 if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
3383 return SCIP_OKAY;
3384
3385 *result = SCIP_DIDNOTRUN;
3386
3387 if( nmipsols == 0 )
3388 return SCIP_OKAY;
3389
3390 SCIP_CALL( SCIPgetVarsData(scip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
3391
3392 if( nbinvars + nintvars == 0 )
3393 return SCIP_OKAY;
3394
3395 SCIP_CALL( SCIPcreateSol(scip, &rootlpsol, NULL) );
3396
3397 /* save root solution LP values in solution */
3398 for( v = 0; v < nbinvars + nintvars; ++v )
3399 {
3400 SCIP_CALL( SCIPsetSolVal(scip, rootlpsol, vars[v], SCIPvarGetRootSol(vars[v])) );
3401 }
3402
3403 /* add the node and the root LP solution */
3404 nsols = nmipsols + 2;
3405
3406 SCIP_CALL( SCIPallocBufferArray(scip, &sols, nsols) );
3407 sols[0] = NULL; /* node LP solution */
3408 sols[1] = rootlpsol;
3409
3410 /* copy the remaining MIP solutions after the LP solutions */
3411 BMScopyMemoryArray(&sols[2], SCIPgetSols(scip), nmipsols); /*lint !e866*/
3412
3413 /* 1. Binary variables are fixed if their values agree in all the solutions */
3414 if( nbinvars > 0 )
3415 {
3416 SCIP_CALL( fixMatchingSolutionValues(scip, sols, nsols, vars, nbinvars, varbuf, valbuf, nfixings) );
3417 }
3418
3419 /* 2. Integer variables are fixed if they have a very low distance between the incumbent and the root LP solution */
3420 for( v = nbinvars; v < nintvars; ++v )
3421 {
3422 SCIP_Real lb;
3423 SCIP_Real ub;
3424 computeIntegerVariableBoundsDins(scip, vars[v], &lb, &ub);
3425
3426 if( ub - lb < 0.5 )
3427 {
3428 assert(SCIPisFeasIntegral(scip, lb));
3429 tryAdd2variableBuffer(scip, vars[v], lb, varbuf, valbuf, nfixings, TRUE);
3430 }
3431 }
3432
3433 *result = SCIP_SUCCESS;
3434
3435 SCIPfreeBufferArray(scip, &sols);
3436
3437 SCIP_CALL( SCIPfreeSol(scip, &rootlpsol) );
3438
3439 return SCIP_OKAY;
3440 }
3441
3442 /** callback for DINS subproblem changes */
3443 static
DECL_CHANGESUBSCIP(changeSubscipDins)3444 DECL_CHANGESUBSCIP(changeSubscipDins)
3445 { /*lint --e{715}*/
3446 SCIP_VAR** vars;
3447 int nintvars;
3448 int nbinvars;
3449 int v;
3450
3451 SCIP_CALL( SCIPgetVarsData(sourcescip, &vars, NULL, &nbinvars, &nintvars, NULL, NULL) );
3452
3453 /* 1. loop over integer variables and tighten the bounds */
3454 for( v = nbinvars; v < nintvars; ++v )
3455 {
3456 SCIP_Real lb;
3457 SCIP_Real ub;
3458
3459 /* skip variables not present in sub-SCIP */
3460 if( subvars[v] == NULL )
3461 continue;
3462
3463 computeIntegerVariableBoundsDins(sourcescip, vars[v], &lb, &ub);
3464
3465 SCIP_CALL( SCIPchgVarLbGlobal(targetscip, subvars[v], lb) );
3466 SCIP_CALL( SCIPchgVarUbGlobal(targetscip, subvars[v], ub) );
3467 ++(*ndomchgs);
3468 }
3469
3470 /* 2. add local branching constraint for binary variables */
3471 SCIP_CALL( addLocalBranchingConstraint(sourcescip, targetscip, subvars, (int)(0.1 * SCIPgetNBinVars(sourcescip)), success, naddedconss) );
3472
3473 *success = TRUE;
3474
3475 return SCIP_OKAY;
3476 }
3477
3478 /** deinitialization callback for DINS before SCIP is freed */
3479 static
DECL_NHFREE(nhFreeDins)3480 DECL_NHFREE(nhFreeDins)
3481 {
3482 assert(neighborhood->data.dins != NULL);
3483
3484 SCIPfreeBlockMemory(scip, &neighborhood->data.dins);
3485
3486 return SCIP_OKAY;
3487 }
3488
3489 /** deinitialization callback for trustregion before SCIP is freed */
3490 static
DECL_NHFREE(nhFreeTrustregion)3491 DECL_NHFREE(nhFreeTrustregion)
3492 {
3493 assert(neighborhood->data.trustregion != NULL);
3494
3495 SCIPfreeBlockMemory(scip, &neighborhood->data.trustregion);
3496
3497 return SCIP_OKAY;
3498 }
3499
3500 /** add trust region neighborhood constraint and auxiliary objective variable */
3501 static
DECL_CHANGESUBSCIP(changeSubscipTrustregion)3502 DECL_CHANGESUBSCIP(changeSubscipTrustregion)
3503 { /*lint --e{715}*/
3504 DATA_TRUSTREGION* data;
3505
3506 data = neighborhood->data.trustregion;
3507
3508 /* adding the neighborhood constraint for the trust region heuristic */
3509 SCIP_CALL( SCIPaddTrustregionNeighborhoodConstraint(sourcescip, targetscip, subvars, data->violpenalty) );
3510
3511 /* incrementing the change in objective since an additional variable is added to the objective to penalize the
3512 * violation of the trust region.
3513 */
3514 ++(*nchgobjs);
3515
3516 return SCIP_OKAY;
3517 }
3518
3519 /** callback that returns the incumbent solution as a reference point */
3520 static
DECL_NHREFSOL(nhRefsolIncumbent)3521 DECL_NHREFSOL(nhRefsolIncumbent)
3522 { /*lint --e{715}*/
3523 assert(scip != NULL);
3524
3525 if( SCIPgetBestSol(scip) != NULL )
3526 {
3527 *result = SCIP_SUCCESS;
3528 *solptr = SCIPgetBestSol(scip);
3529 }
3530 else
3531 {
3532 *result = SCIP_DIDNOTFIND;
3533 }
3534
3535 return SCIP_OKAY;
3536 }
3537
3538
3539 /** callback function that deactivates a neighborhood on problems with no discrete variables */
3540 static
DECL_NHDEACTIVATE(nhDeactivateDiscreteVars)3541 DECL_NHDEACTIVATE(nhDeactivateDiscreteVars)
3542 { /*lint --e{715}*/
3543 assert(scip != NULL);
3544 assert(deactivate != NULL);
3545
3546 /* deactivate if no discrete variables are present */
3547 *deactivate = (SCIPgetNBinVars(scip) + SCIPgetNIntVars(scip) == 0);
3548
3549 return SCIP_OKAY;
3550 }
3551
3552 /** callback function that deactivates a neighborhood on problems with no binary variables */
3553 static
DECL_NHDEACTIVATE(nhDeactivateBinVars)3554 DECL_NHDEACTIVATE(nhDeactivateBinVars)
3555 { /*lint --e{715}*/
3556 assert(scip != NULL);
3557 assert(deactivate != NULL);
3558
3559 /* deactivate if no discrete variables are present */
3560 *deactivate = (SCIPgetNBinVars(scip) == 0);
3561
3562 return SCIP_OKAY;
3563 }
3564
3565 /** callback function that deactivates a neighborhood on problems with no objective variables */
3566 static
DECL_NHDEACTIVATE(nhDeactivateObjVars)3567 DECL_NHDEACTIVATE(nhDeactivateObjVars)
3568 { /*lint --e{715}*/
3569 assert(scip != NULL);
3570 assert(deactivate != NULL);
3571
3572 /* deactivate if no discrete variables are present */
3573 *deactivate = (SCIPgetNObjVars(scip) == 0);
3574
3575 return SCIP_OKAY;
3576 }
3577
3578
3579 /** include all neighborhoods */
3580 static
includeNeighborhoods(SCIP * scip,SCIP_HEURDATA * heurdata)3581 SCIP_RETCODE includeNeighborhoods(
3582 SCIP* scip, /**< SCIP data structure */
3583 SCIP_HEURDATA* heurdata /**< heuristic data of the ALNS heuristic */
3584 )
3585 {
3586 NH* rens;
3587 NH* rins;
3588 NH* mutation;
3589 NH* localbranching;
3590 NH* crossover;
3591 NH* proximity;
3592 NH* zeroobjective;
3593 NH* dins;
3594 NH* trustregion;
3595
3596 heurdata->nneighborhoods = 0;
3597
3598 /* include the RENS neighborhood */
3599 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &rens, "rens",
3600 DEFAULT_MINFIXINGRATE_RENS, DEFAULT_MAXFIXINGRATE_RENS, DEFAULT_ACTIVE_RENS, DEFAULT_PRIORITY_RENS,
3601 varFixingsRens, changeSubscipRens, NULL, NULL, NULL, NULL, nhDeactivateDiscreteVars) );
3602
3603 /* include the RINS neighborhood */
3604 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &rins, "rins",
3605 DEFAULT_MINFIXINGRATE_RINS, DEFAULT_MAXFIXINGRATE_RINS, DEFAULT_ACTIVE_RINS, DEFAULT_PRIORITY_RINS,
3606 varFixingsRins, NULL, NULL, NULL, NULL, nhRefsolIncumbent, nhDeactivateDiscreteVars) );
3607
3608 /* include the mutation neighborhood */
3609 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &mutation, "mutation",
3610 DEFAULT_MINFIXINGRATE_MUTATION, DEFAULT_MAXFIXINGRATE_MUTATION, DEFAULT_ACTIVE_MUTATION, DEFAULT_PRIORITY_MUTATION,
3611 varFixingsMutation, NULL, nhInitMutation, nhExitMutation, NULL, nhRefsolIncumbent, nhDeactivateDiscreteVars) );
3612
3613 /* include the local branching neighborhood */
3614 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &localbranching, "localbranching",
3615 DEFAULT_MINFIXINGRATE_LOCALBRANCHING, DEFAULT_MAXFIXINGRATE_LOCALBRANCHING, DEFAULT_ACTIVE_LOCALBRANCHING, DEFAULT_PRIORITY_LOCALBRANCHING,
3616 NULL, changeSubscipLocalbranching, NULL, NULL, NULL, nhRefsolIncumbent, nhDeactivateBinVars) );
3617
3618 /* include the crossover neighborhood */
3619 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &crossover, "crossover",
3620 DEFAULT_MINFIXINGRATE_CROSSOVER, DEFAULT_MAXFIXINGRATE_CROSSOVER, DEFAULT_ACTIVE_CROSSOVER, DEFAULT_PRIORITY_CROSSOVER,
3621 varFixingsCrossover, NULL,
3622 nhInitCrossover, nhExitCrossover, nhFreeCrossover, nhRefsolCrossover, nhDeactivateDiscreteVars) );
3623
3624 /* allocate data for crossover to include the parameter */
3625 SCIP_CALL( SCIPallocBlockMemory(scip, &crossover->data.crossover) );
3626 crossover->data.crossover->rng = NULL;
3627
3628 /* add crossover neighborhood parameters */
3629 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/alns/crossover/nsols", "the number of solutions that crossover should combine",
3630 &crossover->data.crossover->nsols, TRUE, DEFAULT_NSOLS_CROSSOVER, 2, 10, NULL, NULL) );
3631
3632 /* include the Proximity neighborhood */
3633 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &proximity, "proximity",
3634 DEFAULT_MINFIXINGRATE_PROXIMITY, DEFAULT_MAXFIXINGRATE_PROXIMITY, DEFAULT_ACTIVE_PROXIMITY, DEFAULT_PRIORITY_PROXIMITY,
3635 NULL, changeSubscipProximity, NULL, NULL, NULL, nhRefsolIncumbent, nhDeactivateBinVars) );
3636
3637 /* include the Zeroobjective neighborhood */
3638 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &zeroobjective, "zeroobjective",
3639 DEFAULT_MINFIXINGRATE_ZEROOBJECTIVE, DEFAULT_MAXFIXINGRATE_ZEROOBJECTIVE, DEFAULT_ACTIVE_ZEROOBJECTIVE, DEFAULT_PRIORITY_ZEROOBJECTIVE,
3640 NULL, changeSubscipZeroobjective, NULL, NULL, NULL, nhRefsolIncumbent, nhDeactivateObjVars) );
3641
3642 /* include the DINS neighborhood */
3643 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &dins, "dins",
3644 DEFAULT_MINFIXINGRATE_DINS, DEFAULT_MAXFIXINGRATE_DINS, DEFAULT_ACTIVE_DINS, DEFAULT_PRIORITY_DINS,
3645 varFixingsDins, changeSubscipDins, NULL, NULL, nhFreeDins, nhRefsolIncumbent, nhDeactivateBinVars) );
3646
3647 /* allocate data for DINS to include the parameter */
3648 SCIP_CALL( SCIPallocBlockMemory(scip, &dins->data.dins) );
3649
3650 /* add DINS neighborhood parameters */
3651 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/alns/dins/npoolsols",
3652 "number of pool solutions where binary solution values must agree",
3653 &dins->data.dins->npoolsols, TRUE, DEFAULT_NPOOLSOLS_DINS, 1, 100, NULL, NULL) );
3654
3655 /* include the trustregion neighborhood */
3656 SCIP_CALL( alnsIncludeNeighborhood(scip, heurdata, &trustregion, "trustregion",
3657 DEFAULT_MINFIXINGRATE_TRUSTREGION, DEFAULT_MAXFIXINGRATE_TRUSTREGION, DEFAULT_ACTIVE_TRUSTREGION, DEFAULT_PRIORITY_TRUSTREGION,
3658 NULL, changeSubscipTrustregion, NULL, NULL, nhFreeTrustregion, nhRefsolIncumbent, nhDeactivateBinVars) );
3659
3660 /* allocate data for trustregion to include the parameter */
3661 SCIP_CALL( SCIPallocBlockMemory(scip, &trustregion->data.trustregion) );
3662
3663 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/trustregion/violpenalty",
3664 "the penalty for each change in the binary variables from the candidate solution",
3665 &trustregion->data.trustregion->violpenalty, FALSE, DEFAULT_VIOLPENALTY_TRUSTREGION, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3666
3667 return SCIP_OKAY;
3668 }
3669
3670 /** initialization method of primal heuristic (called after problem was transformed) */
3671 static
SCIP_DECL_HEURINIT(heurInitAlns)3672 SCIP_DECL_HEURINIT(heurInitAlns)
3673 { /*lint --e{715}*/
3674 SCIP_HEURDATA* heurdata;
3675 int i;
3676
3677 assert(scip != NULL);
3678 assert(heur != NULL);
3679
3680 heurdata = SCIPheurGetData(heur);
3681 assert(heurdata != NULL);
3682
3683 /* reactivate all neighborhoods if a new problem is read in */
3684 heurdata->nactiveneighborhoods = heurdata->nneighborhoods;
3685
3686 /* initialize neighborhoods for new problem */
3687 for( i = 0; i < heurdata->nneighborhoods; ++i )
3688 {
3689 NH* neighborhood = heurdata->neighborhoods[i];
3690
3691 SCIP_CALL( neighborhoodInit(scip, neighborhood) );
3692
3693 SCIP_CALL( resetFixingRate(scip, &neighborhood->fixingrate) );
3694
3695 SCIP_CALL( neighborhoodStatsReset(scip, &neighborhood->stats) );
3696 }
3697
3698 /* open reward file for reading */
3699 if( strncasecmp(heurdata->rewardfilename, DEFAULT_REWARDFILENAME, strlen(DEFAULT_REWARDFILENAME)) != 0 )
3700 {
3701 heurdata->rewardfile = fopen(heurdata->rewardfilename, "w");
3702
3703 if( heurdata->rewardfile == NULL )
3704 {
3705 SCIPerrorMessage("Error: Could not open reward file <%s>\n", heurdata->rewardfilename);
3706 return SCIP_FILECREATEERROR;
3707 }
3708
3709 SCIPdebugMsg(scip, "Writing reward information to <%s>\n", heurdata->rewardfilename);
3710 }
3711 else
3712 heurdata->rewardfile = NULL;
3713
3714 return SCIP_OKAY;
3715 }
3716
3717
3718 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
3719 static
SCIP_DECL_HEURINITSOL(heurInitsolAlns)3720 SCIP_DECL_HEURINITSOL(heurInitsolAlns)
3721 { /*lint --e{715}*/
3722 SCIP_HEURDATA* heurdata;
3723 int i;
3724 SCIP_Real* priorities;
3725 unsigned int initseed;
3726
3727 assert(scip != NULL);
3728 assert(heur != NULL);
3729
3730 heurdata = SCIPheurGetData(heur);
3731 assert(heurdata != NULL);
3732 heurdata->nactiveneighborhoods = heurdata->nneighborhoods;
3733
3734 SCIP_CALL( SCIPallocBufferArray(scip, &priorities, heurdata->nactiveneighborhoods) );
3735
3736 /* init neighborhoods for new problem by resetting their statistics and fixing rate */
3737 for( i = heurdata->nneighborhoods - 1; i >= 0; --i )
3738 {
3739 NH* neighborhood = heurdata->neighborhoods[i];
3740 SCIP_Bool deactivate;
3741
3742 SCIP_CALL( neighborhood->nhdeactivate(scip, &deactivate) );
3743
3744 /* disable inactive neighborhoods */
3745 if( deactivate || ! neighborhood->active )
3746 {
3747 if( heurdata->nactiveneighborhoods - 1 > i )
3748 {
3749 assert(heurdata->neighborhoods[heurdata->nactiveneighborhoods - 1]->active);
3750 SCIPswapPointers((void **)&heurdata->neighborhoods[i], (void **)&heurdata->neighborhoods[heurdata->nactiveneighborhoods - 1]);
3751 }
3752 heurdata->nactiveneighborhoods--;
3753 }
3754 }
3755
3756 /* collect neighborhood priorities */
3757 for( i = 0; i < heurdata->nactiveneighborhoods; ++i )
3758 priorities[i] = heurdata->neighborhoods[i]->priority;
3759
3760 initseed = (unsigned int)(heurdata->seed + SCIPgetNVars(scip));
3761
3762 /* active neighborhoods might change between init calls, reset functionality must take this into account */
3763 if( heurdata->bandit != NULL && SCIPbanditGetNActions(heurdata->bandit) != heurdata->nactiveneighborhoods )
3764 {
3765 SCIP_CALL( SCIPfreeBandit(scip, &heurdata->bandit) );
3766
3767 heurdata->bandit = NULL;
3768 }
3769
3770 if( heurdata->nactiveneighborhoods > 0 )
3771 { /* create or reset bandit algorithm */
3772 if( heurdata->bandit == NULL )
3773 {
3774 SCIP_CALL( createBandit(scip, heurdata, priorities, initseed) );
3775
3776 resetMinimumImprovement(heurdata);
3777 resetTargetNodeLimit(heurdata);
3778 }
3779 else if( heurdata->resetweights )
3780 {
3781 SCIP_CALL( SCIPresetBandit(scip, heurdata->bandit, priorities, initseed) );
3782
3783 resetMinimumImprovement(heurdata);
3784 resetTargetNodeLimit(heurdata);
3785 }
3786 }
3787
3788 heurdata->usednodes = 0;
3789 heurdata->ninitneighborhoods = heurdata->nactiveneighborhoods;
3790 resetCurrentNeighborhood(heurdata);
3791
3792 SCIPfreeBufferArray(scip, &priorities);
3793
3794 return SCIP_OKAY;
3795 }
3796
3797
3798 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
3799 static
SCIP_DECL_HEUREXIT(heurExitAlns)3800 SCIP_DECL_HEUREXIT(heurExitAlns)
3801 { /*lint --e{715}*/
3802 SCIP_HEURDATA* heurdata;
3803 int i;
3804
3805 assert(scip != NULL);
3806 assert(heur != NULL);
3807
3808 heurdata = SCIPheurGetData(heur);
3809 assert(heurdata != NULL);
3810
3811 /* free neighborhood specific data */
3812 for( i = 0; i < heurdata->nneighborhoods; ++i )
3813 {
3814 NH* neighborhood = heurdata->neighborhoods[i];
3815
3816 SCIP_CALL( neighborhoodExit(scip, neighborhood) );
3817 }
3818
3819 if( heurdata->rewardfile != NULL )
3820 {
3821 fclose(heurdata->rewardfile);
3822 heurdata->rewardfile = NULL;
3823 }
3824
3825 return SCIP_OKAY;
3826 }
3827
3828 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
3829 static
SCIP_DECL_HEURFREE(heurFreeAlns)3830 SCIP_DECL_HEURFREE(heurFreeAlns)
3831 { /*lint --e{715}*/
3832 SCIP_HEURDATA* heurdata;
3833 int i;
3834
3835 assert(scip != NULL);
3836 assert(heur != NULL);
3837
3838 heurdata = SCIPheurGetData(heur);
3839 assert(heurdata != NULL);
3840
3841 /* bandits are only initialized if a problem has been read */
3842 if( heurdata->bandit != NULL )
3843 {
3844 SCIP_CALL( SCIPfreeBandit(scip, &heurdata->bandit) );
3845 }
3846
3847 /* free neighborhoods */
3848 for( i = 0; i < heurdata->nneighborhoods; ++i )
3849 {
3850 SCIP_CALL( alnsFreeNeighborhood(scip, &(heurdata->neighborhoods[i])) );
3851 }
3852
3853 SCIPfreeBlockMemoryArray(scip, &heurdata->neighborhoods, NNEIGHBORHOODS);
3854
3855 SCIPfreeBlockMemory(scip, &heurdata);
3856
3857 return SCIP_OKAY;
3858 }
3859
3860 /** output method of statistics table to output file stream 'file' */
3861 static
SCIP_DECL_TABLEOUTPUT(tableOutputNeighborhood)3862 SCIP_DECL_TABLEOUTPUT(tableOutputNeighborhood)
3863 { /*lint --e{715}*/
3864 SCIP_HEURDATA* heurdata;
3865
3866 assert(SCIPfindHeur(scip, HEUR_NAME) != NULL);
3867 heurdata = SCIPheurGetData(SCIPfindHeur(scip, HEUR_NAME));
3868 assert(heurdata != NULL);
3869
3870 printNeighborhoodStatistics(scip, heurdata, file);
3871
3872 return SCIP_OKAY;
3873 }
3874
3875 /*
3876 * primal heuristic specific interface methods
3877 */
3878
3879 /** creates the alns primal heuristic and includes it in SCIP */
SCIPincludeHeurAlns(SCIP * scip)3880 SCIP_RETCODE SCIPincludeHeurAlns(
3881 SCIP* scip /**< SCIP data structure */
3882 )
3883 {
3884 SCIP_HEURDATA* heurdata;
3885 SCIP_HEUR* heur;
3886
3887 /* create alns primal heuristic data */
3888 heurdata = NULL;
3889 heur = NULL;
3890
3891 SCIP_CALL( SCIPallocBlockMemory(scip, &heurdata) );
3892 BMSclearMemory(heurdata);
3893
3894 /* TODO make this a user parameter? */
3895 heurdata->lplimfac = LPLIMFAC;
3896
3897 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &heurdata->neighborhoods, NNEIGHBORHOODS) );
3898
3899 /* include primal heuristic */
3900 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
3901 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
3902 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecAlns, heurdata) );
3903
3904 assert(heur != NULL);
3905
3906 /* include all neighborhoods */
3907 SCIP_CALL( includeNeighborhoods(scip, heurdata) );
3908
3909 /* set non fundamental callbacks via setter functions */
3910 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyAlns) );
3911 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeAlns) );
3912 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitAlns) );
3913 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolAlns) );
3914 SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitAlns) );
3915
3916 /* add alns primal heuristic parameters */
3917 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/maxnodes",
3918 "maximum number of nodes to regard in the subproblem",
3919 &heurdata->maxnodes, TRUE,DEFAULT_MAXNODES, 0LL, SCIP_LONGINT_MAX, NULL, NULL) );
3920
3921 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/nodesofs",
3922 "offset added to the nodes budget",
3923 &heurdata->nodesoffset, FALSE, DEFAULT_NODESOFFSET, 0LL, SCIP_LONGINT_MAX, NULL, NULL) );
3924
3925 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/minnodes",
3926 "minimum number of nodes required to start a sub-SCIP",
3927 &heurdata->minnodes, TRUE, DEFAULT_MINNODES, 0LL, SCIP_LONGINT_MAX, NULL, NULL) );
3928
3929 SCIP_CALL( SCIPaddLongintParam(scip, "heuristics/" HEUR_NAME "/waitingnodes",
3930 "number of nodes since last incumbent solution that the heuristic should wait",
3931 &heurdata->waitingnodes, TRUE, DEFAULT_WAITINGNODES, 0LL, SCIP_LONGINT_MAX, NULL, NULL) );
3932
3933 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/nodesquot",
3934 "fraction of nodes compared to the main SCIP for budget computation",
3935 &heurdata->nodesquot, FALSE, DEFAULT_NODESQUOT, 0.0, 1.0, NULL, NULL) );
3936
3937 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/startminimprove",
3938 "initial factor by which ALNS should at least improve the incumbent",
3939 &heurdata->startminimprove, TRUE, DEFAULT_STARTMINIMPROVE, 0.0, 1.0, NULL, NULL) );
3940
3941 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprovelow",
3942 "lower threshold for the minimal improvement over the incumbent",
3943 &heurdata->minimprovelow, TRUE, DEFAULT_MINIMPROVELOW, 0.0, 1.0, NULL, NULL) );
3944
3945 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/minimprovehigh",
3946 "upper bound for the minimal improvement over the incumbent",
3947 &heurdata->minimprovehigh, TRUE, DEFAULT_MINIMPROVEHIGH, 0.0, 1.0, NULL, NULL) );
3948
3949 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/nsolslim",
3950 "limit on the number of improving solutions in a sub-SCIP call",
3951 &heurdata->nsolslim, FALSE, DEFAULT_NSOLSLIM, -1, INT_MAX, NULL, NULL) );
3952
3953 SCIP_CALL( SCIPaddCharParam(scip, "heuristics/" HEUR_NAME "/banditalgo",
3954 "the bandit algorithm: (u)pper confidence bounds, (e)xp.3, epsilon (g)reedy",
3955 &heurdata->banditalgo, TRUE, DEFAULT_BANDITALGO, "ueg", NULL, NULL) );
3956
3957 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/gamma",
3958 "weight between uniform (gamma ~ 1) and weight driven (gamma ~ 0) probability distribution for exp3",
3959 &heurdata->exp3_gamma, TRUE, DEFAULT_GAMMA, 0.0, 1.0, NULL, NULL) );
3960
3961 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/beta",
3962 "reward offset between 0 and 1 at every observation for Exp.3",
3963 &heurdata->exp3_beta, TRUE, DEFAULT_BETA, 0.0, 1.0, NULL, NULL) );
3964
3965 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/alpha",
3966 "parameter to increase the confidence width in UCB",
3967 &heurdata->ucb_alpha, TRUE, DEFAULT_ALPHA, 0.0, 100.0, NULL, NULL) );
3968
3969 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/usedistances",
3970 "distances from fixed variables be used for variable prioritization",
3971 &heurdata->usedistances, TRUE, DEFAULT_USEDISTANCES, NULL, NULL) );
3972
3973 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/useredcost",
3974 "should reduced cost scores be used for variable prioritization?",
3975 &heurdata->useredcost, TRUE, DEFAULT_USEREDCOST, NULL, NULL) );
3976
3977 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/domorefixings",
3978 "should the ALNS heuristic do more fixings by itself based on variable prioritization "
3979 "until the target fixing rate is reached?",
3980 &heurdata->domorefixings, TRUE, DEFAULT_DOMOREFIXINGS, NULL, NULL) );
3981
3982 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/adjustfixingrate",
3983 "should the heuristic adjust the target fixing rate based on the success?",
3984 &heurdata->adjustfixingrate, TRUE, DEFAULT_ADJUSTFIXINGRATE, NULL, NULL) );
3985
3986 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/usesubscipheurs",
3987 "should the heuristic activate other sub-SCIP heuristics during its search?",
3988 &heurdata->usesubscipheurs, TRUE, DEFAULT_USESUBSCIPHEURS, NULL, NULL) );
3989
3990 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/rewardcontrol",
3991 "reward control to increase the weight of the simple solution indicator and decrease the weight of the closed gap reward",
3992 &heurdata->rewardcontrol, TRUE, DEFAULT_REWARDCONTROL, 0.0, 1.0, NULL, NULL) );
3993
3994 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/targetnodefactor",
3995 "factor by which target node number is eventually increased",
3996 &heurdata->targetnodefactor, TRUE, DEFAULT_TARGETNODEFACTOR, 1.0, 1e+5, NULL, NULL) );
3997
3998 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/" HEUR_NAME "/seed",
3999 "initial random seed for bandit algorithms and random decisions by neighborhoods",
4000 &heurdata->seed, FALSE, DEFAULT_SEED, 0, INT_MAX, NULL, NULL) );
4001
4002 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/adjustminimprove",
4003 "should the factor by which the minimum improvement is bound be dynamically updated?",
4004 &heurdata->adjustminimprove, TRUE, DEFAULT_ADJUSTMINIMPROVE, NULL, NULL) );
4005
4006 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/adjusttargetnodes",
4007 "should the target nodes be dynamically adjusted?",
4008 &heurdata->adjusttargetnodes, TRUE, DEFAULT_ADJUSTTARGETNODES, NULL, NULL) );
4009
4010 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/eps",
4011 "increase exploration in epsilon-greedy bandit algorithm",
4012 &heurdata->epsgreedy_eps, TRUE, DEFAULT_EPS, 0.0, 1.0, NULL, NULL) );
4013
4014 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/rewardbaseline",
4015 "the reward baseline to separate successful and failed calls",
4016 &heurdata->rewardbaseline, TRUE, DEFAULT_REWARDBASELINE, 0.0, 0.99, NULL, NULL) );
4017
4018 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/resetweights",
4019 "should the bandit algorithms be reset when a new problem is read?",
4020 &heurdata->resetweights, TRUE, DEFAULT_RESETWEIGHTS, NULL, NULL) );
4021
4022 SCIP_CALL( SCIPaddStringParam(scip, "heuristics/" HEUR_NAME "/rewardfilename", "file name to store all rewards and the selection of the bandit",
4023 &heurdata->rewardfilename, TRUE, DEFAULT_REWARDFILENAME, NULL, NULL) );
4024
4025 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/subsciprandseeds",
4026 "should random seeds of sub-SCIPs be altered to increase diversification?",
4027 &heurdata->subsciprandseeds, TRUE, DEFAULT_SUBSCIPRANDSEEDS, NULL, NULL) );
4028
4029 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/scalebyeffort",
4030 "should the reward be scaled by the effort?",
4031 &heurdata->scalebyeffort, TRUE, DEFAULT_SCALEBYEFFORT, NULL, NULL) );
4032
4033 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/copycuts",
4034 "should cutting planes be copied to the sub-SCIP?",
4035 &heurdata->copycuts, TRUE, DEFAULT_COPYCUTS, NULL, NULL) );
4036
4037 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/fixtol",
4038 "tolerance by which the fixing rate may be missed without generic fixing",
4039 &heurdata->fixtol, TRUE, DEFAULT_FIXTOL, 0.0, 1.0, NULL, NULL) );
4040
4041 SCIP_CALL( SCIPaddRealParam(scip, "heuristics/" HEUR_NAME "/unfixtol",
4042 "tolerance by which the fixing rate may be exceeded without generic unfixing",
4043 &heurdata->unfixtol, TRUE, DEFAULT_UNFIXTOL, 0.0, 1.0, NULL, NULL) );
4044
4045 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/uselocalredcost",
4046 "should local reduced costs be used for generic (un)fixing?",
4047 &heurdata->uselocalredcost, TRUE, DEFAULT_USELOCALREDCOST, NULL, NULL) );
4048
4049 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/" HEUR_NAME "/usepscost",
4050 "should pseudo cost scores be used for variable priorization?",
4051 &heurdata->usepscost, TRUE, DEFAULT_USEPSCOST, NULL, NULL) );
4052
4053 assert(SCIPfindTable(scip, TABLE_NAME_NEIGHBORHOOD) == NULL);
4054 SCIP_CALL( SCIPincludeTable(scip, TABLE_NAME_NEIGHBORHOOD, TABLE_DESC_NEIGHBORHOOD, TRUE,
4055 NULL, NULL, NULL, NULL, NULL, NULL, tableOutputNeighborhood,
4056 NULL, TABLE_POSITION_NEIGHBORHOOD, TABLE_EARLIEST_STAGE_NEIGHBORHOOD) );
4057
4058 return SCIP_OKAY;
4059 }
4060