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_local.c
17 * @brief Improvement heuristic for Steiner problems
18 * @author Daniel Rehfeldt
19 *
20 * This file implements several local heuristics, including vertex insertion, key-path exchange and key-vertex elimination,
21 * ("Fast Local Search for Steiner Trees in Graphs" by Uchoa and Werneck). Other heuristics are for PCSTP and MWCSP.
22 *
23 * A list of all interface methods can be found in heur_local.h.
24 *
25 */
26
27 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
28
29 #include <assert.h>
30 #include <string.h>
31
32 #include "heur_local.h"
33 #include "heur_tm.h"
34 #include "probdata_stp.h"
35
36
37 /* @note if heuristic is running in root node timing is changed there to (SCIP_HEURTIMING_DURINGLPLOOP |
38 * SCIP_HEURTIMING_BEFORENODE), see SCIP_DECL_HEURINITSOL callback
39 */
40
41 #define HEUR_NAME "local"
42 #define HEUR_DESC "improvement heuristic for STP"
43 #define HEUR_DISPCHAR '-'
44 #define HEUR_PRIORITY 100
45 #define HEUR_FREQ 1
46 #define HEUR_FREQOFS 0
47 #define HEUR_MAXDEPTH -1
48 #define HEUR_TIMING (SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_DURINGLPLOOP | SCIP_HEURTIMING_AFTERLPLOOP | SCIP_HEURTIMING_AFTERNODE)
49
50 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
51
52 #define DEFAULT_DURINGROOT TRUE
53 #define DEFAULT_MAXFREQLOC FALSE
54 #define DEFAULT_MAXNBESTSOLS 25
55 #define DEFAULT_NBESTSOLS 10
56 #define DEFAULT_MINNBESTSOLS 6
57 #define LOCAL_MAXRESTARTS 5
58
59 #define GREEDY_MAXRESTARTS 3 /**< Max number of restarts for greedy PC/MW heuristic if improving solution has been found. */
60 #define GREEDY_EXTENSIONS_MW 6 /**< Number of extensions for greedy MW heuristic. MUST BE HIGHER THAN GREEDY_EXTENSIONS */
61 #define GREEDY_EXTENSIONS 5 /**< Number of extensions for greedy PC heuristic. */
62
63
64 /*
65 * Data structures
66 */
67
68 /** primal heuristic data */
69 struct SCIP_HeurData
70 {
71 int nfails; /**< number of fails */
72 int maxnsols; /**< maximal number of best solutions to improve */
73 int nbestsols; /**< number of best solutions to improve */
74 int* lastsolindices; /**< indices of a number of best solutions already tried */
75 SCIP_Bool maxfreq; /**< should the heuristic be called with maximum frequency? */
76 SCIP_Bool duringroot; /**< should the heuristic be called during the root node? */
77 };
78
79
80 /*
81 * Local methods
82 */
83
84 /** recursive methode for a DFS ordering of graph 'g' */
85 static
dfsorder(const GRAPH * graph,const int * edges,const int * node,int * counter,int * dfst)86 void dfsorder(
87 const GRAPH* graph,
88 const int* edges,
89 const int* node,
90 int* counter,
91 int* dfst
92 )
93 {
94 int oedge = graph->outbeg[*node];
95
96 while( oedge >= 0 )
97 {
98 if( edges[oedge] >= 0 )
99 {
100 dfsorder(graph, edges, &(graph->head[oedge]), counter, dfst);
101 }
102 oedge = graph->oeat[oedge];
103 }
104
105 dfst[*counter] = *node;
106 (*counter)++;
107 }
108
109
110 /** computes lowest common ancestors for all pairs {vbase(v), vbase(u)} such that {u,w} is a boundary edge,
111 * first call should be with u := root */
112 static
lca(SCIP * scip,const GRAPH * graph,int u,UF * uf,STP_Bool * nodesmark,int * steineredges,IDX ** lcalists,PHNODE ** boundpaths,int * heapsize,int * vbase)113 SCIP_RETCODE lca(
114 SCIP* scip,
115 const GRAPH* graph,
116 int u,
117 UF* uf,
118 STP_Bool* nodesmark,
119 int* steineredges,
120 IDX** lcalists,
121 PHNODE** boundpaths,
122 int* heapsize,
123 int* vbase
124 )
125 {
126 int* uboundpaths; /* boundary-paths (each one represented by its boundary edge) having node 'u' as an endpoint */
127 int ancestor;
128 int v;
129 int i;
130 int oedge; /* outgoing edge */
131 IDX* curr;
132 uf->parent[u] = u;
133
134 for( oedge = graph->outbeg[u]; oedge != EAT_LAST; oedge = graph->oeat[oedge] )
135 {
136 v = graph->head[oedge];
137 if( steineredges[oedge] == 0 )
138 {
139 SCIP_CALL( lca(scip, graph, v, uf, nodesmark, steineredges, lcalists, boundpaths, heapsize, vbase) );
140 SCIPStpunionfindUnion(uf, u, v, FALSE);
141 uf->parent[SCIPStpunionfindFind(uf, u)] = u;
142 }
143 }
144 nodesmark[u] = TRUE;
145
146 /* iterate through all boundary-paths having one endpoint in the voronoi region of node 'u' */
147 SCIP_CALL( SCIPpairheapBuffarr(scip, boundpaths[u], heapsize[u], &uboundpaths) );
148 for( i = 0; i < heapsize[u]; i++ )
149 {
150 oedge = uboundpaths[i];
151 v = vbase[graph->head[oedge]];
152 if( nodesmark[v] )
153 {
154 ancestor = uf->parent[SCIPStpunionfindFind(uf, v)];
155
156 /* if the ancestor of 'u' and 'v' is one of the two, the boundary-edge is already in boundpaths[u] */
157 if( ancestor != u && ancestor != v)
158 {
159 SCIP_CALL( SCIPallocBlockMemory(scip, &curr) );
160 curr->index = oedge;
161 curr->parent = lcalists[ancestor];
162 lcalists[ancestor] = curr;
163 }
164 }
165 }
166
167 /* free the boundary-paths array */
168 SCIPfreeBufferArray(scip, &uboundpaths);
169
170 return SCIP_OKAY;
171 }
172
173 /** checks whether node is crucial, i.e. a terminal or a vertex with degree at least 3 (w.r.t. the steinertree) */
174 static
nodeIsCrucial(const GRAPH * graph,int * steineredges,int node)175 STP_Bool nodeIsCrucial(
176 const GRAPH* graph,
177 int* steineredges,
178 int node
179 )
180 {
181 assert(graph != NULL);
182 assert(steineredges != NULL);
183
184 if( graph->term[node] == -1 )
185 {
186 int counter = 0;
187 int e = graph->outbeg[node];
188 while( e >= 0 )
189 {
190 /* check if the adjacent node is in the ST */
191 if( steineredges[e] > -1 || steineredges[flipedge(e)] > -1 )
192 {
193 counter++;
194 }
195 e = graph->oeat[e];
196 }
197
198 if( counter < 3 )
199 {
200 return FALSE;
201 }
202 }
203
204 return TRUE;
205 }
206
207 /** perform local heuristics on a given Steiner tree todo delete cost parameter */
SCIPStpHeurLocalRun(SCIP * scip,GRAPH * graph,const SCIP_Real * cost,int * best_result)208 SCIP_RETCODE SCIPStpHeurLocalRun(
209 SCIP* scip, /**< SCIP data structure */
210 GRAPH* graph, /**< graph data structure */
211 const SCIP_Real* cost, /**< arc cost array todo delete this parameter and use graph->cost */
212 int* best_result /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
213 )
214 {
215 NODE* nodes;
216 PATH* vnoi;
217 int* vbase;
218 int* const graphmark = graph->mark;
219 int e;
220 int i;
221 int k;
222 const int root = graph->source;
223 const int nnodes = graph->knots;
224 const int nedges = graph->edges;
225 const int probtype = graph->stp_type;
226 int newnverts;
227 STP_Bool* steinertree;
228 const STP_Bool pc = ((probtype == STP_PCSPG) || (probtype == STP_RPCSPG));
229 const STP_Bool mw = (probtype == STP_MWCSP);
230 const STP_Bool mwpc = (pc || mw);
231 #ifndef NDEBUG
232 const SCIP_Real initialobj = graph_sol_getObj(graph->cost, best_result, 0.0, nedges);
233 #endif
234
235 assert(graph != NULL);
236 assert(cost != NULL);
237 assert(best_result != NULL);
238 assert(graph_valid(graph));
239
240 newnverts = 0;
241
242 if( graph->grad[root] == 0 || graph->terms == 1 )
243 return SCIP_OKAY;
244
245 /* for PC variants test whether solution is trivial */
246 if( mwpc )
247 {
248 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
249 if( !Is_term(graph->term[graph->head[e]]) && best_result[e] )
250 break;
251
252 if( e == EAT_LAST )
253 {
254 SCIPdebugMessage("Local heuristic: return trivial \n");
255 return SCIP_OKAY;
256 }
257 }
258
259 SCIP_CALL( SCIPallocBufferArray(scip, &vnoi, nnodes) );
260 SCIP_CALL( SCIPallocBufferArray(scip, &vbase, nnodes) );
261 SCIP_CALL( SCIPallocBufferArray(scip, &nodes, nnodes) );
262 SCIP_CALL( SCIPallocBufferArray(scip, &steinertree, nnodes) );
263
264 if( mwpc )
265 {
266 SCIP_Bool dummy;
267 SCIP_CALL( SCIPStpHeurLocalExtendPcMw(scip, graph, cost, vnoi, best_result, steinertree, &dummy) );
268 }
269
270 for( i = 0; i < nnodes; i++ )
271 {
272 steinertree[i] = FALSE;
273 SCIPlinkcuttreeInit(&nodes[i]);
274 }
275
276 /* create a link-cut tree representing the current Steiner tree */
277 for( e = 0; e < nedges; e++ )
278 {
279 assert(graph->head[e] == graph->tail[flipedge(e)]);
280
281 /* if edge e is in the tree, so are its incident vertices */
282 if( best_result[e] == CONNECT )
283 {
284 steinertree[graph->tail[e]] = TRUE;
285 steinertree[graph->head[e]] = TRUE;
286 SCIPlinkcuttreeLink(&nodes[graph->head[e]], &nodes[graph->tail[e]], flipedge(e));
287 }
288 }
289
290 assert( nodes[root].edge == -1 );
291 nodes[root].edge = -1;
292
293 /* VERTEX INSERTION */
294 if( probtype == STP_SPG || probtype == STP_RSMT || probtype == STP_OARSMT || probtype == STP_GSTP || (mwpc) )
295 {
296 int newnode;
297 int* insert;
298 int* adds;
299 int* cuts;
300 int* cuts2;
301 int* stdeg;
302
303 SCIP_CALL( SCIPallocBufferArray(scip, &insert, nnodes) );
304 SCIP_CALL( SCIPallocBufferArray(scip, &adds, nnodes) );
305 SCIP_CALL( SCIPallocBufferArray(scip, &cuts, nnodes) );
306
307 if( mw )
308 {
309 SCIP_CALL( SCIPallocBufferArray(scip, &cuts2, nnodes) );
310 SCIP_CALL( SCIPallocBufferArray(scip, &stdeg, nnodes) );
311
312 BMSclearMemoryArray(stdeg, nnodes);
313
314 for( e = 0; e < nedges; e++ )
315 if( best_result[e] == CONNECT )
316 {
317 stdeg[graph->tail[e]]++;
318 stdeg[graph->head[e]]++;
319 }
320 }
321 else
322 {
323 cuts2 = NULL;
324 stdeg = NULL;
325 }
326
327 i = 0;
328 newnode = 0;
329
330 for( ;; )
331 {
332 SCIP_Real diff;
333
334 /* if vertex i is not in the current ST and has at least two adjacent nodes, it might be added */
335 if( !steinertree[i] && graph->grad[i] > 1 && (!mwpc || !Is_term(graph->term[i])) )
336 {
337 NODE* v;
338 int counter;
339 int lastnodeidx;
340 int insertcount = 0;
341
342 /* if an outgoing edge of vertex i points to the current ST, SCIPlinkcuttreeLink the edge to a list */
343 for( int oedge = graph->outbeg[i]; oedge != EAT_LAST; oedge = graph->oeat[oedge])
344 if( steinertree[graph->head[oedge]] && (!mwpc || !Is_term(graph->term[graph->head[oedge]])) )
345 insert[insertcount++] = oedge;
346
347 /* if there are less than two edges connecting node i and the current tree, continue */
348 if( insertcount <= 1 )
349 goto ENDOFLOOP;
350
351 if( mw )
352 SCIPlinkcuttreeInit(&nodes[i]);
353
354 /* the node to insert */
355 v = &nodes[i];
356
357 SCIPlinkcuttreeLink(v, &nodes[graph->head[insert[0]]], insert[0]);
358
359 lastnodeidx = graph->head[insert[0]];
360
361 if( mw )
362 {
363 assert(!SCIPisPositive(scip, graph->prize[i]));
364
365 diff = -1.0;
366 assert(stdeg != NULL);
367 stdeg[lastnodeidx]++;
368 }
369 else
370 diff = graph->cost[v->edge];
371
372 counter = 0;
373
374 /* try to add edges between new vertex and tree */
375 for( k = 1; k < insertcount; k++ )
376 {
377 NODE* firstnode;
378 int firstnodidx;
379 SCIPlinkcuttreeEvert(v);
380
381 /* next vertex in the current Steiner tree adjacent to vertex i resp. v (the one being scrutinized for possible insertion) */
382 firstnodidx = graph->head[insert[k]];
383 firstnode = &nodes[firstnodidx];
384
385 if( mw )
386 {
387 NODE* chainfirst;
388 NODE* chainlast;
389 SCIP_Real minweight;
390
391 assert(stdeg != NULL);
392
393 minweight = SCIPlinkcuttreeFindMinChain(scip, graph->prize, graph->head, stdeg, firstnode, &chainfirst, &chainlast);
394
395 if( SCIPisLT(scip, minweight, graph->prize[i]) )
396 {
397 assert(chainfirst != NULL && chainlast != NULL);
398 for( NODE* mynode = chainfirst; mynode != chainlast; mynode = mynode->parent )
399 {
400 int mynodeidx = graph->head[mynode->edge];
401 steinertree[mynodeidx] = FALSE;
402 stdeg[mynodeidx] = 0;
403 }
404
405 SCIPlinkcuttreeCut(chainfirst);
406 SCIPlinkcuttreeCut(chainlast);
407
408 SCIPlinkcuttreeLink(v, firstnode, insert[k]);
409 stdeg[graph->head[insert[k]]]++;
410
411 diff = graph->prize[i] - minweight;
412 break;
413 }
414 }
415 else
416 {
417 /* if there is an edge with cost greater than that of the current edge... */
418 NODE* max = SCIPlinkcuttreeFindMax(scip, graph->cost, firstnode);
419 if( SCIPisGT(scip, graph->cost[max->edge], graph->cost[insert[k]]) )
420 {
421 diff += graph->cost[insert[k]];
422 diff -= graph->cost[max->edge];
423 cuts[counter] = max->edge;
424 SCIPlinkcuttreeCut(max);
425 SCIPlinkcuttreeLink(v, firstnode, insert[k]);
426 assert(v->edge == insert[k]);
427 adds[counter++] = v->edge;
428 }
429 }
430 }
431
432 if( pc && Is_pterm(graph->term[i]) )
433 diff -= graph->prize[i];
434
435 /* if the new tree is more expensive than the old one, restore the latter */
436 if( mw )
437 {
438 if( SCIPisLT(scip, diff, 0.0) )
439 {
440 assert(stdeg != NULL);
441
442 SCIPlinkcuttreeEvert(v);
443 stdeg[lastnodeidx]--;
444 SCIPlinkcuttreeCut(&nodes[graph->head[insert[0]]]);
445 }
446 else
447 {
448 steinertree[i] = TRUE;
449 newnverts++;
450 }
451 }
452 else
453 {
454 if( !SCIPisNegative(scip, diff) )
455 {
456 SCIPlinkcuttreeEvert(v);
457 for (k = counter - 1; k >= 0; k--)
458 {
459 SCIPlinkcuttreeCut(&nodes[graph->head[adds[k]]]);
460 SCIPlinkcuttreeEvert(&nodes[graph->tail[cuts[k]]]);
461 SCIPlinkcuttreeLink(&nodes[graph->tail[cuts[k]]],
462 &nodes[graph->head[cuts[k]]], cuts[k]);
463 }
464
465 /* finally, cut the edge added first (if it had been cut during the insertion process, it would have been restored above) */
466 SCIPlinkcuttreeEvert(v);
467 SCIPlinkcuttreeCut(&nodes[graph->head[insert[0]]]);
468 }
469 else
470 {
471 SCIPlinkcuttreeEvert(&nodes[root]);
472 adds[counter] = insert[0];
473 newnode = i;
474 steinertree[i] = TRUE;
475 newnverts++;
476 SCIPdebugMessage("ADDED VERTEX \n");
477 }
478 }
479 }
480
481 ENDOFLOOP:
482
483 if( i < nnodes - 1 )
484 i++;
485 else
486 i = 0;
487
488 if( newnode == i )
489 break;
490 }
491
492 /* free buffer memory */
493 if( mw )
494 {
495 SCIPfreeBufferArray(scip, &stdeg);
496 SCIPfreeBufferArray(scip, &cuts2);
497 }
498 SCIPfreeBufferArray(scip, &cuts);
499 SCIPfreeBufferArray(scip, &adds);
500 SCIPfreeBufferArray(scip, &insert);
501
502 for( e = 0; e < nedges; e++ )
503 best_result[e] = UNKNOWN;
504
505 if( newnverts > 0 )
506 {
507 if( mwpc )
508 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, best_result, steinertree) );
509 else
510 SCIP_CALL( SCIPStpHeurTMPrune(scip, graph, graph->cost, 0, best_result, steinertree) );
511
512 for( i = 0; i < nnodes; i++ )
513 SCIPlinkcuttreeInit(&nodes[i]);
514
515 /* create a link-cut tree representing the current Steiner tree */
516 for( e = 0; e < nedges; e++ )
517 {
518 if( best_result[e] == CONNECT )
519 {
520 assert(steinertree[graph->tail[e]]);
521 assert(steinertree[graph->head[e]]);
522 SCIPlinkcuttreeLink(&nodes[graph->head[e]], &nodes[graph->tail[e]], flipedge(e));
523 }
524 }
525 SCIPlinkcuttreeEvert(&nodes[root]);
526 }
527 else
528 {
529 SCIPlinkcuttreeEvert(&nodes[root]);
530 for( i = 0; i < nnodes; i++ )
531 {
532 if( steinertree[i] && nodes[i].edge != -1 )
533 best_result[flipedge(nodes[i].edge)] = 0;
534 }
535 }
536
537 #ifdef SCIP_DEBUG
538 {
539 SCIP_Real obj = 0.0;
540 for( e = 0; e < nedges; e++ )
541 obj += (best_result[e] > -1) ? graph->cost[e] : 0.0;
542 printf("ObjAfterVertexInsertion=%.12e\n", obj);
543
544 }
545 #endif
546 }
547
548 assert(graph_sol_valid(scip, graph, best_result));
549
550 /* Key-Vertex Elimination & Key-Path Exchange */
551 if( !mw )
552 {
553 IDX* blists_curr;
554 IDX** blists_start; /* array [1,..,nnodes],
555 * if node i is in the current ST, blists_start[i] points to a linked list of all nodes having i as their base */
556 PATH* mst; /* minimal spanning tree structure */
557 GRAPH* supergraph;
558 IDX** lvledges_start; /* horizontal edges */
559 IDX* lvledges_curr;
560 PHNODE** boundpaths;
561 UF uf; /* union-find*/
562 SCIP_Real* memdist;
563 SCIP_Real kpcost;
564 SCIP_Real mstcost;
565 SCIP_Real edgecost;
566 SCIP_Real kpathcost;
567 int* state;
568 int* kpedges;
569 int* kpnodes;
570 int* dfstree;
571 int* newedges;
572 int* memvbase;
573 int* heapsize;
574 int* boundedges;
575 int* meminedges;
576 int* supernodes;
577 int* supernodesid;
578 int l;
579 int node;
580 int edge;
581 int count;
582 int nruns;
583 int newedge;
584 int oldedge;
585 int adjnode;
586 int nkpedges;
587 int nstnodes;
588 int nkpnodes;
589 int crucnode; /* current crucial node*/
590 int nresnodes;
591 int kptailnode; /* tail node of the current keypath*/
592 int localmoves = 2;
593 int newpathend = -1;
594 int nsupernodes;
595 int nboundedges;
596 int rootpathstart;
597 STP_Bool* pinned;
598 STP_Bool* scanned;
599 STP_Bool* nodesmark;
600
601 #ifdef SCIP_DEBUG
602 SCIP_Real obj = 0.0;
603 for( e = 0; e < nedges; e++ )
604 obj += (best_result[e] > -1) ? graph->cost[e] : 0.0;
605 printf(" ObjBEFKEYVertexELimination=%.12e\n", obj);
606 #endif
607
608 for( k = 0; k < nnodes; k++ )
609 graphmark[k] = (graph->grad[k] > 0);
610
611 graphmark[root] = TRUE;
612
613 /* allocate memory */
614
615 /* only needed for Key-Path Elimination */
616 SCIP_CALL( SCIPallocBufferArray(scip, &newedges, nedges) );
617 SCIP_CALL( SCIPallocBufferArray(scip, &lvledges_start, nnodes) );
618 SCIP_CALL( SCIPallocBufferArray(scip, &boundedges, nedges) );
619 SCIP_CALL( SCIPallocBufferArray(scip, &supernodesid, nnodes) );
620
621 /* only needed for Key-Path Exchange */
622
623 /* memory needed for both Key-Path Elimination and Exchange */
624 SCIP_CALL( SCIPallocBufferArray(scip, &scanned, nnodes) );
625 SCIP_CALL( SCIPallocBufferArray(scip, &heapsize, nnodes) );
626 SCIP_CALL( SCIPallocBufferArray(scip, &blists_start, nnodes) );
627 SCIP_CALL( SCIPallocBufferArray(scip, &memvbase, nnodes) );
628 SCIP_CALL( SCIPallocBufferArray(scip, &memdist, nnodes) );
629 SCIP_CALL( SCIPallocBufferArray(scip, &meminedges, nnodes) );
630 SCIP_CALL( SCIPallocBufferArray(scip, &boundpaths, nnodes) );
631 SCIP_CALL( SCIPallocBufferArray(scip, &pinned, nnodes) );
632 SCIP_CALL( SCIPallocBufferArray(scip, &dfstree, nnodes) );
633 SCIP_CALL( SCIPallocBufferArray(scip, &nodesmark, nnodes) );
634
635 /* initialize data structures */
636 SCIP_CALL( SCIPStpunionfindInit(scip, &uf, nnodes) );
637
638 for( nruns = 0; nruns < LOCAL_MAXRESTARTS && localmoves > 0; nruns++ )
639 {
640 localmoves = 0;
641
642 BMSclearMemoryArray(blists_start, nnodes);
643
644 /* find a DFS order of the ST nodes */
645 nstnodes = 0;
646 dfsorder(graph, best_result, &(root), &nstnodes, dfstree);
647 assert(root == graph->source);
648
649 /* compute a voronoi diagram with the ST nodes as bases */
650 graph_voronoi(scip, graph, graph->cost, graph->cost, steinertree, vbase, vnoi);
651
652 state = graph->path_state;
653
654 /* initialize data structures */
655 for( k = 0; k < nnodes; k++ )
656 {
657 assert(state[k] == CONNECT || !graphmark[k]);
658
659 pinned[k] = FALSE;
660 scanned[k] = FALSE;
661 nodesmark[k] = FALSE;
662
663 /* initialize pairing heaps */
664 heapsize[k] = 0;
665 boundpaths[k] = NULL;
666
667 lvledges_start[k] = NULL;
668
669 if( !graphmark[k] )
670 continue;
671
672 /* link all nodes to their (respective) voronoi base */
673 SCIP_CALL( SCIPallocBlockMemory(scip, &blists_curr) );
674 blists_curr->index = k;
675 blists_curr->parent = blists_start[vbase[k]];
676 blists_start[vbase[k]] = blists_curr;
677 }
678
679 SCIP_CALL( SCIPallocBufferArray(scip, &supernodes, nstnodes) );
680 SCIP_CALL( SCIPallocBufferArray(scip, &kpnodes, nstnodes) );
681 SCIP_CALL( SCIPallocBufferArray(scip, &kpedges, nstnodes) );
682
683 if( mwpc )
684 {
685 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
686 {
687 k = graph->head[e];
688 if( Is_term(graph->term[k]) )
689 {
690 graphmark[k] = FALSE;
691 for( l = graph->outbeg[k]; l != EAT_LAST; l = graph->oeat[l] )
692 {
693 if( !Is_term(graph->term[graph->head[l]]) )
694 {
695 assert(graph->head[l] != root);
696 pinned[graph->head[l]] = TRUE;
697 }
698 }
699 }
700 }
701 if( probtype != STP_RPCSPG )
702 graphmark[root] = FALSE;
703 }
704
705 /* for each node, store all of its outgoing boundary-edges in a (respective) heap*/
706 for( e = 0; e < nedges; e += 2 )
707 {
708 if( graph->oeat[e] == EAT_FREE )
709 continue;
710
711 node = graph->tail[e];
712 adjnode = graph->head[e];
713 newedges[e] = UNKNOWN;
714 newedges[e + 1] = UNKNOWN;
715
716 /* is edge 'e' a boundary-edge? */
717 if( vbase[node] != vbase[adjnode] && graphmark[node] && graphmark[adjnode] )
718 {
719 edgecost = vnoi[node].dist + graph->cost[e] + vnoi[adjnode].dist;
720
721 /* add the boundary-edge 'e' and its reversed to the corresponding heaps */
722 SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[vbase[node]], e, edgecost, &(heapsize[vbase[node]])) );
723 SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[vbase[adjnode]], flipedge(e), edgecost, &(heapsize[vbase[adjnode]])) );
724 }
725 }
726
727 /* find LCAs for all edges */
728 SCIP_CALL( lca(scip, graph, root, &uf, nodesmark, best_result, lvledges_start, boundpaths, heapsize, vbase) );
729
730 /* henceforth, the union-find structure will be used on the ST */
731 SCIPStpunionfindClear(scip, &uf, nnodes);
732
733 /* henceforth, nodesmark will be used to mark the current supervertices (except for the one representing the root-component) */
734 for( i = 0; dfstree[i] != root; i++ )
735 nodesmark[dfstree[i]] = FALSE;
736 nodesmark[dfstree[i]] = FALSE;
737
738 for( k = 0; k < nnodes; k++ )
739 assert(!nodesmark[k]);
740
741 /* main loop visiting all nodes of the current ST in post-order */
742 for( i = 0; dfstree[i] != root; i++ )
743 {
744 crucnode = dfstree[i];
745 scanned[crucnode] = TRUE;
746
747 SCIPdebugMessage("iteration %d (crucial node: %d) \n", i, crucnode);
748
749 /* has the node been temporarily removed from the ST? */
750 if( !graphmark[crucnode] )
751 continue;
752
753 /* is node 'crucnode' a removable crucial node? (i.e. not pinned or a terminal) */
754 if( !pinned[crucnode] && !Is_term(graph->term[crucnode]) && nodeIsCrucial(graph, best_result, crucnode) )
755 {
756 for( k = 0; k < nnodes; k++ )
757 assert(state[k] == CONNECT || !graphmark[k]);
758
759 /* find all (unique) key-paths starting in node 'crucnode' */
760 k = UNKNOWN;
761 kpcost = 0.0;
762 nkpnodes = 0;
763 nkpedges = 0;
764 nsupernodes = 0;
765 for( edge = graph->outbeg[crucnode]; edge != EAT_LAST; edge = graph->oeat[edge] )
766 {
767 /* check whether the outgoing edge is in the ST */
768 if( (best_result[edge] > -1 && steinertree[graph->head[edge]]) || (best_result[flipedge(edge)] > -1 && steinertree[graph->tail[edge]]) )
769 {
770 kpcost += graph->cost[edge];
771
772 /* check whether the current edge leads to the ST root*/
773 if( best_result[flipedge(edge)] > -1 )
774 {
775 k = flipedge(edge);
776 kpedges[nkpedges++] = k;
777 assert( edge == nodes[crucnode].edge );
778 }
779 else
780 {
781 kpedges[nkpedges++] = edge;
782 adjnode = graph->head[edge];
783 e = edge;
784
785 /* move along the key-path until its end (i.e. a crucial or pinned node) is reached */
786 while( !pinned[adjnode] && !nodeIsCrucial(graph, best_result, adjnode) && steinertree[adjnode] )
787 {
788 /* update the union-find data structure */
789 SCIPStpunionfindUnion(&uf, crucnode, adjnode, FALSE);
790
791 kpnodes[nkpnodes++] = adjnode;
792
793 for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
794 {
795 if( best_result[e] > -1 )
796 {
797 kpcost += graph->cost[e];
798 kpedges[nkpedges++] = e;
799 break;
800 }
801 }
802
803 /* assert that each leaf of the ST is a terminal */
804
805
806 if( e == EAT_LAST )
807 {
808 localmoves = 0;
809
810 goto TERMINATE;
811 }
812 assert(e != EAT_LAST);
813 adjnode = graph->head[e];
814 }
815 /* does the last node on the path belong to a removed component? */
816 if( !steinertree[adjnode] )
817 {
818 kpcost -= graph->cost[e];
819 nkpedges--;
820 adjnode = graph->tail[e];
821 if( adjnode != crucnode )
822 {
823 supernodes[nsupernodes++] = adjnode;
824 nodesmark[adjnode] = TRUE;
825 }
826 }
827 else
828 {
829 supernodes[nsupernodes++] = adjnode;
830 nodesmark[adjnode] = TRUE;
831 }
832 }
833 }
834 }
835
836 /* traverse the key-path leading to the root-component */
837 rootpathstart = nkpnodes;
838 if( k != -1 )
839 {
840 /* begin with the edge starting in the root-component of node 'crucnode' */
841 e = k;
842 adjnode = graph->tail[e];
843 while( !pinned[adjnode] && !nodeIsCrucial(graph, best_result, adjnode) && steinertree[adjnode] )
844 {
845 /* update the union-find data structure */
846 kpnodes[nkpnodes++] = adjnode;
847
848 for( e = graph->inpbeg[adjnode]; e != EAT_LAST; e = graph->ieat[e] )
849 {
850 if( best_result[e] > -1 )
851 {
852 assert(steinertree[graph->tail[e]]);
853 kpcost += graph->cost[e];
854 kpedges[nkpedges++] = e;
855 break;
856 }
857 }
858
859 assert( e != EAT_LAST );
860 adjnode = graph->tail[e];
861 }
862 supernodes[nsupernodes++] = adjnode;
863 }
864
865 /* the last of the key-path nodes to be stored is the current key-node */
866 kpnodes[nkpnodes++] = crucnode;
867
868 /* number of reset nodes */
869 nresnodes = 0;
870
871 /* reset all nodes (referred to as 'C' henceforth) whose bases are internal nodes of the current key-paths */
872 for( k = 0; k < nkpnodes; k++ )
873 {
874 /* reset all nodes having the current (internal) keypath node as their voronoi base */
875 blists_curr = blists_start[kpnodes[k]];
876 while( blists_curr != NULL )
877 {
878 node = blists_curr->index;
879 assert(graphmark[node]);
880
881 /* store all relevant data */
882 memvbase[nresnodes] = vbase[node];
883 memdist[nresnodes] = vnoi[node].dist;
884 meminedges[nresnodes] = vnoi[node].edge;
885 nresnodes++;
886
887 /* reset data */
888 vbase[node] = UNKNOWN;
889 vnoi[node].dist = FARAWAY;
890 vnoi[node].edge = UNKNOWN;
891 state[node] = UNKNOWN;
892 blists_curr = blists_curr->parent;
893 }
894 }
895
896 /* add vertical boundary-paths between the child components and the root-component (wrt node 'crucnode') */
897 nboundedges = 0;
898 for( k = 0; k < nsupernodes - 1; k++ )
899 {
900 l = supernodes[k];
901 edge = UNKNOWN;
902 while( boundpaths[l] != NULL )
903 {
904 SCIP_CALL( SCIPpairheapDeletemin(scip, &edge, &edgecost, &boundpaths[l], &heapsize[l]) );
905
906 node = (vbase[graph->head[edge]] == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(&uf, vbase[graph->head[edge]]);
907
908 /* check whether edge 'edge' represents a boundary-path having an endpoint in the kth-component and in the root-component respectively */
909 if( node != UNKNOWN && !nodesmark[node] && graphmark[node] )
910 {
911 boundedges[nboundedges++] = edge;
912 SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[l], edge, edgecost, &heapsize[l]) );
913 break;
914 }
915 }
916 }
917
918 /* add horizontal boundary-paths (between the child-components) */
919 lvledges_curr = lvledges_start[crucnode];
920 while( lvledges_curr != NULL )
921 {
922 edge = lvledges_curr->index;
923 k = vbase[graph->tail[edge]];
924 l = vbase[graph->head[edge]];
925 node = (l == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(&uf, l);
926 adjnode = (k == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(&uf, k);
927
928 /* check whether the current boundary-path connects two child components */
929 if( node != UNKNOWN && nodesmark[node] && adjnode != UNKNOWN && nodesmark[adjnode] )
930 {
931 assert(graphmark[node]);
932 assert(graphmark[adjnode]);
933 boundedges[nboundedges++] = edge;
934 }
935 lvledges_curr = lvledges_curr->parent;
936 }
937
938 /* try to connect the nodes of C (directly) to COMP(C), as a preprocessing for graph_voronoiRepair */
939 count = 0;
940 for( k = 0; k < nkpnodes; k++ )
941 {
942 blists_curr = blists_start[kpnodes[k]];
943 assert( blists_curr != NULL );
944 while( blists_curr != NULL )
945 {
946 node = blists_curr->index;
947
948 /* iterate through all outgoing edges of 'node' */
949 for( edge = graph->inpbeg[node]; edge != EAT_LAST; edge = graph->ieat[edge] )
950 {
951 adjnode = graph->tail[edge];
952
953 /* check whether the adjacent node is not in C and allows a better voronoi assignment of the current node */
954 if( state[adjnode] == CONNECT && SCIPisGT(scip, vnoi[node].dist, vnoi[adjnode].dist + graph->cost[edge])
955 && graphmark[vbase[adjnode]] && graphmark[adjnode] )
956 {
957 vnoi[node].dist = vnoi[adjnode].dist + graph->cost[edge];
958 vbase[node] = vbase[adjnode];
959 vnoi[node].edge = edge;
960 }
961 }
962 if( vbase[node] != UNKNOWN )
963 {
964 heap_add(graph->path_heap, state, &count, node, vnoi);
965 }
966 blists_curr = blists_curr->parent;
967 }
968 }
969
970 /* if there are no key-path nodes, something has gone wrong */
971 assert(nkpnodes != 0);
972
973 graph_voronoiRepairMult(scip, graph, graph->cost, &count, vbase, boundedges, &nboundedges, nodesmark, &uf, vnoi);
974
975 /* create a supergraph, having the endpoints of the key-paths incident to the current crucial node as (super-) vertices */
976 SCIP_CALL( graph_init(scip, &supergraph, nsupernodes, nboundedges * 2, 1) );
977 supergraph->stp_type = STP_SPG;
978
979 /* add vertices to the supergraph */
980 for( k = 0; k < nsupernodes; k++ )
981 {
982 supernodesid[supernodes[k]] = k;
983 graph_knot_add(supergraph, graph->term[supernodes[k]]);
984 }
985
986 /* the (super-) vertex representing the current root-component of the ST */
987 k = supernodes[nsupernodes - 1];
988
989 /* add edges to the supergraph */
990 for( l = 0; l < nboundedges; l++ )
991 {
992 edge = boundedges[l];
993 node = SCIPStpunionfindFind(&uf, vbase[graph->tail[edge]]);
994 adjnode = SCIPStpunionfindFind(&uf, vbase[graph->head[edge]]);
995
996 /* if node 'node' or 'adjnode' belongs to the root-component, take the (temporary) root-component identifier instead */
997 node = ((nodesmark[node])? node : k);
998 adjnode = ((nodesmark[adjnode])? adjnode : k);
999
1000 /* compute the cost of the boundary-path pertaining to the boundary-edge 'edge' */
1001 edgecost = vnoi[graph->tail[edge]].dist + graph->cost[edge] + vnoi[graph->head[edge]].dist;
1002 graph_edge_add(scip, supergraph, supernodesid[node], supernodesid[adjnode], edgecost, edgecost);
1003 }
1004
1005 /* compute a MST on the supergraph */
1006 SCIP_CALL( SCIPallocBufferArray(scip, &mst, nsupernodes) );
1007 SCIP_CALL( graph_path_init(scip, supergraph) );
1008 graph_path_exec(scip, supergraph, MST_MODE, nsupernodes - 1, supergraph->cost, mst);
1009
1010 /* compute the cost of the MST */
1011 mstcost = 0.0;
1012
1013 /* compute the cost of the MST */
1014 for( l = 0; l < nsupernodes - 1; l++ )
1015 {
1016 /* compute the edge in the original graph corresponding to the current MST edge */
1017 if( mst[l].edge % 2 == 0 )
1018 edge = boundedges[mst[l].edge / 2 ];
1019 else
1020 edge = flipedge(boundedges[mst[l].edge / 2 ]);
1021
1022 mstcost += graph->cost[edge];
1023 assert( newedges[edge] != crucnode && newedges[flipedge(edge)] != crucnode );
1024
1025 /* mark the edge (in the original graph) as visited */
1026 newedges[edge] = crucnode;
1027
1028 /* traverse along the boundary-path belonging to the boundary-edge 'edge' */
1029 for( node = graph->tail[edge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1030 {
1031 e = vnoi[node].edge;
1032
1033 /* if edge 'e' and its reversed have not been visited yet */
1034 if( newedges[e] != crucnode && newedges[flipedge(e)] != crucnode )
1035 {
1036 newedges[e] = crucnode;
1037 mstcost += graph->cost[e];
1038 }
1039 }
1040 for( node = graph->head[edge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1041 {
1042 e = flipedge(vnoi[node].edge);
1043
1044 /* if edge 'e' and its reversed have not been visited yet */
1045 if( newedges[vnoi[node].edge] != crucnode && newedges[e] != crucnode )
1046 {
1047 newedges[e] = crucnode;
1048 mstcost += graph->cost[e];
1049 }
1050 }
1051 }
1052
1053 if( SCIPisLT(scip, mstcost, kpcost) )
1054 {
1055 localmoves++;
1056 SCIPdebugMessage("found improving solution in KEY VERTEX ELIMINATION (round: %d) \n ", nruns);
1057
1058 /* unmark the original edges spanning the supergraph */
1059 for( e = 0; e < nkpedges; e++ )
1060 {
1061 assert(best_result[kpedges[e]] != -1);
1062 best_result[kpedges[e]] = -1;
1063 }
1064
1065 /* mark all ST nodes except for those belonging to the root-component as forbidden */
1066 for( k = rootpathstart; k < nkpnodes; k++ )
1067 {
1068 graphmark[kpnodes[k]] = FALSE;
1069 steinertree[kpnodes[k]] = FALSE;
1070 }
1071
1072 for( k = 0; k < i; k++ )
1073 {
1074 node = SCIPStpunionfindFind(&uf, dfstree[k]);
1075 if( nodesmark[node] || node == crucnode )
1076 {
1077 graphmark[dfstree[k]] = FALSE;
1078 steinertree[dfstree[k]] = FALSE;
1079 }
1080 }
1081
1082 /* add the new edges reconnecting the (super-) components */
1083 for( l = 0; l < nsupernodes - 1; l++ )
1084 {
1085 if( mst[l].edge % 2 == 0 )
1086 edge = boundedges[mst[l].edge / 2 ];
1087 else
1088 edge = flipedge(boundedges[mst[l].edge / 2 ]);
1089
1090 /* change the orientation within the target-component if necessary */
1091 if( !nodesmark[vbase[graph->head[edge]]] )
1092 {
1093 node = vbase[graph->head[edge]];
1094 k = SCIPStpunionfindFind(&uf, node);
1095 assert(nodesmark[k]);
1096 while( node != k )
1097 {
1098 /* the ST edge pointing towards the root */
1099 e = nodes[node].edge;
1100
1101 assert(best_result[e] == -1 && best_result[flipedge(e)] != -1 );
1102 best_result[e] = CONNECT;
1103 best_result[flipedge(e)] = UNKNOWN;
1104 node = graph->head[e];
1105 }
1106 }
1107
1108 /* is the vbase of the current boundary-edge tail in the root-component? */
1109 if( !nodesmark[SCIPStpunionfindFind(&uf, vbase[graph->tail[edge]])] )
1110 {
1111
1112 best_result[edge] = CONNECT;
1113
1114 for( node = graph->tail[edge], adjnode = graph->head[edge]; node != vbase[node]; adjnode = node, node = graph->tail[vnoi[node].edge] )
1115 {
1116 graphmark[node] = FALSE;
1117
1118 if( best_result[flipedge(vnoi[node].edge)] == CONNECT )
1119 best_result[flipedge(vnoi[node].edge)] = UNKNOWN;
1120
1121 best_result[vnoi[node].edge] = CONNECT;
1122 }
1123
1124 assert(!nodesmark[node] && vbase[node] == node);
1125 assert( graphmark[node] == TRUE );
1126
1127 /* is the pinned node its own component identifier? */
1128 if( !Is_term(graph->term[node]) && scanned[node] && !pinned[node] && SCIPStpunionfindFind(&uf, node) == node )
1129 {
1130 graphmark[graph->head[edge]] = FALSE;
1131 oldedge = edge;
1132
1133 for( edge = graph->outbeg[node]; edge != EAT_LAST; edge = graph->oeat[edge] )
1134 {
1135 adjnode = graph->head[edge];
1136 /* check whether edge 'edge' leads to an ancestor of terminal 'node' */
1137 if( best_result[edge] == CONNECT && graphmark[adjnode] && steinertree[adjnode] && SCIPStpunionfindFind(&uf, adjnode) != node )
1138 {
1139
1140 assert(scanned[adjnode]);
1141 /* meld the heaps */
1142 SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &heapsize[node], &heapsize[adjnode]);
1143
1144 /* update the union-find data structure */
1145 SCIPStpunionfindUnion(&uf, node, adjnode, FALSE);
1146
1147 /* move along the key-path until its end (i.e. until a crucial node is reached) */
1148 while( !nodeIsCrucial(graph, best_result, adjnode) && !pinned[adjnode] )
1149 {
1150 for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
1151 {
1152 if( best_result[e] != -1 )
1153 break;
1154 }
1155
1156 /* assert that each leaf of the ST is a terminal */
1157 assert( e != EAT_LAST );
1158 adjnode = graph->head[e];
1159 if( !steinertree[adjnode] )
1160 break;
1161 assert(scanned[adjnode]);
1162 assert(SCIPStpunionfindFind(&uf, adjnode) != node);
1163
1164 /* update the union-find data structure */
1165 SCIPStpunionfindUnion(&uf, node, adjnode, FALSE);
1166
1167 /* meld the heaps */
1168 SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &heapsize[node], &heapsize[adjnode]);
1169 }
1170 }
1171 }
1172 edge = oldedge;
1173 }
1174
1175 /* mark the start node (lying in the root-component of the ST) of the current boundary-path as pinned,
1176 * so that it may not be removed later on */
1177 pinned[node] = TRUE;
1178
1179 for( node = graph->head[edge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1180 {
1181 graphmark[node] = FALSE;
1182 if( best_result[vnoi[node].edge] == CONNECT )
1183 best_result[vnoi[node].edge] = -1;
1184
1185 best_result[flipedge(vnoi[node].edge)] = CONNECT;
1186
1187 }
1188 }
1189 else
1190 {
1191
1192 best_result[edge] = CONNECT;
1193
1194 for( node = graph->tail[edge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1195 {
1196 graphmark[node] = FALSE;
1197 if( best_result[vnoi[node].edge] != CONNECT && best_result[flipedge(vnoi[node].edge)] != CONNECT )
1198 best_result[vnoi[node].edge] = CONNECT;
1199
1200 }
1201
1202 for( node = graph->head[edge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1203 {
1204 graphmark[node] = FALSE;
1205
1206 best_result[flipedge(vnoi[node].edge)] = CONNECT;
1207 best_result[vnoi[node].edge] = UNKNOWN;
1208 }
1209 }
1210 }
1211
1212 for( k = 0; k < nkpnodes; k++ )
1213 {
1214 assert(graphmark[kpnodes[k]] == FALSE);
1215 assert(steinertree[kpnodes[k]] == FALSE);
1216 }
1217 assert(!graphmark[crucnode]);
1218 }
1219 else
1220 {
1221 /* no improving solution has been found during the move */
1222
1223 /* meld the heap pertaining to 'crucnode' and all heaps pertaining to descendant key-paths of node 'crucnode' */
1224 for( k = 0; k < rootpathstart; k++ )
1225 {
1226 SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[kpnodes[k]], &heapsize[crucnode], &heapsize[kpnodes[k]]);
1227 }
1228 for( k = 0; k < nsupernodes - 1; k++ )
1229 {
1230 SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[supernodes[k]], &heapsize[crucnode], &heapsize[supernodes[k]]);
1231
1232 /* update the union-find data structure */
1233 SCIPStpunionfindUnion(&uf, crucnode, supernodes[k], FALSE);
1234 }
1235 }
1236
1237 /* free the supergraph and the MST data structure */
1238 graph_path_exit(scip, supergraph);
1239 graph_free(scip, &supergraph, TRUE);
1240 SCIPfreeBufferArray(scip, &mst);
1241
1242 /* unmark the descendant supervertices */
1243 for( k = 0; k < nsupernodes - 1; k++ )
1244 {
1245 nodesmark[supernodes[k]] = FALSE;
1246 }
1247
1248 /* debug test; to be deleted later on */
1249 for( k = 0; k < nnodes; k++ )
1250 {
1251 assert( !nodesmark[k] );
1252 }
1253
1254 /* restore the original voronoi diagram */
1255 l = 0;
1256 for( k = 0; k < nkpnodes; k++ )
1257 {
1258 /* restore data of all nodes having the current (internal) key-path node as their voronoi base */
1259 blists_curr = blists_start[kpnodes[k]];
1260 while( blists_curr != NULL )
1261 {
1262 node = blists_curr->index;
1263 vbase[node] = memvbase[l];
1264 vnoi[node].dist = memdist[l];
1265 vnoi[node].edge = meminedges[l];
1266 l++;
1267 blists_curr = blists_curr->parent;
1268 }
1269 }
1270
1271 assert(l == nresnodes);
1272 }
1273
1274 /** Key-Path Exchange */
1275 if( probtype != STP_MWCSP )
1276 {
1277 /* if the node has just been eliminated, skip Key-Path Exchange */
1278 if( !graphmark[crucnode] )
1279 continue;
1280
1281 /* is crucnode a crucial or pinned vertex? */
1282 if( (!nodeIsCrucial(graph, best_result, crucnode) && !pinned[crucnode]) )
1283 continue;
1284
1285 if( Is_term(graph->term[crucnode]) || pinned[crucnode] )
1286 {
1287 for( edge = graph->outbeg[crucnode]; edge != EAT_LAST; edge = graph->oeat[edge] )
1288 {
1289 adjnode = graph->head[edge];
1290 /* check whether edge 'edge' leads to an ancestor of terminal 'crucnode' */
1291 if( best_result[edge] == CONNECT && steinertree[adjnode] && graphmark[adjnode] )
1292 {
1293 assert( SCIPStpunionfindFind(&uf, adjnode) != crucnode);
1294 assert(scanned[adjnode]);
1295 /* meld the heaps */
1296 SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[adjnode], &heapsize[crucnode], &heapsize[adjnode]);
1297
1298 /* update the union-find data structure */
1299 SCIPStpunionfindUnion(&uf, crucnode, adjnode, FALSE);
1300
1301 /* move along the key-path until its end (i.e. until a crucial node is reached) */
1302 while( !nodeIsCrucial(graph, best_result, adjnode) && !pinned[adjnode] )
1303 {
1304 for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
1305 {
1306 if( best_result[e] != -1 )
1307 break;
1308 }
1309
1310 /* assert that each leaf of the ST is a terminal */
1311 assert( e != EAT_LAST );
1312 adjnode = graph->head[e];
1313 if( !steinertree[adjnode] || !graphmark[adjnode] )
1314 break;
1315 assert(scanned[adjnode]);
1316 assert(SCIPStpunionfindFind(&uf, adjnode) != crucnode);
1317
1318 /* update the union-find data structure */
1319 SCIPStpunionfindUnion(&uf, crucnode, adjnode, FALSE);
1320
1321 /* meld the heaps */
1322 SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[adjnode], &heapsize[crucnode], &heapsize[adjnode]);
1323 }
1324 }
1325 }
1326
1327 }
1328
1329 /* counts the internal nodes of the keypath */
1330 nkpnodes = 0;
1331
1332 for( k = 0; k < nnodes; k++ )
1333 assert(state[k] == CONNECT || !graphmark[k]);
1334
1335 /* find the (unique) key-path containing the parent of the current crucial node 'crucnode' */
1336 kptailnode = graph->head[nodes[crucnode].edge];
1337 kpathcost = graph->cost[nodes[crucnode].edge];
1338
1339 while( !nodeIsCrucial(graph, best_result, kptailnode) && !pinned[kptailnode] )
1340 {
1341 kpathcost += graph->cost[nodes[kptailnode].edge];
1342
1343 kpnodes[nkpnodes++] = kptailnode;
1344 kptailnode = graph->head[nodes[kptailnode].edge];
1345 }
1346
1347 /* counts the reset nodes during voronoi repair */
1348 nresnodes = 0;
1349
1350 /* reset all nodes (henceforth referred to as 'C') whose bases are internal nodes of the current keypath */
1351 for( k = 0; k < nkpnodes; k++ )
1352 {
1353 /* reset all nodes having the current (internal) keypath node as their voronoi base */
1354 blists_curr = blists_start[kpnodes[k]];
1355 while( blists_curr != NULL )
1356 {
1357 node = blists_curr->index;
1358 memvbase[nresnodes] = vbase[node];
1359 memdist[nresnodes] = vnoi[node].dist;
1360 meminedges[nresnodes] = vnoi[node].edge;
1361 nresnodes++;
1362 vbase[node] = UNKNOWN;
1363 vnoi[node].dist = FARAWAY;
1364 vnoi[node].edge = UNKNOWN;
1365 state[node] = UNKNOWN;
1366 blists_curr = blists_curr->parent;
1367 }
1368 }
1369
1370 edgecost = UNKNOWN;
1371 e = UNKNOWN;
1372 while( boundpaths[crucnode] != NULL )
1373 {
1374 SCIP_CALL( SCIPpairheapDeletemin(scip, &e, &edgecost, &boundpaths[crucnode], &(heapsize[crucnode])) );
1375 assert( e != UNKNOWN );
1376 k = vbase[graph->tail[e]];
1377 l = vbase[graph->head[e]];
1378
1379 assert(graphmark[k]);
1380 node = (l == UNKNOWN || !graphmark[l] )? UNKNOWN : SCIPStpunionfindFind(&uf, l);
1381
1382 /* does the boundary-path end in the root component? */
1383 if( node != UNKNOWN && node != crucnode && graphmark[l] )
1384 {
1385 SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[crucnode], e, edgecost, &(heapsize[crucnode])) );
1386 break;
1387 }
1388 }
1389
1390 if( boundpaths[crucnode] == NULL )
1391 {
1392 oldedge = UNKNOWN;
1393 }
1394 else
1395 {
1396 oldedge = e;
1397 }
1398
1399 /* counts the nodes connected during the following 'preprocessing' */
1400 count = 0;
1401
1402 /* try to connect the nodes of C (directly) to COMP(C), as a preprocessing for voronoi-repair */
1403 for( k = 0; k < nkpnodes; k++ )
1404 {
1405 blists_curr = blists_start[kpnodes[k]];
1406 assert( blists_curr != NULL );
1407 while( blists_curr != NULL )
1408 {
1409 node = blists_curr->index;
1410
1411 /* iterate through all outgoing edges of 'node' */
1412 for( edge = graph->inpbeg[node]; edge != EAT_LAST; edge = graph->ieat[edge] )
1413 {
1414 adjnode = graph->tail[edge];
1415
1416 /* check whether the adjacent node is not in C and allows a better voronoi assignment of the current node */
1417 if( state[adjnode] == CONNECT && SCIPisGT(scip, vnoi[node].dist, vnoi[adjnode].dist + graph->cost[edge])
1418 && graphmark[vbase[adjnode]] && graphmark[adjnode] )
1419 {
1420 vnoi[node].dist = vnoi[adjnode].dist + graph->cost[edge];
1421 vbase[node] = vbase[adjnode];
1422 vnoi[node].edge = edge;
1423 }
1424 }
1425 if( vbase[node] != UNKNOWN )
1426 {
1427 heap_add(graph->path_heap, state, &count, node, vnoi);
1428 }
1429 blists_curr = blists_curr->parent;
1430 }
1431 }
1432 if( nkpnodes > 0 )
1433 assert(count > 0);
1434 newedge = UNKNOWN;
1435
1436 /* if there is no key path, nothing has to be repaired */
1437 if( nkpnodes > 0 )
1438 graph_voronoiRepair(scip, graph, graph->cost, &count, vbase, vnoi, &newedge, crucnode, &uf);
1439 else
1440 newedge = nodes[crucnode].edge;
1441
1442 if( oldedge != UNKNOWN && newedge != UNKNOWN && SCIPisLT(scip, edgecost,
1443 vnoi[graph->tail[newedge]].dist + graph->cost[newedge] + vnoi[graph->head[newedge]].dist) )
1444 newedge = oldedge;
1445 if( oldedge != UNKNOWN && newedge == UNKNOWN )
1446 newedge = oldedge;
1447
1448 assert( newedge != UNKNOWN );
1449 edgecost = vnoi[graph->tail[newedge]].dist + graph->cost[newedge] + vnoi[graph->head[newedge]].dist;
1450 if( SCIPisLT(scip, edgecost, kpathcost) )
1451 {
1452 node = SCIPStpunionfindFind(&uf, vbase[graph->head[newedge]]);
1453
1454 SCIPdebugMessage( "ADDING NEW KEY PATH (%f )\n", edgecost - kpathcost );
1455
1456 localmoves++;
1457
1458 /* remove old keypath */
1459 assert( best_result[flipedge(nodes[crucnode].edge)] != UNKNOWN );
1460 best_result[flipedge(nodes[crucnode].edge)] = UNKNOWN;
1461 steinertree[crucnode] = FALSE;
1462 graphmark[crucnode] = FALSE;
1463
1464 for( k = 0; k < nkpnodes; k++ )
1465 {
1466 assert( best_result[flipedge(nodes[kpnodes[k]].edge)] != UNKNOWN );
1467 best_result[flipedge(nodes[kpnodes[k]].edge)] = UNKNOWN;
1468 steinertree[kpnodes[k]] = FALSE;
1469 graphmark[kpnodes[k]] = FALSE;
1470 }
1471 assert(graphmark[kptailnode]);
1472
1473 if( node == crucnode )
1474 newedge = flipedge(newedge);
1475
1476 for( node = graph->tail[newedge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1477 {
1478 graphmark[node] = FALSE;
1479
1480 best_result[flipedge(vnoi[node].edge)] = CONNECT;
1481 best_result[vnoi[node].edge] = UNKNOWN;
1482 }
1483
1484 for( node = graph->head[newedge]; node != vbase[node]; node = graph->tail[vnoi[node].edge] )
1485 {
1486 graphmark[node] = FALSE;
1487
1488 best_result[vnoi[node].edge] = CONNECT;
1489 }
1490
1491 best_result[flipedge(newedge)] = CONNECT;
1492
1493 newpathend = vbase[graph->tail[newedge]];
1494 assert(node == vbase[graph->head[newedge]] );
1495
1496 /* flip all edges on the ST path between the endnode of the new key-path and the current crucial node */
1497 k = newpathend;
1498 assert(SCIPStpunionfindFind(&uf, newpathend) == crucnode);
1499
1500 while( k != crucnode )
1501 {
1502 assert(graphmark[k]);
1503 assert( best_result[flipedge(nodes[k].edge)] != -1);
1504 best_result[flipedge(nodes[k].edge)] = UNKNOWN;
1505
1506 best_result[nodes[k].edge] = CONNECT;
1507
1508 k = graph->head[nodes[k].edge];
1509 }
1510
1511 for( k = 0; k < i; k++ )
1512 {
1513 if( crucnode == SCIPStpunionfindFind(&uf, dfstree[k]) )
1514 {
1515 graphmark[dfstree[k]] = FALSE;
1516 steinertree[dfstree[k]] = FALSE;
1517 }
1518 }
1519
1520 /* update union find */
1521 if( !Is_term(graph->term[node]) && scanned[node] && !pinned[node] && SCIPStpunionfindFind(&uf, node) == node )
1522 {
1523 for( edge = graph->outbeg[node]; edge != EAT_LAST; edge = graph->oeat[edge] )
1524 {
1525 adjnode = graph->head[edge];
1526 /* check whether edge 'edge' leads to an ancestor of terminal 'node' */
1527 if( best_result[edge] == CONNECT && steinertree[adjnode] && graphmark[adjnode] && SCIPStpunionfindFind(&uf, adjnode) != node )
1528 {
1529 assert(scanned[adjnode]);
1530 /* meld the heaps */
1531 SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &heapsize[node], &heapsize[adjnode]);
1532
1533 /* update the union-find data structure */
1534 SCIPStpunionfindUnion(&uf, node, adjnode, FALSE);
1535
1536 /* move along the key-path until its end (i.e. until a crucial node is reached) */
1537 while( !nodeIsCrucial(graph, best_result, adjnode) && !pinned[adjnode] )
1538 {
1539 for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
1540 {
1541 if( best_result[e] != -1 )
1542 break;
1543 }
1544
1545 /* assert that each leaf of the ST is a terminal */
1546 assert( e != EAT_LAST );
1547 adjnode = graph->head[e];
1548 if( !steinertree[adjnode] )
1549 break;
1550 assert(scanned[adjnode]);
1551 assert(SCIPStpunionfindFind(&uf, adjnode) != node);
1552
1553 /* update the union-find data structure */
1554 SCIPStpunionfindUnion(&uf, node, adjnode, FALSE);
1555
1556 /* meld the heaps */
1557 SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &heapsize[node], &heapsize[adjnode]);
1558 }
1559 }
1560 }
1561
1562 }
1563 pinned[node] = TRUE;
1564 }
1565
1566 /* restore the original voronoi digram */
1567 l = 0;
1568 for( k = 0; k < nkpnodes; k++ )
1569 {
1570 /* reset all nodes having the current (internal) keypath node as their voronoi base */
1571 blists_curr = blists_start[kpnodes[k]];
1572 while( blists_curr != NULL )
1573 {
1574 node = blists_curr->index;
1575 vbase[node] = memvbase[l];
1576 vnoi[node].dist = memdist[l];
1577 vnoi[node].edge = meminedges[l];
1578 l++;
1579 blists_curr = blists_curr->parent;
1580 }
1581 }
1582 assert(l == nresnodes);
1583 }
1584 }
1585
1586
1587 /**********************************************************/
1588
1589 TERMINATE:
1590
1591 SCIPStpunionfindClear(scip, &uf, nnodes);
1592
1593 /* free data structures */
1594 SCIPfreeBufferArray(scip, &kpedges);
1595 SCIPfreeBufferArray(scip, &kpnodes);
1596 SCIPfreeBufferArray(scip, &supernodes);
1597
1598 for( k = nnodes - 1; k >= 0; k-- )
1599 {
1600 if( boundpaths[k] != NULL )
1601 SCIPpairheapFree(scip, &boundpaths[k]);
1602
1603 blists_curr = blists_start[k];
1604 lvledges_curr = lvledges_start[k];
1605 while( lvledges_curr != NULL )
1606 {
1607 lvledges_start[k] = lvledges_curr->parent;
1608 SCIPfreeBlockMemory(scip, &lvledges_curr);
1609 lvledges_curr = lvledges_start[k];
1610 }
1611
1612 while( blists_curr != NULL )
1613 {
1614 blists_start[k] = blists_curr->parent;
1615 SCIPfreeBlockMemory(scip, &blists_curr);
1616 blists_curr = blists_start[k];
1617 }
1618 }
1619
1620 /* has there been a move during this run? */
1621 if( localmoves > 0 )
1622 {
1623 for( i = 0; i < nnodes; i++ )
1624 {
1625 steinertree[i] = FALSE;
1626 graphmark[i] = (graph->grad[i] > 0);
1627 SCIPlinkcuttreeInit(&nodes[i]);
1628 }
1629
1630 graphmark[root] = TRUE;
1631
1632 /* create a link-cut tree representing the current Steiner tree */
1633 for( e = 0; e < nedges; e++ )
1634 {
1635 assert(graph->head[e] == graph->tail[flipedge(e)]);
1636
1637 /* if edge e is in the tree, so are its incident vertices */
1638 if( best_result[e] != -1 )
1639 {
1640 steinertree[graph->tail[e]] = TRUE;
1641 steinertree[graph->head[e]] = TRUE;
1642 SCIPlinkcuttreeLink(&nodes[graph->head[e]], &nodes[graph->tail[e]], flipedge(e));
1643 }
1644 }
1645 assert( nodes[root].edge == -1 );
1646 nodes[root].edge = -1;
1647 }
1648 }
1649
1650 /* free data structures */
1651 SCIPStpunionfindFree(scip, &uf);
1652 SCIPfreeBufferArray(scip, &nodesmark);
1653 SCIPfreeBufferArray(scip, &dfstree);
1654 SCIPfreeBufferArray(scip, &pinned);
1655 SCIPfreeBufferArray(scip, &boundpaths);
1656 SCIPfreeBufferArray(scip, &meminedges);
1657 SCIPfreeBufferArray(scip, &memdist);
1658 SCIPfreeBufferArray(scip, &memvbase);
1659 SCIPfreeBufferArray(scip, &blists_start);
1660 SCIPfreeBufferArray(scip, &heapsize);
1661 SCIPfreeBufferArray(scip, &scanned);
1662 SCIPfreeBufferArray(scip, &supernodesid);
1663 SCIPfreeBufferArray(scip, &boundedges);
1664 SCIPfreeBufferArray(scip, &lvledges_start);
1665 SCIPfreeBufferArray(scip, &newedges);
1666 /******/
1667 }
1668
1669 #ifdef SCIP_DEBUG
1670 {
1671 SCIP_Real obj = 0.0;
1672 for( e = 0; e < nedges; e++ )
1673 obj += (best_result[e] > -1) ? graph->cost[e] : 0.0;
1674
1675 printf(" ObjAfterHeurLocal=%.12e\n", obj);
1676 }
1677 #endif
1678
1679 #ifndef NDEBUG
1680 assert(SCIPisLE(scip, graph_sol_getObj(graph->cost, best_result, 0.0, nedges), initialobj));
1681 #endif
1682
1683 SCIPfreeBufferArray(scip, &steinertree);
1684 SCIPfreeBufferArray(scip, &nodes);
1685 SCIPfreeBufferArray(scip, &vbase);
1686 SCIPfreeBufferArray(scip, &vnoi);
1687
1688 return SCIP_OKAY;
1689 }
1690
1691 /** Greedy Extension local heuristic for (R)PC and MW */
SCIPStpHeurLocalExtendPcMw(SCIP * scip,GRAPH * graph,const SCIP_Real * cost,PATH * path,int * stedge,STP_Bool * stvertex,SCIP_Bool * extensions)1692 SCIP_RETCODE SCIPStpHeurLocalExtendPcMw(
1693 SCIP* scip, /**< SCIP data structure */
1694 GRAPH* graph, /**< graph data structure */
1695 const SCIP_Real* cost, /**< edge cost array*/
1696 PATH* path, /**< shortest data structure array */
1697 int* stedge, /**< initialized array to indicate whether an edge is part of the Steiner tree */
1698 STP_Bool* stvertex, /**< uninitialized array to indicate whether a vertex is part of the Steiner tree */
1699 SCIP_Bool* extensions /**< pointer to store whether extensions have been made */
1700 )
1701 {
1702 GNODE candidates[MAX(GREEDY_EXTENSIONS, GREEDY_EXTENSIONS_MW)];
1703 int candidatesup[MAX(GREEDY_EXTENSIONS, GREEDY_EXTENSIONS_MW)];
1704
1705 PATH* orgpath;
1706 SCIP_PQUEUE* pqueue;
1707 SCIP_Real bestsolval;
1708 int i;
1709 int root;
1710 int nedges;
1711 int nnodes;
1712 int nextensions;
1713 int greedyextensions;
1714 STP_Bool* stvertextmp;
1715
1716 assert(scip != NULL);
1717 assert(path != NULL);
1718 assert(graph != NULL);
1719 assert(stedge != NULL);
1720 assert(cost != NULL);
1721 assert(stvertex != NULL);
1722
1723 root = graph->source;
1724 nnodes = graph->knots;
1725 nedges = graph->edges;
1726
1727 graph_pc_2transcheck(graph);
1728 SCIP_CALL( SCIPallocBufferArray(scip, &stvertextmp, nnodes) );
1729 SCIP_CALL( SCIPallocBufferArray(scip, &orgpath, nnodes) );
1730
1731 /* initialize solution vertex array with FALSE */
1732 BMSclearMemoryArray(stvertex, nnodes);
1733
1734 stvertex[graph->source] = TRUE;
1735 path[graph->source].edge = UNKNOWN;
1736
1737 for( int e = 0; e < nedges; e++ )
1738 if( stedge[e] == CONNECT )
1739 {
1740 i = graph->head[e];
1741 path[i].edge = e;
1742 stvertex[i] = TRUE;
1743 }
1744
1745 graph_path_st_pcmw_extend(scip, graph, cost, path, stvertex, extensions);
1746
1747 BMScopyMemoryArray(orgpath, path, nnodes);
1748
1749 if( graph->stp_type == STP_MWCSP )
1750 greedyextensions = GREEDY_EXTENSIONS_MW;
1751 else
1752 greedyextensions = GREEDY_EXTENSIONS;
1753
1754 /*** compute solution value and save greedyextensions many best unconnected nodes ***/
1755
1756 SCIP_CALL( SCIPpqueueCreate(&pqueue, greedyextensions, 2.0, GNODECmpByDist, NULL) );
1757
1758 bestsolval = 0.0;
1759 nextensions = 0;
1760 for( i = 0; i < nnodes; i++ )
1761 {
1762 if( graph->grad[i] == 0 )
1763 continue;
1764
1765 if( Is_term(graph->term[i]) )
1766 {
1767 if( i != root )
1768 {
1769 int e;
1770 SCIP_Bool connect = FALSE;
1771 SCIP_Real ecost = -1.0;
1772 for( e = graph->inpbeg[i]; e != EAT_LAST; e = graph->ieat[e] )
1773 {
1774 if( root == graph->tail[e] )
1775 ecost = graph->cost[e];
1776 else if( stvertex[graph->tail[e]] )
1777 connect = TRUE;
1778 }
1779
1780 if( !connect )
1781 bestsolval += ecost;
1782 }
1783 }
1784 else if( stvertex[i] )
1785 {
1786 bestsolval += graph->cost[orgpath[i].edge];
1787 }
1788 else if( orgpath[i].edge != UNKNOWN && Is_pterm(graph->term[i]) )
1789 {
1790 if( nextensions < greedyextensions )
1791 {
1792 candidates[nextensions].dist = graph->prize[i] - path[i].dist;
1793 candidates[nextensions].number = i;
1794
1795 SCIP_CALL( SCIPpqueueInsert(pqueue, &(candidates[nextensions++])) );
1796 }
1797 else
1798 {
1799 /* get candidate vertex of minimum value */
1800 GNODE* min = (GNODE*) SCIPpqueueFirst(pqueue);
1801 if( SCIPisLT(scip, min->dist, graph->prize[i] - path[i].dist) )
1802 {
1803 min = (GNODE*) SCIPpqueueRemove(pqueue);
1804 min->dist = graph->prize[i] - path[i].dist;
1805 min->number = i;
1806 SCIP_CALL( SCIPpqueueInsert(pqueue, min) );
1807 }
1808 }
1809 }
1810 }
1811
1812 for( int restartcount = 0; restartcount < GREEDY_MAXRESTARTS; restartcount++ )
1813 {
1814 int l = 0;
1815 SCIP_Bool extensionstmp = FALSE;
1816
1817 i = nextensions;
1818
1819 /* write extension candidates into array, from max to min */
1820 while( SCIPpqueueNElems(pqueue) > 0 )
1821 {
1822 GNODE* min = (GNODE*) SCIPpqueueRemove(pqueue);
1823 assert(i > 0);
1824 candidatesup[--i] = min->number;
1825 }
1826 assert(i == 0);
1827
1828 /* iteratively insert new subpaths and try to improve solution */
1829 for( ; l < nextensions; l++ )
1830 {
1831 i = candidatesup[l];
1832 if( !stvertex[i] )
1833 {
1834 SCIP_Real newsolval = 0.0;
1835 int k = i;
1836
1837 BMScopyMemoryArray(stvertextmp, stvertex, nnodes);
1838 BMScopyMemoryArray(path, orgpath, nnodes);
1839
1840 /* add new extension */
1841 while( !stvertextmp[k] )
1842 {
1843 stvertextmp[k] = TRUE;
1844 assert(orgpath[k].edge != UNKNOWN);
1845 k = graph->tail[orgpath[k].edge];
1846 }
1847
1848 graph_path_st_pcmw_extend(scip, graph, cost, path, stvertextmp, &extensionstmp);
1849
1850 for( int j = 0; j < nnodes; j++ )
1851 {
1852 if( graph->grad[j] == 0 )
1853 continue;
1854
1855 if( Is_term(graph->term[j]) )
1856 {
1857 if( j != root )
1858 {
1859 int e;
1860 SCIP_Bool connect = FALSE;
1861 SCIP_Real ecost = -1.0;
1862 for( e = graph->inpbeg[j]; e != EAT_LAST; e = graph->ieat[e] )
1863 {
1864 if( root == graph->tail[e] )
1865 ecost = graph->cost[e];
1866 else if( stvertextmp[graph->tail[e]] )
1867 connect = TRUE;
1868 }
1869
1870 if( !connect )
1871 newsolval += ecost;
1872 }
1873 }
1874 else if( stvertextmp[j] )
1875 {
1876 newsolval += graph->cost[path[j].edge];
1877 }
1878 }
1879
1880 /* new solution value better than old one? */
1881 if( SCIPisLT(scip, newsolval, bestsolval) )
1882 {
1883 *extensions = TRUE;
1884
1885 bestsolval = newsolval;
1886 BMScopyMemoryArray(stvertex, stvertextmp, nnodes);
1887
1888 /* save greedyextensions many best unconnected nodes */
1889 nextensions = 0;
1890
1891 for( int j = 0; j < nnodes; j++ )
1892 {
1893 orgpath[j].edge = path[j].edge;
1894 orgpath[j].dist = path[j].dist;
1895 if( !stvertex[j] && Is_pterm(graph->term[j]) && path[j].edge != UNKNOWN )
1896 {
1897 if( nextensions < greedyextensions )
1898 {
1899 candidates[nextensions].dist = graph->prize[j] - path[j].dist;
1900 candidates[nextensions].number = j;
1901
1902 SCIP_CALL( SCIPpqueueInsert(pqueue, &(candidates[nextensions++])) );
1903 }
1904 else
1905 {
1906 /* get candidate vertex of minimum value */
1907 GNODE* min = (GNODE*) SCIPpqueueFirst(pqueue);
1908 if( SCIPisLT(scip, min->dist, graph->prize[j] - path[j].dist) )
1909 {
1910 min = (GNODE*) SCIPpqueueRemove(pqueue);
1911 min->dist = graph->prize[j] - path[j].dist;
1912 min->number = j;
1913 SCIP_CALL( SCIPpqueueInsert(pqueue, min) );
1914 }
1915 }
1916 }
1917 }
1918
1919 break;
1920 } /* if new solution value better than old one? */
1921 } /* if !stvertex[i] */
1922 } /* for l < nextension */
1923 /* no more extensions performed? */
1924 if( l == nextensions )
1925 break;
1926 } /* main loop */
1927
1928 /* have vertices been added? */
1929 if( *extensions )
1930 {
1931 SCIPdebugMessage("SCIPStpHeurLocalExtendPcMw found extensions \n");
1932 for( int e = 0; e < nedges; e++ )
1933 stedge[e] = UNKNOWN;
1934 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, cost, stedge, stvertex) );
1935 }
1936
1937 SCIPpqueueFree(&pqueue);
1938 SCIPfreeBufferArray(scip, &orgpath);
1939 SCIPfreeBufferArray(scip, &stvertextmp);
1940
1941 #ifdef SCIP_DEBUG
1942 {
1943 SCIP_Real t = 0.0;
1944 for( int e = 0; e < nedges; e++ )
1945 if( stedge[e] == CONNECT )
1946 t += graph->cost[e];
1947
1948 SCIPdebugMessage("SCIPStpHeurLocalExtendPcMw: exit real cost %f \n", t);
1949 }
1950 #endif
1951
1952 return SCIP_OKAY;
1953 }
1954
1955
1956
1957 /*
1958 * Callback methods of primal heuristic
1959 */
1960
1961 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
1962 static
SCIP_DECL_HEURCOPY(heurCopyLocal)1963 SCIP_DECL_HEURCOPY(heurCopyLocal)
1964 { /*lint --e{715}*/
1965 assert(scip != NULL);
1966 assert(heur != NULL);
1967 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1968
1969 /* call inclusion method of primal heuristic */
1970 SCIP_CALL( SCIPStpIncludeHeurLocal(scip) );
1971
1972 return SCIP_OKAY;
1973 }
1974
1975 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
1976 static
SCIP_DECL_HEURFREE(heurFreeLocal)1977 SCIP_DECL_HEURFREE(heurFreeLocal)
1978 { /*lint --e{715}*/
1979 SCIP_HEURDATA* heurdata;
1980
1981 assert(heur != NULL);
1982 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1983 assert(scip != NULL);
1984
1985 /* free heuristic data */
1986 heurdata = SCIPheurGetData(heur);
1987 assert(heurdata != NULL);
1988 SCIPfreeMemory(scip, &heurdata);
1989 SCIPheurSetData(heur, NULL);
1990
1991 return SCIP_OKAY;
1992 }
1993
1994 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
1995 static
SCIP_DECL_HEURINITSOL(heurInitsolLocal)1996 SCIP_DECL_HEURINITSOL(heurInitsolLocal)
1997 { /*lint --e{715}*/
1998 SCIP_HEURDATA* heurdata;
1999
2000 assert(heur != NULL);
2001 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2002 assert(scip != NULL);
2003
2004 /* free heuristic data */
2005 heurdata = SCIPheurGetData(heur);
2006
2007 heurdata->nfails = 1;
2008 heurdata->nbestsols = DEFAULT_NBESTSOLS;
2009
2010 SCIP_CALL( SCIPallocMemoryArray(scip, &(heurdata->lastsolindices), heurdata->maxnsols) );
2011
2012 for( int i = 0; i < heurdata->maxnsols; i++ )
2013 heurdata->lastsolindices[i] = -1;
2014
2015 return SCIP_OKAY;
2016 }
2017
2018
2019 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */
2020 static
SCIP_DECL_HEUREXITSOL(heurExitsolLocal)2021 SCIP_DECL_HEUREXITSOL(heurExitsolLocal)
2022 { /*lint --e{715}*/
2023 SCIP_HEURDATA* heurdata;
2024
2025 assert(heur != NULL);
2026 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2027 assert(scip != NULL);
2028
2029 /* free heuristic data */
2030 heurdata = SCIPheurGetData(heur);
2031 assert(heurdata != NULL);
2032 assert(heurdata->lastsolindices != NULL);
2033 SCIPfreeMemoryArray(scip, &(heurdata->lastsolindices));
2034
2035 return SCIP_OKAY;
2036 }
2037
2038
2039
2040 /** execution method of primal heuristic */
2041 static
SCIP_DECL_HEUREXEC(heurExecLocal)2042 SCIP_DECL_HEUREXEC(heurExecLocal)
2043 { /*lint --e{715}*/
2044 SCIP_HEURDATA* heurdata;
2045 SCIP_PROBDATA* probdata;
2046 GRAPH* graph; /* graph structure */
2047 SCIP_SOL* newsol; /* new solution */
2048 SCIP_SOL* impsol; /* new improved solution */
2049 SCIP_SOL** sols; /* solutions */
2050 SCIP_VAR** vars; /* SCIP variables */
2051 SCIP_Real pobj;
2052 SCIP_Real* nval;
2053 SCIP_Real* xval;
2054 int i;
2055 int e;
2056 int v;
2057 int min;
2058 int root;
2059 int nvars;
2060 int nsols; /* number of all solutions found so far */
2061 int nedges;
2062 int* results;
2063 int* lastsolindices;
2064 SCIP_Bool feasible;
2065
2066 assert(heur != NULL);
2067 assert(scip != NULL);
2068 assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2069 assert(result != NULL);
2070
2071 /* get heuristic data */
2072 heurdata = SCIPheurGetData(heur);
2073 assert(heurdata != NULL);
2074 lastsolindices = heurdata->lastsolindices;
2075 assert(lastsolindices != NULL);
2076
2077 probdata = SCIPgetProbData(scip);
2078 assert(probdata != NULL);
2079
2080 graph = SCIPprobdataGetGraph(probdata);
2081 assert(graph != NULL);
2082
2083 *result = SCIP_DIDNOTRUN;
2084
2085 /* the local heuristics may not work correctly for several problem variants*/
2086 if( graph->stp_type != STP_SPG && graph->stp_type != STP_RSMT && graph->stp_type != STP_OARSMT &&
2087 graph->stp_type != STP_PCSPG && graph->stp_type != STP_RPCSPG && graph->stp_type != STP_GSTP
2088 && graph->stp_type != STP_MWCSP )
2089 return SCIP_OKAY;
2090
2091 /* don't run local in a Subscip */
2092 if( SCIPgetSubscipDepth(scip) > 0 )
2093 return SCIP_OKAY;
2094
2095 /* no solution available? */
2096 if( SCIPgetBestSol(scip) == NULL )
2097 return SCIP_OKAY;
2098
2099 root = graph->source;
2100 sols = SCIPgetSols(scip);
2101 nsols = SCIPgetNSols(scip);
2102 nedges = graph->edges;
2103
2104 assert(heurdata->maxnsols >= 0);
2105
2106 min = MIN(heurdata->maxnsols, nsols);
2107
2108 /* only process each solution once */
2109 for( v = 0; v < min; v++ )
2110 {
2111 if( SCIPsolGetIndex(sols[v]) != lastsolindices[v] )
2112 {
2113 /* shift all solution indices right of the new solution index */
2114 for( i = min - 1; i >= v + 1; i-- )
2115 lastsolindices[i] = lastsolindices[i - 1];
2116 break;
2117 }
2118 }
2119
2120 /* no new solution available? */
2121 if( v == min )
2122 return SCIP_OKAY;
2123
2124 newsol = sols[v];
2125 lastsolindices[v] = SCIPsolGetIndex(newsol);
2126
2127 /* solution not good enough? */
2128 if( (v > heurdata->nbestsols && !(heurdata->maxfreq)) && graph->stp_type != STP_MWCSP )
2129 return SCIP_OKAY;
2130
2131 /* has the new solution been found by this very heuristic? */
2132 if( SCIPsolGetHeur(newsol) == heur )
2133 return SCIP_OKAY;
2134
2135 *result = SCIP_DIDNOTFIND;
2136
2137 vars = SCIPprobdataGetVars(scip);
2138 nvars = SCIPprobdataGetNVars(scip);
2139 xval = SCIPprobdataGetXval(scip, newsol);
2140
2141 if( vars == NULL )
2142 return SCIP_OKAY;
2143
2144 assert(vars != NULL);
2145 assert(xval != NULL);
2146
2147 /* for PC variants: test whether solution is trivial */
2148 if( graph->stp_type == STP_PCSPG || graph->stp_type == STP_RPCSPG || graph->stp_type == STP_MWCSP )
2149 {
2150 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
2151 if( !Is_term(graph->term[graph->head[e]]) && SCIPisEQ(scip, xval[e], 1.0) )
2152 break;
2153 if( e == EAT_LAST )
2154 return SCIP_OKAY;
2155 }
2156
2157 /* allocate memory */
2158 SCIP_CALL( SCIPallocBufferArray(scip, &results, nedges) );
2159 SCIP_CALL( SCIPallocBufferArray(scip, &nval, nvars) );
2160
2161 /* set solution array */
2162 for( e = 0; e < nedges; e++ )
2163 {
2164 if( SCIPisEQ(scip, xval[e], 1.0) )
2165 results[e] = CONNECT;
2166 else
2167 results[e] = UNKNOWN;
2168 }
2169
2170 if( !graph_sol_valid(scip, graph, results) )
2171 {
2172 SCIPfreeBufferArray(scip, &nval);
2173 SCIPfreeBufferArray(scip, &results);
2174 return SCIP_OKAY;
2175 }
2176
2177 /* pruning necessary? */
2178 if( SCIPsolGetHeur(newsol) == NULL ||
2179 !(strcmp(SCIPheurGetName(SCIPsolGetHeur(newsol)), "rec") == 0 ||
2180 strcmp(SCIPheurGetName(SCIPsolGetHeur(newsol)), "TM") == 0) )
2181 {
2182 const int nnodes = graph->knots;
2183 STP_Bool* steinertree;
2184 SCIP_CALL( SCIPallocBufferArray(scip, &steinertree, nnodes) );
2185 assert(graph_sol_valid(scip, graph, results));
2186
2187 for( v = nnodes - 1; v >= 0; --v )
2188 steinertree[v] = FALSE;
2189
2190 for( e = nedges - 1; e >= 0; --e )
2191 {
2192 if( results[e] == CONNECT )
2193 {
2194 steinertree[graph->tail[e]] = TRUE;
2195 steinertree[graph->head[e]] = TRUE;
2196 }
2197 results[e] = UNKNOWN;
2198 }
2199
2200 if( graph->stp_type == STP_PCSPG || graph->stp_type == STP_RPCSPG || graph->stp_type == STP_MWCSP )
2201 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, results, steinertree) );
2202 else
2203 SCIP_CALL( SCIPStpHeurTMPrune(scip, graph, graph->cost, 0, results, steinertree) );
2204
2205 SCIPfreeBufferArray(scip, &steinertree);
2206 }
2207
2208 /* execute local heuristics */
2209 SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, results) );
2210
2211 /* can we connect the network */
2212 for( v = 0; v < nvars; v++ )
2213 nval[v] = (results[v % nedges] == (v / nedges)) ? 1.0 : 0.0;
2214
2215 SCIP_CALL( SCIPStpValidateSol(scip, graph, nval, &feasible) );
2216
2217 /* solution feasible? */
2218 if( feasible )
2219 {
2220 assert(nedges == nvars);
2221
2222 pobj = 0.0;
2223
2224 for( v = 0; v < nedges; v++ )
2225 pobj += graph->cost[v] * nval[v];
2226
2227 /* has solution been improved? */
2228 if( SCIPisGT(scip, SCIPgetSolOrigObj(scip, newsol) - SCIPprobdataGetOffset(scip), pobj) )
2229 {
2230 SCIP_SOL* bestsol;
2231 SCIP_Bool success;
2232
2233 bestsol = sols[0];
2234 impsol = NULL;
2235 SCIP_CALL( SCIPprobdataAddNewSol(scip, nval, impsol, heur, &success) );
2236
2237 if( success )
2238 {
2239 *result = SCIP_FOUNDSOL;
2240
2241 if( heurdata->nbestsols < heurdata->maxnsols && SCIPisGT(scip, SCIPgetSolOrigObj(scip, bestsol) - SCIPprobdataGetOffset(scip), pobj) )
2242 {
2243 heurdata->nfails = 0;
2244 heurdata->nbestsols++;
2245 }
2246 SCIPdebugMessage("success in local: old: %f new: %f \n", (SCIPgetSolOrigObj(scip, bestsol) - SCIPprobdataGetOffset(scip)), pobj);
2247 }
2248 }
2249 }
2250
2251 if( *result != SCIP_FOUNDSOL )
2252 {
2253 heurdata->nfails++;
2254 if( heurdata->nbestsols > DEFAULT_MINNBESTSOLS && heurdata->nfails > 1 && graph->stp_type != STP_MWCSP )
2255 heurdata->nbestsols--;
2256
2257 SCIPdebugMessage("fail! %d \n", heurdata->nbestsols);
2258 }
2259
2260 SCIPfreeBufferArray(scip, &nval);
2261 SCIPfreeBufferArray(scip, &results);
2262
2263 return SCIP_OKAY;
2264 }
2265
2266 /*
2267 * primal heuristic specific interface methods
2268 */
2269
2270
2271 /** creates the local primal heuristic and includes it in SCIP */
SCIPStpIncludeHeurLocal(SCIP * scip)2272 SCIP_RETCODE SCIPStpIncludeHeurLocal(
2273 SCIP* scip /**< SCIP data structure */
2274 )
2275 {
2276 SCIP_HEURDATA* heurdata;
2277 SCIP_HEUR* heur;
2278
2279 /* create Local primal heuristic data */
2280 SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
2281
2282 /* include primal heuristic */
2283 SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
2284 HEUR_NAME, HEUR_DESC, HEUR_DISPCHAR, HEUR_PRIORITY, HEUR_FREQ, HEUR_FREQOFS,
2285 HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecLocal, heurdata) );
2286
2287 assert(heur != NULL);
2288
2289 /* set non-NULL pointers to callback methods */
2290 SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyLocal) );
2291 SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeLocal) );
2292 SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolLocal) );
2293 SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolLocal) );
2294
2295 /* add local primal heuristic parameters */
2296 SCIP_CALL( SCIPaddBoolParam(scip, "stp/duringroot",
2297 "should the heuristic be called during the root node?",
2298 &heurdata->duringroot, FALSE, DEFAULT_DURINGROOT, NULL, NULL) );
2299
2300 SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/"HEUR_NAME"/maxfreq",
2301 "should the heuristic be executed at maximum frequeny?",
2302 &heurdata->maxfreq, FALSE, DEFAULT_MAXFREQLOC, NULL, NULL) );
2303
2304 SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/maxnsols",
2305 "maximum number of best solutions to improve",
2306 &heurdata->maxnsols, FALSE, DEFAULT_MAXNBESTSOLS, 1, 50, NULL, NULL) );
2307
2308 return SCIP_OKAY;
2309 }
2310