# SCIP

Solving Constraint Integer Programs

shortestpath.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 shortestpath.c
17  * @brief Shortest path based algorithms for Steiner problems
18  * @author Daniel Rehfeldt
19  *
20  * This file encompasses various shortest path based algorithms.
21  * Note: This file is supposed to replace graph_path.c in the long run, as it includes the faster implementations.
22  *
23  */
24
25 //#define SCIP_DEBUG
26 #include "reduce.h"
27 #include "shortestpath.h"
28 #include "portab.h"
29 #include "stpvector.h"
30
31
32 #ifndef NDEBUG
33 /** all terminals reached? */
34 static
36  const GRAPH* g, /**< graph data structure */
37  const STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
38 )
39 {
40  const int nnodes = graph_get_nNodes(g);
41
42  for( int k = 0; k < nnodes; ++k )
43  {
44  if( Is_term(g->term[k]) && !connected[k] )
45  {
46  SCIPdebugMessage("terminal %d not connected! \n", k);
47  return FALSE;
48  }
49  }
50
51  return TRUE;
52 }
53
54
55 /** all pseudo terminals reached? */
56 static
58  const GRAPH* g, /**< graph data structure */
59  const STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
60 )
61 {
62  const int nnodes = graph_get_nNodes(g);
63
64  assert(graph_pc_isPcMw(g));
65
66  for( int k = 0; k < nnodes; ++k )
67  {
68  if( Is_pseudoTerm(g->term[k]) && !connected[k] )
69  {
70  return FALSE;
71  }
72  }
73
74  return TRUE;
75 }
76
77
78 /** all fixed terminals reached? */
79 static
81  const GRAPH* g, /**< graph data structure */
82  const STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
83 )
84 {
85  const int nnodes = graph_get_nNodes(g);
86
87  assert(graph_pc_isPcMw(g));
88
89  for( int k = 0; k < nnodes; ++k )
90  {
91  if( graph_pc_knotIsFixedTerm(g, k) && !connected[k] )
92  {
93  return FALSE;
94  }
95  }
96
97  return TRUE;
98 }
99 #endif
100
101
102 /** initializes */
103 static inline
105  const GRAPH* g, /**< graph data structure */
106  int startnode, /**< start vertex */
107  SPATHS* spaths /**< shortest paths data */
108 )
109 {
110  const int nnodes = graph_get_nNodes(g);
111  DHEAP* dheap = spaths->dheap;
112  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
113  int* RESTRICT nodes_pred = spaths->nodes_pred;
114  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
115
116  assert(startnode >= 0 && startnode < g->knots);
117  assert(nnodes >= 1);
118  assert(spaths->csr->nnodes == nnodes);
119  assert(spaths->csr->start[nnodes] <= g->edges);
120
121  graph_heap_clean(TRUE, dheap);
122
123 #ifndef NDEBUG
124  assert(dheap->size == 0);
125  for( int k = 0; k < nnodes; k++ )
126  {
127  assert(dheap->position[k] == UNKNOWN);
128  nodes_pred[k] = -1;
129  }
130 #endif
131
132  for( int k = 0; k < nnodes; k++ )
133  {
134  nodes_dist[k] = FARAWAY;
135  connected[k] = FALSE;
136  }
137
138  nodes_dist[startnode] = 0.0;
139  nodes_pred[startnode] = -1;
140  graph_heap_correct(startnode, 0.0, dheap);
141  connected[startnode] = TRUE;
142 }
143
144
145 /** connects (also PC/MW potential) terminal to current tree */
146 static inline
148  const GRAPH* g, /**< graph data structure */
149  int k, /**< vertex to connect */
150  const int* nodes_pred, /**< predecessor array (on vertices) */
151  SCIP_Real* RESTRICT nodes_dist, /**< distance array (on vertices) */
152  DHEAP* dheap, /**< Dijkstra heap */
153  STP_Bool* RESTRICT connected /**< array to mark whether a vertex is part of computed Steiner tree */
154 )
155 {
156  assert(k >= 0 && k < g->knots);
157  assert(!connected[k]);
158
159  connected[k] = TRUE;
160  nodes_dist[k] = 0.0;
161
162  SCIPdebugMessage("connect terminal %d \n", k);
163
164  /* connect k to current solution */
165  for( int node = nodes_pred[k]; !connected[node]; node = nodes_pred[node] )
166  {
167  assert(node >= 0 && node < g->knots);
168  assert((!Is_term(g->term[node]) && !Is_pseudoTerm(g->term[node])));
169
170  SCIPdebugMessage("connect path node %d \n", node);
171
172  connected[node] = TRUE;
173  nodes_dist[node] = 0.0;
174  graph_heap_correct(node, 0.0, dheap);
175  }
176 }
177
178
179
180 /** connects node to current tree */
181 static inline
183  const GRAPH* g, /**< graph data structure */
184  int k, /**< vertex to connect */
185  const int* nodes_pred, /**< predecessor array (on vertices) */
186  SCIP_Real* RESTRICT nodes_dist, /**< distance array (on vertices) */
187  DHEAP* dheap, /**< Dijkstra heap */
188  int* termscount, /**< pointer for terminal count */
189  STP_Bool* RESTRICT connected /**< array to mark whether a vertex is part of computed Steiner tree */
190 )
191 {
192  assert(k >= 0 && k < g->knots);
193  assert(!connected[k]);
194
195  SCIPdebugMessage("connect terminal %d (dist=%f) \n", k, nodes_dist[k]);
196
197  connected[k] = TRUE;
198  nodes_dist[k] = 0.0;
199
200  if( Is_term(g->term[k]) )
201  (*termscount)++;
202
203  /* connect k to current solution */
204  for( int node = nodes_pred[k]; !connected[node]; node = nodes_pred[node] )
205  {
206  assert(node >= 0 && node < g->knots);
207
208  SCIPdebugMessage("connect path node %d \n", node);
209
210  if( Is_term(g->term[node]) )
211  (*termscount)++;
212
213  connected[node] = TRUE;
214  nodes_dist[node] = 0.0;
215  graph_heap_correct(node, 0.0, dheap);
216  }
217
218  assert(*termscount <= g->terms);
219 }
220
221
222 /** executes */
223 static inline
225  const GRAPH* g, /**< graph data structure */
226  int startnode, /**< start vertex */
227  SPATHS* spaths /**< shortest paths data */
228 )
229 {
230  const CSR* const csr = spaths->csr;
231  DHEAP* const dheap = spaths->dheap;
232  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
233  int* RESTRICT nodes_pred = spaths->nodes_pred;
234  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
235  int* const state = dheap->position;
236  const SCIP_Real* const cost_csr = csr->cost;
238  const int* const start_csr = csr->start;
239  int termscount = 0;
240
241  assert(dheap->size == 1);
242  assert(connected[startnode]);
243
244  if( Is_term(g->term[startnode]) )
245  termscount++;
246
247  /* main loop */
248  while( dheap->size > 0 )
249  {
250  /* get nearest labelled node */
251  const int k = graph_heap_deleteMinReturnNode(dheap);
252  const int k_start = start_csr[k];
253  const int k_end = start_csr[k + 1];
254  register SCIP_Real k_dist;
255
256  SCIPdebugMessage("take node %d from queue \n", k);
257
258  assert(state[k] == CONNECT);
259  /* NOTE: needs to be set for invariant of heap */
260  state[k] = UNKNOWN;
261
262  if( Is_term(g->term[k]) && k != startnode )
263  {
264  assert(!connected[k]);
265  computeSteinerTree_connectTerminal(g, k, nodes_pred, nodes_dist, dheap, connected);
266
267  assert(termscount < g->terms);
268
269  /* have all terminals been reached? */
270  if( ++termscount == g->terms )
271  {
272  break;
273  }
274  }
275
276  k_dist = nodes_dist[k];
277
278  for( int e = k_start; e != k_end; e++ )
279  {
280  const int m = head_csr[e];
281
282  if( !connected[m] )
283  {
284  const SCIP_Real distnew = k_dist + cost_csr[e];
285
286  /* closer to k than to current predecessor? */
287  if( LT(distnew, nodes_dist[m]) )
288  {
289  nodes_pred[m] = k;
290  nodes_dist[m] = distnew;
291  graph_heap_correct(m, distnew, dheap);
292  }
293  }
294  }
295  }
296 }
297
298
299
300
301 /** Executes directed Steiner tree computation.
302  * Here we always start from the root, but connect the 'startnode' vertex first */
303 static inline
305  SCIP* scip, /**< SCIP data structure */
306  const GRAPH* g, /**< graph data structure */
307  int startnode, /**< start vertex */
308  SPATHS* spaths /**< shortest paths data */
309 )
310 {
311  STP_Vectype(int) lostterms = NULL;
312  const CSR* const csr = spaths->csr;
313  DHEAP* const dheap = spaths->dheap;
314  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
315  int* RESTRICT nodes_pred = spaths->nodes_pred;
316  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
317  int* const state = dheap->position;
318  const SCIP_Real* const cost_csr = csr->cost;
320  const int* const start_csr = csr->start;
321  int termscount = 1;
322
323  assert(dheap->size == 1);
324  assert(graph_typeIsDirected(g));
325  assert(connected[g->source]);
326
327  /* main loop */
328  while( dheap->size > 0 )
329  {
330  /* get nearest labelled node */
331  const int k = graph_heap_deleteMinReturnNode(dheap);
332  const int k_start = start_csr[k];
333  const int k_end = start_csr[k + 1];
334  register SCIP_Real k_dist;
335
336  SCIPdebugMessage("take node %d from queue \n", k);
337
338  assert(state[k] == CONNECT);
339  /* NOTE: needs to be set for invariant of heap */
340  state[k] = UNKNOWN;
341
342  if( (Is_term(g->term[k]) || k == startnode) && !connected[k] )
343  {
344  assert(k != g->source);
345
346  /* NOTE: we don't want to connect terminals until the startnode has been added */
347  if( connected[startnode] || k == startnode )
348  {
349  assert(termscount < g->terms);
350  computeSteinerTree_connectNode(g, k, nodes_pred, nodes_dist, dheap, &termscount, connected);
351
352  if( termscount == g->terms )
353  {
354  break;
355  }
356  }
357  else if( k != startnode )
358  {
359  assert(Is_term(g->term[k]));
360  assert(!connected[k]);
361
362  StpVecPushBack(scip, lostterms, k);
363  }
364  }
365
366  k_dist = nodes_dist[k];
367
368  for( int e = k_start; e != k_end; e++ )
369  {
370  const int m = head_csr[e];
371
372  if( !connected[m] )
373  {
374  const SCIP_Real distnew = k_dist + cost_csr[e];
375
376  /* closer to k than to current predecessor? */
377  if( LT(distnew, nodes_dist[m]) )
378  {
379  nodes_pred[m] = k;
380  nodes_dist[m] = distnew;
381  graph_heap_correct(m, distnew, dheap);
382  }
383  }
384  }
385  }
386
387  for( int i = 0; i < StpVecGetSize(lostterms); i++ )
388  {
389  const int term = lostterms[i];
390  assert(Is_term(g->term[term]));
391
392  if( !connected[term] )
393  {
394  assert(termscount < g->terms);
395  computeSteinerTree_connectNode(g, term, nodes_pred, nodes_dist, dheap, &termscount, connected);
396  }
397  }
398
399  assert(termscount == g->terms);
400
401  StpVecFree(scip, lostterms);
402 }
403
404
405
406 /** executes */
407 static inline
409  const GRAPH* g, /**< graph data structure */
410  const SDPROFIT* sdprofit, /**< implied profit data structure */
411  int startnode, /**< start vertex */
412  SPATHS* spaths /**< shortest paths data */
413 )
414 {
415  const CSR* const csr = spaths->csr;
416  DHEAP* const dheap = spaths->dheap;
417  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
418  int* RESTRICT nodes_pred = spaths->nodes_pred;
419  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
420  int* const state = dheap->position;
421  const SCIP_Real* const cost_csr = csr->cost;
423  const int* const start_csr = csr->start;
424  const int* const nodes_biassource = sdprofit->nodes_biassource;
425  int termscount = 0;
426  const SCIP_Real zeroOffset = spaths->edgecost_zeroOffset;
427
428  assert(dheap->size == 1);
429  assert(connected[startnode]);
430  assert(zeroOffset > 0.0);
431
432  if( Is_term(g->term[startnode]) )
433  termscount++;
434
435  /* main loop */
436  while( dheap->size > 0 )
437  {
438  /* get nearest labelled node */
439  const int k = graph_heap_deleteMinReturnNode(dheap);
440  const int k_start = start_csr[k];
441  const int k_end = start_csr[k + 1];
442  register SCIP_Real k_dist;
443
444  SCIPdebugMessage("take node %d from queue \n", k);
445
446  assert(state[k] == CONNECT);
447  /* NOTE: needs to be set for invariant of heap */
448  state[k] = UNKNOWN;
449
450  if( Is_term(g->term[k]) && k != startnode )
451  {
452  assert(!connected[k]);
453  computeSteinerTree_connectTerminal(g, k, nodes_pred, nodes_dist, dheap, connected);
454
455  assert(termscount < g->terms);
456
457  /* have all terminals been reached? */
458  if( ++termscount == g->terms )
459  {
460  break;
461  }
462  }
463
464  k_dist = nodes_dist[k];
465
466  for( int e = k_start; e != k_end; e++ )
467  {
468  const int m = head_csr[e];
469
470  if( !connected[m] )
471  {
472  SCIP_Real distnew = k_dist + cost_csr[e];
473  const int source = nodes_biassource[k];
474  assert(graph_knot_isInRange(g, source));
475
476  if( source != k )
477  {
478  /* NOTE: with need to be careful with numerical issues here... */
479  if( source != m && !connected[source] && k != nodes_pred[m] )
480  {
481  SCIP_Real profitBias = sdprofit->nodes_bias[k];
482  profitBias = MIN(profitBias, cost_csr[e]);
483  profitBias = MIN(profitBias, k_dist);
484
485  if( GT(profitBias, zeroOffset) )
486  distnew -= profitBias;
487  }
488  }
489  else
490  {
491  assert(EQ(sdprofit->nodes_bias[k], 0.0) || EQ(sdprofit->nodes_bias[k], FARAWAY));
492  }
493
494  /* closer to k than to current predecessor? */
495  if( LT_HARD(distnew, nodes_dist[m]) )
496  {
497  nodes_pred[m] = k;
498  nodes_dist[m] = distnew;
499  graph_heap_correct(m, distnew, dheap);
500  }
501  }
502  }
503  }
504
505  assert(!graph_typeIsSpgLike(g) || termscount == g->terms);
506 }
507
508
509 /** connects node to current tree */
510 static inline
512  const GRAPH* g, /**< graph data structure */
513  int k, /**< vertex to connect */
514  const int* nodes_pred, /**< predecessor array (on vertices) */
515  SCIP_Real* RESTRICT nodes_dist, /**< distance array (on vertices) */
516  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
517  DHEAP* dheap, /**< Dijkstra heap */
518  STP_Bool* RESTRICT connected, /**< array to mark whether a vertex is part of computed Steiner tree */
519  int* RESTRICT nPseudoTerms /**< pointer */
520 )
521 {
522  assert(k >= 0 && k < g->knots);
523  assert(!connected[k]);
524  assert(Is_pseudoTerm(g->term[k]));
525
526  connected[k] = TRUE;
527  nodes_dist[k] = 0.0;
528  shortestpath_pcConnectNode(g, connected, k, spaths_pc);
529  (*nPseudoTerms)++;
530
531  SCIPdebugMessage("connect pseudo terminal %d \n", k);
532
533  /* connect k to current solution */
534  for( int node = nodes_pred[k]; !connected[node]; node = nodes_pred[node] )
535  {
536  assert(node >= 0 && node < g->knots);
537  assert(!connected[node]);
538  assert(!Is_term(g->term[node]));
539
540  SCIPdebugMessage("connect path node %d \n", node);
541
542  if( Is_pseudoTerm(g->term[node]) )
543  {
544  shortestpath_pcConnectNode(g, connected, node, spaths_pc);
545  (*nPseudoTerms)++;
546  }
547
548  connected[node] = TRUE;
549  nodes_dist[node] = 0.0;
550  graph_heap_correct(node, 0.0, dheap);
551  }
552 }
553
554 /** connects node to current tree */
555 static inline
557  const GRAPH* g, /**< graph data structure */
558  int k, /**< vertex to connect */
559  const SCIP_Real* prize, /**< (possibly biased) prize */
560  SCIP_Bool costIsBiased, /**< is cost biased? */
561  const int* nodes_pred, /**< predecessor array (on vertices) */
562  SCIP_Real* RESTRICT nodes_dist, /**< distance array (on vertices) */
563  DHEAP* dheap, /**< Dijkstra heap */
564  STP_Bool* RESTRICT connected, /**< array to mark whether a vertex is part of computed Steiner tree */
565  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
566  SCIP_Bool* isConnected, /**< node connected to tree? */
567  int* RESTRICT nPseudoTerms /**< pointer */
568 )
569 {
570  SCIP_Bool connectK = FALSE;
571
572  /* if k is positive vertex and close enough, connect k to current subtree */
573  if( !connected[k] && Is_pseudoTerm(g->term[k]) )
574  {
575  connectK = (GE(prize[k], nodes_dist[k]));
576
577  /* maybe if we count the prizes on the path, the extension becomes profitable? */
578  if( !connectK )
579  {
580  SCIP_Real prizesum = 0.0;
581  const SCIP_Bool isPc = graph_pc_isPc(g);
582
583  for( int node = nodes_pred[k]; !connected[node]; node = nodes_pred[node] )
584  {
585  if( Is_pseudoTerm(g->term[node]) )
586  prizesum += prize[node];
587  else if( isPc && !costIsBiased && graph_pc_knotIsNonLeafTerm(g, node) )
588  prizesum += prize[node];
589  }
590
591  assert(prizesum >= 0.0 && LT(prizesum, FARAWAY));
592  connectK = GE(prize[k] + prizesum, nodes_dist[k]);
593  }
594
595  if( connectK )
596  {
597  computeSteinerTree_connectPseudoTerm(g, k, nodes_pred, nodes_dist, spaths_pc, dheap, connected, nPseudoTerms);
598  }
599  }
600
601  *isConnected = connectK;
602 }
603
604
605 /** executes */
606 static inline
608  const GRAPH* g, /**< graph data structure */
609  int startnode, /**< start vertex */
610  const SCIP_Real* prize, /**< (possibly biased) prize */
611  SCIP_Bool costIsBiased, /**< is cost biased? */
612  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
613  SPATHS* spaths /**< shortest paths data */
614 )
615 {
616  const CSR* const csr = spaths->csr;
617  DHEAP* const dheap = spaths->dheap;
618  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
619  int* RESTRICT nodes_pred = spaths->nodes_pred;
620  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
621  int* const state = dheap->position;
622  const SCIP_Real* const cost_csr = csr->cost;
624  const int* const start_csr = csr->start;
625  int termsposCount = 0;
626  const int ntermspos = graph_pc_nNonFixedTerms(g);
627
628  assert(g->extended);
629  assert(dheap->size == 1);
630  assert(connected[startnode]);
631  assert(graph_pc_isPcMw(g) && !graph_pc_isRootedPcMw(g));
632
633  shortestpath_pcReset(spaths_pc);
634
635  if( Is_pseudoTerm(g->term[startnode]) )
636  {
637  termsposCount++;
638  shortestpath_pcConnectNode(g, connected, startnode, spaths_pc);
639  }
640
641  /* main loop */
642  while( dheap->size > 0 )
643  {
644  /* get nearest labelled node */
645  const int k = graph_heap_deleteMinReturnNode(dheap);
646  const int k_start = start_csr[k];
647  const int k_end = start_csr[k + 1];
648  register SCIP_Real k_dist;
649  SCIP_Bool k_isConnected = FALSE;
650
651  SCIPdebugMessage("take node %d from queue \n", k);
652
653  assert(state[k] == CONNECT);
654  state[k] = UNKNOWN;
655
656  computeSteinerTree_tryConnectNodePcMw(g, k, prize, costIsBiased, nodes_pred, nodes_dist, dheap, connected,
657  spaths_pc, &k_isConnected, &termsposCount);
658
659  k_dist = nodes_dist[k];
660
661  if( k_isConnected && termsposCount == ntermspos )
662  {
663  SCIPdebugMessage("all potential terminals reached \n");
664  assert(computeSteinerTree_allPseudoTermsAreReached(g, connected));
665  break;
666  }
667  else if( !k_isConnected && k_dist > spaths_pc->maxoutprize )
668  {
669  SCIPdebugMessage("further extension is not profitable \n");
670  break;
671  }
672
673  for( int e = k_start; e != k_end; e++ )
674  {
675  const int m = head_csr[e];
676
677  if( !connected[m] )
678  {
679  const SCIP_Real distnew = k_dist + cost_csr[e];
680
681  /* closer to k than to current predecessor? */
682  if( LT(distnew, nodes_dist[m]) )
683  {
684  nodes_pred[m] = k;
685  nodes_dist[m] = distnew;
686  graph_heap_correct(m, distnew, dheap);
687  }
688  }
689  }
690  }
691
692  assert(termsposCount == ntermspos || !computeSteinerTree_allPseudoTermsAreReached(g, connected));
693 }
694
695
696
697 /** executes */
698 static inline
700  const GRAPH* g, /**< graph data structure */
701  int startnode, /**< start vertex */
702  const SCIP_Real* prize, /**< (possibly biased) prize */
703  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
704  SPATHS* spaths /**< shortest paths data */
705 )
706 {
707  const CSR* const csr = spaths->csr;
708  DHEAP* const dheap = spaths->dheap;
709  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
710  int* RESTRICT nodes_pred = spaths->nodes_pred;
711  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
712  int* const state = dheap->position;
713  const SCIP_Real* const cost_csr = csr->cost;
715  const int* const start_csr = csr->start;
716  const int nfixedterms = graph_pc_nFixedTerms(g);
717  int fixedtermsCount = 0;
718
719  assert(g->extended);
720  assert(dheap->size == 1 && connected[startnode]);
721  assert(graph_pc_isRootedPcMw(g));
722
723  shortestpath_pcReset(spaths_pc);
724
725  if( Is_pseudoTerm(g->term[startnode]) )
726  shortestpath_pcConnectNode(g, connected, startnode, spaths_pc);
727
728  if( Is_term(g->term[startnode]) )
729  {
730  assert(graph_pc_knotIsFixedTerm(g, startnode));
731  fixedtermsCount++;
732  }
733
734  while( dheap->size > 0 )
735  {
736  const int k = graph_heap_deleteMinReturnNode(dheap);
737  const int k_start = start_csr[k];
738  const int k_end = start_csr[k + 1];
739  register SCIP_Real k_dist;
740
741  assert(state[k] == CONNECT);
742  state[k] = UNKNOWN;
743
744  /* if k is fixed terminal positive vertex and close enough, connect its path to current subtree */
745  if( Is_anyTerm(g->term[k]) && (Is_term(g->term[k]) || GE(prize[k], nodes_dist[k])) && !connected[k] )
746  {
747  assert(!graph_pc_knotIsDummyTerm(g, k));
748  assert(graph_pc_knotIsFixedTerm(g, k) || GE(prize[k], nodes_dist[k]));
749  assert(graph_pc_knotIsFixedTerm(g, k) == Is_term(g->term[k]));
750
751  if( Is_term(g->term[k]) )
752  fixedtermsCount++;
753  else if( Is_pseudoTerm(g->term[k]) )
754  shortestpath_pcConnectNode(g, connected, k, spaths_pc);
755
756  connected[k] = TRUE;
757  nodes_dist[k] = 0.0;
758
759  assert(nodes_pred[k] >= 0);
760
761  for( int node = nodes_pred[k]; !connected[node]; node = nodes_pred[node] )
762  {
763  assert(node >= 0 && node < g->knots);
764  SCIPdebugMessage("connect path node %d \n", node);
765
766  if( Is_pseudoTerm(g->term[node]) )
767  shortestpath_pcConnectNode(g, connected, node, spaths_pc);
768
769  connected[node] = TRUE;
770  nodes_dist[node] = 0.0;
771  graph_heap_correct(node, 0.0, dheap);
772  }
773  }
774
775  k_dist = nodes_dist[k];
776  assert(!graph_pc_knotIsFixedTerm(g, k) || EQ(k_dist, 0.0));
777
778  if( fixedtermsCount >= nfixedterms && k_dist > spaths_pc->maxoutprize )
779  {
780  SCIPdebugMessage("all fixed terminals reached \n");
781  assert(fixedtermsCount == nfixedterms);
782  break;
783  }
784
785  for( int e = k_start; e != k_end; e++ )
786  {
787  const int m = head_csr[e];
788
789  if( !connected[m] )
790  {
791  const SCIP_Real distnew = k_dist + cost_csr[e];
792
793  if( distnew < nodes_dist[m] )
794  {
795  nodes_pred[m] = k;
796  nodes_dist[m] = distnew;
797  graph_heap_correct(m, distnew, dheap);
798  }
799  }
800  }
801  }
802
803  assert(computeSteinerTree_allFixedTermsAreReached(g, connected));
804 }
805
806
807 /** executes */
808 static inline
810  const GRAPH* g, /**< graph data structure */
811  int startnode, /**< start vertex */
812  SPATHS* spaths /**< shortest paths data */
813 )
814 {
815  const CSR* const csr = spaths->csr;
816  DHEAP* const dheap = spaths->dheap;
817  SCIP_Real* RESTRICT nodes_dist = spaths->nodes_dist;
818  int* RESTRICT nodes_pred = spaths->nodes_pred;
819  STP_Bool* RESTRICT connected = spaths->nodes_isConnected;
820  int* const state = dheap->position;
821  const SCIP_Real* const cost_csr = csr->cost;
823  const int* const start_csr = csr->start;
824  int termscount = 0;
825  const int nterms = graph_pc_isRootedPcMw(g) ? g->terms : g->terms - 1;
826
827  assert(g->extended && graph_pc_isPcMw(g));
828  assert(dheap->size == 1);
829  assert(connected[startnode]);
830
831  if( Is_term(g->term[startnode]) || Is_pseudoTerm(g->term[startnode]) )
832  termscount++;
833
834  /* main loop */
835  while( dheap->size > 0 )
836  {
837  /* get nearest labelled node */
838  const int k = graph_heap_deleteMinReturnNode(dheap);
839  const int k_start = start_csr[k];
840  const int k_end = start_csr[k + 1];
841  register SCIP_Real k_dist;
842
843  SCIPdebugMessage("take node %d from queue \n", k);
844
845  assert(state[k] == CONNECT);
846  state[k] = UNKNOWN;
847
848  /* connect k to current subtree? */
849  if( !connected[k] && (Is_term(g->term[k]) || Is_pseudoTerm(g->term[k])) )
850  {
851  computeSteinerTree_connectTerminal(g, k, nodes_pred, nodes_dist, dheap, connected);
852
853  assert(termscount < nterms);
854
855  /* have all terminals been reached? */
856  if( ++termscount == nterms )
857  {
858  break;
859  }
860  }
861
862  k_dist = nodes_dist[k];
863
864  for( int e = k_start; e != k_end; e++ )
865  {
866  const int m = head_csr[e];
867
868  if( !connected[m] )
869  {
870  const SCIP_Real distnew = k_dist + cost_csr[e];
871
872  /* closer to k than to current predecessor? */
873  if( LT(distnew, nodes_dist[m]) )
874  {
875  nodes_pred[m] = k;
876  nodes_dist[m] = distnew;
877  graph_heap_correct(m, distnew, dheap);
878  }
879  }
880  }
881  }
882
884 }
885
886
887 /*
888  * Interface methods
889  */
890
891
892 /** initializes */
894  SCIP* scip, /**< SCIP data structure */
895  const GRAPH* graph, /**< graph data structure */
896  const SCIP_Real* costs, /**< cost for edges (might be negative for MWCS or RMWCS) */
897  const SCIP_Real* prizes, /**< prizes for all nodes */
898  SPATHSPC** sppc /**< PC/MW shortest path data */
899  )
900 {
901  SCIP_Real* orderedprizes;
902  int* orderedprizes_id;
903  const int nnodes = graph_get_nNodes(graph);
904  const int nterms = graph->terms;
905  const SCIP_Bool usecosts = (costs && (graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP));
906  int termcount;
907
908  assert(graph_pc_isPcMw(graph));
909  assert(prizes && graph);
910  assert(graph->extended);
911  assert(nterms >= 1);
912
913  SCIP_CALL( SCIPallocMemoryArray(scip, &orderedprizes, nterms + 1) );
914  SCIP_CALL( SCIPallocMemoryArray(scip, &orderedprizes_id, nterms + 1) );
915  SCIP_CALL( SCIPallocMemory(scip, sppc) );
916
917  termcount = 0;
918  for( int k = 0; k < nnodes; k++ )
919  {
920  if( Is_pseudoTerm(graph->term[k]) )
921  {
922  orderedprizes[termcount] = prizes[k];
923
924  /* consider incoming negative arcs */
925  if( usecosts )
926  {
927  SCIP_Real mincost = 0.0;
928  for( int e = graph->inpbeg[k]; e != EAT_LAST; e = graph->ieat[e] )
929  if( costs[e] < mincost )
930  mincost = costs[e];
931
932  if( mincost < 0.0 )
933  orderedprizes[termcount] -= mincost;
934  }
935
936  orderedprizes_id[termcount++] = k;
937  }
938  }
939
940  for( int k = termcount; k < nterms; k++ )
941  {
942  orderedprizes[k] = 0.0;
943  orderedprizes_id[k] = -1;
944  }
945
946  SCIPsortDownRealInt(orderedprizes, orderedprizes_id, nterms);
947
948  /* set sentinel */
949  orderedprizes[nterms] = 0.0;
950  orderedprizes_id[nterms] = -1;
951
952  assert(orderedprizes[0] >= orderedprizes[nterms - 1]);
953
954  (*sppc)->orderedprizes = orderedprizes;
955  (*sppc)->orderedprizes_id = orderedprizes_id;
956  (*sppc)->maxoutprize_idx = -1;
957  (*sppc)->maxoutprize = -FARAWAY;
958
959  return SCIP_OKAY;
960 }
961
962
963 /** frees */
965  SCIP* scip, /**< SCIP data structure */
966  SPATHSPC** sppc /**< PC/MW shortest path data */
967  )
968 {
969  SCIPfreeMemoryArray(scip, &((*sppc)->orderedprizes_id));
970  SCIPfreeMemoryArray(scip, &((*sppc)->orderedprizes));
971
972  SCIPfreeMemory(scip, sppc);
973 }
974
975
976 /** start computation */
978  SPATHSPC* sppc /**< PC/MW shortest path data */
979 )
980 {
981  SCIP_Real* orderedprizes = sppc->orderedprizes;
982
983  assert(orderedprizes);
984
985  sppc->maxoutprize = orderedprizes[0];
986  sppc->maxoutprize_idx = 0;
987 }
988
989
990 /** update maximum prize */
992  const GRAPH* g, /**< graph data structure */
993  const STP_Bool* connected, /**< array to mark whether a vertex is part of computed Steiner tree */
994  int node_connected, /**< node that is removed */
995  SPATHSPC* sppc /**< PC data */
996 )
997 {
998  const SCIP_Real* orderedprizes = sppc->orderedprizes;
999  const int* orderedprizes_id = sppc->orderedprizes_id;
1000  int maxprizeidx = sppc->maxoutprize_idx;
1001  assert(maxprizeidx >= 0 && maxprizeidx <= g->terms);
1002
1003  /* sentinel? */
1004  if( orderedprizes_id[maxprizeidx] < 0 )
1005  {
1006  assert(orderedprizes_id[maxprizeidx] == -1);
1007  assert(EQ(sppc->maxoutprize, 0.0));
1008  return;
1009  }
1010
1011  assert(maxprizeidx < g->terms);
1012
1013  /* is current node at the maximum? */
1014  if( node_connected == orderedprizes_id[maxprizeidx] )
1015  {
1016  while( orderedprizes_id[maxprizeidx] >= 0 && connected[orderedprizes_id[maxprizeidx]] )
1017  {
1018  maxprizeidx++;
1019  assert(maxprizeidx <= g->terms);
1020  }
1021
1022  sppc->maxoutprize_idx = maxprizeidx;
1023
1024  if( orderedprizes_id[maxprizeidx] < 0 )
1025  sppc->maxoutprize = 0.0;
1026  else
1027  sppc->maxoutprize = orderedprizes[maxprizeidx];
1028  }
1029 }
1030
1031
1032 /** shortest path based heuristic for computing a Steiner tree */
1034  const GRAPH* g, /**< graph data structure */
1035  int startnode, /**< start vertex */
1036  SPATHS* spaths /**< shortest paths data */
1037 )
1038 {
1039  assert(g && spaths);
1040  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1041  assert(graph_typeIsSpgLike(g));
1042
1043  computeSteinerTree_init(g, startnode, spaths);
1044
1045  if( g->knots == 1 )
1046  return;
1047
1048  computeSteinerTree_exec(g, startnode, spaths);
1049
1051 }
1052
1053
1054 /** shortest path based heuristic for computing a Steiner tree */
1056  SCIP* scip, /**< SCIP data structure */
1057  const GRAPH* g, /**< graph data structure */
1058  int startnode, /**< start vertex */
1059  SPATHS* spaths /**< shortest paths data */
1060 )
1061 {
1062  assert(g && spaths);
1063  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1064  assert(graph_typeIsDirected(g));
1065
1066  /* NOTE: here we treat g->source as the start */
1067  computeSteinerTree_init(g, g->source, spaths);
1068
1069  if( g->knots == 1 )
1070  return;
1071
1072  SCIPdebugMessage("startnode=%d \n", startnode);
1073
1074  computeSteinerTree_execDirected(scip, g, startnode, spaths);
1075
1077 }
1078
1079
1080 /** shortest path based heuristic for computing a Steiner tree */
1082  const GRAPH* g, /**< graph data structure */
1083  const SDPROFIT* sdprofit, /**< implied profit data structure */
1084  int startnode, /**< start vertex */
1085  SPATHS* spaths /**< shortest paths data */
1086 )
1087 {
1088  assert(g && spaths && sdprofit);
1089  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1090  assert(graph_typeIsSpgLike(g));
1091
1092  computeSteinerTree_init(g, startnode, spaths);
1093
1094  if( g->knots == 1 )
1095  return;
1096
1097  computeSteinerTree_execBiased(g, sdprofit, startnode, spaths);
1098
1100 }
1101
1102
1103 /** shortest path based heuristic for computing a Steiner tree in PC/MW case */
1105  const GRAPH* g, /**< graph data structure */
1106  int startnode, /**< start vertex */
1107  const SCIP_Real* prize, /**< (possibly biased) prize */
1108  SCIP_Bool costIsBiased, /**< is cost biased? */
1109  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
1110  SPATHS* spaths /**< shortest paths data */
1111 )
1112 {
1113  assert(g && spaths);
1114  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1115  assert(graph_pc_isPcMw(g) && !graph_pc_isRootedPcMw(g));
1116  assert(g->extended);
1117  assert(!graph_pc_knotIsDummyTerm(g, startnode));
1118
1119  computeSteinerTree_init(g, startnode, spaths);
1120
1121  if( g->knots == 1 )
1122  return;
1123
1124  computeSteinerTree_execPcMw(g, startnode, prize, costIsBiased, spaths_pc, spaths);
1125 }
1126
1127
1128 /** shortest path based heuristic for computing a Steiner tree in rooted PC/MW case */
1130  const GRAPH* g, /**< graph data structure */
1131  int startnode, /**< start vertex */
1132  const SCIP_Real* prize, /**< (possibly biased) prize */
1133  SPATHSPC* spaths_pc, /**< PC/MW shortest paths data */
1134  SPATHS* spaths /**< shortest paths data */
1135 )
1136 {
1137  assert(g && spaths);
1138  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1139  assert(graph_pc_isRootedPcMw(g) && g->extended);
1140  assert(!graph_pc_knotIsDummyTerm(g, startnode));
1141
1142  computeSteinerTree_init(g, startnode, spaths);
1143
1144  if( g->knots == 1 )
1145  return;
1146
1147  computeSteinerTree_execRpcMw(g, startnode, prize, spaths_pc, spaths);
1148 }
1149
1150
1151 /** shortest path based heuristic for computing a Steiner tree in PC/MW case
1152  * that contains all (potential and fixed) terminals */
1154  const GRAPH* g, /**< graph data structure */
1155  int startnode, /**< start vertex */
1156  SPATHS* spaths /**< shortest paths data */
1157 )
1158 {
1159  assert(g && spaths);
1160  assert(spaths->csr && spaths->nodes_dist && spaths->nodes_pred && spaths->dheap && spaths->nodes_isConnected);
1161  assert(graph_pc_isPcMw(g) && g->extended);
1162  assert(!graph_pc_knotIsDummyTerm(g, startnode));
1163
1164  computeSteinerTree_init(g, startnode, spaths);
1165
1166  if( g->knots == 1 )
1167  return;
1168
1169  computeSteinerTree_execPcMwFull(g, startnode, spaths);
1170 }
#define STP_Vectype(type)
Definition: stpvector.h:44
static volatile int nterms
Definition: interrupt.c:38
const SCIP_Real edgecost_zeroOffset
Definition: shortestpath.h:59
#define StpVecGetSize(vec)
Definition: stpvector.h:139
void shortestpath_computeSteinerTreePcMw(const GRAPH *g, int startnode, const SCIP_Real *prize, SCIP_Bool costIsBiased, SPATHSPC *spaths_pc, SPATHS *spaths)
SCIP_Bool graph_typeIsDirected(const GRAPH *)
Definition: graph_stats.c:83
Definition: graphdefs.h:141
int source
Definition: graphdefs.h:195
static void computeSteinerTree_connectTerminal(const GRAPH *g, int k, const int *nodes_pred, SCIP_Real *RESTRICT nodes_dist, DHEAP *dheap, STP_Bool *RESTRICT connected)
Definition: shortestpath.c:147
SCIP_Bool graph_pc_isRootedPcMw(const GRAPH *)
#define Is_term(a)
Definition: graphdefs.h:102
SCIP_RETCODE shortestpath_pcInit(SCIP *scip, const GRAPH *graph, const SCIP_Real *costs, const SCIP_Real *prizes, SPATHSPC **sppc)
Definition: shortestpath.c:893
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
Shortest path based algorithms for Steiner problems.
int terms
Definition: graphdefs.h:192
int graph_pc_nFixedTerms(const GRAPH *)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
const CSR * csr
Definition: shortestpath.h:50
#define LT_HARD(a, b)
Definition: portab.h:82
#define FALSE
Definition: def.h:87
void shortestpath_computeSteinerTreePcMwFull(const GRAPH *g, int startnode, SPATHS *spaths)
int *RESTRICT inpbeg
Definition: graphdefs.h:202
SCIP_Real * cost
Definition: graphdefs.h:142
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
static void computeSteinerTree_execRpcMw(const GRAPH *g, int startnode, const SCIP_Real *prize, SPATHSPC *spaths_pc, SPATHS *spaths)
Definition: shortestpath.c:699
void graph_heap_clean(SCIP_Bool, DHEAP *)
Definition: graph_util.c:938
static void computeSteinerTree_connectPseudoTerm(const GRAPH *g, int k, const int *nodes_pred, SCIP_Real *RESTRICT nodes_dist, SPATHSPC *spaths_pc, DHEAP *dheap, STP_Bool *RESTRICT connected, int *RESTRICT nPseudoTerms)
Definition: shortestpath.c:511
SCIP_Real *RESTRICT nodes_dist
Definition: shortestpath.h:55
#define EAT_LAST
Definition: graphdefs.h:58
#define SCIPdebugMessage
Definition: pub_message.h:87
static void computeSteinerTree_init(const GRAPH *g, int startnode, SPATHS *spaths)
Definition: shortestpath.c:104
#define FARAWAY
Definition: graphdefs.h:89
static SCIP_Bool computeSteinerTree_allTermsAreReached(const GRAPH *g, const STP_Bool *connected)
Definition: shortestpath.c:35
void SCIPsortDownRealInt(SCIP_Real *realarray, int *intarray, int len)
void shortestpath_computeSteinerTreeDirected(SCIP *scip, const GRAPH *g, int startnode, SPATHS *spaths)
int * start
Definition: graphdefs.h:140
STP_Bool *RESTRICT nodes_isConnected
Definition: shortestpath.h:58
static SCIP_Bool computeSteinerTree_allFixedTermsAreReached(const GRAPH *g, const STP_Bool *connected)
Definition: shortestpath.c:80
header only, simple implementation of an STL like vector
static void computeSteinerTree_connectNode(const GRAPH *g, int k, const int *nodes_pred, SCIP_Real *RESTRICT nodes_dist, DHEAP *dheap, int *termscount, STP_Bool *RESTRICT connected)
Definition: shortestpath.c:182
static void computeSteinerTree_execDirected(SCIP *scip, const GRAPH *g, int startnode, SPATHS *spaths)
Definition: shortestpath.c:304
int * position
Definition: graphdefs.h:305
#define GE(a, b)
Definition: portab.h:85
SCIP_Bool extended
Definition: graphdefs.h:267
int stp_type
Definition: graphdefs.h:264
void shortestpath_pcFree(SCIP *scip, SPATHSPC **sppc)
Definition: shortestpath.c:964
void shortestpath_pcReset(SPATHSPC *sppc)
Definition: shortestpath.c:977
#define NULL
Definition: lpi_spx1.cpp:155
#define STP_RMWCSP
Definition: graphdefs.h:54
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
void shortestpath_computeSteinerTreeBiased(const GRAPH *g, const SDPROFIT *sdprofit, int startnode, SPATHS *spaths)
static SCIP_Bool computeSteinerTree_allPseudoTermsAreReached(const GRAPH *g, const STP_Bool *connected)
Definition: shortestpath.c:57
#define LT(a, b)
Definition: portab.h:81
unsigned char STP_Bool
Definition: portab.h:34
#define GT(a, b)
Definition: portab.h:84
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT ieat
Definition: graphdefs.h:230
static void computeSteinerTree_execPcMwFull(const GRAPH *g, int startnode, SPATHS *spaths)
Definition: shortestpath.c:809
static void computeSteinerTree_exec(const GRAPH *g, int startnode, SPATHS *spaths)
Definition: shortestpath.c:224
int *RESTRICT term
Definition: graphdefs.h:196
SCIP_Bool graph_pc_knotIsNonLeafTerm(const GRAPH *, int)
SCIP_Bool graph_pc_isPc(const GRAPH *)
#define CONNECT
Definition: graphdefs.h:87
Portable definitions.
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
void shortestpath_computeSteinerTreeRpcMw(const GRAPH *g, int startnode, const SCIP_Real *prize, SPATHSPC *spaths_pc, SPATHS *spaths)
#define Is_pseudoTerm(a)
Definition: graphdefs.h:103
SCIP_Bool graph_pc_knotIsDummyTerm(const GRAPH *, int)
int *RESTRICT nodes_pred
Definition: shortestpath.h:56
SCIP_Real *RESTRICT nodes_bias
Definition: reducedefs.h:137
SCIP_Bool graph_pc_knotIsFixedTerm(const GRAPH *, int)
SCIP_Bool graph_pc_isPcMw(const GRAPH *)
#define SCIP_Real
Definition: def.h:177
int graph_heap_deleteMinReturnNode(DHEAP *)
Definition: graph_util.c:1076
void shortestpath_computeSteinerTree(const GRAPH *g, int startnode, SPATHS *spaths)
static void computeSteinerTree_execBiased(const GRAPH *g, const SDPROFIT *sdprofit, int startnode, SPATHS *spaths)
Definition: shortestpath.c:408
#define StpVecFree(scip, vec)
Definition: stpvector.h:153
int edges
Definition: graphdefs.h:219
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
int graph_pc_nNonFixedTerms(const GRAPH *)
void shortestpath_pcConnectNode(const GRAPH *g, const STP_Bool *connected, int node_connected, SPATHSPC *sppc)
Definition: shortestpath.c:991
#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
#define nnodes
Definition: gastrans.c:65
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:167
#define Is_anyTerm(a)
Definition: graphdefs.h:105
static void computeSteinerTree_tryConnectNodePcMw(const GRAPH *g, int k, const SCIP_Real *prize, SCIP_Bool costIsBiased, const int *nodes_pred, SCIP_Real *RESTRICT nodes_dist, DHEAP *dheap, STP_Bool *RESTRICT connected, SPATHSPC *spaths_pc, SCIP_Bool *isConnected, int *RESTRICT nPseudoTerms)
Definition: shortestpath.c:556
SCIP_Bool graph_typeIsSpgLike(const GRAPH *)
Definition: graph_stats.c:57
static void computeSteinerTree_execPcMw(const GRAPH *g, int startnode, const SCIP_Real *prize, SCIP_Bool costIsBiased, SPATHSPC *spaths_pc, SPATHS *spaths)
Definition: shortestpath.c:607
includes various reduction methods for Steiner tree problems
#define STP_MWCSP
Definition: graphdefs.h:51
void graph_heap_correct(int, SCIP_Real, DHEAP *)
Definition: graph_util.c:1166