Scippy

SCIP

Solving Constraint Integer Programs

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