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-2022 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 //#define SCIP_DEBUG
29 #include <assert.h>
30 #include <string.h>
31 #include "heur_local.h"
32 #include "heur_tm.h"
33 #include "probdata_stp.h"
34 #include "cons_stp.h"
35 #include "solstp.h"
36 #include "stpvector.h"
37 
38 
39 /* @note if heuristic is running in root node timing is changed there to (SCIP_HEURTIMING_DURINGLPLOOP |
40  * SCIP_HEURTIMING_BEFORENODE), see SCIP_DECL_HEURINITSOL callback
41  */
42 
43 #define HEUR_NAME "local"
44 #define HEUR_DESC "improvement heuristic for STP"
45 #define HEUR_DISPCHAR '-'
46 #define HEUR_PRIORITY 100
47 #define HEUR_FREQ 1
48 #define HEUR_FREQOFS 0
49 #define HEUR_MAXDEPTH -1
50 #define HEUR_TIMING (SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_DURINGLPLOOP | SCIP_HEURTIMING_AFTERLPLOOP | SCIP_HEURTIMING_AFTERNODE)
51 
52 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
53 
54 #define DEFAULT_DURINGROOT TRUE
55 #define DEFAULT_MAXFREQLOC FALSE
56 #define DEFAULT_MAXNBESTSOLS 50
57 #define DEFAULT_NBESTSOLS 15
58 #define DEFAULT_NBESTSOLS_HARD 50
59 #define DEFAULT_NELITESOLS 3
60 #define DEFAULT_MINNBESTSOLS 10
61 #define DEFAULT_MINNBESTSOLS_HARD 30
62 #define DEFAULT_RANDSEED 1492 /**< random seed */
63 #define LOCAL_MAXRESTARTS 10
64 #define LOCAL_MAXRESTARTS_FAST 3
65 
66 
67 #define GREEDY_MAXRESTARTS 3 /**< Max number of restarts for greedy PC/MW heuristic if improving solution has been found. */
68 #define GREEDY_EXTENSIONS_MW 6 /**< Number of extensions for greedy MW heuristic. MUST BE HIGHER THAN GREEDY_EXTENSIONS */
69 #define GREEDY_EXTENSIONS 5 /**< Number of extensions for greedy PC heuristic. */
70 
71 
72 /*
73  * Data structures
74  */
75 
76 /** primal heuristic data */
77 struct SCIP_HeurData
78 {
79  int nfails; /**< number of fails */
80  int maxnsols; /**< maximal number of best solutions to improve */
81  int nbestsols; /**< number of best solutions to improve */
82  int nbestsols_min; /**< minimum number of best solutions to improve */
83  int* lastsolindices; /**< indices of a number of best solutions already tried */
84  SCIP_Bool maxfreq; /**< should the heuristic be called with maximum frequency? */
85  SCIP_Bool duringroot; /**< should the heuristic be called during the root node? */
86  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
87 };
88 
89 /** Voronoi data for local heuristic */
91 {
92  PATH* vnoi_path; /**< path */
93  int* vnoi_base; /**< base*/
94  SCIP_Real* memvdist; /**< distance */
95  int* memvbase; /**< base*/
96  int* meminedges; /**< in-edge */
97  int* vnoi_nodestate; /**< node state */
98  int nmems; /**< number of memorized elements */
99  int nkpnodes; /**< number of key path nodes */
100 } VNOILOC;
101 
102 
103 /** Connectivity data */
104 typedef struct connectivity_data
105 {
106  STP_Vectype(int)* blists_start; /**< boundary lists starts (on nodes) */
107  STP_Vectype(int)* lvledges_start; /**< horizontal edges starts (on nodes) */
108  PHNODE** pheap_boundpaths; /**< boundary paths (on nodes) */
109  int* pheap_sizes; /**< size (on nodes) */
110  PHNODE* pheap_elements; /**< elements, one per edge! */
111  int* boundaryedges; /**< current boundary edge */
112  UF* uf; /**< union find */
113  int nboundaryedges; /**< number of boundary edges */
114 } CONN;
115 
116 
117 /** Key-paths data */
119 {
120  int* const kpnodes; /**< key path nodes */
121  int* const kpedges; /**< key path edges */
122  SCIP_Real kpcost; /**< cost of key paths */
123  int nkpnodes; /**< number of key path nodes */
124  int nkpedges; /**< number of key path edges */
125  int rootpathstart; /**< start of key path towards root component */
126  int kptailnode; /**< needed for single path */
127 } KPATHS;
128 
129 
130 /** Solution tree data */
131 typedef struct solution_tree_data
132 {
133  STP_Bool* const solNodes; /**< Steiner tree nodes */
134  LCNODE* const linkcutNodes; /**< Steiner tree nodes */
135  int* const solEdges; /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
136  STP_Bool* const nodeIsPinned; /**< of size nodes */
137  STP_Bool* const nodeIsScanned; /**< of size nodes */
138  int* const newedges; /**< marks new edges of the tree */
139 } SOLTREE;
140 
141 
142 /** Super graph data */
143 typedef struct supergraph_data
144 {
145  int* const supernodes; /**< super nodes */
146  STP_Bool* const nodeIsLowSupernode; /**< marks the current super-vertices
147  * (except for the one representing the root-component) */
148  PATH* mst; /**< MST */
149  SCIP_Real mstcost; /**< cost of MST */
150  int nsupernodes; /**< number of super nodes */
151 } SGRAPH;
152 
153 
154 /** Prize-collecting/maximum-weight connected subgraph data */
155 typedef struct pcmw_data
156 {
157  SCIP_Real* const prize_biased; /**< prize */
158  SCIP_Real* const edgecost_biased; /**< fully biased edge costs (PC only) */
159  SCIP_Real* const edgecost_org; /**< original, fully unbiased, edge costs (PC only) */
160  STP_Bool* const prizemark; /**< marked? */
161  int* const prizemarklist; /**< list of all marked */
162  int nprizemarks; /**< size of list */
163  SCIP_Bool isActive; /**< actively used (i.e. Pc/Mw graph)? */
164  int solroot; /**< root of Steiner tree solution */
165 } PCMW;
166 
167 
168 /** local insertion heuristic data */
169 typedef struct insertion_data
170 {
171  int* chainStarts; /**< pointers to starts of current chains (nInsertions many) */
172  int* chainEnds; /**< pointers to ends of current chains (nInsertions many) */
173  const SCIP_Real* edgecosts; /**< the edge costs (original for PC) */
174  int* const solDegreeNonTerm; /**< degree of node [v] in current solution; (pseudo) terminals are marked as UNKNOWN */
175  int* const addedEdges; /**< added edges */
176  int* const cutedgesStart; /**< cut edges for the chains */
177  int* const cutedgesEnd; /**< cut edges for the chains */
178  STP_Bool* const solNodes; /**< solution nodes */
179  SCIP_Bool* const nodeIsBlocked; /**< is node [v] blocked? */
180  int* const blockedList; /**< list of currently blocked nodes */
181  int blockedListSize; /**< size of list */
182  int nInsertions; /**< number of insertions */
183  int insertionVertex; /**< vertex to be inserted */
184 } INSERT;
185 
186 
187 /*
188  * Local methods
189  */
190 
191 
192 /** prune the solution? */
193 static inline
195  SCIP_SOL* initialsol /**< SCIP data structure */
196 )
197 {
198  assert(initialsol);
199 
200  if( SCIPsolGetHeur(initialsol) == NULL )
201  {
202  return TRUE;
203  }
204 
205  if( strcmp(SCIPheurGetName(SCIPsolGetHeur(initialsol)), "rec") == 0 )
206  {
207  return FALSE;
208  }
209 
210  if( strcmp(SCIPheurGetName(SCIPsolGetHeur(initialsol)), "TM") == 0 )
211  {
212  return FALSE;
213  }
214 
215  return TRUE;
216 }
217 
218 /** prune the solution? */
219 static inline
221  SCIP* scip, /**< SCIP data structure */
222  const GRAPH* graph, /**< graph data structure */
223  int* results /**< Steiner tree edges (in/out) */
224 )
225 {
226  STP_Bool* steinertree;
227  const int nnodes = graph->knots;
228 
229 #ifndef NDEBUG
230  const int nedges = graph->edges;
231  const SCIP_Real initialobj = solstp_getObjBounded(graph, results, 0.0, nedges);
232 #endif
233 
234  assert(solstp_isValid(scip, graph, results));
235 
236  SCIP_CALL( SCIPallocBufferArray(scip, &steinertree, nnodes) );
237 
238  solstp_setVertexFromEdge(graph, results, steinertree);
239 
240  SCIP_CALL( solstp_prune(scip, graph, results, steinertree) );
241 
242  SCIPfreeBufferArray(scip, &steinertree);
243 
244 #ifndef NDEBUG
245  {
246  const SCIP_Real initialobj_pruned = solstp_getObjBounded(graph, results, 0.0, nedges);
247 
248  assert(SCIPisLE(scip, initialobj_pruned, initialobj));
249  }
250 #endif
251 
252  return SCIP_OKAY;
253 }
254 
255 
256 /** fills solution array 'results' */
257 static inline
259  SCIP* scip, /**< SCIP data structure */
260  const GRAPH* graph, /**< graph data structure */
261  SCIP_SOL* initialsol, /**< SCIP data structure */
262  int* results /**< Steiner tree edges (out) */
263 )
264 {
265  SCIP_Real* xvals = SCIPprobdataGetXval(scip, initialsol);
266  const int nedges = graph_get_nEdges(graph);
267 
268  assert(results);
269  assert(xvals);
270 
271  /* set solution array */
272  for( int e = 0; e < nedges; e++ )
273  {
274  if( SCIPisEQ(scip, xvals[e], 1.0) )
275  {
276  results[e] = CONNECT;
277  }
278  else
279  {
280  assert(SCIPisEQ(scip, xvals[e], 0.0));
281 
282  results[e] = UNKNOWN;
283  }
284  }
285 }
286 
287 
288 /** initializes */
289 static
291  SCIP* scip, /**< SCIP data structure */
292  SCIP_HEURDATA* heurdata /**< heuristic's data */
293 )
294 {
295  assert(heurdata->nbestsols == -1 && heurdata->nbestsols_min == -1);
296 
298  {
299  heurdata->nbestsols = DEFAULT_NBESTSOLS_HARD;
300  heurdata->nbestsols_min = DEFAULT_MINNBESTSOLS_HARD;
301  }
302  else
303  {
304  heurdata->nbestsols = DEFAULT_NBESTSOLS;
305  heurdata->nbestsols_min = DEFAULT_MINNBESTSOLS;
306  }
307 }
308 
309 
310 /** tries to add solution stored in 'results' to SCIP solution pool */
311 static inline
313  SCIP* scip, /**< SCIP data structure */
314  SCIP_SOL** sols, /**< SCIP solutions */
315  SCIP_HEUR* heur, /**< heuristic */
316  const GRAPH* graph, /**< graph data structure */
317  SCIP_Real initialsol_obj, /**< objective */
318  SCIP_SOL* initialsol, /**< SCIP data structure */
319  int* results, /**< Steiner tree edges (out) */
320  SCIP_RESULT* result /**< pointer */
321 )
322 {
323  SCIP_Real* newsol_vals;
324  const int nedges = graph_get_nEdges(graph);
325  SCIP_Bool feasible = FALSE;
326  SCIP_HEURDATA* heurdata = SCIPheurGetData(heur);
327 
328  assert(heurdata);
329  assert(*result == SCIP_DIDNOTFIND);
330 
331  SCIP_CALL( SCIPallocBufferArray(scip, &newsol_vals, nedges) );
332 
333  for( int v = 0; v < nedges; v++ )
334  newsol_vals[v] = (results[v] == CONNECT) ? 1.0 : 0.0;
335 
336  SCIP_CALL( SCIPStpValidateSol(scip, graph, newsol_vals, FALSE, &feasible) );
337 
338  assert(feasible);
339 
340  /* is new solution feasible? */
341  if( feasible )
342  {
343  const SCIP_Real newsol_obj = solstp_getObjBounded(graph, results, 0.0, nedges);
344  const SCIP_Bool solIsBetter = SCIPisLT(scip, newsol_obj, initialsol_obj);
345  const SCIP_Bool solIsBetter_org
346  = SCIPisLT(scip, newsol_obj, SCIPgetSolOrigObj(scip, initialsol) - SCIPprobdataGetOffset(scip));
347 
348  assert(SCIPisLE(scip, newsol_obj, initialsol_obj));
349 
350 #ifndef NDEBUG
351  {
352  SCIP_Real pobj = 0.0;
353 
354  for( int v = 0; v < nedges; v++ )
355  pobj += graph->cost[v] * newsol_vals[v];
356 
357  assert(SCIPisEQ(scip, pobj, newsol_obj));
358  }
359 #endif
360 
361  // todo if solIsBetter || solIsBetter_org
362 
363  /* has solution been improved? */
364  if( solIsBetter || solIsBetter_org )
365  {
366  SCIP_SOL* const bestsol = sols[0];
367  SCIP_Bool success = FALSE;
368 
369  SCIP_CALL( SCIPprobdataAddNewSol(scip, newsol_vals, heur, &success) );
370 
371  if( success )
372  {
373  *result = SCIP_FOUNDSOL;
374 
375  if( heurdata->nbestsols < heurdata->maxnsols && SCIPisGT(scip, SCIPgetSolOrigObj(scip, bestsol) - SCIPprobdataGetOffset(scip), newsol_obj) )
376  {
377  heurdata->nfails = 0;
378  heurdata->nbestsols++;
379  }
380  SCIPdebugMessage("success in local: old: %f new: %f \n", (SCIPgetSolOrigObj(scip, bestsol) - SCIPprobdataGetOffset(scip)), newsol_obj);
381  }
382  }
383  }
384 
385  SCIPfreeBufferArray(scip, &newsol_vals);
386 
387  if( *result != SCIP_FOUNDSOL )
388  {
389  heurdata->nfails++;
390  if( heurdata->nbestsols > heurdata->nbestsols_min && heurdata->nfails > 1 )
391  heurdata->nbestsols--;
392 
393  SCIPdebugMessage("fail! %d \n", heurdata->nbestsols);
394  }
395 
396  return SCIP_OKAY;
397 }
398 
399 /** can the problem class be used? */
400 static inline
402  const GRAPH* graph /**< graph data structure */
403 )
404 {
405  const int stp_type = graph->stp_type;
406  if( stp_type != STP_SPG && stp_type != STP_RSMT && stp_type != STP_OARSMT && stp_type != STP_PCSPG
407  && stp_type != STP_RPCSPG && stp_type != STP_GSTP && stp_type != STP_MWCSP && stp_type != STP_RMWCSP )
408  return FALSE;
409 
410  return TRUE;
411 }
412 
413 /** submethod for local extend */
414 static
416  SCIP* scip, /**< SCIP data structure */
417  const GRAPH* graph, /**< graph data structure */
418  const PATH* path, /**< shortest data structure array */
419  int i, /**< node */
420  int greedyextensions, /**< greedyextensions */
421  int* nextensions, /**< nextensions */
422  GNODE* candidates, /**< candidates */
423  SCIP_PQUEUE* pqueue /**< pqueue */
424  )
425 {
426  assert(Is_pseudoTerm(graph->term[i]));
427  assert(!graph_pc_knotIsFixedTerm(graph, i));
428 
429  if( *nextensions < greedyextensions )
430  {
431  candidates[*nextensions].dist = graph->prize[i] - path[i].dist;
432  candidates[*nextensions].number = i;
433 
434  SCIP_CALL( SCIPpqueueInsert(pqueue, &(candidates[(*nextensions)++])) );
435  }
436  else
437  {
438  /* get candidate vertex of minimum value */
439  GNODE* min = (GNODE*) SCIPpqueueFirst(pqueue);
440  if( SCIPisLT(scip, min->dist, graph->prize[i] - path[i].dist) )
441  {
442  min = (GNODE*) SCIPpqueueRemove(pqueue);
443  min->dist = graph->prize[i] - path[i].dist;
444  min->number = i;
445  SCIP_CALL( SCIPpqueueInsert(pqueue, min) );
446  }
447  }
448 
449  return SCIP_OKAY;
450 }
451 
452 
453 /** perform local vertex insertion heuristic on given Steiner tree */
454 static
456  SCIP* scip, /**< SCIP data structure */
457  const GRAPH* graph, /**< graph data structure */
458  const int* solEdges, /**< Steiner tree edges */
459  LCNODE* linkcutNodes, /**< Steiner tree nodes */
460  STP_Bool* solNodes /**< Steiner tree nodes */
461  )
462 {
463  const int nnodes = graph->knots;
464  const int nedges = graph->edges;
465 
466  assert(solstp_isValid(scip, graph, solEdges));
467 
468  for( int i = 0; i < nnodes; i++ )
469  {
470  solNodes[i] = FALSE;
471  SCIPlinkcuttreeInitNode(&linkcutNodes[i]);
472  }
473 
474  /* create a link-cut tree representing the current Steiner tree */
475  for( int e = 0; e < nedges; e++ )
476  {
477  if( solEdges[e] == CONNECT )
478  SCIPlinkcuttreeLink(linkcutNodes, graph->head[e], graph->tail[e], flipedge(e));
479  }
480 
481  /* mark current Steiner tree nodes */
482  for( int e = 0; e < nedges; e++ )
483  {
484  if( solEdges[e] == CONNECT )
485  {
486  solNodes[graph->tail[e]] = TRUE;
487  solNodes[graph->head[e]] = TRUE;
488  }
489  }
490 
491  solNodes[graph->source] = TRUE;
492 
493  assert(linkcutNodes[graph->source].edge == -1);
494 }
495 
496 
497 /** is given Steiner tree a trivial solution (i.e. contains only one vertex?) */
498 static
500  const GRAPH* graph, /**< graph data structure */
501  const int* solEdges /**< Steiner tree edges */
502 )
503 {
504  const int root = graph->source;
505  SCIP_Bool isTrivial = TRUE;
506 
507  assert(graph_pc_isPcMw(graph));
508  assert(graph->extended);
509 
510  if( graph_pc_isRootedPcMw(graph) )
511  {
512  for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
513  {
514  if( solEdges[e] )
515  {
516  const int head = graph->head[e];
517  if( graph_pc_knotIsFixedTerm(graph, head) || !Is_term(graph->term[head]) )
518  {
519  isTrivial = FALSE;
520  break;
521  }
522  }
523  }
524  }
525  else
526  {
527  isTrivial = FALSE;
528  }
529 
530 
531  if( isTrivial )
532  {
533  SCIPdebugMessage("trivial solution given \n");
534  }
535 
536  return isTrivial;
537 }
538 
539 
540 /** computes lowest common ancestors for all pairs {vbase(v), vbase(u)} such that {u,w} is a boundary edge,
541  * first call should be with u := root */
542 static
544  SCIP* scip,
545  const GRAPH* graph,
546  int u,
547  UF* uf,
548  STP_Bool* nodesmark,
549  const int* steineredges,
550  STP_Vectype(int)* lcalists,
551  PHNODE** boundpaths,
552  const int* heapsize,
553  const int* vbase
554  )
555 {
556  int* uboundpaths; /* boundary-paths (each one represented by its boundary edge) having node 'u' as an endpoint */
557  uf->parent[u] = u;
558 
559  for( int oedge = graph->outbeg[u]; oedge != EAT_LAST; oedge = graph->oeat[oedge] )
560  {
561  const int v = graph->head[oedge];
562  if( steineredges[oedge] == CONNECT )
563  {
564  SCIP_CALL( lca(scip, graph, v, uf, nodesmark, steineredges, lcalists, boundpaths, heapsize, vbase) );
565  SCIPStpunionfindUnion(uf, u, v, FALSE);
566  uf->parent[SCIPStpunionfindFind(uf, u)] = u;
567  }
568  }
569 
570  nodesmark[u] = TRUE;
571 
572  /* iterate through all boundary-paths having one endpoint in the voronoi region of node 'u' */
573  SCIP_CALL( SCIPpairheapBuffarr(scip, boundpaths[u], heapsize[u], &uboundpaths) );
574 
575  for( int i = 0; i < heapsize[u]; i++ )
576  {
577  const int oedge = uboundpaths[i];
578  const int v = vbase[graph->head[oedge]];
579  if( nodesmark[v] )
580  {
581  const int ancestor = uf->parent[SCIPStpunionfindFind(uf, v)];
582 
583  /* if the ancestor of 'u' and 'v' is one of the two, the boundary-edge is already in boundpaths[u] */
584  if( ancestor != u && ancestor != v)
585  {
586  StpVecPushBack(scip, lcalists[ancestor], oedge);
587  }
588  }
589  }
590 
591  /* free the boundary-paths array */
592  SCIPfreeBufferArrayNull(scip, &uboundpaths);
593 
594  return SCIP_OKAY;
595 }
596 
597 
598 /** computes lowest common ancestors for all pairs {vbase(v), vbase(u)} such that {u,w} is a boundary edge */
599 static
601  SCIP* scip, /**< SCIP data structure */
602  const GRAPH* graph, /**< graph data structure */
603  const VNOILOC* vnoiData, /**< Voronoi data */
604  const SOLTREE* soltreeData, /**< solution tree data */
605  CONN* connectData /**< data */
606 )
607 {
608  PHNODE** const boundpaths = connectData->pheap_boundpaths;
609  STP_Vectype(int)* lvledges_start = connectData->lvledges_start;
610  const int* const solEdges = soltreeData->solEdges;
611  int* const pheapsize = connectData->pheap_sizes;
612  const int* const vnoibase = vnoiData->vnoi_base;
613  STP_Bool* nodesmark;
614  const int nnodes = graph->knots;
615 
616  assert(SCIPStpunionfindIsClear(scip, connectData->uf));
617 
618  SCIP_CALL( SCIPallocBufferArray(scip, &nodesmark, nnodes) );
619 
620  for( int i = 0; i < nnodes; ++i )
621  nodesmark[i] = FALSE;
622 
623  SCIP_CALL( lca(scip, graph, graph->source, connectData->uf, nodesmark, solEdges, lvledges_start, boundpaths, pheapsize, vnoibase) );
624 
625  SCIPfreeBufferArray(scip, &nodesmark);
626 
627  return SCIP_OKAY;
628 }
629 
630 
631 /** recursive method for a DFS ordering of graph 'g' */
632 static
633 void dfsorder(
634  const GRAPH* graph,
635  const int* edges,
636  int node,
637  int* counter,
638  int* dfst
639  )
640 {
641  int oedge = graph->outbeg[node];
642 
643  while( oedge >= 0 )
644  {
645  if( edges[oedge] >= 0 )
646  {
647  dfsorder(graph, edges, graph->head[oedge], counter, dfst);
648  }
649  oedge = graph->oeat[oedge];
650  }
651 
652  dfst[*counter] = node;
653  (*counter)++;
654 }
655 
656 
657 /** helper method for pcmwGetNewEdgePrize */
658 static inline
660  const GRAPH* graph,
661  const STP_Bool* steinertree,
662  const SCIP_Real* prize_biased,
663  const int* graphmark,
664  int node,
665  STP_Bool* prizemark,
666  int* prizemarklist,
667  int* prizemarkcount
668  )
669 {
670  SCIP_Real prizesum = 0.0;
671  assert(graph_pc_isPcMw(graph));
672 
673  if( graphmark[node] && !steinertree[node] && Is_pseudoTerm(graph->term[node]) && !prizemark[node] )
674  {
675  prizesum += prize_biased[node];
676  prizemark[node] = TRUE;
677  prizemarklist[(*prizemarkcount)++] = node;
678  }
679 
680  return prizesum;
681 }
682 
683 
684 /** gets root of solution for Pc/Mw or -1 otherwise */
685 static
687  const GRAPH* graph,
688  const int* solEdges
689 )
690 {
691  int solroot = -1;
692 
693  if( graph_pc_isPcMw(graph) && !graph_pc_isRootedPcMw(graph) )
694  {
695  const int root = graph->source;
696 
697  for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
698  {
699  const int head = graph->head[e];
700  if( solEdges[e] == CONNECT && !Is_term(graph->term[head]) )
701  {
702  assert(graph->cost[e] == 0.0);
703  solroot = head;
704  break;
705  }
706  }
707  assert(solroot >= 0);
708  }
709 
710  return solroot;
711 }
712 
713 
714 /** get prize not already counted that is associated to edge, not including solution nodes or forbidden ones */
715 static
717  const GRAPH* graph, /**< graph */
718  const STP_Bool* steinertree, /**< Steiner tree */
719  int edge, /**< the edge */
720  PCMW* pcmwData /**< data */
721  )
722 {
723  SCIP_Real prizesum = 0.0;
724 
725  assert(graph && steinertree && pcmwData);
726 
727  if( pcmwData->isActive )
728  {
729  const int mhead = graph->head[edge];
730  const int mtail = graph->tail[edge];
731  STP_Bool* const prizemark = pcmwData->prizemark;
732  int* const prizemarklist = pcmwData->prizemarklist;
733 
734  assert(pcmwData->prize_biased);
735  assert(graph_pc_isPcMw(graph));
736 
737  prizesum += getNewPrizeNode(graph, steinertree, pcmwData->prize_biased, graph->mark, mhead,
738  prizemark, prizemarklist, &(pcmwData->nprizemarks));
739  prizesum += getNewPrizeNode(graph, steinertree, pcmwData->prize_biased, graph->mark, mtail,
740  prizemark, prizemarklist, &(pcmwData->nprizemarks));
741  }
742 
743  return prizesum;
744 }
745 
746 
747 #ifndef NDEBUG
748 /** is the data clean? */
749 static
751  const GRAPH* graph, /**< graph data structure */
752  const PCMW* pcmwData /**< data */
753 )
754 {
755  if( !pcmwData->isActive )
756  return TRUE;
757 
758  if( pcmwData->nprizemarks != 0 )
759  {
760  SCIPdebugMessage("nprizemarks != 0 \n");
761  return FALSE;
762  }
763 
764  for( int k = 0; k < graph->knots; k++ )
765  {
766  if( pcmwData->prizemark[k] )
767  {
768  SCIPdebugMessage("entry %d is dirty \n", k);
769  return FALSE;
770  }
771  }
772 
773  return TRUE;
774 }
775 #endif
776 
777 
778 /** clean the data */
779 static
781  const GRAPH* graph, /**< graph data structure */
782  PCMW* pcmwData /**< data */
783 )
784 {
785  if( pcmwData->isActive )
786  {
787  const int* const prizemarklist = pcmwData->prizemarklist;
788  STP_Bool* const prizemark = pcmwData->prizemark;
789  const int prizemarkcount = pcmwData->nprizemarks;
790 
791 
792  for( int pi = 0; pi < prizemarkcount; pi++ )
793  prizemark[prizemarklist[pi]] = FALSE;
794 
795  pcmwData->nprizemarks = 0;
796 
797  assert(pcmwDataIsClean(graph, pcmwData));
798  }
799 }
800 
801 
802 /** updates for PC/MW: graph->mark and pinned */
803 static
805  SCIP* scip, /**< SCIP data structure */
806  GRAPH* graph, /**< graph data structure */
807  SOLTREE* soltreeData, /**< solution tree data */
808  PCMW* pcmwData /**< data */
809  )
810 {
811  SCIP_Bool* reachable_nodes;
812  STP_Bool* const pinned = soltreeData->nodeIsPinned;
813  const int root = graph->source;
814  const int nnodes = graph->knots;
815  int* const graphmark = graph->mark;
816  const int solstartnode = graph_pc_isRootedPcMw(graph)? graph->source : pcmwData->solroot;
817 
818  assert(graph->extended);
819  assert(graph_pc_isPcMw(graph));
820  assert(solstartnode >= 0 && solstartnode < graph->knots);
821  assert(soltreeData->solNodes[solstartnode]);
822 
823  /* remove unconnected vertices */
824  SCIP_CALL( SCIPallocBufferArray(scip, &reachable_nodes, nnodes) );
825 
826  SCIP_CALL( graph_trail_costAware(scip, graph, solstartnode, reachable_nodes) );
827 
828  for( int k = 0; k < nnodes; k++ )
829  graphmark[k] = reachable_nodes[k];
830 
831  SCIPfreeBufferArray(scip, &reachable_nodes);
832 
833  /* unmark and pin pseudo terminals */
834  for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
835  {
836  const int k = graph->head[e];
837  if( Is_term(graph->term[k]) )
838  {
839  if( !graph_pc_knotIsFixedTerm(graph, k) )
840  {
841  const int pterm = graph->head[graph->term2edge[k]];
842 
843  assert(Is_pseudoTerm(graph->term[pterm]));
844 
845  graphmark[k] = FALSE;
846  pinned[pterm] = TRUE;
847  }
848  }
849  }
850 
851  if( !graph_pc_isRootedPcMw(graph) )
852  {
853  graphmark[root] = FALSE;
854  soltreeData->solNodes[root] = FALSE;
855  }
856 
857  return SCIP_OKAY;
858 }
859 
860 
861 /** checks whether node is crucial, i.e. a terminal or a vertex with degree at least 3 (w.r.t. the Steiner tree) */
862 static
864  const GRAPH* graph, /**< graph data structure */
865  const int* steineredges, /**< Steiner edges */
866  int node /**< node to check */
867 )
868 {
869  assert(graph != NULL);
870  assert(steineredges != NULL);
871  assert(node >= 0 && node < graph->knots);
872 
873  if( !Is_term(graph->term[node]) )
874  {
875  int counter = 0;
876 
877  for( int e = graph->outbeg[node]; e != EAT_LAST; e = graph->oeat[e] )
878  {
879  /* check if the adjacent node is in the Steiner tree */
880  if( steineredges[e] == CONNECT || steineredges[flipedge(e)] == CONNECT )
881  counter++;
882  }
883 
884  /* is node of degree smaller 3 in the Steiner tree? */
885  if( counter < 3 )
886  {
887  return FALSE;
888  }
889  }
890 
891  return TRUE;
892 }
893 
894 
895 /** gets unbiased edge cost */
896 static inline
898  const GRAPH* graph, /**< graph data structure */
899  const PCMW* pcmwData, /**< data */
900  int edge /**< edge */
901  )
902 {
903  SCIP_Real cost;
904 
905  assert(graph && pcmwData);
906  assert(edge >= 0 && edge < graph->edges);
907 
908  if( pcmwData->isActive )
909  {
910  assert(graph->extended);
911 
912  if( !graph_pc_isPc(graph) )
913  {
914  assert(graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP);
915 
916  cost = 0.0;
917  }
918  else
919  {
920  assert(pcmwData->edgecost_org);
921 
922  cost = pcmwData->edgecost_org[edge];
923  }
924  }
925  else
926  {
927  cost = graph->cost[edge];
928  }
929 
930  return cost;
931 }
932 
933 
934 /** Get cost of shortest path along boundary edge.
935  * For Pc/Mw: Will consider biased edge costs, but not the actual prizes of proper terminals! */
936 static
938  const GRAPH* graph, /**< graph data structure */
939  const VNOILOC* vnoiData, /**< data */
940  const PCMW* pcmwData, /**< data */
941  int boundaryedge /**< boundary edge*/
942  )
943 {
944  const PATH* const vnoipath = vnoiData->vnoi_path;
945  SCIP_Real pathcost;
946  const int tail = graph->tail[boundaryedge];
947  const int head = graph->head[boundaryedge];
948 
949  assert(boundaryedge >= 0);
950  assert(vnoiData->vnoi_base[tail] != vnoiData->vnoi_base[head]);
951  assert(graph_pc_isPcMw(graph) == pcmwData->isActive);
952 
953  pathcost = vnoipath[tail].dist + vnoipath[head].dist;
954 
955  if( pcmwData->isActive )
956  {
957  /* either vnoipath[head] goes into head, or head is the end of the entire boundary path;
958  * either way we don't want to subtract the prize! */
959  pathcost += getEdgeCostUnbiased(graph, pcmwData, boundaryedge);
960  }
961  else
962  {
963  pathcost += graph->cost[boundaryedge];
964  }
965 
966  assert(pathcost >= 0.0);
967 
968  return pathcost;
969 }
970 
971 
972 /** For Pc/Mw: Consider biased edge costs, and the actual prizes of proper terminals! */
973 static
975  const GRAPH* graph, /**< graph data structure */
976  const VNOILOC* vnoiData, /**< data */
977  const SOLTREE* soltreeData, /**< solution tree data */
978  int boundaryedge, /**< boundary edge */
979  PCMW* pcmwData /**< data */
980  )
981 {
982  SCIP_Real pathcost;
983 
984  assert(boundaryedge >= 0);
985  assert(vnoiData->vnoi_base[graph->tail[boundaryedge]] != vnoiData->vnoi_base[graph->head[boundaryedge]]);
986  assert(graph_pc_isPcMw(graph) == pcmwData->isActive);
987 
988  pathcost = getBoundaryPathCost(graph, vnoiData, pcmwData, boundaryedge);
989 
990  if( pcmwData->isActive )
991  {
992  const STP_Bool* const solNodes = soltreeData->solNodes;
993  const PATH* const vnoipath = vnoiData->vnoi_path;
994  const int* const vnoibase = vnoiData->vnoi_base;
995 
996 #ifdef SCIP_DEBUG
997  printf("boundary path edge ");
998  graph_edge_printInfo(graph, boundaryedge);
999 #endif
1000  pathcost -= pcmwGetNewEdgePrize(graph, solNodes, boundaryedge, pcmwData);
1001 
1002  for( int node = graph->tail[boundaryedge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1003  {
1004 #ifdef SCIP_DEBUG
1005  printf("boundary path edge ");
1006  graph_edge_printInfo(graph, vnoipath[node].edge);
1007 #endif
1008  pathcost -= pcmwGetNewEdgePrize(graph, solNodes, vnoipath[node].edge, pcmwData);
1009  }
1010 
1011  for( int node = graph->head[boundaryedge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1012  {
1013 #ifdef SCIP_DEBUG
1014  printf("boundary path edge ");
1015  graph_edge_printInfo(graph, vnoipath[node].edge);
1016 #endif
1017 
1018  pathcost -= pcmwGetNewEdgePrize(graph, solNodes, vnoipath[node].edge, pcmwData);
1019  }
1020 
1021  if( pathcost < 0.0 )
1022  pathcost = 0.0;
1023  }
1024 
1025  assert(pathcost >= 0.0);
1026 
1027  return pathcost;
1028 }
1029 
1030 
1031 /** returns edge costs (biased for Pc/Mw) */
1032 static inline
1034  const GRAPH* graph, /**< graph data structure */
1035  const PCMW* pcmwData /**< data */
1036 )
1037 {
1038  assert( graph_pc_isPcMw(graph) == pcmwData->isActive);
1039 
1040  return (pcmwData->isActive? pcmwData->edgecost_biased : graph->cost);
1041 }
1042 
1043 
1044 /** update for key-vertex elimination: Save current boundary edges */
1045 static
1047  SCIP* scip, /**< SCIP data structure */
1048  const GRAPH* graph, /**< graph data structure */
1049  const VNOILOC* vnoiData, /**< Voronoi data */
1050  const SGRAPH* supergraphData, /**< super-graph*/
1051  int crucnode, /**< node to eliminate */
1052  CONN* connectData /**< data */
1053 )
1054 {
1055  PHNODE** const boundpaths = connectData->pheap_boundpaths;
1056  STP_Vectype(int)* lvledges_start = connectData->lvledges_start;
1057  STP_Vectype(int) lvledges_curr = NULL;
1058  UF* const uf = connectData->uf;
1059  int* const pheapsize = connectData->pheap_sizes;
1060  const int* const supernodes = supergraphData->supernodes;
1061  const int* const vnoibase = vnoiData->vnoi_base;
1062  const STP_Bool* const isLowSupernode = supergraphData->nodeIsLowSupernode;
1063  int* const boundaryedges = connectData->boundaryedges;
1064  const int* const graphmark = graph->mark;
1065  int nboundaryedges = 0;
1066 
1067  connectData->nboundaryedges = -1;
1068 
1069  /* add vertical boundary-paths between the child components and the root-component (w.r.t. node 'crucnode') */
1070  for( int k = 0; k < supergraphData->nsupernodes - 1; k++ )
1071  {
1072  const int supernode = supernodes[k];
1073  int edge = UNKNOWN;
1074 
1075  while( boundpaths[supernode] != NULL )
1076  {
1077  int node;
1078  SCIP_Real edgecost;
1079 
1080  SCIP_CALL( SCIPpairheapDeletemin(scip, &edge, &edgecost, &boundpaths[supernode], &pheapsize[supernode]) );
1081 
1082  node = (vnoibase[graph->head[edge]] == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(uf, vnoibase[graph->head[edge]]);
1083 
1084  /* check whether edge 'edge' represents a boundary-path having an endpoint in the kth-component and in the root-component respectively */
1085  if( node != UNKNOWN && !isLowSupernode[node] && graphmark[node] )
1086  {
1087  boundaryedges[nboundaryedges++] = edge;
1088  SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[supernode],
1089  connectData->pheap_elements, edge, edgecost, &pheapsize[supernode]) );
1090  break;
1091  }
1092  }
1093  }
1094 
1095  lvledges_curr = lvledges_start[crucnode];
1096 
1097  /* add horizontal boundary-paths (between the child-components) */
1098  for( int i = 0; i < StpVecGetSize(lvledges_curr); i++ )
1099  {
1100  const int edge = lvledges_curr[i];
1101  const int basetail = vnoibase[graph->tail[edge]];
1102  const int basehead = vnoibase[graph->head[edge]];
1103  const int node = (basehead == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(uf, basehead);
1104  const int adjnode = (basetail == UNKNOWN)? UNKNOWN : SCIPStpunionfindFind(uf, basetail);
1105 
1106  /* check whether the current boundary-path connects two child components */
1107  if( node != UNKNOWN && isLowSupernode[node] && adjnode != UNKNOWN && isLowSupernode[adjnode] )
1108  {
1109  assert(graphmark[node] && graphmark[adjnode]);
1110  boundaryedges[nboundaryedges++] = edge;
1111  }
1112  }
1113 
1114  connectData->nboundaryedges = nboundaryedges;
1115 
1116  return SCIP_OKAY;
1117 }
1118 
1119 
1120 /** initialize */
1121 static
1123  SCIP* scip, /**< SCIP data structure */
1124  const GRAPH* graph, /**< graph data structure */
1125  const VNOILOC* vnoiData, /**< Voronoi data */
1126  const SOLTREE* soltreeData, /**< solution tree data */
1127  const PCMW* pcmwData, /**< data */
1128  CONN* connectData /**< data */
1129 )
1130 {
1131  PHNODE** const boundpaths = connectData->pheap_boundpaths;
1132  STP_Vectype(int)* blists_start = connectData->blists_start;
1133  STP_Vectype(int)* lvledges_start = connectData->lvledges_start;
1134  int* const pheapsize = connectData->pheap_sizes;
1135  const int* const vnoibase = vnoiData->vnoi_base;
1136  const int* const graphmark = graph->mark;
1137  const int nnodes = graph->knots;
1138  const int nedges = graph->edges;
1139 
1140  assert(connectData->nboundaryedges == 0);
1141  assert(connectData->boundaryedges);
1142 
1143  for( int k = 0; k < nnodes; k++ )
1144  assert(graph->path_state[k] == CONNECT || !graph->mark[k]);
1145 
1146  for( int k = 0; k < nnodes; ++k )
1147  {
1148  if( blists_start[k] )
1149  StpVecClear(blists_start[k]);
1150  if( lvledges_start[k] )
1151  StpVecClear(lvledges_start[k]);
1152  }
1153 
1154  BMSclearMemoryArray(pheapsize, nnodes);
1155  BMSclearMemoryArray(boundpaths, nnodes);
1156 
1157  /* link all nodes to their (respective) Voronoi base */
1158  for( int k = 0; k < nnodes; ++k )
1159  {
1160  assert(pheapsize[k] == 0 && boundpaths[k] == NULL);
1161 
1162  if( !graphmark[k] )
1163  continue;
1164 
1165  StpVecPushBack(scip, blists_start[vnoibase[k]], k);
1166  }
1167 
1168  /* for each node, store all of its outgoing boundary-edges in a (respective) heap*/
1169  for( int e = 0; e < nedges; e += 2 )
1170  {
1171  if( graph->oeat[e] != EAT_FREE )
1172  {
1173  const int tail = graph->tail[e];
1174  const int head = graph->head[e];
1175 
1176  /* is edge 'e' a boundary-edge? */
1177  if( vnoibase[tail] != vnoibase[head] && graphmark[tail] && graphmark[head] )
1178  {
1179  const SCIP_Real edgecost = getBoundaryPathCost(graph, vnoiData, pcmwData, e);
1180 
1181  assert(vnoibase[tail] != UNKNOWN && vnoibase[head] != UNKNOWN);
1182  assert(SCIPisGE(scip, edgecost, 0.0));
1183 
1184  /* add the boundary-edge 'e' and its reversed to the corresponding heaps */
1185  SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[vnoibase[tail]],
1186  connectData->pheap_elements, e, edgecost, &(pheapsize[vnoibase[tail]])) );
1187  SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[vnoibase[head]],
1188  connectData->pheap_elements, flipedge(e), edgecost, &(pheapsize[vnoibase[head]])) );
1189  }
1190  }
1191  }
1192 
1193  SCIP_CALL( getLowestCommonAncestors(scip, graph, vnoiData, soltreeData, connectData) );
1194 
1195  return SCIP_OKAY;
1196 }
1197 
1198 
1199 /** get key path above given crucial node */
1200 static
1202  SCIP* scip, /**< SCIP data structure */
1203  int crucnode, /**< crucial node to start from */
1204  const GRAPH* graph, /**< graph data structure */
1205  const SOLTREE* soltreeData, /**< solution tree data */
1206  const PCMW* pcmwData, /**< data */
1207  CONN* connectData, /**< data */
1208  KPATHS* keypathsData, /**< key paths */
1209  SCIP_Bool* allGood /**< all good? */
1210 )
1211 {
1212  int* const kpnodes = keypathsData->kpnodes;
1213  const int* const solEdges = soltreeData->solEdges;
1214  const LCNODE* const linkcutNodes = soltreeData->linkcutNodes;
1215  const STP_Bool* const solNodes = soltreeData->solNodes;
1216  const STP_Bool* const pinned = soltreeData->nodeIsPinned;
1217  PHNODE** const boundpaths = connectData->pheap_boundpaths;
1218  const SCIP_Real* const edgecosts = getEdgeCosts(graph, pcmwData);
1219  int* const pheapsize = connectData->pheap_sizes;
1220  const int* const graphmark = graph->mark;
1221  int nkpnodes = 0;
1222  int firstedge2root = -1;
1223  int kptailnode = -1;
1224  SCIP_Real kpcost = -FARAWAY;
1225 
1226  assert(*allGood);
1227 
1228 
1229  if( Is_term(graph->term[crucnode]) || pinned[crucnode] )
1230  {
1231  /* update the data structures; basically merge the whole key path */
1232  for( int edge = graph->outbeg[crucnode]; edge != EAT_LAST; edge = graph->oeat[edge] )
1233  {
1234  int adjnode = graph->head[edge];
1235 
1236  /* check whether edge 'edge' leads to an ancestor of terminal 'crucnode' */
1237  if( solEdges[edge] == CONNECT && solNodes[adjnode] && graphmark[adjnode] )
1238  {
1239  assert( SCIPStpunionfindFind(connectData->uf, adjnode) != crucnode);
1240  assert(soltreeData->nodeIsScanned[adjnode]);
1241 
1242  SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[adjnode], &pheapsize[crucnode], &pheapsize[adjnode]);
1243 
1244  /* update the union-find data structure */
1245  SCIPStpunionfindUnion(connectData->uf, crucnode, adjnode, FALSE);
1246 
1247  /* move along the key-path until its end (i.e. until a crucial node is reached) */
1248  while( !nodeIsCrucial(graph, solEdges, adjnode) && !pinned[adjnode] )
1249  {
1250  int e;
1251  for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
1252  {
1253  if( solEdges[e] == CONNECT )
1254  break;
1255 
1256  assert(UNKNOWN == solEdges[e]);
1257  }
1258 
1259  /* leaf of the Steiner tree is not a terminal? NOTE: should happen very rarely */
1260  if( e == EAT_LAST )
1261  {
1262  *allGood = FALSE;
1263  return;
1264  }
1265 
1266  adjnode = graph->head[e];
1267 
1268  if( !solNodes[adjnode] || !graphmark[adjnode] )
1269  break;
1270 
1271  assert(soltreeData->nodeIsScanned[adjnode]);
1272  assert(SCIPStpunionfindFind(connectData->uf, adjnode) != crucnode);
1273 
1274  /* update the union-find data structure */
1275  SCIPStpunionfindUnion(connectData->uf, crucnode, adjnode, FALSE);
1276 
1277  SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[adjnode], &pheapsize[crucnode], &pheapsize[adjnode]);
1278  }
1279  }
1280  }
1281  }
1282 
1283 #ifndef NDEBUG
1284  if( SCIPisGE(scip, graph->cost[linkcutNodes[crucnode].edge], FARAWAY)
1285  || SCIPisGE(scip, graph->cost[flipedge(linkcutNodes[crucnode].edge)], FARAWAY) )
1286  {
1287  assert(graph_pc_isPcMw(graph));
1288  assert(graph->head[linkcutNodes[crucnode].edge] == graph->source);
1289  }
1290 #endif
1291 
1292  /* traverse the (unique) key-path containing the parent of the current crucial node 'crucnode' */
1293 
1294  firstedge2root = linkcutNodes[crucnode].edge;
1295  kptailnode = graph->head[firstedge2root];
1296 
1297  if( pcmwData->isActive )
1298  {
1299  /* we don't want to count the prize of 'crucnode' */
1300  kpcost = getEdgeCostUnbiased(graph, pcmwData, flipedge(firstedge2root));
1301  }
1302  else
1303  {
1304  kpcost = edgecosts[firstedge2root];
1305  }
1306 
1307  assert(graph->tail[firstedge2root] == crucnode);
1308  assert(soltreeData->solEdges[flipedge(firstedge2root)] == CONNECT);
1309 
1310 #ifdef SCIP_DEBUG
1311  printf("key path edge ");
1312  graph_edge_printInfo(graph, firstedge2root);
1313 #endif
1314 
1315  while( !nodeIsCrucial(graph, solEdges, kptailnode) && !pinned[kptailnode] )
1316  {
1317  const int kpedge2root = linkcutNodes[kptailnode].edge;
1318  kpcost += edgecosts[flipedge(kpedge2root)];
1319 
1320 #ifdef SCIP_DEBUG
1321  printf("key path edge ");
1322  graph_edge_printInfo(graph, kpedge2root);
1323 #endif
1324 
1325  kpnodes[nkpnodes++] = kptailnode;
1326  kptailnode = graph->head[kpedge2root];
1327  }
1328 
1329  keypathsData->kpcost = kpcost;
1330  keypathsData->kptailnode = kptailnode;
1331  keypathsData->nkpnodes = nkpnodes;
1332 }
1333 
1334 
1335 /** unmarks key path nodes */
1336 static
1338  const GRAPH* graph, /**< graph data structure */
1339  const KPATHS* keypathsData, /**< key paths */
1340  SOLTREE* soltreeData /**< solution tree data */
1341 )
1342 {
1343  const int* const kpnodes = keypathsData->kpnodes;
1344  STP_Bool* const solNodes = soltreeData->solNodes;
1345  const int nkpnodes = keypathsData->nkpnodes;
1346 
1347  for( int k = 0; k < nkpnodes; k++ )
1348  {
1349  const int node = kpnodes[k];
1350 
1351  assert(node >= 0 && node < graph->knots);
1352  assert(solNodes[node]);
1353 
1354  solNodes[node] = FALSE;
1355  }
1356 }
1357 
1358 
1359 /** (re-)marks key path nodes */
1360 static
1362  const GRAPH* graph, /**< graph data structure */
1363  const KPATHS* keypathsData, /**< key paths */
1364  SOLTREE* soltreeData /**< solution tree data */
1365 )
1366 {
1367  const int* const kpnodes = keypathsData->kpnodes;
1368  STP_Bool* const solNodes = soltreeData->solNodes;
1369  const int nkpnodes = keypathsData->nkpnodes;
1370 
1371  for( int k = 0; k < nkpnodes; k++ )
1372  {
1373  const int node = kpnodes[k];
1374 
1375  assert(node >= 0 && node < graph->knots);
1376  assert(!solNodes[node]);
1377 
1378  solNodes[node] = TRUE;
1379  }
1380 }
1381 
1382 
1383 /** exchanges key path */
1384 static
1386  SCIP* scip, /**< SCIP data structure */
1387  GRAPH* graph, /**< graph data structure */
1388  const CONN* connectData, /**< data */
1389  const VNOILOC* vnoiData, /**< data */
1390  const KPATHS* keypathsData, /**< key paths */
1391  const int* dfstree, /**< DFS tree */
1392  const STP_Bool* scanned, /**< array to mark which nodes have been scanned */
1393  int dfstree_pos, /**< current position in dfs tree */
1394  int boundedge_new, /**< new Voronoi boundary edge */
1395  SOLTREE* soltreeData /**< solution tree data */
1396 )
1397 {
1398  UF* const uf = connectData->uf;
1399  PHNODE** const boundpaths = connectData->pheap_boundpaths;
1400  int* const pheapsize = connectData->pheap_sizes;
1401  const PATH* const vnoipath = vnoiData->vnoi_path;
1402  const int* const vnoibase = vnoiData->vnoi_base;
1403  const int* const kpnodes = keypathsData->kpnodes;
1404  STP_Bool* const pinned = soltreeData->nodeIsPinned;
1405  const LCNODE* const linkcutNodes = soltreeData->linkcutNodes;
1406  int* const solEdges = soltreeData->solEdges;
1407  STP_Bool* const solNodes = soltreeData->solNodes;
1408  const int nkpnodes = keypathsData->nkpnodes;
1409  const int crucnode = dfstree[dfstree_pos];
1410  int* const graphmark = graph->mark;
1411  int newpathend = -1;
1412  int newedge = boundedge_new;
1413  int node = SCIPStpunionfindFind(uf, vnoibase[graph->head[newedge]]);
1414 
1415  /* remove old keypath */
1416  assert(solEdges[flipedge(linkcutNodes[crucnode].edge)] != UNKNOWN);
1417 
1418  solEdges[flipedge(linkcutNodes[crucnode].edge)] = UNKNOWN;
1419  solNodes[crucnode] = FALSE;
1420  graphmark[crucnode] = FALSE;
1421 
1422  for( int k = 0; k < nkpnodes; k++ )
1423  {
1424  const int keypathnode = kpnodes[k];
1425  assert(solEdges[flipedge(linkcutNodes[keypathnode].edge)] != UNKNOWN);
1426 
1427  solEdges[flipedge(linkcutNodes[keypathnode].edge)] = UNKNOWN;
1428  solNodes[keypathnode] = FALSE;
1429  graphmark[keypathnode] = FALSE;
1430  }
1431 
1432  assert(graphmark[keypathsData->kptailnode]);
1433 
1434  if( node == crucnode )
1435  newedge = flipedge(newedge);
1436 
1437  for( node = graph->tail[newedge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1438  {
1439  graphmark[node] = FALSE;
1440 
1441  solEdges[flipedge(vnoipath[node].edge)] = CONNECT;
1442  solEdges[vnoipath[node].edge] = UNKNOWN;
1443  }
1444 
1445  for( node = graph->head[newedge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1446  {
1447  graphmark[node] = FALSE;
1448 
1449  solEdges[vnoipath[node].edge] = CONNECT;
1450  }
1451 
1452  solEdges[flipedge(newedge)] = CONNECT;
1453 
1454  newpathend = vnoibase[graph->tail[newedge]];
1455  assert(node == vnoibase[graph->head[newedge]] );
1456 
1457  /* flip all edges on the ST path between the endnode of the new key-path and the current crucial node */
1458  assert(SCIPStpunionfindFind(uf, newpathend) == crucnode);
1459 
1460  for( int k = newpathend; k != crucnode; k = graph->head[linkcutNodes[k].edge] )
1461  {
1462  assert(graphmark[k]);
1463  assert( solEdges[flipedge(linkcutNodes[k].edge)] != -1);
1464 
1465  solEdges[flipedge(linkcutNodes[k].edge)] = UNKNOWN;
1466  solEdges[linkcutNodes[k].edge] = CONNECT;
1467  }
1468 
1469  for( int k = 0; k < dfstree_pos; k++ )
1470  {
1471  if( crucnode == SCIPStpunionfindFind(uf, dfstree[k]) )
1472  {
1473  graphmark[dfstree[k]] = FALSE;
1474  solNodes[dfstree[k]] = FALSE;
1475  }
1476  }
1477 
1478  /* update union find */
1479  if( !Is_term(graph->term[node]) && scanned[node] && !pinned[node] && SCIPStpunionfindFind(uf, node) == node )
1480  {
1481  for( int edge = graph->outbeg[node]; edge != EAT_LAST; edge = graph->oeat[edge] )
1482  {
1483  int adjnode = graph->head[edge];
1484 
1485  /* check whether edge 'edge' leads to an ancestor of terminal 'node' */
1486  if( solEdges[edge] == CONNECT && solNodes[adjnode] && graphmark[adjnode] && SCIPStpunionfindFind(uf, adjnode) != node )
1487  {
1488  assert(scanned[adjnode]);
1489  /* meld the heaps */
1490  SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &pheapsize[node], &pheapsize[adjnode]);
1491 
1492  /* update the union-find data structure */
1493  SCIPStpunionfindUnion(uf, node, adjnode, FALSE);
1494 
1495  /* move along the key-path until its end (i.e. until a crucial node is reached) */
1496  while( !nodeIsCrucial(graph, solEdges, adjnode) && !pinned[adjnode] )
1497  {
1498  int e;
1499  for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
1500  {
1501  if( solEdges[e] != -1 )
1502  break;
1503  }
1504 
1505  /* assert that each leaf of the ST is a terminal */
1506  assert( e != EAT_LAST );
1507  adjnode = graph->head[e];
1508 
1509  if( !solNodes[adjnode] )
1510  break;
1511 
1512  assert(scanned[adjnode]);
1513  assert(SCIPStpunionfindFind(uf, adjnode) != node);
1514 
1515  /* update the union-find data structure */
1516  SCIPStpunionfindUnion(uf, node, adjnode, FALSE);
1517 
1518  /* meld the heaps */
1519  SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[adjnode], &pheapsize[node], &pheapsize[adjnode]);
1520  }
1521  }
1522  }
1523 
1524  }
1525  pinned[node] = TRUE;
1526 
1527  return SCIP_OKAY;
1528 }
1529 
1530 
1531 
1532 /** exchanges key-paths star */
1533 static
1535  SCIP* scip, /**< SCIP data structure */
1536  const GRAPH* graph, /**< graph data structure */
1537  const CONN* connectData, /**< data */
1538  const VNOILOC* vnoiData, /**< data */
1539  const KPATHS* keypathsData, /**< key paths */
1540  const SGRAPH* supergraphData, /**< super-graph */
1541  const int* dfstree, /**< DFS tree */
1542  const STP_Bool* scanned, /**< array to mark which nodes have been scanned */
1543  int dfstree_pos, /**< current position in dfs tree */
1544  SOLTREE* soltreeData /**< solution tree data */
1545 )
1546 {
1547  const PATH* mst = supergraphData->mst;
1548  UF* const uf = connectData->uf;
1549  const STP_Bool* const isSupernode = supergraphData->nodeIsLowSupernode;
1550  PHNODE** const boundpaths = connectData->pheap_boundpaths;
1551  int* const pheapsize = connectData->pheap_sizes;
1552  const PATH* const vnoipath = vnoiData->vnoi_path;
1553  const int* const vnoibase = vnoiData->vnoi_base;
1554  const int* const boundedges = connectData->boundaryedges;
1555  const int* const kpnodes = keypathsData->kpnodes;
1556  const int* const kpedges = keypathsData->kpedges;
1557  STP_Bool* const pinned = soltreeData->nodeIsPinned;
1558  const LCNODE* const linkcutNodes = soltreeData->linkcutNodes;
1559  int* const solEdges = soltreeData->solEdges;
1560  int* const graphmark = graph->mark;
1561  STP_Bool* const solNodes = soltreeData->solNodes;
1562  const int nkpnodes = keypathsData->nkpnodes;
1563  const int nkpedges = keypathsData->nkpedges;
1564  const int nsupernodes = supergraphData->nsupernodes;
1565 
1566  /* unmark the original edges spanning the supergraph */
1567  for( int e = 0; e < nkpedges; e++ )
1568  {
1569  assert(solEdges[kpedges[e]] == CONNECT);
1570  solEdges[kpedges[e]] = UNKNOWN;
1571  }
1572 
1573  /* mark all Steiner tree nodes except for those belonging to the root-component as forbidden */
1574  for( int k = keypathsData->rootpathstart; k < nkpnodes; k++ )
1575  {
1576  graphmark[kpnodes[k]] = FALSE;
1577  solNodes[kpnodes[k]] = FALSE;
1578  }
1579 
1580  for( int k = 0; k < dfstree_pos; k++ )
1581  {
1582  const int node = SCIPStpunionfindFind(uf, dfstree[k]);
1583  if( isSupernode[node] || node == dfstree[dfstree_pos] )
1584  {
1585  graphmark[dfstree[k]] = FALSE;
1586  solNodes[dfstree[k]] = FALSE;
1587  }
1588  }
1589 
1590  /* add the new edges reconnecting the (super-) components */
1591  for( int l = 0; l < nsupernodes - 1; l++ )
1592  {
1593  const int edge = (mst[l].edge % 2 == 0)? boundedges[mst[l].edge / 2] : flipedge(boundedges[mst[l].edge / 2]);
1594 
1595  /* change the orientation within the target-component if necessary */
1596  if( !isSupernode[vnoibase[graph->head[edge]]] )
1597  {
1598  int node = vnoibase[graph->head[edge]];
1599  const int nodebase = SCIPStpunionfindFind(uf, node);
1600  assert(isSupernode[nodebase]);
1601 
1602  while( node != nodebase )
1603  {
1604  /* the Steiner tree edge pointing towards the root */
1605  const int e = linkcutNodes[node].edge;
1606 
1607  assert(solEdges[e] == UNKNOWN && solEdges[flipedge(e)] == CONNECT );
1608 
1609  solEdges[e] = CONNECT;
1610  solEdges[flipedge(e)] = UNKNOWN;
1611  node = graph->head[e];
1612  }
1613  }
1614 
1615  /* is the vbase of the current boundary-edge tail in the root-component? */
1616  if( !isSupernode[SCIPStpunionfindFind(uf, vnoibase[graph->tail[edge]])] )
1617  {
1618  int node;
1619  solEdges[edge] = CONNECT;
1620 
1621  for( node = graph->tail[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1622  {
1623  graphmark[node] = FALSE;
1624 
1625  if( solEdges[flipedge(vnoipath[node].edge)] == CONNECT )
1626  solEdges[flipedge(vnoipath[node].edge)] = UNKNOWN;
1627 
1628  solEdges[vnoipath[node].edge] = CONNECT;
1629  }
1630 
1631  assert(!isSupernode[node] && vnoibase[node] == node);
1632  assert(graphmark[node]);
1633 
1634  /* is the pinned node its own component identifier? */
1635  if( !Is_term(graph->term[node]) && scanned[node] && !pinned[node] && SCIPStpunionfindFind(uf, node) == node )
1636  {
1637  graphmark[graph->head[edge]] = FALSE;
1638 
1639  for( int oedge = graph->outbeg[node]; oedge != EAT_LAST; oedge = graph->oeat[oedge] )
1640  {
1641  int head = graph->head[oedge];
1642 
1643  /* check whether edge 'edge' leads to an ancestor of 'node' */
1644  if( solEdges[oedge] == CONNECT && graphmark[head] && solNodes[head] && SCIPStpunionfindFind(uf, head) != node )
1645  {
1646  assert(scanned[head]);
1647 
1648  SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[head], &pheapsize[node], &pheapsize[head]);
1649 
1650  /* update the union-find data structure */
1651  SCIPStpunionfindUnion(uf, node, head, FALSE);
1652 
1653  /* move along the key-path until its end (i.e. until a crucial node is reached) */
1654  while( !nodeIsCrucial(graph, solEdges, head) && !pinned[head] )
1655  {
1656  int e;
1657  for( e = graph->outbeg[head]; e != EAT_LAST; e = graph->oeat[e] )
1658  {
1659  if( solEdges[e] != -1 )
1660  break;
1661  }
1662 
1663  /* assert that each leaf of the ST is a terminal */
1664  assert( e != EAT_LAST );
1665  head = graph->head[e];
1666 
1667  if( !solNodes[head] )
1668  break;
1669 
1670  assert(scanned[head]);
1671  assert(SCIPStpunionfindFind(uf, head) != node);
1672 
1673  /* update the union-find data structure */
1674  SCIPStpunionfindUnion(uf, node, head, FALSE);
1675 
1676  SCIPpairheapMeldheaps(scip, &boundpaths[node], &boundpaths[head], &pheapsize[node], &pheapsize[head]);
1677  }
1678  }
1679  }
1680  }
1681 
1682  /* mark the start node (lying in the root-component of the Steiner tree) of the current boundary-path as pinned,
1683  * so that it may not be removed later on */
1684  pinned[node] = TRUE;
1685 
1686  for( node = graph->head[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1687  {
1688  graphmark[node] = FALSE;
1689  if( solEdges[vnoipath[node].edge] == CONNECT )
1690  solEdges[vnoipath[node].edge] = UNKNOWN;
1691 
1692  solEdges[flipedge(vnoipath[node].edge)] = CONNECT;
1693  }
1694  }
1695  else /* the vbase of the current boundary-edge tail is NOT in the root-component */
1696  {
1697  solEdges[edge] = CONNECT;
1698 
1699  for( int node = graph->tail[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1700  {
1701  graphmark[node] = FALSE;
1702  if( solEdges[vnoipath[node].edge] != CONNECT && solEdges[flipedge(vnoipath[node].edge)] != CONNECT )
1703  {
1704  solEdges[vnoipath[node].edge] = CONNECT;
1705  }
1706  }
1707 
1708  for( int node = graph->head[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1709  {
1710  graphmark[node] = FALSE;
1711 
1712  solEdges[flipedge(vnoipath[node].edge)] = CONNECT;
1713  solEdges[vnoipath[node].edge] = UNKNOWN;
1714  }
1715  }
1716  }
1717 
1718  for( int k = 0; k < nkpnodes; k++ )
1719  {
1720  assert(graphmark[kpnodes[k]] == FALSE);
1721  assert(solNodes[kpnodes[k]] == FALSE);
1722  }
1723 
1724  assert(!graphmark[dfstree[dfstree_pos]]);
1725 
1726  return SCIP_OKAY;
1727 }
1728 
1729 
1730 /** compute cost of alternative key path */
1731 static
1733  SCIP* scip, /**< SCIP data structure */
1734  const GRAPH* graph, /**< graph data structure */
1735  const VNOILOC* vnoiData, /**< data */
1736  const SOLTREE* soltreeData, /**< solution tree data */
1737  SCIP_Real edgecost_old_in, /**< edge cost of old edge */
1738  int boundedge_old, /**< Voronoi boundary edge */
1739  PCMW* pcmwData, /**< data */
1740  int* boundedge_new /**< new Voronoi boundary edge (in/out) */
1741 )
1742 {
1743  SCIP_Real edgecost = FARAWAY;
1744  SCIP_Real edgecost_old = FARAWAY;
1745  SCIP_Real edgecost_new = FARAWAY;
1746  int newedge = *boundedge_new;
1747 
1748  assert(pcmwDataIsClean(graph, pcmwData));
1749 
1750  if( boundedge_old != UNKNOWN )
1751  {
1752  SCIPdebugMessage("get replace path for old edge: \n");
1753 
1754  edgecost_old = getBoundaryPathCostPrized(graph, vnoiData, soltreeData, boundedge_old, pcmwData);
1755 
1756  assert(SCIPisLE(scip, edgecost_old, edgecost_old_in));
1757  assert(SCIPisLT(scip, edgecost_old_in, FARAWAY));
1758 
1759  pcmwDataClean(graph, pcmwData);
1760  }
1761 
1762  if( newedge != UNKNOWN )
1763  {
1764  SCIPdebugMessage("get replace path for new edge: \n");
1765 
1766  edgecost_new = getBoundaryPathCostPrized(graph, vnoiData, soltreeData, newedge, pcmwData);
1767 
1768  pcmwDataClean(graph, pcmwData);
1769  }
1770 
1771  if( SCIPisLT(scip, edgecost_old, edgecost_new) )
1772  {
1773  edgecost = edgecost_old;
1774  *boundedge_new = boundedge_old;
1775  }
1776  else
1777  {
1778  edgecost = edgecost_new;
1779  *boundedge_new = newedge;
1780  }
1781 
1782  assert(SCIPisLT(scip, edgecost, FARAWAY) || (graph_pc_isPcMw(graph) && graph->source == graph->head[*boundedge_new] ) );
1783  assert(*boundedge_new != UNKNOWN);
1784 
1785  return edgecost;
1786 }
1787 
1788 
1789 /** restore supergraph data */
1790 static
1792  const GRAPH* graph, /**< graph data structure */
1793  SGRAPH* supergraphData /**< super-graph*/
1794 )
1795 {
1796  const int* const supernodes = supergraphData->supernodes;
1797  STP_Bool* const nodeIsLowSupernode = supergraphData->nodeIsLowSupernode;
1798 
1799  SCIPfreeMemoryArray(scip, &(supergraphData->mst));
1800 
1801  /* unmark the descendant supervertices */
1802  for( int k = 0; k < supergraphData->nsupernodes - 1; k++ )
1803  nodeIsLowSupernode[supernodes[k]] = FALSE;
1804 
1805 #ifndef NDEBUG
1806  for( int k = 0; k < graph->knots; k++ )
1807  assert(!nodeIsLowSupernode[k]);
1808 #endif
1809 }
1810 
1811 /** compute minimum-spanning tree */
1812 static
1814  SCIP* scip, /**< SCIP data structure */
1815  const GRAPH* graph, /**< graph data structure */
1816  const CONN* connectData, /**< data */
1817  const VNOILOC* vnoiData, /**< data */
1818  const KPATHS* keypathsData, /**< key paths */
1819  int crucnode, /**< node to eliminate */
1820  SOLTREE* soltreeData, /**< solution tree data */
1821  PCMW* pcmwData, /**< data */
1822  SGRAPH* supergraphData /**< super-graph*/
1823 )
1824 {
1825  PATH* mst = NULL;
1826  GRAPH* supergraph = NULL;
1827  UF* const uf = connectData->uf;
1828  const int* const supernodes = supergraphData->supernodes;
1829  int* supernodesid = NULL;
1830  const SCIP_Real* const edgecosts = getEdgeCosts(graph, pcmwData);
1831  const STP_Bool* const isLowerSupernode = supergraphData->nodeIsLowSupernode;
1832  STP_Bool* const solNodes = soltreeData->solNodes;
1833  int* const newedges = soltreeData->newedges;
1834  const PATH* const vnoipath = vnoiData->vnoi_path;
1835  const int* const vnoibase = vnoiData->vnoi_base;
1836  const int* const boundaryedges = connectData->boundaryedges;
1837  SCIP_Real mstcost = 0.0;
1838  const int nboundaryedges = connectData->nboundaryedges;
1839  const int nnodes = graph->knots;
1840  const int nsupernodes = supergraphData->nsupernodes;
1841  /* the (super-) vertex representing the current root-component of the Steiner tree */
1842  const int superroot = supernodes[nsupernodes - 1];
1843 #ifdef SCIP_DEBUG
1844  SCIP_Real mstcost_supergraph = 0.0;
1845  const SCIP_Bool pcmw = graph_pc_isPcMw(graph);
1846 #endif
1847  assert(nboundaryedges > 0);
1848  assert(superroot >= 0);
1849  assert(!supergraphData->mst);
1850 
1851  SCIP_CALL( SCIPallocBufferArray(scip, &supernodesid, nnodes) );
1852 #ifndef NDEBUG
1853  for( int k = 0; k < nnodes; k++ )
1854  supernodesid[k] = -1;
1855 #endif
1856 
1857  /* create a supergraph, having the endpoints of the key-paths incident to the current crucial node as (super-) vertices */
1858  SCIP_CALL( graph_init(scip, &supergraph, nsupernodes, nboundaryedges * 2, 1) );
1859  supergraph->stp_type = STP_SPG;
1860 
1861  for( int k = 0; k < nsupernodes; k++ )
1862  {
1863  supernodesid[supernodes[k]] = k;
1864  graph_knot_add(supergraph, graph->term[supernodes[k]]);
1865  }
1866 
1867  if( pcmwData->isActive )
1868  {
1869  soltreeUnmarkKpNodes(graph, keypathsData, soltreeData);
1870  }
1871 
1872  assert(pcmwDataIsClean(graph, pcmwData));
1873 
1874  SCIPdebugMessage("build super-graph: \n");
1875 
1876  /* add edges to the supergraph */
1877  for( int k = 0; k < nboundaryedges; k++ )
1878  {
1879  SCIP_Real edgecost;
1880  const int edge = boundaryedges[k];
1881  int node = SCIPStpunionfindFind(uf, vnoibase[graph->tail[edge]]);
1882  int adjnode = SCIPStpunionfindFind(uf, vnoibase[graph->head[edge]]);
1883 
1884  /* if node 'node' or 'adjnode' belongs to the root-component, take the (temporary) root-component identifier instead */
1885  node = (isLowerSupernode[node]? node : superroot);
1886  adjnode = (isLowerSupernode[adjnode]? adjnode : superroot);
1887 
1888  /* compute the cost of the boundary-path pertaining to the boundary-edge 'edge' */
1889  edgecost = getBoundaryPathCostPrized(graph, vnoiData, soltreeData, edge, pcmwData);
1890  pcmwDataClean(graph, pcmwData);
1891 
1892  graph_edge_add(scip, supergraph, supernodesid[node], supernodesid[adjnode], edgecost, edgecost);
1893  }
1894 
1895  /* compute a MST on the supergraph */
1896  SCIP_CALL( SCIPallocMemoryArray(scip, &(supergraphData->mst), nsupernodes) );
1897  mst = supergraphData->mst;
1898  SCIP_CALL( graph_path_init(scip, supergraph) );
1899  graph_path_exec(scip, supergraph, MST_MODE, nsupernodes - 1, supergraph->cost, mst);
1900 
1901  assert(pcmwDataIsClean(graph, pcmwData));
1902 
1903  /* compute the cost of the MST */
1904 
1905  mstcost = 0.0;
1906 
1907  for( int l = 0; l < nsupernodes - 1; l++ )
1908  {
1909  /* compute the edge in the original graph corresponding to the current MST edge */
1910  const int edge = (mst[l].edge % 2 == 0)? boundaryedges[mst[l].edge / 2] : flipedge(boundaryedges[mst[l].edge / 2]);
1911 
1912  mstcost += getEdgeCostUnbiased(graph, pcmwData, edge);
1913  mstcost -= pcmwGetNewEdgePrize(graph, solNodes, edge, pcmwData);
1914 
1915 #ifdef SCIP_DEBUG
1916  printf("key vertex mst edge ");
1917  graph_edge_printInfo(graph, edge);
1918  if( pcmw ) printf("nodeweights: %f -> %f \n", graph->prize[graph->tail[edge]], graph->prize[graph->head[edge]]);
1919 
1920  mstcost_supergraph += supergraph->cost[mst[l].edge];
1921 #endif
1922 
1923  assert(newedges[edge] != crucnode && newedges[flipedge(edge)] != crucnode);
1924 
1925  /* mark the edge (in the original graph) as visited */
1926  newedges[edge] = crucnode;
1927 
1928  /* traverse along the boundary-path belonging to the boundary-edge 'edge' */
1929  for( int node = graph->tail[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1930  {
1931  const int e = vnoipath[node].edge;
1932 
1933  /* if edge 'e' and its reversed have not been visited yet */
1934  if( newedges[e] != crucnode && newedges[flipedge(e)] != crucnode )
1935  {
1936  newedges[e] = crucnode;
1937  mstcost += edgecosts[e];
1938  mstcost -= pcmwGetNewEdgePrize(graph, solNodes, e, pcmwData);
1939 
1940 #ifdef SCIP_DEBUG
1941  printf("key vertex mst edge ");
1942  graph_edge_printInfo(graph, e);
1943  if( pcmw ) printf("nodeweights: %f -> %f \n", graph->prize[graph->tail[edge]], graph->prize[graph->head[edge]]);
1944 #endif
1945  }
1946  }
1947 
1948  for( int node = graph->head[edge]; node != vnoibase[node]; node = graph->tail[vnoipath[node].edge] )
1949  {
1950  const int e = vnoipath[node].edge;
1951  const int erev = flipedge(e);
1952 
1953  /* if edge 'e' and its reversed have not been visited yet */
1954  if( newedges[e] != crucnode && newedges[erev] != crucnode )
1955  {
1956  newedges[erev] = crucnode;
1957 
1958  mstcost += edgecosts[e];
1959  mstcost -= pcmwGetNewEdgePrize(graph, solNodes, e, pcmwData);
1960 
1961 #ifdef SCIP_DEBUG
1962  printf("key vertex mst edge ");
1963  graph_edge_printInfo(graph, erev);
1964  if( pcmw ) printf("nodeweights: %f -> %f \n", graph->prize[graph->tail[edge]], graph->prize[graph->head[edge]]);
1965 #endif
1966  }
1967  }
1968  }
1969 #ifdef SCIP_DEBUG
1970  SCIPdebugMessage("mstcost=%f supergraph mstcost=%f \n", mstcost, mstcost_supergraph);
1971 #endif
1972 
1973  pcmwDataClean(graph, pcmwData);
1974 
1975  supergraphData->mstcost = mstcost;
1976 
1977  if( pcmwData->isActive )
1978  {
1979  soltreeMarkKpNodes(graph, keypathsData, soltreeData);
1980  }
1981 
1982  SCIPfreeBufferArray(scip, &supernodesid);
1983  graph_path_exit(scip, supergraph);
1984  graph_free(scip, &supergraph, TRUE);
1985 
1986  return SCIP_OKAY;
1987 }
1988 
1989 
1990 /** get the key-paths star starting from 'keyvertex' */
1991 static
1993  int keyvertex, /**< key vertex to start from */
1994  const GRAPH* graph, /**< graph data structure */
1995  const CONN* connectData, /**< data */
1996  const SOLTREE* soltreeData, /**< solution tree data */
1997  const PCMW* pcmwData, /**< data */
1998  KPATHS* keypathsData, /**< key paths */
1999  SGRAPH* supergraphData, /**< super-graph*/
2000  SCIP_Bool* success /**< success? */
2001 )
2002 {
2003  int* const kpnodes = keypathsData->kpnodes;
2004  int* const kpedges = keypathsData->kpedges;
2005  const int* const solEdges = soltreeData->solEdges;
2006  int* const supernodes = supergraphData->supernodes;
2007  STP_Bool* isLowSupernode = supergraphData->nodeIsLowSupernode;
2008  const STP_Bool* const solNodes = soltreeData->solNodes;
2009  const STP_Bool* const pinned = soltreeData->nodeIsPinned;
2010  int edgeFromRoot = UNKNOWN;
2011  int nkpnodes = 0;
2012  int nkpedges = 0;
2013  int inedgescount = 0;
2014  int nsupernodes = 0;
2015  const SCIP_Real* const edgecosts = getEdgeCosts(graph, pcmwData);
2016 
2017  assert(!pinned[keyvertex] && nodeIsCrucial(graph, solEdges, keyvertex) && !Is_term(graph->term[keyvertex]));
2018 
2019  keypathsData->kpcost = 0.0;
2020  keypathsData->rootpathstart = -1;
2021  keypathsData->nkpedges = -1;
2022  keypathsData->nkpnodes = -1;
2023  supergraphData->nsupernodes = -1;
2024  *success = TRUE;
2025 
2026  /* find all key-paths starting in node 'keyvertex' */
2027  for( int edge = graph->outbeg[keyvertex]; edge != EAT_LAST; edge = graph->oeat[edge] )
2028  {
2029  assert(solNodes[graph->tail[edge]]);
2030 
2031  /* check whether the outgoing edge is in the Steiner tree */
2032  if( (solEdges[edge] == CONNECT && solNodes[graph->head[edge]])
2033  || solEdges[flipedge(edge)] == CONNECT )
2034  {
2035  const int edge_rev = flipedge(edge);
2036 
2037  /* we want to have the edges going into 'keyvertex' to get the prizes right */
2038  keypathsData->kpcost += edgecosts[edge_rev];
2039  inedgescount++;
2040 
2041 #ifdef SCIP_DEBUG
2042  printf("key vertex start edge ");
2043  graph_edge_printInfo(graph, edge_rev);
2044 #endif
2045 
2046  /* does current edge lead to the Steiner tree root? */
2047  if( solEdges[edge_rev] == CONNECT )
2048  {
2049  assert(edgeFromRoot == UNKNOWN);
2050 
2051  edgeFromRoot = edge_rev;
2052  kpedges[nkpedges++] = edgeFromRoot;
2053 
2054  assert(edge == soltreeData->linkcutNodes[keyvertex].edge);
2055  }
2056  else
2057  {
2058  int adjnode = graph->head[edge];
2059  int e = edge;
2060 #ifndef NDEBUG
2061  const int nkpedges_org = nkpedges;
2062  assert(solEdges[flipedge(edge)] == UNKNOWN);
2063 #endif
2064  kpedges[nkpedges++] = e;
2065 
2066  /* move along the key-path until its end (i.e. a crucial or pinned node) is reached */
2067  while( !pinned[adjnode] && !nodeIsCrucial(graph, solEdges, adjnode) && solNodes[adjnode] )
2068  {
2069  /* update the union-find data structure */
2070  SCIPStpunionfindUnion(connectData->uf, keyvertex, adjnode, FALSE);
2071 
2072  kpnodes[nkpnodes++] = adjnode;
2073 
2074  for( e = graph->outbeg[adjnode]; e != EAT_LAST; e = graph->oeat[e] )
2075  {
2076  if( solEdges[e] == CONNECT )
2077  {
2078  /* we sum the cost of the edges pointing toward the key vertex */
2079  keypathsData->kpcost += edgecosts[flipedge(e)];
2080  kpedges[nkpedges++] = e;
2081 #ifdef SCIP_DEBUG
2082  printf("key vertex edge ");
2083  graph_edge_printInfo(graph, flipedge(e));
2084 #endif
2085  break;
2086  }
2087  }
2088 
2089  /* assert that each leaf of the ST is a terminal
2090  * todo check why this fails, the following should never actually happen */
2091  if( e == EAT_LAST )
2092  {
2093  *success = FALSE;
2094  goto TERMINATE;
2095  }
2096 
2097  assert(e != EAT_LAST);
2098  adjnode = graph->head[e];
2099  }
2100 
2101  /* does the last node on the path belong to a removed component? */
2102  if( !solNodes[adjnode] )
2103  {
2104  /* assert that e is not the edge incident to the key vertex */
2105  assert(e != edge && e != edge_rev);
2106  assert(nkpedges - nkpedges_org >= 2);
2107 
2108  keypathsData->kpcost -= edgecosts[flipedge(e)];
2109 #ifdef SCIP_DEBUG
2110  printf("key vertex remove edge ");
2111  graph_edge_printInfo(graph, flipedge(e));
2112 #endif
2113  nkpedges--;
2114  adjnode = graph->tail[e];
2115  if( adjnode != keyvertex )
2116  {
2117  supernodes[nsupernodes++] = adjnode;
2118  isLowSupernode[adjnode] = TRUE;
2119  }
2120  }
2121  else
2122  {
2123  supernodes[nsupernodes++] = adjnode;
2124  isLowSupernode[adjnode] = TRUE;
2125  }
2126 
2127  assert(nkpedges > nkpedges_org);
2128  } /* edge does not go to root */
2129  } /* edge is in solution */
2130  } /* find all (unique) key-paths starting in node 'crucnode' */
2131 
2132  assert(inedgescount >= 0);
2133 
2134  if( inedgescount > 1 && pcmwData->isActive )
2135  {
2136  assert(graph_pc_isPcMw(graph));
2137 
2138  if( !graph_pc_isPc(graph) || graph_pc_knotIsNonLeafTerm(graph, keyvertex) )
2139  {
2140  /* we have implicitly subtracted the prize of the key vertex several times, need to revert */
2141  keypathsData->kpcost += (SCIP_Real)(inedgescount - 1) * graph->prize[keyvertex];
2142  }
2143  else
2144  {
2145  assert(0.0 == graph->prize[keyvertex]);
2146  }
2147  }
2148 
2149  /* traverse the key-path leading to the root-component */
2150  keypathsData->rootpathstart = nkpnodes;
2151  if( edgeFromRoot != UNKNOWN )
2152  {
2153  /* begin with the edge starting in the root-component of key vertex */
2154  int tail = graph->tail[edgeFromRoot];
2155 
2156  while( !pinned[tail] && !nodeIsCrucial(graph, solEdges, tail) && solNodes[tail] )
2157  {
2158  int e;
2159 
2160  kpnodes[nkpnodes++] = tail;
2161 
2162  for( e = graph->inpbeg[tail]; e != EAT_LAST; e = graph->ieat[e] )
2163  {
2164  if( solEdges[e] == CONNECT )
2165  {
2166  assert(solNodes[graph->tail[e]]);
2167  keypathsData->kpcost += edgecosts[e];
2168 #ifdef SCIP_DEBUG
2169  printf("key vertex (root) edge ");
2170  graph_edge_printInfo(graph, e);
2171 #endif
2172  kpedges[nkpedges++] = e;
2173  break;
2174  }
2175  }
2176 
2177  assert(e != EAT_LAST);
2178  tail = graph->tail[e];
2179  }
2180 
2181  supernodes[nsupernodes++] = tail;
2182  }
2183 
2184  /* the last of the key-path nodes to be stored is the current key-node */
2185  kpnodes[nkpnodes++] = keyvertex;
2186 
2187  TERMINATE:
2188 
2189  keypathsData->nkpedges = nkpedges;
2190  keypathsData->nkpnodes = nkpnodes;
2191  supergraphData->nsupernodes = nsupernodes;
2192 }
2193 
2194 
2195 /** preprocessing for Voronoi repair method */
2196 static
2198  SCIP* scip, /**< SCIP data structure */
2199  const GRAPH* graph, /**< graph data structure */
2200  const KPATHS* keypathsData, /**< key paths */
2201  const CONN* connectData, /**< base lists */
2202  const PCMW* pcmwData, /**< data */
2203  VNOILOC* vnoiData, /**< data */
2204  int* nheapelems /**< to store */
2205 )
2206 {
2207  STP_Vectype(int)* blists_start = connectData->blists_start;
2208  PATH* const vnoipath = vnoiData->vnoi_path;
2209  const int* const kpnodes = keypathsData->kpnodes;
2210  int* const vnoibase = vnoiData->vnoi_base;
2211  int* const state = vnoiData->vnoi_nodestate;
2212  const int* const graphmark = graph->mark;
2213  const int nkpnodes = keypathsData->nkpnodes;
2214  const SCIP_Real* const edgecosts = getEdgeCosts(graph, pcmwData);
2215  int count = 0;
2216 
2217  assert(nheapelems);
2218 
2219  for( int k = 0; k < nkpnodes; k++ )
2220  {
2221  STP_Vectype(int) blists_curr = blists_start[kpnodes[k]];
2222  const int size = StpVecGetSize(blists_curr);
2223 
2224  assert(blists_curr != NULL);
2225 
2226  for( int i = 0; i < size; i++ )
2227  {
2228  const int node = blists_curr[i];
2229  assert(graph_knot_isInRange(graph, node));
2230 
2231  for( int edge_in = graph->inpbeg[node]; edge_in != EAT_LAST; edge_in = graph->ieat[edge_in] )
2232  {
2233  const int tail = graph->tail[edge_in];
2234 
2235  /* is tail node not in C and not forbidden? */
2236  if( state[tail] == CONNECT && graphmark[tail] && graphmark[vnoibase[tail]] )
2237  {
2238  const SCIP_Real newdist = vnoipath[tail].dist + edgecosts[edge_in];
2239 
2240  /* is the new path via tail shorter? */
2241  if( SCIPisGT(scip, vnoipath[node].dist, newdist) )
2242  {
2243  vnoipath[node].dist = newdist;
2244  vnoipath[node].edge = edge_in;
2245  vnoibase[node] = vnoibase[tail];
2246  }
2247  }
2248  }
2249 
2250  if( vnoibase[node] != UNKNOWN )
2251  {
2252  graph_pathHeapAdd(vnoipath, node, graph->path_heap, state, &count);
2253  }
2254  }
2255  }
2256 
2257  assert(nkpnodes == 0 || count > 0);
2258 
2259  *nheapelems = count;
2260 }
2261 
2262 /** restore data */
2263 static
2265  const CONN* connectData, /**< base lists */
2266  const KPATHS* keypathsData, /**< key paths */
2267  VNOILOC* vnoiData /**< data */
2268 )
2269 {
2270  STP_Vectype(int)* blists_start = connectData->blists_start;
2271  PATH* vnoipath = vnoiData->vnoi_path;
2272  int* memvbase = vnoiData->memvbase;
2273  int* meminedges = vnoiData->meminedges;
2274  int* vnoibase = vnoiData->vnoi_base;
2275  const int* kpnodes = keypathsData->kpnodes;
2276  SCIP_Real* memvdist = vnoiData->memvdist;
2277  const int nkpnodes = keypathsData->nkpnodes;
2278  int memcount = 0;
2279 
2280  for( int k = 0; k < nkpnodes; k++ )
2281  {
2282  /* restore data of all nodes having the current (internal) key-path node as their voronoi base */
2283  STP_Vectype(int) blists_curr = blists_start[kpnodes[k]];
2284  const int size = StpVecGetSize(blists_curr);
2285 
2286  for( int i = 0; i < size; i++ )
2287  {
2288  const int node = blists_curr[i];
2289  vnoibase[node] = memvbase[memcount];
2290  vnoipath[node].dist = memvdist[memcount];
2291  vnoipath[node].edge = meminedges[memcount];
2292  memcount++;
2293  }
2294  }
2295 
2296  assert(memcount == vnoiData->nmems);
2297  assert(vnoiData->nkpnodes == nkpnodes);
2298 }
2299 
2300 /** reset data */
2301 static
2303  const CONN* connectData, /**< base lists */
2304  const KPATHS* keypathsData, /**< key paths */
2305  const int* graphmark, /**< graph mark */
2306  VNOILOC* vnoiData /**< data */
2307 )
2308 {
2309  STP_Vectype(int)* blists_start = connectData->blists_start;
2310  PATH* vnoipath = vnoiData->vnoi_path;
2311  int* memvbase = vnoiData->memvbase;
2312  int* meminedges = vnoiData->meminedges;
2313  int* state = vnoiData->vnoi_nodestate;
2314  int* vnoibase = vnoiData->vnoi_base;
2315  const int* kpnodes = keypathsData->kpnodes;
2316  SCIP_Real* memvdist = vnoiData->memvdist;
2317  const int nkpnodes = keypathsData->nkpnodes;
2318  int nresnodes = 0;
2319 
2320  /* reset all nodes (referred to as 'C') whose bases are internal nodes of the current key-paths */
2321  for( int k = 0; k < nkpnodes; k++ )
2322  {
2323  /* reset all nodes having the current (internal) key-path node as their Voronoi base */
2324  STP_Vectype(int) blists_curr = blists_start[kpnodes[k]];
2325  const int size = StpVecGetSize(blists_curr);
2326 
2327  for( int i = 0; i < size; i++ )
2328  {
2329  const int node = blists_curr[i];
2330 
2331  assert(graphmark[node]);
2332 
2333  /* store data */
2334  memvbase[nresnodes] = vnoibase[node];
2335  memvdist[nresnodes] = vnoipath[node].dist;
2336  meminedges[nresnodes] = vnoipath[node].edge;
2337  nresnodes++;
2338 
2339  /* reset data */
2340  vnoibase[node] = UNKNOWN;
2341  vnoipath[node].dist = FARAWAY;
2342  vnoipath[node].edge = UNKNOWN;
2343  state[node] = UNKNOWN;
2344  }
2345  }
2346 
2347  vnoiData->nmems = nresnodes;
2348  vnoiData->nkpnodes = nkpnodes;
2349 }
2350 
2351 
2352 #ifndef NDEBUG
2353 static
2355  SCIP* scip, /**< SCIP data structure */
2356  const GRAPH* graph, /**< graph data structure */
2357  const int* solDegree, /**< solution degree */
2358  const LCNODE* linkcutNodes /**< Steiner tree nodes */
2359  )
2360 {
2361  int* solDegreeNew;
2362  SCIP_Bool isValid;
2363  const int nnodes = graph->knots;
2364 
2365  SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &solDegreeNew, nnodes) );
2366 
2367  BMSclearMemoryArray(solDegreeNew, nnodes);
2368 
2369  for( int i = 0; i < nnodes; ++i )
2370  {
2371  const int edge = linkcutNodes[i].edge;
2372 
2373  if( edge == -1 )
2374  continue;
2375 
2376  assert(edge >= 0 && edge < graph->edges);
2377  assert(graph->tail[edge] == i);
2378 
2379  solDegreeNew[graph->tail[edge]]++;
2380  solDegreeNew[graph->head[edge]]++;
2381  }
2382 
2383  isValid = TRUE;
2384 
2385  for( int i = 0; i < nnodes; ++i )
2386  {
2387  if( Is_term(graph->term[i]) || Is_pseudoTerm(graph->term[i]) )
2388  {
2389  assert(UNKNOWN == solDegree[i]);
2390  continue;
2391  }
2392 
2393  if( solDegree[i] != solDegreeNew[i] )
2394  {
2395 #ifdef SCIP_DEBUG
2396  graph_knot_printInfo(graph, i);
2397  SCIPdebugMessage("%d != %d (old/new degree) \n", solDegree[i], solDegreeNew[i]);
2398 #endif
2399  isValid = FALSE;
2400  break;
2401  }
2402  }
2403 
2404  SCIPfreeBufferArray(scip, &solDegreeNew);
2405 
2406  return isValid;
2407 }
2408 
2409 
2410 static
2412  SCIP* scip, /**< SCIP data structure */
2413  const GRAPH* graph, /**< graph data structure */
2414  const STP_Bool* solNodes, /**< Steiner tree nodes */
2415  const LCNODE* linkcutNodes /**< Steiner tree nodes */
2416  )
2417 {
2418  SCIP_Bool isValid;
2419  STP_Bool* solNodes_new;
2420  const int nnodes = graph->knots;
2421 
2422  SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &solNodes_new, nnodes) );
2423 
2424  for( int i = 0; i < nnodes; ++i )
2425  solNodes_new[i] = FALSE;
2426 
2427  for( int i = 0; i < nnodes; ++i )
2428  {
2429  const int edge = linkcutNodes[i].edge;
2430 
2431  if( edge >= 0 )
2432  {
2433  assert(graph->tail[edge] == i);
2434 
2435  solNodes_new[graph->head[edge]] = TRUE;
2436  solNodes_new[graph->tail[edge]] = TRUE;
2437  }
2438  }
2439 
2440  isValid = TRUE;
2441 
2442  for( int i = 0; i < nnodes; ++i )
2443  {
2444  if( solNodes_new[i] != solNodes[i] )
2445  {
2446  isValid = FALSE;
2447  break;
2448  }
2449  }
2450 
2451  SCIPfreeBufferArray(scip, &solNodes_new);
2452 
2453  return isValid;
2454 }
2455 #endif
2456 
2457 
2458 /** reduces solution degree */
2459 static inline
2461  const GRAPH* graph, /**< graph data structure */
2462  int node, /**< replacement edge */
2463  INSERT* insertData /**< insertion data */
2464 )
2465 {
2466  int* const solDegreeNonTerm = insertData->solDegreeNonTerm;
2467 
2468  assert(node >= 0);
2469  assert(insertData->solNodes[node]);
2470 
2471  if( Is_term(graph->term[node]) || Is_pseudoTerm(graph->term[node]) )
2472  {
2473  assert(solDegreeNonTerm[node] == UNKNOWN);
2474  }
2475  else
2476  {
2477  solDegreeNonTerm[node] -= 1;
2478 
2479  assert(solDegreeNonTerm[node] >= 0);
2480  }
2481 }
2482 
2483 
2484 /** increases solution degree */
2485 static inline
2487  const GRAPH* graph, /**< graph data structure */
2488  int node, /**< replacement edge */
2489  INSERT* insertData /**< insertion data */
2490 )
2491 {
2492  int* const solDegreeNonTerm = insertData->solDegreeNonTerm;
2493 
2494  assert(node >= 0);
2495  assert(insertData->solNodes[node]);
2496 
2497  if( Is_term(graph->term[node]) || Is_pseudoTerm(graph->term[node]) )
2498  {
2499  assert(solDegreeNonTerm[node] == UNKNOWN);
2500  }
2501  else
2502  {
2503  solDegreeNonTerm[node] += 1;
2504  }
2505 }
2506 
2507 /** subroutine for insertion heuristic */
2508 static
2510  SCIP_HEURDATA* heurdata, /**< heuristic data */
2511  const GRAPH* graph, /**< graph data structure */
2512  const int* solEdges, /**< Steiner tree edges */
2513  int* solDegreeNonTerm, /**< degree */
2514  SCIP_Bool* nodeIsBlocked, /**< blocked */
2515  int* vertices /**< vertices permuted */
2516  )
2517 {
2518  const int nnodes = graph->knots;
2519  const int nedges = graph->edges;
2520 
2521  assert(solDegreeNonTerm);
2522 
2523  BMSclearMemoryArray(solDegreeNonTerm, nnodes);
2524 
2525  for( int e = 0; e < nedges; e++ )
2526  {
2527  if( solEdges[e] == CONNECT )
2528  {
2529  solDegreeNonTerm[graph->tail[e]]++;
2530  solDegreeNonTerm[graph->head[e]]++;
2531  }
2532  }
2533 
2534  for( int i = 0; i < nnodes; ++i )
2535  {
2536  vertices[i] = i;
2537  nodeIsBlocked[i] = FALSE;
2538 
2539  if( Is_term(graph->term[i]) || Is_pseudoTerm(graph->term[i]) )
2540  solDegreeNonTerm[i] = UNKNOWN;
2541  }
2542 
2543  SCIPrandomPermuteIntArray(heurdata->randnumgen, vertices, 0, nnodes);
2544 }
2545 
2546 /** subroutine for insertion heuristic: prepare insertion of new vertex */
2547 static
2549  SCIP* scip, /**< SCIP data structure */
2550  const GRAPH* graph, /**< graph data structure */
2551  int v_insert, /**< the vertex to add */
2552  int initialEdge, /**< first edge from node to tree */
2553  LCNODE* linkcutNodes, /**< Steiner tree nodes */
2554  INSERT* insertData, /**< insertion data (in/out) */
2555  SCIP_Real* diff /**< the initial diff (out) */
2556  )
2557 {
2558  const int candHeadInitial = graph->head[initialEdge];
2559  const int stp_type = graph->stp_type;
2560  const SCIP_Bool mw = (stp_type == STP_MWCSP || stp_type == STP_RMWCSP);
2561 
2562  assert(graph->tail[initialEdge] == v_insert);
2563  assert(!Is_term(graph->term[v_insert]));
2564 
2565  SCIPlinkcuttreeLink(linkcutNodes, v_insert, candHeadInitial, initialEdge);
2566 
2567  SCIPdebugMessage("try to insert vertex %d \n", v_insert);
2568 
2569  if( mw )
2570  {
2571  assert(!SCIPisPositive(scip, graph->prize[v_insert]));
2572  *diff = -graph->prize[v_insert];
2573  }
2574  else
2575  {
2576  const SCIP_Bool pc = graph_pc_isPc(graph);
2577 
2578  *diff = insertData->edgecosts[initialEdge];
2579 
2580  if( pc )
2581  {
2582  *diff -= graph->prize[v_insert];
2583  }
2584  }
2585 
2586  assert(insertData->solDegreeNonTerm[v_insert] == 0
2587  || (insertData->solDegreeNonTerm[v_insert] == UNKNOWN && graph_pc_isPcMw(graph)));
2588  assert(!insertData->nodeIsBlocked[v_insert]);
2589 
2590  insertData->nInsertions = 0;
2591  insertData->insertionVertex = v_insert;
2592  insertData->solNodes[v_insert] = TRUE;
2593  insertData->nodeIsBlocked[v_insert] = TRUE;
2594  insertionIncrementSolDegree(graph, v_insert, insertData);
2595  insertionIncrementSolDegree(graph, candHeadInitial, insertData);
2596 }
2597 
2598 
2599 /** remove (blocked) chains from tree */
2600 static
2602  const GRAPH* graph, /**< graph data structure */
2603  int v_insert, /**< the vertex to insert */
2604  LCNODE* linkcutNodes, /**< Steiner tree nodes */
2605  INSERT* insertData /**< insertion data */
2606 )
2607 {
2608  SCIP_Bool* const nodeIsBlocked = insertData->nodeIsBlocked;
2609  const int* const blockedList = insertData->blockedList;
2610  STP_Bool* const solNodes = insertData->solNodes;
2611  const int size = insertData->blockedListSize;
2612  int* const solDegreeNonTerm = insertData->solDegreeNonTerm;
2613 
2614  for( int i = 0; i < size; i++ )
2615  {
2616  const int node = blockedList[i];
2617 
2618  assert(node >= 0 && node < graph->knots);
2619  assert(nodeIsBlocked[node]);
2620  assert(solNodes[node]);
2621 
2622  nodeIsBlocked[node] = FALSE;
2623  solNodes[node] = FALSE;
2624 
2625  if( !(Is_term(graph->term[node]) || Is_pseudoTerm(graph->term[node])) )
2626  {
2627  assert(2 == solDegreeNonTerm[node]);
2628 
2629  solDegreeNonTerm[node] = 0;
2630  }
2631  else
2632  {
2633  assert(UNKNOWN == solDegreeNonTerm[node]);
2634  }
2635 
2636  SCIPlinkcuttreeInitNode(&linkcutNodes[node]);
2637  }
2638 
2639 
2640  assert(nodeIsBlocked[v_insert]);
2641  nodeIsBlocked[v_insert] = FALSE;
2642 
2643  for( int k = 0; k < graph->knots; k++ )
2644  assert(!nodeIsBlocked[k]);
2645 
2646  insertData->blockedListSize = 0;
2647 }
2648 
2649 
2650 /** remove (blocked) chains from tree */
2651 static inline
2653  INSERT* insertData /**< insertion data */
2654 )
2655 {
2656  SCIP_Bool* const nodeIsBlocked = insertData->nodeIsBlocked;
2657  const int* const blockedList = insertData->blockedList;
2658  const int size = insertData->blockedListSize;
2659 
2660  for( int i = 0; i < size; i++ )
2661  {
2662  const int node = blockedList[i];
2663 
2664  assert(node >= 0);
2665  assert(nodeIsBlocked[node]);
2666  assert(2 == insertData->solDegreeNonTerm[node]);
2667 
2668  nodeIsBlocked[node] = FALSE;
2669  }
2670 
2671  assert(nodeIsBlocked[insertData->insertionVertex]);
2672  nodeIsBlocked[insertData->insertionVertex] = FALSE;
2673 }
2674 
2675 
2676 /** mark inner chain nodes as blocked */
2677 static inline
2679  const GRAPH* graph, /**< graph data structure */
2680  const LCNODE* lctree, /**< tree */
2681  int chainfirst_index, /**< first chain entry */
2682  int chainlast_index, /**< last chain entry (inside) */
2683  INSERT* insertData /**< insertion data */
2684 )
2685 {
2686  if( chainfirst_index != chainlast_index )
2687  {
2688  SCIP_Bool* const nodeIsBlocked = insertData->nodeIsBlocked;
2689  int* const blockedList = insertData->blockedList;
2690 
2691  for( int node = chainfirst_index; node != chainlast_index; node = lctree[node].parent )
2692  {
2693  int head;
2694 
2695  assert(node != -1);
2696  assert(lctree[node].edge >= 0);
2697 
2698  head = graph->head[lctree[node].edge];
2699 
2700  assert(!nodeIsBlocked[head]);
2701 
2702  nodeIsBlocked[head] = TRUE;
2703  blockedList[insertData->blockedListSize++] = head;
2704  }
2705  }
2706 }
2707 
2708 
2709 /** reinsert (blocked) chains to tree */
2710 static
2712  const GRAPH* graph, /**< graph data structure */
2713  const int* insertCands, /**< insertion candidates */
2714  int lcvertex_insert, /**< the vertex tested for insertion */
2715  LCNODE* linkcutNodes, /**< Steiner tree nodes */
2716  INSERT* insertData /**< insertion data */
2717 )
2718 {
2719  int* chainEnds = insertData->chainEnds;
2720  const int* const addedEdges = insertData->addedEdges;
2721  const int* const cutedgesStart = insertData->cutedgesStart;
2722  const int* const cutedgesEnd = insertData->cutedgesEnd;
2723  const int* const graphTail = graph->tail;
2724  const int* const graphHead = graph->head;
2725  const int v = insertData->insertionVertex;
2726 
2727  SCIPlinkcuttreeEvert(linkcutNodes, lcvertex_insert);
2728 
2729  insertionResetBlockedNodes(insertData);
2730 
2731  /* reinsert each chain and remove the */
2732  for( int k = insertData->nInsertions - 1; k >= 0; k-- )
2733  {
2734  const int cutedgeFirst = cutedgesStart[k];
2735  const int cutedgeLast = cutedgesEnd[k];
2736  const int firstNode = graphTail[cutedgeFirst];
2737  const int secondNode = graphHead[cutedgeFirst];
2738  const int lastNode = chainEnds[k];
2739 
2740  assert(cutedgeLast >= 0);
2741  assert(firstNode == insertData->chainStarts[k]);
2742 
2743  /* remove the newly added edge */
2744  SCIPlinkcuttreeCutNode(&linkcutNodes[graphHead[addedEdges[k]]]);
2745  insertionDecrementSolDegree(graph, graphHead[addedEdges[k]], insertData);
2746 
2747  /* re-link the tail of the chain */
2748  SCIPlinkcuttreeEvert(linkcutNodes, firstNode);
2749  SCIPlinkcuttreeLink(linkcutNodes, firstNode, secondNode, cutedgeFirst);
2750 
2751  /* re-link the head of the chain (if not already done) */
2752  if( lastNode != firstNode )
2753  SCIPlinkcuttreeLink(linkcutNodes, lastNode, graphHead[cutedgeLast], cutedgeLast);
2754 
2755  /* reset solution degrees of chain border vertices */
2756  insertionIncrementSolDegree(graph, graphTail[cutedgeFirst], insertData);
2757  insertionIncrementSolDegree(graph, graphHead[cutedgeLast], insertData);
2758  }
2759 
2760  /* finally, cut the edge added first (if it had been cut during the insertion process, it would have been restored above) */
2761  SCIPlinkcuttreeEvert(linkcutNodes, lcvertex_insert);
2762  SCIPlinkcuttreeCutNode(&linkcutNodes[graphHead[insertCands[0]]]);
2763  insertionDecrementSolDegree(graph, graphHead[insertCands[0]], insertData);
2764 
2765  if( Is_term(graph->term[v]) || Is_pseudoTerm(graph->term[v]) )
2766  {
2767  assert(graph_pc_isPcMw(graph));
2768  insertData->solDegreeNonTerm[v] = UNKNOWN;
2769  }
2770  else
2771  {
2772  insertData->solDegreeNonTerm[v] = 0;
2773  }
2774 
2775  assert(insertData->solNodes[v]);
2776  insertData->solNodes[v] = FALSE;
2777 
2778  for( int k = 0; k < graph->knots; k++ )
2779  assert(!insertData->nodeIsBlocked[k]);
2780 
2781  insertData->blockedListSize = 0;
2782 }
2783 
2784 
2785 /** replaces chain by a single edge */
2786 static
2788  const GRAPH* graph, /**< graph data structure */
2789  int newedge, /**< replacement edge */
2790  LCNODE* lctree, /**< tree */
2791  INSERT* insertData, /**< insertion data */
2792  int v_lc, /**< current vertex */
2793  int headCurr_lc, /**< head of newedge */
2794  int chainfirst_index, /**< first chain entry */
2795  int chainlast_index /**< last chain entry (inside) */
2796 )
2797 {
2798  LCNODE* chainfirst = &lctree[chainfirst_index];
2799  LCNODE* chainlast = &lctree[chainlast_index];
2800  int* const chainStarts = insertData->chainStarts;
2801  int* const chainEnds = insertData->chainEnds;
2802  int* const cutedgesStart = insertData->cutedgesStart;
2803  int* const cutedgesEnd = insertData->cutedgesEnd;
2804  int* const addedEdges = insertData->addedEdges;
2805  const int nInsertions = insertData->nInsertions;
2806  const int newhead = graph->head[newedge];
2807  const int v_insert = insertData->insertionVertex;
2808 
2809  assert(chainlast->edge >= 0);
2810  assert(graph->tail[newedge] == v_insert);
2811 
2812  cutedgesStart[nInsertions] = chainfirst->edge;
2813  cutedgesEnd[nInsertions] = chainlast->edge;
2814  chainStarts[nInsertions] = chainfirst_index;
2815  chainEnds[nInsertions] = chainlast_index;
2816  addedEdges[nInsertions] = newedge;
2817 
2818  insertionIncrementSolDegree(graph, v_insert, insertData);
2819  insertionIncrementSolDegree(graph, newhead, insertData);
2820 
2821  /* decrease the degree of the chain border vertices */
2822  insertionDecrementSolDegree(graph, graph->tail[chainfirst->edge], insertData);
2823  insertionDecrementSolDegree(graph, graph->head[chainlast->edge], insertData);
2824 
2825  SCIPdebugMessage("remove chain %d->%d \n", graph->tail[chainfirst->edge], graph->head[chainlast->edge]);
2826 
2827  insertionBlockChain(graph, lctree, chainfirst_index, chainlast_index, insertData);
2828 
2829  SCIPlinkcuttreeCutNode(chainfirst);
2830  SCIPlinkcuttreeCutNode(chainlast);
2831 
2832  SCIPlinkcuttreeLink(lctree, v_lc, headCurr_lc, newedge);
2833  assert(lctree[v_lc].edge == newedge);
2834 
2835  insertData->nInsertions++;
2836 }
2837 
2838 
2839 /** perform local vertex insertion heuristic on given Steiner tree */
2840 static
2842  SCIP* scip, /**< SCIP data structure */
2843  SCIP_HEURDATA* heurdata, /**< heuristic data */
2844  const GRAPH* graph, /**< graph data structure */
2845  const STP_Bool* solNodes, /**< Steiner tree nodes */
2846  int vertex, /**< vertex to be inserted */
2847  int* insertcands, /**< candidates */
2848  int* ncands /**< number of candidates */
2849  )
2850 {
2851  int insertcount = 0;
2852  const SCIP_Bool mwpc = graph_pc_isPcMw(graph);
2853  const SCIP_Bool rpcmw = graph_pc_isRootedPcMw(graph);
2854 
2855  assert(!Is_term(graph->term[vertex]));
2856  assert(heurdata);
2857 
2858  for( int oedge = graph->outbeg[vertex]; oedge != EAT_LAST; oedge = graph->oeat[oedge] )
2859  {
2860  const int head = graph->head[oedge];
2861 
2862  if( !solNodes[head] )
2863  continue;
2864 
2865  /* skip dummy terminals */
2866  if( mwpc )
2867  {
2868  if( Is_term(graph->term[head]) && !graph_pc_knotIsFixedTerm(graph, head) )
2869  continue;
2870 
2871  /* note: for unrooted problems the (artificial) source is also a fixed terminal */
2872  if( !rpcmw && head == graph->source )
2873  continue;
2874  }
2875 
2876  assert(SCIPisLT(scip, graph->cost[oedge], FARAWAY));
2877 
2878  insertcands[insertcount++] = oedge;
2879  }
2880 
2881  if( insertcount >= 3 )
2882  SCIPrandomPermuteIntArray(heurdata->randnumgen, insertcands, 0, insertcount);
2883 
2884  *ncands = insertcount;
2885 }
2886 
2887 
2888 /** perform local vertex insertion heuristic on given Steiner tree */
2889 static
2891  SCIP* scip, /**< SCIP data structure */
2892  SCIP_HEURDATA* heurdata, /**< heuristic data or NULL */
2893  const GRAPH* graph, /**< graph data structure */
2894  STP_Bool* solNodes, /**< Steiner tree nodes */
2895  LCNODE* linkcutNodes, /**< Steiner tree nodes */
2896  int* solEdges /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
2897  )
2898 {
2899  int* chainStarts;
2900  int* chainEnds;
2901  int* insertCands = NULL;
2902  int* addedEdges = NULL;
2903  int* cutedgesStart = NULL;
2904  int* cutedgesEnd = NULL;
2905  int* solDegree = NULL;
2906  int* vertices = NULL;
2907  int* blockedList = NULL;
2908  SCIP_Bool* nodeIsBlocked = NULL;
2909 
2910  int newnverts = 0;
2911  const int nnodes = graph->knots;
2912  const int nedges = graph->edges;
2913  const int root = graph->source;
2914  const int stp_type = graph->stp_type;
2915  const SCIP_Bool mw = (stp_type == STP_MWCSP || stp_type == STP_RMWCSP);
2916  const SCIP_Bool pc = graph_pc_isPc(graph);
2917  const SCIP_Bool mwpc = graph_pc_isPcMw(graph);
2918  SCIP_Bool solimproved = TRUE;
2919  SCIP_Real* edgecosts;
2920 
2921 
2922 #ifndef NDEBUG
2923  const SCIP_Real initialobj = solstp_getObjBounded(graph, solEdges, 0.0, nedges);
2924  SCIP_Real diff_total = 0.0;
2925 #endif
2926 
2927  if( stp_type != STP_SPG && stp_type != STP_RSMT && stp_type != STP_OARSMT && stp_type != STP_GSTP && !mwpc )
2928  {
2929  SCIPdebugMessage("vertex inclusion does not work for current problem type \n");
2930  return SCIP_OKAY;
2931  }
2932 
2933  if( pc )
2934  {
2935  SCIP_CALL( SCIPallocBufferArray(scip, &edgecosts, nedges) );
2936 
2937  graph_pc_getOrgCosts(scip, graph, edgecosts);
2938  }
2939  else
2940  {
2941  edgecosts = graph->cost;
2942  }
2943 
2944  SCIP_CALL( SCIPallocBufferArray(scip, &chainStarts, nnodes) );
2945  SCIP_CALL( SCIPallocBufferArray(scip, &chainEnds, nnodes) );
2946  SCIP_CALL( SCIPallocBufferArray(scip, &blockedList, nnodes) );
2947  SCIP_CALL( SCIPallocBufferArray(scip, &nodeIsBlocked, nnodes) );
2948  SCIP_CALL( SCIPallocBufferArray(scip, &vertices, nnodes) );
2949  SCIP_CALL( SCIPallocBufferArray(scip, &insertCands, nnodes) );
2950  SCIP_CALL( SCIPallocBufferArray(scip, &addedEdges, nnodes) );
2951  SCIP_CALL( SCIPallocBufferArray(scip, &cutedgesStart, nnodes) );
2952  SCIP_CALL( SCIPallocBufferArray(scip, &cutedgesEnd, nnodes) );
2953  SCIP_CALL( SCIPallocBufferArray(scip, &solDegree, nnodes) );
2954 
2955  insertionInit(heurdata, graph, solEdges, solDegree, nodeIsBlocked, vertices);
2956 
2957  assert(solDegIsValid(scip, graph, solDegree, linkcutNodes));
2958  assert(solNodeIsValid(scip, graph, solNodes, linkcutNodes));
2959  assert(linkcutNodes[graph->source].edge == -1);
2960 
2961  /* main loop */
2962  while( solimproved )
2963  {
2964  INSERT insertData = { .chainStarts = chainStarts, .chainEnds = chainEnds, .edgecosts = edgecosts,
2965  .solDegreeNonTerm = solDegree, .solNodes = solNodes, .nodeIsBlocked = nodeIsBlocked,
2966  .cutedgesStart = cutedgesStart, .cutedgesEnd = cutedgesEnd, .addedEdges = addedEdges,
2967  .blockedList = blockedList, .blockedListSize = 0, .nInsertions = 0, .insertionVertex = -1 };
2968 
2969  solimproved = FALSE;
2970 
2971  for( int i = 0; i < nnodes; i++ )
2972  {
2973  SCIP_Real diff;
2974  int ninsertcands = 0;
2975  const int v = vertices[i];
2976 
2977  assert(v >= 0 && v < graph->knots);
2978 
2979  if( solNodes[v] || graph->grad[v] <= 1 )
2980  continue;
2981 
2982  assert(!Is_term(graph->term[v]));
2983 
2984  insertionGetCandidateEdges(scip, heurdata, graph, solNodes, v, insertCands, &ninsertcands);
2985 
2986  /* if there are less than two edges connecting node i and the current tree, continue */
2987  if( ninsertcands <= 1 )
2988  continue;
2989 
2990  insertionInitInsert(scip, graph, v, insertCands[0], linkcutNodes, &insertData, &diff);
2991 
2992  /* try to add additional edges between new vertex and tree */
2993  for( int k = 1; k < ninsertcands; k++ )
2994  {
2995  int chainfirst = -1;
2996  int chainlast = -1;
2997  const int insertEdgeCurr = insertCands[k];
2998  const int insertHeadCurr = graph->head[insertEdgeCurr];
2999 
3000  if( nodeIsBlocked[insertHeadCurr] )
3001  {
3002  assert(k > 1);
3003 
3004  continue;
3005  }
3006 
3007  SCIPdebugMessage("insertHeadCurr=%d \n", insertHeadCurr);
3008 
3009  SCIPlinkcuttreeEvert(linkcutNodes, v);
3010 
3011  if( mw )
3012  {
3013  const SCIP_Real chainweight = SCIPlinkcuttreeFindMinChainMw(scip, linkcutNodes, graph->prize,
3014  graph->head, solDegree, nodeIsBlocked, insertHeadCurr,
3015  &chainfirst, &chainlast);
3016 
3017  if( SCIPisLT(scip, chainweight, 0.0) )
3018  {
3019  diff += chainweight;
3020  insertionReplaceChain(graph, insertEdgeCurr, linkcutNodes, &insertData, v, insertHeadCurr, chainfirst, chainlast);
3021  }
3022  }
3023  else
3024  {
3025  const SCIP_Real chainweight = SCIPlinkcuttreeFindMaxChain(scip, linkcutNodes, edgecosts, pc ? graph->prize : NULL,
3026  graph->head, solDegree, nodeIsBlocked, insertHeadCurr,
3027  &chainfirst, &chainlast);
3028 
3029  SCIPdebugMessage("chainweight=%f edgecost=%f \n", chainweight, edgecosts[insertEdgeCurr]);
3030 
3031  /* note: comparison needs to be strict to avoid (redundant) removal of current edge */
3032  if( SCIPisGT(scip, chainweight, edgecosts[insertEdgeCurr]) )
3033  {
3034  diff += edgecosts[insertEdgeCurr];
3035  diff -= chainweight;
3036  insertionReplaceChain(graph, insertEdgeCurr, linkcutNodes, &insertData, v, insertHeadCurr, chainfirst, chainlast);
3037  }
3038  }
3039  }
3040 
3041  /* is the new tree better? */
3042  if( SCIPisNegative(scip, diff) )
3043  {
3044  SCIPlinkcuttreeEvert(linkcutNodes, root);
3045  solimproved = TRUE;
3046  newnverts++;
3047  assert(solNodes[v]);
3048 
3049  insertionFinalizeReplacement(graph, v, linkcutNodes, &insertData);
3050 
3051  SCIPdebugMessage("Inclusion: ADDED VERTEX %d \n", v);
3052 #ifndef NDEBUG
3053  diff_total -= diff;
3054 #endif
3055  }
3056  else
3057  {
3058  /* restore the old tree */
3059  insertionRestoreTree(graph, insertCands, v, linkcutNodes, &insertData);
3060  }
3061 
3062  assert(solDegIsValid(scip, graph, solDegree, linkcutNodes));
3063  assert(solNodeIsValid(scip, graph, solNodes, linkcutNodes));
3064 
3065  } /* for i 0:nnodes */
3066  }
3067 
3068  /* free buffer memory */
3069  SCIPfreeBufferArray(scip, &solDegree);
3070  SCIPfreeBufferArray(scip, &cutedgesEnd);
3071  SCIPfreeBufferArray(scip, &cutedgesStart);
3072  SCIPfreeBufferArray(scip, &addedEdges);
3073  SCIPfreeBufferArray(scip, &insertCands);
3074  SCIPfreeBufferArray(scip, &vertices);
3075  SCIPfreeBufferArray(scip, &blockedList);
3076  SCIPfreeBufferArray(scip, &nodeIsBlocked);
3077  SCIPfreeBufferArray(scip, &chainEnds);
3078  SCIPfreeBufferArray(scip, &chainStarts);
3079 
3080  if( pc )
3081  SCIPfreeBufferArray(scip, &edgecosts);
3082 
3083  for( int e = 0; e < nedges; e++ )
3084  solEdges[e] = UNKNOWN;
3085 
3086  if( newnverts > 0 )
3087  {
3088  SCIP_CALL( solstp_prune(scip, graph, solEdges, solNodes) );
3089 
3090  for( int i = 0; i < nnodes; i++ )
3091  SCIPlinkcuttreeInitNode(&linkcutNodes[i]);
3092 
3093  /* create a link-cut tree representing the current Steiner tree */
3094  for( int e = 0; e < nedges; e++ )
3095  {
3096  if( solEdges[e] == CONNECT )
3097  {
3098  assert(solNodes[graph->tail[e]] && solNodes[graph->head[e]]);
3099 
3100  SCIPlinkcuttreeLink(linkcutNodes, graph->head[e], graph->tail[e], flipedge(e));
3101  }
3102  }
3103  SCIPlinkcuttreeEvert(linkcutNodes, root);
3104  }
3105  else
3106  {
3107  SCIPlinkcuttreeEvert(linkcutNodes, root);
3108  for( int i = 0; i < nnodes; i++ )
3109  {
3110  if( solNodes[i] && linkcutNodes[i].edge != -1 )
3111  solEdges[flipedge(linkcutNodes[i].edge)] = CONNECT;
3112  }
3113  }
3114 
3115 #ifndef NDEBUG
3116  {
3117  const SCIP_Real newobj = solstp_getObjBounded(graph, solEdges, 0.0, nedges);
3118  SCIPdebugMessage("vertex inclusion obj before/after: %f/%f \n", initialobj, newobj);
3119  assert(SCIPisLE(scip, newobj + diff_total, initialobj));
3120  assert(SCIPisGE(scip, diff_total, 0.0));
3121  }
3122 #endif
3123 
3124  assert(solstp_isValid(scip, graph, solEdges));
3125 
3126  return SCIP_OKAY;
3127 }
3128 
3129 /** perform local vertex insertion heuristic on given Steiner tree */
3130 static
3132  SCIP* scip, /**< SCIP data structure */
3133  GRAPH* graph, /**< graph data structure */
3134  SCIP_Bool useFast, /**< fast variant? */
3135  STP_Bool* solNodes, /**< Steiner tree nodes */
3136  LCNODE* linkcutNodes, /**< Steiner tree nodes */
3137  int* solEdges, /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
3138  SCIP_Bool* success /**< solution improved? */
3139  )
3140 {
3141  UF uf; /* union-find */
3142  STP_Vectype(int)* blists_start = NULL;
3143  PATH* vnoipath = NULL;
3144  STP_Vectype(int)* lvledges_start = NULL; /* horizontal edges */
3145  PHNODE** boundpaths = NULL;
3146  PHNODE* pheapelements = NULL;
3147  SCIP_Real* memvdist = NULL;
3148  SCIP_Real* edgecostBiased_pc = NULL;
3149  SCIP_Real* edgecostOrg_pc = NULL;
3150  SCIP_Real* prize_pc = NULL;
3151  int* vnoibase = NULL;
3152  int* kpedges = NULL;
3153  int* kpnodes = NULL;
3154  int* dfstree = NULL;
3155  int* newedges = NULL;
3156  int* memvbase = NULL;
3157  int* pheapsize = NULL;
3158  int* boundedges = NULL;
3159  int* meminedges = NULL;
3160  int* supernodes = NULL;
3161  int* prizemarklist = NULL;
3162  STP_Bool* pinned = NULL;
3163  STP_Bool* scanned = NULL;
3164  STP_Bool* supernodesmark = NULL;
3165  STP_Bool* prizemark = NULL;
3166  const int root = graph->source;
3167  const int nnodes = graph->knots;
3168  const int nedges = graph->edges;
3169  const int maxnrestarts = (useFast ? LOCAL_MAXRESTARTS_FAST : LOCAL_MAXRESTARTS);
3170  const int solroot = pcmwGetSolRoot(graph, solEdges);
3171  const STP_Bool mwpc = graph_pc_isPcMw(graph);
3172  SCIP_Bool solimproved = FALSE;
3173 
3174 #ifndef NDEBUG
3175  const SCIP_Real initialobj = solstp_getObjBounded(graph, solEdges, 0.0, graph->edges);
3176  SCIP_Real objimprovement = 0.0;
3177 #endif
3178 
3179  *success = FALSE;
3180 
3181  assert(mwpc || solroot == -1);
3182 
3183  /* memory needed for both Key-Path Elimination and Exchange */
3184  SCIP_CALL( SCIPallocBufferArray(scip, &vnoipath, nnodes) );
3185  SCIP_CALL( SCIPallocBufferArray(scip, &vnoibase, nnodes) );
3186 
3187  /* only needed for Key-Path Elimination */
3188  SCIP_CALL( SCIPallocBufferArray(scip, &newedges, nedges) );
3189  SCIP_CALL( SCIPallocBufferArray(scip, &lvledges_start, nnodes) );
3190  SCIP_CALL( SCIPallocBufferArray(scip, &boundedges, nedges) );
3191 
3192  /* memory needed for both Key-Path Elimination and Exchange */
3193  if( mwpc )
3194  {
3195  if( graph_pc_isPc(graph) )
3196  {
3197  SCIP_CALL(SCIPallocBufferArray(scip, &edgecostOrg_pc, nedges));
3198 
3199  graph_pc_getOrgCosts(scip, graph, edgecostOrg_pc);
3200  }
3201 
3202  SCIP_CALL(SCIPallocBufferArray(scip, &edgecostBiased_pc, nedges));
3203  SCIP_CALL(SCIPallocBufferArray(scip, &prize_pc, nnodes));
3204  SCIP_CALL(SCIPallocBufferArray(scip, &prizemark, nnodes));
3205  SCIP_CALL(SCIPallocBufferArray(scip, &prizemarklist, nnodes));
3206 
3207  for( int k = 0; k < nnodes; k++ )
3208  prizemark[k] = FALSE;
3209 
3210  if( graph_pc_isPc(graph) )
3211  {
3212  graph_pc_getBiased(graph, edgecostBiased_pc, prize_pc);
3213  }
3214  else
3215  {
3216  BMScopyMemoryArray(edgecostBiased_pc, graph->cost, graph->edges);
3217  BMScopyMemoryArray(prize_pc, graph->prize, graph->knots);
3218  }
3219  }
3220 
3221  SCIP_CALL( SCIPallocBufferArray(scip, &scanned, nnodes) );
3222  SCIP_CALL( SCIPallocBufferArray(scip, &pheapsize, nnodes) );
3223  SCIP_CALL( SCIPallocBufferArray(scip, &blists_start, nnodes) );
3224  SCIP_CALL( SCIPallocBufferArray(scip, &memvbase, nnodes) );
3225  SCIP_CALL( SCIPallocBufferArray(scip, &memvdist, nnodes) );
3226  SCIP_CALL( SCIPallocBufferArray(scip, &meminedges, nnodes) );
3227  SCIP_CALL( SCIPallocBufferArray(scip, &boundpaths, nnodes) );
3228  SCIP_CALL( SCIPallocBufferArray(scip, &pinned, nnodes) );
3229  SCIP_CALL( SCIPallocBufferArray(scip, &dfstree, nnodes) );
3230  SCIP_CALL( SCIPallocBufferArray(scip, &supernodesmark, nnodes) );
3231  SCIP_CALL( SCIPallocBufferArray(scip, &supernodes, nnodes) );
3232  SCIP_CALL( SCIPallocBufferArray(scip, &kpnodes, nnodes) );
3233  SCIP_CALL( SCIPallocBufferArray(scip, &kpedges, nnodes) );
3234  SCIP_CALL( SCIPallocMemoryArray(scip, &pheapelements, nedges) );
3235 
3236  for( int k = 0; k < nnodes; k++ )
3237  graph->mark[k] = (graph->grad[k] > 0);
3238 
3239  graph->mark[root] = TRUE;
3240 
3241  SCIP_CALL( SCIPStpunionfindInit(scip, &uf, nnodes) );
3242 
3243  BMSclearMemoryArray(blists_start, nnodes);
3244  BMSclearMemoryArray(lvledges_start, nnodes);
3245 
3246  /* main loop */
3247  for( int nruns = 0, localmoves = 1; nruns < maxnrestarts && localmoves > 0; nruns++ )
3248  {
3249  VNOILOC vnoiData = { .vnoi_path = vnoipath, .vnoi_base = vnoibase, .memvdist = memvdist, .memvbase = memvbase,
3250  .meminedges = meminedges, .vnoi_nodestate = graph->path_state, .nmems = 0, .nkpnodes = -1 };
3251  KPATHS keypathsData = { .kpnodes = kpnodes, .kpedges = kpedges, .kpcost = 0.0, .nkpnodes = 0, .nkpedges = 0,
3252  .kptailnode = -1 };
3253  CONN connectivityData = { .blists_start = blists_start, .pheap_boundpaths = boundpaths, .lvledges_start = lvledges_start,
3254  .pheap_sizes = pheapsize, .uf = &uf, .boundaryedges = boundedges, .nboundaryedges = 0,
3255  .pheap_elements = pheapelements };
3256  SOLTREE soltreeData = { .solNodes = solNodes, .linkcutNodes = linkcutNodes, .solEdges = solEdges, .nodeIsPinned = pinned,
3257  .nodeIsScanned = scanned, .newedges = newedges };
3258  SGRAPH supergraphData = { .supernodes = supernodes, .nodeIsLowSupernode = supernodesmark,
3259  .mst = NULL, .mstcost = 0.0, .nsupernodes = 0 };
3260  PCMW pcmwData = { .prize_biased = prize_pc, .edgecost_biased = edgecostBiased_pc, .edgecost_org = edgecostOrg_pc,
3261  .prizemark = prizemark, .prizemarklist = prizemarklist, .nprizemarks = 0, .isActive = graph_pc_isPcMw(graph),
3262  .solroot = solroot };
3263  int nstnodes = 0;
3264 
3265  localmoves = 0;
3266 
3267  /* find a DFS order of the ST nodes */
3268  dfsorder(graph, solEdges, root, &nstnodes, dfstree);
3269 
3270  /* initialize data structures */
3271  for( int k = 0; k < nnodes; k++ )
3272  {
3273  pinned[k] = FALSE;
3274  scanned[k] = FALSE;
3275  supernodesmark[k] = FALSE;
3276  }
3277 
3278  for( int e = 0; e < nedges; e++ )
3279  pheapelements[e].element = -1;
3280 
3281  for( int e = 0; e < nedges; e++ )
3282  newedges[e] = UNKNOWN;
3283 
3284  if( mwpc )
3285  {
3286  assert(graph->extended);
3287  SCIP_CALL( pcmwUpdate(scip, graph, &soltreeData, &pcmwData) );
3288 
3289  /* compute a Voronoi diagram with the Steiner tree nodes as bases */
3290  graph_voronoi(graph, pcmwData.edgecost_biased, NULL, soltreeData.solNodes, vnoiData.vnoi_base, vnoiData.vnoi_path);
3291  }
3292  else
3293  {
3294  graph_voronoi(graph, graph->cost, NULL, soltreeData.solNodes, vnoiData.vnoi_base, vnoiData.vnoi_path);
3295  }
3296 
3297  SCIP_CALL( connectivityDataInit(scip, graph, &vnoiData, &soltreeData, &pcmwData, &connectivityData) );
3298 
3299  /* henceforth, the union-find structure will be used on the Steiner tree */
3300  assert(uf.nElements == nnodes);
3301  SCIPStpunionfindClear(scip, &uf);
3302 
3303  /* main loop visiting all nodes of the current Steiner tree in post-order */
3304  for( int dfstree_pos = 0; dfstree[dfstree_pos] != root; dfstree_pos++ )
3305  {
3306  /* current crucial node */
3307  const int crucnode = dfstree[dfstree_pos];
3308  int nheapelems = -1;
3309 
3310  assert(!scanned[crucnode]);
3311  scanned[crucnode] = TRUE;
3312 
3313  SCIPdebugMessage("iteration %d (crucial node: %d) \n", dfstree_pos, crucnode);
3314 
3315  if( solroot == crucnode )
3316  {
3317  assert(mwpc && !graph_pc_isRootedPcMw(graph));
3318  continue;
3319  }
3320 
3321  /* has the node been temporarily removed from the Steiner tree? */
3322  if( !graph->mark[crucnode] )
3323  continue;
3324 
3325  /* key vertex elimination: */
3326  /* is node 'crucnode' a removable crucial node? (i.e. not pinned or a terminal) */
3327  if( !pinned[crucnode] && !Is_term(graph->term[crucnode]) && nodeIsCrucial(graph, solEdges, crucnode) )
3328  {
3329  SCIP_Bool allgood;
3330 
3331  getKeyPathsStar(crucnode, graph, &connectivityData, &soltreeData, &pcmwData, &keypathsData, &supergraphData, &allgood);
3332 
3333  /* todo: check this! */
3334  if( !allgood )
3335  {
3336  *success = FALSE;
3337  localmoves = 0;
3338  SCIPdebugMessage("terminate key vertex heuristic \n");
3339  goto TERMINATE;
3340  }
3341 
3342  assert(keypathsData.nkpnodes != 0); /* if there are no key-path nodes, something has gone wrong */
3343 
3344  /* reset all nodes (referred to as 'C' henceforth) whose bases are internal nodes of the current key-paths */
3345  vnoiDataReset(&connectivityData, &keypathsData, graph->mark, &vnoiData);
3346 
3347  SCIP_CALL( connectivityDataKeyElimUpdate(scip, graph, &vnoiData, &supergraphData, crucnode, &connectivityData) );
3348 
3349  /* try to connect the nodes of C (directly) to COMP(C), as a preprocessing for graph_voronoiRepairMult */
3350  vnoiDataRepairPreprocess(scip, graph, &keypathsData, &connectivityData, &pcmwData, &vnoiData, &nheapelems);
3351 
3352  graph_voronoiRepairMult(scip, graph, (pcmwData.isActive? pcmwData.edgecost_biased : graph->cost),
3353  supernodesmark, &nheapelems, vnoibase, connectivityData.boundaryedges, &(connectivityData.nboundaryedges), &uf, vnoipath);
3354 
3355  SCIP_CALL( supergraphComputeMst(scip, graph, &connectivityData, &vnoiData, &keypathsData,
3356  crucnode, &soltreeData, &pcmwData, &supergraphData) );
3357 
3358  assert(crucnode == dfstree[dfstree_pos]);
3359 
3360  SCIPdebugMessage("kpcost=%f mstcost=%f \n", keypathsData.kpcost, supergraphData.mstcost);
3361 
3362  /* improving solution found? */
3363  if( SCIPisLT(scip, supergraphData.mstcost, keypathsData.kpcost) )
3364  {
3365  localmoves++;
3366  solimproved = TRUE;
3367 
3368  SCIPdebugMessage("found improving solution in KEY VERTEX ELIMINATION (round: %d) \n \n ", nruns);
3369 
3370  SCIP_CALL( soltreeElimKeyPathsStar(scip, graph, &connectivityData, &vnoiData, &keypathsData, &supergraphData,
3371  dfstree, scanned, dfstree_pos, &soltreeData) );
3372 
3373 #ifndef NDEBUG
3374  assert((keypathsData.kpcost - supergraphData.mstcost) >= 0.0);
3375  objimprovement += (keypathsData.kpcost - supergraphData.mstcost);
3376 #endif
3377  }
3378  else /* no improving solution has been found during the move */
3379  {
3380  /* meld the heap pertaining to 'crucnode' and all heaps pertaining to descendant key-paths of node 'crucnode' */
3381  for( int k = 0; k < keypathsData.rootpathstart; k++ )
3382  {
3383  SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[kpnodes[k]], &pheapsize[crucnode], &pheapsize[kpnodes[k]]);
3384  }
3385  for( int k = 0; k < supergraphData.nsupernodes - 1; k++ )
3386  {
3387  SCIPpairheapMeldheaps(scip, &boundpaths[crucnode], &boundpaths[supernodes[k]], &pheapsize[crucnode], &pheapsize[supernodes[k]]);
3388  SCIPStpunionfindUnion(&uf, crucnode, supernodes[k], FALSE);
3389  }
3390  }
3391 
3392  supergraphDataRestore(graph, &supergraphData);
3393 
3394  /* restore the original Voronoi diagram */
3395  vnoiDataRestore(&connectivityData, &keypathsData, &vnoiData);
3396  } /* end of key vertex elimination */
3397 
3398  /* Key-Path Exchange:
3399  * If the crucnode has just been eliminated, skip Key-Path Exchange */
3400  if( graph->mark[crucnode] )
3401  {
3402  SCIP_Real edgecost = FARAWAY;
3403  int e = UNKNOWN;
3404  int oldedge = UNKNOWN;
3405  int newedge = UNKNOWN;
3406  SCIP_Bool allgood = TRUE;
3407 
3408  assert(graph->mark[crucnode]);
3409 
3410  /* is crucnode not a crucial node and not a pinned vertex? */
3411  if( (!nodeIsCrucial(graph, solEdges, crucnode) && !pinned[crucnode]) )
3412  continue;
3413 
3414  /* gets key path from crucnode towards tree root */
3415  getKeyPathUpper(scip, crucnode, graph, &soltreeData, &pcmwData, &connectivityData, &keypathsData, &allgood);
3416 
3417  if( !allgood )
3418  {
3419  *success = FALSE;
3420  localmoves = 0;
3421  SCIPdebugMessage("terminate key vertex heuristic \n");
3422  goto TERMINATE;
3423  }
3424 
3425 #ifndef NDEBUG
3426  for( int k = 0; k < nnodes; k++ )
3427  assert(graph->path_state[k] == CONNECT || !graph->mark[k]);
3428 #endif
3429 
3430  /* reset all nodes (henceforth referred to as 'C') whose bases are internal nodes of the current keypath */
3431  vnoiDataReset(&connectivityData, &keypathsData, graph->mark, &vnoiData);
3432 
3433  while( boundpaths[crucnode] != NULL )
3434  {
3435  int base;
3436  int node;
3437 
3438  SCIP_CALL( SCIPpairheapDeletemin(scip, &e, &edgecost, &boundpaths[crucnode], &(pheapsize[crucnode])) );
3439 
3440  assert(e != UNKNOWN);
3441  base = vnoibase[graph->head[e]];
3442 
3443  assert(graph->mark[vnoibase[graph->tail[e]]]);
3444  node = (base == UNKNOWN || !graph->mark[base] )? UNKNOWN : SCIPStpunionfindFind(&uf, base);
3445 
3446  /* does the boundary-path end in the root component? */
3447  if( node != UNKNOWN && node != crucnode && graph->mark[base] )
3448  {
3449  SCIP_CALL( SCIPpairheapInsert(scip, &boundpaths[crucnode],
3450  pheapelements,
3451  e, edgecost, &(pheapsize[crucnode])) );
3452  break;
3453  }
3454  }
3455 
3456  if( boundpaths[crucnode] == NULL )
3457  oldedge = UNKNOWN;
3458  else
3459  oldedge = e;
3460 
3461  /* try to connect the nodes of C (directly) to COMP(C), as a preprocessing for Voronoi-repair */
3462  vnoiDataRepairPreprocess(scip, graph, &keypathsData, &connectivityData, &pcmwData, &vnoiData, &nheapelems);
3463 
3464  newedge = UNKNOWN;
3465 
3466  /* if there is no key path, nothing has to be repaired */
3467  if( keypathsData.nkpnodes > 0 )
3468  {
3469  graph_voronoiRepair(scip, graph, pcmwData.isActive? pcmwData.edgecost_biased : graph->cost,
3470  pcmwData.edgecost_org, &nheapelems, vnoibase, vnoipath, &newedge, crucnode, &uf);
3471  }
3472  else
3473  {
3474  newedge = linkcutNodes[crucnode].edge;
3475  }
3476 
3477  edgecost = getKeyPathReplaceCost(scip, graph, &vnoiData, &soltreeData, edgecost, oldedge, &pcmwData, &newedge);
3478 
3479  if( SCIPisLT(scip, edgecost, keypathsData.kpcost) )
3480  {
3481  localmoves++;
3482  solimproved = TRUE;
3483 
3484  SCIPdebugMessage( "ADDING NEW KEY PATH (%f )\n\n", edgecost - keypathsData.kpcost );
3485 #ifndef NDEBUG
3486  assert((keypathsData.kpcost - edgecost) >= 0.0);
3487  objimprovement += (keypathsData.kpcost - edgecost);
3488  assert(crucnode == dfstree[dfstree_pos]);
3489 #endif
3490 
3491  SCIP_CALL( soltreeExchangeKeyPath(scip, graph, &connectivityData, &vnoiData, &keypathsData,
3492  dfstree, scanned, dfstree_pos, newedge, &soltreeData) );
3493  }
3494 
3495  /* restore the original Voronoi diagram */
3496  vnoiDataRestore(&connectivityData, &keypathsData, &vnoiData);
3497  }
3498  }
3499 
3500 
3501  /**********************************************************/
3502 
3503  TERMINATE:
3504 
3505  assert(uf.nElements == nnodes);
3506  SCIPStpunionfindClear(scip, &uf);
3507 
3508  /* free data structures */
3509 
3510  for( int k = nnodes - 1; k >= 0; k-- )
3511  {
3512  if( boundpaths[k] )
3513  SCIPpairheapFree(scip, &boundpaths[k]);
3514  }
3515 
3516  /* has there been a move during this run? */
3517  if( localmoves > 0 )
3518  {
3519  for( int i = 0; i < nnodes; i++ )
3520  {
3521  solNodes[i] = FALSE;
3522  graph->mark[i] = (graph->grad[i] > 0);
3523  SCIPlinkcuttreeInitNode(&linkcutNodes[i]);
3524  }
3525 
3526  graph->mark[root] = TRUE;
3527 
3528  /* create a link-cut tree representing the current Steiner tree */
3529  for( int e = 0; e < nedges; e++ )
3530  {
3531  assert(graph->head[e] == graph->tail[flipedge(e)]);
3532 
3533  /* if edge e is in the tree, so are its incident vertices */
3534  if( solEdges[e] != -1 )
3535  {
3536  assert(CONNECT == solEdges[e]);
3537 
3538  solNodes[graph->tail[e]] = TRUE;
3539  solNodes[graph->head[e]] = TRUE;
3540  SCIPlinkcuttreeLink(linkcutNodes, graph->head[e], graph->tail[e], flipedge(e));
3541  }
3542  }
3543  assert( linkcutNodes[root].edge == -1 );
3544  linkcutNodes[root].edge = -1;
3545  }
3546  } /* main loop */
3547 
3548  /* free data structures */
3549  SCIPStpunionfindFreeMembers(scip, &uf);
3550  SCIPfreeMemoryArray(scip, &pheapelements);
3551 
3552  for( int k = nnodes - 1; k >= 0; k-- )
3553  {
3554  StpVecFree(scip, lvledges_start[k]);
3555  StpVecFree(scip, blists_start[k]);
3556  }
3557 
3558  SCIPfreeBufferArray(scip, &kpedges);
3559  SCIPfreeBufferArray(scip, &kpnodes);
3560  SCIPfreeBufferArray(scip, &supernodes);
3561  SCIPfreeBufferArray(scip, &supernodesmark);
3562  SCIPfreeBufferArray(scip, &dfstree);
3563  SCIPfreeBufferArray(scip, &pinned);
3564  SCIPfreeBufferArray(scip, &boundpaths);
3565  SCIPfreeBufferArray(scip, &meminedges);
3566  SCIPfreeBufferArray(scip, &memvdist);
3567  SCIPfreeBufferArray(scip, &memvbase);
3568  SCIPfreeBufferArray(scip, &blists_start);
3569  SCIPfreeBufferArray(scip, &pheapsize);
3570  SCIPfreeBufferArray(scip, &scanned);
3571  SCIPfreeBufferArrayNull(scip, &prizemarklist);
3572  SCIPfreeBufferArrayNull(scip, &prizemark);
3573  SCIPfreeBufferArrayNull(scip, &prize_pc);
3574  SCIPfreeBufferArrayNull(scip, &edgecostBiased_pc);
3575  SCIPfreeBufferArrayNull(scip, &edgecostOrg_pc);
3576  SCIPfreeBufferArray(scip, &boundedges);
3577  SCIPfreeBufferArray(scip, &lvledges_start);
3578  SCIPfreeBufferArray(scip, &newedges);
3579  SCIPfreeBufferArray(scip, &vnoibase);
3580  SCIPfreeBufferArray(scip, &vnoipath);
3581  /******/
3582 
3583  if( solimproved )
3584  {
3585  SCIP_CALL( solstp_pruneFromEdges(scip, graph, solEdges) );
3586  *success = TRUE;
3587  }
3588 
3589 #ifndef NDEBUG
3590  {
3591  const SCIP_Real newobj = solstp_getObjBounded(graph, solEdges, 0.0, nedges);
3592  SCIPdebugMessage("key vertex heuristic obj before/after: %f/%f (improvement=%f)\n", initialobj, newobj, objimprovement);
3593  assert(SCIPisLE(scip, newobj + objimprovement, initialobj));
3594  }
3595 #endif
3596 
3597  return SCIP_OKAY;
3598 }
3599 
3600 
3601 /*
3602  * Callback methods of primal heuristic
3603  */
3604 
3605 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
3606 static
3607 SCIP_DECL_HEURCOPY(heurCopyLocal)
3608 { /*lint --e{715}*/
3609  assert(scip != NULL);
3610  assert(heur != NULL);
3611  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
3612 
3613  /* call inclusion method of primal heuristic */
3615 
3616  return SCIP_OKAY;
3617 }
3618 
3619 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
3620 static
3621 SCIP_DECL_HEURFREE(heurFreeLocal)
3622 { /*lint --e{715}*/
3623  SCIP_HEURDATA* heurdata;
3624 
3625  assert(heur != NULL);
3626  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
3627  assert(scip != NULL);
3628 
3629  /* free heuristic data */
3630  heurdata = SCIPheurGetData(heur);
3631 
3632  assert(heurdata != NULL);
3633 
3634  SCIPfreeRandom(scip, &heurdata->randnumgen);
3635  SCIPfreeMemory(scip, &heurdata);
3636  SCIPheurSetData(heur, NULL);
3637 
3638  return SCIP_OKAY;
3639 }
3640 
3641 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
3642 static
3643 SCIP_DECL_HEURINITSOL(heurInitsolLocal)
3644 { /*lint --e{715}*/
3645  SCIP_HEURDATA* heurdata;
3646 
3647  assert(heur != NULL);
3648  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
3649  assert(scip != NULL);
3650 
3651  /* free heuristic data */
3652  heurdata = SCIPheurGetData(heur);
3653 
3654  heurdata->nfails = 1;
3655  heurdata->nbestsols = -1;
3656  heurdata->nbestsols_min = -1;
3657 
3658  SCIP_CALL( SCIPallocMemoryArray(scip, &(heurdata->lastsolindices), heurdata->maxnsols) );
3659 
3660  for( int i = 0; i < heurdata->maxnsols; i++ )
3661  heurdata->lastsolindices[i] = -1;
3662 
3663  return SCIP_OKAY;
3664 }
3665 
3666 
3667 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */
3668 static
3669 SCIP_DECL_HEUREXITSOL(heurExitsolLocal)
3670 { /*lint --e{715}*/
3671  SCIP_HEURDATA* heurdata;
3672 
3673  assert(heur != NULL);
3674  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
3675  assert(scip != NULL);
3676 
3677  /* free heuristic data */
3678  heurdata = SCIPheurGetData(heur);
3679  assert(heurdata != NULL);
3680  assert(heurdata->lastsolindices != NULL);
3681  SCIPfreeMemoryArray(scip, &(heurdata->lastsolindices));
3682 
3683  return SCIP_OKAY;
3684 }
3685 
3686 
3687 /** execution method of primal heuristic */
3688 static
3689 SCIP_DECL_HEUREXEC(heurExecLocal)
3690 { /*lint --e{715}*/
3691  SCIP_HEURDATA* heurdata;
3692  SCIP_PROBDATA* probdata;
3693  GRAPH* graph; /* graph structure */
3694  SCIP_SOL* initialsol; /* initial solution */
3695  SCIP_SOL** sols; /* solutions */
3696  int v;
3697  int nActiveSols;
3698  int nsols; /* number of all solutions found so far */
3699  int nedges;
3700  int* results;
3701  int* lastsolindices;
3702  SCIP_Real initialsol_obj;
3703 
3704  assert(heur != NULL);
3705  assert(scip != NULL);
3706  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
3707  assert(result != NULL);
3708 
3709  /* get heuristic data */
3710  heurdata = SCIPheurGetData(heur);
3711  assert(heurdata != NULL);
3712  lastsolindices = heurdata->lastsolindices;
3713  assert(lastsolindices != NULL);
3714 
3715  probdata = SCIPgetProbData(scip);
3716  assert(probdata != NULL);
3717 
3718  graph = SCIPprobdataGetGraph(probdata);
3719  assert(graph != NULL);
3720 
3721  *result = SCIP_DIDNOTRUN;
3722 
3723  /* the local heuristics may not work correctly for several problem variants */
3724  if( !probtypeIsValidForLocal(graph) )
3725  return SCIP_OKAY;
3726 
3727  /* no solution available? */
3728  if( SCIPgetBestSol(scip) == NULL )
3729  return SCIP_OKAY;
3730 
3731  sols = SCIPgetSols(scip);
3732  nsols = SCIPgetNSols(scip);
3733  nedges = graph->edges;
3734 
3735  assert(heurdata->maxnsols >= 0);
3736 
3737  nActiveSols = MIN(heurdata->maxnsols, nsols);
3738 
3739  /* only process each solution once */
3740  for( v = 0; v < nActiveSols; v++ )
3741  {
3742  if( SCIPsolGetIndex(sols[v]) != lastsolindices[v] )
3743  {
3744  /* shift all solution indices right of the new solution index */
3745  for( int i = nActiveSols - 1; i >= v + 1; i-- )
3746  lastsolindices[i] = lastsolindices[i - 1];
3747  break;
3748  }
3749  }
3750 
3751  /* no new solution available? */
3752  if( v == nActiveSols )
3753  return SCIP_OKAY;
3754 
3755  if( heurdata->nbestsols == -1 )
3756  initSolNumberBounds(scip, heurdata);
3757 
3758  initialsol = sols[v];
3759  lastsolindices[v] = SCIPsolGetIndex(initialsol);
3760 
3761  /* solution not good enough? */
3762  if( (v > heurdata->nbestsols && !(heurdata->maxfreq)) )
3763  return SCIP_OKAY;
3764 
3765  /* has the new solution been found by this very heuristic
3766  * and not among the elite solutions? (note that there can still be an improvement because heuristic is aborted early) */
3767  if( SCIPsolGetHeur(initialsol) == heur && v > DEFAULT_NELITESOLS )
3768  return SCIP_OKAY;
3769 
3770  *result = SCIP_DIDNOTFIND;
3771 
3772  assert(SCIPprobdataGetVars(scip) != NULL);
3773  assert(SCIPprobdataGetNVars(scip) == nedges);
3774 
3775  /* allocate memory */
3776  SCIP_CALL( SCIPallocBufferArray(scip, &results, nedges) );
3777 
3778  /* write solution into 'results' array */
3779  solGetStpSol(scip, graph, initialsol, results);
3780 
3781  if( !solstp_isValid(scip, graph, results) )
3782  {
3783  SCIPfreeBufferArray(scip, &results);
3784 
3785  return SCIP_OKAY;
3786  }
3787 
3788  /* initial pruning necessary? */
3789  if( solNeedsPruning(initialsol) )
3790  {
3791  SCIP_CALL( solPrune(scip, graph, results) );
3792  }
3793 
3794  initialsol_obj = solstp_getObjBounded(graph, results, 0.0, nedges);
3795 
3796  /* execute local heuristics */
3797  SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, results) );
3798 
3799 #if 0
3800  if( graph_pc_isPcMw(graph) )
3801  SCIP_CALL( SCIPStpHeurLocalExtendPcMwImp(scip, graph, results) );
3802 #endif
3803 
3804  /* finally, try to add the solution to the SCIP pool */
3805  SCIP_CALL( solAddTry(scip, sols, heur, graph, initialsol_obj, initialsol, results, result) );
3806 
3807  SCIPfreeBufferArray(scip, &results);
3808 
3809  return SCIP_OKAY;
3810 }
3811 
3812 
3813 /** perform local heuristics on a given Steiner tree */
3814 static
3816  SCIP* scip, /**< SCIP data structure */
3817  GRAPH* graph, /**< graph data structure */
3818  SCIP_Bool useFast, /**< fast variant? */
3819  int* solEdges /**< array indicating whether an arc is part of the solution: CONNECTED/UNKNOWN (in/out) */
3820  )
3821 {
3822  SCIP_HEUR* heur;
3823  SCIP_HEURDATA* heurdata;
3824  LCNODE* linkcutNodes;
3825  const int root = graph->source;
3826  const int nnodes = graph->knots;
3827  STP_Bool* solNodes;
3828  const STP_Bool mwpc = graph_pc_isPcMw(graph);
3829  SCIP_Bool success = FALSE;
3830 #ifndef NDEBUG
3831  const SCIP_Real initialobj = solstp_getObjBounded(graph, solEdges, 0.0, graph->edges);
3832 #endif
3833 
3834  assert(graph && solEdges);
3835  assert(graph_valid(scip, graph));
3836 
3837  if( graph->grad[root] == 0 || graph->terms == 1 )
3838  return SCIP_OKAY;
3839 
3840  if( SCIPisStopped(scip) )
3841  return SCIP_OKAY;
3842 
3843  if( mwpc )
3844  {
3845  assert(graph->extended);
3846 
3847  if( solIsTrivialPcMw(graph, solEdges) )
3848  return SCIP_OKAY;
3849  }
3850 
3851  heur = SCIPfindHeur(scip, "local");
3852  assert(heur != NULL);
3853  heurdata = SCIPheurGetData(heur);
3854  assert(heurdata != NULL);
3855 
3856  SCIP_CALL( SCIPallocBufferArray(scip, &linkcutNodes, nnodes) );
3857  SCIP_CALL( SCIPallocBufferArray(scip, &solNodes, nnodes) );
3858 
3859  /* now call the actual heuristics */
3860 
3861  if( mwpc )
3862  {
3863  SCIP_CALL( SCIPStpHeurLocalExtendPcMw(scip, graph, graph->cost, solEdges, solNodes) );
3864  }
3865 
3866  markSolTreeNodes(scip, graph, solEdges, linkcutNodes, solNodes);
3867 
3868  SCIP_CALL( localVertexInsertion(scip, heurdata, graph, solNodes, linkcutNodes, solEdges) );
3869 
3870  SCIP_CALL( localKeyVertexHeuristics(scip, graph, useFast, solNodes, linkcutNodes, solEdges, &success) );
3871 
3872  if( success && !useFast )
3873  {
3874  markSolTreeNodes(scip, graph, solEdges, linkcutNodes, solNodes);
3875  SCIP_CALL( localVertexInsertion(scip, heurdata, graph, solNodes, linkcutNodes, solEdges) );
3876  }
3877 
3878  if( success && mwpc && !useFast )
3879  {
3880  SCIP_CALL( SCIPStpHeurLocalExtendPcMw(scip, graph, graph->cost, solEdges, solNodes) );
3881  }
3882 
3883 #ifndef NDEBUG
3884  {
3885  const SCIP_Real newobj = solstp_getObjBounded(graph, solEdges, 0.0, graph->edges);
3886  assert(SCIPisLE(scip, newobj, initialobj));
3887  assert(solstp_isValid(scip, graph, solEdges));
3888  }
3889 #endif
3890 
3891  SCIPfreeBufferArray(scip, &solNodes);
3892  SCIPfreeBufferArray(scip, &linkcutNodes);
3893 
3894  return SCIP_OKAY;
3895 }
3896 
3897 
3898 /** perform local heuristics on a given Steiner tree */
3900  SCIP* scip, /**< SCIP data structure */
3901  GRAPH* graph, /**< graph data structure */
3902  int* solEdges /**< array indicating whether an arc is part of the solution: CONNECTED/UNKNOWN (in/out) */
3903  )
3904 {
3905  SCIP_CALL( localRun(scip, graph, FALSE, solEdges) );
3906  return SCIP_OKAY;
3907 }
3908 
3909 
3910 /** perform local heuristics on a given Steiner tree */
3912  SCIP* scip, /**< SCIP data structure */
3913  GRAPH* graph, /**< graph data structure */
3914  int* solEdges /**< array indicating whether an arc is part of the solution: CONNECTED/UNKNOWN (in/out) */
3915  )
3916 {
3917  SCIP_CALL( localRun(scip, graph, TRUE, solEdges) );
3918  return SCIP_OKAY;
3919 }
3920 
3921 
3922 /** Implication based local heuristic for (R)PC and MW */
3924  SCIP* scip, /**< SCIP data structure */
3925  const GRAPH* graph, /**< graph data structure */
3926  int* result /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
3927 )
3928 {
3929  const int* starts = SCIPStpGetPcImplStarts(scip);
3930  const int* verts = SCIPStpGetPcImplVerts(scip);
3931 
3932  assert(graph_pc_isPcMw(graph));
3933  assert(0 && "currently not supported");
3934 
3935  if( starts != NULL )
3936  {
3937  const int nnodes = graph->knots;
3938  STP_Bool* stvertex;
3939  int nfound = 0;
3940  int ptermcount = 0;
3941 
3942  assert(graph->extended);
3943  assert(verts != NULL);
3944 
3945  SCIPallocBufferArray(scip, &stvertex, nnodes);
3946 
3947  solstp_setVertexFromEdge(graph, result, stvertex);
3948 
3949  for( int i = 0; i < nnodes; i++ )
3950  {
3951  if( !(Is_term(graph->term[i]) || Is_pseudoTerm(graph->term[i])) )
3952  continue;
3953 
3954  assert(!graph_pc_knotIsFixedTerm(graph, i));
3955 
3956  ptermcount++;
3957 
3958  if( stvertex[i] )
3959  continue;
3960 
3961  for( int j = starts[ptermcount - 1]; j < starts[ptermcount]; j++ )
3962  {
3963  const int vert = verts[j];
3964  if( stvertex[vert] )
3965  {
3966  /* now connect the vertex */
3967 
3968  graph_knot_printInfo(graph, i);
3969  nfound++;
3970  break;
3971  }
3972  }
3973  }
3974 
3975  assert(ptermcount == graph_pc_nNonFixedTerms(graph));
3976 
3977  if( nfound > 0 )
3978  {
3979  printf("nfound: %d \n\n\n", nfound);
3980  }
3981  else
3982  printf("none %d \n", 0);
3983 
3984  SCIPfreeBufferArray(scip, &stvertex);
3985  }
3986 
3987  return SCIP_OKAY;
3988 }
3989 
3990 /** Greedy Extension local heuristic for (R)PC and MW */
3992  SCIP* scip, /**< SCIP data structure */
3993  GRAPH* graph, /**< graph data structure */
3994  const SCIP_Real* cost, /**< edge cost array */
3995  int* stedge, /**< initialized array to indicate whether an edge is part of the Steiner tree */
3996  STP_Bool* stvertex /**< uninitialized array to indicate whether a vertex is part of the Steiner tree */
3997  )
3998 {
3999  GNODE candidates[GREEDY_EXTENSIONS_MW];
4000  int candidatesup[GREEDY_EXTENSIONS_MW];
4001 
4002  PATH* path;
4003  PATH* orgpath;
4004  SCIP_PQUEUE* pqueue;
4005  SCIP_Real bestsolval;
4006 
4007  int nextensions;
4008  const int greedyextensions = (graph->stp_type == STP_MWCSP) ? GREEDY_EXTENSIONS_MW : GREEDY_EXTENSIONS;
4009  const int nedges = graph->edges;
4010  const int nnodes = graph->knots;
4011  const int root = graph->source;
4012  STP_Bool* stvertextmp;
4013  SCIP_Bool extensions = FALSE;
4014 
4015 #ifndef NDEBUG
4016  SCIP_Real objdiff = 0.0;
4017  const SCIP_Real initialobj = solstp_getObjBounded(graph, stedge, 0.0, nedges);
4018 #endif
4019 
4020  assert(GREEDY_EXTENSIONS_MW >= greedyextensions);
4021  assert(scip != NULL);
4022  assert(graph != NULL);
4023  assert(stedge != NULL);
4024  assert(cost != NULL);
4025  assert(stvertex != NULL);
4026  assert(graph->extended);
4027 
4028  graph_pc_2transcheck(scip, graph);
4029  SCIP_CALL( SCIPallocBufferArray(scip, &stvertextmp, nnodes) );
4030  SCIP_CALL( SCIPallocBufferArray(scip, &orgpath, nnodes) );
4031  SCIP_CALL( SCIPallocBufferArray(scip, &path, nnodes) );
4032 
4033  /* initialize solution vertex array with FALSE */
4034  BMSclearMemoryArray(stvertex, nnodes);
4035 
4036  stvertex[root] = TRUE;
4037 
4038  for( int j = 0; j < nnodes; j++ )
4039  path[j].edge = UNKNOWN;
4040 
4041  for( int e = 0; e < nedges; e++ )
4042  {
4043  if( stedge[e] == CONNECT )
4044  {
4045  path[graph->head[e]].edge = e;
4046  stvertex[graph->head[e]] = TRUE;
4047  }
4048  }
4049 
4050 #ifndef NDEBUG
4051  for( int e = 0; e < nedges; e++ )
4052  if( stedge[e] == CONNECT )
4053  assert(stvertex[graph->tail[e]]);
4054 #endif
4055 
4056  graph_path_st_pcmw_extend(scip, graph, cost, FALSE, path, stvertex, &extensions);
4057 
4058  BMScopyMemoryArray(orgpath, path, nnodes);
4059 
4060  /*** compute solution value and save greedyextensions many best unconnected nodes ***/
4061 
4062  SCIP_CALL( SCIPpqueueCreate(&pqueue, greedyextensions, 2.0, GNODECmpByDist, NULL) );
4063 
4064  assert(orgpath[root].edge == UNKNOWN);
4065 
4066  bestsolval = 0.0;
4067  nextensions = 0;
4068 
4069  for( int i = 0; i < nnodes; i++ )
4070  {
4071  if( graph->grad[i] == 0 || root == i )
4072  continue;
4073 
4074  if( Is_term(graph->term[i]) && !graph_pc_knotIsFixedTerm(graph, i) )
4075  continue;
4076 
4077  if( stvertex[i] )
4078  {
4079  assert(orgpath[i].edge >= 0);
4080 
4081  bestsolval += graph->cost[orgpath[i].edge];
4082 
4083  if( Is_pseudoTerm(graph->term[i]) )
4084  {
4085  bestsolval -= graph->prize[i];
4086  }
4087  }
4088  else if( orgpath[i].edge != UNKNOWN && Is_pseudoTerm(graph->term[i]) )
4089  {
4090  SCIP_CALL( addToCandidates(scip, graph, path, i, greedyextensions, &nextensions, candidates, pqueue) );
4091  }
4092  }
4093 
4094  for( int restartcount = 0; restartcount < GREEDY_MAXRESTARTS && !graph_pc_isRootedPcMw(graph); restartcount++ )
4095  {
4096  int l = 0;
4097  SCIP_Bool extensionstmp = FALSE;
4098  int extcount = nextensions;
4099 
4100  /* write extension candidates into array, from max to min */
4101  while( SCIPpqueueNElems(pqueue) > 0 )
4102  {
4103  GNODE* min = (GNODE*) SCIPpqueueRemove(pqueue);
4104  assert(extcount > 0);
4105  candidatesup[--extcount] = min->number;
4106  }
4107  assert(extcount == 0);
4108 
4109  /* iteratively insert new subpaths and try to improve solution */
4110  for( ; l < nextensions; l++ )
4111  {
4112  const int extensioncand = candidatesup[l];
4113  if( !stvertex[extensioncand] )
4114  {
4115  SCIP_Real newsolval = 0.0;
4116  int k = extensioncand;
4117 
4118  BMScopyMemoryArray(stvertextmp, stvertex, nnodes);
4119  BMScopyMemoryArray(path, orgpath, nnodes);
4120 
4121  /* add new extension */
4122  while( !stvertextmp[k] )
4123  {
4124  stvertextmp[k] = TRUE;
4125  assert(orgpath[k].edge != UNKNOWN);
4126  k = graph->tail[orgpath[k].edge];
4127  assert(k != extensioncand);
4128  }
4129 
4130  graph_path_st_pcmw_extend(scip, graph, cost, TRUE, path, stvertextmp, &extensionstmp);
4131 
4132 
4133  for( int j = 0; j < nnodes; j++ )
4134  {
4135  if( graph->grad[j] == 0 || root == j )
4136  continue;
4137 
4138  if( Is_term(graph->term[j]) && !graph_pc_knotIsFixedTerm(graph, j) )
4139  continue;
4140 
4141  if( stvertextmp[j] )
4142  {
4143  assert(path[j].edge >= 0);
4144 
4145  newsolval += graph->cost[path[j].edge];
4146 
4147  if( Is_pseudoTerm(graph->term[j]) )
4148  {
4149  newsolval -= graph->prize[j];
4150  }
4151  }
4152  }
4153 
4154  /* new solution value better than old one? */
4155  if( SCIPisLT(scip, newsolval, bestsolval) )
4156  {
4157  SCIPdebugMessage("newsolval=%f bestsolval=%f \n", newsolval, bestsolval);
4158 
4159  extensions = TRUE;
4160  bestsolval = newsolval;
4161 
4162 #ifdef SCIP_DEBUG
4163  for( int k2 = 0; k2 < nnodes; k2++ )
4164  {
4165  if( !stvertex[k2] && stvertextmp[k2] )
4166  {
4167  printf("added \n");
4168  graph_knot_printInfo(graph, k2);
4169  graph_edge_printInfo(graph, path[k2].edge);
4170  }
4171  }
4172 #endif
4173 
4174 #ifndef NDEBUG
4175  objdiff += bestsolval - newsolval;
4176 #endif
4177 
4178  BMScopyMemoryArray(stvertex, stvertextmp, nnodes);
4179  BMScopyMemoryArray(orgpath, path, nnodes);
4180 
4181  /* save greedyextensions many best unconnected nodes */
4182  nextensions = 0;
4183 
4184  for( int j = 0; j < nnodes; j++ )
4185  if( !stvertex[j] && Is_pseudoTerm(graph->term[j]) && path[j].edge != UNKNOWN )
4186  SCIP_CALL( addToCandidates(scip, graph, path, j, greedyextensions, &nextensions, candidates, pqueue) );
4187 
4188  break;
4189  } /* if new solution value better than old one? */
4190  } /* if !stvertex[i] */
4191  } /* for l < nextension */
4192 
4193  /* no more extensions performed? */
4194  if( l == nextensions )
4195  break;
4196  } /* main loop */
4197 
4198  /* have vertices been added? */
4199  if( extensions )
4200  {
4201  assert(graph_pc_isPcMw(graph));
4202 
4203  SCIP_CALL( solstp_prune(scip, graph, stedge, stvertex) );
4204  }
4205 
4206  assert(solstp_isValid(scip, graph, stedge));
4207 
4208  SCIPpqueueFree(&pqueue);
4209  SCIPfreeBufferArray(scip, &path);
4210  SCIPfreeBufferArray(scip, &orgpath);
4211  SCIPfreeBufferArray(scip, &stvertextmp);
4212 
4213 #ifndef NDEBUG
4214  {
4215  const SCIP_Real newsolval = solstp_getObjBounded(graph, stedge, 0.0, nedges);
4216  SCIPdebugMessage("pcmw extend obj: initial=%f, new=%f \n", initialobj, newsolval);
4217  assert(SCIPisLE(scip, newsolval + objdiff, initialobj));
4218  }
4219 #endif
4220 
4221  return SCIP_OKAY;
4222 }
4223 
4224 /** Greedy Extension local heuristic for (R)PC and MW */
4226  SCIP* scip, /**< SCIP data structure */
4227  GRAPH* graph, /**< graph data structure */
4228  int* stedge, /**< initialized array to indicate whether an edge is part of the Steiner tree */
4229  STP_Bool* stvertex /**< uninitialized array to indicate whether a vertex is part of the Steiner tree */
4230  )
4231 {
4232  int candidates[GREEDY_EXTENSIONS];
4233  int ncandidates;
4234  DHEAP* dheap;
4235  STP_Bool* stvertextmp;
4236  SCIP_Real* dist;
4237  int* pred;
4238  const int nnodes = graph->knots;
4239  SCIP_Bool extensions = FALSE;
4240  int maxnode;
4241  const SCIP_Bool isexended = graph->extended;
4242 
4243 #ifndef NDEBUG
4244  const int nedges = graph->edges;
4245  const SCIP_Real initialobj = solstp_getObjBounded(graph, stedge, 0.0, nedges);
4246 #endif
4247 
4248  assert(scip && graph && stedge && stvertex);
4249 
4250  graph_pc_2orgcheck(scip, graph);
4251 
4252  solstp_setVertexFromEdge(graph, stedge, stvertex);
4253 
4254  /* compute candidates for extension */
4255 
4256  maxnode = -1;
4257  ncandidates = 0;
4258 
4259  for( int k = 0; k < nnodes; k++ )
4260  if( graph->mark[k] && !stvertex[k] && Is_term(graph->term[k]) && !graph_pc_termIsNonLeafTerm(graph, k) )
4261  {
4262  assert(graph->mark[k]);
4263 
4264  if( maxnode == -1 || graph->prize[k] > graph->prize[maxnode] )
4265  maxnode = k;
4266  }
4267 
4268  if( maxnode != -1 )
4269  {
4270  SCIP_RANDNUMGEN* randnumgen;
4271  int shift;
4272 
4273  SCIP_CALL( SCIPcreateRandom(scip, &randnumgen, 1, TRUE) );
4274 
4275  SCIP_CALL( SCIPallocBufferArray(scip, &dist, nnodes) );
4276  SCIP_CALL( SCIPallocBufferArray(scip, &pred, nnodes) );
4277  SCIP_CALL( SCIPallocBufferArray(scip, &stvertextmp, nnodes) );
4278 
4279  graph_heap_create(scip, nnodes, NULL, NULL, &dheap);
4280  graph_init_csr(scip, graph);
4281 
4282  shift = SCIPrandomGetInt(randnumgen, 0, nnodes - 1);
4283  ncandidates = 1;
4284  candidates[0] = maxnode;
4285 
4286  for( int k = 0; k < nnodes && ncandidates < GREEDY_EXTENSIONS; k++ )
4287  {
4288  const int node = (k + shift) % nnodes;
4289  if( graph->mark[k] && !stvertex[node] && Is_term(graph->term[node])
4290  && !graph_pc_termIsNonLeafTerm(graph, node) && node != maxnode )
4291  {
4292  assert(graph->mark[node]);
4293  candidates[ncandidates++] = node;
4294  }
4295  }
4296 
4297  SCIPfreeRandom(scip, &randnumgen);
4298  }
4299 
4300  /* main loop */
4301  for( int k = 0; k < ncandidates; k++ )
4302  {
4303  const int cand = candidates[k];
4304  SCIP_Bool success = FALSE;
4305 
4306  if( stvertex[cand] )
4307  {
4308  assert(k > 0);
4309  continue;
4310  }
4311 
4312  graph_path_st_pcmw_extendOut(scip, graph, cand, stvertex, dist, pred, stvertextmp, dheap, &success);
4313 
4314  if( success )
4315  extensions = TRUE;
4316  }
4317 
4318  /* have vertices been added? */
4319  if( extensions )
4320  {
4321  graph_pc_2trans(scip, graph);
4322 
4323  assert(graph_pc_isPcMw(graph));
4324 
4325  SCIP_CALL( solstp_prune(scip, graph, stedge, stvertex) );
4326  }
4327 
4328  if( maxnode != -1 )
4329  {
4330  graph_heap_free(scip, TRUE, TRUE, &dheap);
4331  graph_free_csr(scip, graph);
4332 
4333  SCIPfreeBufferArray(scip, &stvertextmp);
4334  SCIPfreeBufferArray(scip, &pred);
4335  SCIPfreeBufferArray(scip, &dist);
4336  }
4337 
4338 #ifndef NDEBUG
4339  assert(SCIPisLE(scip, solstp_getObjBounded(graph, stedge, 0.0, nedges), initialobj));
4340 #endif
4341 
4342  if( isexended && !graph->extended )
4343  graph_pc_2trans(scip, graph);
4344 
4345  if( !isexended && graph->extended )
4346  graph_pc_2org(scip, graph);
4347 
4348  return SCIP_OKAY;
4349 }
4350 
4351 
4352 /** creates the local primal heuristic and includes it in SCIP */
4354  SCIP* scip /**< SCIP data structure */
4355  )
4356 {
4357  SCIP_HEURDATA* heurdata;
4358  SCIP_HEUR* heur;
4359 
4360  /* create Local primal heuristic data */
4361  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
4362 
4363  /* include primal heuristic */
4364  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
4366  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecLocal, heurdata) );
4367 
4368  assert(heur != NULL);
4369 
4370  /* set non-NULL pointers to callback methods */
4371  SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyLocal) );
4372  SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeLocal) );
4373  SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolLocal) );
4374  SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolLocal) );
4375 
4376  /* add local primal heuristic parameters */
4377  SCIP_CALL( SCIPaddBoolParam(scip, "stp/duringroot",
4378  "should the heuristic be called during the root node?",
4379  &heurdata->duringroot, FALSE, DEFAULT_DURINGROOT, NULL, NULL) );
4380 
4381  SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/"HEUR_NAME"/maxfreq",
4382  "should the heuristic be executed at maximum frequeny?",
4383  &heurdata->maxfreq, FALSE, DEFAULT_MAXFREQLOC, NULL, NULL) );
4384 
4385  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/maxnsols",
4386  "maximum number of best solutions to improve",
4387  &heurdata->maxnsols, FALSE, DEFAULT_MAXNBESTSOLS, DEFAULT_NBESTSOLS_HARD, 100, NULL, NULL) );
4388 
4389  SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, DEFAULT_RANDSEED, TRUE) );
4390 
4391  return SCIP_OKAY;
4392 }
enum SCIP_Result SCIP_RESULT
Definition: type_result.h:52
SCIP_Real SCIPlinkcuttreeFindMinChainMw(SCIP *scip, const LCNODE *tree, const SCIP_Real *nodeweight, const int *heads, const int *stdeg, const SCIP_Bool *nodeIsBlocked, int start, int *first, int *last)
Definition: misc_stp.c:381
SCIP_RETCODE SCIPsetHeurExitsol(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEUREXITSOL((*heurexitsol)))
Definition: scip_heur.c:233
#define HEUR_USESSUBSCIP
Definition: heur_local.c:52
#define STP_Vectype(type)
Definition: stpvector.h:44
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
SCIP_Real SCIPlinkcuttreeFindMaxChain(SCIP *scip, const LCNODE *tree, const SCIP_Real *edgecosts, const SCIP_Real *prizes, const int *heads, const int *nonTermDeg, const SCIP_Bool *nodeIsBlocked, int start, int *first, int *last)
Definition: misc_stp.c:445
SCIP_RETCODE graph_init_csr(SCIP *, GRAPH *)
Definition: graph_util.c:1581
static SCIP_RETCODE solPrune(SCIP *scip, const GRAPH *graph, int *results)
Definition: heur_local.c:220
int graph_get_nEdges(const GRAPH *)
Definition: graph_stats.c:548
int SCIPpqueueNElems(SCIP_PQUEUE *pqueue)
Definition: misc.c:1468
int nprizemarks
Definition: heur_local.c:162
#define HEUR_DESC
Definition: heur_local.c:44
static SCIP_Bool probtypeIsValidForLocal(const GRAPH *graph)
Definition: heur_local.c:401
static SCIP_RETCODE solAddTry(SCIP *scip, SCIP_SOL **sols, SCIP_HEUR *heur, const GRAPH *graph, SCIP_Real initialsol_obj, SCIP_SOL *initialsol, int *results, SCIP_RESULT *result)
Definition: heur_local.c:312
#define StpVecGetSize(vec)
Definition: stpvector.h:139
SCIP_RETCODE SCIPStpunionfindInit(SCIP *scip, UF *uf, int length)
Definition: misc_stp.c:817
struct pcmw_data PCMW
void graph_pc_getBiased(const GRAPH *, SCIP_Real *RESTRICT, SCIP_Real *RESTRICT)
#define StpVecClear(vec)
Definition: stpvector.h:134
struct connectivity_data CONN
int *RESTRICT head
Definition: graphdefs.h:224
int source
Definition: graphdefs.h:195
SCIP_Bool graph_pc_isRootedPcMw(const GRAPH *)
static void vnoiDataReset(const CONN *connectData, const KPATHS *keypathsData, const int *graphmark, VNOILOC *vnoiData)
Definition: heur_local.c:2302
#define Is_term(a)
Definition: graphdefs.h:102
static void vnoiDataRestore(const CONN *connectData, const KPATHS *keypathsData, VNOILOC *vnoiData)
Definition: heur_local.c:2264
static SCIP_Bool solDegIsValid(SCIP *scip, const GRAPH *graph, const int *solDegree, const LCNODE *linkcutNodes)
Definition: heur_local.c:2354
SCIP_Real mstcost
Definition: heur_local.c:149
Constraint handler for Steiner problems.
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
signed int edge
Definition: graphdefs.h:287
SCIP_Bool graph_valid(SCIP *, const GRAPH *)
Definition: graph_base.c:1480
void SCIPStpunionfindClear(SCIP *scip, UF *uf)
Definition: misc_stp.c:841
int terms
Definition: graphdefs.h:192
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
int *const solDegreeNonTerm
Definition: heur_local.c:174
void graph_path_exec(SCIP *, const GRAPH *, int, int, const SCIP_Real *, PATH *)
Definition: graph_path.c:541
void SCIPpairheapMeldheaps(SCIP *scip, PHNODE **root1, PHNODE **root2, int *sizeroot1, int *sizeroot2)
Definition: misc_stp.c:708
#define GREEDY_EXTENSIONS
Definition: heur_local.c:69
includes methods for Steiner tree problem solutions
SCIP_RETCODE SCIPStpHeurLocalRunFast(SCIP *scip, GRAPH *graph, int *solEdges)
Definition: heur_local.c:3911
SCIP_RETCODE SCIPStpIncludeHeurLocal(SCIP *scip)
Definition: heur_local.c:4353
#define HEUR_PRIORITY
Definition: heur_local.c:46
void graph_knot_printInfo(const GRAPH *, int)
Definition: graph_stats.c:151
static SCIP_DECL_HEURCOPY(heurCopyLocal)
Definition: heur_local.c:3607
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: graph_base.c:796
void graph_pc_2org(SCIP *, GRAPH *)
SCIP_RETCODE graph_path_init(SCIP *, GRAPH *)
Definition: graph_path.c:480
#define EAT_FREE
Definition: graphdefs.h:57
SCIP_SOL ** SCIPgetSols(SCIP *scip)
Definition: scip_sol.c:2254
static SCIP_Real getEdgeCostUnbiased(const GRAPH *graph, const PCMW *pcmwData, int edge)
Definition: heur_local.c:897
#define FALSE
Definition: def.h:87
int *RESTRICT inpbeg
Definition: graphdefs.h:202
static void supergraphDataRestore(const GRAPH *graph, SGRAPH *supergraphData)
Definition: heur_local.c:1791
static void vnoiDataRepairPreprocess(SCIP *scip, const GRAPH *graph, const KPATHS *keypathsData, const CONN *connectData, const PCMW *pcmwData, VNOILOC *vnoiData, int *nheapelems)
Definition: heur_local.c:2197
int *RESTRICT path_state
Definition: graphdefs.h:256
SCIP_Bool isActive
Definition: heur_local.c:163
void graph_pathHeapAdd(const PATH *, int, int *RESTRICT, int *RESTRICT, int *)
static SCIP_Real getBoundaryPathCost(const GRAPH *graph, const VNOILOC *vnoiData, const PCMW *pcmwData, int boundaryedge)
Definition: heur_local.c:937
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
Problem data for stp problem.
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_RETCODE SCIPStpHeurLocalExtendPcMwOut(SCIP *scip, GRAPH *graph, int *stedge, STP_Bool *stvertex)
Definition: heur_local.c:4225
static void markSolTreeNodes(SCIP *scip, const GRAPH *graph, const int *solEdges, LCNODE *linkcutNodes, STP_Bool *solNodes)
Definition: heur_local.c:455
static void initSolNumberBounds(SCIP *scip, SCIP_HEURDATA *heurdata)
Definition: heur_local.c:290
#define HEUR_MAXDEPTH
Definition: heur_local.c:49
void graph_voronoiRepairMult(SCIP *, const GRAPH *, const SCIP_Real *, const STP_Bool *, int *RESTRICT, int *RESTRICT, int *RESTRICT, int *RESTRICT, UF *RESTRICT, PATH *RESTRICT)
void * SCIPpqueueFirst(SCIP_PQUEUE *pqueue)
Definition: misc.c:1454
static void insertionDecrementSolDegree(const GRAPH *graph, int node, INSERT *insertData)
Definition: heur_local.c:2460
int SCIPprobdataGetNVars(SCIP *scip)
static void insertionGetCandidateEdges(SCIP *scip, SCIP_HEURDATA *heurdata, const GRAPH *graph, const STP_Bool *solNodes, int vertex, int *insertcands, int *ncands)
Definition: heur_local.c:2841
int *const cutedgesEnd
Definition: heur_local.c:177
struct SCIP_HeurData SCIP_HEURDATA
Definition: type_heur.h:67
int *const solEdges
Definition: heur_local.c:135
SCIP_RETCODE SCIPpairheapDeletemin(SCIP *scip, int *element, SCIP_Real *key, PHNODE **root, int *size)
Definition: misc_stp.c:671
static void insertionResetBlockedNodes(INSERT *insertData)
Definition: heur_local.c:2652
SCIP_RETCODE SCIPincludeHeurBasic(SCIP *scip, SCIP_HEUR **heur, const char *name, const char *desc, char dispchar, int priority, int freq, int freqofs, int maxdepth, SCIP_HEURTIMING timingmask, SCIP_Bool usessubscip, SCIP_DECL_HEUREXEC((*heurexec)), SCIP_HEURDATA *heurdata)
Definition: scip_heur.c:108
#define EAT_LAST
Definition: graphdefs.h:58
#define STP_SPG
Definition: graphdefs.h:42
void SCIPpqueueFree(SCIP_PQUEUE **pqueue)
Definition: misc.c:1263
#define SCIPdebugMessage
Definition: pub_message.h:87
int SCIPrandomGetInt(SCIP_RANDNUMGEN *randnumgen, int minrandval, int maxrandval)
Definition: misc.c:10003
int *const supernodes
Definition: heur_local.c:145
struct supergraph_data SGRAPH
#define FARAWAY
Definition: graphdefs.h:89
int *const addedEdges
Definition: heur_local.c:175
static SCIP_Bool pcmwDataIsClean(const GRAPH *graph, const PCMW *pcmwData)
Definition: heur_local.c:750
void SCIPpairheapFree(SCIP *scip, PHNODE **root)
Definition: misc_stp.c:736
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static void insertionFinalizeReplacement(const GRAPH *graph, int v_insert, LCNODE *linkcutNodes, INSERT *insertData)
Definition: heur_local.c:2601
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
int number
Definition: misc_stp.h:55
void SCIPheurSetData(SCIP_HEUR *heur, SCIP_HEURDATA *heurdata)
Definition: heur.c:1362
void graph_voronoiRepair(SCIP *, const GRAPH *, const SCIP_Real *, const SCIP_Real *, int *, int *, PATH *, int *, int, UF *)
Definition: graph_vnoi.c:1100
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:74
header only, simple implementation of an STL like vector
GRAPH * SCIPprobdataGetGraph(SCIP_PROBDATA *probdata)
int *RESTRICT mark
Definition: graphdefs.h:198
SCIP_Bool graph_pc_termIsNonLeafTerm(const GRAPH *, int)
static const SCIP_Real * getEdgeCosts(const GRAPH *graph, const PCMW *pcmwData)
Definition: heur_local.c:1033
void SCIPlinkcuttreeEvert(LCNODE *RESTRICT tree, int root_new)
Definition: misc_stp.c:559
static SCIP_RETCODE localKeyVertexHeuristics(SCIP *scip, GRAPH *graph, SCIP_Bool useFast, STP_Bool *solNodes, LCNODE *linkcutNodes, int *solEdges, SCIP_Bool *success)
Definition: heur_local.c:3131
#define HEUR_FREQ
Definition: heur_local.c:47
static SCIP_RETCODE pcmwUpdate(SCIP *scip, GRAPH *graph, SOLTREE *soltreeData, PCMW *pcmwData)
Definition: heur_local.c:804
#define DEFAULT_DURINGROOT
Definition: heur_local.c:54
#define STP_RSMT
Definition: graphdefs.h:49
int *RESTRICT oeat
Definition: graphdefs.h:231
STP_Bool *const solNodes
Definition: heur_local.c:178
int *const newedges
Definition: heur_local.c:138
void SCIPlinkcuttreeInitNode(LCNODE *v)
Definition: misc_stp.c:347
SCIP_RETCODE SCIPsetHeurInitsol(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURINITSOL((*heurinitsol)))
Definition: scip_heur.c:217
void graph_pc_2trans(SCIP *, GRAPH *)
const char * SCIPheurGetName(SCIP_HEUR *heur)
Definition: heur.c:1441
SCIP_Bool extended
Definition: graphdefs.h:267
SCIP_HEUR * SCIPfindHeur(SCIP *scip, const char *name)
Definition: scip_heur.c:249
int stp_type
Definition: graphdefs.h:264
SCIP_Bool *const nodeIsBlocked
Definition: heur_local.c:179
static SCIP_RETCODE getLowestCommonAncestors(SCIP *scip, const GRAPH *graph, const VNOILOC *vnoiData, const SOLTREE *soltreeData, CONN *connectData)
Definition: heur_local.c:600
SCIP_RETCODE SCIPsetHeurFree(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURFREE((*heurfree)))
Definition: scip_heur.c:169
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_DECL_HEURINITSOL(heurInitsolLocal)
Definition: heur_local.c:3643
void graph_pc_2transcheck(SCIP *, GRAPH *)
static void insertionRestoreTree(const GRAPH *graph, const int *insertCands, int lcvertex_insert, LCNODE *linkcutNodes, INSERT *insertData)
Definition: heur_local.c:2711
static STP_Bool nodeIsCrucial(const GRAPH *graph, const int *steineredges, int node)
Definition: heur_local.c:863
STP_Bool *const nodeIsScanned
Definition: heur_local.c:137
SCIP_Real SCIPprobdataGetOffset(SCIP *scip)
SCIP_Bool SCIPprobdataProbIsAdversarial(SCIP *scip)
static void getKeyPathsStar(int keyvertex, const GRAPH *graph, const CONN *connectData, const SOLTREE *soltreeData, const PCMW *pcmwData, KPATHS *keypathsData, SGRAPH *supergraphData, SCIP_Bool *success)
Definition: heur_local.c:1992
static SCIP_Real getKeyPathReplaceCost(SCIP *scip, const GRAPH *graph, const VNOILOC *vnoiData, const SOLTREE *soltreeData, SCIP_Real edgecost_old_in, int boundedge_old, PCMW *pcmwData, int *boundedge_new)
Definition: heur_local.c:1732
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:128
SCIP_Real dist
Definition: misc_stp.h:56
SCIP_Real * prize
Definition: graphdefs.h:210
SCIP_RETCODE SCIPprobdataAddNewSol(SCIP *scip, SCIP_Real *nval, SCIP_HEUR *heur, SCIP_Bool *success)
#define GREEDY_EXTENSIONS_MW
Definition: heur_local.c:68
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
SCIP_Real dist
Definition: graphdefs.h:286
int *RESTRICT grad
Definition: graphdefs.h:201
SCIP_RETCODE SCIPpairheapInsert(SCIP *scip, PHNODE **root, PHNODE *pheapelems, int element, SCIP_Real key, int *size)
Definition: misc_stp.c:636
struct insertion_data INSERT
void graph_knot_add(GRAPH *, int)
Definition: graph_node.c:64
#define HEUR_TIMING
Definition: heur_local.c:50
static SCIP_RETCODE localVertexInsertion(SCIP *scip, SCIP_HEURDATA *heurdata, const GRAPH *graph, STP_Bool *solNodes, LCNODE *linkcutNodes, int *solEdges)
Definition: heur_local.c:2890
void graph_path_st_pcmw_extend(SCIP *, const GRAPH *, const SCIP_Real *, SCIP_Bool, PATH *, STP_Bool *, SCIP_Bool *)
Definition: graph_path.c:1706
#define DEFAULT_MAXFREQLOC
Definition: heur_local.c:55
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_HEUR * SCIPsolGetHeur(SCIP_SOL *sol)
Definition: sol.c:2604
SCIP_RETCODE SCIPStpHeurLocalRun(SCIP *scip, GRAPH *graph, int *solEdges)
Definition: heur_local.c:3899
#define STP_RMWCSP
Definition: graphdefs.h:54
static SCIP_RETCODE soltreeExchangeKeyPath(SCIP *scip, GRAPH *graph, const CONN *connectData, const VNOILOC *vnoiData, const KPATHS *keypathsData, const int *dfstree, const STP_Bool *scanned, int dfstree_pos, int boundedge_new, SOLTREE *soltreeData)
Definition: heur_local.c:1385
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
void SCIPlinkcuttreeCutNode(LCNODE *v)
Definition: misc_stp.c:372
void SCIPStpunionfindUnion(UF *uf, int p, int q, SCIP_Bool compress)
Definition: misc_stp.c:912
static SCIP_Real pcmwGetNewEdgePrize(const GRAPH *graph, const STP_Bool *steinertree, int edge, PCMW *pcmwData)
Definition: heur_local.c:716
int * term2edge
Definition: graphdefs.h:208
#define MST_MODE
Definition: graphdefs.h:98
int SCIPStpunionfindFind(UF *uf, int element)
Definition: misc_stp.c:885
static void pcmwDataClean(const GRAPH *graph, PCMW *pcmwData)
Definition: heur_local.c:780
#define STP_PCSPG
Definition: graphdefs.h:44
#define STP_OARSMT
Definition: graphdefs.h:50
void * SCIPpqueueRemove(SCIP_PQUEUE *pqueue)
Definition: misc.c:1434
STP_Bool *const nodeIsPinned
Definition: heur_local.c:136
Improvement heuristic for Steiner problems.
SCIP_RETCODE SCIPStpValidateSol(SCIP *, const GRAPH *, const double *, SCIP_Bool, SCIP_Bool *)
Definition: validate.c:202
unsigned char STP_Bool
Definition: portab.h:34
static SCIP_DECL_HEURFREE(heurFreeLocal)
Definition: heur_local.c:3621
int GNODECmpByDist(void *first_arg, void *second_arg)
const SCIP_Real * edgecosts
Definition: heur_local.c:173
SCIP_RETCODE solstp_pruneFromEdges(SCIP *scip, const GRAPH *g, int *result)
Definition: solstp.c:1432
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
SCIP_Real solstp_getObjBounded(const GRAPH *g, const int *soledge, SCIP_Real offset, int nedges)
Definition: solstp.c:1833
static SCIP_Bool solNeedsPruning(SCIP_SOL *initialsol)
Definition: heur_local.c:194
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
void graph_heap_free(SCIP *, SCIP_Bool, SCIP_Bool, DHEAP **)
Definition: graph_util.c:1034
static void solGetStpSol(SCIP *scip, const GRAPH *graph, SCIP_SOL *initialsol, int *results)
Definition: heur_local.c:258
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT ieat
Definition: graphdefs.h:230
int *RESTRICT path_heap
Definition: graphdefs.h:255
#define DEFAULT_RANDSEED
Definition: heur_local.c:62
int *RESTRICT tail
Definition: graphdefs.h:223
#define LOCAL_MAXRESTARTS_FAST
Definition: heur_local.c:64
#define DEFAULT_MINNBESTSOLS
Definition: heur_local.c:60
static void insertData(SCIP *scip, STP_Bitset termsmark, STP_Bitset rootsmark, int64_t nsubsets, int node_pos, int split_pos, int *subrootpos, DPSTREE *dpstree)
Definition: dpterms_util.c:172
void SCIPrandomPermuteIntArray(SCIP_RANDNUMGEN *randnumgen, int *array, int begin, int end)
Definition: misc.c:10044
#define flipedge(edge)
Definition: graphdefs.h:84
static SCIP_RETCODE localRun(SCIP *scip, GRAPH *graph, SCIP_Bool useFast, int *solEdges)
Definition: heur_local.c:3815
static void insertionInit(SCIP_HEURDATA *heurdata, const GRAPH *graph, const int *solEdges, int *solDegreeNonTerm, SCIP_Bool *nodeIsBlocked, int *vertices)
Definition: heur_local.c:2509
static SCIP_DECL_HEUREXEC(heurExecLocal)
Definition: heur_local.c:3689
void SCIPlinkcuttreeLink(LCNODE *tree, int v, int w, int edge)
Definition: misc_stp.c:357
#define LOCAL_MAXRESTARTS
Definition: heur_local.c:63
int SCIPgetNSols(SCIP *scip)
Definition: scip_sol.c:2205
#define DEFAULT_NELITESOLS
Definition: heur_local.c:59
int *RESTRICT term
Definition: graphdefs.h:196
static void insertionIncrementSolDegree(const GRAPH *graph, int node, INSERT *insertData)
Definition: heur_local.c:2486
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:127
void graph_edge_printInfo(const GRAPH *, int)
Definition: graph_stats.c:133
static SCIP_DECL_HEUREXITSOL(heurExitsolLocal)
Definition: heur_local.c:3669
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1435
SCIP_Bool graph_pc_knotIsNonLeafTerm(const GRAPH *, int)
void graph_path_st_pcmw_extendOut(SCIP *, const GRAPH *, int, STP_Bool *, SCIP_Real *, int *, STP_Bool *, DHEAP *, SCIP_Bool *)
Definition: graph_path.c:1051
static void soltreeUnmarkKpNodes(const GRAPH *graph, const KPATHS *keypathsData, SOLTREE *soltreeData)
Definition: heur_local.c:1337
const int * SCIPStpGetPcImplVerts(SCIP *scip)
Definition: cons_stp.c:1269
SCIP_RETCODE graph_trail_costAware(SCIP *, const GRAPH *, int, SCIP_Bool *)
Definition: graph_util.c:353
void graph_path_exit(SCIP *, GRAPH *)
Definition: graph_path.c:515
SCIP_Bool graph_pc_isPc(const GRAPH *)
#define CONNECT
Definition: graphdefs.h:87
int *const cutedgesStart
Definition: heur_local.c:176
int * chainStarts
Definition: heur_local.c:171
struct Voronoi_data_structures VNOILOC
static SCIP_Real getBoundaryPathCostPrized(const GRAPH *graph, const VNOILOC *vnoiData, const SOLTREE *soltreeData, int boundaryedge, PCMW *pcmwData)
Definition: heur_local.c:974
int *const prizemarklist
Definition: heur_local.c:161
static void getKeyPathUpper(SCIP *scip, int crucnode, const GRAPH *graph, const SOLTREE *soltreeData, const PCMW *pcmwData, CONN *connectData, KPATHS *keypathsData, SCIP_Bool *allGood)
Definition: heur_local.c:1201
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
SCIP_RETCODE SCIPpqueueCreate(SCIP_PQUEUE **pqueue, int initsize, SCIP_Real sizefac, SCIP_DECL_SORTPTRCOMP((*ptrcomp)), SCIP_DECL_PQUEUEELEMCHGPOS((*elemchgpos)))
Definition: misc.c:1236
static void insertionBlockChain(const GRAPH *graph, const LCNODE *lctree, int chainfirst_index, int chainlast_index, INSERT *insertData)
Definition: heur_local.c:2678
STP_Bool *const prizemark
Definition: heur_local.c:160
#define Is_pseudoTerm(a)
Definition: graphdefs.h:103
static SCIP_Bool solIsTrivialPcMw(const GRAPH *graph, const int *solEdges)
Definition: heur_local.c:499
static SCIP_RETCODE connectivityDataInit(SCIP *scip, const GRAPH *graph, const VNOILOC *vnoiData, const SOLTREE *soltreeData, const PCMW *pcmwData, CONN *connectData)
Definition: heur_local.c:1122
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2304
#define DEFAULT_NBESTSOLS_HARD
Definition: heur_local.c:58
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_Real getNewPrizeNode(const GRAPH *graph, const STP_Bool *steinertree, const SCIP_Real *prize_biased, const int *graphmark, int node, STP_Bool *prizemark, int *prizemarklist, int *prizemarkcount)
Definition: heur_local.c:659
struct SCIP_ProbData SCIP_PROBDATA
Definition: type_prob.h:44
int solroot
Definition: heur_local.c:164
SCIP_RETCODE SCIPpqueueInsert(SCIP_PQUEUE *pqueue, void *elem)
Definition: misc.c:1335
SCIP_PROBDATA * SCIPgetProbData(SCIP *scip)
Definition: scip_prob.c:962
void graph_pc_2orgcheck(SCIP *, GRAPH *)
static void insertionReplaceChain(const GRAPH *graph, int newedge, LCNODE *lctree, INSERT *insertData, int v_lc, int headCurr_lc, int chainfirst_index, int chainlast_index)
Definition: heur_local.c:2787
SCIP_Bool graph_pc_knotIsFixedTerm(const GRAPH *, int)
SCIP_Real * cost
Definition: graphdefs.h:221
static void dfsorder(const GRAPH *graph, const int *edges, int node, int *counter, int *dfst)
Definition: heur_local.c:633
#define DEFAULT_MINNBESTSOLS_HARD
Definition: heur_local.c:61
static int pcmwGetSolRoot(const GRAPH *graph, const int *solEdges)
Definition: heur_local.c:686
SCIP_RETCODE graph_heap_create(SCIP *, int, int *, DENTRY *, DHEAP **)
Definition: graph_util.c:992
struct solution_tree_data SOLTREE
SCIP_Real *const edgecost_org
Definition: heur_local.c:159
SCIP_Bool graph_pc_isPcMw(const GRAPH *)
#define GREEDY_MAXRESTARTS
Definition: heur_local.c:67
SCIP_Real * memvdist
Definition: heur_local.c:94
struct keypaths_data_structures KPATHS
#define HEUR_NAME
Definition: heur_local.c:43
const int * SCIPStpGetPcImplStarts(SCIP *scip)
Definition: cons_stp.c:1231
#define SCIP_Real
Definition: def.h:177
SCIP_RETCODE SCIPStpHeurLocalExtendPcMwImp(SCIP *scip, const GRAPH *graph, int *result)
Definition: heur_local.c:3923
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:694
#define StpVecFree(scip, vec)
Definition: stpvector.h:149
SCIP_RETCODE SCIPpairheapBuffarr(SCIP *scip, const PHNODE *root, int size, int **elements)
Definition: misc_stp.c:761
void graph_free_csr(SCIP *, GRAPH *)
Definition: graph_util.c:1623
void graph_voronoi(const GRAPH *, const SCIP_Real *, const SCIP_Real *, const STP_Bool *, int *, PATH *)
Definition: graph_vnoi.c:333
SCIP_Real * SCIPprobdataGetXval(SCIP *scip, SCIP_SOL *sol)
int *RESTRICT outbeg
Definition: graphdefs.h:204
#define DEFAULT_NBESTSOLS
Definition: heur_local.c:57
shortest paths based primal heuristics for Steiner problems
STP_Bool *const nodeIsLowSupernode
Definition: heur_local.c:146
int edges
Definition: graphdefs.h:219
SCIP_Real *const prize_biased
Definition: heur_local.c:157
static SCIP_Bool solNodeIsValid(SCIP *scip, const GRAPH *graph, const STP_Bool *solNodes, const LCNODE *linkcutNodes)
Definition: heur_local.c:2411
void graph_pc_getOrgCosts(SCIP *, const GRAPH *, SCIP_Real *)
#define STP_GSTP
Definition: graphdefs.h:53
void solstp_setVertexFromEdge(const GRAPH *g, const int *result, STP_Bool *solnode)
Definition: solstp.c:2078
int graph_pc_nNonFixedTerms(const GRAPH *)
SCIP_RETCODE SCIPStpHeurLocalExtendPcMw(SCIP *scip, GRAPH *graph, const SCIP_Real *cost, int *stedge, STP_Bool *stvertex)
Definition: heur_local.c:3991
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
static void insertionInitInsert(SCIP *scip, const GRAPH *graph, int v_insert, int initialEdge, LCNODE *linkcutNodes, INSERT *insertData, SCIP_Real *diff)
Definition: heur_local.c:2548
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define UNKNOWN
Definition: sepa_mcf.c:4095
static SCIP_RETCODE lca(SCIP *scip, const GRAPH *graph, int u, UF *uf, STP_Bool *nodesmark, const int *steineredges, STP_Vectype(int) *lcalists, PHNODE **boundpaths, const int *heapsize, const int *vbase)
Definition: heur_local.c:543
SCIP_RETCODE SCIPsetHeurCopy(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURCOPY((*heurcopy)))
Definition: scip_heur.c:153
#define nnodes
Definition: gastrans.c:65
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:163
SCIP_RETCODE graph_init(SCIP *, GRAPH **, int, int, int)
Definition: graph_base.c:607
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:123
SCIP_Bool SCIPStpunionfindIsClear(SCIP *scip, const UF *uf)
Definition: misc_stp.c:861
SCIP_Bool solstp_isValid(SCIP *scip, const GRAPH *graph, const int *result)
Definition: solstp.c:1650
STP_Bool *const solNodes
Definition: heur_local.c:133
#define SCIP_CALL_ABORT(x)
Definition: def.h:363
SCIP_HEURDATA * SCIPheurGetData(SCIP_HEUR *heur)
Definition: heur.c:1352
#define HEUR_FREQOFS
Definition: heur_local.c:48
SCIP_Real *const edgecost_biased
Definition: heur_local.c:158
SCIP_RETCODE solstp_prune(SCIP *scip, const GRAPH *g, int *result, STP_Bool *connected)
Definition: solstp.c:1366
LCNODE *const linkcutNodes
Definition: heur_local.c:134
void SCIPStpunionfindFreeMembers(SCIP *scip, UF *uf)
Definition: misc_stp.c:955
static SCIP_RETCODE supergraphComputeMst(SCIP *scip, const GRAPH *graph, const CONN *connectData, const VNOILOC *vnoiData, const KPATHS *keypathsData, int crucnode, SOLTREE *soltreeData, PCMW *pcmwData, SGRAPH *supergraphData)
Definition: heur_local.c:1813
#define HEUR_DISPCHAR
Definition: heur_local.c:45
#define STP_RPCSPG
Definition: graphdefs.h:45
int *const blockedList
Definition: heur_local.c:180
#define DEFAULT_MAXNBESTSOLS
Definition: heur_local.c:56
#define STP_MWCSP
Definition: graphdefs.h:51
static SCIP_RETCODE addToCandidates(SCIP *scip, const GRAPH *graph, const PATH *path, int i, int greedyextensions, int *nextensions, GNODE *candidates, SCIP_PQUEUE *pqueue)
Definition: heur_local.c:415
int SCIPsolGetIndex(SCIP_SOL *sol)
Definition: sol.c:2635
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:48
static SCIP_RETCODE soltreeElimKeyPathsStar(SCIP *scip, const GRAPH *graph, const CONN *connectData, const VNOILOC *vnoiData, const KPATHS *keypathsData, const SGRAPH *supergraphData, const int *dfstree, const STP_Bool *scanned, int dfstree_pos, SOLTREE *soltreeData)
Definition: heur_local.c:1534
static void soltreeMarkKpNodes(const GRAPH *graph, const KPATHS *keypathsData, SOLTREE *soltreeData)
Definition: heur_local.c:1361
static SCIP_RETCODE connectivityDataKeyElimUpdate(SCIP *scip, const GRAPH *graph, const VNOILOC *vnoiData, const SGRAPH *supergraphData, int crucnode, CONN *connectData)
Definition: heur_local.c:1046