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-2015 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file 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 inluding
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 
103  /* Heap shift down
104  * (Oberstes Element runter und korrigieren)
105  */
106  k = heap[1];
107  j = 1;
108  c = 2;
109  heap[1] = heap[(*count)--];
110  state[heap[1]] = 1;
111 
112  if ((*count) > 2)
113  if (LT(pathdist[heap[3]], pathdist[heap[2]]))
114  c++;
115 
116  while((c <= (*count)) && GT(pathdist[heap[j]], pathdist[heap[c]]))
117  {
118  t = heap[c];
119  heap[c] = heap[j];
120  heap[j] = t;
121  state[heap[j]] = j;
122  state[heap[c]] = c;
123  j = c;
124  c += c;
125 
126  if ((c + 1) <= (*count))
127  if (LT(pathdist[heap[c + 1]], pathdist[heap[c]]))
128  c++;
129  }
130  return(k);
131 }
132 
133 
134 /*---------------------------------------------------------------------------*/
135 /*--- Name : CORRECT heap ---*/
136 /*--- Function : Setzt ein neues Element auf den Heap, bzw. korrigiert ---*/
137 /*--- die Position eines vorhandenen Elementes ---*/
138 /*--- Parameter: Derzeitige Entfernungen und benutzten Kanten, ---*/
139 /*--- Neuer Knoten, Vorgaengerknoten, Kante von der man aus ---*/
140 /*--- den neuen Knoten erreicht, Kosten der Kante, ---*/
141 /*--- sowie Betriebsmodus ---*/
142 /*--- Returns : Nichts ---*/
143 /*---------------------------------------------------------------------------*/
144 inline static void correct(
145  SCIP* scip,
146  int* heap,
147  int* state,
148  int* count, /* pointer to store the number of elements on the heap */
149  PATH* path,
150  int l,
151  int k,
152  int e,
153  SCIP_Real cost,
154  int mode)
155 {
156  int t;
157  int c;
158  int j;
159 
160  path[l].dist = (mode == MST_MODE) ? cost : (path[k].dist + cost);
161  path[l].edge = e;
162 
163  /* new node? */
164  if( state[l] == UNKNOWN )
165  {
166  heap[++(*count)] = l;
167  state[l] = (*count);
168  }
169 
170  /* Heap shift up */
171  j = state[l];
172  c = j / 2;
173  while( (j > 1) && SCIPisGT(scip, path[heap[c]].dist, path[heap[j]].dist) )
174  {
175  t = heap[c];
176  heap[c] = heap[j];
177  heap[j] = t;
178  state[heap[j]] = j;
179  state[heap[c]] = c;
180  j = c;
181  c = j / 2;
182  }
183 }
184 
185 
186 /*---------------------------------------------------------------------------*/
187 /*--- Name : CORRECT heap ---*/
188 /*--- Function : Setzt ein neues Element auf den Heap, bzw. korrigiert ---*/
189 /*--- die Position eines vorhandenen Elementes ---*/
190 /*--- Parameter: Derzeitige Entfernungen und benutzten Kanten, ---*/
191 /*--- Neuer Knoten, Vorgaengerknoten, Kante von der man aus ---*/
192 /*--- den neuen Knoten erreicht, Kosten der Kante, ---*/
193 /*--- sowie Betriebsmodus ---*/
194 /*--- Returns : Nichts ---*/
195 /*---------------------------------------------------------------------------*/
196 inline static void correctX(
197  SCIP* scip,
198  int* heap,
199  int* state,
200  int* count, /* pointer to store the number of elements on the heap */
201  SCIP_Real* pathdist,
202  int* pathedge,
203  int l,
204  int k,
205  int e,
206  SCIP_Real cost
207  )
208 {
209  int t;
210  int c;
211  int j;
212 
213  pathdist[l] = (pathdist[k] + cost);
214  pathedge[l] = e;
215 
216  /* Ist der Knoten noch ganz frisch ?
217  */
218  if (state[l] == UNKNOWN)
219  {
220  heap[++(*count)] = l;
221  state[l] = (*count);
222  }
223 
224  /* Heap shift up
225  */
226  j = state[l];
227  c = j / 2;
228 
229  while( (j > 1) && SCIPisGT(scip, pathdist[heap[c]], pathdist[heap[j]]) )
230  {
231  t = heap[c];
232  heap[c] = heap[j];
233  heap[j] = t;
234  state[heap[j]] = j;
235  state[heap[c]] = c;
236  j = c;
237  c = j / 2;
238  }
239 }
240 
241 void heap_add(
242  int* heap, /* heaparray */
243  int* state,
244  int* count, /* pointer to store the number of elements on the heap */
245  int node, /* the node to be added */
246  PATH* path
247  )
248 {
249  int t;
250  int c;
251  int j;
252 
253  heap[++(*count)] = node;
254  state[node] = (*count);
255 
256  /* Heap shift up */
257  j = state[node];
258  c = j / 2;
259 
260  while((j > 1) && GT(path[heap[c]].dist, path[heap[j]].dist))
261  {
262  t = heap[c];
263  heap[c] = heap[j];
264  heap[j] = t;
265  state[heap[j]] = j;
266  state[heap[c]] = c;
267  j = c;
268  c = j / 2;
269  }
270 
271 }
272 inline static void resetX(
273  SCIP* scip,
274  SCIP_Real* pathdist,
275  int* heap,
276  int* state,
277  int* count,
278  int node
279  )
280 {
281  int t;
282  int c;
283  int j;
284 
285  pathdist[node] = 0.0;
286 
287  heap[++(*count)] = node;
288  state[node] = (*count);
289 
290  /* heap shift up */
291  j = state[node];
292  c = j / 2;
293 
294  while( (j > 1) && SCIPisGT(scip, pathdist[heap[c]], pathdist[heap[j]]) )
295  {
296  t = heap[c];
297  heap[c] = heap[j];
298  heap[j] = t;
299  state[heap[j]] = j;
300  state[heap[c]] = c;
301  j = c;
302  c = j / 2;
303  }
304 }
305 
306 inline static void utdist(
307  SCIP* scip,
308  const GRAPH* g,
309  PATH* path,
310  SCIP_Real ecost,
311  int* vbase,
312  int k,
313  int l,
314  int k2,
315  int shift,
316  int nnodes
317  )
318 {
319  SCIP_Real dist;
320  int vbk;
321  int vbk2;
322 
323  if( Is_term(g->term[k]) )
324  vbk = k;
325  else
326  vbk = vbase[k];
327 
328  if( l == 0 )
329  {
330  assert(shift == 0);
331 
332  dist = ecost;
333  if( !Is_term(g->term[k]) )
334  dist += path[k].dist;
335 
336  if( !Is_term(g->term[k2]) )
337  {
338  dist += path[k2].dist;
339  vbk2 = vbase[k2];
340  }
341  else
342  {
343  vbk2 = k2;
344  }
345 
346  if( SCIPisLT(scip, dist, path[vbk].dist) )
347  {
348  path[vbk].dist = dist;
349  vbase[vbk] = vbk2;
350  return;
351  }
352  }
353  else
354  {
355  int max;
356  int pos;
357  int r;
358  int s;
359  int t;
360 
361  pos = vbk + shift;
362  max = MIN((l + 1), 3);
363 
364  for( r = 0; r <= max; r++ )
365  {
366  if( Is_term(g->term[k2]) )
367  {
368  if( r == 0 )
369  t = k2;
370  else
371  break;
372  }
373  else
374  {
375  t = vbase[k2 + (r * nnodes)];
376  }
377  for( s = 0; s < l; s++ )
378  if( vbase[vbk + s * nnodes] == t )
379  break;
380  if( s < l || vbk == t )
381  continue;
382 
383  dist = ecost;
384  if( !Is_term(g->term[k]) )
385  dist += path[k].dist;
386  if( !Is_term(g->term[k2]) )
387  dist += path[k2 + (r * nnodes)].dist;
388 
389  if( SCIPisLT(scip, dist, path[pos].dist) )
390  {
391  path[pos].dist = dist;
392  vbase[pos] = t;
393  return;
394  }
395  }
396  }
397  return;
398 }
399 
400 /*---------------------------------------------------------------------------*/
401 /*--- Name : INIT shortest PATH algorithm ---*/
402 /*--- Function : Initialisiert den benoetigten Speicher fuer die ---*/
403 /*--- kuerzeste Wege Berechnung ---*/
404 /*--- Parameter: Graph in dem berechnet werden soll. ---*/
405 /*--- Returns : Nichts ---*/
406 /*---------------------------------------------------------------------------*/
407 SCIP_RETCODE graph_path_init(
408  SCIP* scip, /**< SCIP data structure */
409  GRAPH* g /**< graph data structure */
410  )
411 {
412  assert(g != NULL);
413  assert(g->path_heap == NULL);
414  assert(g->path_state == NULL);
415 
416  SCIP_CALL( SCIPallocMemoryArray(scip, &(g->path_heap), g->knots + 1) );
417  SCIP_CALL( SCIPallocMemoryArray(scip, &(g->path_state), g->knots) );
418 
419  return SCIP_OKAY;
420 }
421 
422 /*---------------------------------------------------------------------------*/
423 /*--- Name : EXIT shortest PATH algorithm ---*/
424 /*--- Function : Gibt den bei der initialisierung angeforderten Speicher ---*/
425 /*--- wieder frei ---*/
426 /*--- Parameter: Keine ---*/
427 /*--- Returns : Nichts ---*/
428 /*---------------------------------------------------------------------------*/
429 void graph_path_exit(
430  SCIP* scip, /**< SCIP data structure */
431  GRAPH* g /**< graph data structure */
432  )
433 {
434  assert(g->path_heap != NULL);
435  assert(g->path_state != NULL);
436 
437  SCIPfreeMemoryArray(scip, &(g->path_heap));
438  SCIPfreeMemoryArray(scip, &(g->path_state));
439 }
440 
441 /*---------------------------------------------------------------------------*/
442 /*--- Name : FIND shortest PATH / minimum spanning tree ---*/
443 /*--- Function : Dijkstras Algorithmus zur bestimmung eines kuerzesten ---*/
444 /*--- Weges in einem gerichteten Graphen. ---*/
445 /*--- Parameter: Graph, Nummer des Startknotens und Betriebsmodus ---*/
446 /*--- (Kuerzeste Wege oder minimal spannender Baum) ---*/
447 /*--- Returns : Liefert einen Vektor mit einem PATH Element je Knotem. ---*/
448 /*---- Setzt das .dist Feld zu jedem Knoten auf die Entfernung ---*/
449 /*--- zum Startknoten, sowie das .prev Feld auf den Knoten von ---*/
450 /*--- dem man zum Knoten gekommen ist. Im MST Modus steht im ---*/
451 /*--- .dist Feld die Entfernung zum naechsten Knoten. ---*/
452 /*---------------------------------------------------------------------------*/
453 void graph_path_exec(
454  SCIP* scip, /**< SCIP data structure */
455  const GRAPH* p, /**< graph data structure */
456  int mode, /**< shortest path (FSP_MODE) or minimum spanning tree (MST_MODE)? */
457  int start, /**< start vertex */
458  SCIP_Real* cost, /**< edge costs */
459  PATH* path /**< shortest paths data structure */
460  )
461 {
462  int k;
463  int m;
464  int i;
465  int* heap;
466  int* state;
467  int count;
468 
469  assert(scip != NULL);
470  assert(p != NULL);
471  assert(start >= 0);
472  assert(start < p->knots);
473  assert((mode == FSP_MODE) || (mode == MST_MODE));
474  assert(p->path_heap != NULL);
475  assert(p->path_state != NULL);
476  assert(path != NULL);
477  assert(cost != NULL);
478 
479  /* no nodes?, return*/
480  if (p->knots == 0)
481  return;
482 
483  heap = p->path_heap;
484  state = p->path_state;
485 
486  /* Erstmal alles auf null, unbekannt und weit weg
487  */
488  for(i = 0; i < p->knots; i++)
489  {
490  state[i] = UNKNOWN;
491  path[i].dist = FARAWAY + 1;
492  path[i].edge = -1;
493  }
494  /* Startknoten in den Heap
495  */
496  k = start;
497  path[k].dist = 0.0;
498 
499  /* Wenn nur ein Knoten drin ist funktioniert der Heap nicht sehr gut,
500  * weil dann fuer genau 0 Elemente Platz ist.
501  */
502  if( p->knots > 1 )
503  {
504  count = 1;
505  heap[count] = k;
506  state[k] = count;
507 
508  /* Wenn nichts mehr auf dem Heap ist, sind wir fertig
509  * und jetzt erstmal Hula Loop
510  */
511  while( count > 0 )
512  {
513  /* Na, wer ist der Naechste ?
514  */
515  k = nearest(heap, state, &count, path);
516 
517  /* Wieder einen erledigt
518  */
519  state[k] = CONNECT;
520 
521  /* Verbunden Knoten berichtigen ...
522  *
523  * Wenn ein Knoten noch nicht erledigt ist
524  * werden wir dann auf diesem Wege besser ?
525  */
526  for( i = p->outbeg[k]; i != EAT_LAST; i = p->oeat[i] )
527  {
528  m = p->head[i];
529  /* node not scanned and valid? */
530  if( (state[m]) && (p->mark[m]) )
531  {
532  /* closer than previously? */
533  if( SCIPisGT(scip, path[m].dist, (mode == MST_MODE) ? cost[i] : (path[k].dist + cost[i])) )
534  correct(scip, heap, state, &count, path, m, k, i, cost[i], mode);
535  }
536  }
537  }
538  }
539 }
540 
541 /** limited Dijkstra, stopping at terminals */
542 void sdpaths(
543  SCIP* scip, /**< SCIP data structure */
544  const GRAPH* g, /**< graph data structure */
545  PATH* path, /**< shortest paths data structure */
546  SCIP_Real* cost, /**< edge costs */
547  SCIP_Real distlimit, /**< distance limit of the search */
548  int* heap, /**< array representing a heap */
549  int* state, /**< array to indicate whether a node has been scanned during SP calculation */
550  int* memlbl, /**< array to save labelled nodes */
551  int* nlbl, /**< number of labelled nodes */
552  int tail, /**< tail of the edge */
553  int head, /**< head of the edge */
554  int limit /**< maximum number of edges to consider during execution */
555  )
556 {
557  int k;
558  int m;
559  int e;
560  int limit1;
561  int count;
562  int nchecks;
563 
564  assert(g != NULL);
565  assert(heap != NULL);
566  assert(path != NULL);
567  assert(cost != NULL);
568  assert(nlbl != NULL);
569  assert(memlbl != NULL);
570 
571  limit1 = limit - limit / 3;
572  *nlbl = 0;
573  nchecks = 0;
574 
575  if( g->grad[tail] == 0 || g->grad[head] == 0 )
576  return;
577 
578  assert(g->mark[head] && g->mark[tail]);
579 
580  count = 0;
581  path[tail].dist = 0.0;
582  state[tail] = CONNECT;
583  memlbl[(*nlbl)++] = tail;
584 
585  if( g->stp_type != STP_MAX_NODE_WEIGHT )
586  g->mark[head] = FALSE;
587 
588  for( e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
589  {
590  m = g->head[e];
591 
592  if( g->mark[m] )
593  {
594  assert(SCIPisGT(scip, path[m].dist, path[tail].dist + cost[e]));
595  /* m labelled the first time */
596  memlbl[(*nlbl)++] = m;
597  correct(scip, heap, state, &count, path, m, tail, e, cost[e], FSP_MODE);
598  }
599  if( nchecks++ > limit1 )
600  break;
601  }
602  g->mark[head] = TRUE;
603 
604  while( count > 0 && nchecks <= limit )
605  {
606  /* get nearest labelled node */
607  k = nearest(heap, state, &count, path);
608 
609  /* scanned */
610  state[k] = CONNECT;
611 
612  /* distance limit reached? */
613  if( SCIPisGT(scip, path[k].dist, distlimit) )
614  break;
615 
616  /* stop at terminals */
617  if( Is_term(g->term[k]) || k == head )
618  continue;
619 
620  /* correct incident nodes */
621  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
622  {
623  m = g->head[e];
624  if( state[m] && g->mark[m] && SCIPisGT(scip, path[m].dist, path[k].dist + cost[e]) )
625  {
626  /* m labelled for the first time? */
627  if( state[m] == UNKNOWN )
628  memlbl[(*nlbl)++] = m;
629  correct(scip, heap, state, &count, path, m, k, e, cost[e], FSP_MODE);
630  }
631  if( nchecks++ > limit )
632  break;
633  }
634  }
635 }
636 
637 
638 /** Dijkstra's algorithm starting from node 'start' */
639 void graph_path_execX(
640  SCIP* scip, /**< SCIP data structure */
641  const GRAPH* g, /**< graph data structure */
642  int start, /**< start vertex */
643  SCIP_Real* cost, /**< edgecosts */
644  SCIP_Real* pathdist, /**< distance array (on vertices) */
645  int* pathedge /**< predecessor edge array (on vertices) */
646  )
647 {
648  int k;
649  int m;
650  int i;
651  int count;
652  int* heap;
653  int* state;
654 
655  assert(g != NULL);
656  assert(start >= 0);
657  assert(start < g->knots);
658  assert(g->path_heap != NULL);
659  assert(g->path_state != NULL);
660  assert(pathdist != NULL);
661  assert(pathedge != NULL);
662  assert(cost != NULL);
663 
664  if (g->knots == 0)
665  return;
666 
667  heap = g->path_heap;
668  state = g->path_state;
669 
670  for( i = 0; i < g->knots; i++ )
671  {
672  state[i] = UNKNOWN;
673  pathdist[i] = FARAWAY;
674  pathedge[i] = -1;
675  }
676 
677  k = start;
678  pathdist[k] = 0.0;
679 
680  if( g->knots > 1 )
681  {
682  count = 1;
683  heap[count] = k;
684  state[k] = count;
685 
686  while( count > 0 )
687  {
688  k = nearestX(heap, state, &count, pathdist);
689 
690  state[k] = CONNECT;
691 #if 0
692  if( Is_term(g->term[k]) && k != start )
693  continue;
694 #endif
695  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
696  {
697  m = g->head[i];
698 
699  if( state[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[i])) )
700  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, i, cost[i]);
701  }
702  }
703  }
704 }
705 
706 /** Find a directed tree rooted in node 'start' and spanning all terminals */
707 void graph_path_st(
708  SCIP* scip, /**< SCIP data structure */
709  const GRAPH* g, /**< graph data structure */
710  SCIP_Real* cost, /**< edgecosts */
711  SCIP_Real* pathdist, /**< distance array (on vertices) */
712  int* pathedge, /**< predecessor edge array (on vertices) */
713  int start, /**< start vertex */
714  unsigned int* randseed, /**< random seed */
715  char* connected /**< array to mark whether a vertex is part of computed Steiner tree */
716  )
717 {
718  int z;
719  int k;
720  int m;
721  int e;
722  int count;
723  int nnodes;
724  int* heap;
725  int* state;
726  char cgt;
727 
728  assert(randseed != NULL);
729  assert(pathdist != NULL);
730  assert(pathedge != NULL);
731  assert(g != NULL);
732  assert(start >= 0);
733  assert(start < g->knots);
734  assert(cost != NULL);
735  assert(connected != NULL);
736 
737  count = 0;
738  nnodes = g->knots;
739  heap = g->path_heap;
740  state = g->path_state;
741 
742  /* initialize */
743  for( k = 0; k < nnodes; k++ )
744  {
745  state[k] = UNKNOWN;
746  pathdist[k] = FARAWAY;
747  pathedge[k] = -1;
748  connected[k] = FALSE;
749  }
750 
751  /* add start vertex to heap */
752  k = start;
753  pathdist[k] = 0.0;
754  connected[k] = TRUE;
755 
756  if( nnodes > 1 )
757  {
758  int node;
759  int nterms = 0;
760  int hnnodes = nnodes / 2;
761  z = SCIPgetRandomInt(0, nnodes - 1, randseed);
762  count = 1;
763  heap[count] = k;
764  state[k] = count;
765 
766  /* repeat until heap is empty */
767  while( count > 0 )
768  {
769  /* get closest node */
770  k = nearestX(heap, state, &count, pathdist);
771  state[k] = UNKNOWN;
772 
773  /* if k is terminal, connect its path to current subtree */
774  if( Is_term(g->term[k]) )
775  {
776  assert(k == start || !connected[k]);
777  connected[k] = TRUE;
778  pathdist[k] = 0.0;
779  node = k;
780  z = SCIPgetRandomInt(0, nnodes - 1, randseed);
781 
782  if( k != start )
783  {
784  assert(pathedge[k] != - 1);
785 
786  while( !connected[node = g->tail[pathedge[node]]] )
787  {
788  assert(pathedge[node] != - 1);
789  connected[node] = TRUE;
790  resetX(scip, pathdist, heap, state, &count, node);
791  }
792  }
793  /* have all terminals been reached? */
794  if( ++nterms == g->terms )
795  break;
796  }
797 
798  z = (k + z) % nnodes;
799  cgt = (z > hnnodes);
800 
801  /* update adjacent vertices */
802  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
803  {
804  m = g->head[e];
805 
806  /* is m not connected, allowed and closer (as close)? */
807  if( cgt || (state[m]) )
808  {
809  if( !connected[m] && g->mark[m] && SCIPisGT(scip, pathdist[m], (pathdist[k] + cost[e])) )
810  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
811  }
812  else
813  {
814  if( !connected[m] && g->mark[m] && SCIPisGE(scip, pathdist[m], (pathdist[k] + cost[e])) )
815  correctX(scip, heap, state, &count, pathdist, pathedge, m, k, e, cost[e]);
816  }
817  }
818  }
819  }
820 }
821 #if 0
822 /* computes the shortest path from each terminal to every other vertex */
823 void calculate_distances(
824  SCIP* scip, /**< SCIP data structure */
825  const GRAPH* g, /**< graph data structure */
826  PATH** path,
827  double* cost,
828  int mode
829  )
830 {
831  int i;
832 
833  assert(mode == FSP_MODE || mode == BSP_MODE);
834 
835  SCIPdebug(fputc('C', stdout));
836  SCIPdebug(fflush(stdout));
837 
838  for(i = 0; i < g->knots; i++)
839  {
840  if (Is_term(g->term[i]) && (g->grad[i] > 0))
841  {
842  if (path[i] == NULL)
843  path[i] = malloc((size_t)g->knots * sizeof(PATH));
844 
845  assert(path[i] != NULL);
846 
847  graph_path_exec(scip, g, mode, i, cost, path[i]);
848  }
849  else
850  {
851  if (path[i] != NULL)
852  {
853  free(path[i]);
854 
855  path[i] = NULL;
856  }
857  }
858  }
859 }
860 #endif
861 /** extend a voronoi region until all neighbouring terminals are spanned */
862 SCIP_RETCODE voronoi_extend(
863  SCIP* scip, /**< SCIP data structure */
864  const GRAPH* g, /**< graph data structure */
865  SCIP_Real* cost, /**< edgecosts */
866  PATH* path, /**< shortest paths data structure */
867  VLIST** adjterms, /**< structure for adjacent terminal data */
868  char* termsmark, /**< array to mark terminal */
869  int* reachednodes, /**< array to mark reached nodes */
870  int* nreachednodes, /**< pointer to number of reached nodes */
871  int* nodenterms, /**< array to store number of terimals to each node */
872  int nneighbterms, /**< number of neighbouring terminals */
873  int base, /**< voronoi base */
874  int countex /**< number of heap elements */
875  )
876 {
877  int* heap;
878  int* state;
879  int k;
880  int m;
881  int i;
882  int count = countex;
883  VLIST* curr;
884 
885  assert(g != NULL);
886  assert(g->path_heap != NULL);
887  assert(g->path_state != NULL);
888  assert(path != NULL);
889  assert(cost != NULL);
890 
891 
892  if( g->knots == 0 || nneighbterms <= 0 )
893  return SCIP_OKAY;
894 
895  heap = g->path_heap;
896  state = g->path_state;
897 
898 
899  if( g->knots > 1 )
900  {
901  /*
902  */
903  while( count > 0 && nneighbterms > 0 )
904  {
905 
906  k = nearest(heap, state, &count, path);
907 
908  state[k] = CONNECT;
909  reachednodes[(*nreachednodes)++] = k;
910  SCIP_CALL( SCIPallocMemory(scip, &curr) );
911  curr->dist = path[k].dist;
912  curr->edge = path[k].edge;
913  curr->base = base;
914  curr->next = adjterms[k];
915  adjterms[k] = curr;
916  nodenterms[k]++;
917 
918  if( termsmark[k] == TRUE )
919  {
920  termsmark[k] = FALSE;
921  if( --nneighbterms == 0 )
922  {
923  while( count > 0 )
924  {
925  reachednodes[(*nreachednodes)++] = heap[count--];
926  }
927  return SCIP_OKAY;
928  }
929 
930  }
931  else
932  {
933  /* iterate over all outgoing edges of vertex k */
934  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
935  {
936  m = g->head[i];
937 
938 
939  /* check whether the path (to m) including (k, m) is shorter than the so far best known */
940  if( (state[m]) && (g->mark[m]) && (GT(path[m].dist, (path[k].dist + cost[i]))) )
941  {
942  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
943  }
944  }
945  }
946  }
947  assert(nneighbterms == 0);
948  }
949  return SCIP_OKAY;
950 }
951 
952 
953 /** extend a voronoi region until all neighbouring terminals are spanned */
954 SCIP_RETCODE voronoi_extend2(
955  SCIP* scip, /**< SCIP data structure */
956  const GRAPH* g, /**< graph data structure */
957  SCIP_Real* cost, /**< edgecosts */
958  PATH* path, /**< shortest paths data structure */
959  SCIP_Real** distarr, /**< array to store distance from each node to its base */
960  int** basearr, /**< array to store the bases */
961  int** edgearr, /**< array to store the ancestor edge */
962  char* termsmark, /**< array to mark terminal */
963  int* reachednodes, /**< array to mark reached nodes */
964  int* nreachednodes, /**< pointer to number of reached nodes */
965  int* nodenterms, /**< array to store number of terimals to each node */
966  int nneighbterms, /**< number of neighbouring terminals */
967  int base, /**< voronoi base */
968  int countex /**< count of heap elements */
969  )
970 {
971  int k;
972  int m;
973  int i;
974  int* heap;
975  int* state;
976  int count = countex;
977 
978  assert(g != NULL);
979  assert(g->path_heap != NULL);
980  assert(g->path_state != NULL);
981  assert(path != NULL);
982  assert(cost != NULL);
983 
984  if( g->knots == 0 || nneighbterms <= 0 )
985  return SCIP_OKAY;
986 
987  heap = g->path_heap;
988  state = g->path_state;
989 
990  if( g->knots > 1 )
991  {
992  while( count > 0 && nneighbterms > 0 )
993  {
994  k = nearest(heap, state, &count, path);
995  state[k] = CONNECT;
996  reachednodes[(*nreachednodes)++] = k;
997  distarr[k][nodenterms[k]] = path[k].dist;
998  edgearr[k][nodenterms[k]] = path[k].edge;
999  basearr[k][nodenterms[k]] = base;
1000 
1001  nodenterms[k]++;
1002 
1003  if( termsmark[k] == TRUE )
1004  {
1005  termsmark[k] = FALSE;
1006  if( --nneighbterms == 0 )
1007  {
1008  while( count > 0 )
1009  reachednodes[(*nreachednodes)++] = heap[count--];
1010 
1011  return SCIP_OKAY;
1012  }
1013  }
1014  else if( !Is_term(g->term[k]) )
1015  {
1016  /* iterate over all outgoing edges of vertex k */
1017  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1018  {
1019  m = g->head[i];
1020 
1021  /* check whether the path (to m) including (k, m) is shorter than the so far best known */
1022  if( (state[m]) && (g->mark[m]) && (GT(path[m].dist, (path[k].dist + cost[i]))) )
1023  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
1024  }
1025  }
1026  }
1027  assert(nneighbterms == 0);
1028  }
1029  return SCIP_OKAY;
1030 }
1031 
1032 
1033 /** build a voronoi region, w.r.t. shortest paths, for a given set of bases */
1034 void voronoi(
1035  SCIP* scip, /**< SCIP data structure */
1036  const GRAPH* g, /**< graph data structure */
1037  SCIP_Real* cost, /**< edge costs */
1038  SCIP_Real* costrev, /**< reversed edge costs */
1039  char* base, /**< array to indicate whether a vertex is a Voronoi base */
1040  int* vbase, /**< voronoi base to each vertex */
1041  PATH* path /**< path data struture (leading to respective Voronoi base) */
1042  )
1043 {
1044  int k;
1045  int m;
1046  int i;
1047  int* heap;
1048  int* state;
1049  int count = 0;
1050  int root;
1051  int nbases = 0;
1052 
1053  assert(scip != NULL);
1054  assert(g != NULL);
1055  assert(g->path_heap != NULL);
1056  assert(g->path_state != NULL);
1057  assert(path != NULL);
1058  assert(cost != NULL);
1059  assert(costrev != NULL);
1060  if( g->knots == 0 )
1061  return;
1062 
1063  root = g->source[0];
1064  heap = g->path_heap;
1065  state = g->path_state;
1066 
1067  /* initialize */
1068  for( i = 0; i < g->knots; i++ )
1069  {
1070  /* set the base of vertex i */
1071  if( base[i] )
1072  {
1073  nbases++;
1074  if( g->knots > 1 )
1075  heap[++count] = i;
1076  vbase[i] = i;
1077  path[i].dist = 0.0;
1078  path[i].edge = UNKNOWN;
1079  state[i] = count;
1080  }
1081  else
1082  {
1083  vbase[i] = UNKNOWN;
1084  path[i].dist = FARAWAY + 1;
1085  path[i].edge = UNKNOWN;
1086  state[i] = UNKNOWN;
1087  }
1088  }
1089  assert(nbases > 0);
1090 
1091  if( g->knots > 1 )
1092  {
1093  /* until the heap is empty */
1094  while( count > 0 )
1095  {
1096  /* get the next (i.e. a nearest) vertex of the heap */
1097  k = nearest(heap, state, &count, path);
1098 
1099  /* mark vertex k as scanned */
1100  state[k] = CONNECT;
1101 
1102  /* iterate over all outgoing edges of vertex k */
1103  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1104  {
1105  m = g->head[i];
1106 
1107  /* check whether the path (to m) including k is shorter than the so far best known */
1108  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + ((vbase[k] == root)? cost[i] : costrev[i])) )
1109  {
1110  assert(!base[m]);
1111  correct(scip, heap, state, &count, path, m, k, i, ((vbase[k] == root)? cost[i] : costrev[i]), FSP_MODE);
1112  vbase[m] = vbase[k];
1113  }
1114  }
1115  }
1116  }
1117 }
1118 
1119 
1120 /** Voronoi Steiner tree extension for RPC, PC and MW */
1121 void voronoiSteinerTreeExt(
1122  SCIP* scip, /**< SCIP data structure */
1123  const GRAPH* g, /**< graph data structure */
1124  SCIP_Real* costrev, /**< reversed edge costs */
1125  int* vbase, /**< voronoi base to each vertex */
1126  char* stvertex, /**< array to indicate whether a vertex is a Voronoi base */
1127  PATH* path /**< path data struture (leading to respective Voronoi base) */
1128  )
1129 {
1130  int k;
1131  int m;
1132  int i;
1133  int count;
1134  int nnodes;
1135  int nbases;
1136  int* heap;
1137  int* state;
1138 
1139  assert(scip != NULL);
1140  assert(g != NULL);
1141  assert(g->path_heap != NULL);
1142  assert(g->path_state != NULL);
1143  assert(path != NULL);
1144  assert(costrev != NULL);
1145  assert(stvertex != NULL);
1146 
1147  heap = g->path_heap;
1148  count = 0;
1149  state = g->path_state;
1150  nbases = 0;
1151  nnodes = g->knots;
1152 
1153  /* initialize */
1154  for( i = 0; i < nnodes; i++ )
1155  {
1156  /* set the base of vertex i */
1157  if( !stvertex[i] && Is_pterm(g->term[i]) )
1158  {
1159  nbases++;
1160  if( nnodes > 1 )
1161  heap[++count] = i;
1162  vbase[i] = i;
1163  path[i].dist = -(g->prize[i]);
1164  path[i].edge = UNKNOWN;
1165  state[i] = count;
1166  }
1167  else
1168  {
1169  vbase[i] = UNKNOWN;
1170  path[i].dist = FARAWAY;
1171  path[i].edge = UNKNOWN;
1172  state[i] = UNKNOWN;
1173  }
1174  }
1175 
1176  if( nbases == 0)
1177  return;
1178 
1179  if( nnodes > 1 )
1180  {
1181  /* until the heap is empty */
1182  while( count > 0 )
1183  {
1184  /* get the next (i.e. a nearest) vertex of the heap */
1185  k = nearest(heap, state, &count, path);
1186 
1187  /* mark vertex k as scanned */
1188  state[k] = CONNECT;
1189 
1190  if( SCIPisGT(scip, path[k].dist, 0.0) )
1191  break;
1192 
1193  /* stop at Steiner tree */
1194  if( stvertex[k] )
1195  continue;
1196 
1197  /* iterate over all outgoing edges of vertex k */
1198  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1199  {
1200  m = g->head[i];
1201 
1202  /* check whether the path (to m) including k is shorter than the so far best known */
1203  if( vbase[m] != m && (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + costrev[i]) && !Is_term(g->term[m]) )
1204  {
1205  correct(scip, heap, state, &count, path, m, k, i, costrev[i], FSP_MODE);
1206  vbase[m] = vbase[k];
1207  }
1208  }
1209  }
1210  }
1211 }
1212 
1213 /** 2th next terminal to all non terminal nodes */
1214 void get2next(
1215  SCIP* scip, /**< SCIP data structure */
1216  const GRAPH* g, /**< graph data structure */
1217  SCIP_Real* cost, /**< edge costs */
1218  SCIP_Real* costrev, /**< reversed edge costs */
1219  PATH* path, /**< path data struture (leading to first and second nearest terminal) */
1220  int* vbase, /**< first and second nearest terminal to each non terminal */
1221  int* heap, /**< array representing a heap */
1222  int* state /**< array to mark the state of each node during calculation */
1223  )
1224 {
1225  int k;
1226  int j;
1227  int i;
1228  int e;
1229  int root;
1230  int count;
1231  int nnodes;
1232 
1233  assert(g != NULL);
1234  assert(path != NULL);
1235  assert(cost != NULL);
1236  assert(heap != NULL);
1237  assert(state != NULL);
1238  assert(costrev != NULL);
1239 
1240  root = g->source[0];
1241  count = 0;
1242  nnodes = g->knots;
1243 
1244  /* initialize */
1245  for( i = 0; i < nnodes; i++ )
1246  {
1247  state[i] = CONNECT;
1248 
1249  /* copy of node i */
1250  k = i + nnodes;
1251  vbase[k] = UNKNOWN;
1252  state[k] = UNKNOWN;
1253  path[k].edge = UNKNOWN;
1254  path[k].dist = FARAWAY;
1255  }
1256 
1257  /* scan original nodes */
1258  for( i = 0; i < nnodes; i++ )
1259  {
1260  if( !g->mark[i] )
1261  continue;
1262 
1263  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
1264  {
1265  j = g->head[e];
1266  k = j + nnodes;
1267  if( !Is_term(g->term[j]) && SCIPisGT(scip, path[k].dist, path[i].dist + ((root == vbase[i])? cost[e] : costrev[e])) &&
1268  vbase[i] != vbase[j] && g->mark[j] )
1269  {
1270  correct(scip, heap, state, &count, path, k, i, e, ((root == vbase[i])? cost[e] : costrev[e]), FSP_MODE);
1271  vbase[k] = vbase[i];
1272  }
1273  }
1274  }
1275 
1276  if( nnodes > 1 )
1277  {
1278  int jc;
1279  /* until the heap is empty */
1280  while( count > 0 )
1281  {
1282  /* get the next (i.e. a nearest) vertex of the heap */
1283  k = nearest(heap, state, &count, path);
1284 
1285  /* mark vertex k as removed from heap */
1286  state[k] = UNKNOWN;
1287 
1288  assert(k - nnodes >= 0);
1289  /* iterate over all outgoing edges of vertex k */
1290  for( e = g->outbeg[k - nnodes]; e != EAT_LAST; e = g->oeat[e] )
1291  {
1292  j = g->head[e];
1293  if( Is_term(g->term[j]) || !g->mark[j] )
1294  continue;
1295  jc = j + nnodes;
1296 
1297  /* check whether the path (to j) including k is shorter than the so far best known */
1298  if( vbase[j] != vbase[k] && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) )
1299  {
1300  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
1301  vbase[jc] = vbase[k];
1302  }
1303  }
1304  }
1305  }
1306  return;
1307 }
1308 
1309 /* 3th next terminal to all non terminal nodes */
1310 void get3next(
1311  SCIP* scip, /**< SCIP data structure */
1312  const GRAPH* g, /**< graph data structure */
1313  SCIP_Real* cost, /**< edge costs */
1314  SCIP_Real* costrev, /**< reversed edge costs */
1315  PATH* path, /**< path data struture (leading to first, second and third nearest terminal) */
1316  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1317  int* heap, /**< array representing a heap */
1318  int* state /**< array to mark the state of each node during calculation */
1319  )
1320 {
1321  int k;
1322  int j;
1323  int i;
1324  int l;
1325  int v;
1326  int e;
1327  int root;
1328  int count = 0;
1329  int nnodes;
1330  int dnnodes;
1331 
1332  assert(g != NULL);
1333  assert(path != NULL);
1334  assert(cost != NULL);
1335  assert(heap != NULL);
1336  assert(state != NULL);
1337  assert(costrev != NULL);
1338 
1339  root = g->source[0];
1340  nnodes = g->knots;
1341  dnnodes = 2 * nnodes;
1342 
1343  /* initialize */
1344  for( i = 0; i < nnodes; i++ )
1345  {
1346  /* copy of node i */
1347  k = i + dnnodes;
1348  vbase[k] = UNKNOWN;
1349  state[k] = UNKNOWN;
1350  path[k].edge = UNKNOWN;
1351  path[k].dist = FARAWAY;
1352  }
1353 
1354  /* scan original nodes */
1355  for( i = 0; i < nnodes; i++ )
1356  {
1357  state[i] = CONNECT;
1358  state[i + nnodes] = CONNECT;
1359  if( !g->mark[i] )
1360  continue;
1361 
1362  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
1363  {
1364  j = g->head[e];
1365  k = j + dnnodes;
1366  if( !Is_term(g->term[j]) && g->mark[j] )
1367  {
1368  v = i;
1369  for( l = 0; l < 2; l++ )
1370  {
1371  if( SCIPisGT(scip, path[k].dist, path[v].dist + ((root == vbase[v])? cost[e] : costrev[e])) &&
1372  vbase[v] != vbase[j] && vbase[v] != vbase[j + nnodes] )
1373  {
1374  correct(scip, heap, state, &count, path, k, v, e, ((root == vbase[v])? cost[e] : costrev[e]), FSP_MODE);
1375  vbase[k] = vbase[v];
1376  }
1377  v += nnodes;
1378  }
1379  }
1380  }
1381  }
1382  if( nnodes > 1 )
1383  {
1384  int jc;
1385  /* until the heap is empty */
1386  while( count > 0 )
1387  {
1388  /* get the next (i.e. a nearest) vertex of the heap */
1389  k = nearest(heap, state, &count, path);
1390 
1391  /* mark vertex k as removed from heap */
1392  state[k] = UNKNOWN;
1393 
1394  assert(k - dnnodes >= 0);
1395  /* iterate over all outgoing edges of vertex k */
1396  for( e = g->outbeg[k - dnnodes]; e != EAT_LAST; e = g->oeat[e] )
1397  {
1398  j = g->head[e];
1399  if( Is_term(g->term[j]) || !g->mark[j] )
1400  continue;
1401  jc = j + dnnodes;
1402  /* check whether the path (to j) including k is shorter than the so far best known */
1403  if( vbase[j] != vbase[k] && vbase[j + nnodes] != vbase[k]
1404  && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) ) /*TODO(state[jc])??*/
1405  {
1406  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
1407  vbase[jc] = vbase[k];
1408  }
1409  }
1410  }
1411  }
1412  return;
1413 }
1414 
1415 
1416 /* 4th next terminal to all non terminal nodes */
1417 void get4next(
1418  SCIP* scip, /**< SCIP data structure */
1419  const GRAPH* g, /**< graph data structure */
1420  SCIP_Real* cost, /**< edge costs */
1421  SCIP_Real* costrev, /**< reversed edge costs */
1422  PATH* path, /**< path data struture (leading to first, second and third nearest terminal) */
1423  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1424  int* heap, /**< array representing a heap */
1425  int* state /**< array to mark the state of each node during calculation */
1426  )
1427 {
1428  int k;
1429  int j;
1430  int i;
1431  int l;
1432  int v;
1433  int e;
1434  int root;
1435  int count = 0;
1436  int nnodes;
1437  int dnnodes;
1438  int tnnodes;
1439 
1440  assert(g != NULL);
1441  assert(path != NULL);
1442  assert(cost != NULL);
1443  assert(heap != NULL);
1444  assert(state != NULL);
1445  assert(costrev != NULL);
1446 
1447  root = g->source[0];
1448  nnodes = g->knots;
1449  dnnodes = 2 * nnodes;
1450  tnnodes = 3 * nnodes;
1451 
1452  /* initialize */
1453  for( i = 0; i < nnodes; i++ )
1454  {
1455  /* copy of node i */
1456  k = i + tnnodes;
1457  vbase[k] = UNKNOWN;
1458  state[k] = UNKNOWN;
1459  path[k].edge = UNKNOWN;
1460  path[k].dist = FARAWAY;
1461  }
1462 
1463  /* scan original nodes */
1464  for( i = 0; i < nnodes; i++ )
1465  {
1466  state[i] = CONNECT;
1467  state[i + nnodes] = CONNECT;
1468  state[i + dnnodes] = CONNECT;
1469  if( !g->mark[i] )
1470  continue;
1471 
1472  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
1473  {
1474  j = g->head[e];
1475  k = j + tnnodes;
1476  if( !Is_term(g->term[j]) && g->mark[j] )
1477  {
1478  v = i;
1479  for( l = 0; l < 3; l++ )
1480  {
1481  if( SCIPisGT(scip, path[k].dist, path[v].dist + ((root == vbase[v])? cost[e] : costrev[e])) &&
1482  vbase[v] != vbase[j] && vbase[v] != vbase[j + nnodes] && vbase[v] != vbase[j + dnnodes] )
1483  {
1484  correct(scip, heap, state, &count, path, k, v, e, ((root == vbase[v])? cost[e] : costrev[e]), FSP_MODE);
1485  vbase[k] = vbase[v];
1486  }
1487  v += nnodes;
1488  }
1489  }
1490  }
1491  }
1492  if( nnodes > 1 )
1493  {
1494  int jc;
1495  /* until the heap is empty */
1496  while( count > 0 )
1497  {
1498  /* get the next (i.e. a nearest) vertex of the heap */
1499  k = nearest(heap, state, &count, path);
1500 
1501  /* mark vertex k as removed from heap */
1502  state[k] = UNKNOWN;
1503 
1504  assert(k - tnnodes >= 0);
1505  /* iterate over all outgoing edges of vertex k */
1506  for( e = g->outbeg[k - tnnodes]; e != EAT_LAST; e = g->oeat[e] )
1507  {
1508  j = g->head[e];
1509  if( Is_term(g->term[j]) || !g->mark[j] )
1510  continue;
1511  jc = j + tnnodes;
1512  /* check whether the path (to j) including k is shorter than the so far best known */
1513  if( vbase[j] != vbase[k] && vbase[j + nnodes] != vbase[k] && vbase[j + dnnodes] != vbase[k]
1514  && SCIPisGT(scip, path[jc].dist, path[k].dist + ((root == vbase[k])? cost[e] : costrev[e])) ) /*TODO(state[jc])??*/
1515  {
1516  correct(scip, heap, state, &count, path, jc, k, e, (root == vbase[k])? cost[e] : costrev[e], FSP_MODE);
1517  vbase[jc] = vbase[k];
1518  }
1519  }
1520  }
1521  }
1522  return;
1523 }
1524 
1525 /** build a voronoi region in presolving, w.r.t. shortest paths, for all terminals*/
1526 void getnext3terms(
1527  SCIP* scip, /**< SCIP data structure */
1528  const GRAPH* g, /**< graph data structure */
1529  SCIP_Real* cost, /**< edge costs */
1530  SCIP_Real* costrev, /**< reversed edge costs */
1531  PATH* path3, /**< path data struture (leading to first, second and third nearest terminal) */
1532  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1533  int* heap, /**< array representing a heap */
1534  int* state /**< array to mark the state of each node during calculation */
1535  )
1536 {
1537  int k;
1538  assert(g != NULL);
1539  assert(path3 != NULL);
1540  assert(cost != NULL);
1541  assert(costrev != NULL);
1542  assert(heap != NULL);
1543  assert(state != NULL);
1544 
1545  if( g->stp_type != STP_PRIZE_COLLECTING && g->stp_type != STP_ROOTED_PRIZE_COLLECTING )
1546  for( k = 0; k < g->knots; k++ )
1547  g->mark[k] = (g->grad[k] > 0);
1548 
1549  /* build voronoi diagram */
1550  voronoi_terms(scip, g, cost, path3, vbase, heap, state);
1551 
1552  /* get 2nd nearest terms */
1553  get2next(scip, g, cost, costrev, path3, vbase, heap, state);
1554 
1555  /* get 3th nearest terms */
1556  get3next(scip, g, cost, costrev, path3, vbase, heap, state);
1557 
1558  return;
1559 }
1560 
1561 /** build a voronoi region in presolving, w.r.t. shortest paths, for all terminals*/
1562 void getnext4terms(
1563  SCIP* scip, /**< SCIP data structure */
1564  const GRAPH* g, /**< graph data structure */
1565  SCIP_Real* cost, /**< edge costs */
1566  SCIP_Real* costrev, /**< reversed edge costs */
1567  PATH* path, /**< path data struture (leading to first, second, third and fouth nearest terminal) */
1568  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1569  int* heap, /**< array representing a heap */
1570  int* state /**< array to mark the state of each node during calculation */
1571  )
1572 {
1573  int k;
1574 
1575  assert(g != NULL);
1576  assert(path != NULL);
1577  assert(cost != NULL);
1578  assert(heap != NULL);
1579  assert(state != NULL);
1580  assert(costrev != NULL);
1581 
1582  if( g->stp_type != STP_PRIZE_COLLECTING && g->stp_type != STP_ROOTED_PRIZE_COLLECTING )
1583  for( k = 0; k < g->knots; k++ )
1584  g->mark[k] = (g->grad[k] > 0);
1585 
1586  /* build voronoi diagram */
1587  voronoi_terms(scip, g, cost, path, vbase, heap, state);
1588 
1589  /* get 2nd nearest terms */
1590  get2next(scip, g, cost, costrev, path, vbase, heap, state);
1591 
1592  /* get 3th nearest terms */
1593  get3next(scip, g, cost, costrev, path, vbase, heap, state);
1594 
1595  /* get 4th nearest terms */
1596  get4next(scip, g, cost, costrev, path, vbase, heap, state);
1597 
1598  return;
1599 }
1600 
1601 
1602 /** get 4 close terminals to each terminal */
1603 void getnext4tterms(
1604  SCIP* scip, /**< SCIP data structure */
1605  const GRAPH* g, /**< graph data structure */
1606  SCIP_Real* cost, /**< edge costs */
1607  SCIP_Real* boundedges, /**< array to store boundary edges */
1608  PATH* path, /**< path data struture (leading to first, second, third and fouth nearest terminal) */
1609  int* vbase, /**< first, second and third nearest terminal to each non terminal */
1610  int* heap, /**< array representing a heap */
1611  int* state /**< array to mark the state of each node during calculation */
1612  )
1613 {
1614  int k;
1615  int e;
1616  int l;
1617  int k2;
1618  int bedge;
1619  int shift;
1620  int nnodes;
1621  int nboundedges;
1622 
1623  assert(g != NULL);
1624  assert(path != NULL);
1625  assert(cost != NULL);
1626  assert(heap != NULL);
1627  assert(state != NULL);
1628  assert(boundedges != NULL);
1629 
1630  shift = 0;
1631  nnodes = g->knots;
1632 
1633  if( g->stp_type != STP_PRIZE_COLLECTING && g->stp_type != STP_ROOTED_PRIZE_COLLECTING )
1634  for( k = 0; k < g->knots; k++ )
1635  g->mark[k] = (g->grad[k] > 0);
1636 
1637  nboundedges = 0;
1638  for( k = 0; k < nnodes; k ++ )
1639  {
1640  if( !g->mark[k] )
1641  continue;
1642 
1643  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1644  {
1645  k2 = g->head[e];
1646  if( !g->mark[k2] || k2 < k )
1647  continue;
1648 
1649  /* is e a boundary edge? */
1650  if( vbase[k] != vbase[k2] )
1651  {
1652  boundedges[nboundedges++] = e;
1653  }
1654  }
1655  if( Is_term(g->term[k]) )
1656  {
1657  path[k].dist = FARAWAY;
1658  vbase[k] = UNKNOWN;
1659  }
1660 
1661  }
1662 
1663  for( l = 0; l < 4; l++ )
1664  {
1665  for( e = 0; e < nboundedges; e++ )
1666  {
1667  bedge = (int) boundedges[e];
1668  k = g->tail[bedge];
1669  k2 = g->head[bedge];
1670  utdist(scip, g, path, cost[bedge], vbase, k, l, k2, shift, nnodes);
1671  utdist(scip, g, path, cost[bedge], vbase, k2, l, k, shift, nnodes);
1672  }
1673  shift += nnodes;
1674  }
1675  return;
1676 }
1677 
1678 /** build a Voronoi region in presolving, w.r.t. shortest paths, for all terminals */
1679 void voronoi_terms(
1680  SCIP* scip, /**< SCIP data structure */
1681  const GRAPH* g, /**< graph data structure */
1682  SCIP_Real* cost, /**< edge costs */
1683  PATH* path, /**< path data struture (leading to respective Voronoi base) */
1684  int* vbase, /**< Voronoi base to each vertex */
1685  int* heap, /**< array representing a heap */
1686  int* state /**< array to mark the state of each node during calculation */
1687  )
1688 {
1689  int k;
1690  int m;
1691  int i;
1692  int count = 0;
1693  int nbases = 0;
1694 
1695  assert(g != NULL);
1696  assert(path != NULL);
1697  assert(cost != NULL);
1698  assert(heap != NULL);
1699  assert(state != NULL);
1700 
1701  /* initialize */
1702  for( i = 0; i < g->knots; i++ )
1703  {
1704  /* set the base of vertex i */
1705  if( Is_term(g->term[i]) && g->mark[i] )
1706  {
1707  nbases++;
1708  if( g->knots > 1 )
1709  heap[++count] = i;
1710  vbase[i] = i;
1711  path[i].dist = 0.0;
1712  path[i].edge = UNKNOWN;
1713  state[i] = count;
1714  }
1715  else
1716  {
1717  vbase[i] = UNKNOWN;
1718  path[i].dist = FARAWAY;
1719  path[i].edge = UNKNOWN;
1720  state[i] = UNKNOWN;
1721  }
1722  }
1723  if( nbases == 0 )
1724  return;
1725  if( g->knots > 1 )
1726  {
1727  /* until the heap is empty */
1728  while( count > 0 )
1729  {
1730  /* get the next (i.e. a nearest) vertex of the heap */
1731  k = nearest(heap, state, &count, path);
1732 
1733  /* mark vertex k as scanned */
1734  state[k] = CONNECT;
1735 
1736  /* iterate over all outgoing edges of vertex k */
1737  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1738  {
1739  m = g->head[i];
1740 
1741  /* check whether the path (to m) including k is shorter than the so far best known */
1742  if( (state[m]) && SCIPisGT(scip, path[m].dist, path[k].dist + cost[i]) && g->mark[m] )
1743  {
1744  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
1745  vbase[m] = vbase[k];
1746  }
1747  }
1748  }
1749  }
1750  return;
1751 }
1752 
1753 
1754 /** build a voronoi region, w.r.t. shortest paths, for all terminal and the distance */
1755 SCIP_RETCODE voronoi_dist(
1756  SCIP* scip, /**< SCIP data structure */
1757  const GRAPH* g, /**< graph data structure */
1758  SCIP_Real* cost, /**< edge costs */
1759  SCIP_Real* distance, /**< array storing path from a terminal over shortest
1760  incident edge to nearest terminal */
1761  int* minedgepred, /**< shortest edge predecessor array */
1762  int* vbase, /**< array containing Voronoi base to each node */
1763  int* minarc, /**< array to mark whether an edge is one a path corresponding to 'distance' */
1764  int* heap, /**< array representing a heap */
1765  int* state, /**< array indicating state of each vertex during calculation of Voronoi regions */
1766  int* distnode, /**< array to store terminal corresponding to distance stored in distance array */
1767  PATH* path /**< array containing Voronoi paths data */
1768  )
1769 {
1770  SCIP_Real new;
1771  int e;
1772  int k;
1773  int m;
1774  int i;
1775  int pred;
1776  int count;
1777  int nbases;
1778  int nnodes;
1779  int prededge;
1780 
1781  assert(g != NULL);
1782  assert(path != NULL);
1783  assert(cost != NULL);
1784  assert(heap != NULL);
1785  assert(state != NULL);
1786  assert(distance != NULL);
1787 
1788  count = 0;
1789  nbases = 0;
1790  nnodes = g->knots;
1791 
1792  /* initialize */
1793  for( i = 0; i < nnodes; i++ )
1794  {
1795  distance[i] = FARAWAY;
1796  if( distnode != NULL )
1797  distnode[i] = UNKNOWN;
1798 
1799  /* set the base of vertex i */
1800  if( Is_term(g->term[i]) && g->mark[i] )
1801  {
1802  nbases++;
1803  if( nnodes > 1 )
1804  heap[++count] = i;
1805  vbase[i] = i;
1806  state[i] = count;
1807  path[i].dist = 0.0;
1808  }
1809  else
1810  {
1811  vbase[i] = UNKNOWN;
1812  state[i] = UNKNOWN;
1813  path[i].dist = FARAWAY;
1814  }
1815  path[i].edge = UNKNOWN;
1816  }
1817 
1818  /* initialize predecessor array */
1819 
1820  for( e = 0; e < g->edges; e++ )
1821  minedgepred[e] = FALSE;
1822 
1823  for( k = 0; k < nbases; k++ )
1824  if( minarc[k] != UNKNOWN )
1825  minedgepred[minarc[k]] = TRUE;
1826 
1827  /* if graph contains more than one vertices, start main loop */
1828  if( nnodes > 1 )
1829  {
1830  /* until the heap is empty */
1831  while( count > 0 )
1832  {
1833  /* get the next (i.e. a nearest) vertex of the heap */
1834  k = nearest(heap, state, &count, path);
1835 
1836  /* mark vertex k as scanned */
1837  state[k] = CONNECT;
1838 
1839  prededge = path[k].edge;
1840 
1841  if( prededge != UNKNOWN )
1842  {
1843  pred = g->tail[prededge];
1844 
1845  assert(vbase[k] != UNKNOWN);
1846  assert(vbase[pred] != UNKNOWN);
1847  assert(vbase[pred] == vbase[k]);
1848  assert(g->head[prededge] == k);
1849 
1850  if( !Is_term(g->term[pred]) && g->mark[pred] )
1851  {
1852  assert(path[pred].edge != UNKNOWN);
1853  minedgepred[prededge] = minedgepred[path[pred].edge];
1854  }
1855 
1856  }
1857 
1858  /* iterate over all outgoing edges of vertex k */
1859  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
1860  {
1861  m = g->head[i];
1862 
1863  if( state[m] == CONNECT && vbase[m] != vbase[k] && g->mark[m] )
1864  {
1865  if( minedgepred[i] || (prededge != UNKNOWN && minedgepred[prededge] ) )
1866  {
1867  new = path[k].dist + g->cost[i] + path[m].dist;
1868  if( SCIPisLT(scip, new, distance[vbase[k]]) )
1869  {
1870  if( distnode != NULL )
1871  distnode[vbase[k]] = vbase[m];
1872 #if 0
1873  if( vbase[k] == 0 )
1874  printf("update dist to %f k: %d m: %d \n", new , k, m);
1875 #endif
1876  distance[vbase[k]] = new;
1877  }
1878  }
1879  if( minedgepred[Edge_anti(i)] || (path[m].edge != UNKNOWN && minedgepred[path[m].edge] ) )
1880  {
1881  new = path[m].dist + g->cost[i] + path[k].dist;
1882  if( SCIPisLT(scip, new, distance[vbase[m]]) )
1883  {
1884  if( distnode != NULL )
1885  distnode[vbase[m]] = vbase[k];
1886 #if 0
1887  if( vbase[m] == 0 )
1888  printf("update dist to %f k: %d m: %d \n", new , k, m);
1889 #endif
1890  distance[vbase[m]] = new;
1891  }
1892  }
1893  }
1894 
1895  /* Check whether the path (to m) including k is shorter than the so far best known.
1896  In case of equality, also update if k is sucessor of minedge. */
1897  if( state[m] && g->mark[m] &&
1898  (SCIPisGT(scip, path[m].dist, path[k].dist + cost[i]) ||
1899  (prededge != UNKNOWN && minedgepred[prededge] && SCIPisEQ(scip, path[m].dist, path[k].dist + cost[i]))) )
1900  {
1901  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
1902  vbase[m] = vbase[k];
1903  }
1904  }
1905  }
1906  }
1907 
1908  return SCIP_OKAY;
1909 }
1910 
1911 /** build voronoi regions, w.r.t. shortest paths, for all terminals and compute the radii */
1912 SCIP_RETCODE voronoi_radius(
1913  SCIP* scip, /**< SCIP data structure */
1914  const GRAPH* graph, /**< graph data structure */
1915  GRAPH* adjgraph, /**< graph data structure */
1916  PATH* path, /**< array containing Voronoi paths data */
1917  SCIP_Real* rad, /**< array storing shortest way from a terminal out of its Voronoi region */
1918  SCIP_Real* cost, /**< edge costs */
1919  SCIP_Real* costrev, /**< reversed edge costs */
1920  int* vbase, /**< array containing Voronoi base of each node */
1921  int* heap, /**< array representing a heap */
1922  int* state /**< array to mark state of each node during calculation */
1923  )
1924 {
1925  int* nodesid;
1926  int i;
1927  int k;
1928  int m;
1929  int vbm;
1930  int vbk;
1931  int root;
1932  int count = 0;
1933  int nnodes;
1934  int nterms = 0;
1935  char pc;
1936  char mw;
1937 
1938  assert(graph != NULL);
1939  assert(heap != NULL);
1940  assert(state != NULL);
1941  assert(path != NULL);
1942  assert(cost != NULL);
1943  assert(costrev != NULL);
1944  assert(rad != NULL);
1945  assert(vbase != NULL);
1946 
1947  nnodes = graph->knots;
1948  if( nnodes == 0 || graph->terms == 0 )
1949  return SCIP_OKAY;
1950  root = graph->source[0];
1951  mw = (graph->stp_type == STP_MAX_NODE_WEIGHT);
1952  pc = ((graph->stp_type == STP_PRIZE_COLLECTING) || (graph->stp_type == STP_ROOTED_PRIZE_COLLECTING));
1953  SCIP_CALL( SCIPallocBufferArray(scip, &nodesid, nnodes) );
1954 
1955  /* initialize */
1956  for( i = 0; i < nnodes; i++ )
1957  {
1958  rad[i] = FARAWAY;
1959 
1960  /* set the base of vertex i */
1961  if( Is_term(graph->term[i]) && graph->mark[i] )
1962  {
1963  if( nnodes > 1 )
1964  heap[++count] = i;
1965 
1966  if( !mw )
1967  {
1968  adjgraph->mark[adjgraph->knots] = TRUE;
1969  graph_knot_add(adjgraph, -1);
1970  }
1971 
1972  nodesid[i] = nterms++;
1973  vbase[i] = i;
1974 
1975  if( mw )
1976  assert(SCIPisGE(scip, graph->prize[i], 0.0));
1977  if( mw )
1978  path[i].dist = -graph->prize[i];
1979  else
1980  path[i].dist = 0.0;
1981 
1982  path[i].edge = UNKNOWN;
1983  state[i] = count;
1984  }
1985  else
1986  {
1987  vbase[i] = UNKNOWN;
1988  nodesid[i] = UNKNOWN;
1989  path[i].dist = FARAWAY;
1990  path[i].edge = UNKNOWN;
1991  state[i] = UNKNOWN;
1992  }
1993  }
1994 
1995  if( nnodes > 1 )
1996  {
1997  SCIP_Real c1;
1998  SCIP_Real c2;
1999  SCIP_Real ecost;
2000  int ne;
2001 
2002  /* until the heap is empty */
2003  while( count > 0 )
2004  {
2005  /* get the next (i.e. a nearest) vertex of the heap */
2006  k = nearest(heap, state, &count, path);
2007 
2008  /* mark vertex k as scanned */
2009  state[k] = CONNECT;
2010 
2011  /* iterate over all outgoing edges of vertex k */
2012  for( i = graph->outbeg[k]; i != EAT_LAST; i = graph->oeat[i] )
2013  {
2014  m = graph->head[i];
2015  vbm = vbase[m];
2016  vbk = vbase[k];
2017 
2018  if( state[m] == CONNECT && vbm != vbk && graph->mark[m] )
2019  {
2020  assert(graph->mark[vbm]);
2021  assert(graph->mark[vbk]);
2022  if( mw )
2023  {
2024  if( SCIPisGT(scip, rad[vbk], path[k].dist) )
2025  rad[vbk] = path[k].dist;
2026  if( SCIPisGT(scip, rad[vbm], path[m].dist) )
2027  rad[vbm] = path[m].dist;
2028 #if 0
2029  if( SCIPisGT(scip, path[m].dist + graph->prize[vbm], path[k].dist + graph->prize[vbk]) )
2030  ecost = graph->prize[vbm] - path[k].dist;
2031  else
2032  ecost = graph->prize[vbk] - path[m].dist;
2033  if( SCIPisLT(scip, ecost, 0.0) )
2034  ecost = 0.0;
2035  /* find edge in adjgraph */
2036  for( ne = adjgraph->outbeg[nodesid[vbk]]; ne != EAT_LAST; ne = adjgraph->oeat[ne] )
2037  if( adjgraph->head[ne] == nodesid[vbm] )
2038  break;
2039 
2040  /* edge exists? */
2041  if( ne != EAT_LAST )
2042  {
2043  assert(ne >= 0);
2044  assert(adjgraph->head[ne] == nodesid[vbm]);
2045  assert(adjgraph->tail[ne] == nodesid[vbk]);
2046  if( SCIPisLT(scip, adjgraph->cost[ne], ecost) )
2047  {
2048  adjgraph->cost[ne] = ecost;
2049  adjgraph->cost[Edge_anti(ne)] = ecost;
2050  }
2051  }
2052  else
2053  {
2054  graph_edge_add(scip, adjgraph, nodesid[vbm], nodesid[vbk], ecost, ecost);
2055  }
2056 #endif
2057  }
2058  else
2059  {
2060  if( graph->stp_type == STP_HOP_CONS )
2061  {
2062  if( m == root )
2063  c1 = path[m].dist + costrev[i];
2064  else
2065  c1 = path[m].dist + cost[i];
2066  if( k == root )
2067  c2 = path[k].dist + cost[i];
2068  else
2069  c2 = path[k].dist + costrev[i];
2070 
2071  if( SCIPisGT(scip, c1, c2) )
2072  ecost = c2;
2073  else
2074  ecost = c1;
2075  }
2076  else
2077  {
2078  if( SCIPisGT(scip, path[m].dist, path[k].dist) )
2079  ecost = path[k].dist + cost[i];
2080  else
2081  ecost = path[m].dist + cost[i];
2082  }
2083 
2084  if( pc )
2085  {
2086  if( SCIPisGT(scip, ecost, graph->prize[vbm]) && root != vbm )
2087  ecost = graph->prize[vbm];
2088  if( SCIPisGT(scip, ecost, graph->prize[vbk]) && root != vbk )
2089  ecost = graph->prize[vbk];
2090  }
2091 
2092  /* find edge in adjgraph */
2093  for( ne = adjgraph->outbeg[nodesid[vbk]]; ne != EAT_LAST; ne = adjgraph->oeat[ne] )
2094  if( adjgraph->head[ne] == nodesid[vbm] )
2095  break;
2096 
2097  /* edge exists? */
2098  if( ne != EAT_LAST )
2099  {
2100  assert(ne >= 0);
2101  assert(adjgraph->head[ne] == nodesid[vbm]);
2102  assert(adjgraph->tail[ne] == nodesid[vbk]);
2103  if( SCIPisGT(scip, adjgraph->cost[ne], ecost) )
2104  {
2105  adjgraph->cost[ne] = ecost;
2106  adjgraph->cost[Edge_anti(ne)] = ecost;
2107  }
2108  }
2109  else
2110  {
2111  graph_edge_add(scip, adjgraph, nodesid[vbm], nodesid[vbk], ecost, ecost);
2112  }
2113 
2114  if( SCIPisGT(scip, rad[vbk], path[k].dist + ((vbk == root)? cost[i] : costrev[i])) )
2115  rad[vbk] = path[k].dist + ((vbk == root)? cost[i] : costrev[i]);
2116 
2117  if( SCIPisGT(scip, rad[vbm], path[m].dist + ((vbm == root)? costrev[i] : cost[i])) )
2118  rad[vbm] = path[m].dist + ((vbm == root)? costrev[i] : cost[i]);
2119  }
2120  }
2121 
2122  /* check whether the path (to m) including k is shorter than the so far best known */
2123  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])) )
2124  {
2125  correct(scip, heap, state, &count, path, m, k, i, ((vbk == root)? cost[i] : costrev[i]), FSP_MODE);
2126  vbase[m] = vbk;
2127  }
2128  }
2129  }
2130  }
2131  SCIPfreeBufferArray(scip, &nodesid);
2132  return SCIP_OKAY;
2133 }
2134 
2135 /** repair the voronoi diagram for SL reduction */
2136 void voronoi_slrepair(
2137  SCIP* scip, /**< SCIP data structure */
2138  const GRAPH* g, /**< graph data structure */
2139  SCIP_Real* cost, /**< edge costs */
2140  PATH* path, /**< array to store */
2141  int* vbase, /**< array containing Voronoi base of each node*/
2142  int* heap, /**< array representing a heap */
2143  int* state, /**< array indicating state of each vertex during calculation of Voronoi regions */
2144  int start, /**< node to start repair from */
2145  int adjnode /**< adjacent node */
2146  )
2147 {
2148  int k;
2149  int m;
2150  int i;
2151  int count = 0;
2152 
2153  assert(g != NULL);
2154  assert(path != NULL);
2155  assert(cost != NULL);
2156  assert(heap != NULL);
2157  assert(state != NULL);
2158 
2159  /* initialize */
2160  for( i = 0; i < g->knots; i++ )
2161  state[i] = UNKNOWN;
2162 
2163  if( SCIPisGT(scip, path[start].dist, path[adjnode].dist) )
2164  {
2165  if( vbase[adjnode] == adjnode )
2166  {
2167  assert(Is_term(g->term[start]));
2168  vbase[start] = start;
2169  }
2170  else
2171  {
2172  vbase[start] = vbase[adjnode];
2173  }
2174  path[start].dist = path[adjnode].dist;
2175  path[start].edge = path[adjnode].edge;
2176  }
2177  k = start;
2178 
2179  if( g->knots > 1 )
2180  {
2181  count = 1;
2182  heap[count] = k;
2183  state[k] = count;
2184 
2185  /* until the heap is empty */
2186  while( count > 0 )
2187  {
2188  /* get the next (i.e. a nearest) vertex of the heap */
2189  k = nearest(heap, state, &count, path);
2190 
2191  /* mark vertex k as scanned */
2192  state[k] = CONNECT;
2193 
2194  /* iterate over all outgoing edges of vertex k */
2195  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2196  {
2197  m = g->head[i];
2198  /* check whether the path (to m) including k is shorter than the so far best known */
2199  if( (state[m] != CONNECT) && SCIPisGE(scip, path[m].dist, path[k].dist + cost[i]) && g->mark[m] )
2200  {
2201  correct(scip, heap, state, &count, path, m, k, i, cost[i], FSP_MODE);
2202  vbase[m] = vbase[k];
2203  }
2204  }
2205  }
2206  }
2207  return;
2208 }
2209 
2210 /** repair the voronoi diagram for a given set nodes */
2211 void voronoi_repair(
2212  SCIP* scip, /**< SCIP data structure */
2213  const GRAPH* g, /**< graph data structure */
2214  SCIP_Real* cost, /**< edge costs */
2215  int* count, /**< pointer to number of heap elements */
2216  int* vbase, /**< array containing Voronoi base of each node */
2217  PATH* path, /**< Voronoi paths data struture */
2218  int* newedge, /**< the new edge */
2219  int crucnode, /**< the current crucial node */
2220  UF* uf /**< union find data structure */
2221  )
2222 {
2223  int k;
2224  int m;
2225  int i;
2226  int e;
2227  int* heap;
2228  int* state;
2229  int node1;
2230  int node2;
2231 
2232  *newedge = UNKNOWN;
2233  e = UNKNOWN;
2234  assert(g != NULL);
2235  assert(g->path_heap != NULL);
2236  assert(g->path_state != NULL);
2237  assert(path != NULL);
2238  assert(cost != NULL);
2239 
2240  if( g->knots == 0 )
2241  return;
2242 
2243  heap = g->path_heap;
2244  state = g->path_state;
2245 
2246  if( g->knots > 1 )
2247  {
2248  /* until the heap is empty */
2249  while( *count > 0 )
2250  {
2251  /* get the next (i.e. a nearest) vertex from the heap */
2252  k = nearest(heap, state, count, path);
2253 
2254  /* mark vertex k as scanned */
2255  state[k] = CONNECT;
2256  assert(g->mark[k]);
2257  assert(g->mark[vbase[k]]);
2258  /* iterate over all outgoing edges of vertex k */
2259  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2260  {
2261  m = g->head[i];
2262 
2263  /* check whether the path (to m) including k is shorter than the so far best known */
2264  if( (state[m]) && (SCIPisGT(scip, path[m].dist, path[k].dist + cost[i])) )/*&& g->mark[m] ) */
2265  {
2266  assert(g->mark[m]);
2267  correct(scip, heap, state, count, path, m, k, i, cost[i], FSP_MODE);
2268  vbase[m] = vbase[k];
2269  }
2270 
2271  /* check whether there is a better new boundary edge adjacent to vertex k */
2272  else
2273  {
2274  node1 = SCIPunionfindFind(uf, vbase[m]);
2275  node2 = SCIPunionfindFind(uf, vbase[k]);
2276  if( state[m] == CONNECT && ((node1 == crucnode) != (node2 == crucnode)) && g->mark[m] && g->mark[vbase[m]]
2277  && g->mark[node1] && g->mark[node2] && ((e == UNKNOWN) || SCIPisGT(scip,
2278  (path[g->tail[e]].dist + cost[e] + path[g->head[e]].dist), path[k].dist + cost[i] + path[m].dist)) )
2279  e = i;
2280  }
2281  }
2282  }
2283  }
2284  *newedge = e;
2285 }
2286 
2287 
2288 /** repair the voronoi diagram for a given set nodes */
2289 void voronoi_repair_mult(
2290  SCIP* scip, /**< SCIP data structure */
2291  const GRAPH* g, /**< graph data structure */
2292  SCIP_Real* cost, /**< edge costs */
2293  int* count, /**< pointer to number of heap elements */
2294  int* vbase, /**< array containing Voronoi base of each node */
2295  int* boundedges, /**< boundary edges */
2296  int* nboundedges, /**< number of boundary edges */
2297  char* nodesmark, /**< array to mark temporarily discarded nodes */
2298  UF* uf, /**< union find data structure */
2299  PATH* path /**< Voronoi paths data struture */
2300  )
2301 {
2302  int k;
2303  int m;
2304  int i;
2305  int* heap;
2306  int* state;
2307  int node1;
2308  int node2;
2309 
2310  assert(g != NULL);
2311  assert(g->path_heap != NULL);
2312  assert(g->path_state != NULL);
2313  assert(path != NULL);
2314  assert(cost != NULL);
2315 
2316  if( g->knots == 0 )
2317  return;
2318 
2319  heap = g->path_heap;
2320  state = g->path_state;
2321 
2322  if( g->knots > 1 )
2323  {
2324  /* until the heap is empty */
2325  while( (*count) > 0 )
2326  {
2327  /* get the next (i.e. a nearest) vertex from the heap */
2328  k = nearest(heap, state, count, path);
2329 
2330  /* mark vertex k as scanned */
2331  state[k] = CONNECT;
2332 
2333  /* iterate all outgoing edges of vertex k */
2334  for( i = g->outbeg[k]; i != EAT_LAST; i = g->oeat[i] )
2335  {
2336 
2337  m = g->head[i];
2338  /* printf("scrut edge %d->%d vbases: %d %d \n ", k,m, vbase[k], vbase[m]);
2339  printf(" uf : %d %d \n ", SCIPunionfindFind(uf, vbase[k]), SCIPunionfindFind(uf, vbase[m]) ); */
2340  /* check whether the path (to m) including k is shorter than the so far best known */
2341  if( state[m] && (SCIPisGT(scip,path[m].dist, path[k].dist + cost[i])) && g->mark )
2342  {
2343  correct(scip, heap, state, count, path, m, k, i, cost[i], FSP_MODE);
2344  vbase[m] = vbase[k];
2345  }
2346  /* check whether there is a new boundary edge adjacent to vertex k */
2347  else if( (state[m] == CONNECT) && ((node1 = SCIPunionfindFind(uf, vbase[m])) != (node2 = SCIPunionfindFind(uf, vbase[k])))
2348  && g->mark[node1] && g->mark[node2] && (nodesmark[node1] || nodesmark[node2]) )
2349  {
2350  boundedges[(*nboundedges)++] = i;
2351  /*printf("adding new boundaryedge [%d] %d_%d \n", (*nboundedges) - 1, k,m);*/
2352  }
2353  }
2354  }
2355  }
2356 }
2357 
2358 /* ARGSUSED */
2359 void graph_path_length(
2360  const GRAPH* g, /**< graph data structure */
2361  const PATH* p /**< path data structure */
2362  )
2363 {
2364 #ifndef NDEBUG
2365  int len = 0;
2366  SCIP_Real dist = 0.0;
2367  int e;
2368  int i;
2369 
2370  for(i = 0; i < g->knots; i++)
2371  {
2372  if ((e = p[i].edge) < 0)
2373  continue;
2374 
2375  len++;
2376  dist += g->cost[e];
2377  }
2378  (void)printf("Graph path length: %d Distance=%g\n", len, dist);
2379 
2380 #endif /* NDEBUG */
2381 }
Is_term
#define Is_term(a)
Definition: grph.h:158
voronoi_repair
void voronoi_repair(SCIP *scip, const GRAPH *g, SCIP_Real *cost, int *count, int *vbase, PATH *path, int *newedge, int crucnode, UF *uf)
Definition: grphpath.c:2211
LT
#define LT(a, b)
Definition: portab.h:64
voronoiSteinerTreeExt
void voronoiSteinerTreeExt(SCIP *scip, const GRAPH *g, SCIP_Real *costrev, int *vbase, char *stvertex, PATH *path)
Definition: grphpath.c:1121
getnext4terms
void getnext4terms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1562
GRAPH
Definition: grph.h:55
TRUE
#define TRUE
Definition: portab.h:34
shortest_path::edge
signed int edge
Definition: grph.h:140
MST_MODE
#define MST_MODE
Definition: grph.h:152
GRAPH::path_heap
int * path_heap
Definition: grph.h:114
GRAPH::terms
int terms
Definition: grph.h:63
getnext4tterms
void getnext4tterms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *boundedges, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1603
SCIPunionfindFind
int SCIPunionfindFind(UF *uf, int element)
Definition: misc_stp.c:620
GT
#define GT(a, b)
Definition: portab.h:66
STP_HOP_CONS
#define STP_HOP_CONS
Definition: grph.h:48
UnionFind_Structure
Definition: misc_stp.h:42
EAT_LAST
#define EAT_LAST
Definition: grph.h:31
Edge_anti
#define Edge_anti(a)
Definition: grph.h:161
GRAPH::path_state
int * path_state
Definition: grph.h:115
shortest_path
Definition: grph.h:137
sdpaths
void 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:542
get4next
void get4next(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1417
nearestX
static int nearestX(int *heap, int *state, int *count, const SCIP_Real *pathdist)
Definition: grphpath.c:92
graph_path_length
void graph_path_length(const GRAPH *g, const PATH *p)
Definition: grphpath.c:2359
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:144
GRAPH::oeat
int * oeat
Definition: grph.h:100
voronoi_extend2
SCIP_RETCODE voronoi_extend2(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, SCIP_Real **distarr, int **basearr, int **edgearr, char *termsmark, int *reachednodes, int *nreachednodes, int *nodenterms, int nneighbterms, int base, int countex)
Definition: grphpath.c:954
voronoi_radius
SCIP_RETCODE voronoi_radius(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:1912
GRAPH::head
int * head
Definition: grph.h:94
Vnoi_List_Node
Definition: misc_stp.h:50
FALSE
#define FALSE
Definition: portab.h:37
GRAPH::mark
int * mark
Definition: grph.h:71
Vnoi_List_Node::edge
signed int edge
Definition: misc_stp.h:53
CONNECT
#define CONNECT
Definition: grph.h:147
GRAPH::stp_type
int stp_type
Definition: grph.h:123
Vnoi_List_Node::next
struct Vnoi_List_Node * next
Definition: misc_stp.h:55
GRAPH::outbeg
int * outbeg
Definition: grph.h:77
Is_pterm
#define Is_pterm(a)
Definition: grph.h:159
GRAPH::prize
SCIP_Real * prize
Definition: grph.h:92
graph_path_exit
void graph_path_exit(SCIP *scip, GRAPH *g)
Definition: grphpath.c:429
GRAPH::tail
int * tail
Definition: grph.h:93
voronoi_repair_mult
void voronoi_repair_mult(SCIP *scip, const GRAPH *g, SCIP_Real *cost, int *count, int *vbase, int *boundedges, int *nboundedges, char *nodesmark, UF *uf, PATH *path)
Definition: grphpath.c:2289
FSP_MODE
#define FSP_MODE
Definition: grph.h:153
GRAPH::term
int * term
Definition: grph.h:69
GRAPH::knots
int knots
Definition: grph.h:61
Vnoi_List_Node::base
signed int base
Definition: misc_stp.h:54
voronoi
void voronoi(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, char *base, int *vbase, PATH *path)
Definition: grphpath.c:1034
resetX
static void resetX(SCIP *scip, SCIP_Real *pathdist, int *heap, int *state, int *count, int node)
Definition: grphpath.c:272
STP_MAX_NODE_WEIGHT
#define STP_MAX_NODE_WEIGHT
Definition: grph.h:47
voronoi_dist
SCIP_RETCODE voronoi_dist(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:1755
FARAWAY
#define FARAWAY
Definition: grph.h:149
voronoi_terms
void voronoi_terms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1679
BSP_MODE
#define BSP_MODE
Definition: grph.h:154
UNKNOWN
#define UNKNOWN
Definition: grph.h:148
get2next
void get2next(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1214
STP_PRIZE_COLLECTING
#define STP_PRIZE_COLLECTING
Definition: grph.h:40
get3next
void get3next(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path, int *vbase, int *heap, int *state)
Definition: grphpath.c:1310
grph.h
includes various files containing graph methods used for Steiner problems
graph_path_init
SCIP_RETCODE graph_path_init(SCIP *scip, GRAPH *g)
Definition: grphpath.c:407
portab.h
Portable defintions.
GRAPH::grad
int * grad
Definition: grph.h:74
graph_path_exec
void graph_path_exec(SCIP *scip, const GRAPH *p, int mode, int start, SCIP_Real *cost, PATH *path)
Definition: grphpath.c:453
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:196
GRAPH::cost
SCIP_Real * cost
Definition: grph.h:91
voronoi_extend
SCIP_RETCODE voronoi_extend(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, VLIST **adjterms, char *termsmark, int *reachednodes, int *nreachednodes, int *nodenterms, int nneighbterms, int base, int countex)
Definition: grphpath.c:862
getnext3terms
void getnext3terms(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, PATH *path3, int *vbase, int *heap, int *state)
Definition: grphpath.c:1526
STP_ROOTED_PRIZE_COLLECTING
#define STP_ROOTED_PRIZE_COLLECTING
Definition: grph.h:41
heap_add
void heap_add(int *heap, int *state, int *count, int node, PATH *path)
Definition: grphpath.c:241
GRAPH::source
int * source
Definition: grph.h:67
GRAPH::edges
int edges
Definition: grph.h:89
graph_knot_add
void graph_knot_add(GRAPH *, int)
Definition: grphbase.c:1362
Vnoi_List_Node::dist
double dist
Definition: misc_stp.h:52
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:306
graph_path_execX
void graph_path_execX(SCIP *scip, const GRAPH *g, int start, SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge)
Definition: grphpath.c:639
graph_path_st
void graph_path_st(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *pathdist, int *pathedge, int start, unsigned int *randseed, char *connected)
Definition: grphpath.c:707
nearest
static int nearest(int *heap, int *state, int *count, const PATH *path)
Definition: grphpath.c:43
shortest_path::dist
double dist
Definition: grph.h:139
graph_edge_add
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
voronoi_slrepair
void voronoi_slrepair(SCIP *scip, const GRAPH *g, SCIP_Real *cost, PATH *path, int *vbase, int *heap, int *state, int start, int adjnode)
Definition: grphpath.c:2136