# SCIP

Solving Constraint Integer Programs

graph_tpath.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 /* */
7 /* fuer Informationstechnik Berlin */
8 /* */
10 /* */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file graph_tpath.c
17  * @brief Graph algorithms for Steiner problems related to shortest paths to terminals.
18  * @author Thorsten Koch
19  * @author Daniel Rehfeldt
20  *
21  * This file encompasses various algorithms that are used to compute shortest paths from
22  * non-terminals to terminals.
23  *
24  */
25
26 #include <stdlib.h>
27 #include <stdio.h>
28 #include <stddef.h>
29 #include <assert.h>
30 #include "graph.h"
31 #include "graphdefs.h"
32 #include "stpvector.h"
33 #include "reduce.h"
34
35 #define STP_TPATHS_NTERMBASES 4
36 #define STP_TPATHS_RESERVESIZE 16
37
38 /** reset single path struct */
39 #define pathResetSingle(entry) \
40  do \
41  { \
42  entry.dist = FARAWAY; \
43  entry.edge = UNKNOWN; \
44  } while( 0 )
45
46
47 /** resets node */
48 #define tpathsRepairResetNode(scip, node, resetnodes, stack, nodes_isvisited) \
49  do \
50  { \
51  assert(scip && resetnodes && stack && nodes_isvisited); \
52  assert(!nodes_isvisited[node]); \
53  StpVecPushBack(scip, resetnodes, node); \
54  StpVecPushBack(scip, stack, node); \
55  nodes_isvisited[node] = TRUE; \
56  } while( 0 )
57
58
59 /** Steiner nodes to terminal paths
60  * NOTE: all arrays are of size STP_TPATHS_NTERMBASES * nnodes */
62 {
63  PATH* termpaths; /**< path data (leading to first, second, ... terminal) */
64  int* termbases; /**< terminals to each non terminal */
65  int* termbases_org; /**< original terminals to each non terminal; only needed for repair */
66  int* state; /**< array to mark the state of each node during calculation */
67  int nnodes; /**< number of nodes of underlying graph */
68 };
69
70
71 /** data needed to repair node to terminal paths for edge-deletion */
72 typedef struct tpaths_repair
73 {
74  STP_Vectype(int)* resetnodes; /**< nodes to be reseted; for different levels (1st to 4th) */
75  STP_Vectype(int) stack; /**< temporary vector */
76  TPATHS* tpaths; /**< the terminal paths */
77  SCIP_Bool* nodes_isvisited; /**< visited nodes during repair process */
78  int edge; /**< edge about to be eliminated */
79  int nHeapElems; /**< number of heap elements */
80  SCIP_Bool withEdgeReplacement;/**< with edge replacement? */
81 } TREPAIR;
82
83
84
85
86 /*
87  * Local methods
88  */
89
90
91
92 /*---------------------------------------------------------------------------*/
93 /*--- Name : get NEAREST knot ---*/
94 /*--- Function : Holt das oberste Element vom Heap (den Knoten mit der ---*/
95 /*--- geringsten Entfernung zu den bereits Verbundenen) ---*/
96 /*--- Parameter: Derzeitige Entfernungen und benutzte Kanten ---*/
97 /*--- Returns : Nummer des bewussten Knotens ---*/
98 /*---------------------------------------------------------------------------*/
99 inline static
101  int* RESTRICT heap,
102  int* RESTRICT state,
103  int* RESTRICT count, /* pointer to store the number of elements on the heap */
104  const PATH* path)
105 {
106  int k;
107  int t;
108  int c;
109  int j;
110
111  /* Heap shift down
112  * (Oberstes Element runter und korrigieren)
113  */
114  k = heap[1];
115  j = 1;
116  c = 2;
117  heap[1] = heap[(*count)--];
118  state[heap[1]] = 1;
119
120  if ((*count) > 2)
121  if (LT(path[heap[3]].dist, path[heap[2]].dist))
122  c++;
123
124  while((c <= (*count)) && GT(path[heap[j]].dist, path[heap[c]].dist))
125  {
126  t = heap[c];
127  heap[c] = heap[j];
128  heap[j] = t;
129  state[heap[j]] = j;
130  state[heap[c]] = c;
131  j = c;
132  c += c;
133
134  if ((c + 1) <= (*count))
135  if (LT(path[heap[c + 1]].dist, path[heap[c]].dist))
136  c++;
137  }
138  return(k);
139 }
140
141
142
143 /** set new distance and predeccesor */
144 inline static
146  int* RESTRICT heap,
147  int* RESTRICT state,
148  int* RESTRICT count, /* pointer to store the number of elements on the heap */
149  PATH* RESTRICT path,
150  int l,
151  int k,
152  int edge,
153  SCIP_Real newdist
154  )
155 {
156  int t;
157  int c;
158  int j;
159
160  path[l].dist = newdist;
161  path[l].edge = edge;
162
163  /* new node? */
164  if( state[l] == UNKNOWN )
165  {
166  heap[++(*count)] = l;
167  state[l] = (*count);
168  }
169
170  /* Heap shift up */
171  j = state[l];
172  c = j / 2;
173  while( (j > 1) && path[heap[c]].dist > path[heap[j]].dist )
174  {
175  t = heap[c];
176  heap[c] = heap[j];
177  heap[j] = t;
178  state[heap[j]] = j;
179  state[heap[c]] = c;
180  j = c;
181  c = j / 2;
182  }
183 }
184
185
186 /** terminal to terminal distance */
187 inline
188 static void utdist(
189  SCIP* scip,
190  const GRAPH* g,
191  PATH* path,
192  SCIP_Real ecost,
193  int* vbase,
194  int k,
195  int l,
196  int k2,
197  int shift,
198  int nnodes
199  )
200 {
201  SCIP_Real dist;
202  int vbk;
203  int vbk2;
204
205  if( Is_term(g->term[k]) )
206  vbk = k;
207  else
208  vbk = vbase[k];
209
210  if( l == 0 )
211  {
212  assert(shift == 0);
213
214  dist = ecost;
215  if( !Is_term(g->term[k]) )
216  dist += path[k].dist;
217
218  if( !Is_term(g->term[k2]) )
219  {
220  dist += path[k2].dist;
221  vbk2 = vbase[k2];
222  }
223  else
224  {
225  vbk2 = k2;
226  }
227
228  if( SCIPisLT(scip, dist, path[vbk].dist) )
229  {
230  path[vbk].dist = dist;
231  vbase[vbk] = vbk2;
232  return;
233  }
234  }
235  else
236  {
237  int max;
238  int pos;
239  int r;
240  int s;
241  int t;
242
243  pos = vbk + shift;
244  max = MIN((l + 1), 3);
245
246  for( r = 0; r <= max; r++ )
247  {
248  if( Is_term(g->term[k2]) )
249  {
250  if( r == 0 )
251  t = k2;
252  else
253  break;
254  }
255  else
256  {
257  t = vbase[k2 + (r * nnodes)];
258  }
259  for( s = 0; s < l; s++ )
260  if( vbase[vbk + s * nnodes] == t )
261  break;
262  if( s < l || vbk == t )
263  continue;
264
265  dist = ecost;
266  if( !Is_term(g->term[k]) )
267  dist += path[k].dist;
268  if( !Is_term(g->term[k2]) )
269  dist += path[k2 + (r * nnodes)].dist;
270
271  if( SCIPisLT(scip, dist, path[pos].dist) )
272  {
273  path[pos].dist = dist;
274  vbase[pos] = t;
275  return;
276  }
277  }
278  }
279  return;
280 }
281
282
283 /** prints terminal path from given node */
284 static
286  const GRAPH* g, /**< graph data structure */
287  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
288  const TPATHS* tpaths, /**< storage for terminal paths */
289  int node, /**< node to start from */
290  int level /**< between 1 and 4 */
291 )
292 {
293  const int* const termbases = tpaths->termbases;
294  const int nnodes = graph_get_nNodes(g);
295  const int offset = (level - 1) * nnodes;
296  const int base = termbases[node + offset];
297  const PATH* const termpaths = tpaths->termpaths;
298  int node_pred;
299  int node_curr;
300  int node_next;
301  int edge_curr;
302  assert(1 <= level && level <= STP_TPATHS_NTERMBASES);
303  assert(graph_knot_isInRange(g, base));
304  assert(base != UNKNOWN);
305
306  node_pred = -1;
307  node_curr = node;
308  edge_curr = termpaths[node_curr + offset].edge;
309  node_next = (edge_curr >= 0) ? g->tail[edge_curr] : -1;
310
311  while( node_curr != base )
312  {
313  assert(graph_knot_isInRange(g, node_curr));
314
315  if( !graph_edge_isInRange(g, edge_curr) )
316  {
317  printf("found deleted edge, break \n");
318  break;
319  }
320
321  if( sdprofit && node_curr != node )
322  {
323  const SCIP_Real profit = reduce_sdprofitGetProfit(sdprofit, node_curr, node_pred, node_next);
324  printf("node %d profit: %f \n", node_curr, profit);
325  }
326
327  graph_edge_printInfo(g, edge_curr);
328
329  if( graph_edge_isDeleted(g, edge_curr) )
330  {
331  printf("...edge killed! \n");
332  }
333
334  node_pred = node_curr;
335  node_curr = node_next;
336  edge_curr = termpaths[node_curr + offset].edge;
337  node_next = (edge_curr >= 0) ? g->tail[edge_curr] : -1;
338  }
339 }
340
341
342 /** allocates TPATHS data */
343 static
345  SCIP* scip, /**< SCIP */
346  const GRAPH* g, /**< graph */
347  TPATHS** tpaths /**< the terminal paths */
348 )
349 {
350
351  TPATHS* tp;
352  const int nnodes = graph_get_nNodes(g);
353  assert(nnodes >= 1);
354
355  SCIP_CALL( SCIPallocMemory(scip, tpaths) );
356
357  tp = *tpaths;
358  tp->nnodes = nnodes;
359
361  SCIP_CALL( SCIPallocMemoryArray(scip, &(tp->termbases), nnodes * STP_TPATHS_NTERMBASES) );
362  SCIP_CALL( SCIPallocMemoryArray(scip, &(tp->state), nnodes * STP_TPATHS_NTERMBASES) );
363  tp->termbases_org = NULL;
364
365  return SCIP_OKAY;
366 }
367
368
369 /** resets TPATHS members */
370 static
372  const GRAPH* g, /**< graph */
373  int level, /**< level between 2 and 4 */
374  TPATHS* tpaths /**< the terminal paths */
375 )
376 {
377  PATH* RESTRICT path = tpaths->termpaths;
378  int* RESTRICT vbase = tpaths->termbases;
379  int* RESTRICT state = tpaths->state;
380  const int nnodes = graph_get_nNodes(g);
381  const int offset = (level - 1) * nnodes;
382
383  assert(2 <= level && level <= 4);
384  assert(offset == nnodes || offset == 2 * nnodes || offset == 3 * nnodes);
385
386  for( int i = 0; i < nnodes; i++ )
387  {
388  /* copy of node i */
389  const int k = i + offset;
390
391  vbase[k] = UNKNOWN;
392  state[k] = UNKNOWN;
393  pathResetSingle(path[k]);
394  }
395
396  for( int i = 0; i < offset; i++ )
397  {
398  state[i] = CONNECT;
399  }
400 }
401
402
403 /** gets (up to) four close terminals to given node i */
404 static inline
406  const GRAPH* g, /**< graph */
407  const TPATHS* tpaths, /**< the terminal paths */
408  int maxnumber, /**< number of close terms */
409  int node, /**< node */
410  SCIP_Real maxdist, /**< maximum valid distance */
411  SCIP_Bool distIsStrict, /**< is 'maxdist' strict? */
412  int* RESTRICT closeterms, /**< four close terminals */
413  int* RESTRICT firstedges, /**< corresponding first edge (can be NULL) */
414  SCIP_Real* RESTRICT closeterms_dist, /**< four close terminal distance */
415  int* RESTRICT ncloseterms /**< number of close terminals found */
416 )
417 {
418  const PATH* const termpaths = tpaths->termpaths;
419  const int* const termbases = tpaths->termbases;
420  const int nnodes = tpaths->nnodes;
421  int pos = node;
422  int nnterms = 0;
423  const SCIP_Bool withFirstEdges = (firstedges != NULL);
424
425  assert(closeterms && closeterms_dist && ncloseterms && g);
426  assert(maxnumber <= STP_TPATHS_NTERMBASES);
427  assert(0 <= node && node < nnodes);
428  assert(nnodes == g->knots);
429  assert(GE(maxdist, 0.0));
430
431  if( Is_term(g->term[node]))
432  {
433  *ncloseterms = 1;
434  closeterms[0] = node;
435  closeterms_dist[0] = 0.0;
436
437  if( withFirstEdges )
438  firstedges[0] = UNKNOWN;
439
440  return;
441  }
442
443  for( int k = 0; k < maxnumber; k++ )
444  {
446  distIsStrict ?
447  LT(termpaths[pos].dist, maxdist)
448  : LE(termpaths[pos].dist, maxdist);
449
451  {
452  closeterms[nnterms] = termbases[pos];
453
454  if( withFirstEdges )
455  {
456  firstedges[nnterms] = termpaths[pos].edge;
457  }
458  closeterms_dist[nnterms++] = termpaths[pos].dist;
459  }
460  else
461  {
462  assert(distIsStrict || !LE(termpaths[pos].dist, maxdist));
463  }
464
465  pos += nnodes;
466  }
467
468
469 #ifndef NDEBUG
470  for( int k = 0; k < nnterms; k++ )
471  {
472  const int term = closeterms[k];
473  assert(graph_knot_isInRange(g, term));
474  assert(Is_term(g->term[term]) || g->grad[term] == 0);
475  }
476 #endif
477
478  *ncloseterms = nnterms;
479 }
480
481
482 /** gets new distance for extension */
483 inline static
485  const GRAPH* g, /**< graph */
486  const SCIP_Real* cost, /**< edge costs */
487  const SCIP_Real* costrev, /**< reversed edge costs */
488  const int* vbase, /**< bases */
489  int node, /**< node to extend from */
490  int nextnode, /**< node to extend to */
491  int edge, /**< the edge */
492  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
493  const PATH* path /**< the actual terminal paths */
494 )
495 {
496  SCIP_Real distnew;
497  const SCIP_Real ecost = ((g->source == vbase[node]) ? cost[edge] : costrev[edge]);
498
499 #ifndef NDEBUG
500  assert(node % g->knots == g->tail[edge]);
501  assert(nextnode % g->knots == g->head[edge]);
502 #endif
503
504  if( sdprofit && path[node].edge >= 0 )
505  {
506  const int nnodes = graph_get_nNodes(g);
507  const int node_pred = g->tail[path[node].edge];
508
509  assert(!Is_term(g->term[node % g->knots]));
510
511  if( node_pred == nextnode % nnodes )
512  {
513  distnew = path[node].dist + ecost;
514  }
515  else
516  {
517  distnew = reduce_sdprofitGetBiasedDist(sdprofit, node % nnodes,
518  ecost, path[node].dist, nextnode % nnodes, node_pred);
519  }
520  }
521  else
522  {
523  distnew = path[node].dist + ecost;
524  }
525
526  assert(graph_pc_isMw(g) || GE(distnew, 0.0));
527  assert(LE(distnew, path[node].dist + ecost));
528
529  return distnew;
530 }
531
532
533 /** allocates TPATHS data */
534 static inline
536  SCIP* scip, /**< SCIP */
537  GRAPH* g, /**< graph */
538  TPATHS* tpaths /**< the terminal paths */
539 )
540 {
541  assert(STP_TPATHS_NTERMBASES == 4);
542
543  graph_tpathsSetAll4(g, g->cost, g->cost, NULL, tpaths);
544 }
545
546 /** allocates TPATHS data */
547 static inline
549  const SDPROFIT* sdprofit, /**< SD profit */
550  GRAPH* g, /**< graph */
551  TPATHS* tpaths /**< the terminal paths */
552 )
553 {
554  assert(STP_TPATHS_NTERMBASES == 4);
555
556  graph_tpathsSetAll4(g, g->cost, g->cost, sdprofit, tpaths);
557 }
558
559
560 /** sub-method to find closest terminal to each non terminal (Voronoi diagram) */
561 static
563  SCIP* scip, /**< SCIP or NULL */
564  const GRAPH* g, /**< graph data structure */
565  const SCIP_Real* cost, /**< edge costs */
566  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
567  int heapsize, /**< size of heap */
568  TPATHS* tpaths, /**< storage for terminal paths */
569  TREPAIR* repair /**< data for repairing or NULL */
570 )
571 {
572  int nHeapElems = heapsize;
573  int* RESTRICT heap = g->path_heap;
574  PATH* RESTRICT path = tpaths->termpaths;
575  int* RESTRICT vbase = tpaths->termbases;
576  int* RESTRICT state = tpaths->state;
577  const SCIP_Bool withProfit = (sdprofit != NULL);
578  const SCIP_Bool extend = (repair != NULL);
579  SCIP_Bool* nodes_isvisited = extend ? repair->nodes_isvisited : NULL;
580
581  assert(nHeapElems >= 0);
582  assert(repair == NULL || repair->withEdgeReplacement);
583
584  while( nHeapElems > 0 )
585  {
586  const int k = pathheapGetNearest(heap, state, &nHeapElems, path);
587  const SCIP_Real k_dist = path[k].dist;
588
589
590
591  /* mark vertex k as scanned */
592  state[k] = CONNECT;
593
594  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
595  {
596  const int m = g->head[e];
597
598  if( state[m] != CONNECT && g->mark[m] )
599  {
600  SCIP_Real distnew;
601
602  if( withProfit && vbase[k] != k )
603  {
604  int k_pred;
605  assert(path[k].edge >= 0);
606  k_pred = g->tail[path[k].edge];
607  distnew = reduce_sdprofitGetBiasedDist(sdprofit, k, cost[e], k_dist, m, k_pred);
608  }
609  else
610  {
611  assert(!withProfit || Is_term(g->term[k]));
612  distnew = k_dist + cost[e];
613  }
614
615  assert(GE(distnew, 0.0) && LE(distnew, k_dist + cost[e]));
616
617  /* check whether the path (to m) including k is shorter than the so far best known */
618  if( GT(path[m].dist, distnew) )
619  {
620  if( extend && !nodes_isvisited[m] )
621  {
622  StpVecPushBack(scip, repair->resetnodes[0], m);
623  nodes_isvisited[m] = TRUE;
624  }
625
626  pathheapCorrectDist(heap, state, &nHeapElems, path, m, k, e, distnew);
627  vbase[m] = vbase[k];
628  }
629  }
630  }
631  }
632 }
633
634
635 /** sub-method to find 2nd closest terminal to each non terminal */
636 static
638  SCIP* scip, /**< SCIP or NULL */
639  const GRAPH* g, /**< graph data structure */
640  const SCIP_Real* cost, /**< edge costs */
641  const SCIP_Real* costrev, /**< reverse edge costs */
642  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
643  int heapsize, /**< size of heap */
644  TPATHS* tpaths, /**< storage for terminal paths */
645  TREPAIR* repair /**< data for repairing or NULL */
646 )
647 {
648  const int nnodes = graph_get_nNodes(g);
649  int nheapElems = heapsize;
650  int* RESTRICT heap = g->path_heap;
651  PATH* RESTRICT path = tpaths->termpaths;
652  int* RESTRICT vbase = tpaths->termbases;
653  int* RESTRICT state = tpaths->state;
654  const SCIP_Bool extend = (repair != NULL);
655  SCIP_Bool* nodes_isvisited = extend ? repair->nodes_isvisited : NULL;
656
657  assert(nheapElems >= 0);
658  assert(repair == NULL || repair->withEdgeReplacement);
659
660  while( nheapElems > 0 )
661  {
662  const int k2 = pathheapGetNearest(heap, state, &nheapElems, path);
663  const int k = k2 - nnodes;
664
665  assert(graph_knot_isInRange(g, k));
666
667  /* mark vertex k as removed from heap */
668  state[k2] = UNKNOWN;
669
670  /* iterate over all outgoing edges of vertex (ancestor of) k */
671  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
672  {
674
676  {
678  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, k2, head2, e, sdprofit, path);
679
680  /* check whether the path (to j) including k is shorter than the so far best known */
682  {
683  if( extend && !nodes_isvisited[head] )
684  {
687  }
688
689  pathheapCorrectDist(heap, state, &nheapElems, path, head2, k2, e, distnew);
691  }
692  }
693  }
694  }
695 }
696
697
698 /** sub-method to find 3rd closest terminal to each non terminal */
699 static
701  SCIP* scip, /**< SCIP or NULL */
702  const GRAPH* g, /**< graph data structure */
703  const SCIP_Real* cost, /**< edge costs */
704  const SCIP_Real* costrev, /**< reverse edge costs */
705  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
706  int heapsize, /**< size of heap */
707  TPATHS* tpaths, /**< storage for terminal paths */
708  TREPAIR* repair /**< data for repairing or NULL */
709 )
710 {
711  const int nnodes = graph_get_nNodes(g);
712  const int dnnodes = 2 * nnodes;
713  int nheapElems = heapsize;
714  int* RESTRICT heap = g->path_heap;
715  PATH* RESTRICT path = tpaths->termpaths;
716  int* RESTRICT vbase = tpaths->termbases;
717  int* RESTRICT state = tpaths->state;
718  const SCIP_Bool extend = (repair != NULL);
719  SCIP_Bool* nodes_isvisited = extend ? repair->nodes_isvisited : NULL;
720
721  assert(nheapElems >= 0);
722  assert(repair == NULL || repair->withEdgeReplacement);
723
724  while( nheapElems > 0 )
725  {
726  /* get the next (i.e. a nearest) vertex of the heap */
727  const int k3 = pathheapGetNearest(heap, state, &nheapElems, path);
728  const int k = k3 - dnnodes;
729  assert(0 <= k && k < nnodes);
730
731  /* mark vertex k as removed from heap */
732  state[k3] = UNKNOWN;
733
734  /* iterate over all outgoing edges of vertex k */
735  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
736  {
738
740  {
742  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, k3, head3, e, sdprofit, path);
743
744  /* check whether the path (to j) including k is shorter than the so far best known */
746  {
747  if( extend && !nodes_isvisited[head] )
748  {
751  }
752
753  pathheapCorrectDist(heap, state, &nheapElems, path, head3, k3, e, distnew);
755  }
756  }
757  }
758  }
759 }
760
761
762 /** sub-method to find 4th closest terminal to each non terminal */
763 static
765  SCIP* scip, /**< SCIP or NULL */
766  const GRAPH* g, /**< graph data structure */
767  const SCIP_Real* cost, /**< edge costs */
768  const SCIP_Real* costrev, /**< reverse edge costs */
769  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
770  int heapsize, /**< size of heap */
771  TPATHS* tpaths, /**< storage for terminal paths */
772  TREPAIR* repair /**< data for repairing or NULL */
773 )
774 {
775  const int nnodes = graph_get_nNodes(g);
776  const int dnnodes = 2 * nnodes;
777  const int tnnodes = 3 * nnodes;
778  int nHeapElems = heapsize;
779  int* RESTRICT heap = g->path_heap;
780  PATH* RESTRICT path = tpaths->termpaths;
781  int* RESTRICT vbase = tpaths->termbases;
782  int* RESTRICT state = tpaths->state;
783  const SCIP_Bool extend = (repair != NULL);
784  SCIP_Bool* nodes_isvisited = extend ? repair->nodes_isvisited : NULL;
785
786  assert(nHeapElems >= 0);
787  assert(repair == NULL || repair->withEdgeReplacement);
788
789  /* until the heap is empty */
790  while( nHeapElems > 0 )
791  {
792  /* get the next (i.e. a nearest) vertex of the heap */
793  const int k4 = pathheapGetNearest(heap, state, &nHeapElems, path);
794  const int k = k4 - tnnodes;
795  assert(0 <= k && k < nnodes);
796
797  /* mark vertex k as removed from heap */
798  state[k4] = UNKNOWN;
799
800  /* iterate over all outgoing edges of vertex k */
801  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
802  {
804
807  {
809  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, k4, head4, e, sdprofit, path);
810
811  /* check whether the path4 (to j) including k is shorter than the so far best known */
813  {
814  if( extend && !nodes_isvisited[head] )
815  {
818  }
819
820  pathheapCorrectDist(heap, state, &nHeapElems, path, head4, k4, e, distnew);
822  }
823  }
824  }
825  }
826 }
827
828
829 /** initializes */
830 static inline
832  SCIP* scip, /**< SCIP */
833  int level, /**< 0 - 3 */
834  const GRAPH* g, /**< graph */
835  TREPAIR* repair /**< data for repairing */
836  )
837 {
838  const int nnodes = graph_get_nNodes(g);
839
840 #ifndef NDEBUG
841  const int* const state = repair->tpaths->state;
842  for( int i = 0; i < nnodes; ++i )
843  {
844  assert(state[i] == UNKNOWN);
845  }
846 #endif
847  assert(0 <= level && level <= 3);
848  assert(!repair->nodes_isvisited);
849  SCIP_CALL( SCIPallocCleanBufferArray(scip, &(repair->nodes_isvisited), nnodes) );
850
851  repair->nHeapElems = -1;
852
853  return SCIP_OKAY;
854 }
855
856
857 /** cleans and frees */
858 static inline
860  SCIP* scip, /**< SCIP */
861  int level, /**< 0 - 3 */
862  const GRAPH* g, /**< graph */
863  TREPAIR* repair /**< data for repairing */
864  )
865 {
866  const STP_Vectype(int) resetnodes_level = repair->resetnodes[level];
867  int* RESTRICT state = repair->tpaths->state;
868  SCIP_Bool* RESTRICT nodes_isvisited = repair->nodes_isvisited;
869  const int size = StpVecGetSize(resetnodes_level);
870  const int shift = g->knots * level;
871
872  assert(0 <= level && level <= 3);
873
874  /* clear visited list */
875  for( int i = 0; i < size; i++ )
876  {
877  const int node = resetnodes_level[i];
878  assert(graph_knot_isInRange(g, node));
879
880  nodes_isvisited[node] = FALSE;
881  state[node + shift] = UNKNOWN;
882  }
883
884 #ifndef NDBEUG
885  for( int i = 0; i < g->knots; ++i )
886  {
887  assert(!nodes_isvisited[i]);
888  }
889
890  for( int i = 0; i < g->knots * STP_TPATHS_NTERMBASES; ++i )
891  {
892  assert(UNKNOWN == state[i]);
893  }
894 #endif
895
896  SCIPfreeCleanBufferArray(scip, &(repair->nodes_isvisited));
897 }
898
899
900 /** traverses graph to find nodes that need to be reseted */
901 static
903  SCIP* scip, /**< SCIP */
904  int start, /**< node to start from */
905  int pred, /**< predecessor node */
906  const GRAPH* g, /**< graph */
907  TREPAIR* repair /**< data for repairing */
908 )
909 {
910  const TPATHS* tpaths = repair->tpaths;
911  const PATH* termpaths = tpaths->termpaths;
912  SCIP_Bool* RESTRICT nodes_isvisited = repair->nodes_isvisited;
913
914  if( termpaths[start].edge >= 0 && g->tail[termpaths[start].edge] == pred )
915  {
916  STP_Vectype(int) resetnodes1st = repair->resetnodes[0];
917  STP_Vectype(int) stack = repair->stack;
918
919  assert(!Is_term(g->term[start]));
920  assert(resetnodes1st);
921
922  tpathsRepairResetNode(scip, start, resetnodes1st, stack, nodes_isvisited);
923
924  SCIPdebugMessage("starting DFS from %d ... \n", start);
925
926  /* DFS loop */
927  while( !StpVecIsEmpty(stack) )
928  {
929  const int node = stack[StpVecGetSize(stack) - 1];
930  StpVecPopBack(stack);
931
932  assert(graph_knot_isInRange(g, node));
933
934  for( int a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
935  {
937
939  continue;
940
942  {
945
947
948  tpathsRepairResetNode(scip, head, resetnodes1st, stack, nodes_isvisited);
949  }
950  }
951  }
952
953  repair->stack = stack;
954  repair->resetnodes[0] = resetnodes1st;
955  }
956
957  return SCIP_OKAY;
958 }
959
960
961 /** DFS on level */
962 static
964  SCIP* scip, /**< SCIP */
965  const GRAPH* g, /**< graph */
966  int level, /**< level from 1-3 */
967  TREPAIR* repair /**< data for repairing */
968 )
969 {
970  TPATHS* tpaths = repair->tpaths;
971  const PATH* termpaths = tpaths->termpaths;
972  const int* vbase_org = tpaths->termbases_org;
973  SCIP_Bool* RESTRICT nodes_isvisited = repair->nodes_isvisited;
974  const int nnodes = graph_get_nNodes(g);
975  const int shift = level * nnodes;
976  STP_Vectype(int) stack = repair->stack;
977  STP_Vectype(int) resetnodes_level = repair->resetnodes[level];
978
979  assert(scip && nodes_isvisited && resetnodes_level);
980  assert(1 <= level && level <= 3);
981
982  SCIPdebugMessage("starting DFS for level %d ... \n", level);
983
984  /* DFS loop */
985  while( !StpVecIsEmpty(stack) )
986  {
987  const int node = stack[StpVecGetSize(stack) - 1];
988  const int nodebaseorg_top = vbase_org[node + shift];
989  StpVecPopBack(stack);
990
991  assert(graph_knot_isInRange(g, node));
992
993  for( int a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
994  {
996
998  {
1000
1002  {
1003  if( nodebaseorg_top != vbase_org[head_shifted] )
1004  continue;
1005
1007
1008  tpathsRepairResetNode(scip, head, resetnodes_level, stack, nodes_isvisited);
1009  }
1010  }
1011  }
1012  }
1013
1014  repair->stack = stack;
1015  repair->resetnodes[level] = resetnodes_level;
1016 }
1017
1018
1019 /** adds outdated nodes whose predecessor points along edge to be deleted */
1020 static
1022  SCIP* scip, /**< SCIP */
1023  const GRAPH* g, /**< graph */
1024  int level, /**< level (1-3, level 0 is handled separately) */
1025  TREPAIR* repair /**< data for repairing */
1026 )
1027 {
1028  const TPATHS* tpaths = repair->tpaths;
1029  const PATH* termpaths = tpaths->termpaths;
1030  const int nnodes = graph_get_nNodes(g);
1031  const int shift = level * nnodes;
1032  const int edge = repair->edge;
1033  const int tail = g->tail[edge];
1035
1036  if( termpaths[tail + shift].edge >= 0 && g->tail[termpaths[tail + shift].edge] == head )
1037  {
1038  tpathsRepairResetNode(scip, tail, repair->resetnodes[level], repair->stack, repair->nodes_isvisited);
1039  }
1040
1041  if( termpaths[head + shift].edge >= 0 && g->tail[termpaths[head + shift].edge] == tail )
1042  {
1043  tpathsRepairResetNode(scip, head, repair->resetnodes[level], repair->stack, repair->nodes_isvisited);
1044  }
1045 }
1046
1047
1048 /** adds outdated nodes whose predecessor point to lower level */
1049 static
1051  SCIP* scip, /**< SCIP */
1052  const GRAPH* g, /**< graph */
1053  int level, /**< level (1-3, level 0 is handled separately) */
1054  TREPAIR* repair /**< data for repairing */
1055 )
1056 {
1057  const TPATHS* tpaths = repair->tpaths;
1058  const PATH* termpaths = tpaths->termpaths;
1059  const int* vbase = tpaths->termbases;
1060  const int* vbase_org = tpaths->termbases_org;
1061  SCIP_Bool* RESTRICT nodes_isvisited = repair->nodes_isvisited;
1062  const int nnodes = graph_get_nNodes(g);
1063  const int shift = level * nnodes;
1064
1065  assert(repair->resetnodes[level] && vbase && nodes_isvisited);
1066  assert(1 <= level && level <= 3);
1067
1068  for( int d = 0; d < level; d++ )
1069  {
1070  const STP_Vectype(int) resetnodes_d = repair->resetnodes[d];
1071  const int size = StpVecGetSize(resetnodes_d);
1072
1073  assert(size >= 0);
1074
1075  for( int i = 0; i < size; i++ )
1076  {
1077  const int node = resetnodes_d[i];
1078  const int nodebase_d = vbase[node + d * nnodes];
1079  const int nodebaseorg_d = vbase_org[node + d * nnodes];
1080  assert(graph_knot_isInRange(g, node));
1081
1082  /* duplicate? */
1083  if( !nodes_isvisited[node] && nodebase_d == vbase[node + shift] )
1084  {
1085  tpathsRepairResetNode(scip, node, repair->resetnodes[level], repair->stack, nodes_isvisited);
1086  }
1087
1088  for( int a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
1089  {
1091
1093  {
1095
1097  {
1098  if( nodebaseorg_d != vbase[head_shifted] )
1099  continue;
1100
1102
1104
1105  tpathsRepairResetNode(scip, head, repair->resetnodes[level], repair->stack, nodes_isvisited);
1106  }
1107  }
1108  }
1109  }
1110  }
1111 }
1112
1113
1114 /** Traverses graph to find nodes on given level that need to be reseted,
1115  * by updating from previous levels. Also considers duplicates. */
1116 static
1118  SCIP* scip, /**< SCIP */
1119  const GRAPH* g, /**< graph */
1120  int level, /**< level (1-3, level 0 is handled separately) */
1121  TREPAIR* repair /**< data for repairing */
1122 )
1123 {
1124  assert(1 <= level && level <= 3);
1125
1128
1129  tpathsRepairTraverseLevelWithStack(scip, g, level, repair);
1130 }
1131
1132
1133 /** updates reset nodes from non-reset nodes */
1134 static
1136  SCIP* scip, /**< SCIP */
1137  const GRAPH* g, /**< graph */
1138  TREPAIR* repair /**< data for repairing */
1139 )
1140 {
1141  const STP_Vectype(int) resetnodes1st = repair->resetnodes[0];
1142  const SCIP_Bool* nodes_isvisited = repair->nodes_isvisited;
1143  TPATHS* tpaths = repair->tpaths;
1144  PATH* RESTRICT termpaths = tpaths->termpaths;
1145  int* RESTRICT termbases = tpaths->termbases;
1146  int* RESTRICT heap = g->path_heap;
1147  int* RESTRICT state = tpaths->state;
1148  int nHeapElems = 0;
1149  const int edge = repair->edge;
1150  const int edge_rev = flipedge(edge);
1151  const int size = StpVecGetSize(resetnodes1st);
1152
1153  for( int i = 0; i < size; i++ )
1154  {
1155  const int node = resetnodes1st[i];
1156  assert(graph_knot_isInRange(g, node));
1157  assert(!Is_term(g->term[node]));
1158
1159  pathResetSingle(termpaths[node]);
1160  termbases[node] = UNKNOWN;
1161
1162  for( int e = g->inpbeg[node]; e != EAT_LAST; e = g->ieat[e] )
1163  {
1164  SCIP_Real dist;
1165  const int tail = g->tail[e];
1166
1167  if( nodes_isvisited[tail] || e == edge || e == edge_rev )
1168  continue;
1169
1170  dist = g->cost[e] + termpaths[tail].dist;
1171
1172  if( dist < termpaths[node].dist )
1173  {
1174  termpaths[node].dist = dist;
1175  termpaths[node].edge = e;
1176  termbases[node] = termbases[tail];
1177
1178  assert(graph_knot_isInRange(g, termbases[node]));
1179  }
1180  }
1181
1182  if( termbases[node] != UNKNOWN )
1183  {
1184  SCIPdebugMessage("update node %d: base=%d, pred=%d dist=%f \n", node,
1185  termbases[node], g->tail[termpaths[node].edge], termpaths[node].dist);
1186
1187  graph_pathHeapAdd(termpaths, node, heap, state, &nHeapElems);
1188  }
1189  }
1190
1191  repair->nHeapElems = nHeapElems;
1192 }
1193
1194
1195
1196
1197 /** updates reset nodes from non-reset nodes of given level */
1198 static
1200  SCIP* scip, /**< SCIP */
1201  const GRAPH* g, /**< graph */
1202  int level, /**< level (1-3, 0 is handled separately) */
1203  TREPAIR* repair /**< data for repairing */
1204 )
1205 {
1206  const STP_Vectype(int) resetnodes_top = repair->resetnodes[level];
1207  const SCIP_Bool* nodes_isvisited = repair->nodes_isvisited;
1208  TPATHS* tpaths = repair->tpaths;
1209  PATH* RESTRICT termpaths = tpaths->termpaths;
1210  int* RESTRICT vbase = tpaths->termbases;
1211  int* RESTRICT heap = g->path_heap;
1212  int* RESTRICT state = tpaths->state;
1213  int nHeapElems = 0;
1214  const int nnodes = graph_get_nNodes(g);
1215  const int shift = level * nnodes;
1216  const int edge = repair->edge;
1217  const int edge_rev = flipedge(edge);
1218  const int size = StpVecGetSize(resetnodes_top);
1219
1220  assert(1 <= level && level <= 3);
1221
1222  for( int i = 0; i < size; i++ )
1223  {
1224  const int node = resetnodes_top[i];
1225  const int node_shift = node + shift;
1226  assert(graph_knot_isInRange(g, node));
1227  assert(!Is_term(g->term[node]));
1228  assert(nodes_isvisited[node]);
1229
1230  pathResetSingle(termpaths[node_shift]);
1231  vbase[node_shift] = UNKNOWN;
1232
1233  for( int e = g->inpbeg[node]; e != EAT_LAST; e = g->ieat[e] )
1234  {
1235  const int tail = g->tail[e];
1236
1237  if( e == edge || e == edge_rev )
1238  continue;
1239
1240  /* loop over all levels of tail */
1241  for( int d = 0; d <= level; d++ )
1242  {
1243  SCIP_Real dist;
1244  const int tail_d = tail + d * nnodes;
1245  const int tailbase_d = vbase[tail_d];
1246  int h;
1247
1248  if( d == level && nodes_isvisited[tail] )
1249  break;
1250
1251  for( h = 0; h < level; h++ )
1252  if( vbase[node + h * nnodes] == tailbase_d )
1253  break;
1254
1255  /* has tail_d a base from node? */
1256  if( h != level )
1257  continue;
1258
1259  dist = g->cost[e] + termpaths[tail_d].dist;
1260
1261  if( dist < termpaths[node_shift].dist )
1262  {
1263  termpaths[node_shift].dist = dist;
1264  termpaths[node_shift].edge = e;
1265  vbase[node_shift] = tailbase_d;
1266
1267  assert(graph_knot_isInRange(g, vbase[node_shift]));
1268  }
1269  }
1270  }
1271
1272  if( vbase[node_shift] != UNKNOWN )
1273  {
1274  SCIPdebugMessage("update node %d: base=%d, pred=%d dist=%f \n", node,
1275  vbase[node_shift], g->tail[termpaths[node_shift].edge], termpaths[node_shift].dist);
1276
1277  graph_pathHeapAdd(termpaths, node_shift, heap, state, &nHeapElems);
1278  }
1279  }
1280
1281  repair->nHeapElems = nHeapElems;
1282 }
1283
1284
1285 /** repairs TPATHS structure for imminent edge deletion (1st level) */
1286 static
1288  SCIP* scip, /**< SCIP */
1289  const GRAPH* g, /**< graph */
1290  TREPAIR* repair /**< data for repairing */
1291 )
1292 {
1293  const int edge = repair->edge;
1294  const int tail = g->tail[edge];
1296
1297  SCIP_CALL( tpathsRepairInitLevel(scip, 0, g, repair) );
1298
1299  /* find nodes that need to be reseted */
1300  SCIP_CALL( tpathsRepairTraverse1st(scip, head, tail, g, repair) );
1301  SCIP_CALL( tpathsRepairTraverse1st(scip, tail, head, g, repair) );
1302
1303  /* update newly found nodes from remaining nodes */
1304  tpathsRepairUpdate1st(scip, g, repair);
1305
1306  /* complete the repair process for this level */
1307  tpathsScan1st(scip, g, g->cost, NULL, repair->nHeapElems, repair->tpaths,
1308  repair->withEdgeReplacement ? repair : NULL);
1309
1310  tpathsRepairExitLevel(scip, 0, g, repair);
1311
1312  return SCIP_OKAY;
1313 }
1314
1315
1316 /** repairs TPATHS structure for imminent edge deletion (2nd level) */
1317 static
1319  SCIP* scip, /**< SCIP */
1320  const GRAPH* g, /**< graph */
1321  TREPAIR* repair /**< data for repairing */
1322 )
1323 {
1324  const int level = 1; /* NOTE: level is 0-indexed */
1325
1326  SCIP_CALL( tpathsRepairInitLevel(scip, level, g, repair) );
1327
1328  tpathsRepairTraverseLevel(scip, g, level, repair);
1329  tpathsRepairUpdateLevel(scip, g, level, repair);
1330  tpathsScan2nd(scip, g, g->cost, g->cost, NULL, repair->nHeapElems, repair->tpaths,
1331  repair->withEdgeReplacement ? repair : NULL);
1332
1333  tpathsRepairExitLevel(scip, level, g, repair);
1334
1335  return SCIP_OKAY;
1336 }
1337
1338 /** repairs TPATHS structure for imminent edge deletion (3rd level) */
1339 static
1341  SCIP* scip, /**< SCIP */
1342  const GRAPH* g, /**< graph */
1343  TREPAIR* repair /**< data for repairing */
1344 )
1345 {
1346  const int level = 2; /* NOTE: level is 0-indexed */
1347
1348  SCIP_CALL( tpathsRepairInitLevel(scip, level, g, repair) );
1349
1350  tpathsRepairTraverseLevel(scip, g, level, repair);
1351  tpathsRepairUpdateLevel(scip, g, level, repair);
1352  tpathsScan3rd(scip, g, g->cost, g->cost, NULL, repair->nHeapElems, repair->tpaths,
1353  repair->withEdgeReplacement ? repair : NULL);
1354
1355  tpathsRepairExitLevel(scip, level, g, repair);
1356
1357  return SCIP_OKAY;
1358 }
1359
1360
1361 /** repairs TPATHS structure for imminent edge deletion (4th level) */
1362 static
1364  SCIP* scip, /**< SCIP */
1365  const GRAPH* g, /**< graph */
1366  TREPAIR* repair /**< data for repairing */
1367 )
1368 {
1369  const int level = 3; /* NOTE: level is 0-indexed */
1370
1371  SCIP_CALL( tpathsRepairInitLevel(scip, level, g, repair) );
1372
1373  tpathsRepairTraverseLevel(scip, g, level, repair);
1374  tpathsRepairUpdateLevel(scip, g, level, repair);
1375  tpathsScan4th(scip, g, g->cost, g->cost, NULL, repair->nHeapElems, repair->tpaths,
1376  repair->withEdgeReplacement ? repair : NULL);
1377
1378  tpathsRepairExitLevel(scip, level, g, repair);
1379
1380  return SCIP_OKAY;
1381 }
1382
1383
1384 /** initializes */
1385 static
1387  SCIP* scip, /**< SCIP */
1388  const GRAPH* g, /**< graph */
1389  TREPAIR* repair /**< data for repairing */
1390 )
1391 {
1392
1393 #ifndef NDEBUG
1394  TPATHS* tpaths = repair->tpaths;
1395  for( int i = 0; i < STP_TPATHS_NTERMBASES * g->knots; i++ )
1396  {
1397  assert(tpaths->termbases_org[i] == tpaths->termbases[i]);
1398  }
1399 #endif
1400
1401  StpVecReserve(scip, repair->stack, STP_TPATHS_RESERVESIZE);
1402
1403  for( int i = 0; i < STP_TPATHS_NTERMBASES; i++ )
1404  {
1405  assert(!repair->resetnodes[i]);
1406  StpVecReserve(scip, repair->resetnodes[i], STP_TPATHS_RESERVESIZE);
1407  assert(repair->resetnodes[i]);
1408  }
1409 }
1410
1411
1412 /** frees and resets */
1413 static
1415  SCIP* scip, /**< SCIP */
1416  const GRAPH* g, /**< graph */
1417  TREPAIR* repair /**< data for repairing */
1418 )
1419 {
1420  TPATHS* tpaths = repair->tpaths;
1421  int* RESTRICT vbase_org = tpaths->termbases_org;
1422  const int* const vbase = tpaths->termbases;
1423  const int nnodes = graph_get_nNodes(g);
1424
1425  for( int i = STP_TPATHS_NTERMBASES - 1; i >= 0; i-- )
1426  {
1427  STP_Vectype(int) resetnodes_i = repair->resetnodes[i];
1428  const int size = StpVecGetSize(resetnodes_i);
1429  const int shift = nnodes * i;
1430
1431  assert(resetnodes_i);
1432
1433  for( int j = 0; j < size; j++ )
1434  {
1435  const int node = resetnodes_i[j];
1436  const int node_shift = node + shift;
1437  assert(graph_knot_isInRange(g, node));
1438
1439  vbase_org[node_shift] = vbase[node_shift];
1440  }
1441
1442  StpVecFree(scip, repair->resetnodes[i]);
1443  }
1444
1445  StpVecFree(scip, repair->stack);
1446
1447 #ifndef NDEBUG
1448  for( int i = 0; i < STP_TPATHS_NTERMBASES * nnodes; i++ )
1449  {
1450  assert(vbase_org[i] == vbase[i]);
1451
1452  }
1453 #endif
1454 }
1455
1456
1457 /*
1458  * Interface methods
1459  */
1460
1461
1462 /** builds a Voronoi region w.r.t. shortest paths, for all terminals
1463  * NOTE: serves as a wrapper */
1465  const GRAPH* g, /**< graph data structure */
1466  const SCIP_Real* cost, /**< edge costs */
1467  PATH* path, /**< path data structure (leading to respective Voronoi base) */
1468  int* vbase, /**< Voronoi base to each vertex */
1469  int* state /**< array to mark the state of each node during calculation */
1470  )
1471 {
1472  const int nnodes = graph_get_nNodes(g);
1473  TPATHS tpaths = { .termpaths = path, .termbases = vbase, .state = state, .nnodes = nnodes };
1474
1475  assert(path && vbase && state);
1476
1478 }
1479
1480 /** 2nd next terminal to all non terminal nodes
1481  * NOTE: legacy wrapper */
1483  const GRAPH* g, /**< graph data structure */
1484  const SCIP_Real* cost, /**< edge costs */
1485  const SCIP_Real* costrev, /**< reversed edge costs */
1486  PATH* path2, /**< path data structure (leading to first and second nearest terminal) */
1487  int* vbase2, /**< first and second nearest terminal to each non terminal */
1488  int* state2 /**< array to mark the state of each node during calculation */
1489  )
1490 {
1491  const int nnodes = graph_get_nNodes(g);
1492  TPATHS tpaths = { .termpaths = path2, .termbases = vbase2, .state = state2, .nnodes = nnodes };
1493
1494  assert(path2 && vbase2 && state2);
1495
1496  graph_tpathsAdd2nd(g, cost, costrev, NULL, &tpaths);
1497 }
1498
1499 /** 3rd next terminal to all non terminal nodes
1500  * NOTE: legacy wrapper */
1502  const GRAPH* g, /**< graph data structure */
1503  const SCIP_Real* cost, /**< edge costs */
1504  const SCIP_Real* costrev, /**< reversed edge costs */
1505  PATH* path3, /**< path data structure (leading to first, second and third nearest terminal) */
1506  int* vbase3, /**< first, second and third nearest terminal to each non terminal */
1507  int* state3 /**< array to mark the state of each node during calculation */
1508  )
1509 {
1510  const int nnodes = graph_get_nNodes(g);
1511  TPATHS tpaths = { .termpaths = path3, .termbases = vbase3, .state = state3, .nnodes = nnodes };
1512
1513  assert(path3 && vbase3 && state3);
1514
1515  graph_tpathsAdd3rd(g, cost, costrev, NULL, &tpaths);
1516 }
1517
1518
1519 /* 4th next terminal to all non terminal nodes
1520  * NOTE: legacy wrapper */
1522  const GRAPH* g, /**< graph data structure */
1523  const SCIP_Real* cost, /**< edge costs */
1524  const SCIP_Real* costrev, /**< reversed edge costs */
1525  PATH* path4, /**< path data structure (leading to first, second, third and fourth nearest terminal) */
1526  int* vbase4, /**< first, second, third, and fourth nearest terminal to each non terminal */
1527  int* state4 /**< array to mark the state of each node during calculation */
1528  )
1529 {
1530  const int nnodes = graph_get_nNodes(g);
1531  TPATHS tpaths = { .termpaths = path4, .termbases = vbase4, .state = state4, .nnodes = nnodes };
1532
1533  assert(path4 && vbase4 && state4);
1534
1535  graph_tpathsAdd4th(g, cost, costrev, NULL, &tpaths);
1536 }
1537
1538
1539
1540 /** gets non-terminal shortest paths to 2 closest terminal for each non-terminal
1541  * NOTE: legacy wrapper */
1543  GRAPH* g, /**< graph data structure */
1544  const SCIP_Real* cost, /**< edge costs */
1545  const SCIP_Real* costrev, /**< reversed edge costs */
1546  PATH* path2, /**< path data structure (leading to first, second terminal) */
1547  int* vbase2, /**< first, second and third nearest terminal to each non terminal */
1548  int* state2 /**< array to mark the state of each node during calculation */
1549  )
1550 {
1551  const int nnodes = graph_get_nNodes(g);
1552  TPATHS tpaths = { .termpaths = path2, .termbases = vbase2, .state = state2, .nnodes = nnodes };
1553
1554  assert(path2 && vbase2 && state2);
1555
1556  graph_tpathsSetAll2(g, cost, costrev, NULL, &tpaths);
1557 }
1558
1559
1560 /** gets non-terminal shortest paths to 3 closest terminal for each non-terminal
1561  * NOTE: legacy wrapper */
1563  GRAPH* g, /**< graph data structure */
1564  const SCIP_Real* cost, /**< edge costs */
1565  const SCIP_Real* costrev, /**< reversed edge costs */
1566  PATH* path3, /**< path data structure (leading to first, second and third nearest terminal) */
1567  int* vbase3, /**< first, second and third nearest terminal to each non terminal */
1568  int* state3 /**< array to mark the state of each node during calculation */
1569  )
1570 {
1571  const int nnodes = graph_get_nNodes(g);
1572  TPATHS tpaths = { .termpaths = path3, .termbases = vbase3, .state = state3, .nnodes = nnodes };
1573
1574  assert(path3 && vbase3 && state3);
1575
1576  graph_tpathsSetAll3(g, cost, costrev, NULL, &tpaths);
1577 }
1578
1579
1580 /** gets non-terminal shortest paths to 4 closest terminal for each non-terminal
1581  * NOTE: legacy wrapper */
1583  GRAPH* g, /**< graph data structure */
1584  const SCIP_Real* cost, /**< edge costs */
1585  const SCIP_Real* costrev, /**< reversed edge costs */
1586  PATH* path4, /**< path data struture (leading to first, second, third and fouth nearest terminal) */
1587  int* vbase4, /**< first, second and third nearest terminal to each non terminal */
1588  int* state4 /**< array to mark the state of each node during calculation */
1589  )
1590 {
1591  const int nnodes = graph_get_nNodes(g);
1592  TPATHS tpaths = { .termpaths = path4, .termbases = vbase4, .state = state4, .nnodes = nnodes };
1593
1594  assert(path4 && vbase4 && state4);
1595
1596  graph_tpathsSetAll4(g, cost, costrev, NULL, &tpaths);
1597 }
1598
1599
1600 /** get 4 close terminals to each terminal */
1602  SCIP* scip, /**< SCIP data structure */
1603  GRAPH* g, /**< graph data structure */
1604  const SCIP_Real* cost, /**< edge costs */
1605  PATH* path, /**< path data structure (leading to first, second, third and fourth nearest terminal) */
1606  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1607  int* heap, /**< array representing a heap */
1608  int* state /**< array to mark the state of each node during calculation */
1609  )
1610 {
1611  int* boundedges;
1612  int k;
1613  int e;
1614  int l;
1615  int k2;
1616  int bedge;
1617  int shift;
1618  int nnodes;
1619  int nboundedges;
1620
1621  assert(g != NULL);
1622  assert(path != NULL);
1623  assert(cost != NULL);
1624  assert(heap != NULL);
1625  assert(state != NULL);
1626
1627  SCIP_CALL( SCIPallocBufferArray(scip, &boundedges, g->edges) );
1628
1629  shift = 0;
1630  nnodes = g->knots;
1631
1632  if( !graph_pc_isPcMw(g) )
1633  graph_mark(g);
1634
1635  nboundedges = 0;
1636  for( k = 0; k < nnodes; k++ )
1637  {
1638  if( !g->mark[k] )
1639  continue;
1640
1641  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1642  {
1644  if( !g->mark[k2] || k2 < k )
1645  continue;
1646
1647  /* is e a boundary edge? */
1648  if( vbase[k] != vbase[k2] )
1649  {
1650  boundedges[nboundedges++] = e;
1651  }
1652  }
1653  if( Is_term(g->term[k]) )
1654  {
1655  path[k].dist = FARAWAY;
1656  vbase[k] = UNKNOWN;
1657  }
1658
1659  }
1660
1661  for( l = 0; l < 4; l++ )
1662  {
1663  for( e = 0; e < nboundedges; e++ )
1664  {
1665  bedge = boundedges[e];
1666  k = g->tail[bedge];
1668  utdist(scip, g, path, cost[bedge], vbase, k, l, k2, shift, nnodes);
1669  utdist(scip, g, path, cost[bedge], vbase, k2, l, k, shift, nnodes);
1670  }
1671  shift += nnodes;
1672  }
1673
1674  SCIPfreeBufferArray(scip, &boundedges);
1675
1676  return SCIP_OKAY;
1677 }
1678
1679
1680
1681 /** initializes TPATHS structure */
1683  SCIP* scip, /**< SCIP */
1684  GRAPH* g, /**< graph NOTE: will mark the graph, thus not const :(
1685  terrible design */
1686  TPATHS** tpaths /**< the terminal paths */
1687 )
1688 {
1689  assert(scip && g);
1690  assert(STP_TPATHS_NTERMBASES == 4);
1691
1692  SCIP_CALL( tpathsAlloc(scip, g, tpaths) );
1693  tpathsBuild(scip, g, *tpaths);
1694
1695  return SCIP_OKAY;
1696 }
1697
1698
1699 /** initializes biased TPATHS structure */
1701  SCIP* scip, /**< SCIP */
1702  const SDPROFIT* sdprofit, /**< SD profit */
1703  GRAPH* g, /**< graph NOTE: will mark the graph, thus not const :(
1704  terrible design */
1705  TPATHS** tpaths /**< the terminal paths */
1706 )
1707 {
1708  assert(scip && g && sdprofit);
1709  assert(STP_TPATHS_NTERMBASES == 4);
1710
1711  SCIP_CALL( tpathsAlloc(scip, g, tpaths) );
1712  tpathsBuildBiased(sdprofit, g, *tpaths);
1713
1714  return SCIP_OKAY;
1715 }
1716
1717
1718 /** sets up TPATHS for subsequent edge deletions */
1720  const GRAPH* g, /**< graph */
1721  TPATHS* tpaths /**< the terminal paths */
1722 )
1723 {
1724  const int nnodes = graph_get_nNodes(g);
1725  int* RESTRICT state = tpaths->state;
1726
1727  if( !tpaths->termbases_org )
1728  {
1730  }
1731
1733
1734  for( int i = 0; i < nnodes * STP_TPATHS_NTERMBASES; i++ )
1735  {
1736  state[i] = UNKNOWN;
1737  }
1738
1739  return SCIP_OKAY;
1740 }
1741
1742
1743 /** repairs TPATHS structure for imminent edge deletion */
1745  SCIP* scip, /**< SCIP */
1746  int edge, /**< edge about to be eliminated */
1747  SCIP_Bool withEdgeReplacement,/**< with edge replacement? */
1748  const GRAPH* g, /**< graph */
1749  TPATHS* tpaths /**< the terminal paths */
1750 )
1751 {
1752  STP_Vectype(int) resetnodes[STP_TPATHS_NTERMBASES] = { NULL, NULL, NULL, NULL };
1753  TREPAIR repair = { .resetnodes = resetnodes, .stack = NULL, .tpaths = tpaths, .nodes_isvisited = NULL,
1754  .edge = edge, .nHeapElems = 0, .withEdgeReplacement = withEdgeReplacement };
1755
1756  assert(scip && g && tpaths);
1757  assert(STP_TPATHS_NTERMBASES == 4);
1758  assert(graph_edge_isInRange(g, edge));
1759  assert(graph_isMarked(g));
1760  assert(tpaths->termbases_org);
1761
1762  tpathsRepairInit(scip, g, &repair);
1763
1764  SCIP_CALL( tpathsRepair1st(scip, g, &repair) );
1765  SCIP_CALL( tpathsRepair2nd(scip, g, &repair) );
1766  SCIP_CALL( tpathsRepair3rd(scip, g, &repair) );
1767  SCIP_CALL( tpathsRepair4th(scip, g, &repair) );
1768
1769  tpathsRepairExit(scip, g, &repair);
1770
1771  return SCIP_OKAY;
1772 }
1773
1774
1775 /** recomputes biased TPATHS structure */
1777  const SDPROFIT* sdprofit, /**< SD profit */
1778  GRAPH* g, /**< graph NOTE: will mark the graph, thus not const
1779  :( terrible design */
1780  TPATHS* tpaths /**< the terminal paths */
1781 )
1782 {
1783  assert(tpaths && g && sdprofit);
1784  assert(STP_TPATHS_NTERMBASES == 4);
1785
1786  tpathsBuildBiased(sdprofit, g, tpaths);
1787
1788  return SCIP_OKAY;
1789 }
1790
1791
1792 /** prints terminal paths from given node */
1794  const GRAPH* g, /**< graph data structure */
1795  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
1796  const TPATHS* tpaths, /**< storage for terminal paths */
1797  int node /**< node to start from */
1798 )
1799 {
1800  const int* termbases;
1801  const int nnodes = graph_get_nNodes(g);
1802  int pos;
1803
1804  assert(0 && "currently does not work, need to use vbase to find correct ancestor");
1805  assert(tpaths);
1806  assert(graph_knot_isInRange(g, node));
1807
1808  printf("printing terminal paths for ");
1809  graph_knot_printInfo(g, node);
1810
1811  if( Is_term(g->term[node]))
1812  {
1813  printf("single node path only: \n");
1814  graph_knot_printInfo(g, node);
1815
1816  return;
1817  }
1818
1819  termbases = tpaths->termbases;
1820  pos = node;
1821
1822  for( int k = 0; k < STP_TPATHS_NTERMBASES; k++ )
1823  {
1824  const int base = termbases[pos];
1825
1826  if( base != UNKNOWN )
1827  {
1828  assert(graph_knot_isInRange(g, base));
1829
1830  printf("Path to terminal (distance=%f) ", tpaths->termpaths[pos].dist);
1831  graph_knot_printInfo(g, base);
1832
1833  tpathsPrintPath(g, sdprofit, tpaths, node, k + 1);
1834  }
1835  pos += nnodes;
1836  }
1837
1838 }
1839
1840
1841 /** gets nodes with positive profit on path */
1843  SCIP* scip, /**< SCIP */
1844  const GRAPH* g, /**< graph data structure */
1845  const TPATHS* tpaths, /**< storage for terminal paths */
1846  const SDPROFIT* sdprofit, /**< SD bias for nodes */
1847  int start, /**< start vertex */
1848  int base, /**< base vertex (terminal) */
1849  STP_Vectype(int) profitnodes /**< NOTE: needs to be of capacity at least g->knots! */
1850 )
1851 {
1852  const int nnodes = graph_get_nNodes(g);
1853  int level;
1854
1855  assert(tpaths && profitnodes && sdprofit);
1856  assert(graph_knot_isInRange(g, start));
1857  assert(graph_knot_isInRange(g, base));
1858  assert(Is_term(g->term[base]));
1859  assert(StpVecGetcapacity(profitnodes) >= g->knots);
1860
1861  StpVecClear(profitnodes);
1862
1863  if( base != start )
1864  {
1865  const int* const termbases = tpaths->termbases;
1866  const PATH* const termpaths = tpaths->termpaths;
1867  int node_pred;
1868  int node_curr;
1869  int node_next;
1870  int edge_curr;
1871
1872  for( level = 0; termbases[start + nnodes * level] != base; level++ )
1873  {
1874  assert(level <= STP_TPATHS_NTERMBASES);
1875  }
1876  assert(termbases[start + nnodes * level] == base);
1877  assert(0 <= level && level <= STP_TPATHS_NTERMBASES);
1878
1879  node_pred = -1;
1880  node_curr = start;
1881  edge_curr = termpaths[node_curr + level * nnodes].edge;
1882  assert(edge_curr >= 0);
1883  node_next = g->tail[edge_curr];
1884
1885  // printf("%d->%d \n", start, base);
1886
1887  while( node_curr != base )
1888  {
1889  assert(graph_knot_isInRange(g, node_curr));
1890  assert(graph_edge_isInRange(g, edge_curr));
1891  assert(termbases[node_curr + level * nnodes] == base);
1892
1893  // printf("%d \n", node_curr);
1894
1895  if( node_curr != start )
1896  {
1897  const SCIP_Real profit = reduce_sdprofitGetProfit(sdprofit, node_curr, node_pred, node_next);
1898  if( GT(profit, 0.0) )
1899  {
1900  StpVecPushBack(scip, profitnodes, node_curr);
1901  }
1902  }
1903
1904  node_pred = node_curr;
1905  node_curr = node_next;
1906
1907  while( termbases[node_curr + level * nnodes] != base )
1908  {
1909  level--;
1910  assert(level >= 0);
1911  }
1912
1913  edge_curr = termpaths[node_curr + level * nnodes].edge;
1914  node_next = (edge_curr >= 0) ? g->tail[edge_curr] : -1;
1915  }
1916  }
1917
1918 #ifndef NDEBUG
1919  for( int k = 0; k < StpVecGetSize(profitnodes); k++ )
1920  {
1921  const int node = profitnodes[k];
1922  assert(graph_knot_isInRange(g, node));
1923  assert(!Is_term(g->term[node]));
1924  }
1925 #endif
1926 }
1927
1928
1929
1930 /** Computes next terminal to all non-terminal nodes.
1931  * Essentially a Voronoi diagram */
1933  const GRAPH* g, /**< graph data structure */
1934  const SCIP_Real* cost, /**< edge costs */
1935  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
1936  TPATHS* tpaths /**< storage for terminal paths */
1937 )
1938 {
1939  int nHeapElems = 0;
1940  int nbases = 0;
1941  const int nnodes = graph_get_nNodes(g);
1942  int* RESTRICT heap = g->path_heap;
1943  PATH* RESTRICT path = tpaths->termpaths;
1944  int* RESTRICT vbase = tpaths->termbases;
1945  int* RESTRICT state = tpaths->state;
1946
1947  assert(path != NULL);
1948  assert(cost != NULL);
1949  assert(heap != NULL);
1950  assert(state != NULL);
1951  assert(nnodes == tpaths->nnodes);
1952
1953  /* initialize */
1954  for( int i = 0; i < nnodes; i++ )
1955  {
1956  /* set the base of vertex i */
1957  if( Is_term(g->term[i]) && g->mark[i] )
1958  {
1959  nbases++;
1960  if( nnodes > 1 )
1961  {
1962  heap[++nHeapElems] = i;
1963  }
1964
1965  vbase[i] = i;
1966  path[i].dist = 0.0;
1967  path[i].edge = UNKNOWN;
1968  state[i] = nHeapElems;
1969  }
1970  else
1971  {
1972  vbase[i] = UNKNOWN;
1973  path[i].dist = FARAWAY;
1974  path[i].edge = UNKNOWN;
1975  state[i] = UNKNOWN;
1976  }
1977  }
1978
1979  if( nbases == 0 || nnodes <= 1 )
1980  return;
1981
1982  tpathsScan1st(NULL, g, cost, sdprofit, nHeapElems, tpaths, NULL);
1983 }
1984
1985
1986 /** computes 2nd next terminal to all non-terminal nodes */
1988  const GRAPH* g, /**< graph data structure */
1989  const SCIP_Real* cost, /**< edge costs */
1990  const SCIP_Real* costrev, /**< reversed edge costs */
1991  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
1992  TPATHS* tpaths /**< storage for terminal paths */
1993 )
1994 {
1995  int nheapElems;
1996  const int nnodes = graph_get_nNodes(g);
1997  int* RESTRICT heap = g->path_heap;
1998  PATH* RESTRICT path = tpaths->termpaths;
1999  int* RESTRICT vbase = tpaths->termbases;
2000  int* RESTRICT state = tpaths->state;
2001
2002  assert(nnodes == tpaths->nnodes);
2003  assert(path != NULL);
2004  assert(cost != NULL);
2005  assert(heap != NULL);
2006  assert(state != NULL);
2007  assert(costrev != NULL);
2008
2009  nheapElems = 0;
2010  tpathsResetMembers(g, 2, tpaths);
2011
2012  /* scan original nodes */
2013  for( int k = 0; k < nnodes; k++ )
2014  {
2015  if( !g->mark[k] )
2016  continue;
2017
2018  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2019  {
2021
2023  {
2025  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, k, head2, e, sdprofit, path);
2026
2028  {
2029  pathheapCorrectDist(heap, state, &nheapElems, path, head2, k, e, distnew);
2031  }
2032  }
2033  }
2034  }
2035
2036  if( nnodes <= 1 )
2037  return;
2038
2039  tpathsScan2nd(NULL, g, cost, costrev, sdprofit, nheapElems, tpaths, NULL);
2040 }
2041
2042
2043 /** computes 3rd next terminal to all non-terminal nodes */
2045  const GRAPH* g, /**< graph data structure */
2046  const SCIP_Real* cost, /**< edge costs */
2047  const SCIP_Real* costrev, /**< reversed edge costs */
2048  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
2049  TPATHS* tpaths /**< storage for terminal paths */
2050 )
2051 {
2052  int nheapElems = 0;
2053  const int nnodes = graph_get_nNodes(g);
2054  const int dnnodes = 2 * nnodes;
2055  int* RESTRICT heap = g->path_heap;
2056  PATH* RESTRICT path = tpaths->termpaths;
2057  int* RESTRICT vbase = tpaths->termbases;
2058  int* RESTRICT state = tpaths->state;
2059
2060  assert(nnodes == tpaths->nnodes);
2061  assert(path != NULL);
2062  assert(cost != NULL);
2063  assert(heap != NULL);
2064  assert(state != NULL);
2065  assert(costrev != NULL);
2066
2067  tpathsResetMembers(g, 3, tpaths);
2068
2069  /* scan original nodes */
2070  for( int k = 0; k < nnodes; k++ )
2071  {
2072  if( !g->mark[k] )
2073  continue;
2074
2075  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2076  {
2079
2081  {
2082  int kn = k;
2083
2084  for( int level = 0; level < 2; level++ )
2085  {
2086  if( vbase[kn] != vbase[head] && vbase[kn] != vbase[head + nnodes] )
2087  {
2088  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, kn, head3, e, sdprofit, path);
2089
2091  {
2092  pathheapCorrectDist(heap, state, &nheapElems, path, head3, kn, e, distnew);
2094  }
2095  }
2096
2097  kn += nnodes;
2098  }
2099  }
2100  }
2101  }
2102
2103  if( nnodes <= 1 )
2104  return;
2105
2106  tpathsScan3rd(NULL, g, cost, costrev, sdprofit, nheapElems, tpaths, NULL);
2107 }
2108
2109
2110
2111 /** computes 4th next terminal to all non-terminal nodes */
2113  const GRAPH* g, /**< graph data structure */
2114  const SCIP_Real* cost, /**< edge costs */
2115  const SCIP_Real* costrev, /**< reversed edge costs */
2116  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
2117  TPATHS* tpaths /**< storage for terminal paths */
2118 )
2119 {
2120  int nHeapElems = 0;
2121  const int nnodes = graph_get_nNodes(g);
2122  const int dnnodes = 2 * nnodes;
2123  const int tnnodes = 3 * nnodes;
2124  int* RESTRICT heap = g->path_heap;
2125  PATH* RESTRICT path = tpaths->termpaths;
2126  int* RESTRICT vbase = tpaths->termbases;
2127  int* RESTRICT state = tpaths->state;
2128
2129  assert(nnodes == tpaths->nnodes);
2130  assert(path != NULL);
2131  assert(cost != NULL);
2132  assert(costrev != NULL);
2133  assert(heap != NULL);
2134  assert(state != NULL);
2135  assert(costrev != NULL);
2136
2137  tpathsResetMembers(g, 4, tpaths);
2138
2139  /* scan original nodes */
2140  for( int k = 0; k < nnodes; k++ )
2141  {
2142  if( !g->mark[k] )
2143  continue;
2144
2145  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2146  {
2149
2151  {
2152  int kn = k;
2153
2154  for( int level = 0; level < 3; level++ )
2155  {
2156  if( vbase[kn] != vbase[head] && vbase[kn] != vbase[head + nnodes] && vbase[kn] != vbase[head + dnnodes] )
2157  {
2158  const SCIP_Real distnew = tpathsGetDistNew(g, cost, costrev, vbase, kn, head4, e, sdprofit, path);
2159
2161  {
2162  pathheapCorrectDist(heap, state, &nHeapElems, path, head4, kn, e, distnew);
2164  }
2165  }
2166  kn += nnodes;
2167  }
2168  }
2169  }
2170  }
2171
2172  if( nnodes <= 1 )
2173  return;
2174
2175  tpathsScan4th(NULL, g, cost, costrev, sdprofit, nHeapElems, tpaths, NULL);
2176 }
2177
2178
2179 /** computes 2 next terminal to all non-terminal nodes */
2181  GRAPH* g, /**< graph data structure */
2182  const SCIP_Real* cost, /**< edge costs */
2183  const SCIP_Real* costrev, /**< reversed edge costs */
2184  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
2185  TPATHS* tpaths /**< storage for terminal paths */
2186 )
2187 {
2188  assert(g != NULL);
2189  assert(cost != NULL);
2190  assert(costrev != NULL);
2191  assert(tpaths != NULL);
2192
2193  if( !graph_pc_isPcMw(g) )
2194  graph_mark(g);
2195
2197  graph_tpathsAdd2nd(g, cost, costrev, sdprofit, tpaths);
2198
2199 #ifndef NDEBUG
2200  if( !sdprofit )
2201  {
2202  const PATH* RESTRICT path2 = tpaths->termpaths;
2203  const int nnodes = graph_get_nNodes(g);
2204
2205  for( int level = 0; level < 1; level++ )
2206  {
2207  for( int k = 0; k < nnodes; ++k )
2208  {
2209  assert(LE_FEAS_EPS(path2[level * nnodes + k].dist, path2[(level + 1) * nnodes + k].dist, EPSILON));
2210  }
2211  }
2212  }
2213 #endif
2214
2215  return;
2216 }
2217
2218 /** computes 3 next terminal to all non-terminal nodes */
2220  GRAPH* g, /**< graph data structure */
2221  const SCIP_Real* cost, /**< edge costs */
2222  const SCIP_Real* costrev, /**< reversed edge costs */
2223  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
2224  TPATHS* tpaths /**< storage for terminal paths */
2225 )
2226 {
2227  assert(g != NULL);
2228  assert(cost != NULL);
2229  assert(costrev != NULL);
2230  assert(tpaths != NULL);
2231
2232  if( !graph_pc_isPcMw(g) )
2233  graph_mark(g);
2234
2236  graph_tpathsAdd2nd(g, cost, costrev, sdprofit, tpaths);
2237  graph_tpathsAdd3rd(g, cost, costrev, sdprofit, tpaths);
2238
2239 #ifndef NDEBUG
2240  if( !sdprofit )
2241  {
2242  const PATH* RESTRICT path3 = tpaths->termpaths;
2243  const int nnodes = graph_get_nNodes(g);
2244
2245  for( int level = 0; level < 2; level++ )
2246  {
2247  for( int k = 0; k < nnodes; ++k )
2248  {
2249  assert(LE_FEAS_EPS(path3[level * nnodes + k].dist, path3[(level + 1) * nnodes + k].dist, EPSILON));
2250  }
2251  }
2252  }
2253 #endif
2254
2255  return;
2256 }
2257
2258
2259 /** computes 4 next terminal to all non-terminal nodes */
2261  GRAPH* g, /**< graph data structure */
2262  const SCIP_Real* cost, /**< edge costs */
2263  const SCIP_Real* costrev, /**< reversed edge costs */
2264  const SDPROFIT* sdprofit, /**< SD bias for nodes or NULL */
2265  TPATHS* tpaths /**< storage for terminal paths */
2266 )
2267 {
2268  assert(g != NULL);
2269  assert(cost != NULL);
2270  assert(costrev != NULL);
2271
2272  graph_mark(g);
2273
2275  graph_tpathsAdd2nd(g, cost, costrev, sdprofit, tpaths);
2276  graph_tpathsAdd3rd(g, cost, costrev, sdprofit, tpaths);
2277  graph_tpathsAdd4th(g, cost, costrev, sdprofit, tpaths);
2278
2279 #ifndef NDEBUG
2280  if( !sdprofit )
2281  {
2282  PATH* RESTRICT path4 = tpaths->termpaths;
2283  const int nnodes = graph_get_nNodes(g);
2284
2285  for( int level = 0; level < 3; level++ )
2286  {
2287  for( int k = 0; k < nnodes; ++k )
2288  {
2289  assert(LE(path4[level * nnodes + k].dist, path4[(level + 1) * nnodes + k].dist));
2290  }
2291  }
2292  }
2293 #endif
2294
2295  return;
2296 }
2297
2298
2299 /** frees TPATHS structure */
2301  SCIP* scip, /**< SCIP */
2302  TPATHS** tpaths /**< the terminal paths */
2303 )
2304 {
2305  TPATHS* tp;
2306  assert(scip && tpaths);
2307
2308  tp = *tpaths;
2309  assert(tp);
2310
2311  SCIPfreeMemoryArrayNull(scip, &(tp->termbases_org));
2312  SCIPfreeMemoryArray(scip, &(tp->state));
2313  SCIPfreeMemoryArray(scip, &(tp->termbases));
2314  SCIPfreeMemoryArray(scip, &(tp->termpaths));
2315
2316  SCIPfreeMemory(scip, tpaths);
2317 }
2318
2319
2320 /** gets (up to) four close terminals to given node i;
2321  * with strict upper bound on allowed distances */
2323  const GRAPH* g, /**< graph */
2324  const TPATHS* tpaths, /**< the terminal paths */
2325  int node, /**< node */
2326  int* RESTRICT closeterm, /**< terminal */
2327  int* RESTRICT firstedge, /**< corresponding first edge (can be NULL) */
2328  SCIP_Real* RESTRICT closeterm_dist /**< terminal distance */
2329 )
2330 {
2331  assert(g && tpaths && closeterm && closeterm_dist);
2332  assert(graph_knot_isInRange(g, node));
2333
2334  if( Is_term(g->term[node]) )
2335  {
2336  *closeterm = node;
2337  *closeterm_dist = 0.0;
2338  if( firstedge )
2339  *firstedge = UNKNOWN;
2340  }
2341  else
2342  {
2343  const PATH* const termpaths = tpaths->termpaths;
2344  const int* const termbases = tpaths->termbases;
2345
2346  *closeterm = termbases[node];
2347  *closeterm_dist = termpaths[node].dist;
2348  if( firstedge )
2349  *firstedge = termpaths[node].edge;
2350  }
2351 }
2352
2353
2354 /** gets (up to) three close terminals to given node i;
2355  * with strict upper bound on allowed distances */
2357  const GRAPH* g, /**< graph */
2358  const TPATHS* tpaths, /**< the terminal paths */
2359  int node, /**< node */
2360  SCIP_Real maxdist_strict, /**< maximum valid distance (strict) */
2361  int* RESTRICT closeterms, /**< three close terminals */
2362  int* RESTRICT firstedges, /**< corresponding first edge (can be NULL) */
2363  SCIP_Real* RESTRICT closeterms_dist, /**< three close terminal distance */
2364  int* RESTRICT ncloseterms /**< number of close terminals found */
2365 )
2366 {
2367  const SCIP_Bool isStrict = TRUE;
2368  tpathsGetKCloseTerms(g, tpaths, 3, node, maxdist_strict, isStrict, closeterms, firstedges, closeterms_dist, ncloseterms);
2369 }
2370
2371
2372 /** gets (up to) four close terminals to given node i;
2373  * with strict upper bound on allowed distances */
2375  const GRAPH* g, /**< graph */
2376  const TPATHS* tpaths, /**< the terminal paths */
2377  int node, /**< node */
2378  SCIP_Real maxdist_strict, /**< maximum valid distance (strict) */
2379  int* RESTRICT closeterms, /**< four close terminals */
2380  int* RESTRICT firstedges, /**< corresponding first edge (can be NULL) */
2381  SCIP_Real* RESTRICT closeterms_dist, /**< four close terminal distance */
2382  int* RESTRICT ncloseterms /**< number of close terminals found */
2383 )
2384 {
2385  const SCIP_Bool isStrict = TRUE;
2386  tpathsGetKCloseTerms(g, tpaths, 4, node, maxdist_strict, isStrict, closeterms, firstedges, closeterms_dist, ncloseterms);
2387 }
2388
2389
2390 /** gets (up to) four close terminals to given node i
2391  * with non-strict upper bound on allowed distances */
2393  const GRAPH* g, /**< graph */
2394  const TPATHS* tpaths, /**< the terminal paths */
2395  int node, /**< node */
2396  SCIP_Real maxdist_nonstrict, /**< maximum valid distance (not strict) */
2397  int* RESTRICT closeterms, /**< four close terminals */
2398  int* RESTRICT firstedges, /**< corresponding first edge (can be NULL) */
2399  SCIP_Real* RESTRICT closeterms_dist, /**< four close terminal distance */
2400  int* RESTRICT ncloseterms /**< number of close terminals found */
2401 )
2402 {
2403  const SCIP_Bool isStrict = FALSE;
2404  tpathsGetKCloseTerms(g, tpaths, 4, node, maxdist_nonstrict, isStrict, closeterms, firstedges, closeterms_dist, ncloseterms);
2405 }
#define STP_Vectype(type)
Definition: stpvector.h:44
static void tpathsRepairExitLevel(SCIP *scip, int level, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:859
#define STP_TPATHS_RESERVESIZE
Definition: graph_tpath.c:36
void graph_get3nextTermPaths(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path3, int *vbase3, int *state3)
Definition: graph_tpath.c:1562
#define StpVecGetSize(vec)
Definition: stpvector.h:139
static void tpathsRepairExit(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1414
#define StpVecClear(vec)
Definition: stpvector.h:134
SCIP_Bool graph_pc_isMw(const GRAPH *)
Definition: graphdefs.h:224
static void tpathsScan4th(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, int heapsize, TPATHS *tpaths, TREPAIR *repair)
Definition: graph_tpath.c:764
int source
Definition: graphdefs.h:195
void graph_tpathsGet4CloseTerms(const GRAPH *g, const TPATHS *tpaths, int node, SCIP_Real maxdist_strict, int *RESTRICT closeterms, int *RESTRICT firstedges, SCIP_Real *RESTRICT closeterms_dist, int *RESTRICT ncloseterms)
Definition: graph_tpath.c:2374
#define SCIPfreeMemoryArrayNull(scip, ptr)
Definition: scip_mem.h:72
#define Is_term(a)
Definition: graphdefs.h:102
static void tpathsScan3rd(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, int heapsize, TPATHS *tpaths, TREPAIR *repair)
Definition: graph_tpath.c:700
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
#define EPSILON
Definition: lpi_qso.c:93
signed int edge
Definition: graphdefs.h:287
void graph_tpathsGet3CloseTerms(const GRAPH *g, const TPATHS *tpaths, int node, SCIP_Real maxdist_strict, int *RESTRICT closeterms, int *RESTRICT firstedges, SCIP_Real *RESTRICT closeterms_dist, int *RESTRICT ncloseterms)
Definition: graph_tpath.c:2356
void graph_tpathsAdd2nd(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:1987
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
void graph_add3rdTermPaths(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path3, int *vbase3, int *state3)
Definition: graph_tpath.c:1501
void graph_mark(GRAPH *)
Definition: graph_base.c:1130
static void tpathsResetMembers(const GRAPH *g, int level, TPATHS *tpaths)
Definition: graph_tpath.c:371
void graph_knot_printInfo(const GRAPH *, int)
Definition: graph_stats.c:151
#define FALSE
Definition: def.h:87
void graph_add4thTermPaths(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path4, int *vbase4, int *state4)
Definition: graph_tpath.c:1521
int *RESTRICT inpbeg
Definition: graphdefs.h:202
static SCIP_RETCODE tpathsRepair3rd(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1340
void graph_pathHeapAdd(const PATH *, int, int *RESTRICT, int *RESTRICT, int *)
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
includes various files containing graph methods used for Steiner tree problems
static void tpathsPrintPath(const GRAPH *g, const SDPROFIT *sdprofit, const TPATHS *tpaths, int node, int level)
Definition: graph_tpath.c:285
SCIP_RETCODE graph_tpathsRecomputeBiased(const SDPROFIT *sdprofit, GRAPH *g, TPATHS *tpaths)
Definition: graph_tpath.c:1776
void graph_tpathsGetClosestTerm(const GRAPH *g, const TPATHS *tpaths, int node, int *RESTRICT closeterm, int *RESTRICT firstedge, SCIP_Real *RESTRICT closeterm_dist)
Definition: graph_tpath.c:2322
#define StpVecPopBack(vec)
Definition: stpvector.h:182
#define EAT_LAST
Definition: graphdefs.h:58
#define StpVecReserve(scip, vec, _size_)
Definition: stpvector.h:186
#define SCIPdebugMessage
Definition: pub_message.h:87
#define FARAWAY
Definition: graphdefs.h:89
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
SCIP_RETCODE graph_tpathsInit(SCIP *scip, GRAPH *g, TPATHS **tpaths)
Definition: graph_tpath.c:1682
void graph_tpathsAdd3rd(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:2044
header only, simple implementation of an STL like vector
static SCIP_RETCODE tpathsAlloc(SCIP *scip, const GRAPH *g, TPATHS **tpaths)
Definition: graph_tpath.c:344
int *RESTRICT mark
Definition: graphdefs.h:198
static void utdist(SCIP *scip, const GRAPH *g, PATH *path, SCIP_Real ecost, int *vbase, int k, int l, int k2, int shift, int nnodes)
Definition: graph_tpath.c:188
static void tpathsRepairTraverseLevel(SCIP *scip, const GRAPH *g, int level, TREPAIR *repair)
Definition: graph_tpath.c:1117
static void tpathsRepairUpdate1st(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1135
#define SCIPallocCleanBufferArray(scip, ptr, num)
Definition: scip_mem.h:133
static void tpathsRepairUpdateLevel(SCIP *scip, const GRAPH *g, int level, TREPAIR *repair)
Definition: graph_tpath.c:1199
void graph_tpathsSetAll2(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:2180
int *RESTRICT oeat
Definition: graphdefs.h:231
SCIP_RETCODE graph_tpathsRepair(SCIP *scip, int edge, SCIP_Bool withEdgeReplacement, const GRAPH *g, TPATHS *tpaths)
Definition: graph_tpath.c:1744
#define LE(a, b)
Definition: portab.h:83
static void tpathsRepairTraverseLevelWithStack(SCIP *scip, const GRAPH *g, int level, TREPAIR *repair)
Definition: graph_tpath.c:963
#define GE(a, b)
Definition: portab.h:85
#define STP_TPATHS_NTERMBASES
Definition: graph_tpath.c:35
void graph_get2nextTermPaths(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path2, int *vbase2, int *state2)
Definition: graph_tpath.c:1542
#define StpVecGetcapacity(vec)
Definition: stpvector.h:129
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static void pathheapCorrectDist(int *RESTRICT heap, int *RESTRICT state, int *RESTRICT count, PATH *RESTRICT path, int l, int k, int edge, SCIP_Real newdist)
Definition: graph_tpath.c:145
static void tpathsBuild(SCIP *scip, GRAPH *g, TPATHS *tpaths)
Definition: graph_tpath.c:535
SCIP_Real dist
Definition: graphdefs.h:286
Definition: graphdefs.h:201
void graph_add1stTermPaths(const GRAPH *g, const SCIP_Real *cost, PATH *path, int *vbase, int *state)
Definition: graph_tpath.c:1464
#define NULL
Definition: lpi_spx1.cpp:155
static void tpathsScan1st(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SDPROFIT *sdprofit, int heapsize, TPATHS *tpaths, TREPAIR *repair)
Definition: graph_tpath.c:562
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
SCIP_VAR * h
Definition: circlepacking.c:59
#define LT(a, b)
Definition: portab.h:81
static SCIP_Real reduce_sdprofitGetProfit(const SDPROFIT *sdprofit, int node, int nonsource1, int nonsource2)
Definition: reducedefs.h:178
static int pathheapGetNearest(int *RESTRICT heap, int *RESTRICT state, int *RESTRICT count, const PATH *path)
Definition: graph_tpath.c:100
void graph_tpathsPrintForNode(const GRAPH *g, const SDPROFIT *sdprofit, const TPATHS *tpaths, int node)
Definition: graph_tpath.c:1793
static SCIP_RETCODE tpathsRepair2nd(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1318
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define GT(a, b)
Definition: portab.h:84
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT ieat
Definition: graphdefs.h:230
int *RESTRICT path_heap
Definition: graphdefs.h:255
SCIP_Bool graph_edge_isDeleted(const GRAPH *, int)
Definition: graph_stats.c:121
static SCIP_RETCODE tpathsRepair1st(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1287
includes graph definitions used for Steiner tree problems
SCIP_Bool graph_isMarked(const GRAPH *)
Definition: graph_base.c:1165
int *RESTRICT tail
Definition: graphdefs.h:223
void graph_add2ndTermPaths(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path2, int *vbase2, int *state2)
Definition: graph_tpath.c:1482
#define pathResetSingle(entry)
Definition: graph_tpath.c:39
#define flipedge(edge)
Definition: graphdefs.h:84
void graph_tpathsSetAll4(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:2260
#define tpathsRepairResetNode(scip, node, resetnodes, stack, nodes_isvisited)
Definition: graph_tpath.c:48
SCIP_RETCODE graph_get4nextTTerms(SCIP *scip, GRAPH *g, const SCIP_Real *cost, PATH *path, int *vbase, int *heap, int *state)
Definition: graph_tpath.c:1601
static void tpathsGetKCloseTerms(const GRAPH *g, const TPATHS *tpaths, int maxnumber, int node, SCIP_Real maxdist, SCIP_Bool distIsStrict, int *RESTRICT closeterms, int *RESTRICT firstedges, SCIP_Real *RESTRICT closeterms_dist, int *RESTRICT ncloseterms)
Definition: graph_tpath.c:405
int *RESTRICT term
Definition: graphdefs.h:196
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:127
void graph_edge_printInfo(const GRAPH *, int)
Definition: graph_stats.c:133
static void tpathsScan2nd(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, int heapsize, TPATHS *tpaths, TREPAIR *repair)
Definition: graph_tpath.c:637
#define StpVecIsEmpty(vec)
Definition: stpvector.h:148
#define CONNECT
Definition: graphdefs.h:87
void graph_tpathsGet4CloseTermsLE(const GRAPH *g, const TPATHS *tpaths, int node, SCIP_Real maxdist_nonstrict, int *RESTRICT closeterms, int *RESTRICT firstedges, SCIP_Real *RESTRICT closeterms_dist, int *RESTRICT ncloseterms)
Definition: graph_tpath.c:2392
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
SCIP_Real * r
Definition: circlepacking.c:50
static void tpathsRepairTraverseStackAddEdge(SCIP *scip, const GRAPH *g, int level, TREPAIR *repair)
Definition: graph_tpath.c:1021
void graph_tpathsAdd1st(const GRAPH *g, const SCIP_Real *cost, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:1932
static void tpathsRepairInit(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1386
static SCIP_RETCODE tpathsRepairInitLevel(SCIP *scip, int level, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:831
SCIP_Real * cost
Definition: graphdefs.h:221
SCIP_Bool graph_edge_isInRange(const GRAPH *, int)
Definition: graph_stats.c:110
void graph_tpathsGetProfitNodes(SCIP *scip, const GRAPH *g, const TPATHS *tpaths, const SDPROFIT *sdprofit, int start, int base, STP_Vectype(int) profitnodes)
Definition: graph_tpath.c:1842
SCIP_Bool graph_pc_isPcMw(const GRAPH *)
SCIP_VAR * a
Definition: circlepacking.c:57
static void tpathsBuildBiased(const SDPROFIT *sdprofit, GRAPH *g, TPATHS *tpaths)
Definition: graph_tpath.c:548
#define SCIP_Real
Definition: def.h:177
#define SCIPfreeCleanBufferArray(scip, ptr)
Definition: scip_mem.h:137
void graph_get4nextTermPaths(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path4, int *vbase4, int *state4)
Definition: graph_tpath.c:1582
#define StpVecFree(scip, vec)
Definition: stpvector.h:153
void graph_tpathsFree(SCIP *scip, TPATHS **tpaths)
Definition: graph_tpath.c:2300
static SCIP_Real tpathsGetDistNew(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const int *vbase, int node, int nextnode, int edge, const SDPROFIT *sdprofit, const PATH *path)
Definition: graph_tpath.c:484
struct tpaths_repair TREPAIR
int *RESTRICT outbeg
Definition: graphdefs.h:204
static SCIP_Real reduce_sdprofitGetBiasedDist(const SDPROFIT *sdprofit, int node, SCIP_Real edgecost, SCIP_Real nodedist, int nonsource1, int nonsource2)
Definition: reducedefs.h:214
int edges
Definition: graphdefs.h:219
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
#define UNKNOWN
Definition: sepa_mcf.c:4095
SCIP_RETCODE graph_tpathsRepairSetUp(const GRAPH *g, TPATHS *tpaths)
Definition: graph_tpath.c:1719
static void tpathsRepairTraverseStackAddBelow(SCIP *scip, const GRAPH *g, int level, TREPAIR *repair)
Definition: graph_tpath.c:1050
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:167
void graph_tpathsAdd4th(const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:2112
#define LE_FEAS_EPS(a, b, eps)
Definition: portab.h:75
static SCIP_RETCODE tpathsRepair4th(SCIP *scip, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:1363
includes various reduction methods for Steiner tree problems
SCIP_RETCODE graph_tpathsInitBiased(SCIP *scip, const SDPROFIT *sdprofit, GRAPH *g, TPATHS **tpaths)
Definition: graph_tpath.c:1700
static SCIP_RETCODE tpathsRepairTraverse1st(SCIP *scip, int start, int pred, const GRAPH *g, TREPAIR *repair)
Definition: graph_tpath.c:902
void graph_tpathsSetAll3(GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, const SDPROFIT *sdprofit, TPATHS *tpaths)
Definition: graph_tpath.c:2219