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_ascendprune.c
17 * @brief reduction-based primal heuristic for Steiner problems
18 * @author Daniel Rehfeldt
19 *
20 * This file implements a reduction and dual-cost based heuristic for Steiner problems. See
21 * "SCIP-Jack - A solver for STP and variants with parallelization extensions" (2016) by
22 * Gamrath, Koch, Maher, Rehfeldt and Shinano
23 *
24 * A list of all interface methods can be found in heur_ascendprune.h
25 *
26 */
27
28 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
29
30 #include <assert.h>
31 #include <string.h>
32 #include <stdio.h>
33 #include "scip/scip.h"
34 #include "scip/scipdefplugins.h"
35 #include "scip/cons_linear.h"
36 #include "heur_ascendprune.h"
37 #include "heur_local.h"
38 #include "heur_prune.h"
39 #include "grph.h"
40 #include "heur_tm.h"
41 #include "cons_stp.h"
42 #include "scip/pub_misc.h"
43 #include "probdata_stp.h"
44
45 #define HEUR_NAME "ascendprune"
46 #define HEUR_DESC "Dual-cost reduction heuristic for Steiner problems"
47 #define HEUR_DISPCHAR 'A'
48 #define HEUR_PRIORITY 2
49 #define HEUR_FREQ -1
50 #define HEUR_FREQOFS 0
51 #define HEUR_MAXDEPTH -1
52 #define HEUR_TIMING (SCIP_HEURTIMING_DURINGLPLOOP | SCIP_HEURTIMING_AFTERLPLOOP | SCIP_HEURTIMING_AFTERNODE)
53 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
54
55 #define DEFAULT_MAXFREQPRUNE FALSE /**< executions of the heuristic at maximum frequency? */
56 #define ASCENPRUNE_MINLPIMPROVE 0.05 /**< minimum percentual improvement of dual bound (wrt to gap) mandatory to execute heuristic */
57
58 #ifdef WITH_UG
59 int getUgRank(void);
60 #endif
61
62 /*
63 * Data structures
64 */
65
66 /** primal heuristic data */
67 struct SCIP_HeurData
68 {
69 SCIP_Real lastdualbound; /**< dual bound after the previous run */
70 int bestsolindex; /**< best solution during the previous run */
71 int nfailures; /**< number of failures since last successful call */
72 SCIP_Bool maxfreq; /**< should the heuristic be called at maximum frequency? */
73 };
74
75
76 /*
77 * Local methods
78 */
79
80 /* put your local methods here, and declare them static */
81
82
83 /*
84 * Callback methods of primal heuristic
85 */
86
87
88 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
89 static
SCIP_DECL_HEURCOPY(heurCopyAscendPrune)90 SCIP_DECL_HEURCOPY(heurCopyAscendPrune)
91 { /*lint --e{715}*/
92 assert(scip != NULL);
93 assert(heur != NULL);
94 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
95
96 /* call inclusion method of primal heuristic */
97 SCIP_CALL( SCIPStpIncludeHeurAscendPrune(scip) );
98
99 return SCIP_OKAY;
100 }
101
102 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
103 static
SCIP_DECL_HEURFREE(heurFreeAscendPrune)104 SCIP_DECL_HEURFREE(heurFreeAscendPrune)
105 { /*lint --e{715}*/
106 SCIP_HEURDATA* heurdata;
107
108 assert(heur != NULL);
109 assert(scip != NULL);
110
111 /* get heuristic data */
112 heurdata = SCIPheurGetData(heur);
113 assert(heurdata != NULL);
114
115 /* free heuristic data */
116 SCIPfreeMemory(scip, &heurdata);
117 SCIPheurSetData(heur, NULL);
118
119 return SCIP_OKAY;
120 }
121
122
123 /** initialization method of primal heuristic (called after problem was transformed) */
124
125 static
SCIP_DECL_HEURINIT(heurInitAscendPrune)126 SCIP_DECL_HEURINIT(heurInitAscendPrune)
127 { /*lint --e{715}*/
128 SCIP_HEURDATA* heurdata;
129
130 assert(heur != NULL);
131 assert(scip != NULL);
132
133 /* get heuristic's data */
134 heurdata = SCIPheurGetData(heur);
135
136 assert(heurdata != NULL);
137
138 /* initialize data */
139 heurdata->nfailures = 0;
140 heurdata->bestsolindex = -1;
141 heurdata->lastdualbound = 0.0;
142
143 return SCIP_OKAY;
144 }
145
146 /** execution method of primal heuristic */
147 static
SCIP_DECL_HEUREXEC(heurExecAscendPrune)148 SCIP_DECL_HEUREXEC(heurExecAscendPrune)
149 { /*lint --e{715}*/
150 SCIP_HEURDATA* heurdata;
151 SCIP_PROBDATA* probdata;
152 SCIP_VAR** vars;
153 SCIP_SOL* bestsol; /* best known solution */
154 GRAPH* graph;
155 SCIP_Real dualbound;
156 SCIP_Real gap;
157 SCIP_Real* redcosts;
158 SCIP_Bool success;
159 int e;
160 int nnodes;
161 int nedges;
162 int* edgearrint;
163 int* nodearrint;
164 STP_Bool* nodearrchar;
165
166 assert(heur != NULL);
167 assert(scip != NULL);
168 assert(result != NULL);
169 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
170
171 /* get heuristic data */
172 heurdata = SCIPheurGetData(heur);
173 assert(heurdata != NULL);
174
175 /* get problem data */
176 probdata = SCIPgetProbData(scip);
177 assert(probdata != NULL);
178
179 /* get graph */
180 graph = SCIPprobdataGetGraph(probdata);
181 assert(graph != NULL);
182
183 vars = SCIPprobdataGetVars(scip);
184 assert(vars != NULL);
185
186 *result = SCIP_DIDNOTRUN;
187
188 /* todo: delete this file and move to slack-prune */
189
190 nedges = graph->edges;
191 nnodes = graph->knots;
192 success = FALSE;
193
194 /* get best current solution */
195 bestsol = SCIPgetBestSol(scip);
196
197 /* no solution available? */
198 if( bestsol == NULL )
199 return SCIP_OKAY;
200
201 /* get dual bound */
202 dualbound = SCIPgetDualbound(scip);
203
204 /* no new best solution available? */
205 if( heurdata->bestsolindex == SCIPsolGetIndex(SCIPgetBestSol(scip)) && !(heurdata->maxfreq) )
206 {
207 /* current optimality gap */
208 gap = SCIPgetSolOrigObj(scip, bestsol) - dualbound;
209
210 if( SCIPisLT(scip, dualbound - heurdata->lastdualbound, gap * ASCENPRUNE_MINLPIMPROVE ) )
211 return SCIP_OKAY;
212 }
213
214 heurdata->lastdualbound = dualbound;
215
216 /* allocate memory for ascent and prune */
217 SCIP_CALL( SCIPallocBufferArray(scip, &redcosts, nedges) );
218 SCIP_CALL( SCIPallocBufferArray(scip, &edgearrint, nedges) );
219 SCIP_CALL( SCIPallocBufferArray(scip, &nodearrint, nnodes ) );
220 SCIP_CALL( SCIPallocBufferArray(scip, &nodearrchar, nnodes) );
221
222 for( e = 0; e < nedges; e++ )
223 {
224 assert(SCIPvarIsBinary(vars[e]));
225
226 /* variable is already fixed, we must not trust the reduced cost */
227 if( SCIPvarGetLbLocal(vars[e]) + 0.5 > SCIPvarGetUbLocal(vars[e]) )
228 {
229 if( SCIPvarGetLbLocal(vars[e]) > 0.5 )
230 redcosts[e] = 0.0;
231 else
232 {
233 assert(SCIPvarGetUbLocal(vars[e]) < 0.5);
234 redcosts[e] = FARAWAY;
235 }
236 }
237 else
238 {
239 if( SCIPisFeasZero(scip, SCIPgetSolVal(scip, NULL, vars[e])) )
240 {
241 assert(!SCIPisDualfeasNegative(scip, SCIPgetVarRedcost(scip, vars[e])));
242 redcosts[e] = SCIPgetVarRedcost(scip, vars[e]);
243 }
244 else
245 {
246 assert(!SCIPisDualfeasPositive(scip, SCIPgetVarRedcost(scip, vars[e])));
247 assert(SCIPisFeasEQ(scip, SCIPgetSolVal(scip, NULL, vars[e]), 1.0) || SCIPisDualfeasZero(scip, SCIPgetVarRedcost(scip, vars[e])));
248 redcosts[e] = 0.0;
249 }
250 }
251
252 if( SCIPisLT(scip, redcosts[e], 0.0) )
253 redcosts[e] = 0.0;
254
255 assert(SCIPisGE(scip, redcosts[e], 0.0));
256 assert(!SCIPisEQ(scip, redcosts[e], SCIP_INVALID));
257 }
258
259 /* perform ascent and prune */
260 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, heur, graph, redcosts, edgearrint, nodearrint, graph->source, nodearrchar, &success, TRUE) );
261
262 if( success )
263 {
264 heurdata->nfailures = 0;
265 *result = SCIP_FOUNDSOL;
266 }
267 else
268 {
269 heurdata->nfailures++;
270 }
271
272 heurdata->bestsolindex = SCIPsolGetIndex(SCIPgetBestSol(scip));
273
274 /* free memory */
275 SCIPfreeBufferArray(scip, &nodearrchar);
276 SCIPfreeBufferArray(scip, &nodearrint);
277 SCIPfreeBufferArray(scip, &edgearrint);
278 SCIPfreeBufferArray(scip, &redcosts);
279
280 return SCIP_OKAY;
281 }
282
283
284 /*
285 * primal heuristic specific interface methods
286 */
287
288
289 /** ascent and prune */
SCIPStpHeurAscendPruneRun(SCIP * scip,SCIP_HEUR * heur,const GRAPH * g,const SCIP_Real * redcosts,int * edgearrint,int * nodearrint,int root,STP_Bool * nodearrchar,SCIP_Bool * solfound,SCIP_Bool addsol)290 SCIP_RETCODE SCIPStpHeurAscendPruneRun(
291 SCIP* scip, /**< SCIP data structure */
292 SCIP_HEUR* heur, /**< heuristic data structure or NULL */
293 const GRAPH* g, /**< the graph */
294 const SCIP_Real* redcosts, /**< the reduced costs */
295 int* edgearrint, /**< int edges array to store solution */
296 int* nodearrint, /**< int vertices array for internal computations */
297 int root, /**< the root (used for dual ascent) */
298 STP_Bool* nodearrchar, /**< STP_Bool vertices array for internal computations */
299 SCIP_Bool* solfound, /**< has a solution been found? */
300 SCIP_Bool addsol /**< should the solution be added to SCIP by this method? */
301 )
302 {
303 GRAPH* newgraph;
304 SCIP_Real* nval;
305 int* const mark = g->mark;
306 int* const newedges = edgearrint;
307 int* const nodechild = nodearrint;
308 int* edgeancestor;
309 const int nnodes = g->knots;
310 const int nedges = g->edges;
311 const int probtype = g->stp_type;
312 int nnewnodes = 0;
313 int nnewedges = 0;
314 const SCIP_Bool pcmw = (probtype == STP_PCSPG || probtype == STP_MWCSP || probtype == STP_RPCSPG || probtype == STP_RMWCSP);
315 SCIP_Bool success;
316
317 assert(g != NULL);
318 assert(scip != NULL);
319 assert(redcosts != NULL);
320 assert(edgearrint != NULL);
321 assert(nodearrint != NULL);
322 assert(nodearrchar != NULL);
323
324 if( root < 0 )
325 root = g->source;
326
327 assert(Is_term(g->term[root]));
328 assert(graph_valid(g));
329
330 if( addsol )
331 {
332 const int nvars = SCIPprobdataGetNVars(scip);
333 SCIP_CALL( SCIPallocBufferArray(scip, &nval, nvars) );
334 }
335 else
336 {
337 nval = NULL;
338 }
339
340 /* DFS to identify 0-redcost subgraph */
341 {
342 int* const queue = nodearrint;
343 STP_Bool* const scanned = nodearrchar;
344 int qsize;
345
346 /*
347 * construct new graph corresponding to zero cost paths from the root to all terminals
348 */
349
350 BMSclearMemoryArray(mark, nnodes);
351 BMSclearMemoryArray(scanned, nnodes);
352
353 qsize = 0;
354 mark[root] = TRUE;
355 queue[qsize++] = root;
356 nnewnodes++;
357
358 /* DFS */
359 while( qsize )
360 {
361 const int k = queue[--qsize];
362 scanned[k] = TRUE;
363
364 /* traverse outgoing arcs */
365 for( int a = g->outbeg[k]; a != EAT_LAST; a = g->oeat[a] )
366 {
367 const int head = g->head[a];
368
369 if( SCIPisZero(scip, redcosts[a]) )
370 {
371 if( pcmw && k == root && Is_term(g->term[head]) )
372 continue;
373
374 /* vertex not labeled yet? */
375 if( !mark[head] )
376 {
377 mark[head] = TRUE;
378 nnewnodes++;
379 queue[qsize++] = head;
380 }
381
382 if( (!scanned[head] || !SCIPisZero(scip, redcosts[flipedge(a)])) && k != root )
383 {
384 assert(g->tail[a] != root);
385 assert(g->head[a] != root);
386
387 newedges[nnewedges++] = a;
388 }
389 }
390 }
391 }
392
393 #ifndef NDEBUG
394 for( int k = 0; k < nnewedges && pcmw; k++ )
395 {
396 const int e = newedges[k];
397 assert(!(g->tail[e] == root && Is_pterm(g->term[g->head[e]])));
398 assert(!(g->head[e] == root && Is_pterm(g->term[g->tail[e]])));
399 }
400 #endif
401
402 for( int a = g->outbeg[root]; a != EAT_LAST; a = g->oeat[a] )
403 {
404 const int head = g->head[a];
405 if( mark[head] )
406 newedges[nnewedges++] = a;
407 }
408 }
409
410 SCIP_CALL( SCIPallocBufferArray(scip, &edgeancestor, 2 * nnewedges) );
411
412 /* initialize new graph */
413 SCIP_CALL( graph_init(scip, &newgraph, nnewnodes, 2 * nnewedges, 1) );
414
415 if( probtype == STP_RSMT || probtype == STP_OARSMT || probtype == STP_GSTP )
416 newgraph->stp_type = STP_SPG;
417 else
418 newgraph->stp_type = probtype;
419
420 if( pcmw )
421 SCIP_CALL( graph_pc_init(scip, newgraph, nnewnodes, nnewnodes) );
422
423 for( int k = 0; k < nnodes; k++ )
424 {
425 if( mark[k] )
426 {
427 if( pcmw )
428 {
429 if( (!Is_term(g->term[k])) )
430 newgraph->prize[newgraph->knots] = g->prize[k];
431 else
432 newgraph->prize[newgraph->knots] = 0.0;
433 }
434
435 nodechild[k] = newgraph->knots;
436 graph_knot_add(newgraph, g->term[k]);
437 }
438 }
439
440 if( pcmw )
441 {
442 newgraph->norgmodelknots = nnewnodes;
443 newgraph->extended = TRUE;
444 }
445
446 assert(nnewnodes == newgraph->knots);
447
448 /* set root of new graph */
449 newgraph->source = nodechild[root];
450 assert(newgraph->source >= 0);
451
452 if( g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP )
453 newgraph->prize[newgraph->source] = FARAWAY;
454
455 /* add edges to new graph */
456 for( int a = 0; a < nnewedges; a++ )
457 {
458 int i;
459 const int e = newedges[a];
460 const int tail = nodechild[g->tail[e]];
461 const int head = nodechild[g->head[e]];
462
463 assert(tail >= 0);
464 assert(head >= 0);
465
466 for( i = newgraph->outbeg[tail]; i != EAT_LAST; i = newgraph->oeat[i] )
467 if( newgraph->head[i] == head )
468 break;
469
470 if( i == EAT_LAST )
471 {
472 edgeancestor[newgraph->edges] = e;
473 edgeancestor[newgraph->edges + 1] = flipedge(e);
474
475 if( pcmw )
476 graph_pc_updateTerm2edge(newgraph, g, tail, head, g->tail[e], g->head[e]);
477
478 graph_edge_add(scip, newgraph, tail, head, g->cost[e], g->cost[flipedge(e)]);
479 }
480 }
481
482 nnewedges = newgraph->edges;
483 newgraph->norgmodeledges = nnewedges;
484
485 assert(!pcmw || -1 == newgraph->term2edge[newgraph->source]);
486
487 /* initialize ancestors of new graph edges */
488 SCIP_CALL( graph_init_history(scip, newgraph) );
489
490 /* initialize shortest path algorithm */
491 SCIP_CALL( graph_path_init(scip, newgraph) );
492
493 SCIP_CALL( level0(scip, newgraph) );
494
495 #ifdef DEBUG_ASCENDPRUNE
496 for( int k = 0; k < nnodes && !pcmw; k++ )
497 {
498 if( Is_term(g->term[k]) )
499 {
500 const int i = nodechild[k];
501 if( i < 0 )
502 {
503 printf("k %d root %d \n", k, root);
504 printf("FAIL in AP \n\n\n");
505 return SCIP_ERROR;
506 }
507
508 if( newgraph->grad[i] == 0 && newgraph->knots > 1 )
509 {
510 printf("FAIL GRAD \n\n\n");
511 return SCIP_ERROR;
512
513 }
514 }
515 }
516 #endif
517 assert(graph_valid(newgraph));
518
519 /* get solution on new graph by PRUNE heuristic */
520 SCIP_CALL( SCIPStpHeurPruneRun(scip, NULL, newgraph, newedges, &success, FALSE, TRUE) );
521
522 #ifdef DEBUG_ASCENDPRUNE
523 for( int k = 0; k < newgraph->knots; ++k )
524 {
525 if( Is_term(newgraph->term[k]) && newgraph->grad[k] == 0 && k != newgraph->source )
526 {
527 printf("after i %d r %d \n", k, root);
528 return SCIP_ERROR;
529 }
530 }
531
532 if( !graph_sol_valid(scip, newgraph, newedges) )
533 {
534 printf("not valid %d \n", 0);
535 return SCIP_ERROR;
536 }
537 #endif
538 if( !success )
539 {
540 SCIPdebugMessage("failed to build tree in ascend-prune (by prune) \n");
541 goto TERMINATE;
542 }
543
544 assert(success && graph_sol_valid(scip, newgraph, newedges));
545
546 SCIPdebugMessage("obj after prune %f \n", graph_sol_getObj(newgraph->cost, newedges, 0.0, newgraph->edges));
547
548 SCIP_CALL( SCIPStpHeurLocalRun(scip, newgraph, newgraph->cost, newedges) );
549
550 SCIPdebugMessage("obj after local %f \n", graph_sol_getObj(newgraph->cost, newedges, 0.0, newgraph->edges));
551
552 assert(graph_sol_valid(scip, newgraph, newedges));
553 graph_path_exit(scip, newgraph);
554
555
556 /*
557 * prune solution (in the original graph)
558 */
559
560 BMSclearMemoryArray(nodearrchar, nnodes);
561
562 for( int e = 0; e < nnewedges; e++ )
563 if( newedges[e] == CONNECT )
564 {
565 const int eorg = edgeancestor[e];
566 nodearrchar[g->tail[eorg]] = TRUE;
567 nodearrchar[g->head[eorg]] = TRUE;
568 }
569
570 for( int e = 0; e < nedges; e++ )
571 newedges[e] = UNKNOWN;
572
573 if( pcmw )
574 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, g, g->cost, newedges, nodearrchar) );
575 else
576 SCIP_CALL( SCIPStpHeurTMPrune(scip, g, g->cost, 0, newedges, nodearrchar) );
577
578 assert(graph_sol_valid(scip, g, newedges));
579
580 if( addsol )
581 {
582 assert(nval != NULL);
583 for( int e = 0; e < nedges; e++ )
584 {
585 if( newedges[e] == CONNECT )
586 nval[e] = 1.0;
587 else
588 nval[e] = 0.0;
589 }
590 }
591
592 success = graph_sol_valid(scip, g, newedges);
593
594 if( success && addsol )
595 {
596 /* add solution */
597 SCIP_SOL* sol = NULL;
598 SCIP_CALL( SCIPprobdataAddNewSol(scip, nval, sol, heur, &success) );
599 SCIPdebugMessage("Ascend-and-prune added solution \n");
600 }
601
602 *solfound = success;
603
604 TERMINATE:
605
606 for( int k = 0; k < nnodes; k++ )
607 mark[k] = (g->grad[k] > 0);
608
609 /* free memory */
610 graph_free(scip, &newgraph, TRUE);
611 SCIPfreeBufferArray(scip, &edgeancestor);
612 SCIPfreeBufferArrayNull(scip, &nval);
613
614 return SCIP_OKAY;
615 }
616
617
618 /** creates the prune primal heuristic and includes it in SCIP */
SCIPStpIncludeHeurAscendPrune(SCIP * scip)619 SCIP_RETCODE SCIPStpIncludeHeurAscendPrune(
620 SCIP* scip /**< SCIP data structure */
621 )
622 {
623 SCIP_HEURDATA* heurdata;
624 SCIP_HEUR* heur;
625
626 /* create prune primal heuristic data */
627 SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
628
629 /* include primal heuristic */
630 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
631 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
632 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecAscendPrune, heurdata) );
633
634 assert(heur != NULL);
635
636 /* set non fundamental callbacks via setter functions */
637 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyAscendPrune) );
638 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeAscendPrune) );
639 SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitAscendPrune) );
640
641 /* add ascend and prune primal heuristic parameters */
642 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/"HEUR_NAME"/maxfreq",
643 "should the heuristic be executed at maximum frequeny?",
644 &heurdata->maxfreq, FALSE, DEFAULT_MAXFREQPRUNE, NULL, NULL) );
645
646 return SCIP_OKAY;
647 }
648