Scippy

SCIP

Solving Constraint Integer Programs

grphpath.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file grphpath.c
17  * @brief Shortest path based graph algorithms for Steiner problems
18  * @author Thorsten Koch
19  * @author Daniel Rehfeldt
20  *
21  * This file encompasses various (heap-based) shortest path based algorithms including
22  * Dijkstra's algorithm and Voronoi diagram algorithms
23  *
24  * The underlying heap routines can be found in Jon Bentley, Programming Pearls, Addison-Wesley 1989
25  *
26  * The heap array is initialized with n elements (nodes), but only at most n-1 nodes can be included
27  * in the array, since element 0 is not used for storing.
28  */
29 #include <stdlib.h>
30 #include <stdio.h>
31 #include <stddef.h>
32 #include <assert.h>
33 #include "portab.h"
34 #include "grph.h"
35 
36 /*---------------------------------------------------------------------------*/
37 /*--- Name : get NEAREST knot ---*/
38 /*--- Function : Holt das oberste Element vom Heap (den Knoten mit der ---*/
39 /*--- geringsten Entfernung zu den bereits Verbundenen) ---*/
40 /*--- Parameter: Derzeitige Entfernungen und benutzte Kanten ---*/
41 /*--- Returns : Nummer des bewussten Knotens ---*/
42 /*---------------------------------------------------------------------------*/
43 inline static int nearest(
44  int* heap,
45  int* state,
46  int* count, /* pointer to store the number of elements on the heap */
47  const PATH* path)
48 {
49  int k;
50  int t;
51  int c;
52  int j;
53 
54  /* Heap shift down
55  * (Oberstes Element runter und korrigieren)
56  */
57  k = heap[1];
58  j = 1;
59  c = 2;
60  heap[1] = heap[(*count)--];
61  state[heap[1]] = 1;
62 
63  if ((*count) > 2)
64  if (LT(path[heap[3]].dist, path[heap[2]].dist))
65  c++;
66 
67  while((c <= (*count)) && GT(path[heap[j]].dist, path[heap[c]].dist))
68  {
69  t = heap[c];
70  heap[c] = heap[j];
71  heap[j] = t;
72  state[heap[j]] = j;
73  state[heap[c]] = c;
74  j = c;
75  c += c;
76 
77  if ((c + 1) <= (*count))
78  if (LT(path[heap[c + 1]].dist, path[heap[c]].dist))
79  c++;
80  }
81  return(k);
82 }
83 
84 
85 /*---------------------------------------------------------------------------*/
86 /*--- Name : get NEAREST knot ---*/
87 /*--- Function : Holt das oberste Element vom Heap (den Knoten mit der ---*/
88 /*--- geringsten Entfernung zu den bereits Verbundenen) ---*/
89 /*--- Parameter: Derzeitige Entfernungen und benutzte Kanten ---*/
90 /*--- Returns : Nummer des bewussten Knotens ---*/
91 /*---------------------------------------------------------------------------*/
92 inline static int nearestX(
93  int* heap,
94  int* state,
95  int* count, /* pointer to store the number of elements on the heap */
96  const SCIP_Real* pathdist)
97 {
98  int k;
99  int t;
100  int c;
101  int j;
102  int dcount;
103 
104  /* Heap shift down
105  * (Oberstes Element runter und korrigieren)
106  */
107  k = heap[1];
108  j = 1;
109  c = 2;
110  heap[1] = heap[(*count)--];
111  state[heap[1]] = 1;
112 
113  dcount = *count;
114 
115  if (dcount > 2)
116  if (LT(pathdist[heap[3]], pathdist[heap[2]]))
117  c++;
118 
119  while((c <= dcount) && GT(pathdist[heap[j]], pathdist[heap[c]]))
120  {
121  t = heap[c];
122  heap[c] = heap[j];
123  heap[j] = t;
124  state[heap[j]] = j;
125  state[heap[c]] = c;
126  j = c;
127  c += c;
128 
129  if ((c + 1) <= dcount)
130  if (LT(pathdist[heap[c + 1]], pathdist[heap[c]]))
131  c++;
132  }
133  return(k);
134 }
135 
136 
137 /*---------------------------------------------------------------------------*/
138 /*--- Name : CORRECT heap ---*/
139 /*--- Function : Setzt ein neues Element auf den Heap, bzw. korrigiert ---*/
140 /*--- die Position eines vorhandenen Elementes ---*/
141 /*--- Parameter: Derzeitige Entfernungen und benutzten Kanten, ---*/
142 /*--- Neuer Knoten, Vorgaengerknoten, Kante von der man aus ---*/
143 /*--- den neuen Knoten erreicht, Kosten der Kante, ---*/
144 /*--- sowie Betriebsmodus ---*/
145 /*--- Returns : Nichts ---*/
146 /*---------------------------------------------------------------------------*/
147 inline static void correct(
148  SCIP* scip,
149  int* heap,
150  int* state,
151  int* count, /* pointer to store the number of elements on the heap */
152  PATH* path,
153  int l,
154  int k,
155  int e,
156  SCIP_Real cost,
157  int mode)
158 {
159  int t;
160  int c;
161  int j;
162 
163  path[l].dist = (mode == MST_MODE) ? cost : (path[k].dist + cost);
164  path[l].edge = e;
165 
166  /* new node? */
167  if( state[l] == UNKNOWN )
168  {
169  heap[++(*count)] = l;
170  state[l] = (*count);
171  }
172 
173  /* Heap shift up */
174  j = state[l];
175  c = j / 2;
176  while( (j > 1) && SCIPisGT(scip, path[heap[c]].dist, path[heap[j]].dist) )
177  {
178  t = heap[c];
179  heap[c] = heap[j];
180  heap[j] = t;
181  state[heap[j]] = j;
182  state[heap[c]] = c;
183  j = c;
184  c = j / 2;
185  }
186 }
187 
188 
189 /*---------------------------------------------------------------------------*/
190 /*--- Name : CORRECT heap ---*/
191 /*--- Function : Setzt ein neues Element auf den Heap, bzw. korrigiert ---*/
192 /*--- die Position eines vorhandenen Elementes ---*/
193 /*--- Parameter: Derzeitige Entfernungen und benutzten Kanten, ---*/
194 /*--- Neuer Knoten, Vorgaengerknoten, Kante von der man aus ---*/
195 /*--- den neuen Knoten erreicht, Kosten der Kante, ---*/
196 /*--- sowie Betriebsmodus ---*/
197 /*--- Returns : Nichts ---*/
198 /*---------------------------------------------------------------------------*/
199 inline static void correctX(
200  SCIP* scip,
201  int* heap,
202  int* state,
203  int* count, /* pointer to store the number of elements on the heap */
204  SCIP_Real* pathdist,
205  int* pathedge,
206  int l,
207  int k,
208  int e,
209  SCIP_Real cost
210  )
211 {
212  int t;
213  int c;
214  int j;
215 
216  pathdist[l] = (pathdist[k] + cost);
217  pathedge[l] = e;
218 
219  /* Ist der Knoten noch ganz frisch ?
220  */
221  if (state[l] == UNKNOWN)
222  {
223  heap[++(*count)] = l;
224  state[l] = (*count);
225  }
226 
227  /* Heap shift up
228  */
229  j = state[l];
230  c = j / 2;
231 
232  while( (j > 1) && SCIPisGT(scip, pathdist[heap[c]], pathdist[heap[j]]) )
233  {
234  t = heap[c];
235  heap[c] = heap[j];
236  heap[j] = t;
237  state[heap[j]] = j;
238  state[heap[c]] = c;
239  j = c;
240  c = j / 2;
241  }
242 }
243 
244 void heap_add(
245  int* heap, /* heaparray */
246  int* state,
247  int* count, /* pointer to store the number of elements on the heap */
248  int node, /* the node to be added */
249  PATH* path
250  )
251 {
252  int t;
253  int c;
254  int j;
255 
256  heap[++(*count)] = node;
257  state[node] = (*count);
258 
259  /* Heap shift up */
260  j = state[node];
261  c = j / 2;
262 
263  while((j > 1) && GT(path[heap[c]].dist, path[heap[j]].dist))
264  {
265  t = heap[c];
266  heap[c] = heap[j];
267  heap[j] = t;
268  state[heap[j]] = j;
269  state[heap[c]] = c;
270  j = c;
271  c = j / 2;
272  }
273 
274 }
275 inline static void resetX(
276  SCIP* scip,
277  SCIP_Real* pathdist,
278  int* heap,
279  int* state,
280  int* count,
281  int node
282  )
283 {
284  int t;
285  int c;
286  int j;
287 
288  pathdist[node] = 0.0;
289 
290  heap[++(*count)] = node;
291  state[node] = (*count);
292 
293  /* heap shift up */
294  j = state[node];
295  c = j / 2;
296 
297  while( (j > 1) && SCIPisGT(scip, pathdist[heap[c]], pathdist[heap[j]]) )
298  {
299  t = heap[c];
300  heap[c] = heap[j];
301  heap[j] = t;
302  state[heap[j]] = j;
303  state[heap[c]] = c;
304  j = c;
305  c = j / 2;
306  }
307 }
308 
309 inline static void reset(
310  SCIP* scip,
311  PATH* path,
312  int* heap,
313  int* state,
314  int* count,
315  int node
316  )
317 {
318  int t;
319  int c;
320  int j;
321 
322  path[node].dist = 0.0;
323 
324  heap[++(*count)] = node;
325  state[node] = (*count);
326 
327  /* heap shift up */
328  j = state[node];
329  c = j / 2;
330 
331  while( (j > 1) && SCIPisGT(scip, path[heap[c]].dist, path[heap[j]].dist) )
332  {
333  t = heap[c];
334  heap[c] = heap[j];
335  heap[j] = t;
336  state[heap[j]] = j;
337  state[heap[c]] = c;
338  j = c;
339  c = j / 2;
340  }
341 }
342 
343 inline static void utdist(
344  SCIP* scip,
345  const GRAPH* g,
346  PATH* path,
347  SCIP_Real ecost,
348  int* vbase,
349  int k,
350  int l,
351  int k2,
352  int shift,
353  int nnodes
354  )
355 {
356  SCIP_Real dist;
357  int vbk;
358  int vbk2;
359 
360  if( Is_term(g->term[k]) )
361  vbk = k;
362  else
363  vbk = vbase[k];
364 
365  if( l == 0 )
366  {
367  assert(shift == 0);
368 
369  dist = ecost;
370  if( !Is_term(g->term[k]) )
371  dist += path[k].dist;
372 
373  if( !Is_term(g->term[k2]) )
374  {
375  dist += path[k2].dist;
376  vbk2 = vbase[k2];
377  }
378  else
379  {
380  vbk2 = k2;
381  }
382 
383  if( SCIPisLT(scip, dist, path[vbk].dist) )
384  {
385  path[vbk].dist = dist;
386  vbase[vbk] = vbk2;
387  return;
388  }
389  }
390  else
391  {
392  int max;
393  int pos;
394  int r;
395  int s;
396  int t;
397 
398  pos = vbk + shift;
399  max = MIN((l + 1), 3);
400 
401  for( r = 0; r <= max; r++ )
402  {
403  if( Is_term(g->term[k2]) )
404  {
405  if( r == 0 )
406  t = k2;
407  else
408  break;
409  }
410  else
411  {
412  t = vbase[k2 + (r * nnodes)];
413  }
414  for( s = 0; s < l; s++ )
415  if( vbase[vbk + s * nnodes] == t )
416  break;
417  if( s < l || vbk == t )
418  continue;
419 
420  dist = ecost;
421  if( !Is_term(g->term[k]) )
422  dist += path[k].dist;
423  if( !Is_term(g->term[k2]) )
424  dist += path[k2 + (r * nnodes)].dist;
425 
426  if( SCIPisLT(scip, dist, path[pos].dist) )
427  {
428  path[pos].dist = dist;
429  vbase[pos] = t;
430  return;
431  }
432  }
433  }
434  return;
435 }
436 
437 /*---------------------------------------------------------------------------*/
438 /*--- Name : INIT shortest PATH algorithm ---*/
439 /*--- Function : Initialisiert den benoetigten Speicher fuer die ---*/
440 /*--- kuerzeste Wege Berechnung ---*/
441 /*--- Parameter: Graph in dem berechnet werden soll. ---*/
442 /*--- Returns : Nichts ---*/
443 /*---------------------------------------------------------------------------*/
444 SCIP_RETCODE graph_path_init(
445  SCIP* scip, /**< SCIP data structure */
446  GRAPH* g /**< graph data structure */
447  )
448 {
449  assert(g != NULL);
450  assert(g->path_heap == NULL);
451  assert(g->path_state == NULL);
452 
453  SCIP_CALL( SCIPallocMemoryArray(scip, &(g->path_heap), g->knots + 1) );
454  SCIP_CALL( SCIPallocMemoryArray(scip, &(g->path_state), g->knots) );
455 
456  return SCIP_OKAY;
457 }
458 
459 /*---------------------------------------------------------------------------*/
460 /*--- Name : EXIT shortest PATH algorithm ---*/
461 /*--- Function : Gibt den bei der initialisierung angeforderten Speicher ---*/
462 /*--- wieder frei ---*/
463 /*--- Parameter: Keine ---*/
464 /*--- Returns : Nichts ---*/
465 /*---------------------------------------------------------------------------*/
466 void graph_path_exit(
467  SCIP* scip, /**< SCIP data structure */
468  GRAPH* g /**< graph data structure */
469  )
470 {
471  assert(g != NULL);
472  assert(g->path_heap != NULL);
473  assert(g->path_state != NULL);
474 
475  SCIPfreeMemoryArray(scip, &(g->path_state));
476  SCIPfreeMemoryArray(scip, &(g->path_heap));
477 }
478 
479 /*---------------------------------------------------------------------------*/
480 /*--- Name : FIND shortest PATH / minimum spanning tree ---*/
481 /*--- Function : Dijkstras Algorithmus zur bestimmung eines kuerzesten ---*/
482 /*--- Weges in einem gerichteten Graphen. ---*/
483 /*--- Parameter: Graph, Nummer des Startknotens und Betriebsmodus ---*/
484 /*--- (Kuerzeste Wege oder minimal spannender Baum) ---*/
485 /*--- Returns : Liefert einen Vektor mit einem PATH Element je Knotem. ---*/
486 /*---- Setzt das .dist Feld zu jedem Knoten auf die Entfernung ---*/
487 /*--- zum Startknoten, sowie das .prev Feld auf den Knoten von ---*/
488 /*--- dem man zum Knoten gekommen ist. Im MST Modus steht im ---*/
489 /*--- .dist Feld die Entfernung zum naechsten Knoten. ---*/
490 /*---------------------------------------------------------------------------*/
491 void graph_path_exec(
492  SCIP* scip, /**< SCIP data structure */
493  const GRAPH* p, /**< graph data structure */
494  const int mode, /**< shortest path (FSP_MODE) or minimum spanning tree (MST_MODE)? */
495  int start, /**< start vertex */
496  const SCIP_Real* cost, /**< edge costs */
497  PATH* path /**< shortest paths data structure */
498  )
499 {
500  int* heap;
501  int* state;
502  int k;
503  const int nnodes = p->knots;
504  int count;
505 
506  assert(scip != NULL);
507  assert(p != NULL);
508  assert(start >= 0);
509  assert(start < p->knots);
510  assert((mode == FSP_MODE) || (mode == MST_MODE));
511  assert(p->path_heap != NULL);
512  assert(p->path_state != NULL);
513  assert(path != NULL);
514  assert(cost != NULL);
515 
516  /* no nodes?, return*/
517  if( nnodes == 0 )
518  return;
519 
520  heap = p->path_heap;
521  state = p->path_state;
522 
523  /* initialize */
524  for( int i = 0; i < nnodes; i++ )
525  {
526  state[i] = UNKNOWN;
527  path[i].dist = FARAWAY + 1;
528  path[i].edge = UNKNOWN;
529  }
530  /* add first node to heap */
531  k = start;
532  path[k].dist = 0.0;
533 
534  if( nnodes > 1 )
535  {
536  count = 1;
537  heap[count] = k;
538  state[k] = count;
539 
540  while( count > 0 )
541  {
542  /* get nearest labeled node */
543  k = nearest(heap, state, &count, path);
544 
545  /* mark as scanned */
546  state[k] = CONNECT;
547 
548  for( int i = p->outbeg[k]; i >= 0; i = p->oeat[i] )
549  {
550  const int m = p->head[i];
551 
552  assert(i != EAT_LAST);
553 
554  /* node not scanned and valid? */
555  if( state[m] )
556  {
557  /* closer than previously and valid? */
558  if( path[m].dist > ((mode == MST_MODE) ? cost[i] : (path[k].dist + cost[i])) && p->mark[m] )
559  correct(scip, heap, state, &count, path, m, k, i, cost[i], mode);
560  }
561  }
562  }
563  }
564 }
565 
566 /** limited Dijkstra, stopping at terminals */
567 void graph_sdPaths(
568  SCIP* scip, /**< SCIP data structure */
569  const GRAPH* g, /**< graph data structure */
570  PATH* path, /**< shortest paths data structure */
571  SCIP_Real* cost, /**< edge costs */
572  SCIP_Real distlimit, /**< distance limit of the search */
573  int* heap, /**< array representing a heap */
574  int* state, /**< array to indicate whether a node has been scanned during SP calculation */
575  int* memlbl, /**< array to save labelled nodes */
576  int* nlbl, /**< number of labelled nodes */
577  int tail, /**< tail of the edge */
578  int head, /**< head of the edge */
579  int limit /**< maximum number of edges to consider during execution */
580  )
581 {
582 
583  int count;
584  int nchecks;
585  const int limit1 = limit / 2;
586 
587  assert(g != NULL);
588  assert(heap != NULL);
589  assert(path != NULL);
590  assert(cost != NULL);
591  assert(nlbl != NULL);
592  assert(memlbl != NULL);
593  assert(limit1 >= 0);
594 
595  *nlbl = 0;
596 
597  if( g->grad[tail] == 0 || g->grad[head] == 0 )
598  return;
599 
600  assert(g->mark[head] && g->mark[tail]);
601 
602  count = 0;
603  nchecks = 0;
604  path[tail].dist = 0.0;
605  state[tail] = CONNECT;
606  memlbl[(*nlbl)++] = tail;
607 
608  if( g->stp_type != STP_MWCSP )
609  g->mark[head] = FALSE;
610 
611  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
612  {
613  const int m = g->head[e];
614 
615  if( g->mark[m] && (distlimit >= cost[e]) )
616  {
617  assert(SCIPisGT(scip, path[m].dist, path[tail].dist + cost[e]));
618 
619  /* m labelled the first time */
620  memlbl[(*nlbl)++] = m;
621  correct(scip, heap, state, &count, path, m, tail, e, cost[e], FSP_MODE);
622 
623  if( nchecks++ > limit1 )
624  break;
625  }
626  }
627  g->mark[head] = TRUE;
628 
629  while( count > 0 && nchecks <= limit )
630  {
631  /* get nearest labelled node */
632  const int k = nearest(heap, state, &count, path);
633 
634  /* scanned */
635  state[k] = CONNECT;
636 
637  /* distance limit reached? */
638  if( SCIPisGT(scip, path[k].dist, distlimit) )
639  break;
640 
641  /* stop at terminals */
642  if( Is_term(g->term[k]) || k == head )
643  continue;
644 
645  /* correct incident nodes */
646  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
647  {
648  const int m = g->head[e];
649  if( state[m] && g->mark[m] && (distlimit >= cost[e]) && (path[m].dist > path[k].dist + cost[e]) )
650  {
651  /* m labelled for the first time? */
652  if( state[m] == UNKNOWN )
653  memlbl[(*nlbl)++] = m;
654  correct(scip, heap, state, &count, path, m, k, e, cost[e], FSP_MODE);
655  }
656  if( nchecks++ > limit )
657  break;
658  }
659  }
660 }
661 
662 
663 /** limited Dijkstra for PCSTP, taking terminals into account */
664 void graph_path_PcMwSd(
665  SCIP* scip, /**< SCIP data structure */
666  const GRAPH* g, /**< graph data structure */
667  PATH* path, /**< shortest paths data structure */
668  SCIP_Real* cost, /**< edge costs */
669  SCIP_Real distlimit, /**< distance limit of the search */
670  int* pathmaxnode, /**< maximum weight node on path */
671  int* heap, /**< array representing a heap */
672  int* state, /**< array to indicate whether a node has been scanned during SP calculation */
673  int* stateblock, /**< array to indicate whether a node has been scanned during previous SP calculation */
674  int* memlbl, /**< array to save labelled nodes */
675  int* nlbl, /**< number of labelled nodes */
676  int tail, /**< tail of the edge */
677  int head, /**< head of the edge */
678  int limit /**< maximum number of edges to consider during execution */
679  )
680 {
681  const int limit1 = limit / 2;
682  int count;
683  int nchecks;
684  const SCIP_Bool block = (stateblock != NULL);
685 
686  assert(limit > 0);
687  assert(g != NULL);
688  assert(heap != NULL);
689  assert(path != NULL);
690  assert(cost != NULL);
691  assert(nlbl != NULL);
692  assert(memlbl != NULL);
693  assert(g->prize != NULL);
694  assert(pathmaxnode != NULL);
695 
696  *nlbl = 0;
697 
698  if( g->grad[tail] == 0 || g->grad[head] == 0 )
699  return;
700 
701  assert(g->mark[head] && g->mark[tail]);
702 
703  nchecks = 0;
704  count = 0;
705  path[tail].dist = 0.0;
706  state[tail] = CONNECT;
707  memlbl[(*nlbl)++] = tail;
708 
709  if( g->stp_type != STP_MWCSP )
710  g->mark[head] = FALSE;
711 
712  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
713  {
714  const int m = g->head[e];
715 
716  if( g->mark[m] && SCIPisLE(scip, cost[e], distlimit) )
717  {
718  assert(SCIPisGT(scip, path[m].dist, path[tail].dist + cost[e]));
719 
720  /* m labelled the first time */
721  memlbl[(*nlbl)++] = m;
722  correct(scip, heap, state, &count, path, m, tail, e, cost[e], FSP_MODE);
723 
724  if( nchecks++ > limit1 )
725  break;
726  }
727  }
728 
729  g->mark[head] = TRUE;
730 
731  /* main loop */
732  while( count > 0 )
733  {
734  const int k = nearest(heap, state, &count, path);
735  SCIP_Real maxweight = pathmaxnode[k] >= 0 ? g->prize[pathmaxnode[k]] : 0.0;
736 
737  assert(k != tail);
738  assert(maxweight >= 0);
739  assert(SCIPisLE(scip, path[k].dist - maxweight, distlimit));
740 
741  /* scanned */
742  state[k] = CONNECT;
743 
744  /* stop at other end */
745  if( k == head )
746  continue;
747 
748  if( Is_term(g->term[k]) && g->prize[k] > maxweight && distlimit >= path[k].dist )
749  {
750  pathmaxnode[k] = k;
751  maxweight = g->prize[k];
752  }
753 
754  /* stop at node scanned in first run */
755  if( block && stateblock[k] == CONNECT )
756  continue;
757 
758  /* correct incident nodes */
759  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
760  {
761  const int m = g->head[e];
762 
763  if( state[m] && g->mark[m] && path[m].dist > (path[k].dist + cost[e])
764  && distlimit >= (path[k].dist + cost[e] - maxweight) )
765  {
766  if( state[m] == UNKNOWN ) /* m labeled for the first time? */
767  memlbl[(*nlbl)++] = m;
768 
769  pathmaxnode[m] = pathmaxnode[k];
770  correct(scip, heap, state, &count, path, m, k, e, cost[e], FSP_MODE);
771  }
772  if( nchecks++ > limit )
773  break;
774  }
775  }
776 }
777 
778 
779 
780 /** Dijkstra's algorithm starting from node 'start' */
781 void graph_path_execX(
782  SCIP* scip, /**< SCIP data structure */
783  const GRAPH* g, /**< graph data structure */
784  int start, /**< start vertex */
785  const SCIP_Real* cost, /**< edge costs */
786  SCIP_Real* pathdist, /**< distance array (on vertices) */
787  int* pathedge /**< predecessor edge array (on vertices) */
788  )
789 {
790  int k;
791  int m;
792  int i;
793  int count;
794  int nnodes;
795  int* heap;
796  int* state;
797 
798  assert(g != NULL);
799  assert(start >= 0);
800  assert(start < g->knots);
801  assert(g->path_heap != NULL);
802  assert(g->path_state != NULL);
803  assert(pathdist != NULL);
804  assert(pathedge != NULL);
805  assert(cost != NULL);
806 
807  nnodes = g->knots;
808 
809  if( nnodes == 0 )
810  return;
811 
812  heap = g->path_heap;
813  state = g->path_state;
814 
815  for( i = nnodes - 1; i >= 0; --i )
816  {
817  state[i] = UNKNOWN;
818  pathdist[i] = FARAWAY;
819  pathedge[i] = -1;
820  }
821 
822  k = start;
823  pathdist[k] = 0.0;
824 
825  if( nnodes > 1 )
826  {
827  count = 1;
828  heap[count] = k;
829  state[k] = count;
830 
831  while( count > 0 )
832  {
833  k = nearestX(heap, state, &count, pathdist);
834 
835  state[k] = CONNECT;
836 
837  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
838  {
839  m = g->head[i];
840 
841  if( state[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[i])) )
842  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, i, cost[i]);
843  }
844  }
845  }
846 }
847 
848 /** Dijkstra on incoming edges until root is reached */
849 void graph_path_invroot(
850  SCIP* scip, /**< SCIP data structure */
851  const GRAPH* g, /**< graph data structure */
852  int start, /**< start vertex */
853  const SCIP_Real* cost, /**< edge costs */
854  SCIP_Real* pathdist, /**< distance array (on vertices) */
855  int* pathedge /**< predecessor edge array (on vertices) */
856  )
857 {
858  SCIP_Real rootdist;
859  int k;
860  int m;
861  int i;
862  int count;
863  int nnodes;
864  int* heap;
865  int* state;
866 
867  assert(g != NULL);
868  assert(start >= 0);
869  assert(start < g->knots);
870  assert(g->path_heap != NULL);
871  assert(g->path_state != NULL);
872  assert(pathdist != NULL);
873  assert(pathedge != NULL);
874  assert(cost != NULL);
875 
876  nnodes = g->knots;
877 
878  if( nnodes == 0 )
879  return;
880 
881  heap = g->path_heap;
882  state = g->path_state;
883  rootdist = FARAWAY;
884 
885  for( i = nnodes - 1; i >= 0; --i )
886  {
887  state[i] = UNKNOWN;
888  pathdist[i] = FARAWAY;
889  pathedge[i] = -1;
890  }
891 
892  k = start;
893  pathdist[k] = 0.0;
894 
895  if( nnodes > 1 )
896  {
897  int root = g->source;
898 
899  count = 1;
900  heap[count] = k;
901  state[k] = count;
902 
903  while( count > 0 )
904  {
905  k = nearestX(heap, state, &count, pathdist);
906 
907  state[k] = CONNECT;
908 
909  if( k == root )
910  rootdist = pathdist[k];
911  else if( SCIPisGT(scip, pathdist[k], rootdist) )
912  break;
913 
914  for( i = g->inpbeg[k]; i != EAT_LAST; i = g->ieat[i] )
915  {
916  m = g->tail[i];
917 
918  if( state[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[i])) )
919  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, i, cost[i]);
920  }
921  }
922  }
923 }
924 
925 /** Find a directed tree rooted in node 'start' and spanning all terminals */
926 void graph_path_st(
927  SCIP* scip, /**< SCIP data structure */
928  const GRAPH* g, /**< graph data structure */
929  SCIP_Real* cost, /**< edgecosts */
930  SCIP_Real* pathdist, /**< distance array (on vertices) */
931  int* pathedge, /**< predecessor edge array (on vertices) */
932  int start, /**< start vertex */
933  SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
934  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
935  )
936 {
937  int k;
938  int m;
939  int e;
940  int count;
941  int nnodes;
942  int* heap;
943  int* state;
944 
945  assert(randnumgen != NULL);
946  assert(pathdist != NULL);
947  assert(pathedge != NULL);
948  assert(g != NULL);
949  assert(start >= 0);
950  assert(start < g->knots);
951  assert(cost != NULL);
952  assert(connected != NULL);
953 
954  count = 0;
955  nnodes = g->knots;
956  heap = g->path_heap;
957  state = g->path_state;
958 
959  /* initialize */
960  for( k = 0; k < nnodes; k++ )
961  {
962  state[k] = UNKNOWN;
963  pathdist[k] = FARAWAY;
964  pathedge[k] = -1;
965  connected[k] = FALSE;
966  }
967 
968  /* add start vertex to heap */
969  k = start;
970  pathdist[k] = 0.0;
971  connected[k] = TRUE;
972 
973  if( nnodes > 1 )
974  {
975  int node;
976  int nterms = 0;
977 
978  count = 1;
979  heap[count] = k;
980  state[k] = count;
981 
982  /* repeat until heap is empty */
983  while( count > 0 )
984  {
985  /* get closest node */
986  k = nearestX(heap, state, &count, pathdist);
987  state[k] = UNKNOWN;
988 
989  /* if k is terminal, connect its path to current subtree */
990  if( Is_term(g->term[k]) )
991  {
992  assert(k == start || !connected[k]);
993  connected[k] = TRUE;
994  pathdist[k] = 0.0;
995  node = k;
996 
997  if( k != start )
998  {
999  assert(pathedge[k] != - 1);
1000 
1001  while( !connected[node = g->tail[pathedge[node]]] )
1002  {
1003  assert(pathedge[node] != - 1);
1004  connected[node] = TRUE;
1005  resetX(scip, pathdist, heap, state, &count, node);
1006  }
1007  }
1008  /* have all terminals been reached? */
1009  if( ++nterms == g->terms )
1010  break;
1011  }
1012 
1013  /* update adjacent vertices */
1014  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1015  {
1016  m = g->head[e];
1017 
1018  assert(state[m]);
1019 
1020  /* is m not connected, allowed and closer (as close)? */
1021  if( !connected[m] && pathdist[m] > (pathdist[k] + cost[e]) && g->mark[m] )
1022  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
1023  }
1024  }
1025  }
1026 }
1027 
1028 
1029 /** For rooted prize-collecting problem find a tree rooted in node 'start' and connecting
1030  * positive vertices as long as this is profitable.
1031  * Note that this function overwrites g->mark.
1032  * */
1033 void graph_path_st_rpc(
1034  SCIP* scip, /**< SCIP data structure */
1035  const GRAPH* g, /**< graph data structure */
1036  const SCIP_Real* cost, /**< edge costs */
1037  SCIP_Real* pathdist, /**< distance array (on vertices) */
1038  int* pathedge, /**< predecessor edge array (on vertices) */
1039  int start, /**< start vertex */
1040  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
1041  )
1042 {
1043  SCIP_Real maxprize;
1044  int k;
1045  int m;
1046  int e;
1047  int root;
1048  int count;
1049  int nnodes;
1050  int* heap;
1051  int* state;
1052 
1053  assert(pathdist != NULL);
1054  assert(pathedge != NULL);
1055  assert(g != NULL);
1056  assert(start >= 0);
1057  assert(start < g->knots);
1058  assert(cost != NULL);
1059  assert(connected != NULL);
1060  assert(g->stp_type == STP_RPCSPG);
1061 
1062  root = g->source;
1063  heap = g->path_heap;
1064  count = 0;
1065  state = g->path_state;
1066  nnodes = g->knots;
1067  maxprize = 0.0;
1068 
1069  /* initialize */
1070  for( k = nnodes - 1; k >= 0; --k )
1071  {
1072  g->mark[k] = ((g->grad[k] > 0) && !Is_term(g->term[k]));
1073  state[k] = UNKNOWN;
1074  pathdist[k] = FARAWAY;
1075  pathedge[k] = -1;
1076  connected[k] = FALSE;
1077 
1078  if( Is_pterm(g->term[k]) && SCIPisGT(scip, g->prize[k], maxprize) && k != start )
1079  {
1080  maxprize = g->prize[k];
1081  }
1082  }
1083 
1084  g->mark[root] = TRUE;
1085 
1086  /* add start vertex to heap */
1087  k = start;
1088  pathdist[k] = 0.0;
1089  connected[k] = TRUE;
1090 
1091  if( nnodes > 1 )
1092  {
1093  int node;
1094  int nterms = 0;
1095 
1096  count = 1;
1097  heap[count] = k;
1098  state[k] = count;
1099 
1100  /* repeat until heap is empty */
1101  while( count > 0 )
1102  {
1103  /* get closest node */
1104  k = nearestX(heap, state, &count, pathdist);
1105  state[k] = UNKNOWN;
1106 
1107  /* if k is positive vertex and close enough, connect its path to current subtree */
1108  if( ((Is_pterm(g->term[k]) && SCIPisGE(scip, g->prize[k], pathdist[k])) || k == root) && !connected[k] )
1109  {
1110  connected[k] = TRUE;
1111  pathdist[k] = 0.0;
1112  node = k;
1113 
1114  if( k != start )
1115  {
1116  assert(pathedge[k] != - 1);
1117 
1118  while( !connected[node = g->tail[pathedge[node]]] )
1119  {
1120  assert(pathedge[node] != - 1);
1121  connected[node] = TRUE;
1122  resetX(scip, pathdist, heap, state, &count, node);
1123  }
1124  }
1125  /* have all terminals been reached? */
1126  if( ++nterms == g->terms - 1 )
1127  break;
1128  }
1129  else if( SCIPisGT(scip, pathdist[k], maxprize) && connected[root] )
1130  {
1131  break;
1132  }
1133 
1134  /* update adjacent vertices */
1135  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1136  {
1137  m = g->head[e];
1138 
1139  assert(state[m]);
1140 
1141  /* is m not connected, allowed and closer (as close)? */
1142  if( !connected[m] && pathdist[m] > (pathdist[k] + cost[e]) && g->mark[m] )
1143  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
1144  }
1145  }
1146  }
1147 }
1148 
1149 
1150 /** Find a tree rooted in node 'start' and connecting
1151  * positive vertices as long as this is profitable.
1152  * Note that this function overwrites g->mark.
1153  * */
1154 void graph_path_st_pcmw(
1155  SCIP* scip, /**< SCIP data structure */
1156  const GRAPH* g, /**< graph data structure */
1157  const SCIP_Real* cost, /**< edge costs */
1158  SCIP_Real* pathdist, /**< distance array (on vertices) */
1159  int* pathedge, /**< predecessor edge array (on vertices) */
1160  int start, /**< start vertex */
1161  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
1162  )
1163 {
1164  SCIP_Real maxprize;
1165  int k;
1166  int count;
1167  const int nnodes = g->knots;
1168  int* const heap = g->path_heap;
1169  int* const state = g->path_state;
1170 
1171  assert(pathdist != NULL);
1172  assert(pathedge != NULL);
1173  assert(g != NULL);
1174  assert(start >= 0);
1175  assert(start < g->knots);
1176  assert(cost != NULL);
1177  assert(connected != NULL);
1178 
1179  maxprize = 0.0;
1180  count = 0;
1181 
1182  /* initialize */
1183  for( k = nnodes - 1; k >= 0; --k )
1184  {
1185  g->mark[k] = ((g->grad[k] > 0) && !Is_term(g->term[k]));
1186  state[k] = UNKNOWN;
1187  pathdist[k] = FARAWAY;
1188  pathedge[k] = -1;
1189  connected[k] = FALSE;
1190 
1191  if( Is_pterm(g->term[k]) && (g->prize[k] > maxprize) && k != start )
1192  maxprize = g->prize[k];
1193  }
1194 
1195  /* add start vertex to heap */
1196  k = start;
1197  pathdist[k] = 0.0;
1198  connected[k] = TRUE;
1199 
1200  if( nnodes > 1 )
1201  {
1202  int nterms = 0;
1203 
1204  count = 1;
1205  heap[count] = k;
1206  state[k] = count;
1207 
1208  /* repeat until heap is empty */
1209  while( count > 0 )
1210  {
1211  /* get closest node */
1212  k = nearestX(heap, state, &count, pathdist);
1213  state[k] = UNKNOWN;
1214 
1215  /* if k is positive vertex and close enough, connect its path to current subtree */
1216  if( !connected[k] && Is_pterm(g->term[k]) && (g->prize[k] >= pathdist[k]) )
1217  {
1218  int node = k;
1219 
1220  connected[k] = TRUE;
1221  pathdist[k] = 0.0;
1222 
1223  if( k != start )
1224  {
1225  assert(pathedge[k] != - 1);
1226 
1227  while( !connected[node = g->tail[pathedge[node]]] )
1228  {
1229  assert(pathedge[node] != - 1);
1230  connected[node] = TRUE;
1231  resetX(scip, pathdist, heap, state, &count, node);
1232  assert(state[node]);
1233  }
1234  }
1235  /* have all terminals been reached? */
1236  if( ++nterms == g->terms - 1 )
1237  break;
1238  }
1239  else if( SCIPisGT(scip, pathdist[k], maxprize) )
1240  {
1241  break;
1242  }
1243 
1244  /* update adjacent vertices */
1245  for( int e = g->outbeg[k]; e >= 0; e = g->oeat[e] )
1246  {
1247  const int m = g->head[e];
1248 
1249  assert(state[m] && e != EAT_LAST);
1250  /* is m not connected, allowed and closer (as close)? */
1251 
1252  if( !connected[m] && pathdist[m] > (pathdist[k] + cost[e]) && g->mark[m] )
1253  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
1254  }
1255  }
1256  }
1257 }
1258 
1259 #if 1
1260 
1261 /** Reduce given solution
1262  * Note that this function overwrites g->mark.
1263  * */
1264 void graph_path_st_pcmw_reduce(
1265  SCIP* scip, /**< SCIP data structure */
1266  const GRAPH* g, /**< graph data structure */
1267  const SCIP_Real* cost, /**< edge costs */
1268  SCIP_Real* tmpnodeweight, /**< node weight array */
1269  int* result, /**< incoming arc array */
1270  int start, /**< start vertex */
1271  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
1272  )
1273 {
1274  assert(tmpnodeweight != NULL);
1275  assert(result != NULL);
1276  assert(g != NULL);
1277  assert(start >= 0);
1278  assert(start < g->knots);
1279  assert(cost != NULL);
1280  assert(connected != NULL);
1281 
1282  for( int e = g->outbeg[start]; e != EAT_LAST; e = g->oeat[e] )
1283  {
1284  if( result[e] == CONNECT )
1285  {
1286  const int head = g->head[e];
1287 
1288  if( Is_term(g->term[head]) )
1289  continue;
1290 
1291  graph_path_st_pcmw_reduce(scip, g, cost, tmpnodeweight, result, head, connected);
1292 
1293  assert(connected[head]);
1294 
1295  if( SCIPisGE(scip, cost[e], tmpnodeweight[head]) )
1296  {
1297  connected[head] = FALSE;
1298  result[e] = UNKNOWN;
1299  }
1300  else
1301  {
1302  tmpnodeweight[start] += tmpnodeweight[head] - cost[e];
1303  }
1304  }
1305  }
1306 
1307  SCIPfreeBufferArray(scip, &tmpnodeweight);
1308 }
1309 #endif
1310 
1311 
1312 /** Find a tree rooted in node 'start' and connecting
1313  * all positive vertices.
1314  * Note that this function overwrites g->mark.
1315  * */
1316 void graph_path_st_pcmw_full(
1317  SCIP* scip, /**< SCIP data structure */
1318  const GRAPH* g, /**< graph data structure */
1319  const SCIP_Real* cost, /**< edge costs */
1320  SCIP_Real* pathdist, /**< distance array (on vertices) */
1321  int* pathedge, /**< predecessor edge array (on vertices) */
1322  int start, /**< start vertex */
1323  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
1324  )
1325 {
1326  int* const heap = g->path_heap;
1327  int* const state = g->path_state;
1328  int k;
1329  const int nnodes = g->knots;
1330 
1331  assert(pathdist != NULL);
1332  assert(pathedge != NULL);
1333  assert(g != NULL);
1334  assert(start >= 0);
1335  assert(start < g->knots);
1336  assert(cost != NULL);
1337  assert(connected != NULL);
1338 
1339  /* initialize */
1340  for( k = nnodes - 1; k >= 0; --k )
1341  {
1342  g->mark[k] = ((g->grad[k] > 0) && !Is_term(g->term[k]));
1343  state[k] = UNKNOWN;
1344  pathdist[k] = FARAWAY;
1345  pathedge[k] = -1;
1346  connected[k] = FALSE;
1347  }
1348 
1349  /* add start vertex to heap */
1350  k = start;
1351  pathdist[k] = 0.0;
1352  connected[k] = TRUE;
1353 
1354  if( nnodes > 1 )
1355  {
1356  int count = 1;
1357  int nterms = 0;
1358 
1359  heap[count] = k;
1360  state[k] = count;
1361 
1362  /* repeat until heap is empty */
1363  while( count > 0 )
1364  {
1365  /* get closest node */
1366  k = nearestX(heap, state, &count, pathdist);
1367  state[k] = UNKNOWN;
1368 
1369  /* if k is positive vertex, connect its path to current subtree */
1370  if( !connected[k] && Is_pterm(g->term[k]) )
1371  {
1372  connected[k] = TRUE;
1373  pathdist[k] = 0.0;
1374 
1375  if( k != start )
1376  {
1377  int node = k;
1378  assert(pathedge[k] != - 1);
1379 
1380  while( !connected[node = g->tail[pathedge[node]]] )
1381  {
1382  assert(pathedge[node] != - 1);
1383  connected[node] = TRUE;
1384  resetX(scip, pathdist, heap, state, &count, node);
1385  assert(state[node]);
1386  }
1387  }
1388  /* have all terminals been reached? */
1389  if( ++nterms == g->terms - 1 )
1390  break;
1391  }
1392 
1393  /* update adjacent vertices */
1394  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1395  {
1396  int m = g->head[e];
1397 
1398  assert(state[m]);
1399  /* is m not connected, allowed and closer (as close)? */
1400 
1401  if( !connected[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[e])) )
1402  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
1403  }
1404  }
1405  }
1406 }
1407 
1408 
1409 /** greedy extension of a given tree for PC or MW problems */
1410 void graph_path_st_pcmw_extend(
1411  SCIP* scip, /**< SCIP data structure */
1412  const GRAPH* g, /**< graph data structure */
1413  const SCIP_Real* cost, /**< edge costs */
1414  PATH* path, /**< shortest paths data structure */
1415  STP_Bool* connected, /**< array to mark whether a vertex is part of computed Steiner tree */
1416  SCIP_Bool* extensions /**< extensions performed? */
1417  )
1418 {
1419  SCIP_Real maxprize;
1420  int k;
1421  int count;
1422  int nnodes;
1423  int nstnodes;
1424  int* heap;
1425  int* state;
1426 
1427  assert(path != NULL);
1428  assert(g != NULL);
1429  assert(cost != NULL);
1430  assert(connected != NULL);
1431 
1432  maxprize = 0.0;
1433  count = 0;
1434  nstnodes = 0;
1435  nnodes = g->knots;
1436  heap = g->path_heap;
1437  state = g->path_state;
1438  *extensions = FALSE;
1439 
1440  /* initialize */
1441  for( k = 0; k < nnodes; k++ )
1442  {
1443  g->mark[k] = ((g->grad[k] > 0) && !Is_term(g->term[k]));
1444  if( connected[k] && g->mark[k] )
1445  {
1446  /* add node to heap */
1447  nstnodes++;
1448  if( nnodes > 1 )
1449  heap[++count] = k;
1450 
1451  state[k] = count;
1452  path[k].dist = 0.0;
1453  assert(path[k].edge != UNKNOWN || k == g->source);
1454  }
1455  else
1456  {
1457  state[k] = UNKNOWN;
1458  path[k].dist = FARAWAY;
1459 
1460  if( Is_pterm(g->term[k]) && SCIPisGT(scip, g->prize[k], maxprize) && g->mark[k] )
1461  {
1462  maxprize = g->prize[k];
1463  }
1464  }
1465 
1466  if( !connected[k] )
1467  path[k].edge = UNKNOWN;
1468  }
1469 
1470  /* nothing to extend? */
1471  if( nstnodes == 0 )
1472  return;
1473 
1474  if( nnodes > 1 )
1475  {
1476  int node;
1477  int nterms = 0;
1478 
1479  /* repeat until heap is empty */
1480  while( count > 0 )
1481  {
1482  /* get closest node */
1483  k = nearest(heap, state, &count, path);
1484  state[k] = UNKNOWN;
1485 
1486  /* if k is positive vertex and close enough (or fixnode), connect its path to current subtree */
1487  if( !connected[k] && Is_pterm(g->term[k]) &&
1488  SCIPisGE(scip, g->prize[k], path[k].dist) )
1489  {
1490  connected[k] = TRUE;
1491  *extensions = TRUE;
1492  path[k].dist = 0.0;
1493  node = k;
1494 
1495  assert(path[k].edge != UNKNOWN);
1496 
1497  while( !connected[node = g->tail[path[node].edge]] )
1498  {
1499  assert(path[node].edge != UNKNOWN);
1500  connected[node] = TRUE;
1501  reset(scip, path, heap, state, &count, node);
1502  assert(state[node]);
1503  }
1504 
1505  /* have all terminals been reached? */
1506  if( ++nterms == g->terms - 1 )
1507  break;
1508  }
1509  else if( SCIPisGT(scip, path[k].dist, maxprize) )
1510  {
1511  break;
1512  }
1513 
1514  /* update adjacent vertices */
1515  for( int e = g->outbeg[k]; e >= 0; e = g->oeat[e] )
1516  {
1517  const int m = g->head[e];
1518 
1519  assert(state[m]);
1520  assert(e != EAT_LAST);
1521 
1522  /* is m not connected, allowed and closer (as close)? */
1523 
1524  if( !connected[m] && path[m].dist > (path[k].dist + cost[e]) && g->mark[m] )
1525  correct(scip, heap, state, &count, path, m, k, e, cost[e], FSP_MODE);
1526  }
1527  }
1528  }
1529 }
1530 
1531 
1532 /** Shortest path heuristic for the RMWCSP
1533  * Find a directed tree rooted in node 'start' and connecting all terminals as well as all
1534  * positive vertices (as long as this is profitable).
1535  * */
1536 void graph_path_st_rmw(
1537  SCIP* scip, /**< SCIP data structure */
1538  const GRAPH* g, /**< graph data structure */
1539  const SCIP_Real* cost, /**< edge costs */
1540  SCIP_Real* pathdist, /**< distance array (on vertices) */
1541  int* pathedge, /**< predecessor edge array (on vertices) */
1542  int start, /**< start vertex */
1543  STP_Bool* connected /**< array to mark whether a vertex is part of computed Steiner tree */
1544  )
1545 {
1546  SCIP_Real maxprize;
1547  int k;
1548  int e;
1549  int root;
1550  int count;
1551  int nrterms;
1552  int nnodes;
1553  int* heap;
1554  int* state;
1555 
1556  assert(pathdist != NULL);
1557  assert(pathedge != NULL);
1558  assert(g != NULL);
1559  assert(start >= 0);
1560  assert(start < g->knots);
1561  assert(cost != NULL);
1562  assert(connected != NULL);
1563 
1564  root = g->source;
1565  count = 0;
1566  nrterms = 0;
1567  state = g->path_state;
1568  heap = g->path_heap;
1569  nnodes = g->knots;
1570  maxprize = 0.0;
1571 
1572  for( k = 0; k < nnodes; k++ )
1573  g->mark[k] = (g->grad[k] > 0);
1574 
1575  for( e = g->outbeg[root]; e != EAT_LAST; e = g->oeat[e] )
1576  {
1577  if( SCIPisGT(scip, g->cost[e], 0.0) && Is_term(g->term[g->head[e]]) )
1578  {
1579  k = g->head[e];
1580  g->mark[k] = FALSE;
1581  assert(g->grad[k] == 2);
1582  }
1583  }
1584 
1585  for( k = 0; k < nnodes; k++ )
1586  {
1587  state[k] = UNKNOWN;
1588  pathdist[k] = FARAWAY;
1589  pathedge[k] = -1;
1590  connected[k] = FALSE;
1591 
1592  if( Is_term(g->term[k]) && g->mark[k] )
1593  nrterms++;
1594 
1595  if( Is_pterm(g->term[k]) )
1596  if( SCIPisGT(scip, g->prize[k], maxprize) && k != start )
1597  maxprize = g->prize[k];
1598  }
1599 
1600  /* add start vertex to heap */
1601  k = start;
1602  pathdist[k] = 0.0;
1603  connected[k] = TRUE;
1604 
1605  if( nnodes > 1 )
1606  {
1607  int node;
1608  int nterms = g->terms;
1609  int termscount = 0;
1610  int rtermscount = 0;
1611 
1612  count = 1;
1613  heap[count] = k;
1614  state[k] = count;
1615 
1616  /* repeat until heap is empty */
1617  while( count > 0 )
1618  {
1619  /* get closest node */
1620  k = nearestX(heap, state, &count, pathdist);
1621  state[k] = UNKNOWN;
1622  /* if k is positive vertex and close enough, connect its path to current subtree */
1623  if( Is_gterm(g->term[k]) && (Is_term(g->term[k]) || SCIPisGE(scip, g->prize[k], pathdist[k])) && !connected[k] )
1624  {
1625  if( Is_term(g->term[k]) )
1626  rtermscount++;
1627  connected[k] = TRUE;
1628  pathdist[k] = 0.0;
1629  node = k;
1630 
1631  if( k != start )
1632  {
1633  assert(pathedge[k] != - 1);
1634  while( !connected[node = g->tail[pathedge[node]]] )
1635  {
1636  assert(pathedge[node] != - 1);
1637 
1638  connected[node] = TRUE;
1639  resetX(scip, pathdist, heap, state, &count, node);
1640  }
1641  }
1642 
1643  /* have all terminals been reached? */
1644  if( ++termscount == nterms )
1645  {
1646  break;
1647  }
1648  }
1649  else if( SCIPisGT(scip, pathdist[k], maxprize) && rtermscount >= nrterms )
1650  {
1651  break;
1652  }
1653 
1654  /* update adjacent vertices */
1655  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1656  {
1657  int m = g->head[e];
1658 
1659  assert(state[m]);
1660 
1661  /* is m not connected, allowed and closer (as close)? */
1662  if( !connected[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[e])) )
1663  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
1664  }
1665  }
1666  }
1667 }
1668 
1669 
1670 /** extend a voronoi region until all neighbouring terminals are spanned */
1671 SCIP_RETCODE graph_voronoiExtend(
1672  SCIP* scip, /**< SCIP data structure */
1673  const GRAPH* g, /**< graph data structure */
1674  SCIP_Real* cost, /**< edgecosts */
1675  PATH* path, /**< shortest paths data structure */
1676  SCIP_Real** distarr, /**< array to store distance from each node to its base */
1677  int** basearr, /**< array to store the bases */
1678  int** edgearr, /**< array to store the ancestor edge */
1679  STP_Bool* termsmark, /**< array to mark terminal */
1680  int* reachednodes, /**< array to mark reached nodes */
1681  int* nreachednodes, /**< pointer to number of reached nodes */
1682  int* nodenterms, /**< array to store number of terimals to each node */
1683  int nneighbterms, /**< number of neighbouring terminals */
1684  int base, /**< voronoi base */
1685  int countex /**< count of heap elements */
1686  )
1687 {
1688  int k;
1689  int m;
1690  int i;
1691  int* heap;
1692  int* state;
1693  int count = countex;
1694 
1695  assert(g != NULL);
1696  assert(g->path_heap != NULL);
1697  assert(g->path_state != NULL);
1698  assert(path != NULL);
1699  assert(cost != NULL);
1700 
1701  if( g->knots == 0 || nneighbterms <= 0 )
1702  return SCIP_OKAY;
1703 
1704  heap = g->path_heap;
1705  state = g->path_state;
1706 
1707  if( g->knots > 1 )
1708  {
1709  while( count > 0 && nneighbterms > 0 )
1710  {
1711  k = nearest(heap, state, &count, path);
1712  state[k] = CONNECT;
1713  reachednodes[(*nreachednodes)++] = k;
1714  distarr[k][nodenterms[k]] = path[k].dist;
1715  edgearr[k][nodenterms[k]] = path[k].edge;
1716  basearr[k][nodenterms[k]] = base;
1717 
1718  nodenterms[k]++;
1719 
1720  if( termsmark[k] == TRUE )
1721  {
1722  termsmark[k] = FALSE;
1723  if( --nneighbterms == 0 )
1724  {
1725  while( count > 0 )
1726  reachednodes[(*nreachednodes)++] = heap[count--];
1727 
1728  return SCIP_OKAY;
1729  }
1730  }
1731  else if( !Is_term(g->term[k]) )
1732  {
1733  /* iterate over all outgoing edges of vertex k */
1734  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1735  {
1736  m = g->head[i];
1737 
1738  /* check whether the path (to m) including (k, m) is shorter than the so far best known */
1739  if( (state[m]) && (g->mark[m]) && (GT(path[m].dist, (path[k].dist + cost[i]))) )
1740  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
1741  }
1742  }
1743  }
1744  assert(nneighbterms == 0);
1745  }
1746  return SCIP_OKAY;
1747 }
1748 
1749 
1750 /** build a voronoi region, w.r.t. shortest paths, for a given set of bases */
1751 void graph_voronoi(
1752  SCIP* scip, /**< SCIP data structure */
1753  const GRAPH* g, /**< graph data structure */
1754  SCIP_Real* cost, /**< edge costs */
1755  SCIP_Real* costrev, /**< reversed edge costs */
1756  STP_Bool* base, /**< array to indicate whether a vertex is a Voronoi base */
1757  int* vbase, /**< voronoi base to each vertex */
1758  PATH* path /**< path data struture (leading to respective Voronoi base) */
1759  )
1760 {
1761  int k;
1762  int m;
1763  int i;
1764  int* heap;
1765  int* state;
1766  int count = 0;
1767  int root;
1768  int nbases = 0;
1769 
1770  assert(g != NULL);
1771  assert(scip != NULL);
1772  assert(path != NULL);
1773  assert(cost != NULL);
1774  assert(costrev != NULL);
1775  assert(g->path_heap != NULL);
1776  assert(g->path_state != NULL);
1777 
1778  if( g->knots == 0 )
1779  return;
1780 
1781  root = g->source;
1782  heap = g->path_heap;
1783  state = g->path_state;
1784 
1785  /* initialize */
1786  for( i = 0; i < g->knots; i++ )
1787  {
1788  /* set the base of vertex i */
1789  if( base[i] )
1790  {
1791  nbases++;
1792  if( g->knots > 1 )
1793  heap[++count] = i;
1794  vbase[i] = i;
1795  path[i].dist = 0.0;
1796  path[i].edge = UNKNOWN;
1797  state[i] = count;
1798  }
1799  else
1800  {
1801  vbase[i] = UNKNOWN;
1802  path[i].dist = FARAWAY + 1;
1803  path[i].edge = UNKNOWN;
1804  state[i] = UNKNOWN;
1805  }
1806  }
1807  assert(nbases > 0);
1808 
1809  if( g->knots > 1 )
1810  {
1811  /* until the heap is empty */
1812  while( count > 0 )
1813  {
1814  /* get the next (i.e. a nearest) vertex of the heap */
1815  k = nearest(heap, state, &count, path);
1816 
1817  /* mark vertex k as scanned */
1818  state[k] = CONNECT;
1819 
1820  /* iterate over all outgoing edges of vertex k */
1821  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1822  {
1823  m = g->head[i];
1824 
1825  /* check whether the path (to m) including k is shorter than the so far best known */
1826  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + ((vbase[k] == root)? cost[i] : costrev[i])) )
1827  {
1828  assert(!base[m]);
1829  correct(scip, heap, state, &count, path, m, k, i, ((vbase[k] == root)? cost[i] : costrev[i]), FSP_MODE);
1830  vbase[m] = vbase[k];
1831  }
1832  }
1833  }
1834  }
1835 }
1836 
1837 /** 2th next terminal to all non terminal nodes */
1838 void graph_get2next(
1839  SCIP* scip, /**< SCIP data structure */
1840  const GRAPH* g, /**< graph data structure */
1841  const SCIP_Real* cost, /**< edge costs */
1842  const SCIP_Real* costrev, /**< reversed edge costs */
1843  PATH* path, /**< path data structure (leading to first and second nearest terminal) */
1844  int* vbase, /**< first and second nearest terminal to each non terminal */
1845  int* heap, /**< array representing a heap */
1846  int* state /**< array to mark the state of each node during calculation */
1847  )
1848 {
1849  int k;
1850  int j;
1851  int i;
1852  int e;
1853  int root;
1854  int count;
1855  int nnodes;
1856 
1857  assert(g != NULL);
1858  assert(path != NULL);
1859  assert(cost != NULL);
1860  assert(heap != NULL);
1861  assert(state != NULL);
1862  assert(costrev != NULL);
1863 
1864  root = g->source;
1865  count = 0;
1866  nnodes = g->knots;
1867 
1868  /* initialize */
1869  for( i = 0; i < nnodes; i++ )
1870  {
1871  state[i] = CONNECT;
1872 
1873  /* copy of node i */
1874  k = i + nnodes;
1875  vbase[k] = UNKNOWN;
1876  state[k] = UNKNOWN;
1877  path[k].edge = UNKNOWN;
1878  path[k].dist = FARAWAY;
1879  }
1880 
1881  /* scan original nodes */
1882  for( i = 0; i < nnodes; i++ )
1883  {
1884  if( !g->mark[i] )
1885  continue;
1886 
1887  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
1888  {
1889  j = g->head[e];
1890  k = j + nnodes;
1891  if( !Is_term(g->term[j]) && SCIPisGT(scip, path[k].dist, path[i].dist + ((root == vbase[i])? cost[e] : costrev[e])) &&
1892  vbase[i] != vbase[j] && g->mark[j] )
1893  {
1894  correct(scip, heap, state, &count, path, k, i, e, ((root == vbase[i])? cost[e] : costrev[e]), FSP_MODE);
1895  vbase[k] = vbase[i];
1896  }
1897  }
1898  }
1899 
1900  if( nnodes > 1 )
1901  {
1902  int jc;
1903  /* until the heap is empty */
1904  while( count > 0 )
1905  {
1906  /* get the next (i.e. a nearest) vertex of the heap */
1907  k = nearest(heap, state, &count, path);
1908 
1909  /* mark vertex k as removed from heap */
1910  state[k] = UNKNOWN;
1911 
1912  assert(k - nnodes >= 0);
1913  /* iterate over all outgoing edges of vertex k */
1914  for( e = g->outbeg[k - nnodes]; e != EAT_LAST; e = g->oeat[e] )
1915  {
1916  j = g->head[e];
1917  if( Is_term(g->term[j]) || !g->mark[j] )
1918  continue;
1919  jc = j + nnodes;
1920 
1921  /* check whether the path (to j) including k is shorter than the so far best known */
1922  if( vbase[j] != vbase[k] && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) )
1923  {
1924  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
1925  vbase[jc] = vbase[k];
1926  }
1927  }
1928  }
1929  }
1930  return;
1931 }
1932 
1933 /* 3th next terminal to all non terminal nodes */
1934 void graph_get3next(
1935  SCIP* scip, /**< SCIP data structure */
1936  const GRAPH* g, /**< graph data structure */
1937  const SCIP_Real* cost, /**< edge costs */
1938  const SCIP_Real* costrev, /**< reversed edge costs */
1939  PATH* path, /**< path data structure (leading to first, second and third nearest terminal) */
1940  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1941  int* heap, /**< array representing a heap */
1942  int* state /**< array to mark the state of each node during calculation */
1943  )
1944 {
1945  int k;
1946  int j;
1947  int i;
1948  int l;
1949  int v;
1950  int e;
1951  int root;
1952  int count = 0;
1953  int nnodes;
1954  int dnnodes;
1955 
1956  assert(g != NULL);
1957  assert(path != NULL);
1958  assert(cost != NULL);
1959  assert(heap != NULL);
1960  assert(state != NULL);
1961  assert(costrev != NULL);
1962 
1963  root = g->source;
1964  nnodes = g->knots;
1965  dnnodes = 2 * nnodes;
1966 
1967  /* initialize */
1968  for( i = 0; i < nnodes; i++ )
1969  {
1970  /* copy of node i */
1971  k = i + dnnodes;
1972  vbase[k] = UNKNOWN;
1973  state[k] = UNKNOWN;
1974  path[k].edge = UNKNOWN;
1975  path[k].dist = FARAWAY;
1976  }
1977 
1978  /* scan original nodes */
1979  for( i = 0; i < nnodes; i++ )
1980  {
1981  state[i] = CONNECT;
1982  state[i + nnodes] = CONNECT;
1983  if( !g->mark[i] )
1984  continue;
1985 
1986  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
1987  {
1988  j = g->head[e];
1989  k = j + dnnodes;
1990  if( !Is_term(g->term[j]) && g->mark[j] )
1991  {
1992  v = i;
1993  for( l = 0; l < 2; l++ )
1994  {
1995  if( SCIPisGT(scip, path[k].dist, path[v].dist + ((root == vbase[v])? cost[e] : costrev[e])) &&
1996  vbase[v] != vbase[j] && vbase[v] != vbase[j + nnodes] )
1997  {
1998  correct(scip, heap, state, &count, path, k, v, e, ((root == vbase[v])? cost[e] : costrev[e]), FSP_MODE);
1999  vbase[k] = vbase[v];
2000  }
2001  v += nnodes;
2002  }
2003  }
2004  }
2005  }
2006  if( nnodes > 1 )
2007  {
2008  int jc;
2009  /* until the heap is empty */
2010  while( count > 0 )
2011  {
2012  /* get the next (i.e. a nearest) vertex of the heap */
2013  k = nearest(heap, state, &count, path);
2014 
2015  /* mark vertex k as removed from heap */
2016  state[k] = UNKNOWN;
2017 
2018  assert(k - dnnodes >= 0);
2019  /* iterate over all outgoing edges of vertex k */
2020  for( e = g->outbeg[k - dnnodes]; e != EAT_LAST; e = g->oeat[e] )
2021  {
2022  j = g->head[e];
2023  if( Is_term(g->term[j]) || !g->mark[j] )
2024  continue;
2025  jc = j + dnnodes;
2026  /* check whether the path (to j) including k is shorter than the so far best known */
2027  if( vbase[j] != vbase[k] && vbase[j + nnodes] != vbase[k]
2028  && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) ) /*TODO(state[jc])??*/
2029  {
2030  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
2031  vbase[jc] = vbase[k];
2032  }
2033  }
2034  }
2035  }
2036  return;
2037 }
2038 
2039 
2040 /* 4th next terminal to all non terminal nodes */
2041 void graph_get4next(
2042  SCIP* scip, /**< SCIP data structure */
2043  const GRAPH* g, /**< graph data structure */
2044  const SCIP_Real* cost, /**< edge costs */
2045  const SCIP_Real* costrev, /**< reversed edge costs */
2046  PATH* path, /**< path data struture (leading to first, second and third nearest terminal) */
2047  int* vbase, /**< first, second and third nearest terminal to each non terminal */
2048  int* heap, /**< array representing a heap */
2049  int* state /**< array to mark the state of each node during calculation */
2050  )
2051 {
2052  int k;
2053  int j;
2054  int i;
2055  int l;
2056  int v;
2057  int e;
2058  int root;
2059  int count = 0;
2060  int nnodes;
2061  int dnnodes;
2062  int tnnodes;
2063 
2064  assert(g != NULL);
2065  assert(path != NULL);
2066  assert(cost != NULL);
2067  assert(heap != NULL);
2068  assert(state != NULL);
2069  assert(costrev != NULL);
2070 
2071  root = g->source;
2072  nnodes = g->knots;
2073  dnnodes = 2 * nnodes;
2074  tnnodes = 3 * nnodes;
2075 
2076  /* initialize */
2077  for( i = 0; i < nnodes; i++ )
2078  {
2079  /* copy of node i */
2080  k = i + tnnodes;
2081  vbase[k] = UNKNOWN;
2082  state[k] = UNKNOWN;
2083  path[k].edge = UNKNOWN;
2084  path[k].dist = FARAWAY;
2085  }
2086 
2087  /* scan original nodes */
2088  for( i = 0; i < nnodes; i++ )
2089  {
2090  state[i] = CONNECT;
2091  state[i + nnodes] = CONNECT;
2092  state[i + dnnodes] = CONNECT;
2093  if( !g->mark[i] )
2094  continue;
2095 
2096  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
2097  {
2098  j = g->head[e];
2099  k = j + tnnodes;
2100  if( !Is_term(g->term[j]) && g->mark[j] )
2101  {
2102  v = i;
2103  for( l = 0; l < 3; l++ )
2104  {
2105  if( SCIPisGT(scip, path[k].dist, path[v].dist + ((root == vbase[v])? cost[e] : costrev[e])) &&
2106  vbase[v] != vbase[j] && vbase[v] != vbase[j + nnodes] && vbase[v] != vbase[j + dnnodes] )
2107  {
2108  correct(scip, heap, state, &count, path, k, v, e, ((root == vbase[v])? cost[e] : costrev[e]), FSP_MODE);
2109  vbase[k] = vbase[v];
2110  }
2111  v += nnodes;
2112  }
2113  }
2114  }
2115  }
2116  if( nnodes > 1 )
2117  {
2118  int jc;
2119  /* until the heap is empty */
2120  while( count > 0 )
2121  {
2122  /* get the next (i.e. a nearest) vertex of the heap */
2123  k = nearest(heap, state, &count, path);
2124 
2125  /* mark vertex k as removed from heap */
2126  state[k] = UNKNOWN;
2127 
2128  assert(k - tnnodes >= 0);
2129  /* iterate over all outgoing edges of vertex k */
2130  for( e = g->outbeg[k - tnnodes]; e != EAT_LAST; e = g->oeat[e] )
2131  {
2132  j = g->head[e];
2133  if( Is_term(g->term[j]) || !g->mark[j] )
2134  continue;
2135  jc = j + tnnodes;
2136  /* check whether the path (to j) including k is shorter than the so far best known */
2137  if( vbase[j] != vbase[k] && vbase[j + nnodes] != vbase[k] && vbase[j + dnnodes] != vbase[k]
2138  && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) ) /*TODO(state[jc])??*/
2139  {
2140  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
2141  vbase[jc] = vbase[k];
2142  }
2143  }
2144  }
2145  }
2146  return;
2147 }
2148 
2149 /** build a voronoi region in presolving, w.r.t. shortest paths, for all terminals*/
2150 void graph_get3nextTerms(
2151  SCIP* scip, /**< SCIP data structure */
2152  const GRAPH* g, /**< graph data structure */
2153  SCIP_Real* cost, /**< edge costs */
2154  SCIP_Real* costrev, /**< reversed edge costs */
2155  PATH* path3, /**< path data struture (leading to first, second and third nearest terminal) */
2156  int* vbase, /**< first, second and third nearest terminal to each non terminal */
2157  int* heap, /**< array representing a heap */
2158  int* state /**< array to mark the state of each node during calculation */
2159  )
2160 {
2161  int k;
2162  assert(g != NULL);
2163  assert(path3 != NULL);
2164  assert(cost != NULL);
2165  assert(costrev != NULL);
2166  assert(heap != NULL);
2167  assert(state != NULL);
2168 
2169  if( g->stp_type != STP_PCSPG && g->stp_type != STP_RPCSPG )
2170  for( k = 0; k < g->knots; k++ )
2171  g->mark[k] = (g->grad[k] > 0);
2172 
2173  /* build voronoi diagram */
2174  graph_voronoiTerms(scip, g, cost, path3, vbase, heap, state);
2175 
2176  /* get 2nd nearest terms */
2177  graph_get2next(scip, g, cost, costrev, path3, vbase, heap, state);
2178 
2179  /* get 3th nearest terms */
2180  graph_get3next(scip, g, cost, costrev, path3, vbase, heap, state);
2181 
2182  return;
2183 }
2184 
2185 /** build a voronoi region in presolving, w.r.t. shortest paths, for all terminals*/
2186 void graph_get4nextTerms(
2187  SCIP* scip, /**< SCIP data structure */
2188  const GRAPH* g, /**< graph data structure */
2189  SCIP_Real* cost, /**< edge costs */
2190  SCIP_Real* costrev, /**< reversed edge costs */
2191  PATH* path, /**< path data struture (leading to first, second, third and fouth nearest terminal) */
2192  int* vbase, /**< first, second and third nearest terminal to each non terminal */
2193  int* heap, /**< array representing a heap */
2194  int* state /**< array to mark the state of each node during calculation */
2195  )
2196 {
2197  int k;
2198 
2199  assert(g != NULL);
2200  assert(path != NULL);
2201  assert(cost != NULL);
2202  assert(heap != NULL);
2203  assert(state != NULL);
2204  assert(costrev != NULL);
2205 
2206  if( g->stp_type != STP_PCSPG && g->stp_type != STP_RPCSPG )
2207  for( k = 0; k < g->knots; k++ )
2208  g->mark[k] = (g->grad[k] > 0);
2209 
2210  /* build voronoi diagram */
2211  graph_voronoiTerms(scip, g, cost, path, vbase, heap, state);
2212 
2213  /* get 2nd nearest terms */
2214  graph_get2next(scip, g, cost, costrev, path, vbase, heap, state);
2215 
2216  /* get 3th nearest terms */
2217  graph_get3next(scip, g, cost, costrev, path, vbase, heap, state);
2218 
2219  /* get 4th nearest terms */
2220  graph_get4next(scip, g, cost, costrev, path, vbase, heap, state);
2221 
2222  return;
2223 }
2224 
2225 
2226 /** get 4 close terminals to each terminal */
2227 SCIP_RETCODE graph_get4nextTTerms(
2228  SCIP* scip, /**< SCIP data structure */
2229  const GRAPH* g, /**< graph data structure */
2230  SCIP_Real* cost, /**< edge costs */
2231  PATH* path, /**< path data structure (leading to first, second, third and fourth nearest terminal) */
2232  int* vbase, /**< first, second and third nearest terminal to each non terminal */
2233  int* heap, /**< array representing a heap */
2234  int* state /**< array to mark the state of each node during calculation */
2235  )
2236 {
2237  int* boundedges;
2238  int k;
2239  int e;
2240  int l;
2241  int k2;
2242  int bedge;
2243  int shift;
2244  int nnodes;
2245  int nboundedges;
2246 
2247  assert(g != NULL);
2248  assert(path != NULL);
2249  assert(cost != NULL);
2250  assert(heap != NULL);
2251  assert(state != NULL);
2252 
2253  SCIP_CALL( SCIPallocBufferArray(scip, &boundedges, g->edges) );
2254 
2255  shift = 0;
2256  nnodes = g->knots;
2257 
2258  if( g->stp_type != STP_PCSPG && g->stp_type != STP_RPCSPG )
2259  for( k = 0; k < g->knots; k++ )
2260  g->mark[k] = (g->grad[k] > 0);
2261 
2262  nboundedges = 0;
2263  for( k = 0; k < nnodes; k ++ )
2264  {
2265  if( !g->mark[k] )
2266  continue;
2267 
2268  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2269  {
2270  k2 = g->head[e];
2271  if( !g->mark[k2] || k2 < k )
2272  continue;
2273 
2274  /* is e a boundary edge? */
2275  if( vbase[k] != vbase[k2] )
2276  {
2277  boundedges[nboundedges++] = e;
2278  }
2279  }
2280  if( Is_term(g->term[k]) )
2281  {
2282  path[k].dist = FARAWAY;
2283  vbase[k] = UNKNOWN;
2284  }
2285 
2286  }
2287 
2288  for( l = 0; l < 4; l++ )
2289  {
2290  for( e = 0; e < nboundedges; e++ )
2291  {
2292  bedge = boundedges[e];
2293  k = g->tail[bedge];
2294  k2 = g->head[bedge];
2295  utdist(scip, g, path, cost[bedge], vbase, k, l, k2, shift, nnodes);
2296  utdist(scip, g, path, cost[bedge], vbase, k2, l, k, shift, nnodes);
2297  }
2298  shift += nnodes;
2299  }
2300 
2301  SCIPfreeBufferArray(scip, &boundedges);
2302 
2303  return SCIP_OKAY;
2304 }
2305 
2306 /** build a Voronoi region in presolving, w.r.t. shortest paths, for all terminals */
2307 void graph_voronoiTerms(
2308  SCIP* scip, /**< SCIP data structure */
2309  const GRAPH* g, /**< graph data structure */
2310  const SCIP_Real* cost, /**< edge costs */
2311  PATH* path, /**< path data structure (leading to respective Voronoi base) */
2312  int* vbase, /**< Voronoi base to each vertex */
2313  int* heap, /**< array representing a heap */
2314  int* state /**< array to mark the state of each node during calculation */
2315  )
2316 {
2317  int k;
2318  int m;
2319  int i;
2320  int count = 0;
2321  int nbases = 0;
2322 
2323  assert(g != NULL);
2324  assert(path != NULL);
2325  assert(cost != NULL);
2326  assert(heap != NULL);
2327  assert(state != NULL);
2328 
2329  /* initialize */
2330  for( i = 0; i < g->knots; i++ )
2331  {
2332  /* set the base of vertex i */
2333  if( Is_term(g->term[i]) && g->mark[i] )
2334  {
2335  nbases++;
2336  if( g->knots > 1 )
2337  heap[++count] = i;
2338  vbase[i] = i;
2339  path[i].dist = 0.0;
2340  path[i].edge = UNKNOWN;
2341  state[i] = count;
2342  }
2343  else
2344  {
2345  vbase[i] = UNKNOWN;
2346  path[i].dist = FARAWAY;
2347  path[i].edge = UNKNOWN;
2348  state[i] = UNKNOWN;
2349  }
2350  }
2351  if( nbases == 0 )
2352  return;
2353  if( g->knots > 1 )
2354  {
2355  /* until the heap is empty */
2356  while( count > 0 )
2357  {
2358  /* get the next (i.e. a nearest) vertex of the heap */
2359  k = nearest(heap, state, &count, path);
2360 
2361  /* mark vertex k as scanned */
2362  state[k] = CONNECT;
2363 
2364  /* iterate over all outgoing edges of vertex k */
2365  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2366  {
2367  m = g->head[i];
2368 
2369  /* check whether the path (to m) including k is shorter than the so far best known */
2370  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + cost[i]) && g->mark[m] )
2371  {
2372  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
2373  vbase[m] = vbase[k];
2374  }
2375  }
2376  }
2377  }
2378  return;
2379 }
2380 
2381 
2382 /** build a Voronoi region, w.r.t. shortest paths, for all positive vertices */
2383 void graph_voronoiMw(
2384  SCIP* scip, /**< SCIP data structure */
2385  const GRAPH* g, /**< graph data structure */
2386  const SCIP_Real* costrev, /**< reversed edge costs */
2387  PATH* path, /**< path data structure (leading to respective Voronoi base) */
2388  int* vbase, /**< Voronoi base to each vertex */
2389  int* heap, /**< array representing a heap */
2390  int* state /**< array to mark the state of each node during calculation */
2391  )
2392 {
2393  int k;
2394  int m;
2395  int i;
2396  int count = 0;
2397  int nbases = 0;
2398  int nnodes;
2399 
2400  assert(g != NULL);
2401  assert(path != NULL);
2402  assert(costrev != NULL);
2403  assert(heap != NULL);
2404  assert(state != NULL);
2405 
2406  nnodes = g->knots;
2407 
2408  /* initialize */
2409  for( i = 0; i < nnodes; i++ )
2410  {
2411  /* set the base of vertex i */
2412  if( Is_term(g->term[i]) && g->mark[i] )
2413  {
2414  nbases++;
2415  if( nnodes > 1 )
2416  heap[++count] = i;
2417  vbase[i] = i;
2418  path[i].dist = -g->prize[i];
2419  path[i].edge = UNKNOWN;
2420  state[i] = count;
2421  }
2422  else
2423  {
2424  vbase[i] = UNKNOWN;
2425  path[i].dist = FARAWAY;
2426  path[i].edge = UNKNOWN;
2427  state[i] = UNKNOWN;
2428  }
2429  }
2430  if( nbases == 0 )
2431  return;
2432  if( nnodes > 1 )
2433  {
2434  /* until the heap is empty */
2435  while( count > 0 )
2436  {
2437  /* get the next (i.e. a nearest) vertex of the heap */
2438  k = nearest(heap, state, &count, path);
2439 
2440  /* mark vertex k as scanned */
2441  state[k] = CONNECT;
2442 
2443  /* iterate over all outgoing edges of vertex k */
2444  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2445  {
2446  m = g->head[i];
2447 
2448  /* check whether the path (to m) including k is shorter than the so far best known */
2449  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + costrev[i]) && g->mark[m] && !Is_term(g->term[m]) )
2450  {
2451  assert(vbase[m] != m);
2452  correct(scip, heap, state, &count, path, m, k, i, costrev[i], FSP_MODE);
2453  vbase[m] = vbase[k];
2454  }
2455  }
2456  }
2457  }
2458  return;
2459 }
2460 
2461 
2462 /** build a voronoi region, w.r.t. shortest paths, for all terminal and the distance */
2463 SCIP_RETCODE graph_voronoiWithDist(
2464  SCIP* scip, /**< SCIP data structure */
2465  const GRAPH* g, /**< graph data structure */
2466  SCIP_Real* cost, /**< edge costs */
2467  SCIP_Real* distance, /**< array storing path from a terminal over shortest
2468  incident edge to nearest terminal */
2469  int* minedgepred, /**< shortest edge predecessor array */
2470  int* vbase, /**< array containing Voronoi base to each node */
2471  int* minarc, /**< array to mark whether an edge is one a path corresponding to 'distance' */
2472  int* heap, /**< array representing a heap */
2473  int* state, /**< array indicating state of each vertex during calculation of Voronoi regions */
2474  int* distnode, /**< array to store terminal corresponding to distance stored in distance array */
2475  PATH* path /**< array containing Voronoi paths data */
2476  )
2477 {
2478  SCIP_Real new;
2479  int e;
2480  int k;
2481  int m;
2482  int i;
2483  int pred;
2484  int count;
2485  int nbases;
2486  int nnodes;
2487  int prededge;
2488 
2489  assert(g != NULL);
2490  assert(path != NULL);
2491  assert(cost != NULL);
2492  assert(heap != NULL);
2493  assert(state != NULL);
2494  assert(distance != NULL);
2495 
2496  count = 0;
2497  nbases = 0;
2498  nnodes = g->knots;
2499 
2500  /* initialize */
2501  for( i = 0; i < nnodes; i++ )
2502  {
2503  distance[i] = FARAWAY;
2504  if( distnode != NULL )
2505  distnode[i] = UNKNOWN;
2506 
2507  /* set the base of vertex i */
2508  if( Is_term(g->term[i]) && g->mark[i] )
2509  {
2510  nbases++;
2511  if( nnodes > 1 )
2512  heap[++count] = i;
2513  vbase[i] = i;
2514  state[i] = count;
2515  path[i].dist = 0.0;
2516  }
2517  else
2518  {
2519  vbase[i] = UNKNOWN;
2520  state[i] = UNKNOWN;
2521  path[i].dist = FARAWAY;
2522  }
2523  path[i].edge = UNKNOWN;
2524  }
2525 
2526  /* initialize predecessor array */
2527 
2528  for( e = 0; e < g->edges; e++ )
2529  minedgepred[e] = FALSE;
2530 
2531  for( k = 0; k < nbases; k++ )
2532  if( minarc[k] != UNKNOWN )
2533  minedgepred[minarc[k]] = TRUE;
2534 
2535  /* if graph contains more than one vertices, start main loop */
2536  if( nnodes > 1 )
2537  {
2538  /* until the heap is empty */
2539  while( count > 0 )
2540  {
2541  /* get the next (i.e. a nearest) vertex of the heap */
2542  k = nearest(heap, state, &count, path);
2543 
2544  /* mark vertex k as scanned */
2545  state[k] = CONNECT;
2546 
2547  prededge = path[k].edge;
2548 
2549  if( prededge != UNKNOWN )
2550  {
2551  pred = g->tail[prededge];
2552 
2553  assert(vbase[k] != UNKNOWN);
2554  assert(vbase[pred] != UNKNOWN);
2555  assert(vbase[pred] == vbase[k]);
2556  assert(g->head[prededge] == k);
2557 
2558  if( !Is_term(g->term[pred]) && g->mark[pred] )
2559  {
2560  assert(path[pred].edge != UNKNOWN);
2561  minedgepred[prededge] = minedgepred[path[pred].edge];
2562  }
2563 
2564  }
2565 
2566  /* iterate over all outgoing edges of vertex k */
2567  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2568  {
2569  m = g->head[i];
2570 
2571  if( state[m] == CONNECT && vbase[m] != vbase[k] && g->mark[m] )
2572  {
2573  if( minedgepred[i] || (prededge != UNKNOWN && minedgepred[prededge] ) )
2574  {
2575  new = path[k].dist + g->cost[i] + path[m].dist;
2576  if( SCIPisLT(scip, new, distance[vbase[k]]) )
2577  {
2578  if( distnode != NULL )
2579  distnode[vbase[k]] = vbase[m];
2580 
2581  distance[vbase[k]] = new;
2582  }
2583  }
2584  if( minedgepred[Edge_anti(i)] || (path[m].edge != UNKNOWN && minedgepred[path[m].edge] ) )
2585  {
2586  new = path[m].dist + g->cost[i] + path[k].dist;
2587  if( SCIPisLT(scip, new, distance[vbase[m]]) )
2588  {
2589  if( distnode != NULL )
2590  distnode[vbase[m]] = vbase[k];
2591 
2592  distance[vbase[m]] = new;
2593  }
2594  }
2595  }
2596 
2597  /* Check whether the path (to m) including k is shorter than the so far best known.
2598  In case of equality, also update if k is sucessor of minedge. */
2599  if( state[m] && g->mark[m] &&
2600  (SCIPisGT(scip, path[m].dist, path[k].dist + cost[i]) ||
2601  (prededge != UNKNOWN && minedgepred[prededge] && SCIPisEQ(scip, path[m].dist, path[k].dist + cost[i]))) )
2602  {
2603  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
2604  vbase[m] = vbase[k];
2605  }
2606  }
2607  }
2608  }
2609 
2610  return SCIP_OKAY;
2611 }
2612 
2613 /** build voronoi regions, w.r.t. shortest paths, for all terminals and compute the radii */
2614 SCIP_RETCODE graph_voronoiWithRadius(
2615  SCIP* scip, /**< SCIP data structure */
2616  const GRAPH* graph, /**< graph data structure */
2617  GRAPH* adjgraph, /**< graph data structure */
2618  PATH* path, /**< array containing Voronoi paths data */
2619  SCIP_Real* rad, /**< array storing shortest way from a terminal out of its Voronoi region */
2620  SCIP_Real* cost, /**< edge costs */
2621  SCIP_Real* costrev, /**< reversed edge costs */
2622  int* vbase, /**< array containing Voronoi base of each node */
2623  int* heap, /**< array representing a heap */
2624  int* state /**< array to mark state of each node during calculation */
2625  )
2626 {
2627  int* nodesid;
2628  int i;
2629  int k;
2630  int m;
2631  int vbm;
2632  int vbk;
2633  int root;
2634  int count = 0;
2635  int nnodes;
2636  int nterms = 0;
2637  STP_Bool pc;
2638  STP_Bool mw;
2639 
2640  assert(graph != NULL);
2641  assert(heap != NULL);
2642  assert(state != NULL);
2643  assert(path != NULL);
2644  assert(cost != NULL);
2645  assert(costrev != NULL);
2646  assert(rad != NULL);
2647  assert(vbase != NULL);
2648 
2649  nnodes = graph->knots;
2650 
2651  if( graph->terms == 0 )
2652  return SCIP_OKAY;
2653 
2654  root = graph->source;
2655  mw = (graph->stp_type == STP_MWCSP);
2656  pc = ((graph->stp_type == STP_PCSPG) || (graph->stp_type == STP_RPCSPG));
2657 
2658  SCIP_CALL( SCIPallocBufferArray(scip, &nodesid, nnodes) );
2659 
2660  /* initialize */
2661  for( i = 0; i < nnodes; i++ )
2662  {
2663  rad[i] = FARAWAY;
2664 
2665  /* set the base of vertex i */
2666  if( Is_term(graph->term[i]) && graph->mark[i] )
2667  {
2668  if( nnodes > 1 )
2669  heap[++count] = i;
2670 
2671  if( !mw )
2672  {
2673  adjgraph->mark[adjgraph->knots] = TRUE;
2674  graph_knot_add(adjgraph, -1);
2675  }
2676 
2677  nodesid[i] = nterms++;
2678  vbase[i] = i;
2679 
2680  if( mw )
2681  assert(SCIPisGE(scip, graph->prize[i], 0.0));
2682  if( mw )
2683  path[i].dist = -graph->prize[i];
2684  else
2685  path[i].dist = 0.0;
2686 
2687  path[i].edge = UNKNOWN;
2688  state[i] = count;
2689  }
2690  else
2691  {
2692  vbase[i] = UNKNOWN;
2693  nodesid[i] = UNKNOWN;
2694  path[i].dist = FARAWAY;
2695  path[i].edge = UNKNOWN;
2696  state[i] = UNKNOWN;
2697  }
2698  }
2699 
2700  if( nnodes > 1 )
2701  {
2702  SCIP_Real c1;
2703  SCIP_Real c2;
2704  SCIP_Real ecost;
2705  int ne;
2706 
2707  /* until the heap is empty */
2708  while( count > 0 )
2709  {
2710  /* get the next (i.e. a nearest) vertex of the heap */
2711  k = nearest(heap, state, &count, path);
2712 
2713  /* mark vertex k as scanned */
2714  state[k] = CONNECT;
2715 
2716  /* iterate over all outgoing edges of vertex k */
2717  for( i = graph->outbeg[k]; i != EAT_LAST; i = graph->oeat[i] )
2718  {
2719  m = graph->head[i];
2720  vbm = vbase[m];
2721  vbk = vbase[k];
2722 
2723  if( state[m] == CONNECT && vbm != vbk && graph->mark[m] )
2724  {
2725  assert(graph->mark[vbm]);
2726  assert(graph->mark[vbk]);
2727  if( mw )
2728  {
2729  if( SCIPisGT(scip, rad[vbk], path[k].dist) )
2730  rad[vbk] = path[k].dist;
2731  if( SCIPisGT(scip, rad[vbm], path[m].dist) )
2732  rad[vbm] = path[m].dist;
2733 #if 0
2734  if( SCIPisGT(scip, path[m].dist + graph->prize[vbm], path[k].dist + graph->prize[vbk]) )
2735  ecost = graph->prize[vbm] - path[k].dist;
2736  else
2737  ecost = graph->prize[vbk] - path[m].dist;
2738  if( SCIPisLT(scip, ecost, 0.0) )
2739  ecost = 0.0;
2740  /* find edge in adjgraph */
2741  for( ne = adjgraph->outbeg[nodesid[vbk]]; ne != EAT_LAST; ne = adjgraph->oeat[ne] )
2742  if( adjgraph->head[ne] == nodesid[vbm] )
2743  break;
2744 
2745  /* edge exists? */
2746  if( ne != EAT_LAST )
2747  {
2748  assert(ne >= 0);
2749  assert(adjgraph->head[ne] == nodesid[vbm]);
2750  assert(adjgraph->tail[ne] == nodesid[vbk]);
2751  if( SCIPisLT(scip, adjgraph->cost[ne], ecost) )
2752  {
2753  adjgraph->cost[ne] = ecost;
2754  adjgraph->cost[Edge_anti(ne)] = ecost;
2755  }
2756  }
2757  else
2758  {
2759  graph_edge_add(scip, adjgraph, nodesid[vbm], nodesid[vbk], ecost, ecost);
2760  }
2761 #endif
2762  }
2763  else
2764  {
2765  if( graph->stp_type == STP_DHCSTP )
2766  {
2767  if( m == root )
2768  c1 = path[m].dist + costrev[i];
2769  else
2770  c1 = path[m].dist + cost[i];
2771  if( k == root )
2772  c2 = path[k].dist + cost[i];
2773  else
2774  c2 = path[k].dist + costrev[i];
2775 
2776  if( SCIPisGT(scip, c1, c2) )
2777  ecost = c2;
2778  else
2779  ecost = c1;
2780  }
2781  else
2782  {
2783  if( SCIPisGT(scip, path[m].dist, path[k].dist) )
2784  ecost = path[k].dist + cost[i];
2785  else
2786  ecost = path[m].dist + cost[i];
2787  }
2788 
2789  if( pc )
2790  {
2791  if( SCIPisGT(scip, ecost, graph->prize[vbm])
2792  && root != vbm )
2793  ecost = graph->prize[vbm];
2794  if( SCIPisGT(scip, ecost, graph->prize[vbk])
2795  && root != vbk )
2796  ecost = graph->prize[vbk];
2797  }
2798 
2799  /* find edge in adjgraph */
2800  for( ne = adjgraph->outbeg[nodesid[vbk]]; ne != EAT_LAST; ne = adjgraph->oeat[ne] )
2801  if( adjgraph->head[ne] == nodesid[vbm] )
2802  break;
2803 
2804  /* edge exists? */
2805  if( ne != EAT_LAST )
2806  {
2807  assert(ne >= 0);
2808  assert(adjgraph->head[ne] == nodesid[vbm]);
2809  assert(adjgraph->tail[ne] == nodesid[vbk]);
2810  if( SCIPisGT(scip, adjgraph->cost[ne], ecost) )
2811  {
2812  adjgraph->cost[ne] = ecost;
2813  adjgraph->cost[Edge_anti(ne)] = ecost;
2814  }
2815  }
2816  else
2817  {
2818  graph_edge_add(scip, adjgraph, nodesid[vbm], nodesid[vbk],
2819  ecost, ecost);
2820  }
2821 
2822  if( SCIPisGT(scip, rad[vbk],
2823  path[k].dist + ((vbk == root) ? cost[i] : costrev[i])) )
2824  rad[vbk] = path[k].dist
2825  + ((vbk == root) ? cost[i] : costrev[i]);
2826 
2827  if( SCIPisGT(scip, rad[vbm],
2828  path[m].dist + ((vbm == root) ? costrev[i] : cost[i])) )
2829  rad[vbm] = path[m].dist
2830  + ((vbm == root) ? costrev[i] : cost[i]);
2831  }
2832  }
2833 
2834  /* check whether the path (to m) including k is shorter than the so far best known */
2835  if( state[m] && graph->mark[m] && !Is_term(graph->term[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + ((vbk == root)? cost[i] : costrev[i])) )
2836  {
2837  correct(scip, heap, state, &count, path, m, k, i, ((vbk == root)? cost[i] : costrev[i]), FSP_MODE);
2838  vbase[m] = vbk;
2839  }
2840  }
2841  }
2842  }
2843  SCIPfreeBufferArray(scip, &nodesid);
2844 
2845  return SCIP_OKAY;
2846 }
2847 
2848 /** Build partition of graph for MWCSP into regions around the positive vertices.
2849  * Store the weight of a minimum weight center-boundary path for each region
2850  * in the radius array (has to be reverted to obtain the final r() value).
2851  */
2852 void graph_voronoiWithRadiusMw(
2853  SCIP* scip, /**< SCIP data structure */
2854  const GRAPH* g, /**< graph data structure */
2855  PATH* path, /**< array containing graph decomposition data */
2856  const SCIP_Real* cost, /**< edge costs */
2857  SCIP_Real* radius, /**< array storing shortest way from a positive vertex out of its region */
2858  int* vbase, /**< array containing base of each node */
2859  int* heap, /**< array representing a heap */
2860  int* state /**< array to mark state of each node during calculation */
2861  )
2862 {
2863  SCIP_Real* nodeweight;
2864  int i;
2865  int k;
2866  int m;
2867  int count;
2868  int nbases;
2869  int nnodes;
2870  int nterms;
2871 
2872  assert(g != NULL);
2873  assert(heap != NULL);
2874  assert(state != NULL);
2875  assert(path != NULL);
2876  assert(cost != NULL);
2877  assert(radius != NULL);
2878  assert(vbase != NULL);
2879 
2880  count = 0;
2881  nbases = 0;
2882  nnodes = g->knots;
2883  nterms = g->terms - 1;
2884  nodeweight = g->prize;
2885 
2886  assert(nodeweight != NULL);
2887 
2888  if( nnodes == 0 || nterms <= 0 )
2889  return;
2890 
2891  assert(!g->mark[g->source]);
2892 
2893  /* initialize data */
2894  for( i = 0; i < nnodes; i++ )
2895  {
2896  radius[i] = FARAWAY;
2897 
2898  /* set the base of vertex i */
2899  if( Is_term(g->term[i]) && g->mark[i] )
2900  {
2901  nbases++;
2902  if( nnodes > 1 )
2903  heap[++count] = i;
2904  vbase[i] = i;
2905  path[i].dist = 0.0;
2906  path[i].edge = UNKNOWN;
2907  state[i] = count;
2908  }
2909  else
2910  {
2911  vbase[i] = UNKNOWN;
2912  path[i].dist = FARAWAY;
2913  path[i].edge = UNKNOWN;
2914  state[i] = UNKNOWN;
2915  }
2916  }
2917 
2918  if( nbases == 0 )
2919  return;
2920 
2921  if( nnodes > 1 )
2922  {
2923  int basem;
2924  int basek;
2925 
2926  /* until the heap is empty */
2927  while( count > 0 )
2928  {
2929  /* get the next (i.e. a nearest) vertex of the heap */
2930  k = nearest(heap, state, &count, path);
2931 
2932  /* mark vertex k as scanned */
2933  state[k] = CONNECT;
2934  basek = vbase[k];
2935 
2936  assert(g->mark[basek]);
2937 
2938  /* iterate over all outgoing edges of vertex k */
2939  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2940  {
2941  m = g->head[i];
2942 
2943  if( !(g->mark[m]) )
2944  continue;
2945 
2946  basem = vbase[m];
2947 
2948  /* are both m and k finally in different regions? */
2949  if( basem != basek && (state[m] == CONNECT || vbase[m] == m ||
2950  SCIPisGE(scip, path[k].dist, nodeweight[basek])) )
2951  {
2952  if( SCIPisGT(scip, radius[basek], path[k].dist) )
2953  radius[basek] = path[k].dist;
2954  if( (state[m] == CONNECT || vbase[m] == m) && SCIPisGT(scip, radius[basem], path[m].dist) )
2955  {
2956  assert(g->mark[basem]);
2957  radius[basem] = path[m].dist;
2958  }
2959  }
2960 
2961  /* is the distance of vertex k smaller than the weight of its base node and smaller than the temp. radius? */
2962  if( SCIPisLT(scip, path[k].dist, nodeweight[basek]) && SCIPisLT(scip, path[k].dist, radius[basek]) )
2963  {
2964  /* check whether the path (to m) including k is shorter than the so far best known */
2965  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + cost[i]) )
2966  {
2967  assert(vbase[m] != m);
2968  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
2969  vbase[m] = basek;
2970  }
2971  }
2972  }
2973  }
2974  }
2975  return;
2976 }
2977 
2978 /** repair a Voronoi diagram for a given set of base nodes */
2979 void graph_voronoiRepair(
2980  SCIP* scip, /**< SCIP data structure */
2981  const GRAPH* g, /**< graph data structure */
2982  SCIP_Real* cost, /**< edge costs */
2983  int* count, /**< pointer to number of heap elements */
2984  int* vbase, /**< array containing Voronoi base of each node */
2985  PATH* path, /**< Voronoi paths data struture */
2986  int* newedge, /**< the new edge */
2987  int crucnode, /**< the current crucial node */
2988  UF* uf /**< union find data structure */
2989  )
2990 {
2991  int k;
2992  int m;
2993  int i;
2994  int e;
2995  int* heap;
2996  int* state;
2997  int node1;
2998  int node2;
2999 
3000  *newedge = UNKNOWN;
3001  e = UNKNOWN;
3002  assert(g != NULL);
3003  assert(g->path_heap != NULL);
3004  assert(g->path_state != NULL);
3005  assert(path != NULL);
3006  assert(cost != NULL);
3007 
3008  if( g->knots == 0 )
3009  return;
3010 
3011  heap = g->path_heap;
3012  state = g->path_state;
3013 
3014  if( g->knots > 1 )
3015  {
3016  /* until the heap is empty */
3017  while( *count > 0 )
3018  {
3019  /* get the next (i.e. a nearest) vertex from the heap */
3020  k = nearest(heap, state, count, path);
3021 
3022  /* mark vertex k as scanned */
3023  state[k] = CONNECT;
3024  assert(g->mark[k]);
3025  assert(g->mark[vbase[k]]);
3026  /* iterate over all outgoing edges of vertex k */
3027  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
3028  {
3029  m = g->head[i];
3030 
3031  /* check whether the path (to m) including k is shorter than the so far best known */
3032  if( (state[m]) && (SCIPisGT(scip, path[m].dist, path[k].dist + cost[i])) )/*&& g->mark[m] ) */
3033  {
3034  assert(g->mark[m]);
3035  correct(scip, heap, state, count, path, m, k, i, cost[i], FSP_MODE);
3036  vbase[m] = vbase[k];
3037  }
3038 
3039  /* check whether there is a better new boundary edge adjacent to vertex k */
3040  else
3041  {
3042  node1 = SCIPStpunionfindFind(uf, vbase[m]);
3043  node2 = SCIPStpunionfindFind(uf, vbase[k]);
3044  if( state[m] == CONNECT && ((node1 == crucnode) != (node2 == crucnode)) && g->mark[m] && g->mark[vbase[m]]
3045  && g->mark[node1] && g->mark[node2] && ((e == UNKNOWN) || SCIPisGT(scip,
3046  (path[g->tail[e]].dist + cost[e] + path[g->head[e]].dist), path[k].dist + cost[i] + path[m].dist)) )
3047  e = i;
3048  }
3049  }
3050  }
3051  }
3052  *newedge = e;
3053 }
3054 
3055 
3056 /** repair the Voronoi diagram for a given set nodes */
3057 void graph_voronoiRepairMult(
3058  SCIP* scip, /**< SCIP data structure */
3059  const GRAPH* g, /**< graph data structure */
3060  SCIP_Real* cost, /**< edge costs */
3061  int* count, /**< pointer to number of heap elements */
3062  int* vbase, /**< array containing Voronoi base of each node */
3063  int* boundedges, /**< boundary edges */
3064  int* nboundedges, /**< number of boundary edges */
3065  STP_Bool* nodesmark, /**< array to mark temporarily discarded nodes */
3066  UF* uf, /**< union find data structure */
3067  PATH* path /**< Voronoi paths data structure */
3068  )
3069 {
3070  int k;
3071  int m;
3072  int i;
3073  int* heap;
3074  int* state;
3075  int node1;
3076  int node2;
3077 
3078  assert(g != NULL);
3079  assert(g->mark != NULL);
3080  assert(g->path_heap != NULL);
3081  assert(g->path_state != NULL);
3082  assert(path != NULL);
3083  assert(cost != NULL);
3084 
3085  if( g->knots == 0 )
3086  return;
3087 
3088  heap = g->path_heap;
3089  state = g->path_state;
3090 
3091  assert(heap != NULL);
3092  assert(state != NULL);
3093 
3094  if( g->knots > 1 )
3095  {
3096  /* until the heap is empty */
3097  while( (*count) > 0 )
3098  {
3099  /* get the next (i.e. a nearest) vertex from the heap */
3100  k = nearest(heap, state, count, path);
3101 
3102  /* mark vertex k as scanned */
3103  state[k] = CONNECT;
3104 
3105  /* iterate all outgoing edges of vertex k */
3106  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
3107  {
3108  m = g->head[i];
3109 #if 1
3110  if( !(g->mark[m]) )
3111  continue;
3112 #endif
3113  /* check whether the path (to m) including k is shorter than the so far best known */
3114  if( state[m] && (SCIPisGT(scip,path[m].dist, path[k].dist + cost[i])) )
3115  {
3116  correct(scip, heap, state, count, path, m, k, i, cost[i], FSP_MODE);
3117  vbase[m] = vbase[k];
3118  }
3119  /* check whether there is a new boundary edge adjacent to vertex k */
3120  else if( (state[m] == CONNECT) && ((node1 = SCIPStpunionfindFind(uf, vbase[m])) != (node2 = SCIPStpunionfindFind(uf, vbase[k])))
3121  && g->mark[node1] && g->mark[node2] && (nodesmark[node1] || nodesmark[node2]) )
3122  {
3123  boundedges[(*nboundedges)++] = i;
3124  }
3125  }
3126  }
3127  }
3128 }
graph_path_st_pcmw_full
void graph_path_st_pcmw_full(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, STP_Bool *connected)
Definition: grphpath.c:1316
SCIPisEQ
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:438
graph_get3nextTerms
void graph_get3nextTerms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path3, int *vbase, int *heap, int *state)
Definition: grphpath.c:2150
nterms
static volatile int nterms
Definition: interrupt.c:38
graph_voronoiWithDist
SCIP_RETCODE graph_voronoiWithDist(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *distance, int *minedgepred, int *vbase, int *minarc, int *heap, int *state, int *distnode, PATH *path)
Definition: grphpath.c:2463
graph_path_invroot
void graph_path_invroot(SCIP *scip, const GRAPH *g, int start, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge)
Definition: grphpath.c:849
GRAPH::head
int *RESTRICT head
Definition: grph.h:96
GRAPH
Definition: grph.h:57
GRAPH::source
int source
Definition: grph.h:67
SCIPisGE
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:490
SCIPfreeMemoryArray
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:69
shortest_path::edge
signed int edge
Definition: grph.h:146
MST_MODE
#define MST_MODE
Definition: grph.h:162
GRAPH::terms
int terms
Definition: grph.h:64
graph_path_st_rmw
void graph_path_st_rmw(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, STP_Bool *connected)
Definition: grphpath.c:1536
SCIPallocMemoryArray
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:53
graph_path_execX
void graph_path_execX(SCIP *scip, const GRAPH *g, int start, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge)
Definition: grphpath.c:781
UnionFind_Structure
Definition: misc_stp.h:65
EAT_LAST
#define EAT_LAST
Definition: grph.h:31
Edge_anti
#define Edge_anti(a)
Definition: grph.h:171
graph_sdPaths
void graph_sdPaths(SCIP *scip, const GRAPH *g, PATH *path, SCIP_Real *cost, SCIP_Real distlimit, int *heap, int *state, int *memlbl, int *nlbl, int tail, int head, int limit)
Definition: grphpath.c:567
FALSE
#define FALSE
Definition: def.h:73
GRAPH::inpbeg
int *RESTRICT inpbeg
Definition: grph.h:74
GRAPH::path_state
int *RESTRICT path_state
Definition: grph.h:119
TRUE
#define TRUE
Definition: def.h:72
SCIP_RETCODE
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
STP_DHCSTP
#define STP_DHCSTP
Definition: grph.h:48
SCIP_RandNumGen
Definition: struct_misc.h:259
shortest_path
Definition: grph.h:143
nearestX
static int nearestX(int *heap, int *state, int *count, const SCIP_Real *pathdist)
Definition: grphpath.c:92
correct
static void correct(SCIP *scip, int *heap, int *state, int *count, PATH *path, int l, int k, int e, SCIP_Real cost, int mode)
Definition: grphpath.c:147
graph_get3next
void graph_get3next(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1934
STP_PCSPG
#define STP_PCSPG
Definition: grph.h:40
graph_path_exec
void graph_path_exec(SCIP *scip, const GRAPH *p, const int mode, int start, const SCIP_Real *cost, PATH *path)
Definition: grphpath.c:491
graph_voronoi
void graph_voronoi(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, STP_Bool *base, int *vbase, PATH *path)
Definition: grphpath.c:1751
graph_voronoiExtend
SCIP_RETCODE graph_voronoiExtend(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, SCIP_Real **distarr, int **basearr, int **edgearr, STP_Bool *termsmark, int *reachednodes, int *nreachednodes, int *nodenterms, int nneighbterms, int base, int countex)
Definition: grphpath.c:1671
graph_path_st_pcmw
void graph_path_st_pcmw(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, STP_Bool *connected)
Definition: grphpath.c:1154
SCIPfreeBufferArray
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
graph_path_st_pcmw_extend
void graph_path_st_pcmw_extend(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, PATH *path, STP_Bool *connected, SCIP_Bool *extensions)
Definition: grphpath.c:1410
graph_voronoiTerms
void graph_voronoiTerms(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:2307
GRAPH::mark
int *RESTRICT mark
Definition: grph.h:70
SCIPisLE
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:464
graph_voronoiRepairMult
void graph_voronoiRepairMult(SCIP *scip, const GRAPH *g, SCIP_Real *cost, int *count, int *vbase, int *boundedges, int *nboundedges, STP_Bool *nodesmark, UF *uf, PATH *path)
Definition: grphpath.c:3057
GRAPH::oeat
int *RESTRICT oeat
Definition: grph.h:103
CONNECT
#define CONNECT
Definition: grph.h:154
graph_path_st
void graph_path_st(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, SCIP_RANDNUMGEN *randnumgen, STP_Bool *connected)
Definition: grphpath.c:926
GRAPH::stp_type
int stp_type
Definition: grph.h:127
Is_pterm
#define Is_pterm(a)
Definition: grph.h:169
STP_Bool
unsigned char STP_Bool
Definition: grph.h:52
GRAPH::prize
SCIP_Real * prize
Definition: grph.h:82
graph_get4next
void graph_get4next(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:2041
graph_path_exit
void graph_path_exit(SCIP *scip, GRAPH *g)
Definition: grphpath.c:466
shortest_path::dist
SCIP_Real dist
Definition: grph.h:145
GRAPH::grad
int *RESTRICT grad
Definition: grph.h:73
FSP_MODE
#define FSP_MODE
Definition: grph.h:163
NULL
#define NULL
Definition: lpi_spx1.cpp:155
graph_get4nextTerms
void graph_get4nextTerms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:2186
SCIP_OKAY
Definition: type_retcode.h:33
graph_voronoiMw
void graph_voronoiMw(SCIP *scip, const GRAPH *g, const SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:2383
graph_get2next
void graph_get2next(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1838
GRAPH::knots
int knots
Definition: grph.h:62
SCIP_CALL
#define SCIP_CALL(x)
Definition: def.h:370
resetX
static void resetX(SCIP *scip, SCIP_Real *pathdist, int *heap, int *state, int *count, int node)
Definition: grphpath.c:275
SCIPStpunionfindFind
int SCIPStpunionfindFind(UF *uf, int element)
Definition: misc_stp.c:728
LT
#define LT(a, b)
Definition: portab.h:64
Is_gterm
#define Is_gterm(a)
Definition: grph.h:170
FARAWAY
#define FARAWAY
Definition: grph.h:156
SCIPallocBufferArray
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
GT
#define GT(a, b)
Definition: portab.h:66
SCIP_Bool
#define SCIP_Bool
Definition: def.h:70
GRAPH::ieat
int *RESTRICT ieat
Definition: grph.h:102
GRAPH::path_heap
int *RESTRICT path_heap
Definition: grph.h:118
graph_voronoiWithRadius
SCIP_RETCODE graph_voronoiWithRadius(SCIP *scip, const GRAPH *graph, GRAPH *adjgraph, PATH *path, SCIP_Real *rad, SCIP_Real *cost, SCIP_Real *costrev, int *vbase, int *heap, int *state)
Definition: grphpath.c:2614
STP_MWCSP
#define STP_MWCSP
Definition: grph.h:47
GRAPH::tail
int *RESTRICT tail
Definition: grph.h:95
reset
static void reset(SCIP *scip, PATH *path, int *heap, int *state, int *count, int node)
Definition: grphpath.c:309
graph_voronoiRepair
void graph_voronoiRepair(SCIP *scip, const GRAPH *g, SCIP_Real *cost, int *count, int *vbase, PATH *path, int *newedge, int crucnode, UF *uf)
Definition: grphpath.c:2979
GRAPH::term
int *RESTRICT term
Definition: grph.h:68
graph_voronoiWithRadiusMw
void graph_voronoiWithRadiusMw(SCIP *scip, const GRAPH *g, PATH *path, const SCIP_Real *cost, SCIP_Real *radius, int *vbase, int *heap, int *state)
Definition: grphpath.c:2852
grph.h
includes various files containing graph methods used for Steiner tree problems
graph_path_st_rpc
void graph_path_st_rpc(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, STP_Bool *connected)
Definition: grphpath.c:1033
graph_path_init
SCIP_RETCODE graph_path_init(SCIP *scip, GRAPH *g)
Definition: grphpath.c:444
portab.h
Portable defintions.
r
SCIP_Real * r
Definition: circlepacking.c:50
Is_term
#define Is_term(a)
Definition: grph.h:168
correctX
static void correctX(SCIP *scip, int *heap, int *state, int *count, SCIP_Real *pathdist, int *pathedge, int l, int k, int e, SCIP_Real cost)
Definition: grphpath.c:199
GRAPH::cost
SCIP_Real * cost
Definition: grph.h:94
SCIP_Real
#define SCIP_Real
Definition: def.h:163
heap_add
void heap_add(int *heap, int *state, int *count, int node, PATH *path)
Definition: grphpath.c:244
SCIPisLT
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:451
SCIPisGT
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:477
GRAPH::outbeg
int *RESTRICT outbeg
Definition: grph.h:76
GRAPH::edges
int edges
Definition: grph.h:92
graph_knot_add
void graph_knot_add(GRAPH *, int)
Definition: grphbase.c:2200
graph_path_PcMwSd
void graph_path_PcMwSd(SCIP *scip, const GRAPH *g, PATH *path, SCIP_Real *cost, SCIP_Real distlimit, int *pathmaxnode, int *heap, int *state, int *stateblock, int *memlbl, int *nlbl, int tail, int head, int limit)
Definition: grphpath.c:664
UNKNOWN
#define UNKNOWN
Definition: sepa_mcf.c:4095
nnodes
#define nnodes
Definition: gastrans.c:65
graph_get4nextTTerms
SCIP_RETCODE graph_get4nextTTerms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:2227
utdist
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: grphpath.c:343
nearest
static int nearest(int *heap, int *state, int *count, const PATH *path)
Definition: grphpath.c:43
graph_edge_add
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
STP_RPCSPG
#define STP_RPCSPG
Definition: grph.h:41
graph_path_st_pcmw_reduce
void graph_path_st_pcmw_reduce(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *tmpnodeweight, int *result, int start, STP_Bool *connected)
Definition: grphpath.c:1264