Scippy

SCIP

Solving Constraint Integer Programs

reduce_alt.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-2019 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 scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file reduce_alt.c
17  * @brief Altenative based reduction tests for Steiner problems
18  * @author Thorsten Koch
19  * @author Stephen Maher
20  * @author Daniel Rehfeldt
21  *
22  * This file implements alternative-based reduction techniques for several Steiner problems.
23  * All tests can be found in "Combining NP-Hard Reduction Techniques and Strong Heuristics in an Exact Algorithm for the
24  * Maximum-Weight Connected Subgraph Problem" by Daniel Rehfeldt and Thorsten Koch,
25  * or in "Reduction Techniques for the Prize-Collecting Steiner Tree Problem and the Maximum-Weight Connected Subgraph Problem"
26  * by Daniel Rehfeldt et al.
27  *
28  *
29  * A list of all interface methods can be found in grph.h.
30  *
31  */
32 
33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include "grph.h"
39 #include "misc_stp.h"
40 #include "probdata_stp.h"
41 #include "scip/scip.h"
42 
43 #define VERTEX_CONNECT 0
44 #define VERTEX_TEMPNEIGHBOR 1
45 #define VERTEX_NEIGHBOR 2
46 #define VERTEX_OTHER 3
47 #define STP_RED_CNSNN 25
48 #define STP_RED_ANSMAXCANDS 500
49 #define STP_RED_ANSEXMAXCANDS 50
50 #define STP_RED_ANSMAXNEIGHBORS 25
51 
52 
53 /* can edge be deleted in SD test in case of equality? If so, 'forbidden' array is adapted */
54 static
55 SCIP_Bool sddeltable(
56  SCIP* scip,
57  GRAPH* g,
58  PATH* path,
59  int* vbase,
60  int* forbidden,
61  int tpos,
62  int hpos,
63  int tail,
64  int head,
65  int edge,
66  int nnodes
67  )
68 {
69  int e;
70 
71  assert(g != NULL);
72  assert(path != NULL);
73  assert(scip != NULL);
74  assert(vbase != NULL);
75  assert(forbidden != NULL);
76 
77  assert(tpos >= 0);
78  assert(hpos >= 0);
79  assert(tail >= 0);
80  assert(head >= 0);
81  assert(edge >= 0);
82  assert(nnodes >= 0);
83 
84  /* edge between non-terminals */
85  if( !Is_term(g->term[tail]) && !Is_term(g->term[head]) )
86  return TRUE;
87 
88  /* has edge been marked as forbidden? */
89  if( forbidden[edge] )
90  return FALSE;
91 
92  /* edge between a terminal and a non terminal */
93  if( !Is_term(g->term[tail]) || !Is_term(g->term[head]) )
94  {
95 
96  /* check whether edge is used in shortest path */
97 
98  if( !Is_term(g->term[tail]) && Is_term(g->term[head]) )
99  {
100  e = path[tail + tpos * nnodes].edge;
101 
102  if( g->ieat[e] == EAT_FREE )
103  return TRUE;
104 
105  assert(g->head[e] == tail);
106 
107  if( g->tail[e] == head )
108  return FALSE;
109  }
110  else if( Is_term(g->term[tail]) && !Is_term(g->term[head]) )
111  {
112  e = path[head + hpos * nnodes].edge;
113 
114  if( g->ieat[e] == EAT_FREE )
115  return TRUE;
116 
117  assert(g->head[e] == head);
118 
119  if( g->tail[e] == tail )
120  return FALSE;
121  }
122  }
123 
124  /* update forbidden edges todo check bottleneck distance between terminals, and don't delete if distance is higher than ecost */
125 
126  if( Is_term(g->term[head]) )
127  {
128  SCIP_Real ecost = g->cost[edge];
129  for( e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
130  {
131  assert(e >= 0);
132 
133  if( SCIPisEQ(scip, g->cost[e], ecost) )
134  {
135  if( !(forbidden[e]) && Is_term(g->term[g->head[e]]) )
136  {
137  forbidden[e] = TRUE;
138  forbidden[flipedge(e)] = TRUE;
139  }
140  }
141  }
142  }
143 
144  if( Is_term(g->term[tail]) )
145  {
146  SCIP_Real ecost = g->cost[edge];
147  for( e = g->outbeg[head]; e != EAT_LAST; e = g->oeat[e] )
148  {
149  assert(e >= 0);
150 
151  if( SCIPisEQ(scip, g->cost[e], ecost) )
152  {
153  if( !(forbidden[e]) && Is_term(g->term[g->head[e]]) )
154  {
155  forbidden[e] = TRUE;
156  forbidden[flipedge(e)] = TRUE;
157  }
158  }
159  }
160  }
161 
162  return TRUE;
163 }
164 
165 
166 static
167 int getcloseterms(
168  SCIP* scip,
169  const PATH* vnoi,
170  SCIP_Real* termdist,
171  SCIP_Real ecost,
172  const int* vbase,
173  int* neighbterms,
174  int i,
175  int nnodes
176  )
177 {
178  int k;
179  int pos = i;
180  int nnterms = 0;
181 
182  assert(scip != NULL);
183  assert(vnoi != NULL);
184  assert(vbase != NULL);
185  assert(termdist != NULL);
186  assert(neighbterms != NULL);
187 
188  for( k = 0; k < 4; k++ )
189  {
190  if( SCIPisLT(scip, vnoi[pos].dist, ecost) )
191  {
192  neighbterms[nnterms] = vbase[pos];
193  termdist[nnterms++] = vnoi[pos].dist;
194  }
195  else
196  {
197  break;
198  }
199  pos += nnodes;
200  }
201 
202  return nnterms;
203 }
204 
205 static
206 int getcloseterms2term(
207  SCIP* scip,
208  const GRAPH* g,
209  const GRAPH* netgraph,
210  SCIP_Real* termdist,
211  SCIP_Real ecost,
212  int* neighbterms,
213  int* nodesid,
214  int* nodesorg,
215  int i
216 )
217 {
218  int nnterms = 0;
219 
220  for( int k = 0; k < 4; k++ )
221  neighbterms[k] = UNKNOWN;
222 
223  /* get three nearest terms */
224  for( int ne = netgraph->outbeg[nodesid[i]]; ne != EAT_LAST; ne = netgraph->oeat[ne] )
225  {
226  if( SCIPisLT(scip, netgraph->cost[ne], ecost) )
227  {
228  const SCIP_Real necost = netgraph->cost[ne];
229  int j = nodesorg[netgraph->head[ne]];
230 
231  assert(Is_term(g->term[j]));
232  assert(j != i);
233 
234  if( nnterms < 4 )
235  nnterms++;
236  for( int k = 0; k < 4; k++ )
237  {
238  if( neighbterms[k] == UNKNOWN || SCIPisGT(scip, termdist[k], necost) )
239  {
240  for( int l = 3; l > k; l-- )
241  {
242  neighbterms[l] = neighbterms[l - 1];
243  termdist[l] = termdist[l - 1];
244  }
245  neighbterms[k] = j;
246  termdist[k] = necost;
247  break;
248  }
249  }
250  }
251  }
252 
253  return nnterms;
254 }
255 
256 static
257 int getlecloseterms(
258  SCIP* scip,
259  PATH* vnoi,
260  SCIP_Real* termdist,
261  SCIP_Real ecost,
262  int* vbase,
263  int* neighbterms,
264  int i,
265  int nnodes
266  )
267 {
268  int k;
269  int pos = i;
270  int nnterms = 0;
271 
272  assert(scip != NULL);
273  assert(vnoi != NULL);
274  assert(vbase != NULL);
275  assert(termdist != NULL);
276  assert(neighbterms != NULL);
277 
278  for( k = 0; k < 4; k++ )
279  {
280  if( SCIPisLE(scip, vnoi[pos].dist, ecost) )
281  {
282  neighbterms[nnterms] = vbase[pos];
283  termdist[nnterms++] = vnoi[pos].dist;
284  }
285  else
286  {
287  break;
288  }
289  pos += nnodes;
290  }
291 
292  return nnterms;
293 }
294 
295 #if 0
296 static int issmaller(
297  SCIP* scip,
298  const double* pathdist,
299  const double* pathtran,
300  int a,
301  int b)
302 {
303  return (SCIPisLT(scip, pathdist[a], pathdist[b]) || (!SCIPisGT(scip, pathdist[a], pathdist[b]) && SCIPisLT(scip, pathtran[a], pathtran[b])));
304 }
305 static int islarger(
306  SCIP* scip,
307  const double* pathdist,
308  const double* pathtran,
309  int a,
310  int b)
311 {
312  return (SCIPisGT(scip, pathdist[a], pathdist[b]) || (!SCIPisLT(scip, pathdist[a], pathdist[b]) && SCIPisGT(scip, pathtran[a], pathtran[b])));
313 }
314 /*---------------------------------------------------------------------------*/
315 /*--- Name : get NEAREST knot ---*/
316 /*--- Function : Holt das oberste Element vom Heap (den Knoten mit der ---*/
317 /*--- geringsten Entfernung zu den bereits Verbundenen) ---*/
318 /*--- Parameter: Derzeitige Entfernungen und benutzte Kanten ---*/
319 /*--- Returns : Nummer des bewussten Knotens ---*/
320 /*---------------------------------------------------------------------------*/
321 inline static int nearest(
322  SCIP* scip,
323  int* heap,
324  int* state,
325  int* count, /* pointer to store number of elements of heap */
326  const double* pathdist,
327  const double* pathtran)
328 {
329  int k;
330  int t;
331  int c;
332  int j;
333 
334  /* Heap shift down
335  * (Oberstes Element runter und korrigieren)
336  */
337  k = heap[1];
338  j = 1;
339  c = 2;
340  heap[1] = heap[(*count)--];
341  state[heap[1]] = 1;
342 
343  if ((*count) > 2)
344  if (islarger(scip, pathdist, pathtran, heap[2], heap[3]))
345  c++;
346 
347  while((c <= (*count)) && islarger(scip, pathdist, pathtran, heap[j], heap[c]))
348  {
349  t = heap[c];
350  heap[c] = heap[j];
351  heap[j] = t;
352  state[heap[j]] = j;
353  state[heap[c]] = c;
354  j = c;
355  c += c;
356 
357  if ((c + 1) <= (*count))
358  if (issmaller(scip, pathdist, pathtran, heap[c + 1], heap[c]))
359  c++;
360  }
361  return(k);
362 }
363 
364 
365 /** insert respectively change element in heap */
366 inline static void correct(
367  SCIP* scip,
368  int* heap,
369  int* state,
370  int* count, /* pointer to store number of elements of heap */
371  double* pathdist,
372  double* pathtran,
373  int l)
374 {
375  int t;
376  int c;
377  int j;
378 
379  /* Ist der Knoten noch ganz frisch ?
380  */
381  if (state[l] == UNKNOWN)
382  {
383  heap[++(*count)] = l;
384  state[l] = (*count);
385  }
386 
387  /* Heap shift up
388  */
389  j = state[l];
390  c = j / 2;
391 
392  while((j > 1) && (islarger(scip, pathdist, pathtran, heap[c], heap[j])))
393  {
394  t = heap[c];
395  heap[c] = heap[j];
396  heap[j] = t;
397  state[heap[j]] = j;
398  state[heap[c]] = c;
399  j = c;
400  c = j / 2;
401  }
402 
403 }
404 
405 /** Dijkstra's algorithm for shortest path or minimum spanning tree */
406 static
407 void compute_sd(
408  SCIP* scip,
409  const GRAPH* p,
410  int start,
411  const double* cost,
412  const double* randarr,
413  int* heap,
414  int* state,
415  int* count, /* pointer to store number of elements of heap */
416  double* pathdist,
417  double* pathtran,
418  double* pathrand
419  )
420 {
421  int k;
422  int m;
423  int i;
424  int done = 0;
425  double tran;
426  double dist;
427  double temprand;
428 
429  assert(scip != NULL);
430  assert(p != NULL);
431  assert(start >= 0);
432  assert(start < p->knots);
433  assert(heap != NULL);
434  assert(state != NULL);
435  assert(pathdist != NULL);
436  assert(pathtran != NULL);
437  assert(cost != NULL);
438  assert(count != NULL);
439  assert(*count >= 0);
440 
441  /* Kein Baum ohne Knoten
442  */
443  if (p->knots == 0)
444  return;
445 
446  (*count) = 0;
447 
448  /* Erstmal alles auf null, unbekannt und weit weg
449  */
450  for(i = 0; i < p->knots; i++)
451  {
452  state[i] = UNKNOWN;
453  pathdist[i] = FARAWAY;
454  pathtran[i] = FARAWAY;
455  }
456  /* Startknoten in den Heap
457  */
458  k = start;
459  pathdist[k] = 0.0;
460  pathtran[k] = 0.0;
461  pathrand[k] = 0.0;
462 
463  /* Wenn nur ein Knoten drin ist funktioniert der Heap nicht sehr gut,
464  * weil dann fuer genau 0 Elemente Platz ist.
465  */
466  if (p->knots > 1)
467  {
468  (*count) = 1;
469  heap[(*count)] = k;
470  state[k] = (*count);
471 
472  /* Wenn nichts mehr auf dem Heap ist, sind wir fertig
473  * und jetzt erstmal Hula Loop
474  */
475  while((*count) > 0)
476  {
477  /* Na, wer ist der Naechste ?
478  */
479  k = nearest(scip, heap, state, count, pathdist, pathtran);
480 
481  /* Wieder einen erledigt
482  */
483  state[k] = CONNECT;
484 
485  if (p->mark[k] == 2)
486  if (++done >= p->grad[start])
487  break;
488 
489  /* Verbunden Knoten berichtigen ...
490  *
491  * Wenn ein Knoten noch nicht erledigt ist
492  * werden wir dann auf diesem Wege besser ?
493  */
494  for(i = p->outbeg[k]; i != EAT_LAST; i = p->oeat[i])
495  {
496  m = p->head[i];
497 
498  /* 1. Ist der Knoten noch nicht festgelegt ?
499  * Ist der wohlmoeglich tabu ?
500  */
501  if ((state[m]) && (p->mark[m]))
502  {
503  /* 2. Ist es ueberhaupt eine Verbesserung diesen Weg zu nehmen ?
504  * Wenn ja, dann muss die Entferung kuerzer sein.
505  */
506  /* The special distance is the length of the longest path between two terminals (elementary path)
507  * contained in the path between knots i and j.
508  * - tran measures the distance between two terminals.
509  * - dist stores the current longest elementary path.
510  */
511  if( Is_term(p->term[m]) )
512  {
513  tran = 0.0;
514  temprand = 0.0;
515  }
516  else
517  {
518  tran = pathtran[k] + cost[i];
519  temprand = pathrand[k] + randarr[i];
520  }
521 
522  if( SCIPisGE(scip, pathdist[k], pathtran[k] + cost[i]) )
523  dist = pathdist[k];
524  else
525  dist = pathtran[k] + cost[i];
526 
527  if( SCIPisLT(scip, dist, pathdist[m])
528  || (SCIPisEQ(scip, dist, pathdist[m]) && SCIPisLT(scip, tran, pathtran[m])) )
529  {
530  pathdist[m] = dist;
531  pathtran[m] = tran;
532  pathrand[m] = temprand;
533 
534  correct(scip, heap, state, count, pathdist, pathtran, m);
535  }
536  }
537  }
538  }
539  }
540 }
541 #endif
542 
543 
544 
545 #if 0
546 /* C. W. Duin and A. Volganant
547  *
548  * "Reduction Tests for the Steiner Problem in Graphs"
549  *
550  * Networks, Volume 19 (1989), 549-567
551  *
552  * Nearest Special Vertex 3 Test
553  *
554  * Taken from:
555  *
556  * Maculan et. al.
557  *
558  * "An approach for the Steiner problem in directed graphs"
559  *
560  * Annals of Operations Research, Volume 33 (1991), 471-480
561  *
562  * Nearest Vertex Test (for optimal arcs)
563  *
564  * and
565  *
566  * T. Polzin
567  *
568  * "Algorithms for the Steiner problem in networks"
569  *
570  * Section 3.3.3 pp. 54-55
571  *
572  * This is only called for the directed Steiner tree problem
573  */
574 SCIP_RETCODE reduce_nv_optimal(
575  SCIP* scip,
576  GRAPH* g,
577  double* fixed,
578  int* elimins,
579  int runnum)
580 {
581  PATH** path;
582  PATH* pathfromterm;
583  PATH* pathfromsource;
584  PATH* pathhops;
585  double* distance;
586  double* radius;
587  double* hopscost;
588  int* vregion;
589  int* heap;
590  int* state;
591  int* pred;
592  int* minArc1;
593  int* minArc2;
594  int* terms;
595  int termcount;
596  int i;
597  int j;
598  int e;
599  double min1;
600  double min2;
601  int minhops;
602  int shortarc;
603  int shortarctail;
604  STP_Bool antiedgeexists;
605  int knotoffset;
606 
607  SCIPdebugMessage("NSV-Reduction: ");
608  fflush(stdout);
609  /*
610  graph_show(g);
611  */
612  path = malloc((size_t)g->knots * sizeof(PATH*));
613 
614  assert(path != NULL);
615 
616  pathfromterm = malloc((size_t)g->knots * sizeof(PATH));
617  pathfromsource = malloc((size_t)g->knots * sizeof(PATH));
618 
619  assert(pathfromterm != NULL);
620  assert(pathfromsource != NULL);
621 
622  pathhops = malloc((size_t)g->knots * sizeof(PATH));
623 
624  assert(pathhops != NULL);
625 
626  distance = malloc((size_t)g->knots * sizeof(double));
627  radius = malloc((size_t)g->knots * sizeof(double));
628 
629  assert(distance != NULL);
630  assert(radius != NULL);
631 
632  hopscost = malloc((size_t)g->edges * sizeof(double));
633 
634  assert(hopscost != NULL);
635 
636  vregion = malloc((size_t)g->knots * sizeof(int));
637 
638  assert(vregion != NULL);
639 
640  heap = malloc((size_t)g->knots * sizeof(int));
641  state = malloc((size_t)g->knots * sizeof(int));
642 
643  assert(heap != NULL);
644  assert(state != NULL);
645 
646  *elimins = 0;
647  pred = malloc((size_t)g->edges * sizeof(int));
648  minArc1 = malloc((size_t)g->knots * sizeof(int));
649  minArc2 = malloc((size_t)g->knots * sizeof(int));
650  terms = malloc((size_t)g->terms * sizeof(int));
651 
652  termcount = 0;
653  for(i = 0; i < g->knots; i++)
654  {
655  if( Is_term(g->term[i]) )
656  {
657  terms[termcount] = i;
658  termcount++;
659  }
660  g->mark[i] = (g->grad[i] > 0);
661  minArc1[i] = -1;
662  minArc2[i] = -1;
663  path[i] = NULL;
664 
665  for(e = g->inpbeg[i]; e != EAT_LAST; e = g->ieat[e])
666  {
667  if( LT(g->cost[e], FARAWAY) )
668  hopscost[e] = 1;
669  else
670  hopscost[e] = FARAWAY;
671 
672  if( LT(g->cost[Edge_anti(e)], FARAWAY) )
673  hopscost[Edge_anti(e)] = 1;
674  else
675  hopscost[Edge_anti(e)] = FARAWAY;
676  }
677  }
678 
679  assert(g->source >= 0);
680 
681  /* computing the voronoi regions inward to a node */
682  voronoi_term(g, g->cost, distance, radius, pathfromterm, vregion, heap, state, pred, 1);
683 
684  /* computing the shortest paths from the source node */
685  graph_path_exec(scip, g, FSP_MODE, g->source, g->cost, pathfromsource);
686 
687  /* computing the shortest hops paths from the source node */
688  graph_path_exec(scip, g, FSP_MODE, g->source, hopscost, pathhops);
689 
690  /* this is the offset used to minimise the number of knots to examine in large graphs. */
691  srand(runnum*100);
692  knotoffset = rand() % KNOTFREQ;
693 
694  for(i = 0; i < g->knots; i++)
695  {
696  /* For the prize collecting variants all edges from the "dummy" root node must be retained. */
697  if ((g->stp_type == STP_PCSPG || g->stp_type == STP_MWCSP) && i == g->source )
698  continue;
699 
700  if( g->stp_type == STP_SAP && i != g->source )
701  continue;
702 
703 
704  if( g->knots > KNOTLIMIT && i % KNOTFREQ != knotoffset )
705  continue;
706 
707  if (Is_term(g->term[i]) && g->grad[i] >= 3)
708  {
709 
710  min1 = FARAWAY;
711  min2 = FARAWAY;
712  minhops = g->hoplimit;
713  shortarctail = -1;
714  shortarc = -1;
715  antiedgeexists = FALSE;
716  for(e = g->inpbeg[i]; e != EAT_LAST; e = g->ieat[e])
717  {
718  if( g->stp_type == STP_SAP && i == g->source )
719  {
720  if( LT(g->cost[e], FARAWAY) )
721  {
722  g->cost[e] = FARAWAY;
723  (*elimins)++;
724  }
725 
726  continue;
727  }
728 
729  if (g->cost[e] < min1)
730  {
731  shortarc = e;
732  shortarctail = g->tail[e];
733 
734  min2 = min1;
735  min1 = g->cost[e];
736  }
737 
738  if( LT(pathfromsource[g->tail[e]].hops, minhops) )
739  minhops = pathfromsource[g->tail[e]].hops;
740 
741  if( LT(g->cost[Edge_anti(e)], FARAWAY) )
742  antiedgeexists = TRUE;
743 
744  }
745 
746  if( g->stp_type == STP_SAP )
747  continue;
748 
749  if (LT(min1, FARAWAY) && LE(pathfromsource[shortarctail].dist + min1, min2))
750  {
751  assert(shortarc >= 0);
752  assert(shortarctail >= 0);
753 
754  if ((g->stp_type == STP_PCSPG || g->stp_type == STP_MWCSP
755  || g->stp_type == STP_SAP) && shortarctail == g->source )
756  continue;
757 
758  if( g->stp_type == STP_DHCSTP && GT(pathfromsource[shortarctail].hops, pathhops[i].dist - 1) )
759  continue;
760 
761  if( antiedgeexists == TRUE )
762  {
763  if( LT(min2, FARAWAY) )
764  {
765  for(e = g->inpbeg[i]; e != EAT_LAST; e = j)
766  {
767  j = g->ieat[e];
768 
769  if( e != shortarc )
770  {
771  if( LT(g->cost[Edge_anti(e)], FARAWAY) )
772  g->cost[e] = FARAWAY;
773  else
774  graph_edge_del(scip, g, e, TRUE);
775  }
776  }
777  (*elimins)++;
778  }
779  }
780  else
781  {
782  *fixed += min1;
783  SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->fixedges), g->ancestors[shortarc], NULL) ); /* I think that this should be
784  shortarc instead of shortarctail */
785  SCIP_CALL( graph_knot_contract(scip, g, i, shortarctail) );
786 
787  if( g->stp_type == STP_DHCSTP )
788  {
789  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
790  g->cost[e] = FARAWAY;
791  }
792 
793  (*elimins)++;
794  }
795 
796  /* computing the shortest paths from the source node */
797  graph_path_exec(scip, g, FSP_MODE, g->source, g->cost, pathfromsource);
798  }
799  }
800  /* The knot is not a terminal so we can perform the short link test */
801 #if 0
802  else
803  {
804  for(e = g->inpbeg[i]; e != EAT_LAST; e = g->ieat[e])
805  {
806  j = g->tail[e];
807  if( vregion[i] != vregion[j] )
808  {
809  if( minArc1[vregion[i]] < 0 )
810  minArc1[vregion[i]] = e;
811  else if( g->cost[e] < g->cost[minArc1[vregion[i]]] )
812  {
813  minArc2[vregion[i]] = minArc1[vregion[i]];
814  minArc1[vregion[i]] = e;
815  }
816  }
817  }
818  }
819 #endif
820  }
821 
822 #if 0
823  for( k = 0; k < termcount; k++ )
824  {
825  assert(terms[k] >= 0 && terms[k] < g->knots);
826 
827  if( minArc1[terms[k]] >= 0 && minArc2[terms[k]] >= 0 && pathfromsource[g->tail[minArc1[terms[k]]]].dist
828  + g->cost[minArc1[terms[k]]] + pathfromterm[g->head[minArc1[terms[k]]]].dist < g->cost[minArc2[terms[k]]] )
829  {
830  e = minArc1[terms[k]];
831  i = g->head[e];
832  j = g->tail[e];
833 
834  if ((g->stp_type == STP_PCSPG || g->stp_type == STP_MWCSP) && (i == g->source || j == g->source) )
835  continue;
836 
837 
838  if( Is_term(g->term[i]) )
839  {
840  SCIPintListNodeAppendCopy(&(g->fixedges), g->ancestors[e], NULL);
841  *fixed += g->cost[e];
842  }
843  graph_knot_contract(g, j, i);
844 
845 
846  elimins++;
847  }
848  }
849 #endif
850 
851  free(terms);
852  free(minArc2);
853  free(minArc1);
854  free(pred);
855  free(state);
856  free(heap);
857  free(hopscost);
858  free(radius);
859  free(distance);
860  free(vregion);
861  free(pathhops);
862  free(pathfromsource);
863  free(pathfromterm);
864  free(path);
865 
866  assert(graph_valid(g));
867  SCIPdebugMessage(" %d Knots deleted\n", *elimins);
868 
869  return SCIP_OKAY;
870 }
871 
872 * C. W. Duin and A. Volganant
873  *
874  * "An Edge Elimination Test for the Steiner Problem in Graphs"
875  *
876  * Operations Research Letters 8, (1989), 79-83
877  *
878  * Special Distance Test
879  */
880 int reduce_sduction(
881  SCIP* scip,
882  GRAPH* g,
883  double* sddist,
884  double* sdtrans,
885  double* sdrand,
886  double* cost,
887  double* random,
888  int* heap,
889  int* state,
890  int* knotexamined,
891  int runnum
892  )
893 {
894  SCIP_Real redstarttime;
895  SCIP_Real timelimit;
896  SCIP_Real stalltime;
897  int count = 0;
898  int i;
899  int e;
900  int j;
901  int elimins = 0;
902  int knotoffset = 0;
903 
904  SCIPdebugMessage("SD-Reduktion: ");
905  fflush(stdout);
906 
907  /*
908  heap = malloc((size_t)g->knots * sizeof(int));
909  state = malloc((size_t)g->knots * sizeof(int));
910  */
911  assert(heap != NULL);
912  assert(state != NULL);
913  /*
914  sd = malloc((size_t)g->knots * sizeof(SDPTH));
915  */
916  assert(sddist != NULL);
917  assert(sdtrans != NULL);
918  /*
919  cost = malloc((size_t)g->edges * sizeof(double));
920  */
921  assert(cost != NULL);
922 
923  assert(knotexamined != NULL);
924 
925  redstarttime = SCIPgetTotalTime(scip);
926  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
927  stalltime = timelimit*0.1; /* this should be set as a parameter */
928 
929  for(i = 0; i < g->knots; i++)
930  {
931  g->mark[i] = (g->grad[i] > 0);
932  for(e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e])
933  {
934  random[e] = (double)(rand() % 512);
935  cost[e] = g->cost[e] * 1000.0 + random[e];
936  }
937  }
938 
939  /* this is the offset used to minimise the number of knots to examine in large graphs. */
940  if( g->knots > KNOTLIMIT )
941  {
942  srand(runnum*100);
943  i = 0;
944  do
945  {
946  knotoffset = rand() % KNOTFREQ;
947  i++;
948  } while( g->knots > KNOTLIMIT && knotexamined[knotoffset] >= 0 && i < 50 );
949  knotexamined[knotoffset]++;
950  }
951 
952 
953  for(i = 0; i < g->knots; i++)
954  {
955  if( i % 100 == 0 && elimins == 0 && SCIPgetTotalTime(scip) - redstarttime > stalltime)
956  break;
957 
958  if (!(i % 100))
959  {
960  SCIPdebug(fputc('.', stdout));
961  SCIPdebug(fflush(stdout));
962  }
963  if (g->grad[i] == 0)
964  continue;
965 
966 
967  if( g->knots > KNOTLIMIT && i % KNOTFREQ != knotoffset )
968  continue;
969 
970  /* For the prize collecting variants all edges from the "dummy" root node must be retained. */
971  if ( (g->stp_type == STP_PRIZE_COLLECTING || g->stp_type == STP_ROOTED_PRIZE_COLLECTING
972  || g->stp_type == STP_MAX_NODE_WEIGHT) && i == g->source )
973  continue;
974 
975  for(e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e])
976  {
977  assert(g->mark[g->head[e]] == 1);
978 
979  g->mark[g->head[e]] = 2;
980  }
981 
982  compute_sd(g, i, cost, random, heap, state, &count, sddist, sdtrans, sdrand);
983 
984  for(e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e])
985  {
986  assert(g->mark[g->head[e]] == 2);
987  /* assert(sd[g->head[e]].dist < FARAWAY); */
988 
989  g->mark[g->head[e]] = 1;
990  }
991 
992  for(e = g->outbeg[i]; e != EAT_LAST; e = j)
993  {
994  assert(g->tail[e] == i);
995 
996  j = g->oeat[e];
997 
998  if (LT(g->cost[e], FARAWAY) && LT(sddist[g->head[e]], cost[e])
999  && LT(sddist[g->head[e]] - sdrand[g->head[e]], cost[e] - random[e]))
1000  {
1001  SCIPindexListNodeFree(&((g->ancestors)[e]));
1002  assert(g->ancestors[e] == NULL);
1003  graph_edge_del(g, e);
1004  elimins++;
1005  }
1006  }
1007  }
1008 #if 0
1009  free(heap);
1010  free(state);
1011 
1012  heap = NULL;
1013  state = NULL;
1014 
1015  free(sd);
1016  free(cost);
1017 #endif
1018  assert(graph_valid(g));
1019 
1020  SCIPdebugMessage("%d Edges deleted\n", elimins * 2);
1021  /*printf("%d SD: Edges deleted\n", elimins * 2);*/
1022  return(elimins);
1023 }
1024 #endif
1025 
1026 /** ans subtest */
1027 static
1028 void ansProcessCandidate(
1029  SCIP* scip, /**< SCIP data structure */
1030  GRAPH* g, /**< graph data structure */
1031  int* marked, /**< nodes array */
1032  int* count, /**< pointer to number of reductions */
1033  SCIP_Real min, /**< value to not surpass */
1034  int candvertex /**< candidate */
1035 )
1036 {
1037  int e2;
1038  int bridgeedge = -1;
1039  unsigned misses = 0;
1040 
1041  assert(g->mark[candvertex]);
1042  assert(candvertex != g->source);
1043  assert(!Is_pterm(g->term[candvertex]));
1044 
1045  for( e2 = g->outbeg[candvertex]; e2 != EAT_LAST; e2 = g->oeat[e2] )
1046  {
1047  if( !marked[g->head[e2]] )
1048  {
1049  misses++;
1050  if( misses >= 2 )
1051  break;
1052  bridgeedge = e2;
1053  }
1054  }
1055 
1056  /* neighbors of candvertex subset of those of k? */
1057  if( misses == 0 && SCIPisLE(scip, g->prize[candvertex], min) )
1058  {
1059  (*count) += g->grad[candvertex];
1060  while( g->outbeg[candvertex] != EAT_LAST )
1061  graph_edge_del(scip, g, g->outbeg[candvertex], TRUE);
1062 
1063  g->mark[candvertex] = FALSE;
1064  marked[candvertex] = FALSE;
1065  }
1066  else if( misses == 1 )
1067  {
1068  int e3;
1069  const int neighbor = g->head[bridgeedge];
1070 
1071  if( SCIPisGT(scip, g->prize[neighbor] + g->prize[candvertex], min) )
1072  return;
1073 
1074  for( e3 = g->outbeg[neighbor]; e3 != EAT_LAST; e3 = g->oeat[e3] )
1075  {
1076  const int head = g->head[e3];
1077  if( !marked[head] )
1078  break;
1079  }
1080  if( e3 == EAT_LAST )
1081  {
1082  // delete both vertices?
1083  if( SCIPisLE(scip, g->prize[neighbor], min) && SCIPisLE(scip, g->prize[candvertex], min) )
1084  {
1085  (*count) += g->grad[candvertex] + g->grad[neighbor] - 1;
1086  while( g->outbeg[candvertex] != EAT_LAST )
1087  graph_edge_del(scip, g, g->outbeg[candvertex], TRUE);
1088  while( g->outbeg[neighbor] != EAT_LAST )
1089  graph_edge_del(scip, g, g->outbeg[neighbor], TRUE);
1090 
1091  g->mark[candvertex] = FALSE;
1092  g->mark[neighbor] = FALSE;
1093  marked[candvertex] = FALSE;
1094  marked[neighbor] = FALSE;
1095  }
1096  else
1097  {
1098  graph_edge_del(scip, g, bridgeedge, TRUE);
1099  }
1100  }
1101  }
1102 }
1103 
1104 /** Special distance test */
1105 SCIP_RETCODE reduce_sd(
1106  SCIP* scip, /**< SCIP data structure */
1107  GRAPH* g, /**< graph data structure */
1108  PATH* vnoi, /**< Voronoi data structure */
1109  SCIP_Real* edgepreds, /**< array to store edge predecessors of auxiliary graph */
1110  SCIP_Real* mstsdist, /**< array to store mst distances in auxiliary graph */
1111  int* heap, /**< array representing a heap */
1112  int* state, /**< array to indicate whether a node has been scanned during SP calculation */
1113  int* vbase, /**< Voronoi base to each vertex */
1114  int* nodesid, /**< array to map nodes in auxiliary graph to original ones */
1115  int* nodesorg, /**< array to map terminals of original graph vertices of auxiliary graph */
1116  int* forbidden, /**< array to mark whether an edge may be eliminated */
1117  int* nelims, /**< point to store number of deleted edges */
1118  SCIP_Bool nodereplacing, /**< should node replacement (by edges) be performed? */
1119  int* edgestate /**< array to store status of (directed) edge (for propagation, can otherwise be set to NULL) */
1120  )
1121 {
1122  GRAPH* netgraph;
1123  PATH* mst;
1124  SCIP_Real termdist1[4];
1125  SCIP_Real termdist2[4];
1126  SCIP_Real ecost;
1127  SCIP_Real dist;
1128  int neighbterms1[4];
1129  int neighbterms2[4];
1130  int e;
1131  int i;
1132  int j;
1133  int k;
1134  int l;
1135  int i2;
1136  int tj;
1137  int tk;
1138  int ne;
1139  int nj;
1140  int nk;
1141  int id1;
1142  int id2;
1143  int enext;
1144  int nnodes;
1145  int nterms;
1146  int nedges;
1147  int nnterms1;
1148  int nnterms2;
1149  int maxnedges;
1150  SCIP_Bool checkstate = (edgestate != NULL);
1151 
1152  assert(g != NULL);
1153  assert(scip != NULL);
1154  assert(heap != NULL);
1155  assert(vnoi != NULL);
1156  assert(state != NULL);
1157  assert(vbase != NULL);
1158  assert(nelims != NULL);
1159  assert(nodesid != NULL);
1160  assert(mstsdist != NULL);
1161  assert(nodesorg != NULL);
1162  assert(edgepreds != NULL);
1163  assert(forbidden != NULL);
1164 
1165  nnodes = g->knots;
1166  nterms = g->terms;
1167  nedges = g->edges;
1168  *nelims = 0;
1169  maxnedges = MIN(nedges, (nterms - 1) * nterms);
1170 
1171  /* only one terminal left? */
1172  if( nterms == 1 )
1173  return SCIP_OKAY;
1174 
1175  /* compute nearest four terminals to all non-terminals */
1176  graph_get4nextTerms(scip, g, g->cost, g->cost, vnoi, vbase, heap, state);
1177 
1178  /* construct auxiliary graph to compute paths between terminals */
1179 
1180  /* initialize the new graph */
1181  SCIP_CALL( graph_init(scip, &netgraph, nterms, maxnedges, 1) );
1182 
1183  j = 0;
1184  for( k = 0; k < nnodes; k++ )
1185  {
1186  if( Is_term(g->term[k]) && g->grad[k] > 0 )
1187  {
1188  graph_knot_add(netgraph, -1);
1189  nodesid[k] = j;
1190  mstsdist[j] = -1.0;
1191  netgraph->mark[j] = TRUE;
1192  nodesorg[j++] = k;
1193  }
1194  else
1195  {
1196  nodesid[k] = UNKNOWN;
1197  }
1198  }
1199 
1200  for( k = 0; k < nedges; k++ )
1201  {
1202  forbidden[k] = FALSE;
1203  edgepreds[k] = -1;
1204  }
1205 
1206  assert(netgraph->knots == j);
1207  assert(netgraph->knots == nterms);
1208 
1209  for( k = 0; k < nnodes; k++ )
1210  {
1211  if( g->grad[k] == 0 )
1212  continue;
1213 
1214  i = vbase[k];
1215  id1 = nodesid[i];
1216 
1217  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1218  {
1219  if( i != vbase[g->head[e]] )
1220  {
1221  i2 = vbase[g->head[e]];
1222  id2 = nodesid[i2];
1223 
1224  assert(id1 >= 0);
1225  assert(id2 >= 0);
1226  assert(Is_term(g->term[i]));
1227  assert(Is_term(g->term[i2]));
1228 
1229  for( ne = netgraph->outbeg[id1]; ne != EAT_LAST; ne = netgraph->oeat[ne] )
1230  if( netgraph->head[ne] == id2 )
1231  break;
1232 
1233  /* cost of the edge in the auxiliary graph */
1234  ecost = g->cost[e] + vnoi[g->head[e]].dist + vnoi[g->tail[e]].dist;
1235 
1236  /* does edge already exist? */
1237  if( ne != EAT_LAST )
1238  {
1239  assert(ne >= 0);
1240  assert(netgraph->tail[ne] == id1);
1241  assert(netgraph->head[ne] == id2);
1242 
1243  /* is the new edge better than the existing one? */
1244  if( SCIPisGT(scip, netgraph->cost[ne], ecost) )
1245  {
1246  netgraph->cost[ne] = ecost;
1247  netgraph->cost[flipedge(ne)] = ecost;
1248 
1249  edgepreds[ne] = e;
1250  edgepreds[flipedge(ne)] = flipedge(e);
1251 
1252  assert(ne <= maxnedges);
1253  }
1254  }
1255  else
1256  {
1257  edgepreds[netgraph->edges] = e;
1258  edgepreds[netgraph->edges + 1] = flipedge(e);
1259 
1260  graph_edge_add(scip, netgraph, id1, id2, ecost, ecost);
1261 
1262  assert(netgraph->edges <= maxnedges);
1263  }
1264  }
1265  }
1266  }
1267 
1268  /* compute MST on netgraph */
1269  graph_knot_chg(netgraph, 0, 0);
1270  netgraph->source = 0;
1271 
1272  SCIP_CALL( graph_path_init(scip, netgraph) );
1273  SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
1274 
1275  graph_path_exec(scip, netgraph, MST_MODE, 0, netgraph->cost, mst);
1276 
1277  /* mark (original) edges of MST */
1278  for( k = 1; k < netgraph->knots; k++ )
1279  {
1280  assert(mst[k].edge != -1);
1281  assert((int) edgepreds[mst[k].edge] != -1);
1282  assert((int) edgepreds[flipedge(mst[k].edge)] != -1);
1283 
1284  e = (int) edgepreds[mst[k].edge];
1285 
1286  assert(vbase[g->tail[e]] == nodesorg[k] || vbase[g->head[e]] == nodesorg[k]);
1287 
1288  if( Is_term(g->tail[e]) && Is_term(g->head[e]) )
1289  {
1290  forbidden[e] = TRUE;
1291  forbidden[flipedge(e)] = TRUE;
1292  }
1293  }
1294 
1295  /* traverse all edges */
1296  for( i = 0; i < nnodes; i++ )
1297  {
1298  if( g->grad[i] <= 0 )
1299  continue;
1300 
1301  nnterms1 = 1;
1302  if( Is_term(g->term[i]) )
1303  {
1304  termdist1[0] = 0.0;
1305  neighbterms1[0] = i;
1306  }
1307 
1308  enext = g->outbeg[i];
1309  while( enext != EAT_LAST )
1310  {
1311  e = enext;
1312  i2 = g->head[e];
1313  enext = g->oeat[e];
1314 
1315  if( i2 < i || !g->mark[i2] )
1316  continue;
1317 
1318  ecost = g->cost[e];
1319 
1320  /* is i a terminal? If not, get three closest terminals of distance <= ecost */
1321  if( !Is_term(g->term[i]) )
1322  {
1323  nnterms1 = getlecloseterms(scip, vnoi, termdist1, ecost, vbase, neighbterms1, i, nnodes);
1324 
1325  if( nnterms1 == 0 )
1326  continue;
1327  }
1328 
1329  if( Is_term(g->term[i2]) )
1330  {
1331  nnterms2 = 1;
1332  termdist2[0] = 0.0;
1333  neighbterms2[0] = i2;
1334  }
1335  else
1336  {
1337  /* get closest terminals of distance <= ecost */
1338  nnterms2 = getlecloseterms(scip, vnoi, termdist2, ecost, vbase, neighbterms2, i2, nnodes);
1339 
1340  if( nnterms2 == 0 )
1341  continue;
1342  }
1343 
1344  for( j = 0; j < nnterms1; j++ )
1345  {
1346  /* has edge already been deleted? */
1347  if( g->oeat[e] == EAT_FREE )
1348  break;
1349 
1350  tj = neighbterms1[j];
1351 
1352  assert(tj >= 0);
1353 
1354  for( k = 0; k < nnterms2; k++ )
1355  {
1356  tk = neighbterms2[k];
1357 
1358  assert(tk >= 0);
1359  assert(Is_term(g->term[tk]));
1360  assert(Is_term(g->term[tj]));
1361 
1362  if( tj == tk )
1363  {
1364  if( SCIPisGE(scip, termdist1[j], termdist2[k] ) )
1365  dist = termdist1[j];
1366  else
1367  dist = termdist2[k];
1368 
1369  assert(SCIPisGE(scip, ecost, dist));
1370 
1371  if( SCIPisEQ(scip, dist, ecost) )
1372  if( !sddeltable(scip, g, vnoi, vbase, forbidden, j, k, i, i2, e, nnodes ) )
1373  continue;
1374 
1375  if( checkstate && (edgestate[e] == EDGE_BLOCKED) )
1376  continue;
1377 
1378  graph_edge_del(scip, g, e, TRUE);
1379  (*nelims)++;
1380  break;
1381  }
1382  else
1383  {
1384  /* get sd between (terminals) tj and tk */
1385  nj = nodesid[tj];
1386  nk = nodesid[tk];
1387 
1388  assert(nj != nk);
1389 
1390  l = nj;
1391  dist = 0.0;
1392  mstsdist[l] = 0.0;
1393 
1394  /* go down to the root */
1395  while( l != 0 )
1396  {
1397  ne = mst[l].edge;
1398 
1399  assert(netgraph->head[ne] == l);
1400 
1401  l = netgraph->tail[ne];
1402 
1403  if( SCIPisGT(scip, netgraph->cost[ne], dist) )
1404  dist = netgraph->cost[ne];
1405 
1406  mstsdist[l] = dist;
1407 
1408  if( l == nk )
1409  break;
1410  }
1411 
1412  if( l == nk )
1413  {
1414  l = 0;
1415  }
1416  else
1417  {
1418  l = nk;
1419  dist = 0.0;
1420  }
1421 
1422  while( l != 0 )
1423  {
1424  ne = mst[l].edge;
1425  l = netgraph->tail[ne];
1426 
1427  if( SCIPisGT(scip, netgraph->cost[ne], dist) )
1428  dist = netgraph->cost[ne];
1429 
1430  if( mstsdist[l] >= 0 )
1431  {
1432  if( SCIPisGT(scip, mstsdist[l], dist) )
1433  dist = mstsdist[l];
1434  break;
1435  }
1436  if( l == 0 )
1437  assert(nj == 0);
1438  }
1439 
1440  /* restore */
1441  l = nj;
1442  mstsdist[l] = -1.0;
1443 
1444  while( l != 0 )
1445  {
1446  ne = mst[l].edge;
1447  l = netgraph->tail[ne];
1448  mstsdist[l] = -1.0;
1449  if( l == nk )
1450  break;
1451  }
1452 
1453  assert(SCIPisGT(scip, dist, 0.0));
1454 
1455  if( SCIPisLT(scip, dist, termdist1[j]) )
1456  dist = termdist1[j];
1457 
1458  if( SCIPisLT(scip, dist, termdist2[k]) )
1459  dist = termdist2[k];
1460 
1461  if( SCIPisGE(scip, ecost, dist) )
1462  {
1463  if( SCIPisEQ(scip, ecost, dist) )
1464  if( !(sddeltable(scip, g, vnoi, vbase, forbidden, j, k, i, i2, e, nnodes)) )
1465  continue;
1466 
1467  assert(SCIPisGE(scip, ecost, termdist1[j]));
1468  assert(SCIPisGE(scip, ecost, termdist2[k]));
1469 
1470  if( checkstate && (edgestate[e] == EDGE_BLOCKED) )
1471  continue;
1472 
1473  graph_edge_del(scip, g, e, TRUE);
1474  (*nelims)++;
1475  break;
1476  }
1477  } /* tj != tk (else) */
1478  } /* k < nnterms2 */
1479  } /* j < nnterms1 */
1480 
1481  } /* while( enext != EAT_LAST ) */
1482  }
1483 
1484  if( nodereplacing )
1485  {
1486  SCIP_CALL( reduce_bdr(scip, g, netgraph, mst, vnoi, mstsdist, termdist1, termdist2, vbase, nodesid, neighbterms1, neighbterms2, nelims) );
1487  }
1488 
1489  /* free memory*/
1490  SCIPfreeBufferArray(scip, &mst);
1491  graph_path_exit(scip, netgraph);
1492  graph_free(scip, &netgraph, TRUE);
1493 
1494  return SCIP_OKAY;
1495 }
1496 
1497 
1498 /** SD test for PC */
1499 SCIP_RETCODE reduce_sdPc(
1500  SCIP* scip, /**< SCIP data structure */
1501  GRAPH* g, /**< graph data structure */
1502  PATH* vnoi, /**< Voronoi data structure */
1503  int* heap, /**< heap array */
1504  int* state, /**< array to store state of a node during Voronoi computation*/
1505  int* vbase, /**< Voronoi base to each node */
1506  int* nodesid, /**< array */
1507  int* nodesorg, /**< array */
1508  int* nelims /**< pointer to store number of eliminated edges */
1509  )
1510 {
1511  GRAPH* netgraph;
1512  SCIP_Real termdist1[4];
1513  SCIP_Real termdist2[4];
1514  int neighbterms1[4];
1515  int neighbterms2[4];
1516 
1517  int j;
1518  int maxnedges;
1519  const int root = g->source;
1520  const int nnodes = g->knots;
1521  const int nterms = g->terms;
1522  const int nedges = g->edges;
1523 
1524  assert(g != NULL);
1525  assert(heap != NULL);
1526  assert(vnoi != NULL);
1527  assert(state != NULL);
1528  assert(vbase != NULL);
1529  assert(scip != NULL);
1530  assert(nelims != NULL);
1531  assert(nodesid != NULL);
1532  assert(nodesorg != NULL);
1533 
1534  *nelims = 0;
1535 
1536  if( nterms <= 1 )
1537  return SCIP_OKAY;
1538  else
1539  {
1540  const SCIP_Longint longedges = (SCIP_Longint) nedges;
1541  const SCIP_Longint longtermsq = (SCIP_Longint) (nterms - 1) * nterms;
1542 
1543  if( longedges <= longtermsq )
1544  maxnedges = nedges;
1545  else
1546  maxnedges = ((nterms - 1) * nterms);
1547  }
1548 
1549  /* compute nearest four terminals to each non-terminal */
1550  graph_get4nextTerms(scip, g, g->cost, g->cost, vnoi, vbase, heap, state);
1551 
1552  /*
1553  * construct auxiliary graph to compute paths between terminals
1554  */
1555 
1556  SCIP_CALL( graph_init(scip, &netgraph, nterms, maxnedges, 1) );
1557 
1558  for( int k = 0; k < 4; k++ )
1559  {
1560  termdist1[k] = FARAWAY;
1561  termdist2[k] = FARAWAY;
1562  }
1563 
1564  j = 0;
1565  for( int k = 0; k < nnodes; k++ )
1566  {
1567  if( Is_term(g->term[k]) )
1568  {
1569  assert(g->grad[k] > 0);
1570  graph_knot_add(netgraph, -1);
1571  nodesid[k] = j;
1572  nodesorg[j++] = k;
1573  }
1574  else
1575  {
1576  nodesid[k] = UNKNOWN;
1577  }
1578  }
1579 
1580  assert(netgraph->knots == j);
1581  assert(netgraph->knots == nterms);
1582 
1583  /* insert Voronoi boundary paths as edges into netgraph */
1584  for( int k = 0; k < nnodes; k++ )
1585  {
1586  if( !g->mark[k] )
1587  continue;
1588 
1589  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1590  {
1591  int i;
1592  if( !g->mark[g->head[e]] )
1593  continue;
1594  i = vbase[k];
1595  assert(i != UNKNOWN);
1596  if( i != vbase[g->head[e]] )
1597  {
1598  SCIP_Real ecost;
1599  int ne;
1600  const int i2 = vbase[g->head[e]];
1601 
1602  assert(i2 != UNKNOWN);
1603  assert(Is_term(g->term[i]));
1604  assert(Is_term(g->term[i2]));
1605  assert(nodesid[i] >= 0);
1606  assert(nodesid[i2] >= 0);
1607 
1608  for( ne = netgraph->outbeg[nodesid[i]]; ne != EAT_LAST; ne = netgraph->oeat[ne] )
1609  if( netgraph->head[ne] == nodesid[i2] )
1610  break;
1611 
1612  ecost = g->cost[e] + vnoi[g->head[e]].dist + vnoi[g->tail[e]].dist;
1613 
1614  /* edge exists? */
1615  if( ne != EAT_LAST )
1616  {
1617  assert(ne >= 0);
1618  assert(netgraph->head[ne] == nodesid[i2]);
1619  assert(netgraph->tail[ne] == nodesid[i]);
1620 
1621  if( SCIPisGT(scip, netgraph->cost[ne], ecost) )
1622  {
1623  netgraph->cost[ne] = ecost;
1624  netgraph->cost[flipedge(ne)] = ecost;
1625  assert(ne <= maxnedges);
1626  }
1627  }
1628  else
1629  {
1630  graph_edge_add(scip, netgraph, nodesid[i], nodesid[i2], ecost, ecost);
1631  assert(netgraph->edges <= maxnedges);
1632  }
1633  }
1634  }
1635  }
1636 
1637  /* compute four close terminals to each terminal */
1638  SCIP_CALL( graph_get4nextTTerms(scip, g, g->cost, vnoi, vbase, heap, state) );
1639 
1640  /* traverse all edges */
1641  for( int i = 0; i < nnodes; i++ )
1642  {
1643  int enext;
1644  if( !g->mark[i] )
1645  continue;
1646 
1647  enext = g->outbeg[i];
1648  while( enext != EAT_LAST )
1649  {
1650  SCIP_Real ecost;
1651  int e = enext;
1652  int nnterms1;
1653  int nnterms2;
1654  const int i2 = g->head[e];
1655  enext = g->oeat[e];
1656 
1657  if( i2 < i || Is_term(g->term[i2]) || !g->mark[i2] )
1658  continue;
1659  ecost = g->cost[e];
1660 
1661  /* @todo: fix */
1662 #if 1
1663  if( Is_term(g->term[i]) )
1664  nnterms1 = getcloseterms2term(scip, g, netgraph, termdist1, ecost, neighbterms1, nodesid, nodesorg, i);
1665  else
1666  nnterms1 = getcloseterms(scip, vnoi, termdist1, ecost, vbase, neighbterms1, i, nnodes);
1667 #else
1668  if( Is_term(g->term[i]) )
1669  nnterms1 = 0;
1670  else
1671  nnterms1 = getcloseterms(scip, vnoi, termdist1, ecost, vbase, neighbterms1, i, nnodes);
1672 #endif
1673 
1674  if( nnterms1 == 0 )
1675  continue;
1676 
1677  /* @todo: fix */
1678 #if 1
1679  if( Is_term(g->term[i2]) )
1680  nnterms2 = getcloseterms2term(scip, g, netgraph, termdist2, ecost, neighbterms2, nodesid, nodesorg, i2);
1681  else
1682  nnterms2 = getcloseterms(scip, vnoi, termdist2, ecost, vbase, neighbterms2, i2, nnodes);
1683 #else
1684  if( Is_term(g->term[i2]) )
1685  nnterms2 = 0;
1686  else
1687  nnterms2 = getcloseterms(scip, vnoi, termdist2, ecost, vbase, neighbterms2, i2, nnodes);
1688 #endif
1689 
1690  if( nnterms2 == 0 )
1691  continue;
1692 
1693  /* todo: mark nearest terminals! */
1694  for( j = 0; j < nnterms1; j++ )
1695  {
1696  int tj;
1697 
1698  /* has edge already been deleted? */
1699  if( g->oeat[e] == EAT_FREE )
1700  break;
1701 
1702  tj = neighbterms1[j];
1703 
1704  assert(tj >= 0);
1705  assert(Is_term(g->term[tj]));
1706 
1707  for( int k = 0; k < nnterms2; k++ )
1708  {
1709  const int tk = neighbterms2[k];
1710 
1711  assert(tk >= 0);
1712  assert(Is_term(g->term[tk]));
1713 
1714  if( tj == tk )
1715  {
1716  if( SCIPisGT(scip, ecost, termdist1[j] + termdist2[k] - g->prize[tj]) || tj == root )
1717  {
1718  graph_edge_del(scip, g, e, TRUE);
1719  (*nelims)++;
1720  break;
1721  }
1722  }
1723  else
1724  {
1725  SCIP_Real necost = FARAWAY;
1726  int e2;
1727  int pos;
1728 
1729  /* get distance between terminals */
1730  for( e2 = netgraph->outbeg[nodesid[tj]]; e2 != EAT_LAST; e2 = netgraph->oeat[e2] )
1731  {
1732  if( netgraph->head[e2] == nodesid[tk] )
1733  {
1734  necost = netgraph->cost[e2];
1735  break;
1736  }
1737  }
1738  pos = tj;
1739  for( int l = 0; l < 4; l++ )
1740  {
1741  if( vbase[pos] == UNKNOWN )
1742  break;
1743  if( vbase[pos] == tk && SCIPisLT(scip, vnoi[pos].dist, necost) )
1744  {
1745  necost = vnoi[pos].dist;
1746  break;
1747  }
1748  pos += nnodes;
1749  }
1750 
1751  if( SCIPisGT(scip, ecost, necost)
1752  && SCIPisGT(scip, ecost, necost + termdist1[j] - g->prize[tj])
1753  && SCIPisGT(scip, ecost, necost + termdist2[k] - g->prize[tk])
1754  && SCIPisGT(scip, ecost, necost + termdist1[j] + termdist2[k] - g->prize[tj] - g->prize[tk]) )
1755  {
1756  SCIPdebugMessage("SDPC delete: %d %d (%d)\n", g->tail[e], g->head[e], e);
1757  graph_edge_del(scip, g, e, TRUE);
1758  (*nelims)++;
1759  break;
1760  }
1761  }
1762  }
1763  }
1764  }
1765  }
1766 
1767  SCIPdebugMessage("SDPC eliminations: %d \n", *nelims);
1768  graph_free(scip, &netgraph, TRUE);
1769 
1770  assert(graph_valid(g));
1771 
1772  return SCIP_OKAY;
1773 }
1774 
1775 /** get RSD to a single edge*/
1776 static
1777 SCIP_Real getRSD(
1778  SCIP* scip, /**< SCIP data structure */
1779  GRAPH* g, /**< graph structure */
1780  GRAPH* netgraph, /**< auxiliary graph structure */
1781  PATH* mst, /**< MST structure */
1782  PATH* vnoi, /**< path structure */
1783  SCIP_Real* mstsdist, /**< MST distance in aux-graph */
1784  SCIP_Real* termdist1, /**< dist array */
1785  SCIP_Real* termdist2, /**< second dist array */
1786  SCIP_Real sd_initial, /**< initial sd or -1.0 */
1787  int* vbase, /**< bases for nearest terminals */
1788  int* nodesid, /**< nodes identification array */
1789  int* neighbterms1, /**< neighbour terminals array */
1790  int* neighbterms2, /**< second neighbour terminals array */
1791  int i, /**< first vertex */
1792  int i2, /**< second vertex */
1793  int limit /**< limit for incident edges to consider */
1794  )
1795 {
1796  SCIP_Real sd;
1797  SCIP_Real max;
1798  SCIP_Real dist;
1799  int l;
1800  int e;
1801  int j;
1802  int k;
1803  int ne;
1804  int tj;
1805  int tk;
1806  int nj;
1807  int nk;
1808  int nnodes;
1809  int nnterms1;
1810  int nnterms2;
1811 
1812  assert(scip != NULL);
1813  assert(g != NULL);
1814  assert(netgraph != NULL);
1815  assert(mst != NULL);
1816  assert(mstsdist != NULL);
1817  assert(termdist1 != NULL);
1818  assert(termdist2 != NULL);
1819  assert(neighbterms1 != NULL);
1820  assert(neighbterms2 != NULL);
1821 
1822  nnodes = g->knots;
1823  l = 0;
1824 
1825  if( sd_initial >= 0.0 )
1826  sd = sd_initial;
1827  else
1828  sd = FARAWAY;
1829 
1830  /* compare restricted sd with edge cost (if existing) */
1831  for( e = g->outbeg[i]; (l++ <= limit) && (e != EAT_LAST); e = g->oeat[e] )
1832  {
1833  if( g->head[e] == i2 )
1834  {
1835  if( g->cost[e] < sd )
1836  sd = g->cost[e];
1837  break;
1838  }
1839  }
1840 
1841  /* is i a terminal? If not, get three closest terminals of distance smaller sd */
1842  if( Is_term(g->term[i]) )
1843  {
1844  nnterms1 = 1;
1845  termdist1[0] = 0.0;
1846  neighbterms1[0] = i;
1847  }
1848  else
1849  {
1850  nnterms1 = getcloseterms(scip, vnoi, termdist1, sd, vbase, neighbterms1, i, nnodes);
1851 
1852  if( nnterms1 == 0 )
1853  return sd;
1854  }
1855 
1856  /* is i2 a terminal? If not, get three closest terminals of distance smaller sd */
1857  if( Is_term(g->term[i2]) )
1858  {
1859  nnterms2 = 1;
1860  termdist2[0] = 0.0;
1861  neighbterms2[0] = i2;
1862  }
1863  else
1864  {
1865  /* get closest terminals of distance smaller sd */
1866  nnterms2 = getcloseterms(scip, vnoi, termdist2, sd, vbase, neighbterms2, i2, nnodes);
1867 
1868  if( nnterms2 == 0 )
1869  return sd;
1870  }
1871 
1872  for( j = 0; j < nnterms1; j++ )
1873  {
1874  tj = neighbterms1[j];
1875  assert(tj >= 0);
1876  for( k = 0; k < nnterms2; k++ )
1877  {
1878  tk = neighbterms2[k];
1879  assert(Is_term(g->term[tk]));
1880  assert(Is_term(g->term[tj]));
1881  assert(tk >= 0);
1882 
1883  if( SCIPisGT(scip, termdist1[j], termdist2[k]) )
1884  max = termdist1[j];
1885  else
1886  max = termdist2[k];
1887 
1888  if( tj == tk )
1889  {
1890  if( SCIPisLT(scip, max, sd) )
1891  sd = max;
1892  }
1893  else
1894  {
1895  /* get sd between (terminals) tj and tk */
1896  nj = nodesid[tj];
1897  nk = nodesid[tk];
1898  assert(nj != nk);
1899 
1900  l = nj;
1901  dist = 0.0;
1902  mstsdist[l] = 0.0;
1903 
1904  while( l != 0 )
1905  {
1906  ne = mst[l].edge;
1907 
1908  assert(netgraph->head[ne] == l);
1909  l = netgraph->tail[ne];
1910  if( SCIPisGT(scip, netgraph->cost[ne], dist) )
1911  dist = netgraph->cost[ne];
1912 
1913  mstsdist[l] = dist;
1914  if( l == nk )
1915  break;
1916  }
1917 
1918  if( l == nk )
1919  {
1920  l = 0;
1921  }
1922  else
1923  {
1924  l = nk;
1925  dist = 0.0;
1926  }
1927  while( l != 0 )
1928  {
1929  ne = mst[l].edge;
1930  l = netgraph->tail[ne];
1931  if( SCIPisGT(scip, netgraph->cost[ne], dist) )
1932  dist = netgraph->cost[ne];
1933 
1934  if( mstsdist[l] >= 0 )
1935  {
1936  if( SCIPisGT(scip, mstsdist[l], dist) )
1937  dist = mstsdist[l];
1938  break;
1939  }
1940  if( l == 0 )
1941  assert( nj == 0);
1942  }
1943 
1944  /* restore */
1945  l = nj;
1946  mstsdist[l] = -1.0;
1947  while( l != 0 )
1948  {
1949  ne = mst[l].edge;
1950  l = netgraph->tail[ne];
1951  mstsdist[l] = -1.0;
1952  if( l == nk )
1953  break;
1954  }
1955 
1956  assert(SCIPisGT(scip, dist, 0.0));
1957  if( SCIPisGT(scip, dist, max) )
1958  max = dist;
1959  if( SCIPisLT(scip, max, sd) )
1960  sd = max;
1961  }
1962  } /* k < nnterms2 */
1963  } /* j < nnterms1 */
1964  return sd;
1965 }
1966 
1967 
1968 
1969 /** get SD to a single edge*/
1970 SCIP_RETCODE reduce_getSd(
1971  SCIP* scip,
1972  GRAPH* g,
1973  PATH* pathtail,
1974  PATH* pathhead,
1975  SCIP_Real* sdist,
1976  SCIP_Real distlimit,
1977  int* heap,
1978  int* statetail,
1979  int* statehead,
1980  int* memlbltail,
1981  int* memlblhead,
1982  int i,
1983  int i2,
1984  int limit,
1985  SCIP_Bool pc,
1986  SCIP_Bool mw
1987  )
1988 {
1989  SCIP_Real sd;
1990  SCIP_Real dist;
1991  int k;
1992  int l;
1993  int e;
1994  int nnodes;
1995  int nlbltail;
1996  int nlblhead;
1997  const SCIP_Bool pcmw = (pc || mw);
1998 
1999  assert(g != NULL);
2000  assert(scip != NULL);
2001  assert(pathtail != NULL);
2002  assert(pathhead != NULL);
2003  assert(heap != NULL);
2004  assert(statetail != NULL);
2005  assert(statehead != NULL);
2006  assert(memlbltail != NULL);
2007  assert(memlblhead != NULL);
2008  assert(sdist != NULL);
2009 
2010  nnodes = g->knots;
2011 
2012  /* start from tail */
2013  graph_sdPaths(scip, g, pathtail, g->cost, distlimit, heap, statetail, memlbltail, &nlbltail, i, i2, limit);
2014 
2015  /* test whether edge e can be eliminated */
2016  graph_sdPaths(scip, g, pathhead, g->cost, distlimit, heap, statehead, memlblhead, &nlblhead, i2, i, limit);
2017 
2018  sd = FARAWAY;
2019 #if 0
2020  if( statetail[i2] != UNKNOWN )
2021  {
2022  sd = pathtail[i2].dist;
2023  assert(SCIPisGT(scip, FARAWAY, sd));
2024  }
2025  if( statehead[i] != UNKNOWN && SCIPisGT(scip, sd, pathhead[i].dist) )
2026  {
2027  sd = pathhead[i].dist;
2028  assert(SCIPisGT(scip, FARAWAY, sd));
2029  }
2030 #endif
2031  /* get restore state and path of tail and head */
2032  for( k = 0; k < nlbltail; k++ )
2033  {
2034  l = memlbltail[k];
2035  assert(statetail[l] != UNKNOWN);
2036  if( statehead[l] != UNKNOWN )
2037  {
2038  assert(SCIPisGT(scip, FARAWAY, pathtail[l].dist));
2039  assert(SCIPisGT(scip, FARAWAY, pathhead[l].dist));
2040  if( Is_term(g->term[l]) )
2041  {
2042  dist = 0.0;
2043  if( SCIPisLT(scip, dist, pathhead[l].dist) )
2044  dist = pathhead[l].dist;
2045  if( SCIPisLT(scip, dist, pathtail[l].dist) )
2046  dist = pathtail[l].dist;
2047  if( pcmw && SCIPisLT(scip, dist, pathhead[l].dist + pathtail[l].dist - g->prize[l]) )
2048  dist = pathhead[l].dist + pathtail[l].dist - g->prize[l];
2049  if( SCIPisGT(scip, sd, dist) )
2050  sd = dist;
2051  }
2052  else
2053  {
2054  if( mw && l != i && l != i2 )
2055  assert(SCIPisLE(scip, g->prize[l], 0.0));
2056  if( mw && SCIPisLT(scip, g->prize[l], 0.0) )
2057  dist = pathhead[l].dist + pathtail[l].dist + g->prize[l];
2058  else
2059  dist = pathhead[l].dist + pathtail[l].dist;
2060  if( SCIPisGT(scip, sd, dist) )
2061  sd = dist;
2062  }
2063  }
2064 
2065  statetail[l] = UNKNOWN;
2066  pathtail[l].dist = FARAWAY;
2067  pathtail[l].edge = UNKNOWN;
2068  }
2069  /* restore state and path of tail and head */
2070  for( k = 0; k < nlblhead; k++ )
2071  {
2072  l = memlblhead[k];
2073  statehead[l] = UNKNOWN;
2074  pathhead[l].dist = FARAWAY;
2075  pathhead[l].edge = UNKNOWN;
2076  }
2077 
2078 
2079  for( k = 0; k < nnodes; k++ )
2080  {
2081  assert(statetail[k] == UNKNOWN);
2082  assert(pathtail[k].dist == FARAWAY);
2083  assert(pathtail[k].edge == UNKNOWN);
2084  assert(statehead[k] == UNKNOWN);
2085  assert(pathhead[k].dist == FARAWAY);
2086  assert(pathhead[k].edge == UNKNOWN);
2087  }
2088 
2089 
2090  l = 0;
2091  /* compare restricted sd with edge cost (if existing) */
2092  for( e = g->outbeg[i]; (l++ <= limit) && (e != EAT_LAST); e = g->oeat[e] )
2093  {
2094  if( g->head[e] == i2 )
2095  {
2096  if( mw )
2097  sd = 0.0;
2098  else if( SCIPisGT(scip, sd, g->cost[e]) )
2099  sd = g->cost[e];
2100  break;
2101  }
2102  }
2103 
2104  *sdist = sd;
2105  return SCIP_OKAY;
2106 }
2107 
2108 
2109 /** get SD to a single edge*/
2110 SCIP_RETCODE reduce_getSdPcMw(
2111  SCIP* scip,
2112  const GRAPH* g,
2113  PATH* pathtail,
2114  PATH* pathhead,
2115  SCIP_Real* sdist,
2116  SCIP_Real distlimit,
2117  int* heap,
2118  int* statetail,
2119  int* statehead,
2120  int* memlbltail,
2121  int* memlblhead,
2122  int* pathmaxnodetail,
2123  int* pathmaxnodehead,
2124  int i,
2125  int i2,
2126  int limit
2127  )
2128 {
2129  SCIP_Real sd;
2130  int nlbltail;
2131  int nlblhead;
2132  const int nnodes = g->knots;
2133  const SCIP_Bool mw = g->stp_type == STP_MWCSP;
2134 
2135  assert((g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG) || mw);
2136  assert(g != NULL);
2137  assert(scip != NULL);
2138  assert(pathtail != NULL);
2139  assert(pathhead != NULL);
2140  assert(heap != NULL);
2141  assert(statetail != NULL);
2142  assert(statehead != NULL);
2143  assert(memlbltail != NULL);
2144  assert(memlblhead != NULL);
2145  assert(pathmaxnodetail != NULL);
2146  assert(pathmaxnodehead != NULL);
2147  assert(sdist != NULL);
2148 
2149  graph_path_PcMwSd(scip, g, pathtail, g->cost, distlimit, pathmaxnodetail, heap, statetail, NULL, memlbltail, &nlbltail, i, i2, limit);
2150  graph_path_PcMwSd(scip, g, pathhead, g->cost, distlimit, pathmaxnodehead, heap, statehead, statetail, memlblhead, &nlblhead, i2, i, limit);
2151 
2152  sd = FARAWAY;
2153 
2154  /* get restore state and path of tail and head */
2155  for( int k = 0; k < nlbltail; k++ )
2156  {
2157  const int l = memlbltail[k];
2158  assert(statetail[l] != UNKNOWN);
2159 
2160  if( statehead[l] != UNKNOWN )
2161  {
2162  SCIP_Real dist = FARAWAY;
2163  const int tailmaxterm = pathmaxnodetail[l];
2164  const int headmaxterm = pathmaxnodehead[l];
2165 
2166  assert(SCIPisGT(scip, FARAWAY, pathtail[l].dist));
2167  assert(SCIPisGT(scip, FARAWAY, pathhead[l].dist));
2168  assert(tailmaxterm != i && headmaxterm != i);
2169  assert(tailmaxterm != i2 && headmaxterm != i2);
2170 
2171  /* any terminal on the path? */
2172  if( tailmaxterm >= 0 || headmaxterm >= 0 )
2173  {
2174  if( tailmaxterm == headmaxterm )
2175  {
2176  assert(tailmaxterm == l);
2177  assert(SCIPisPositive(scip, g->prize[tailmaxterm]));
2178 
2179  dist = misc_stp_maxReal((SCIP_Real []) {
2180  pathhead[headmaxterm].dist,
2181  pathtail[tailmaxterm].dist,
2182  pathhead[l].dist + pathtail[l].dist - g->prize[l]
2183  }, 3);
2184  SCIPdebugMessage("sd1 %f \n", dist);
2185  }
2186  else if( tailmaxterm >= 0 && headmaxterm >= 0 )
2187  {
2188  const SCIP_Real distl2tailmax = pathtail[l].dist - pathtail[tailmaxterm].dist;
2189  const SCIP_Real distl2headmax = pathhead[l].dist - pathhead[headmaxterm].dist;
2190 
2191  assert(tailmaxterm != headmaxterm);
2192  assert(!SCIPisNegative(scip, distl2tailmax));
2193  assert(!SCIPisNegative(scip, distl2headmax));
2194  assert(SCIPisPositive(scip, g->prize[tailmaxterm]) && SCIPisPositive(scip, g->prize[headmaxterm]));
2195 
2196  dist = misc_stp_maxReal((SCIP_Real []) {
2197  pathhead[headmaxterm].dist,
2198  pathtail[tailmaxterm].dist,
2199  distl2tailmax + distl2headmax,
2200  distl2tailmax + pathhead[l].dist - g->prize[headmaxterm],
2201  distl2headmax + pathtail[l].dist - g->prize[tailmaxterm],
2202  pathhead[l].dist + pathtail[l].dist - g->prize[tailmaxterm] - g->prize[headmaxterm]
2203  }, 6);
2204  SCIPdebugMessage("sd2 %f \n", dist);
2205  }
2206  else if( tailmaxterm >= 0 )
2207  {
2208  const SCIP_Real distl2tailmax = pathtail[l].dist - pathtail[tailmaxterm].dist;
2209 
2210  assert(headmaxterm < 0);
2211  assert(SCIPisPositive(scip, g->prize[tailmaxterm]));
2212 
2213  dist = misc_stp_maxReal((SCIP_Real []) {
2214  pathtail[tailmaxterm].dist,
2215  distl2tailmax + pathhead[l].dist,
2216  pathhead[l].dist + pathtail[l].dist - g->prize[tailmaxterm]
2217  }, 3);
2218  SCIPdebugMessage("sd3 %f \n", dist);
2219  }
2220  else if( headmaxterm >= 0 )
2221  {
2222  const SCIP_Real distl2headmax = pathhead[l].dist - pathhead[headmaxterm].dist;
2223 
2224  assert(tailmaxterm < 0);
2225  assert(SCIPisPositive(scip, g->prize[headmaxterm]));
2226 
2227  dist = misc_stp_maxReal((SCIP_Real []) {
2228  pathhead[headmaxterm].dist,
2229  distl2headmax + pathtail[l].dist,
2230  pathhead[l].dist + pathtail[l].dist - g->prize[headmaxterm]
2231  }, 3);
2232  SCIPdebugMessage("sd4 %f \n", dist);
2233  }
2234  }
2235  else
2236  {
2237  dist = pathhead[l].dist + pathtail[l].dist;
2238  }
2239 
2240  if( dist < sd )
2241  sd = dist;
2242 
2243  if( Is_term(g->term[l]) )
2244  {
2245  dist = misc_stp_maxReal((SCIP_Real []) {
2246  pathhead[l].dist,
2247  pathtail[l].dist,
2248  pathhead[l].dist + pathtail[l].dist - g->prize[l]
2249  }, 3);
2250  if( dist < sd )
2251  sd = dist;
2252  }
2253  }
2254  }
2255 
2256  /* restore state and path of tail and head */
2257 
2258  for( int k = 0; k < nlbltail; k++ )
2259  {
2260  const int l = memlbltail[k];
2261  statetail[l] = UNKNOWN;
2262  pathtail[l].dist = FARAWAY;
2263  pathtail[l].edge = UNKNOWN;
2264  pathmaxnodetail[l] = -1;
2265  }
2266 
2267  for( int k = 0; k < nlblhead; k++ )
2268  {
2269  const int l = memlblhead[k];
2270  statehead[l] = UNKNOWN;
2271  pathhead[l].dist = FARAWAY;
2272  pathhead[l].edge = UNKNOWN;
2273  pathmaxnodehead[l] = -1;
2274  }
2275 
2276  for( int k = 0; k < nnodes; k++ )
2277  {
2278  assert(statetail[k] == UNKNOWN);
2279  assert(pathtail[k].dist == FARAWAY);
2280  assert(pathtail[k].edge == UNKNOWN);
2281  assert(statehead[k] == UNKNOWN);
2282  assert(pathhead[k].dist == FARAWAY);
2283  assert(pathhead[k].edge == UNKNOWN);
2284  assert(pathmaxnodehead[k] == -1);
2285  assert(pathmaxnodetail[k] == -1);
2286  }
2287 
2288  /* compare restricted sd with edge cost (if existing) */
2289  for( int e = g->outbeg[i], count = 0; (count++ <= limit) && (e != EAT_LAST); e = g->oeat[e] )
2290  {
2291  if( g->head[e] == i2 )
2292  {
2293  if( mw )
2294  sd = 0.0;
2295  else if( sd > g->cost[e] )
2296  sd = g->cost[e];
2297  break;
2298  }
2299  }
2300 
2301  *sdist = sd;
2302 
2303  return SCIP_OKAY;
2304 }
2305 
2306 
2307 /** SDC test for the SAP using a limited version of Dijkstra's algorithm from both endpoints of an arc */
2308 SCIP_RETCODE reduce_sdspSap(
2309  SCIP* scip,
2310  GRAPH* g,
2311  PATH* pathtail,
2312  PATH* pathhead,
2313  int* heap,
2314  int* statetail,
2315  int* statehead,
2316  int* memlbltail,
2317  int* memlblhead,
2318  int* nelims,
2319  int limit
2320  )
2321 {
2322  SCIP_Real sdist;
2323  SCIP_Real* costrev;
2324  int i;
2325  int k;
2326  int l;
2327  int e;
2328  int i2;
2329  int enext;
2330  int nnodes;
2331  int nedges;
2332  int nlbltail;
2333  int nlblhead;
2334 
2335  assert(g != NULL);
2336  assert(scip != NULL);
2337  assert(pathtail != NULL);
2338  assert(pathhead != NULL);
2339  assert(heap != NULL);
2340  assert(statetail != NULL);
2341  assert(statehead != NULL);
2342  assert(memlbltail != NULL);
2343  assert(memlblhead != NULL);
2344  assert(nelims != NULL);
2345 
2346  nedges = g->edges;
2347  nnodes = g->knots;
2348  *nelims = 0;
2349 
2350  for( i = 0; i < nnodes; i++ )
2351  {
2352  g->mark[i] = (g->grad[i] > 0);
2353  statetail[i] = UNKNOWN;
2354  pathtail[i].dist = FARAWAY;
2355  pathtail[i].edge = UNKNOWN;
2356  statehead[i] = UNKNOWN;
2357  pathhead[i].dist = FARAWAY;
2358  pathhead[i].edge = UNKNOWN;
2359  }
2360 
2361  SCIP_CALL( SCIPallocBufferArray(scip, &costrev, nedges) );
2362 
2363  for( e = 0; e < nedges; e++ )
2364  costrev[e] = g->cost[flipedge(e)];
2365 
2366  /* iterate through all nodes */
2367  for( i = 0; i < nnodes; i++ )
2368  {
2369  if( !g->mark[i] )
2370  continue;
2371  /* traverse neighbours */
2372  e = g->outbeg[i];
2373  while( e != EAT_LAST )
2374  {
2375 
2376  enext = g->oeat[e];
2377  i2 = g->head[e];
2378 
2379  assert(g->mark[i2]);
2380 
2381  /* start limited dijkstra from i, marking all reached vertices */
2382  graph_sdPaths(scip, g, pathtail, g->cost, g->cost[e], heap, statetail, memlbltail, &nlbltail, i, i2, limit);
2383 
2384  /* start limited dijkstra from i2, marking all reached vertices */
2385  graph_sdPaths(scip, g, pathhead, costrev, g->cost[e], heap, statehead, memlblhead, &nlblhead, i2, i, limit);
2386 
2387  sdist = FARAWAY;
2388 #if 0
2389  if( statetail[i2] != UNKNOWN )
2390  {
2391  sdist = pathtail[i2].dist;
2392  assert(SCIPisGT(scip, FARAWAY, sdist));
2393  }
2394  if( statehead[i] != UNKNOWN && SCIPisGT(scip, sdist, pathhead[i].dist) )
2395  {
2396  sdist = pathhead[i].dist;
2397  assert(SCIPisGT(scip, FARAWAY, sdist));
2398  }
2399 #endif
2400  /* get restore state and path of tail and head */
2401  for( k = 0; k < nlbltail; k++ )
2402  {
2403  l = memlbltail[k];
2404  assert(g->mark[l]);
2405  assert(statetail[l] != UNKNOWN);
2406  if( statehead[l] != UNKNOWN )
2407  {
2408  assert(SCIPisGT(scip, FARAWAY, pathtail[l].dist));
2409  assert(SCIPisGT(scip, FARAWAY, pathhead[l].dist));
2410 
2411  if( SCIPisGT(scip, sdist, pathhead[l].dist + pathtail[l].dist) )
2412  sdist = pathhead[l].dist + pathtail[l].dist;
2413  }
2414 
2415  statetail[l] = UNKNOWN;
2416  pathtail[l].dist = FARAWAY;
2417  pathtail[l].edge = UNKNOWN;
2418  }
2419  /* restore state and path of tail and head */
2420  for( k = 0; k < nlblhead; k++ )
2421  {
2422  l = memlblhead[k];
2423  statehead[l] = UNKNOWN;
2424  pathhead[l].dist = FARAWAY;
2425  pathhead[l].edge = UNKNOWN;
2426  }
2427 
2428 #if 1
2429  for( k = 0; k < nnodes; k++ )
2430  {
2431  assert(statetail[k] == UNKNOWN);
2432  assert(pathtail[k].dist == FARAWAY);
2433  assert(pathtail[k].edge == UNKNOWN);
2434  assert(statehead[k] == UNKNOWN);
2435  assert(pathhead[k].dist == FARAWAY);
2436  assert(pathhead[k].edge == UNKNOWN);
2437  }
2438 #endif
2439  /* can edge be deleted? */
2440  if( SCIPisGE(scip, g->cost[e], sdist) )
2441  {
2442  if( SCIPisGE(scip, costrev[e], FARAWAY) )
2443  graph_edge_del(scip, g, e, TRUE);
2444  else
2445  g->cost[e] = FARAWAY;
2446  (*nelims)++;
2447  }
2448 
2449  e = enext;
2450  }
2451  }
2452 
2453  SCIPfreeBufferArray(scip, &costrev);
2454 
2455  return SCIP_OKAY;
2456 }
2457 
2458 /** SD test using only limited dijkstra from both endpoints of an edge */
2459 SCIP_RETCODE reduce_sdsp(
2460  SCIP* scip,
2461  GRAPH* g,
2462  PATH* pathtail,
2463  PATH* pathhead,
2464  int* heap,
2465  int* statetail,
2466  int* statehead,
2467  int* memlbltail,
2468  int* memlblhead,
2469  int* nelims,
2470  int limit,
2471  int* edgestate
2472 
2473  )
2474 {
2475  int* pathmaxnodetail = NULL;
2476  int* pathmaxnodehead = NULL;
2477  const int nnodes = g->knots;
2478  const SCIP_Bool pc = (g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG);
2479  const SCIP_Bool checkstate = (edgestate != NULL);
2480 
2481  assert(g != NULL);
2482  assert(scip != NULL);
2483  assert(pathtail != NULL);
2484  assert(pathhead != NULL);
2485  assert(heap != NULL);
2486  assert(statetail != NULL);
2487  assert(statehead != NULL);
2488  assert(memlbltail != NULL);
2489  assert(memlblhead != NULL);
2490  assert(nelims != NULL);
2491 
2492  *nelims = 0;
2493 
2494  if( pc )
2495  {
2496  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodetail, nnodes) );
2497  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodehead, nnodes) );
2498 
2499  for( int i = 0; i < nnodes; i++ )
2500  {
2501  pathmaxnodetail[i] = -1;
2502  pathmaxnodehead[i] = -1;
2503  }
2504  }
2505  else
2506  {
2507  for( int i = 0; i < nnodes; i++ )
2508  g->mark[i] = (g->grad[i] > 0);
2509  }
2510 
2511  for( int i = 0; i < nnodes; i++ )
2512  {
2513  statetail[i] = UNKNOWN;
2514  pathtail[i].dist = FARAWAY;
2515  pathtail[i].edge = UNKNOWN;
2516  statehead[i] = UNKNOWN;
2517  pathhead[i].dist = FARAWAY;
2518  pathhead[i].edge = UNKNOWN;
2519  }
2520 
2521  /* iterate through all nodes */
2522  for( int i = 0; i < nnodes; i++ )
2523  {
2524  int e;
2525  if( !g->mark[i] )
2526  continue;
2527 
2528  /* traverse neighbours */
2529  e = g->outbeg[i];
2530  while( e != EAT_LAST )
2531  {
2532  const SCIP_Real ecost = g->cost[e];
2533  const int i2 = g->head[e];
2534  const int enext = g->oeat[e];
2535  int nlbltail;
2536  int nlblhead;
2537  SCIP_Bool deletable;
2538 
2539  /* avoid double checking */
2540  if( i2 < i || !g->mark[i2] )
2541  {
2542  e = enext;
2543  continue;
2544  }
2545 
2546  /* execute limited Dijkstra from both sides */
2547 
2548  if( pc )
2549  {
2550  graph_path_PcMwSd(scip, g, pathtail, g->cost, ecost, pathmaxnodetail, heap, statetail, NULL, memlbltail, &nlbltail, i, i2, limit);
2551  graph_path_PcMwSd(scip, g, pathhead, g->cost, ecost, pathmaxnodehead, heap, statehead, statetail, memlblhead, &nlblhead, i2, i, limit);
2552  }
2553  else
2554  {
2555  graph_sdPaths(scip, g, pathtail, g->cost, ecost, heap, statetail, memlbltail, &nlbltail, i, i2, limit);
2556  graph_sdPaths(scip, g, pathhead, g->cost, ecost, heap, statehead, memlblhead, &nlblhead, i2, i, limit);
2557  }
2558 
2559  deletable = FALSE;
2560 
2561  /* check whether edge e can be deleted and restore data structures */
2562  for( int k = 0; k < nlbltail && !deletable; k++ )
2563  {
2564  const int l = memlbltail[k];
2565 
2566  assert(g->mark[l]);
2567  assert(statetail[l] != UNKNOWN);
2568 
2569  if( statehead[l] != UNKNOWN )
2570  {
2571  assert(SCIPisGT(scip, FARAWAY, pathtail[l].dist));
2572  assert(SCIPisGT(scip, FARAWAY, pathhead[l].dist));
2573 
2574  if( pc )
2575  {
2576  const int tailmaxterm = pathmaxnodetail[l];
2577  const int headmaxterm = pathmaxnodehead[l];
2578 
2579  assert(tailmaxterm != i && headmaxterm != i);
2580  assert(tailmaxterm != i2 && headmaxterm != i2);
2581 
2582  /* any terminal on the path? */
2583  if( tailmaxterm >= 0 || headmaxterm >= 0 )
2584  {
2585  if( tailmaxterm == headmaxterm )
2586  {
2587  assert(tailmaxterm == l);
2588  assert(SCIPisPositive(scip, g->prize[tailmaxterm]));
2589  assert(SCIPisGE(scip, ecost, pathhead[headmaxterm].dist) && SCIPisGE(scip, ecost, pathtail[tailmaxterm].dist));
2590 
2591  if( SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[l]) )
2592  {
2593  deletable = TRUE;
2594  SCIPdebugMessage("delete1Term \n");
2595  }
2596  }
2597  else if( tailmaxterm >= 0 && headmaxterm >= 0 )
2598  {
2599  const SCIP_Real distl2tailmax = pathtail[l].dist - pathtail[tailmaxterm].dist;
2600  const SCIP_Real distl2headmax = pathhead[l].dist - pathhead[headmaxterm].dist;
2601 
2602  assert(tailmaxterm != headmaxterm);
2603  assert(!SCIPisNegative(scip, distl2tailmax));
2604  assert(!SCIPisNegative(scip, distl2headmax));
2605  assert(SCIPisPositive(scip, g->prize[tailmaxterm]) && SCIPisPositive(scip, g->prize[headmaxterm]));
2606  assert(SCIPisGE(scip, ecost, pathhead[headmaxterm].dist) && SCIPisGE(scip, ecost, pathtail[tailmaxterm].dist));
2607 
2608  if( SCIPisGE(scip, ecost, distl2tailmax + distl2headmax)
2609  && SCIPisGE(scip, ecost, distl2tailmax + pathhead[l].dist - g->prize[headmaxterm])
2610  && SCIPisGE(scip, ecost, distl2headmax + pathtail[l].dist - g->prize[tailmaxterm])
2611  && SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[tailmaxterm] - g->prize[headmaxterm]) )
2612  {
2613  deletable = TRUE;
2614  SCIPdebugMessage("delete2Term \n");
2615  }
2616  }
2617  else if( tailmaxterm >= 0 )
2618  {
2619  const SCIP_Real distl2tailmax = pathtail[l].dist - pathtail[tailmaxterm].dist;
2620  // todo consider l == term?
2621  assert(headmaxterm < 0);
2622  assert(SCIPisGE(scip, ecost, pathtail[tailmaxterm].dist));
2623  assert(SCIPisPositive(scip, g->prize[tailmaxterm]));
2624 
2625  if( SCIPisGE(scip, ecost, distl2tailmax + pathhead[l].dist)
2626  && SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[tailmaxterm]) )
2627  {
2628  deletable = TRUE;
2629  SCIPdebugMessage("deleteHalfTerm1 \n");
2630  }
2631  }
2632  else if( headmaxterm >= 0 )
2633  {
2634  const SCIP_Real distl2headmax = pathhead[l].dist - pathhead[headmaxterm].dist;
2635  // todo consider l == term?
2636  assert(tailmaxterm < 0);
2637  assert(SCIPisGE(scip, ecost, pathhead[headmaxterm].dist));
2638  assert(SCIPisPositive(scip, g->prize[headmaxterm]));
2639 
2640  if( SCIPisGE(scip, ecost, distl2headmax + pathtail[l].dist)
2641  && SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[headmaxterm]) )
2642  {
2643  deletable = TRUE;
2644  SCIPdebugMessage("deleteHalfTerm2 \n");
2645  }
2646  }
2647  }
2648  else if( SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist) )
2649  {
2650  deletable = TRUE;
2651  }
2652 
2653  if( Is_term(g->term[l]) )
2654  {
2655  if( SCIPisGE(scip, ecost, pathhead[l].dist) && SCIPisGE(scip, ecost, pathtail[l].dist)
2656  && SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[l]) )
2657  deletable = TRUE;
2658  }
2659  }
2660  else
2661  {
2662  if( Is_term(g->term[l]) )
2663  {
2664  if( SCIPisGE(scip, ecost, pathhead[l].dist) && SCIPisGE(scip, ecost, pathtail[l].dist) )
2665  {
2666  deletable = TRUE;
2667 #if 0
2668  if( pc && SCIPisLT(scip, ecost, pathhead[l].dist + pathtail[l].dist - g->prize[l]) )
2669  deletable = FALSE;
2670 #endif
2671  }
2672  }
2673  else
2674  {
2675  if( SCIPisGE(scip, ecost, pathhead[l].dist + pathtail[l].dist) )
2676  deletable = TRUE;
2677  }
2678  }
2679  }
2680  }
2681 
2682  /* restore data */
2683 
2684  for( int k = 0; k < nlbltail; k++ )
2685  {
2686  const int l = memlbltail[k];
2687  statetail[l] = UNKNOWN;
2688  pathtail[l].dist = FARAWAY;
2689  pathtail[l].edge = UNKNOWN;
2690  if( pc )
2691  pathmaxnodetail[l] = -1;
2692  }
2693 
2694  for( int k = 0; k < nlblhead; k++ )
2695  {
2696  const int l = memlblhead[k];
2697  statehead[l] = UNKNOWN;
2698  pathhead[l].dist = FARAWAY;
2699  pathhead[l].edge = UNKNOWN;
2700  if( pc )
2701  pathmaxnodehead[l] = -1;
2702  }
2703 
2704 #ifndef NDEBUG
2705  for( int k = 0; k < nnodes; k++ )
2706  {
2707  assert(statetail[k] == UNKNOWN);
2708  assert(pathtail[k].dist == FARAWAY);
2709  assert(pathtail[k].edge == UNKNOWN);
2710  assert(statehead[k] == UNKNOWN);
2711  assert(pathhead[k].dist == FARAWAY);
2712  assert(pathhead[k].edge == UNKNOWN);
2713  if( pc )
2714  {
2715  assert(pathmaxnodetail[k] == -1);
2716  assert(pathmaxnodehead[k] == -1);
2717  }
2718  }
2719 #endif
2720  /* can edge be deleted? */
2721  if( deletable )
2722  {
2723  if( !checkstate || (edgestate[e] != EDGE_BLOCKED) )
2724  {
2725  graph_edge_del(scip, g, e, TRUE);
2726  (*nelims)++;
2727  }
2728  }
2729 
2730  e = enext;
2731  }
2732  }
2733 
2734  SCIPfreeBufferArrayNull(scip, &pathmaxnodehead);
2735  SCIPfreeBufferArrayNull(scip, &pathmaxnodetail);
2736 
2737  assert(graph_valid(g));
2738 
2739  return SCIP_OKAY;
2740 }
2741 
2742 
2743 #define STP_BDR_MAXDEGREE 4
2744 #define STP_BDR_MAXDNEDGES 6
2745 
2746 /** bd_k test for given Steiner bottleneck distances todo bd5 */
2747 SCIP_RETCODE reduce_bdr(
2748  SCIP* scip, /**< SCIP data structure */
2749  GRAPH* g, /**< graph structure */
2750  GRAPH* netgraph, /**< auxiliary graph structure */
2751  PATH* netmst, /**< MST structure */
2752  PATH* vnoi, /**< path structure */
2753  SCIP_Real* mstsdist, /**< MST distance in aux-graph */
2754  SCIP_Real* termdist1, /**< dist array */
2755  SCIP_Real* termdist2, /**< second dist array */
2756  int* vbase, /**< bases for nearest terminals */
2757  int* nodesid, /**< nodes identification array */
2758  int* neighbterms1, /**< neighbour terminals array */
2759  int* neighbterms2, /**< second neighbour terminals array */
2760  int* nelims /**< number of eliminations */
2761  )
2762 {
2763  SCIP_Real cutoffs[STP_BDR_MAXDNEDGES];
2764  SCIP_Real sd[STP_BDR_MAXDEGREE];
2765  SCIP_Real ecost[STP_BDR_MAXDEGREE];
2766  int adjvert[STP_BDR_MAXDEGREE];
2767  GRAPH* auxg;
2768  PATH* mst;
2769  SCIP_Real csum;
2770  SCIP_Real mstcost;
2771  const int nnodes = g->knots;
2772  const SCIP_Bool pc = g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG;
2773 
2774  assert(g != NULL);
2775  assert(netgraph != NULL);
2776  assert(netmst != NULL);
2777 
2778  /* initialize mst struct and new graph for bd4 test */
2779  SCIP_CALL( SCIPallocBufferArray(scip, &mst, STP_BDR_MAXDEGREE) );
2780  SCIP_CALL( graph_init(scip, &auxg, STP_BDR_MAXDEGREE, 2 * STP_BDR_MAXDNEDGES, 1) );
2781 
2782  for( int k = 0; k < STP_BDR_MAXDEGREE; k++ )
2783  graph_knot_add(auxg, -1);
2784 
2785  for( int k = 0; k < STP_BDR_MAXDEGREE - 1; k++ )
2786  for( int k2 = STP_BDR_MAXDEGREE - 1; k2 >= k + 1; k2-- )
2787  graph_edge_add(scip, auxg, k, k2, 1.0, 1.0);
2788 
2789  assert(auxg->edges == 2 * STP_BDR_MAXDNEDGES);
2790 
2791  /* init graph for mst computation */
2792  SCIP_CALL( graph_path_init(scip, auxg) );
2793 
2794  SCIPdebugMessage("BD3-R Reduction: ");
2795 
2796  for( int i = 0; i < STP_BDR_MAXDEGREE; i++ )
2797  sd[i] = 0.0;
2798 
2799  if( !pc )
2800  for( int i = 0; i < nnodes; i++ )
2801  g->mark[i] = (g->grad[i] > 0);
2802 
2803  for( int degree = 3; degree <= STP_BDR_MAXDEGREE; degree ++ )
2804  {
2805  for( int i = 0; i < nnodes; i++ )
2806  {
2807  int k;
2808 
2809  if( Is_term(g->term[i]) || g->grad[i] != degree )
2810  continue;
2811 
2812  assert(g->mark[i]);
2813 
2814  k = 0;
2815  for( int e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
2816  {
2817  ecost[k] = g->cost[e];
2818  adjvert[k++] = g->head[e];
2819  }
2820 
2821  assert(k == degree);
2822 
2823  /* vertex of degree 3? */
2824  if( degree == 3 )
2825  {
2826  csum = ecost[0] + ecost[1] + ecost[2];
2827 
2828  sd[0] = getRSD(scip, g, netgraph, netmst, vnoi, mstsdist, termdist1, termdist2, ecost[0] + ecost[1], vbase, nodesid, neighbterms1, neighbterms2, adjvert[0],
2829  adjvert[1], 300);
2830  sd[1] = getRSD(scip, g, netgraph, netmst, vnoi, mstsdist, termdist1, termdist2, ecost[1] + ecost[2], vbase, nodesid, neighbterms1, neighbterms2, adjvert[1],
2831  adjvert[2], 300);
2832  sd[2] = getRSD(scip, g, netgraph, netmst, vnoi, mstsdist, termdist1, termdist2, ecost[2] + ecost[0], vbase, nodesid, neighbterms1, neighbterms2, adjvert[2],
2833  adjvert[0], 300);
2834 
2835  if( sd[0] + sd[1] <= csum || sd[0] + sd[2] <= csum || sd[1] + sd[2] <= csum )
2836  {
2837  SCIP_Bool success;
2838 
2839  cutoffs[0] = sd[0];
2840  cutoffs[1] = sd[2];
2841  cutoffs[2] = sd[1];
2842 
2843  SCIP_CALL(graph_knot_delPseudo(scip, g, g->cost, cutoffs, NULL, i, &success));
2844 
2845  assert(success);
2846  assert(g->grad[i] == 0);
2847 
2848  SCIPdebugMessage("BD3-R Reduction: %f %f %f csum: %f\n ", sd[0], sd[1], sd[2], csum);
2849  (*nelims)++;
2850  }
2851  }
2852  /* vertex of degree 4? */
2853  else if( degree == 4 )
2854  {
2855  for( k = 0; k < 4; k++ )
2856  {
2857  auxg->mark[k] = TRUE;
2858  for( int e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
2859  {
2860  const int k2 = auxg->head[e];
2861  if( k2 > k )
2862  {
2863  auxg->cost[e] = getRSD(scip, g, netgraph, netmst, vnoi, mstsdist, termdist1, termdist2, ecost[k] + ecost[k2], vbase, nodesid, neighbterms1,
2864  neighbterms2, adjvert[k], adjvert[k2], 200);
2865  auxg->cost[flipedge(e)] = auxg->cost[e];
2866  }
2867  }
2868  }
2869 
2870  for( int l = 0; l < 4; l++ )
2871  mst[l].dist = UNKNOWN;
2872 
2873  /* compute mst on all neighbours */
2874  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
2875  mstcost = 0.0;
2876 
2877 #ifndef NDEBUG
2878  for( int l = 1; l < 4; l++ )
2879  assert(mst[l].dist != UNKNOWN);
2880 #endif
2881 
2882  for( int l = 1; l < 4; l++ )
2883  mstcost += mst[l].dist;
2884 
2885  k = UNKNOWN;
2886  csum = ecost[0] + ecost[1] + ecost[2] + ecost[3];
2887 
2888  if( csum >= mstcost )
2889  {
2890  /* compute mst on all 3-subsets of all neigbours */
2891  for( k = 0; k < 4; k++ )
2892  {
2893  auxg->mark[k] = FALSE;
2894  mstcost = 0.0;
2895 
2896  if( k != 0 )
2897  {
2898  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
2899  for( int l = 1; l < 4; l++ )
2900  if( auxg->mark[l] )
2901  mstcost += mst[l].dist;
2902  }
2903  else
2904  {
2905  graph_path_exec(scip, auxg, MST_MODE, 1, auxg->cost, mst);
2906  for( int l = 2; l < 4; l++ )
2907  mstcost += mst[l].dist;
2908  }
2909 
2910  auxg->mark[k] = TRUE;
2911  csum -= ecost[k];
2912 
2913  if( csum < mstcost )
2914  break;
2915 
2916  csum += ecost[k];
2917  }
2918  }
2919 
2920  if( k == 4 )
2921  {
2922 
2923  int edgecount = 0;
2924  SCIP_Bool success;
2925  for( k = 0; k < 3; k++ )
2926  {
2927  for( int e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
2928  {
2929  const int k2 = auxg->head[e];
2930  if( k2 > k )
2931  cutoffs[edgecount++] = auxg->cost[e];
2932  }
2933  }
2934 
2935  SCIP_CALL(graph_knot_delPseudo(scip, g, g->cost, cutoffs, NULL, i, &success));
2936 
2937  if( success )
2938  (*nelims)++;
2939  }
2940  }
2941  }
2942  }
2943 
2944  graph_path_exit(scip, auxg);
2945  graph_free(scip, &auxg, TRUE);
2946  SCIPfreeBufferArray(scip, &mst);
2947 
2948  return SCIP_OKAY;
2949 }
2950 
2951 
2952 #define STP_BD_MAXDEGREE 4
2953 #define STP_BD_MAXDNEDGES 6
2954 
2955 
2956 /* C. W. Duin and A. Volganant
2957  *
2958  * "Reduction Tests for the Steiner Problem in Graphs"
2959  *
2960  * Networks, Volume 19 (1989), 549-567
2961  *
2962  * Bottleneck Degree 3,4 Test
2963  */
2964 SCIP_RETCODE reduce_bd34(
2965  SCIP* scip, /**< SCIP data structure */
2966  GRAPH* g, /**< graph structure */
2967  PATH* pathtail, /**< array for internal use */
2968  PATH* pathhead, /**< array for internal use */
2969  int* heap, /**< array for internal use */
2970  int* statetail, /**< array for internal use */
2971  int* statehead, /**< array for internal use */
2972  int* memlbltail, /**< array for internal use */
2973  int* memlblhead, /**< array for internal use */
2974  int* nelims, /**< point to return number of eliminations */
2975  int limit /**< limit for edges to consider for each vertex */
2976  )
2977 {
2978  SCIP_Real cutoffs[STP_BD_MAXDNEDGES];
2979  PATH mst[STP_BD_MAXDEGREE];
2980  SCIP_Real sd[STP_BD_MAXDEGREE];
2981  SCIP_Real ecost[STP_BD_MAXDEGREE];
2982  GRAPH* auxg;
2983  int* pathmaxnodetail = NULL;
2984  int* pathmaxnodehead = NULL;
2985  int adjvert[STP_BD_MAXDEGREE];
2986  const int nnodes = g->knots;
2987  const SCIP_Bool pc = g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG;
2988 
2989  SCIPdebugMessage("BD34-Reduction: ");
2990 
2991  assert(scip != NULL);
2992  assert(g != NULL);
2993  assert(heap != NULL);
2994  assert(nelims != NULL);
2995 
2996  /* initialize new graph for bd4 tests */
2997  SCIP_CALL( graph_init(scip, &auxg, STP_BD_MAXDEGREE, 2 * STP_BD_MAXDNEDGES, 1) );
2998 
2999  for( int k = 0; k < STP_BD_MAXDEGREE; k++ )
3000  graph_knot_add(auxg, -1);
3001 
3002  for( int k = 0; k < STP_BD_MAXDEGREE; k++ )
3003  for( int k2 = STP_BD_MAXDEGREE - 1; k2 >= k + 1; k2-- )
3004  graph_edge_add(scip, auxg, k, k2, 1.0, 1.0);
3005 
3006  /* init graph for mst computation */
3007  SCIP_CALL( graph_path_init(scip, auxg) );
3008 
3009  *nelims = 0;
3010 
3011  for( int i = 0; i < STP_BD_MAXDEGREE; i++ )
3012  sd[i] = 0.0;
3013 
3014  if( pc )
3015  {
3016  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodetail, nnodes) );
3017  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodehead, nnodes) );
3018  }
3019  else
3020  {
3021  for( int i = 0; i < nnodes; i++ )
3022  g->mark[i] = (g->grad[i] > 0);
3023  }
3024 
3025  for( int i = 0; i < nnodes; i++ )
3026  {
3027  statetail[i] = UNKNOWN;
3028  pathtail[i].dist = FARAWAY;
3029  pathtail[i].edge = UNKNOWN;
3030  statehead[i] = UNKNOWN;
3031  pathhead[i].dist = FARAWAY;
3032  pathhead[i].edge = UNKNOWN;
3033  if( pc )
3034  {
3035  pathmaxnodetail[i] = -1;
3036  pathmaxnodehead[i] = -1;
3037  }
3038  }
3039 
3040  for( int degree = 3; degree <= STP_BD_MAXDEGREE; degree++ )
3041  {
3042  for( int i = 0; i < nnodes; i++ )
3043  {
3044  int edgecount;
3045  if( Is_term(g->term[i]) || g->grad[i] != degree )
3046  continue;
3047 
3048  edgecount = 0;
3049  for( int e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
3050  {
3051  ecost[edgecount] = g->cost[e];
3052  adjvert[edgecount++] = g->head[e];
3053  }
3054 
3055  /* vertex of degree 3? */
3056  if( degree == 3 )
3057  {
3058  const SCIP_Real csum = ecost[0] + ecost[1] + ecost[2];
3059 
3060  assert(edgecount == 3);
3061 
3062  if( pc )
3063  {
3064  SCIP_CALL(
3065  reduce_getSdPcMw(scip, g, pathtail, pathhead, &(sd[0]), csum, heap, statetail, statehead, memlbltail, memlblhead,
3066  pathmaxnodetail, pathmaxnodehead, adjvert[0], adjvert[1], limit));
3067  SCIP_CALL(
3068  reduce_getSdPcMw(scip, g, pathtail, pathhead, &(sd[1]), csum, heap, statetail, statehead, memlbltail, memlblhead,
3069  pathmaxnodetail, pathmaxnodehead, adjvert[1], adjvert[2], limit));
3070  SCIP_CALL(
3071  reduce_getSdPcMw(scip, g, pathtail, pathhead, &(sd[2]), csum, heap, statetail, statehead, memlbltail, memlblhead,
3072  pathmaxnodetail, pathmaxnodehead, adjvert[2], adjvert[0], limit));
3073  }
3074  else
3075  {
3076  SCIP_CALL(
3077  reduce_getSd(scip, g, pathtail, pathhead, &(sd[0]), csum, heap, statetail, statehead, memlbltail, memlblhead, adjvert[0], adjvert[1], limit, pc, FALSE));
3078  SCIP_CALL(
3079  reduce_getSd(scip, g, pathtail, pathhead, &(sd[1]), csum, heap, statetail, statehead, memlbltail, memlblhead, adjvert[1], adjvert[2], limit, pc, FALSE));
3080  SCIP_CALL(
3081  reduce_getSd(scip, g, pathtail, pathhead, &(sd[2]), csum, heap, statetail, statehead, memlbltail, memlblhead, adjvert[2], adjvert[0], limit, pc, FALSE));
3082  }
3083 
3084  if( SCIPisLE(scip, sd[0] + sd[1], csum) || SCIPisLE(scip, sd[0] + sd[2], csum) || SCIPisLE(scip, sd[1] + sd[2], csum) )
3085  {
3086  SCIP_Bool success;
3087 
3088  cutoffs[0] = sd[0];
3089  cutoffs[1] = sd[2];
3090  cutoffs[2] = sd[1];
3091 
3092  SCIP_CALL(graph_knot_delPseudo(scip, g, g->cost, cutoffs, NULL, i, &success));
3093  assert(success);
3094  assert(g->grad[i] == 0);
3095 
3096  SCIPdebugMessage("BD3 Reduction: %f %f %f csum: %f\n ", sd[0], sd[1], sd[2], csum);
3097  (*nelims)++;
3098  }
3099  }
3100  /* vertex of degree 4? */
3101  else if( degree == 4 )
3102  {
3103  int k;
3104  SCIP_Real csum = ecost[0] + ecost[1] + ecost[2] + ecost[3];
3105  SCIP_Real mstcost;
3106 
3107  assert(edgecount == 4);
3108 
3109  for( k = 0; k < 4; k++ )
3110  {
3111  auxg->mark[k] = TRUE;
3112  for( int e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
3113  {
3114  const int k2 = auxg->head[e];
3115  if( k2 > k )
3116  {
3117  SCIP_Real s1 = -1.0;
3118  if( pc )
3119  SCIP_CALL(
3120  reduce_getSdPcMw(scip, g, pathtail, pathhead, &(s1), csum, heap, statetail, statehead, memlbltail, memlblhead,
3121  pathmaxnodetail, pathmaxnodehead, adjvert[k], adjvert[k2], limit));
3122  else
3123  SCIP_CALL(
3124  reduce_getSd(scip, g, pathtail, pathhead, &(s1), csum, heap, statetail, statehead, memlbltail, memlblhead, adjvert[k], adjvert[k2], limit, pc, FALSE));
3125  assert(s1 >= 0);
3126  auxg->cost[e] = s1;
3127  auxg->cost[flipedge(e)] = s1;
3128  }
3129  }
3130  }
3131 
3132  for( int l = 0; l < 4; l++ )
3133  mst[l].dist = UNKNOWN;
3134 
3135  /* compute mst on all neighbours */
3136  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
3137  mstcost = 0.0;
3138 
3139  /* sum cost of (root == 0) */
3140  for( int l = 1; l < 4; l++ )
3141  {
3142  assert(mst[l].dist != UNKNOWN);
3143  mstcost += mst[l].dist;
3144  }
3145 
3146  k = 0;
3147  if( SCIPisGE(scip, csum, mstcost) )
3148  {
3149  /* compute mst on all 3-subsets of all neighbors */
3150  for( k = 0; k < 4; k++ )
3151  {
3152  auxg->mark[k] = FALSE;
3153  mstcost = 0.0;
3154 
3155  if( k != 0 )
3156  {
3157  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
3158  for( int l = 1; l < 4; l++ )
3159  if( auxg->mark[l] )
3160  mstcost += mst[l].dist;
3161  }
3162  else
3163  {
3164  graph_path_exec(scip, auxg, MST_MODE, 1, auxg->cost, mst);
3165  for( int l = 2; l < 4; l++ )
3166  mstcost += mst[l].dist;
3167  }
3168 
3169  auxg->mark[k] = TRUE;
3170 
3171  if( SCIPisLT(scip, csum - ecost[k], mstcost) )
3172  break;
3173  }
3174  }
3175 
3176  if( k == 4 )
3177  {
3178  SCIP_Bool success;
3179  edgecount = 0;
3180 
3181  for( k = 0; k < 3; k++ )
3182  {
3183  for( int e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
3184  {
3185  const int k2 = auxg->head[e];
3186  if( k2 > k )
3187  cutoffs[edgecount++] = auxg->cost[e];
3188  }
3189  }
3190 
3191  SCIP_CALL(graph_knot_delPseudo(scip, g, g->cost, cutoffs, NULL, i, &success));
3192 
3193  if( success )
3194  (*nelims)++;
3195  }
3196  }
3197  }
3198  }
3199 
3200  /* free memory */
3201 
3202  SCIPfreeBufferArrayNull(scip, &pathmaxnodehead);
3203  SCIPfreeBufferArrayNull(scip, &pathmaxnodetail);
3204 
3205  graph_path_exit(scip, auxg);
3206  graph_free(scip, &auxg, TRUE);
3207 
3208  SCIPdebugMessage("bd34: %d nodes deleted\n", *nelims);
3209 
3210  assert(graph_valid(g));
3211 
3212  return SCIP_OKAY;
3213 }
3214 
3215 
3216 /** Non-Terminal-Set test*/
3217 SCIP_RETCODE reduce_nts(
3218  SCIP* scip, /**< SCIP data structure */
3219  GRAPH* g, /**< graph structure */
3220  PATH* pathtail, /**< array for internal use */
3221  PATH* pathhead, /**< array for internal use */
3222  int* heap, /**< array for internal use */
3223  int* statetail, /**< array for internal use */
3224  int* statehead, /**< array for internal use */
3225  int* memlbltail, /**< array for internal use */
3226  int* memlblhead, /**< array for internal use */
3227  int* nelims, /**< point to return number of eliminations */
3228  int limit /**< limit for edges to consider for each vertex */
3229  )
3230 {
3231  PATH mst[STP_BD_MAXDEGREE];
3232  SCIP_Real ecost[STP_BD_MAXDEGREE];
3233  int outedge[STP_BD_MAXDEGREE];
3234  int adjvert[STP_BD_MAXDEGREE];
3235  GRAPH* auxg;
3236  int* pathmaxnodetail = NULL;
3237  int* pathmaxnodehead = NULL;
3238  const int nnodes = g->knots;
3239  const int nedges = g->edges;
3240  const SCIP_Bool pc = g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG;
3241 
3242  SCIPdebugMessage("NTS-Reduction: ");
3243 
3244  assert(scip != NULL);
3245  assert(g != NULL);
3246  assert(heap != NULL);
3247  assert(nelims != NULL);
3248  assert(limit > 0);
3249 
3250  // todo extra method, also called from bd34
3251 
3252  /* initialize new graph for nts tests */
3253  SCIP_CALL( graph_init(scip, &auxg, STP_BD_MAXDEGREE, 2 * STP_BD_MAXDNEDGES, 1) );
3254 
3255  for( int k = 0; k < STP_BD_MAXDEGREE; k++ )
3256  {
3257  graph_knot_add(auxg, -1);
3258  auxg->mark[k] = TRUE;
3259  }
3260 
3261  for( int k = 0; k < STP_BD_MAXDEGREE; k++ )
3262  for( int k2 = STP_BD_MAXDEGREE - 1; k2 >= k + 1; k2-- )
3263  graph_edge_add(scip, auxg, k, k2, 1.0, 1.0);
3264 
3265  /* init graph for mst computation */
3266  SCIP_CALL( graph_path_init(scip, auxg) );
3267 
3268  *nelims = 0;
3269 
3270  if( pc )
3271  {
3272  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodetail, nnodes) );
3273  SCIP_CALL( SCIPallocBufferArray(scip, &pathmaxnodehead, nnodes) );
3274  }
3275  else
3276  {
3277  for( int i = 0; i < nnodes; i++ )
3278  g->mark[i] = (g->grad[i] > 0);
3279  }
3280 
3281  for( int i = 0; i < nnodes; i++ )
3282  {
3283  statetail[i] = UNKNOWN;
3284  pathtail[i].dist = FARAWAY;
3285  pathtail[i].edge = UNKNOWN;
3286  statehead[i] = UNKNOWN;
3287  pathhead[i].dist = FARAWAY;
3288  pathhead[i].edge = UNKNOWN;
3289  if( pc )
3290  {
3291  pathmaxnodetail[i] = -1;
3292  pathmaxnodehead[i] = -1;
3293  }
3294  }
3295 
3296  for( int degree = 4; degree <= STP_BD_MAXDEGREE; degree++ )
3297  {
3298  for( int edge = 0; edge < nedges; edge += 2 )
3299  {
3300  /* edge not deleted? */
3301  if( g->oeat[edge] != EAT_FREE )
3302  {
3303  SCIP_Real csum;
3304  SCIP_Real mstcost;
3305  const SCIP_Real edgecost = g->cost[edge];
3306  int edgecount;
3307  const int tail = g->tail[edge];
3308  const int head = g->head[edge];
3309  SCIP_Bool success = TRUE;
3310  SCIP_Bool duplicate;
3311 
3312  assert(edge >= 0);
3313  assert(tail >= 0 && head >= 0);
3314 
3315  // todo
3316  if( g->grad[tail] != 3 || g->grad[head] != 3 )
3317  continue;
3318 
3319  if( Is_term(g->term[tail]) || Is_term(g->term[head]) )
3320  continue;
3321 
3322  edgecount = 0;
3323  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
3324  {
3325  if( g->head[e] == head )
3326  continue;
3327  outedge[edgecount] = e;
3328  ecost[edgecount] = g->cost[e];
3329  adjvert[edgecount++] = g->head[e];
3330  }
3331 
3332  for( int e = g->outbeg[head]; e != EAT_LAST; e = g->oeat[e] )
3333  {
3334  if( g->head[e] == tail )
3335  continue;
3336  outedge[edgecount] = e;
3337  ecost[edgecount] = g->cost[e];
3338  adjvert[edgecount++] = g->head[e];
3339  }
3340 
3341  assert(edgecount == degree);
3342 
3343  duplicate = FALSE;
3344 
3345  /* check for shared neighbor */
3346  for( int i = 0; i < degree - 1; i++ )
3347  for( int j = i + 1; j < degree; j++ )
3348  if( adjvert[i] == adjvert[j] && adjvert[j] >= 0 )
3349  {
3350  adjvert[j] = -adjvert[j] - 1;
3351  duplicate = TRUE;
3352  }
3353 
3354  // todo reshuffle
3355  if( duplicate )
3356  {
3357  continue;
3358  }
3359 
3360  assert(g->mark[head] && g->mark[tail]);
3361 
3362  g->mark[head] = FALSE;
3363  g->mark[tail] = FALSE;
3364 
3365  csum = ecost[0] + ecost[1] + ecost[2] + ecost[3];
3366 
3367  /* compute sd values and check 2-subsets of neighbors */
3368  for( int k = 0; k < degree - 1 && success; k++ )
3369  {
3370  for( int e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
3371  {
3372  const int k2 = auxg->head[e];
3373  if( k2 > k )
3374  {
3375  SCIP_Real s1 = -1.0;
3376 
3377  if( pc )
3378  SCIP_CALL( reduce_getSdPcMw(scip, g, pathtail, pathhead, &(s1), csum, heap, statetail, statehead, memlbltail, memlblhead,
3379  pathmaxnodetail, pathmaxnodehead, adjvert[k], adjvert[k2], limit) );
3380  else
3381  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &(s1), csum, heap, statetail, statehead, memlbltail, memlblhead,
3382  adjvert[k], adjvert[k2], limit, pc, FALSE) );
3383 
3384  assert(s1 >= 0);
3385 
3386  if( g->tail[outedge[k]] != g->tail[outedge[k2]] )
3387  {
3388  const SCIP_Real innercost = ecost[k] + ecost[k2] + edgecost;
3389 
3390  if( SCIPisGT(scip, s1, innercost) )
3391  {
3392  // success = FALSE;
3393  // break;
3394  }
3395  }
3396 
3397  auxg->cost[e] = s1;
3398  auxg->cost[flipedge(e)] = s1;
3399  }
3400  }
3401  }
3402 
3403  g->mark[head] = TRUE;
3404  g->mark[tail] = TRUE;
3405 
3406  for( int l = 0; l < 4; l++ )
3407  mst[l].dist = UNKNOWN;
3408 
3409  /* compute MST on all neighbors */
3410  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
3411  mstcost = 0.0;
3412 
3413  /* sum cost of (root == 0) */
3414  for( int l = 1; l < 4; l++ )
3415  {
3416  assert(mst[l].dist != UNKNOWN);
3417  mstcost += mst[l].dist;
3418  }
3419 
3420  success = (success && SCIPisGE(scip, csum, mstcost));
3421 
3422 
3423  if( success && 0)
3424  {
3425  /* compute MST on all 3-subsets of all neighbors */
3426  for( int k = 0; k < degree; k++ )
3427  {
3428  assert(auxg->mark[k]);
3429 
3430  auxg->mark[k] = FALSE;
3431  mstcost = 0.0;
3432 
3433  if( k != 0 )
3434  {
3435  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
3436  for( int l = 1; l < 4; l++ )
3437  if( auxg->mark[l] )
3438  mstcost += mst[l].dist;
3439  }
3440  else
3441  {
3442  graph_path_exec(scip, auxg, MST_MODE, 1, auxg->cost, mst);
3443  for( int l = 2; l < 4; l++ )
3444  mstcost += mst[l].dist;
3445  }
3446 
3447  auxg->mark[k] = TRUE;
3448 
3449  if( SCIPisLT(scip, csum - ecost[k], mstcost) )
3450  {
3451  success = FALSE;
3452  break;
3453  }
3454  }
3455  }
3456 
3457  if( success )
3458  {
3459  // graph_edge_del(scip, g, edge, TRUE);
3460 
3461  (*nelims)++;
3462  }
3463  } /* check single edge */
3464  } /* for all edges */
3465  } /* for all degrees */
3466 
3467 
3468  /* free memory */
3469 
3470  SCIPfreeBufferArrayNull(scip, &pathmaxnodehead);
3471  SCIPfreeBufferArrayNull(scip, &pathmaxnodetail);
3472 
3473  graph_path_exit(scip, auxg);
3474  graph_free(scip, &auxg, TRUE);
3475 
3476  assert(graph_valid(g));
3477 
3478  return SCIP_OKAY;
3479 }
3480 
3481 
3482 /* shortest link reduction */
3483 SCIP_RETCODE reduce_sl(
3484  SCIP* scip, /**< SCIP data structure */
3485  GRAPH* g, /**< graph data structure */
3486  PATH* vnoi, /**< Voronoi data structure */
3487  double* fixed, /**< offset pointer */
3488  int* heap, /**< heap array */
3489  int* state, /**< shortest path array */
3490  int* vbase, /**< Voronoi/shortest path bases array */
3491  int* vrnodes, /**< Voronoi/shortest path array */
3492  STP_Bool* visited, /**< Voronoi/shortest path array */
3493  int* solnode, /**< node array to mark whether an node is part of a given solution (CONNECT),
3494  or NULL */
3495  int* nelims /**< pointer to store number of eliminations */
3496  )
3497 {
3498  SCIP_QUEUE* queue;
3499  SCIP_Real cost;
3500  SCIP_Real mincost2;
3501  SCIP_Real mincost3;
3502  int i;
3503  int k;
3504  int e;
3505  int j;
3506  int t;
3507  int old;
3508  int head;
3509  int tail;
3510  int root;
3511  int nnodes;
3512  int vrcount;
3513  int minedge;
3514  int* qnode;
3515  STP_Bool contract;
3516  STP_Bool foundterms;
3517  STP_Bool* forbidden;
3518  STP_Bool* newterm;
3519  SCIP_Bool pc;
3520 
3521  assert(g != NULL);
3522  assert(vnoi != NULL);
3523  assert(heap != NULL);
3524  assert(state != NULL);
3525  assert(vbase != NULL);
3526  assert(vrnodes != NULL);
3527  assert(visited != NULL);
3528 
3529  *nelims = 0;
3530  foundterms = FALSE;
3531  nnodes = g->knots;
3532  root = g->source;
3533  pc = (g->stp_type == STP_PCSPG) || (g->stp_type == STP_RPCSPG);
3534 
3535  if( g->terms <= 1 )
3536  return SCIP_OKAY;
3537 
3538  SCIP_CALL( SCIPallocBufferArray(scip, &forbidden, nnodes) );
3539  SCIP_CALL( SCIPallocBufferArray(scip, &newterm, nnodes) );
3540 
3541  if( !pc )
3542  for( i = 0; i < nnodes; i++ )
3543  g->mark[i] = (g->grad[i] > 0);
3544 
3545  graph_voronoiTerms(scip, g, g->cost, vnoi, vbase, heap, state);
3546 
3547  SCIP_CALL( SCIPqueueCreate(&queue, nnodes, 2.0) );
3548  for( j = 0; j < nnodes; j++ )
3549  {
3550  newterm[j] = FALSE;
3551  forbidden[j] = FALSE;
3552  visited[j] = FALSE;
3553  }
3554  for( i = 0; i < nnodes; i++ )
3555  {
3556  /* is i terminal and not disabled? */
3557  if( Is_term(g->term[i]) && g->mark[i] && !forbidden[i] )
3558  {
3559  /* traverse voronoi-region of (terminal) i */
3560  assert(SCIPqueueIsEmpty(queue));
3561  t = i;
3562  SCIP_CALL( SCIPqueueInsert(queue, &t) );
3563  vrcount = 1;
3564  vrnodes[0] = i;
3565  visited[i] = TRUE;
3566  minedge = UNKNOWN;
3567  mincost2 = FARAWAY;
3568  mincost3 = FARAWAY;
3569 
3570  while( !SCIPqueueIsEmpty(queue) )
3571  {
3572  qnode = (SCIPqueueRemove(queue));
3573  /* traverse all adjacent edges */
3574  for( e = g->outbeg[*qnode]; e != EAT_LAST; e = g->oeat[e] )
3575  {
3576  j = g->head[e];
3577 
3578  if( !g->mark[j] )
3579  continue;
3580 
3581  k = vbase[j];
3582  assert(k != UNKNOWN);
3583  if( !visited[j] && k == i )
3584  {
3585  visited[j] = TRUE;
3586  vrnodes[vrcount++] = j;
3587  SCIP_CALL( SCIPqueueInsert(queue, &(g->head[e])) );
3588  }
3589  else if( k != i )
3590  /* update shortest and second shortest edge (cost) leaving the voronoi region */
3591  {
3592  cost = g->cost[e];
3593  if( minedge == UNKNOWN )
3594  {
3595  minedge = e;
3596  }
3597  else if( SCIPisLE(scip, cost, g->cost[minedge]) )
3598  {
3599  mincost3 = mincost2;
3600  mincost2 = g->cost[minedge];
3601  minedge = e;
3602  }
3603  else if( SCIPisLE(scip, cost, mincost2) )
3604  {
3605  mincost3 = mincost2;
3606  mincost2 = g->cost[e];
3607  }
3608  else if( SCIPisLT(scip, cost, mincost3) )
3609  {
3610  mincost3 = g->cost[e];
3611  }
3612  }
3613  }
3614  }
3615  for( j = 0; j < vrcount; j++ )
3616  visited[vrnodes[j]] = FALSE;
3617  if( minedge == UNKNOWN )
3618  continue;
3619  e = minedge;
3620  tail = g->tail[e];
3621  head = g->head[e];
3622  assert(vbase[tail] == i);
3623 
3624  contract = FALSE;
3625  cost = vnoi[tail].dist + g->cost[e] + vnoi[head].dist;
3626  if( SCIPisGE(scip, mincost2, cost) )
3627  {
3628  contract = TRUE;
3629  }
3630 
3631  /* check whether minedge can be removed */
3632  if( contract )
3633  {
3634  if( pc )
3635  {
3636  if( root != vbase[head] && !SCIPisLE(scip, g->cost[e] + vnoi[tail].dist + vnoi[head].dist, g->prize[vbase[head]]) )
3637  continue;
3638  if( i == tail )
3639  {
3640  if( root != i && !SCIPisLE(scip, vnoi[tail].dist + g->cost[e], g->prize[i]) )
3641  continue;
3642  }
3643  else
3644  {
3645  if( root != i && !SCIPisLT(scip, vnoi[tail].dist + g->cost[e], g->prize[i]) )
3646  continue;
3647  }
3648  if( Is_term(g->term[head]) && Is_term(g->term[tail]) )
3649  continue;
3650  }
3651 
3652  *fixed += g->cost[e];
3653  assert(g->mark[tail] && g->mark[head]);
3654  assert(!Is_pterm(g->term[tail]) && !Is_pterm(g->term[head]));
3655 
3656  if( Is_term(g->term[head]) )
3657  {
3658  j = head;
3659  k = tail;
3660  }
3661  else
3662  {
3663  j = tail;
3664  k = head;
3665  }
3666 
3667  old = g->grad[j] + g->grad[k] - 1;
3668 
3669  if( pc )
3670  {
3671  SCIP_CALL( graph_pc_contractEdge(scip, g, solnode, j, k, i) );
3672  }
3673  else
3674  {
3675  SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->fixedges), g->ancestors[e], NULL) );
3676  SCIP_CALL( graph_knot_contract(scip, g, solnode, j, k) );
3677 
3678  assert(g->grad[k] == 0 && g->grad[j] >= 0);
3679 
3680  if( !Is_term(g->term[j]) )
3681  {
3682  newterm[j] = TRUE;
3683  foundterms = TRUE;
3684  }
3685  }
3686 
3687  assert(old - g->grad[j] - g->grad[k] > 0);
3688  (*nelims) += old - g->grad[j] - g->grad[k];
3689  forbidden[vbase[j]] = TRUE;
3690  forbidden[vbase[k]] = TRUE;
3691  }
3692  }
3693  }
3694 
3695  for( i = 0; i < nnodes && foundterms; i++ )
3696  if( newterm[i] && !Is_term(g->term[i]) && g->grad[i] > 0 )
3697  graph_knot_chg(g, i, 0);
3698 
3699 
3700  /* free memory */
3701  SCIPqueueFree(&queue);
3702  SCIPfreeBufferArray(scip, &newterm);
3703  SCIPfreeBufferArray(scip, &forbidden);
3704 
3705  return SCIP_OKAY;
3706 }
3707 
3708 
3709 /* NV reduction from T. Polzin's "Algorithms for the Steiner problem in networks" */
3710 SCIP_RETCODE reduce_nv(
3711  SCIP* scip, /**< SCIP data structure */
3712  GRAPH* g, /**< graph data structure */
3713  PATH* vnoi, /**< Voronoi data structure */
3714  double* fixed, /**< offset pointer */
3715  int* edgearrint, /**< edge int array for internal computations */
3716  int* heap, /**< heap array */
3717  int* state, /**< array for internal computations */
3718  int* vbase, /**< array for internal computations */
3719  int* solnode, /**< node array to mark whether an node is part of a given solution (CONNECT),
3720  or NULL */
3721  int* nelims /**< pointer to store number of eliminations */
3722  )
3723 {
3724  SCIP_Real* distance;
3725  SCIP_Real min1;
3726  SCIP_Real min2;
3727  SCIP_Real pi;
3728  SCIP_Real pt;
3729  SCIP_Real ttdist;
3730  int* term;
3731  int* minedge1;
3732  int* distnode;
3733  int i;
3734  int l;
3735  int k;
3736  int e;
3737  int t;
3738  int old;
3739  int edge1;
3740  int nnodes;
3741  int nterms;
3742  int mingrad;
3743  int termcount;
3744  SCIP_Bool pc;
3745  assert(g != NULL);
3746  assert(vnoi != NULL);
3747  assert(heap != NULL);
3748  assert(state != NULL);
3749  assert(vbase != NULL);
3750 
3751  t = 0;
3752  termcount = 0;
3753  *nelims = 0;
3754  pi = 0;
3755  pt = 0;
3756 
3757  nnodes = g->knots;
3758  nterms = g->terms;
3759  pc = (g->stp_type == STP_PCSPG) || (g->stp_type == STP_RPCSPG);
3760  SCIP_CALL( SCIPallocBufferArray(scip, &term, nterms) );
3761  SCIP_CALL( SCIPallocBufferArray(scip, &minedge1, nterms) );
3762  SCIP_CALL( SCIPallocBufferArray(scip, &distance, nnodes) );
3763 
3764  /* minimal grad of a vertex to be scrutinized */
3765  if( pc )
3766  {
3767  if( g->stp_type == STP_RPCSPG )
3768  mingrad = 3;
3769  else
3770  mingrad = 4;
3771 
3772  SCIP_CALL( SCIPallocBufferArray(scip, &distnode, nnodes) );
3773  }
3774  else
3775  {
3776  mingrad = 2;
3777  distnode = NULL;
3778  for( i = 0; i < nnodes; i++ )
3779  g->mark[i] = (g->grad[i] > 0);
3780  }
3781 
3782  for( i = 0; i < nnodes; i++ )
3783  {
3784  if( Is_term(g->term[i]) && g->mark[i] && g->grad[i] > 0 )
3785  {
3786  /* compute shortest incident edge */
3787  edge1 = UNKNOWN;
3788  if( g->grad[i] >= 1 )
3789  {
3790  min1 = FARAWAY;
3791 
3792  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
3793  {
3794  if( !g->mark[g->head[e]] )
3795  continue;
3796  if( SCIPisLE(scip, g->cost[e], min1) )
3797  {
3798  edge1 = e;
3799  min1 = g->cost[e];
3800  }
3801  }
3802  }
3803  minedge1[termcount] = edge1;
3804  term[termcount++] = i;
3805  }
3806  }
3807 
3808  /* compute Voronoi regions and distances */
3809  SCIP_CALL( graph_voronoiWithDist(scip, g, g->cost, distance, edgearrint, vbase, minedge1, heap, state, distnode, vnoi) );
3810 
3811  for( l = 0; l < termcount; l++ )
3812  {
3813  /* get l'th terminal */
3814  i = term[l];
3815 
3816  if( g->grad[i] < mingrad )
3817  continue;
3818 
3819  assert(minedge1[l] != UNKNOWN);
3820  /* get shortest two edges */
3821  edge1 = UNKNOWN;
3822  min2 = FARAWAY;
3823  min1 = FARAWAY;
3824  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
3825  {
3826  if( !g->mark[g->head[e]] )
3827  continue;
3828  if( SCIPisLE(scip, g->cost[e], min1) )
3829  {
3830  edge1 = e;
3831  min2 = min1;
3832  min1 = g->cost[e];
3833  }
3834  else if( SCIPisLE(scip, g->cost[e], min2) )
3835  {
3836  min2 = g->cost[e];
3837  }
3838  }
3839 
3840  assert(edge1 != UNKNOWN);
3841  assert(i == g->tail[edge1]);
3842  k = g->head[edge1];
3843 
3844  /* covered in degree test */
3845  if( Is_term(g->term[k]) )
3846  continue;
3847 
3848  if( vbase[k] != i )
3849  {
3850  if( pc )
3851  t = vbase[k];
3852  ttdist = g->cost[edge1] + vnoi[k].dist;
3853  }
3854  else
3855  {
3856  if( distnode != NULL )
3857  t = distnode[i];
3858  ttdist = distance[i];
3859  }
3860  if( pc )
3861  {
3862  if( i != g->source )
3863  pi = g->prize[i];
3864  else
3865  pi = FARAWAY;
3866 
3867  if( t != g->source )
3868  pt = g->prize[t];
3869  else
3870  pt = FARAWAY;
3871  }
3872 
3873  if( SCIPisGE(scip, min2, ttdist)
3874  && (!pc || (SCIPisLE(scip, g->cost[edge1], pi) && SCIPisLE(scip, ttdist, pt))) )
3875  {
3876  old = g->grad[i] + g->grad[k] - 1;
3877  *fixed += g->cost[edge1];
3878 
3879  if( pc )
3880  {
3881  SCIP_CALL( graph_pc_contractEdge(scip, g, solnode, i, k, i) );
3882  }
3883  else
3884  {
3885  SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->fixedges), g->ancestors[edge1], NULL) );
3886  SCIP_CALL( graph_knot_contract(scip, g, solnode, i, k) );
3887  }
3888  assert(old - g->grad[i] - g->grad[k] > 0);
3889  (*nelims) += old - g->grad[i] - g->grad[k];
3890  }
3891  }
3892 
3893  SCIPfreeBufferArrayNull(scip, &distnode);
3894  SCIPfreeBufferArray(scip, &distance);
3895  SCIPfreeBufferArray(scip, &minedge1);
3896  SCIPfreeBufferArray(scip, &term);
3897 
3898  return SCIP_OKAY;
3899 }
3900 
3901 
3902 /* advanced NV reduction */
3903 SCIP_RETCODE reduce_nvAdv(
3904  SCIP* scip, /**< SCIP data structure */
3905  GRAPH* g, /**< graph data structure */
3906  PATH* vnoi, /**< Voronoi data structure */
3907  SCIP_Real* distance, /**< nodes-sized distance array */
3908  double* fixed, /**< offset pointer */
3909  int* edgearrint, /**< edges-sized array */
3910  int* heap, /**< heap array */
3911  int* state, /**< shortest path array */
3912  int* vbase, /**< Voronoi base array */
3913  int* neighb, /**< nodes-sized neighborhood array */
3914  int* distnode, /**< nodes-sized distance array */
3915  int* solnode, /**< node array to mark whether an node is part of a given solution (CONNECT),
3916  or NULL */
3917  int* nelims /**< pointer to store number of eliminations */
3918  )
3919 {
3920  SCIP_Real min1;
3921  SCIP_Real min2;
3922  SCIP_Real min3;
3923  SCIP_Real pi;
3924  SCIP_Real pt;
3925  SCIP_Real ttdist;
3926  int* term;
3927  int* minedge1;
3928  int i;
3929  int l;
3930  int k;
3931  int e;
3932  int t;
3933  int edge1;
3934  int edge2;
3935  int nnodes;
3936  int nterms;
3937  int mingrad;
3938  int termcount;
3939  SCIP_Bool pc;
3940  SCIP_Bool contract;
3941 
3942  assert(g != NULL);
3943  assert(neighb != NULL);
3944  assert(vnoi != NULL);
3945  assert(heap != NULL);
3946  assert(state != NULL);
3947  assert(vbase != NULL);
3948 
3949  t = 0;
3950  pi = 0;
3951  pt = 0;
3952  pc = (g->stp_type == STP_PCSPG) || (g->stp_type == STP_RPCSPG);
3953  nnodes = g->knots;
3954  nterms = g->terms;
3955  *nelims = 0;
3956  termcount = 0;
3957 
3958  if( nterms <= 1 )
3959  return SCIP_OKAY;
3960 
3961  /* allocate memory */
3962  SCIP_CALL( SCIPallocBufferArray(scip, &term, nterms) );
3963  SCIP_CALL( SCIPallocBufferArray(scip, &minedge1, nterms) );
3964 
3965  /* set minimal grad of a vertex to be scrutinized */
3966  if( pc )
3967  {
3968  if( g->stp_type == STP_RPCSPG )
3969  mingrad = 3;
3970  else
3971  mingrad = 4;
3972 
3973  assert(distnode != NULL);
3974  }
3975  else
3976  {
3977  mingrad = 2;
3978  for( i = 0; i < nnodes; i++ )
3979  g->mark[i] = (g->grad[i] > 0);
3980  }
3981 
3982  /* compute shortest incident edge to each terminal */
3983  for( i = 0; i < nnodes; i++ )
3984  {
3985  neighb[i] = FALSE;
3986  if( Is_term(g->term[i]) && g->mark[i] && g->grad[i] > 0 )
3987  {
3988  /* compute shortest incident edge */
3989  edge1 = UNKNOWN;
3990  if( g->grad[i] >= 1 )
3991  {
3992  min1 = FARAWAY;
3993 
3994  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
3995  {
3996  if( g->mark[g->head[e]] && SCIPisLE(scip, g->cost[e], min1) )
3997  {
3998  edge1 = e;
3999  min1 = g->cost[e];
4000  }
4001  }
4002  }
4003 
4004  minedge1[termcount] = edge1;
4005  term[termcount++] = i;
4006  }
4007  }
4008 
4009  /* compute Voronoi regions and distances */
4010  SCIP_CALL( graph_voronoiWithDist(scip, g, g->cost, distance, edgearrint, vbase, minedge1, heap, state, distnode, vnoi) );
4011 
4012  /* main loop: try to contract (shortest) edges into terminals */
4013  for( l = 0; l < termcount; l++ )
4014  {
4015  /* get l'th terminal */
4016  i = term[l];
4017 
4018  if( g->grad[i] < mingrad )
4019  continue;
4020 
4021  assert(minedge1[l] != UNKNOWN);
4022 
4023  /* get shortest two edges */
4024 
4025  min3 = FARAWAY;
4026  min2 = FARAWAY;
4027  min1 = FARAWAY;
4028  edge1 = UNKNOWN;
4029  edge2 = UNKNOWN;
4030 
4031  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
4032  {
4033  if( !g->mark[g->head[e]] )
4034  continue;
4035  neighb[g->head[e]] = TRUE;
4036 
4037  if( SCIPisLE(scip, g->cost[e], min1) )
4038  {
4039  edge2 = edge1;
4040  edge1 = e;
4041 
4042  min3 = min2;
4043  min2 = min1;
4044  min1 = g->cost[e];
4045  }
4046  else if( SCIPisLE(scip, g->cost[e], min2) )
4047  {
4048  edge2 = e;
4049 
4050  min3 = min2;
4051  min2 = g->cost[e];
4052  }
4053  else if( SCIPisLE(scip, g->cost[e], min3) )
4054  {
4055  min3 = g->cost[e];
4056  }
4057  }
4058 
4059  assert(edge1 != UNKNOWN);
4060  assert(i == g->tail[edge1]);
4061 
4062  k = g->head[edge1];
4063 
4064  /* covered in degree test */
4065  if( Is_term(g->term[k]) )
4066  continue;
4067 
4068  if( vbase[k] != i )
4069  {
4070  if( pc )
4071  t = vbase[k];
4072  ttdist = g->cost[edge1] + vnoi[k].dist;
4073  }
4074  else
4075  {
4076  if( pc )
4077  {
4078  assert(distnode != NULL);
4079  t = distnode[i];
4080  }
4081  ttdist = distance[i];
4082  }
4083  if( pc )
4084  {
4085  if( i != g->source )
4086  pi = g->prize[i];
4087  else
4088  pi = FARAWAY;
4089 
4090  if( t == UNKNOWN )
4091  pt = -1.0;
4092  else if( t != g->source )
4093  pt = g->prize[t];
4094  else
4095  pt = FARAWAY;
4096  }
4097 
4098  contract = FALSE;
4099 
4100  if( SCIPisGE(scip, min2, ttdist) )
4101  {
4102  contract = TRUE;
4103  }
4104  else if( edge2 != UNKNOWN && !Is_term(g->term[g->head[edge2]]) && SCIPisGE(scip, min3, ttdist) )
4105  {
4106  t = g->head[edge2];
4107  for( e = g->outbeg[t]; e != EAT_LAST; e = g->oeat[e] )
4108  if( e != flipedge(edge2) && SCIPisLT(scip, g->cost[e], ttdist) )/*&& !neighb[g->head[e]] ) */
4109  break;
4110 
4111  if( e == EAT_LAST )
4112  contract = TRUE;
4113  }
4114 
4115  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
4116  neighb[g->head[e]] = FALSE;
4117 
4118  if( contract && (!pc || (SCIPisLE(scip, g->cost[edge1], pi) && SCIPisLE(scip, ttdist, pt))) )
4119  {
4120  (*nelims)++;
4121  *fixed += g->cost[edge1];
4122 
4123  if( pc )
4124  {
4125  SCIP_CALL( graph_pc_contractEdge(scip, g, solnode, i, k, i) );
4126  }
4127  else
4128  {
4129  SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->fixedges), g->ancestors[edge1], NULL) );
4130  SCIP_CALL( graph_knot_contract(scip, g, solnode, i, k) );
4131  }
4132  }
4133  }
4134 
4135  SCIPfreeBufferArray(scip, &minedge1);
4136  SCIPfreeBufferArray(scip, &term);
4137 
4138  return SCIP_OKAY;
4139 }
4140 
4141 
4142 /* longest edge reduction test from T. Polzin's "Algorithms for the Steiner problem in networks" (Lemma 20) */
4143 SCIP_RETCODE reduce_ledge(
4144  SCIP* scip,
4145  GRAPH* g,
4146  PATH* vnoi,
4147  int* heap,
4148  int* state,
4149  int* vbase,
4150  int* nelims,
4151  int* edgestate
4152 )
4153 {
4154  GRAPH* netgraph;
4155  PATH* mst;
4156  SCIP_Real cost;
4157  SCIP_Real maxcost;
4158  int v1;
4159  int v2;
4160  int k;
4161  int e;
4162  int ne;
4163  int nedges;
4164  int nnodes;
4165  int nterms;
4166  int maxnedges;
4167  int netnnodes;
4168  int* nodesid;
4169  int* edgeorg;
4170  SCIP_Bool checkstate = (edgestate != NULL);
4171  STP_Bool* blocked;
4172 
4173  assert(g != NULL);
4174  assert(vnoi != NULL);
4175  assert(heap != NULL);
4176  assert(state != NULL);
4177  assert(vbase != NULL);
4178 
4179  *nelims = 0;
4180  nedges = g->edges;
4181  nnodes = g->knots;
4182  assert(graph_valid(g));
4183 
4184  nterms = 0;
4185  for( k = 0; k < nnodes; k++ )
4186  {
4187  g->mark[k] = (g->grad[k] > 0);
4188  if( Is_term(g->term[k]) && g->mark[k] )
4189  nterms++;
4190  }
4191 
4192  if( nterms <= 1 )
4193  return SCIP_OKAY;
4194 
4195  SCIP_CALL( SCIPallocBufferArray(scip, &blocked, nedges / 2) );
4196  SCIP_CALL( SCIPallocBufferArray(scip, &edgeorg, nedges / 2) );
4197 
4198  graph_voronoiTerms(scip, g, g->cost, vnoi, vbase, heap, state);
4199 
4200  if( nedges >= (nterms - 1) * nterms )
4201  maxnedges = (nterms - 1) * nterms;
4202  else
4203  maxnedges = nedges;
4204 
4205  SCIP_CALL( SCIPallocBufferArray(scip, &nodesid, nnodes) );
4206 
4207  /* initialize the new graph */
4208  SCIP_CALL( graph_init(scip, &netgraph, nterms, maxnedges, 1) );
4209 
4210  e = 0;
4211  for( k = 0; k < nnodes; k++ )
4212  {
4213  if( Is_term(g->term[k]) && g->grad[k] > 0 )
4214  {
4215  if( e == 0 )
4216  graph_knot_add(netgraph, 0);
4217  else
4218  graph_knot_add(netgraph, -1);
4219  nodesid[k] = e++;
4220  }
4221  else
4222  {
4223  nodesid[k] = UNKNOWN;
4224  }
4225  }
4226 
4227  netnnodes = netgraph->knots;
4228  assert(netnnodes == e);
4229  assert(netnnodes == nterms);
4230 
4231  for( e = 0; e < nedges / 2; e++ )
4232  {
4233  blocked[e] = FALSE;
4234  edgeorg[e] = UNKNOWN;
4235  }
4236 
4237  for( k = 0; k < nnodes; k++ )
4238  {
4239  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4240  {
4241  v1 = vbase[k];
4242  assert(k == g->tail[e]);
4243 
4244  if( v1 != vbase[g->head[e]] )
4245  {
4246  v2 = vbase[g->head[e]];
4247  assert(Is_term(g->term[v1]));
4248  assert(Is_term(g->term[v2]));
4249  assert(nodesid[v1] >= 0);
4250  assert(nodesid[v2] >= 0);
4251 
4252  for( ne = netgraph->outbeg[nodesid[v1]]; ne != EAT_LAST; ne = netgraph->oeat[ne] )
4253  if( netgraph->head[ne] == nodesid[v2] )
4254  break;
4255 
4256  cost = g->cost[e] + vnoi[g->head[e]].dist + vnoi[g->tail[e]].dist;
4257  /* edge exists? */
4258  if( ne != EAT_LAST )
4259  {
4260  assert(ne >= 0);
4261  assert(netgraph->head[ne] == nodesid[v2]);
4262  assert(netgraph->tail[ne] == nodesid[v1]);
4263  if( SCIPisGT(scip, netgraph->cost[ne], cost) )
4264  {
4265  netgraph->cost[ne] = cost;
4266  netgraph->cost[Edge_anti(ne)] = cost;
4267  edgeorg[ne / 2] = e;
4268  assert(ne <= maxnedges);
4269  }
4270  }
4271  else
4272  {
4273  edgeorg[netgraph->edges / 2] = e;
4274  graph_edge_add(scip, netgraph, nodesid[v1], nodesid[v2], cost, cost);
4275  assert(netgraph->edges <= maxnedges);
4276  }
4277  }
4278  }
4279  }
4280  netgraph->source = 0;
4281 
4282  assert(graph_valid(netgraph));
4283 
4284  for( k = 0; k < netnnodes; k++ )
4285  netgraph->mark[k] = TRUE;
4286 
4287  /* compute a MST on netgraph */
4288  SCIP_CALL( SCIPallocBufferArray(scip, &mst, netnnodes) );
4289  SCIP_CALL( graph_path_init(scip, netgraph) );
4290  graph_path_exec(scip, netgraph, MST_MODE, 0, netgraph->cost, mst);
4291 
4292  maxcost = -1;
4293  assert(mst[0].edge == -1);
4294 
4295  for( k = 1; k < netnnodes; k++ )
4296  {
4297  assert(netgraph->path_state[k] == CONNECT);
4298  e = mst[k].edge;
4299  assert(e >= 0);
4300  cost = netgraph->cost[e];
4301  if( SCIPisGT(scip, cost, maxcost) )
4302  maxcost = cost;
4303 
4304  ne = edgeorg[e / 2];
4305  blocked[ne / 2] = TRUE;
4306  for( v1 = g->head[ne]; v1 != vbase[v1]; v1 = g->tail[vnoi[v1].edge] )
4307  blocked[vnoi[v1].edge / 2] = TRUE;
4308 
4309  for( v1 = g->tail[ne]; v1 != vbase[v1]; v1 = g->tail[vnoi[v1].edge] )
4310  blocked[vnoi[v1].edge / 2] = TRUE;
4311  assert(e != EAT_LAST);
4312  }
4313 
4314  for( k = 0; k < nnodes; k++ )
4315  {
4316  e = g->outbeg[k];
4317  while( e != EAT_LAST )
4318  {
4319  assert(e >= 0);
4320  if( checkstate && (edgestate[e] == EDGE_BLOCKED) )
4321  {
4322  e = g->oeat[e];
4323  continue;
4324  }
4325 
4326  if( SCIPisGE(scip, g->cost[e], maxcost) && !blocked[e / 2] )
4327  {
4328  (*nelims)++;
4329  v1 = g->oeat[e];
4330  graph_edge_del(scip, g, e, TRUE);
4331  e = v1;
4332  }
4333  else
4334  {
4335  e = g->oeat[e];
4336  }
4337  }
4338  }
4339 
4340  /* graph might have become disconnected */
4341  if( *nelims > 0 )
4342  {
4343  SCIP_CALL( level0(scip, g) );
4344  }
4345 
4346  /* free netgraph and MST data structure */
4347  graph_path_exit(scip, netgraph);
4348  graph_free(scip, &netgraph, TRUE);
4349  SCIPfreeBufferArray(scip, &mst);
4350  SCIPfreeBufferArray(scip, &nodesid);
4351  SCIPfreeBufferArray(scip, &edgeorg);
4352  SCIPfreeBufferArray(scip, &blocked);
4353 
4354  assert(graph_valid(g));
4355  return SCIP_OKAY;
4356 }
4357 
4358 
4359 #if 0
4360 /** domination vertex reduction for the SPG */
4361 void reduce_alt_dv(
4362  SCIP* scip, /**< SCIP data structure */
4363  GRAPH* g, /**< graph data structure */
4364  STP_Bool* marked, /**< nodes array */
4365  int* count /**< pointer to number of reductions */
4366  )
4367 {
4368  const int nnodes = g->knots;
4369  int nreds = 0;
4370 
4371  assert(scip != NULL);
4372  assert(g != NULL);
4373  assert(count != NULL);
4374  assert(marked != NULL);
4375  assert(g->stp_type == STP_SPG);
4376 
4377  BMSclearMemoryArray(marked, nnodes);
4378 
4379  /* main loop */
4380  for( int k = 0; k < nnodes; k++ )
4381  {
4382  const SCIP_Real maxgrad = g->grad[k];
4383 
4384  if( maxgrad < 3 )
4385  continue;
4386 
4387  assert(g->mark[k]);
4388 
4389  /* mark adjacent vertices and k*/
4390  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4391  marked[g->head[e]] = TRUE;
4392 
4393  marked[k] = TRUE;
4394 
4395  /* check all other nodes */
4396  for( int i = 0; i < nnodes; i++ )
4397  {
4398  /* valid candidate? */
4399  if( !Is_term(g->term[i]) && g->grad[i] <= maxgrad && g->mark[i] && k != i )
4400  {
4401  SCIP_Real min;
4402  int e2;
4403  for( e2 = g->outbeg[i]; e2 != EAT_LAST; e2 = g->oeat[e2] )
4404  if( !marked[g->head[e2]] )
4405  break;
4406 
4407  /* neighbors of j subset of those of k? */
4408  if( e2 == EAT_LAST )
4409  {
4410 #if 0
4411  int maxe = g->outbeg[i];
4412  for( e2 = g->outbeg[i]; e2 != EAT_LAST; e2 = g->oeat[e2] )
4413  if( g->cost[e] > g->cost[maxe] )
4414  maxe = e;
4415 
4416  min = 0.0;
4417 
4418 
4419 
4420 
4421  (*count) += g->grad[i];
4422  while( g->outbeg[i] != EAT_LAST )
4423  graph_edge_del(scip, g, g->outbeg[j] TRUE);
4424 
4425  g->mark[i] = FALSE;
4426  marked[i] = FALSE;
4427 #endif
4428  }
4429  }
4430  }
4431 
4432  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4433  marked[g->head[e]] = FALSE;
4434 
4435  marked[k] = FALSE;
4436 
4437  for( int j = 0; j < nnodes; j++ )
4438  assert(marked[j] == FALSE);
4439  }
4440 
4441  *count += nreds;
4442 }
4443 
4444 #endif
4445 
4446 /** adjacent neighbourhood reduction for the MWCSP */
4447 void reduce_ans(
4448  SCIP* scip, /**< SCIP data structure */
4449  GRAPH* g, /**< graph data structure */
4450  int* marked, /**< nodes array */
4451  int* count /**< pointer to number of reductions */
4452  )
4453 {
4454  const int nnodes = g->knots;
4455 
4456  assert(scip != NULL);
4457  assert(g != NULL);
4458  assert(count != NULL);
4459  assert(marked != NULL);
4460  assert(g->stp_type == STP_MWCSP);
4461  assert(graph_valid(g));
4462 
4463  *count = 0;
4464 
4465  /* unmark all nodes */
4466  for( int k = 0; k < nnodes; k++ )
4467  marked[k] = FALSE;
4468 
4469  /* check neighborhood of all nodes */
4470  for( int k = 0; k < nnodes; k++ )
4471  {
4472  SCIP_Real min;
4473  int e;
4474  int neighborcount = 0;
4475 
4476  if( !g->mark[k] )
4477  continue;
4478 
4479  /* mark adjacent vertices and k*/
4480  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4481  marked[g->head[e]] = TRUE;
4482 
4483  marked[k] = TRUE;
4484 
4485  if( SCIPisLT(scip, g->prize[k], 0.0) )
4486  min = g->prize[k];
4487  else
4488  min = 0.0;
4489 
4490  /* check all neighbors of k */
4491  e = g->outbeg[k];
4492  while( e != EAT_LAST && neighborcount++ < STP_RED_ANSMAXNEIGHBORS )
4493  {
4494  const int j = g->head[e];
4495  e = g->oeat[e];
4496 
4497  /* valid candidate? */
4498  if( g->grad[j] <= g->grad[k] && !Is_term(g->term[j]) && g->mark[j] )
4499  ansProcessCandidate(scip, g, marked, count, min, j);
4500  }
4501 
4502  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4503  marked[g->head[e]] = FALSE;
4504 
4505  marked[k] = FALSE;
4506 
4507  for( int j = 0; j < nnodes; j++ )
4508  assert(marked[j] == FALSE);
4509  }
4510 
4511  assert(graph_valid(g));
4512 }
4513 
4514 /** advanced adjacent neighbourhood reduction for the MWCSP */
4515 void reduce_ansAdv(
4516  SCIP* scip, /**< SCIP data structure */
4517  GRAPH* g, /**< graph data structure */
4518  int* marked, /**< nodes array */
4519  int* count, /**< pointer to number of performed reductions */
4520  SCIP_Bool extneigbhood /**< use extended neighbour hood */
4521  )
4522 {
4523  int candidates[MAX(STP_RED_ANSMAXCANDS, STP_RED_ANSEXMAXCANDS)];
4524  int neighbarr[STP_RED_CNSNN];
4525  const int nnodes = g->knots;
4526 
4527  assert(scip != NULL);
4528  assert(g != NULL);
4529  assert(count != NULL);
4530  assert(marked != NULL);
4531  assert(g->stp_type == STP_MWCSP);
4532 
4533  *count = 0;
4534 
4535  /* unmark all nodes */
4536  for( int k = 0; k < nnodes; k++ )
4537  marked[k] = FALSE;
4538 
4539  /* check neighborhood of all nodes */
4540  for( int k = 0; k < nnodes; k++ )
4541  {
4542  SCIP_Real min;
4543  int nn;
4544  int ncands;
4545  int maxgrad;
4546 
4547  if( (!(g->mark[k])) || (g->grad[k] < 2) || Is_term(g->term[k]) )
4548  continue;
4549 
4550  nn = 0;
4551 
4552  /* mark adjacent vertices and k */
4553  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4554  {
4555  const int j = g->head[e];
4556  marked[j] = TRUE;
4557  if( SCIPisGT(scip, g->prize[j], 0.0) && nn < STP_RED_CNSNN )
4558  neighbarr[nn++] = j;
4559  }
4560 
4561  marked[k] = TRUE;
4562  maxgrad = g->grad[k];
4563  ncands = 0;
4564 
4565  for( int l = 0; l < nn; l++ )
4566  {
4567  for( int e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e] )
4568  {
4569  marked[g->head[e]] = TRUE;
4570  }
4571  maxgrad += g->grad[neighbarr[l]];
4572  }
4573 
4574  assert(SCIPisLE(scip, g->prize[k], 0.0));
4575 
4576  min = g->prize[k];
4577 
4578  if( extneigbhood )
4579  {
4580  assert(0 && "implement me");
4581  }
4582  else
4583  {
4584  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4585  {
4586  const int j = g->head[e];
4587  if( g->grad[j] <= maxgrad && g->mark[j] && !Is_term(g->term[j]) )
4588  {
4589  candidates[ncands++] = j;
4590  if( ncands >= STP_RED_ANSMAXCANDS )
4591  {
4592  SCIPdebugMessage("REACHED ANS LIMIT %d \n", ncands);
4593  break;
4594  }
4595  }
4596  }
4597  }
4598 
4599  /* check all neighbors of k */
4600  for( int l = 0; l < ncands; l++ )
4601  ansProcessCandidate(scip, g, marked, count, min, candidates[l]);
4602 
4603  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4604  marked[g->head[e]] = FALSE;
4605 
4606  for( int l = 0; l < nn; l++ )
4607  for( int e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e] )
4608  marked[g->head[e]] = FALSE;
4609 
4610  marked[k] = FALSE;
4611 
4612  for( int k2 = 0; k2 < nnodes; k2++ )
4613  assert(marked[k2] == FALSE);
4614  }
4615 
4616  assert(graph_valid(g));
4617 }
4618 
4619 
4620 /** alternative advanced adjacent neighbourhood reduction for the MWCSP */
4621 void reduce_ansAdv2(
4622  SCIP* scip, /**< SCIP data structure */
4623  GRAPH* g, /**< graph data structure */
4624  int* marked, /**< nodes array */
4625  int* count /**< pointer to number of reductions */
4626  )
4627 {
4628  int neighbarr[STP_RED_CNSNN + 1];
4629  SCIP_Real min;
4630  const int nnodes = g->knots;
4631 
4632  assert(scip != NULL);
4633  assert(g != NULL);
4634  assert(count != NULL);
4635  assert(marked != NULL);
4636  assert(g->stp_type == STP_MWCSP);
4637 
4638  *count = 0;
4639 
4640  /* unmark all nodes */
4641  for( int k = 0; k < nnodes; k++ )
4642  marked[k] = FALSE;
4643 
4644  /* check neighbourhood of all nodes */
4645  for( int k = 0; k < nnodes; k++ )
4646  {
4647  SCIP_Real maxprize;
4648  int neighbor0;
4649 
4650  if( (!(g->mark[k])) || (g->grad[k] < 2) || Is_term(g->term[k]) )
4651  continue;
4652 
4653  maxprize = g->prize[k];
4654  neighbor0 = -1;
4655 
4656  for( int run = 0; run < 2; run++ )
4657  {
4658  int e;
4659  int nn = 0;
4660  int k2 = UNKNOWN;
4661  int maxgrad = g->grad[k];
4662 
4663  /* mark adjacent vertices and k */
4664  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4665  {
4666  const int j = g->head[e];
4667  marked[j] = TRUE;
4668  if( SCIPisGE(scip, g->prize[j], 0.0) && nn < STP_RED_CNSNN - 1 )
4669  {
4670  neighbarr[nn++] = j;
4671  }
4672  else if( SCIPisGT(scip, g->prize[j], maxprize) && g->mark[j] && j != neighbor0 )
4673  {
4674  maxprize = g->prize[j];
4675  k2 = j;
4676  }
4677  }
4678 
4679  marked[k] = TRUE;
4680 
4681  if( run == 0 && k2 != UNKNOWN )
4682  neighbarr[nn++] = k2;
4683 
4684  for( int l = 0; l < nn; l++ )
4685  {
4686  for( e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e] )
4687  {
4688  const int j = g->head[e];
4689  if( run == 1 && g->mark[j] && !Is_term(g->term[j]) && SCIPisGT(scip, g->prize[j], maxprize) && j != neighbor0 )
4690  {
4691  maxprize = g->prize[j];
4692  k2 = j;
4693  }
4694  marked[j] = TRUE;
4695  }
4696  maxgrad += g->grad[neighbarr[l]];
4697  }
4698 
4699  if( run == 1 && k2 != UNKNOWN )
4700  {
4701  maxgrad += g->grad[k2];
4702  neighbarr[nn++] = k2;
4703 
4704  for( e = g->outbeg[k2]; e != EAT_LAST; e = g->oeat[e] )
4705  marked[g->head[e]] = TRUE;
4706  }
4707 
4708  assert(SCIPisLE(scip, g->prize[k], 0.0));
4709 
4710  min = g->prize[k];
4711  if( k2 != UNKNOWN )
4712  {
4713  neighbor0 = k2;
4714  min += g->prize[k2];
4715  }
4716 
4717  /* check all neighbours of k */
4718  e = g->outbeg[k];
4719  while( e != EAT_LAST )
4720  {
4721  const int j = g->head[e];
4722  e = g->oeat[e];
4723 
4724  /* valid candidate? */
4725  if( g->grad[j] <= maxgrad && g->mark[j] && SCIPisLE(scip, g->prize[j], min) && k2 != j )
4726  ansProcessCandidate(scip, g, marked, count, min, j);
4727  }
4728 
4729  for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4730  marked[g->head[e]] = FALSE;
4731 
4732  for( int l = 0; l < nn; l++ )
4733  for( e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e] )
4734  marked[g->head[e]] = FALSE;
4735 
4736  marked[k] = FALSE;
4737 
4738  for( k2 = 0; k2 < nnodes; k2++ )
4739  assert(marked[k2] == FALSE);
4740  }
4741  }
4742 
4743  assert(graph_valid(g));
4744 }
4745 
4746 
4747 /** advanced connected neighborhood subset reduction test for the MWCSP */
4748 SCIP_RETCODE reduce_cnsAdv(
4749  SCIP* scip, /**< SCIP data structure */
4750  GRAPH* g, /**< graph data structure */
4751  int* marked, /**< nodes array for internal use */
4752  int* count /**< pointer to number of reductions */
4753  )
4754 {
4755  SCIP_Real min;
4756  SCIP_Real kprize;
4757  int neighbarr[STP_RED_CNSNN + 1];
4758  int neighbarr2[STP_RED_CNSNN + 1];
4759  int k;
4760  int j;
4761  int e;
4762  int k2;
4763  int e2;
4764  int l;
4765  int nn;
4766  int nn2;
4767  int k2grad;
4768  int nnodes;
4769  int maxgrad;
4770 
4771  assert(scip != NULL);
4772  assert(g != NULL);
4773  assert(count != NULL);
4774  assert(marked != NULL);
4775  assert(g->stp_type == STP_MWCSP);
4776 
4777  k2grad = 0;
4778  *count = 0;
4779  nnodes = g->knots;
4780 
4781  /* unmark all nodes */
4782  for( k = 0; k < nnodes; k++ )
4783  marked[k] = VERTEX_OTHER;
4784 
4785  /* first run: consider node plus adjacent terminals */
4786 
4787  /* check neighborhood of all nodes */
4788  for( k = 0; k < nnodes; k++ )
4789  {
4790  if( !(g->mark[k]) || (g->grad[k] < 2) )
4791  continue;
4792 
4793  nn = 0;
4794  k2 = UNKNOWN;
4795  nn2 = 0;
4796  kprize = g->prize[k];
4797  maxgrad = g->grad[k];
4798 
4799  /* mark adjacent vertices and k */
4800  for (e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e])
4801  {
4802  j = g->head[e];
4803 
4804  if( !g->mark[j] )
4805  continue;
4806 
4807  if( SCIPisGE(scip, g->prize[j], 0.0) && nn < STP_RED_CNSNN - 1 )
4808  {
4809  neighbarr[nn++] = j;
4810  marked[j] = VERTEX_CONNECT;
4811  }
4812  else
4813  {
4814  marked[j] = VERTEX_NEIGHBOR;
4815  }
4816  }
4817 
4818  marked[k] = VERTEX_CONNECT;
4819 
4820  /* traverse all connected non-negative nodes and mark their neighbors */
4821  for (l = 0; l < nn; l++)
4822  {
4823  for (e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e])
4824  {
4825  j = g->head[e];
4826  if( !g->mark[j] )
4827  continue;
4828 
4829  if( marked[j] == VERTEX_OTHER )
4830  marked[j] = VERTEX_NEIGHBOR;
4831  }
4832  maxgrad += g->grad[neighbarr[l]] - 1;
4833  }
4834 
4835  if( Is_term(g->term[k]) )
4836  min = 0.0;
4837  else
4838  min = g->prize[k];
4839 
4840  /* traverse all vertices (main loop) */
4841  for (j = 0; j < nnodes; j++)
4842  {
4843  /* vertex part of the current connected subset? Or terminal? Or belonging to the extension of the graph? */
4844  if( marked[j] != VERTEX_CONNECT && g->mark[j] && !Is_term(g->term[j])
4845  &&
4846  /* valid candidate? */
4847  g->grad[j] <= maxgrad && SCIPisLE(scip, g->prize[j], min) )
4848  {
4849  for (e2 = g->outbeg[j]; e2 != EAT_LAST; e2 = g->oeat[e2])
4850  if( marked[g->head[e2]] == VERTEX_OTHER )
4851  break;
4852 
4853  /* neighbors of j subset of those of k? */
4854  if( e2 == EAT_LAST )
4855  {
4856  /* yes, delete vertex */
4857  while (g->outbeg[j] != EAT_LAST)
4858  {
4859  e2 = g->outbeg[j];
4860  (*count)++;
4861  graph_edge_del(scip, g, e2, TRUE);
4862  }
4863  g->mark[j] = FALSE;
4864  marked[j] = VERTEX_OTHER;
4865  }
4866  }
4867  }
4868 
4869  for (e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e])
4870  marked[g->head[e]] = VERTEX_OTHER;
4871 
4872  for (l = 0; l < nn; l++)
4873  for (e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e])
4874  marked[g->head[e]] = VERTEX_OTHER;
4875 
4876  marked[k] = VERTEX_OTHER;
4877 
4878 #ifdef DEBUG
4879  for( l = 0; l < nnodes; l++ )
4880  assert(marked[l] == VERTEX_OTHER);
4881 #endif
4882 
4883  }
4884  /* second run: consider the same plus an additional (non-positive) vertex */
4885 
4886  for (k = 0; k < nnodes; k++)
4887  {
4888  if( !(g->mark[k]) || g->grad[k] < 2 || Is_term(g->term[k]) )
4889  continue;
4890 
4891  nn = 0;
4892  k2 = UNKNOWN;
4893  nn2 = 0;
4894  kprize = g->prize[k];
4895  maxgrad = g->grad[k];
4896 
4897  /* mark adjacent vertices and k */
4898  for (e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e])
4899  {
4900  j = g->head[e];
4901 
4902  if( !g->mark[j] )
4903  continue;
4904 
4905  if( SCIPisGE(scip, g->prize[j], 0.0) && nn < STP_RED_CNSNN - 1 )
4906  {
4907  neighbarr[nn++] = j;
4908  marked[j] = VERTEX_CONNECT;
4909  }
4910  else if( (SCIPisGT(scip, g->prize[j], kprize) && nn2 < STP_RED_CNSNN)
4911  || (SCIPisGE(scip, g->prize[j], kprize) && j > k && nn2 < 3) )
4912  {
4913  neighbarr2[nn2++] = j;
4914  marked[j] = VERTEX_NEIGHBOR;
4915  }
4916  else
4917  {
4918  marked[j] = VERTEX_NEIGHBOR;
4919  }
4920  }
4921 
4922  marked[k] = VERTEX_CONNECT;
4923 
4924  /* traverse all connected non-negative nodes and mark their neighbors */
4925  for (l = 0; l < nn; l++)
4926  {
4927  for (e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e])
4928  {
4929  j = g->head[e];
4930  if( !g->mark[j] )
4931  continue;
4932 
4933  if( marked[j] == VERTEX_OTHER && nn2 < STP_RED_CNSNN
4934  && (SCIPisGT(scip, g->prize[j], kprize)
4935  || (SCIPisGE(scip, g->prize[j], kprize) && j > k
4936  && nn2 < 3)) )
4937  {
4938  neighbarr2[nn2++] = j;
4939  marked[j] = VERTEX_NEIGHBOR;
4940  }
4941  else if( marked[j] == VERTEX_OTHER )
4942  {
4943  marked[j] = VERTEX_NEIGHBOR;
4944  }
4945  }
4946  maxgrad += g->grad[neighbarr[l]] - 1;
4947  }
4948 
4949  if( Is_term(g->term[k]) )
4950  min = 0.0;
4951  else
4952  min = g->prize[k];
4953 
4954  for (l = 0; l < nn2; l++)
4955  {
4956  k2 = neighbarr2[l];
4957 
4958  if( !g->mark[k2] )
4959  continue;
4960 
4961  for (e = g->outbeg[k2]; e != EAT_LAST; e = g->oeat[e])
4962  if( marked[g->head[e]] == VERTEX_OTHER && g->mark[g->head[e]] )
4963  marked[g->head[e]] = VERTEX_TEMPNEIGHBOR;
4964  min += g->prize[k2];
4965  k2grad = g->grad[k2];
4966  maxgrad += k2grad - 1;
4967  assert(SCIPisLE(scip, g->prize[k2], 0.0));
4968 
4969  /* traverse all vertices (main loop) */
4970  for (j = 0; j < nnodes; j++)
4971  {
4972  /* vertex part of the current connected subset? Or terminal? Or belonging to the extension of the graph? */
4973  if( marked[j] != VERTEX_CONNECT && g->mark[j]
4974  && !Is_term(g->term[j]) &&
4975  /* valid candidate? */
4976  g->grad[j] <= maxgrad && SCIPisLE(scip, g->prize[j], min) )
4977  {
4978  for (e2 = g->outbeg[j]; e2 != EAT_LAST; e2 = g->oeat[e2])
4979  if( marked[g->head[e2]] == VERTEX_OTHER )
4980  break;
4981 
4982  /* neighbors of j subset of those of k? */
4983  if( e2 == EAT_LAST )
4984  {
4985  /* yes, delete vertex */
4986  while (g->outbeg[j] != EAT_LAST)
4987  {
4988  e2 = g->outbeg[j];
4989  (*count)++;
4990  graph_edge_del(scip, g, e2, TRUE);
4991  }
4992  g->mark[j] = FALSE;
4993  marked[j] = VERTEX_OTHER;
4994  }
4995  }
4996  }
4997  for (e = g->outbeg[k2]; e != EAT_LAST; e = g->oeat[e])
4998  if( marked[g->head[e]] == VERTEX_TEMPNEIGHBOR
4999  && g->mark[g->head[e]] )
5000  marked[g->head[e]] = VERTEX_OTHER;
5001  min -= g->prize[k2];
5002  maxgrad -= k2grad - 1;
5003 
5004  }
5005  for (e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e])
5006  marked[g->head[e]] = VERTEX_OTHER;
5007 
5008  for (l = 0; l < nn; l++)
5009  for (e = g->outbeg[neighbarr[l]]; e != EAT_LAST; e = g->oeat[e])
5010  marked[g->head[e]] = VERTEX_OTHER;
5011 
5012  marked[k] = VERTEX_OTHER;
5013 #ifdef DEBUG
5014  for( l = 0; l < nnodes; l++ )
5015  assert(marked[l] == VERTEX_OTHER);
5016 #endif
5017  }
5018 
5019  return SCIP_OKAY;
5020 }
5021 
5022 
5023 /** non-positive vertex reduction for the MWCSP */
5024 SCIP_RETCODE reduce_npv(
5025  SCIP* scip,
5026  GRAPH* g,
5027  PATH* pathtail,
5028  PATH* pathhead,
5029  int* heap,
5030  int* statetail,
5031  int* statehead,
5032  int* memlbltail,
5033  int* memlblhead,
5034  int* nelims,
5035  int limit
5036  )
5037 {
5038  GRAPH* auxg;
5039  PATH mst[5];
5040  int adjverts[5];
5041 
5042  SCIP_Real prize;
5043  SCIP_Real sdist0;
5044  SCIP_Real sdist1;
5045  SCIP_Real sdist2;
5046 
5047  const int nnodes = g->knots;;
5048 
5049  assert(g != NULL);
5050  assert(scip != NULL);
5051  assert(pathtail != NULL);
5052  assert(pathhead != NULL);
5053  assert(heap != NULL);
5054  assert(statetail != NULL);
5055  assert(statehead != NULL);
5056  assert(memlbltail != NULL);
5057  assert(memlblhead != NULL);
5058  assert(nelims != NULL);
5059 
5060  *nelims = 0;
5061 
5062  /* initialize arrays */
5063  for( int i = 0; i < nnodes; i++ )
5064  {
5065  statetail[i] = UNKNOWN;
5066  pathtail[i].dist = FARAWAY;
5067  pathtail[i].edge = UNKNOWN;
5068  statehead[i] = UNKNOWN;
5069  pathhead[i].dist = FARAWAY;
5070  pathhead[i].edge = UNKNOWN;
5071  }
5072 
5073 
5074  /* --- NPV3 test --- */
5075 
5076  /* try to eliminate non-positive vertices of degree 3 */
5077  for( int i = 0; i < nnodes; i++ )
5078  {
5079  int k;
5080 
5081  assert(g->grad[i] >= 0);
5082 
5083  /* only non-positive vertices of degree 3 */
5084  if( !g->mark[i] || g->grad[i] != 3 || Is_term(g->term[i]) )
5085  continue;
5086 
5087  k = 0;
5088  for( int e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
5089  {
5090  assert(g->head[e] != g->source);
5091  assert(k < 3);
5092 
5093  adjverts[k++] = g->head[e];
5094  }
5095 
5096  assert(k == 3);
5097 
5098  g->mark[i] = FALSE;
5099 
5100  prize = g->prize[i];
5101  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &sdist0, -prize, heap, statetail, statehead, memlbltail, memlblhead, adjverts[0], adjverts[1], limit, FALSE, TRUE) );
5102  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &sdist1, -prize, heap, statetail, statehead, memlbltail, memlblhead, adjverts[1], adjverts[2], limit, FALSE, TRUE) );
5103  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &sdist2, -prize, heap, statetail, statehead, memlbltail, memlblhead, adjverts[2], adjverts[0], limit, FALSE, TRUE) );
5104 
5105  /* can vertex be deleted? */
5106  if( (SCIPisGE(scip, -sdist0 - sdist1, prize) && SCIPisGE(scip, -sdist2, prize))
5107  || (SCIPisGE(scip, -sdist1 - sdist2, prize) && SCIPisGE(scip, -sdist0, prize))
5108  || (SCIPisGE(scip, -sdist2 - sdist0, prize) && SCIPisGE(scip, -sdist1, prize))
5109  )
5110  {
5111  SCIPdebugMessage("npv3Reduction delete: %d (prize: %f) \n", i, g->prize[i]);
5112  (*nelims) += g->grad[i];
5113  while( g->outbeg[i] != EAT_LAST )
5114  graph_edge_del(scip, g, g->outbeg[i], TRUE);
5115  }
5116  else
5117  {
5118  g->mark[i] = TRUE;
5119  }
5120  }
5121 
5122  /* --- NPV4 test --- */
5123 
5124  /* initialize mst struct and new graph for further tests */
5125  SCIP_CALL( graph_init(scip, &auxg, 5, 40, 1) );
5126 
5127  for( int k = 0; k < 4; k++ )
5128  graph_knot_add(auxg, -1);
5129  for( int k = 0; k < 4; k++ )
5130  for( int k2 = k + 1; k2 < 4; k2++ )
5131  graph_edge_add(scip, auxg, k, k2, 1.0, 1.0);
5132 
5133  SCIP_CALL( graph_path_init(scip, auxg) );
5134 
5135  /* try to eliminate non-positive vertices of degree 4 */
5136  for( int i = 0; i < nnodes; i++ )
5137  {
5138  int e;
5139  int k;
5140 
5141  /* only non-positive vertices of degree 4 */
5142  if( !g->mark[i] || g->grad[i] != 4 || Is_term(g->term[i]) )
5143  continue;
5144  k = 0;
5145 
5146  /* store neighbours */
5147  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
5148  adjverts[k++] = g->head[e];
5149 
5150  assert(k == 4);
5151 
5152  g->mark[i] = FALSE;
5153  prize = g->prize[i];
5154 
5155  /* compute mw bottleneck distance to each pair of neighbours */
5156  for( k = 0; k < 4; k++ )
5157  {
5158  auxg->mark[k] = TRUE;
5159  for( e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
5160  {
5161  const int k2 = auxg->head[e];
5162  if( k2 > k )
5163  {
5164  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &(sdist0), -prize, heap, statetail, statehead, memlbltail, memlblhead, adjverts[k], adjverts[k2], limit, FALSE, TRUE) );
5165  auxg->cost[e] = sdist0;
5166  if( SCIPisGT(scip, prize, -auxg->cost[e]) )
5167  break;
5168 
5169  auxg->cost[flipedge(e)] = auxg->cost[e];
5170  }
5171  }
5172  if( e != EAT_LAST )
5173  break;
5174  }
5175 
5176  k = UNKNOWN;
5177  if( e == EAT_LAST )
5178  {
5179  /* compute mst on all neighbours */
5180  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
5181 
5182  /* calculate mst cost */
5183  sdist0 = 0.0;
5184  for( int l = 1; l < 4; l++ )
5185  sdist0 += mst[l].dist;
5186 
5187  if( SCIPisLE(scip, prize, -sdist0) )
5188  {
5189  /* compute subset msts on all neighbours */
5190  for( k = 0; k < 4; k++ )
5191  {
5192  auxg->mark[k] = FALSE;
5193  sdist0 = 0.0;
5194  if( k != 0 )
5195  {
5196  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
5197  for( int l = 1; l < 4; l++ )
5198  if( auxg->mark[l] )
5199  sdist0 += mst[l].dist;
5200  }
5201  else
5202  {
5203  graph_path_exec(scip, auxg, MST_MODE, 1, auxg->cost, mst);
5204  for( int l = 2; l < 4; l++ )
5205  sdist0 += mst[l].dist;
5206  }
5207  auxg->mark[k] = TRUE;
5208  if( SCIPisGT(scip, prize, -sdist0) )
5209  break;
5210  }
5211  }
5212  }
5213 
5214  /* can node be eliminated? */
5215  if( k == 4 )
5216  {
5217  SCIPdebugMessage("npv4Reduction delete: %d (prize: %f) \n", i, g->prize[i]);
5218  (*nelims) += g->grad[i];
5219  while( g->outbeg[i] != EAT_LAST )
5220  graph_edge_del(scip, g, g->outbeg[i], TRUE);
5221  }
5222  else
5223  {
5224  g->mark[i] = TRUE;
5225  }
5226  }
5227 
5228  /* --- NPV5 test --- */
5229 
5230  /* enlarge graph for NPV5 test*/
5231  graph_knot_add(auxg, -1);
5232  for( int k = 0; k < 4; k++ )
5233  graph_edge_add(scip, auxg, k, 4, 1.0, 1.0);
5234  graph_path_exit(scip, auxg);
5235  SCIP_CALL( graph_path_init(scip, auxg) );
5236 
5237  /* try to eliminate non-positive vertices of degree 5 */
5238  for( int i = 0; i < nnodes; i++ )
5239  {
5240  int e;
5241  int k;
5242 
5243  /* only non-positive vertices of degree 5 */
5244  if( !g->mark[i] || g->grad[i] != 5 || Is_term(g->term[i]) )
5245  continue;
5246  k = 0;
5247 
5248  /* store neighbours */
5249  for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
5250  adjverts[k++] = g->head[e];
5251 
5252  assert(k == 5);
5253 
5254  g->mark[i] = FALSE;
5255  prize = g->prize[i];
5256 
5257  for( k = 0; k < 5; k++ )
5258  {
5259  auxg->mark[k] = TRUE;
5260  for( e = auxg->outbeg[k]; e != EAT_LAST; e = auxg->oeat[e] )
5261  {
5262  const int k2 = auxg->head[e];
5263  if( k2 > k )
5264  {
5265  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &(sdist0), -prize, heap, statetail, statehead, memlbltail, memlblhead, adjverts[k], adjverts[k2], limit, FALSE, TRUE) );
5266  auxg->cost[e] = sdist0;
5267  if( SCIPisGT(scip, prize, -auxg->cost[e]) )
5268  break;
5269 
5270  auxg->cost[flipedge(e)] = auxg->cost[e];
5271  }
5272  }
5273  if( e != EAT_LAST )
5274  break;
5275  }
5276  k = UNKNOWN;
5277  if( e == EAT_LAST )
5278  {
5279  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
5280  sdist0 = 0.0;
5281  for( int l = 1; l < 5; l++ )
5282  sdist0 += mst[l].dist;
5283 
5284  if( SCIPisLE(scip, prize, -sdist0) )
5285  {
5286  for( k = 0; k < 5; k++ )
5287  {
5288  int k2 = UNKNOWN;
5289  auxg->mark[k] = FALSE;
5290  sdist0 = 0.0;
5291  if( k != 0 )
5292  {
5293  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
5294  for( int l = 1; l < 5; l++ )
5295  if( auxg->mark[l] )
5296  sdist0 += mst[l].dist;
5297  }
5298  else
5299  {
5300  graph_path_exec(scip, auxg, MST_MODE, 1, auxg->cost, mst);
5301  for( int l = 2; l < 5; l++ )
5302  sdist0 += mst[l].dist;
5303  }
5304 
5305  if( SCIPisLE(scip, prize, -sdist0) )
5306  {
5307  for( k2 = k + 1; k2 < 5; k2++ )
5308  {
5309  if( k2 == k )
5310  continue;
5311  auxg->mark[k2] = FALSE;
5312  sdist0 = 0.0;
5313  if( k2 != 0 && k != 0)
5314  {
5315  graph_path_exec(scip, auxg, MST_MODE, 0, auxg->cost, mst);
5316  for( int l = 1; l < 5; l++ )
5317  if( auxg->mark[l] )
5318  sdist0 += mst[l].dist;
5319  }
5320  else
5321  {
5322  int s;
5323  if( k != 1 && k2 != 1 )
5324  s = 1;
5325  else
5326  s = 2;
5327  graph_path_exec(scip, auxg, MST_MODE, s, auxg->cost, mst);
5328  for( int l = 0; l < 5; l++ )
5329  if( auxg->mark[l] && l != s )
5330  sdist0 += mst[l].dist;
5331  }
5332  auxg->mark[k2] = TRUE;
5333  if( SCIPisGT(scip, prize, -sdist0) )
5334  break;
5335  }
5336  }
5337  auxg->mark[k] = TRUE;
5338  if( k2 != 5 )
5339  break;
5340  }
5341  }
5342  }
5343 
5344  if( k == 5 )
5345  {
5346  SCIPdebugMessage(" \n npv5Reduction delete: %d (prize: %f) \n", i, g->prize[i]);
5347  (*nelims) += g->grad[i];
5348  while( g->outbeg[i] != EAT_LAST )
5349  graph_edge_del(scip, g, g->outbeg[i], TRUE);
5350 
5351  }
5352  else
5353  {
5354  g->mark[i] = TRUE;
5355  }
5356  }
5357 
5358  /* free memory*/
5359  graph_path_exit(scip, auxg);
5360  graph_free(scip, &auxg, TRUE);
5361 
5362  return SCIP_OKAY;
5363 }
5364 
5365 
5366 /** chain reduction (NPV_2) for the MWCSP */
5367 SCIP_RETCODE reduce_chain2(
5368  SCIP* scip,
5369  GRAPH* g,
5370  PATH* pathtail,
5371  PATH* pathhead,
5372  int* heap,
5373  int* statetail,
5374  int* statehead,
5375  int* memlbltail,
5376  int* memlblhead,
5377  int* nelims,
5378  int limit
5379  )
5380 {
5381  SCIP_Real sdist;
5382  int i;
5383  int i1;
5384  int i2;
5385  int e1;
5386  int e2;
5387  int nnodes;
5388 
5389  assert(g != NULL);
5390  assert(scip != NULL);
5391  assert(pathtail != NULL);
5392  assert(pathhead != NULL);
5393  assert(heap != NULL);
5394  assert(statetail != NULL);
5395  assert(statehead != NULL);
5396  assert(memlbltail != NULL);
5397  assert(memlblhead != NULL);
5398  assert(nelims != NULL);
5399 
5400  *nelims = 0;
5401  nnodes = g->knots;
5402 
5403  /* initialize arrays */
5404  for( i = 0; i < nnodes; i++ )
5405  {
5406  statetail[i] = UNKNOWN;
5407  pathtail[i].dist = FARAWAY;
5408  pathtail[i].edge = UNKNOWN;
5409  statehead[i] = UNKNOWN;
5410  pathhead[i].dist = FARAWAY;
5411  pathhead[i].edge = UNKNOWN;
5412  }
5413 
5414  for( i = 0; i < nnodes; i++ )
5415  {
5416  assert(g->grad[i] >= 0);
5417  if( !g->mark[i] || g->grad[i] == 0 || Is_term(g->term[i]) || g->grad[i] != 2 )
5418  continue;
5419 
5420  /* non-positive chains */
5421  e1 = g->outbeg[i];
5422  e2 = g->oeat[e1];
5423  i1 = g->head[e1];
5424  i2 = g->head[e2];
5425 
5426  assert(e1 >= 0);
5427  assert(e2 >= 0);
5428  assert(g->mark[i1]);
5429  assert(g->mark[i2]);
5430  g->mark[i] = FALSE;
5431  SCIP_CALL( reduce_getSd(scip, g, pathtail, pathhead, &sdist, -(g->prize[i]), heap, statetail, statehead, memlbltail, memlblhead, i1, i2, limit, FALSE, TRUE) );
5432  if( SCIPisGE(scip, -sdist, g->prize[i]) )
5433  {
5434  SCIPdebugMessage("delete : %d prize: %f sd: %f \n", i, g->prize[i], -sdist );
5435  graph_edge_del(scip, g, e1, TRUE);
5436  graph_edge_del(scip, g, e2, TRUE);
5437  (*nelims) += 2;
5438  }
5439  else
5440  {
5441  g->mark[i] = TRUE;
5442  }
5443  }
5444  return SCIP_OKAY;
5445 }
5446 
5447 
5448 
5449 /** non-negative path reduction for the MWCSP */
5450 void reduce_nnp(
5451  SCIP* scip, /**< SCIP data structure */
5452  GRAPH* g, /**< graph data structure */
5453  int* marked, /**< nodes array */
5454  int* count /**< pointer to number of reductions */
5455  )
5456 {
5457  int localcount = 0;
5458  const int nnodes = g->knots;
5459 
5460  assert(scip != NULL);
5461  assert(g != NULL);
5462  assert(count != NULL);
5463  assert(marked != NULL);
5464  assert(g->stp_type == STP_MWCSP);
5465 
5466  /* unmark all nodes */
5467  for( int k = 0; k < nnodes; k++ )
5468  marked[k] = FALSE;
5469 
5470  /* check neighborhood of terminals */
5471  for( int k = 0; k < nnodes; k++ )
5472  {
5473  if( !g->mark[k] || g->prize[k] < 0.0 )
5474  continue;
5475 
5476  /* mark adjacent vertices of k */
5477  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
5478  if( g->mark[g->head[e]] )
5479  marked[g->head[e]] = TRUE;
5480 
5481  /* ... and traverse them */
5482  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
5483  {
5484  const int j = g->head[e];
5485 
5486  if( marked[j] )
5487  {
5488  int e2 = g->outbeg[j];
5489 
5490  while( e2 != EAT_LAST )
5491  {
5492  const int candedge = e2;
5493  e2 = g->oeat[e2];
5494  if( marked[g->head[candedge]] )
5495  {
5496  graph_edge_del(scip, g, candedge, TRUE);
5497  localcount++;
5498  }
5499  }
5500  }
5501  }
5502 
5503  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
5504  marked[g->head[e]] = FALSE;
5505 
5506  for( int j = 0; j < nnodes; j++ )
5507  assert(marked[j] == FALSE);
5508  }
5509 
5510  *count = localcount;
5511 
5512  assert(graph_valid(g));
5513 }
graph_path_init
SCIP_RETCODE graph_path_init(SCIP *, GRAPH *)
Definition: grphpath.c:444
SCIPisEQ
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:436
nterms
static volatile int nterms
Definition: interrupt.c:37
SCIPisPositive
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:561
getlecloseterms
static int getlecloseterms(SCIP *scip, PATH *vnoi, SCIP_Real *termdist, SCIP_Real ecost, int *vbase, int *neighbterms, int i, int nnodes)
Definition: reduce_alt.c:257
NULL
#define NULL
Definition: def.h:253
misc_stp_maxReal
SCIP_Real misc_stp_maxReal(SCIP_Real *realarr, unsigned nreals)
Definition: misc_stp.c:41
STP_BDR_MAXDEGREE
#define STP_BDR_MAXDEGREE
Definition: reduce_alt.c:2743
reduce_getSd
SCIP_RETCODE reduce_getSd(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, SCIP_Real *sdist, SCIP_Real distlimit, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int i, int i2, int limit, SCIP_Bool pc, SCIP_Bool mw)
Definition: reduce_alt.c:1970
GRAPH::head
int *RESTRICT head
Definition: grph.h:96
STP_BD_MAXDEGREE
#define STP_BD_MAXDEGREE
Definition: reduce_alt.c:2952
GRAPH
Definition: grph.h:57
reduce_npv
SCIP_RETCODE reduce_npv(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit)
Definition: reduce_alt.c:5024
GRAPH::source
int source
Definition: grph.h:67
SCIPisGE
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:488
graph_voronoiWithDist
SCIP_RETCODE graph_voronoiWithDist(SCIP *, const GRAPH *, SCIP_Real *, double *, int *, int *, int *, int *, int *, int *, PATH *)
shortest_path::edge
signed int edge
Definition: grph.h:146
MST_MODE
#define MST_MODE
Definition: grph.h:162
SCIPqueueInsert
SCIP_RETCODE SCIPqueueInsert(SCIP_QUEUE *queue, void *elem)
Definition: misc.c:1018
VERTEX_CONNECT
#define VERTEX_CONNECT
Definition: reduce_alt.c:43
graph_valid
SCIP_Bool graph_valid(const GRAPH *)
Definition: grphbase.c:4324
GRAPH::terms
int terms
Definition: grph.h:64
reduce_sdsp
SCIP_RETCODE reduce_sdsp(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit, int *edgestate)
Definition: reduce_alt.c:2459
reduce_getSdPcMw
SCIP_RETCODE reduce_getSdPcMw(SCIP *scip, const GRAPH *g, PATH *pathtail, PATH *pathhead, SCIP_Real *sdist, SCIP_Real distlimit, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *pathmaxnodetail, int *pathmaxnodehead, int i, int i2, int limit)
Definition: reduce_alt.c:2110
reduce_nts
SCIP_RETCODE reduce_nts(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit)
Definition: reduce_alt.c:3217
EAT_LAST
#define EAT_LAST
Definition: grph.h:31
Edge_anti
#define Edge_anti(a)
Definition: grph.h:171
SCIPqueueRemove
void * SCIPqueueRemove(SCIP_QUEUE *queue)
Definition: misc.c:1069
FALSE
#define FALSE
Definition: def.h:73
VERTEX_TEMPNEIGHBOR
#define VERTEX_TEMPNEIGHBOR
Definition: reduce_alt.c:44
GRAPH::inpbeg
int *RESTRICT inpbeg
Definition: grph.h:74
GRAPH::path_state
int *RESTRICT path_state
Definition: grph.h:119
graph_free
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: grphbase.c:3674
probdata_stp.h
Problem data for stp problem.
TRUE
#define TRUE
Definition: def.h:72
SCIPdebug
#define SCIPdebug(x)
Definition: pub_message.h:74
SCIP_RETCODE
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
sddeltable
static SCIP_Bool sddeltable(SCIP *scip, GRAPH *g, PATH *path, int *vbase, int *forbidden, int tpos, int hpos, int tail, int head, int edge, int nnodes)
Definition: reduce_alt.c:55
graph_path_exit
void graph_path_exit(SCIP *, GRAPH *)
Definition: grphpath.c:466
STP_DHCSTP
#define STP_DHCSTP
Definition: grph.h:48
shortest_path
Definition: grph.h:143
reduce_cnsAdv
SCIP_RETCODE reduce_cnsAdv(SCIP *scip, GRAPH *g, int *marked, int *count)
Definition: reduce_alt.c:4748
graph_knot_delPseudo
SCIP_RETCODE graph_knot_delPseudo(SCIP *, GRAPH *, const SCIP_Real *, const SCIP_Real *, const SCIP_Real *, int, SCIP_Bool *)
Definition: grphbase.c:2258
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
STP_PCSPG
#define STP_PCSPG
Definition: grph.h:40
SCIPdebugMessage
#define SCIPdebugMessage
Definition: pub_message.h:77
STP_RED_CNSNN
#define STP_RED_CNSNN
Definition: reduce_alt.c:47
SCIPfreeBufferArray
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
graph_knot_chg
void graph_knot_chg(GRAPH *, int, int)
Definition: grphbase.c:2218
reduce_sdspSap
SCIP_RETCODE reduce_sdspSap(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit)
Definition: reduce_alt.c:2308
getRSD
static SCIP_Real getRSD(SCIP *scip, GRAPH *g, GRAPH *netgraph, PATH *mst, PATH *vnoi, SCIP_Real *mstsdist, SCIP_Real *termdist1, SCIP_Real *termdist2, SCIP_Real sd_initial, int *vbase, int *nodesid, int *neighbterms1, int *neighbterms2, int i, int i2, int limit)
Definition: reduce_alt.c:1777
STP_BDR_MAXDNEDGES
#define STP_BDR_MAXDNEDGES
Definition: reduce_alt.c:2744
GRAPH::mark
int *RESTRICT mark
Definition: grph.h:70
VERTEX_NEIGHBOR
#define VERTEX_NEIGHBOR
Definition: reduce_alt.c:45
GRAPH::fixedges
IDX * fixedges
Definition: grph.h:85
graph_pc_contractEdge
SCIP_RETCODE graph_pc_contractEdge(SCIP *, GRAPH *, int *, int, int, int)
Definition: grphbase.c:2101
SCIPisLE
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:462
graph_sdPaths
void graph_sdPaths(SCIP *, const GRAPH *, PATH *, SCIP_Real *, SCIP_Real, int *, int *, int *, int *, int, int, int)
Definition: grphpath.c:567
STP_RED_ANSMAXNEIGHBORS
#define STP_RED_ANSMAXNEIGHBORS
Definition: reduce_alt.c:50
reduce_ansAdv2
void reduce_ansAdv2(SCIP *scip, GRAPH *g, int *marked, int *count)
Definition: reduce_alt.c:4621
SCIPqueueIsEmpty
SCIP_Bool SCIPqueueIsEmpty(SCIP_QUEUE *queue)
Definition: misc.c:1173
GRAPH::oeat
int *RESTRICT oeat
Definition: grph.h:103
LE
#define LE(a, b)
Definition: portab.h:65
CONNECT
#define CONNECT
Definition: grph.h:154
reduce_ans
void reduce_ans(SCIP *scip, GRAPH *g, int *marked, int *count)
Definition: reduce_alt.c:4447
misc_stp.h
miscellaneous methods used for solving Steiner problems
SCIPintListNodeAppendCopy
SCIP_RETCODE SCIPintListNodeAppendCopy(SCIP *scip, IDX **node1, IDX *node2, SCIP_Bool *conflict)
Definition: misc_stp.c:94
STP_SAP
#define STP_SAP
Definition: grph.h:39
GRAPH::stp_type
int stp_type
Definition: grph.h:127
GRAPH::ancestors
IDX ** ancestors
Definition: grph.h:86
level0
SCIP_RETCODE level0(SCIP *, GRAPH *)
Definition: reduce.c:153
reduce_bdr
SCIP_RETCODE reduce_bdr(SCIP *scip, GRAPH *g, GRAPH *netgraph, PATH *netmst, PATH *vnoi, SCIP_Real *mstsdist, SCIP_Real *termdist1, SCIP_Real *termdist2, int *vbase, int *nodesid, int *neighbterms1, int *neighbterms2, int *nelims)
Definition: reduce_alt.c:2747
Is_pterm
#define Is_pterm(a)
Definition: grph.h:169
STP_Bool
unsigned char STP_Bool
Definition: grph.h:52
SCIPfreeBufferArrayNull
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:124
GRAPH::prize
SCIP_Real * prize
Definition: grph.h:82
SCIP_Queue
Definition: struct_misc.h:51
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
graph_path_exec
void graph_path_exec(SCIP *, const GRAPH *, const int, int, const SCIP_Real *, PATH *)
Definition: grphpath.c:491
SCIP_OKAY
Definition: type_retcode.h:33
GRAPH::knots
int knots
Definition: grph.h:62
reduce_sdPc
SCIP_RETCODE reduce_sdPc(SCIP *scip, GRAPH *g, PATH *vnoi, int *heap, int *state, int *vbase, int *nodesid, int *nodesorg, int *nelims)
Definition: reduce_alt.c:1499
SCIP_CALL
#define SCIP_CALL(x)
Definition: def.h:365
SCIPqueueCreate
SCIP_RETCODE SCIPqueueCreate(SCIP_QUEUE **queue, int initsize, SCIP_Real sizefac)
Definition: misc.c:932
LT
#define LT(a, b)
Definition: portab.h:64
graph_path_PcMwSd
void graph_path_PcMwSd(SCIP *, const GRAPH *, PATH *, SCIP_Real *, SCIP_Real, int *, int *, int *, int *, int *, int *, int, int, int)
Definition: grphpath.c:664
reduce_nvAdv
SCIP_RETCODE reduce_nvAdv(SCIP *scip, GRAPH *g, PATH *vnoi, SCIP_Real *distance, double *fixed, int *edgearrint, int *heap, int *state, int *vbase, int *neighb, int *distnode, int *solnode, int *nelims)
Definition: reduce_alt.c:3903
FARAWAY
#define FARAWAY
Definition: grph.h:156
SCIPgetRealParam
SCIP_RETCODE SCIPgetRealParam(SCIP *scip, const char *name, SCIP_Real *value)
Definition: scip_param.c:297
voronoi_term
void voronoi_term(const GRAPH *, double *, double *, double *, PATH *, int *, int *, int *, int *, int)
STP_SPG
#define STP_SPG
Definition: grph.h:38
SCIPallocBufferArray
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
GT
#define GT(a, b)
Definition: portab.h:66
reduce_nnp
void reduce_nnp(SCIP *scip, GRAPH *g, int *marked, int *count)
Definition: reduce_alt.c:5450
SCIP_Bool
#define SCIP_Bool
Definition: def.h:70
SCIPgetTotalTime
SCIP_Real SCIPgetTotalTime(SCIP *scip)
Definition: scip_timing.c:332
GRAPH::ieat
int *RESTRICT ieat
Definition: grph.h:102
STP_MWCSP
#define STP_MWCSP
Definition: grph.h:47
GRAPH::tail
int *RESTRICT tail
Definition: grph.h:95
ansProcessCandidate
static void ansProcessCandidate(SCIP *scip, GRAPH *g, int *marked, int *count, SCIP_Real min, int candvertex)
Definition: reduce_alt.c:1028
getcloseterms
static int getcloseterms(SCIP *scip, const PATH *vnoi, SCIP_Real *termdist, SCIP_Real ecost, const int *vbase, int *neighbterms, int i, int nnodes)
Definition: reduce_alt.c:167
MIN
#define MIN(x, y)
Definition: def.h:223
SCIPisNegative
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
Definition: scip_numerics.c:573
STP_BD_MAXDNEDGES
#define STP_BD_MAXDNEDGES
Definition: reduce_alt.c:2953
GRAPH::term
int *RESTRICT term
Definition: grph.h:68
graph_edge_del
void graph_edge_del(SCIP *, GRAPH *, int, SCIP_Bool)
Definition: grphbase.c:2837
reduce_nv
SCIP_RETCODE reduce_nv(SCIP *scip, GRAPH *g, PATH *vnoi, double *fixed, int *edgearrint, int *heap, int *state, int *vbase, int *solnode, int *nelims)
Definition: reduce_alt.c:3710
reduce_sl
SCIP_RETCODE reduce_sl(SCIP *scip, GRAPH *g, PATH *vnoi, double *fixed, int *heap, int *state, int *vbase, int *vrnodes, STP_Bool *visited, int *solnode, int *nelims)
Definition: reduce_alt.c:3483
getcloseterms2term
static int getcloseterms2term(SCIP *scip, const GRAPH *g, const GRAPH *netgraph, SCIP_Real *termdist, SCIP_Real ecost, int *neighbterms, int *nodesid, int *nodesorg, int i)
Definition: reduce_alt.c:206
grph.h
includes various files containing graph methods used for Steiner tree problems
graph_get4nextTerms
void graph_get4nextTerms(SCIP *, const GRAPH *, SCIP_Real *, SCIP_Real *, PATH *, int *, int *, int *)
Definition: grphpath.c:2186
b
SCIP_VAR ** b
Definition: circlepacking.c:56
reduce_ansAdv
void reduce_ansAdv(SCIP *scip, GRAPH *g, int *marked, int *count, SCIP_Bool extneigbhood)
Definition: reduce_alt.c:4515
Is_term
#define Is_term(a)
Definition: grph.h:168
MAX
#define MAX(x, y)
Definition: def.h:222
EDGE_BLOCKED
#define EDGE_BLOCKED
Definition: grph.h:159
EAT_FREE
#define EAT_FREE
Definition: grph.h:30
reduce_bd34
SCIP_RETCODE reduce_bd34(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit)
Definition: reduce_alt.c:2964
GRAPH::cost
SCIP_Real * cost
Definition: grph.h:94
a
SCIP_VAR * a
Definition: circlepacking.c:57
SCIP_Real
#define SCIP_Real
Definition: def.h:164
reduce_sd
SCIP_RETCODE reduce_sd(SCIP *scip, GRAPH *g, PATH *vnoi, SCIP_Real *edgepreds, SCIP_Real *mstsdist, int *heap, int *state, int *vbase, int *nodesid, int *nodesorg, int *forbidden, int *nelims, SCIP_Bool nodereplacing, int *edgestate)
Definition: reduce_alt.c:1105
SCIPisLT
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:449
SCIPisGT
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
Definition: scip_numerics.c:475
GRAPH::outbeg
int *RESTRICT outbeg
Definition: grph.h:76
SCIP_Longint
#define SCIP_Longint
Definition: def.h:149
GRAPH::edges
int edges
Definition: grph.h:92
flipedge
#define flipedge(edge)
Definition: grph.h:150
graph_knot_add
void graph_knot_add(GRAPH *, int)
Definition: grphbase.c:2196
graph_get4nextTTerms
SCIP_RETCODE graph_get4nextTTerms(SCIP *, const GRAPH *, SCIP_Real *, PATH *, int *, int *, int *)
Definition: grphpath.c:2227
graph_knot_contract
SCIP_RETCODE graph_knot_contract(SCIP *, GRAPH *, int *, int, int)
Definition: grphbase.c:2429
UNKNOWN
#define UNKNOWN
Definition: sepa_mcf.c:4081
STP_RED_ANSEXMAXCANDS
#define STP_RED_ANSEXMAXCANDS
Definition: reduce_alt.c:49
nnodes
#define nnodes
Definition: gastrans.c:65
reduce_chain2
SCIP_RETCODE reduce_chain2(SCIP *scip, GRAPH *g, PATH *pathtail, PATH *pathhead, int *heap, int *statetail, int *statehead, int *memlbltail, int *memlblhead, int *nelims, int limit)
Definition: reduce_alt.c:5367
VERTEX_OTHER
#define VERTEX_OTHER
Definition: reduce_alt.c:46
SCIPqueueFree
void SCIPqueueFree(SCIP_QUEUE **queue)
Definition: misc.c:956
BMSclearMemoryArray
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:120
reduce_ledge
SCIP_RETCODE reduce_ledge(SCIP *scip, GRAPH *g, PATH *vnoi, int *heap, int *state, int *vbase, int *nelims, int *edgestate)
Definition: reduce_alt.c:4143
nearest
static int nearest(int *heap, int *state, int *count, const PATH *path)
Definition: grphpath.c:43
GRAPH::hoplimit
int hoplimit
Definition: grph.h:89
graph_edge_add
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
STP_RPCSPG
#define STP_RPCSPG
Definition: grph.h:41
graph_voronoiTerms
void graph_voronoiTerms(SCIP *, const GRAPH *, const SCIP_Real *, PATH *, int *, int *, int *)
Definition: grphpath.c:2307
scip.h
SCIP callable library.
graph_init
SCIP_RETCODE graph_init(SCIP *, GRAPH **, int, int, int)
Definition: grphbase.c:3491
STP_RED_ANSMAXCANDS
#define STP_RED_ANSMAXCANDS
Definition: reduce_alt.c:48