Scippy

SCIP

Solving Constraint Integer Programs

reduce_bnd.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-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file reduce_bnd.c
17  * @brief bound based reductions for Steiner tree problems
18  * @author Daniel Rehfeldt
19  *
20  * This file implements bound-based reduction techniques for several Steiner problems.
21  * All tests can be found in "A Generic Approach to Solving the Steiner Tree Problem and Variants" by Daniel Rehfeldt.
22  *
23  * A list of all interface methods can be found in grph.h.
24  *
25  */
26 
27 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31 #include <assert.h>
32 #include "scip/scip.h"
33 #include "grph.h"
34 #include "heur_tm.h"
35 #include "heur_ascendprune.h"
36 #include "cons_stp.h"
37 #include "heur_local.h"
38 #include "misc_stp.h"
39 #include "prop_stp.h"
40 #include "probdata_stp.h"
41 #include "heur_rec.h"
42 #include "heur_slackprune.h"
43 
44 #define DEFAULT_HEURRUNS 100 /**< number of runs of constructive heuristic */
45 #define DEFAULT_DARUNS 7 /**< number of runs for dual ascent heuristic */
46 #define DEFAULT_NMAXROOTS 8 /**< max number of roots to use for new graph in dual ascent heuristic */
47 #define PERTUBATION_RATIO 0.05 /**< pertubation ratio for dual-ascent primal bound computation */
48 #define PERTUBATION_RATIO_PC 0.005 /**< pertubation ratio for dual-ascent primal bound computation */
49 #define SOLPOOL_SIZE 20 /**< size of presolving solution pool */
50 #define STP_RED_MINBNDTERMS 750
51 #define STP_DABD_MAXDEGREE 5
52 #define STP_DABD_MAXDNEDGES 10
53 #define STP_DAEX_MAXDFSDEPTH 6
54 #define STP_DAEX_MINDFSDEPTH 4
55 #define STP_DAEX_MAXGRAD 8
56 #define STP_DAEX_MINGRAD 6
57 #define STP_DAEX_EDGELIMIT 50000
58 #define DAMAXDEVIATION_RANDOM_LOWER 0.15 /**< random upper bound for max deviation for dual ascent */
59 #define DAMAXDEVIATION_RANDOM_UPPER 0.30 /**< random upper bound for max deviation for dual ascent */
60 #define DAMAXDEVIATION_FAST 0.75
61 
62 #define EXEDGE_FREE 0
63 #define EXEDGE_FIXED 1
64 #define EXEDGE_KILLED 2
65 
66 
67 /** mark ancestors of given edge */
68 static
70  const GRAPH* graph, /**< graph data structure */
71  int edge, /**< edge to use */
72  int* ancestormark /**< ancestor mark array */
73  )
74 {
75  int count = 0;
76  assert(edge >= 0);
77 
78  for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
79  {
80  const unsigned idx = ((unsigned) curr->index) / 2;
81 
82  assert(curr->index >= 0 && idx < (unsigned) (MAX(graph->edges, graph->orgedges) / 2));
83 
84  if( ancestormark[idx] )
85  return TRUE;
86 
87  ancestormark[idx] = 1;
88  }
89 
90  return FALSE;
91 }
92 
93 
94 /** mark ancestors of given edge */
95 static
97  const GRAPH* graph, /**< graph data structure */
98  int edge, /**< edge to use */
99  int* ancestormark /**< ancestor mark array */
100  )
101 {
102  int count = 0;
103  assert(edge >= 0);
104 
105  for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
106  {
107  assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
108  assert((ancestormark[((unsigned) curr->index) / 2]) == 0);
109 
110  ancestormark[((unsigned) curr->index) / 2] = 1;
111  }
112 }
113 
114 /** unmark ancestors of given edge */
115 static
117  const GRAPH* graph, /**< graph data structure */
118  int edge, /**< edge to use */
119  int* ancestormark /**< ancestor mark array */
120  )
121 {
122  int count = 0;
123  assert(edge >= 0);
124 
125  for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
126  {
127  assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
128  ancestormark[((unsigned) curr->index) / 2] = 0;
129  }
130 }
131 
132 /** unmark ancestors of given edge */
133 static
135  const GRAPH* graph, /**< graph data structure */
136  int edge, /**< edge to use */
137  int* ancestormark /**< ancestor mark array */
138  )
139 {
140  int count = 0;
141  assert(edge >= 0);
142 
143  for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
144  {
145  const unsigned idx = ((unsigned) curr->index) / 2;
146 
147  assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
148  assert(ancestormark[idx] == 1);
149 
150  ancestormark[idx] = 0;
151  }
152 }
153 
154 
155 /** finalize subtree computations (clean up, update global bound) */
156 static
158  const GRAPH* graph, /**< graph data structure */
159  const int* edgeends, /**< heads or tail of edge */
160  const int* treeedges, /**< tree edges */
161  int dfsdepth, /**< dfs depth */
162  SCIP_Bool stopped, /**< has extension been stopped? */
163  SCIP_Real localbound, /**< local bound */
164  SCIP_Real* globalbound, /**< points to global bound */
165  SCIP_Bool* ruleout, /**< rule out? */
166  int* nodepos, /**< node position in tree */
167  int* ancestormark /**< ancestor mark array */
168 )
169 {
170 
171  if( localbound > *globalbound )
172  *globalbound = localbound;
173 
174  if( !stopped )
175  *ruleout = TRUE;
176 
177  if( !(*ruleout) )
178  {
179  for( int etree = 0; etree < dfsdepth; etree++ )
180  {
181  const int node = edgeends[treeedges[etree]];
182  assert(nodepos[node]);
183  nodepos[node] = 0;
184  unmarkAncestors(graph, treeedges[etree], ancestormark);
185  }
186  }
187  else
188  {
189  assert(dfsdepth == 0);
190  }
191 }
192 
193 /** should we truncate subtree? */
194 static
196  const GRAPH* graph, /**< graph data structure */
197  SCIP_Real extendedcost, /**< cost of subtree */
198  int root, /**< DA root */
199  int currhead, /**< latest node */
200  int maxgrad, /**< maximum allowed degree */
201  int dfsdepth, /**< dfs depth */
202  int maxdfsdepth, /**< max dfs depth */
203  SCIP_Real* minbound, /**< bound */
204  SCIP_Bool* stopped /**< real truncation? */
205 )
206 {
207  if( Is_term(graph->term[currhead]) || graph->grad[currhead] > maxgrad || dfsdepth >= maxdfsdepth )
208  {
209  /* real truncation? */
210  if( root != currhead )
211  {
212  *stopped = TRUE;
213  if( extendedcost < *minbound )
214  *minbound = extendedcost;
215  }
216  return TRUE;
217  }
218  return FALSE;
219 }
220 #if 0
221 
222 /** extend subtree and return new nadded_edges */
223 static
224 int reduceAndExtendSubtree(
225  SCIP* scip, /**< SCIP */
226  const GRAPH* graph, /**< graph data structure */
227  const int* graph_start, /**< start */
228  const int* graph_next, /**< next */
229  const int* graph_endp, /**< next */
230  int currnode, /**< latest node */
231  int nadded_edges, /**< added edges */
232  int* nodepos, /**< node position in tree */
233  int* edgestack, /**< stack */
234  unsigned char* extensionmark /**< mark array for extension */
235 )
236 {
237  int n = nadded_edges;
238  const int start = n;
239  SCIP_Bool cleanup = FALSE;
240 
241  /* extend from currnode */
242  for( int e = graph_start[currnode]; e != EAT_LAST; e = graph_next[e] )
243  if( !nodepos[graph_endp[e]] )
244  {
245  /* check whether edge needs to be added */
246 
247  const int vertex_exnew = graph_endp[e];
248 
249  assert((n - start) <= STP_DAEX_MAXGRAD);
250  extensionmark[n - start] = EXEDGE_FREE;
251 
252  for( int e2 = graph->outbeg[vertex_exnew]; e2 != EAT_LAST; e2 = graph->oeat[e2] )
253  {
254  const int head = graph->head[e2];
255 
256  /* have we already extended to head? */
257  if( nodepos[head] < 0 )
258  {
259  const int stackposition = (-nodepos[head]) - 1;
260  const int edge_exold = edgestack[stackposition];
261  const SCIP_Real cost_exnew = graph->cost[e];
262  const SCIP_Real cost_exold = graph->cost[edge_exold];
263  const SCIP_Real cost_alt = graph->cost[e2];
264 
265  assert(edge_exold >= 0);
266  assert(e != e2 && e != edge_exold && edge_exold != e2);
267 
268  if( (SCIPisLT(scip, cost_alt, cost_exold) || SCIPisLT(scip, cost_alt, cost_exnew)) )
269  {
270  /* can we eliminate new edge? */
271  if( cost_exold <= cost_exnew && extensionmark[n - start] == EXEDGE_FREE )
272  {
273  assert(stackposition >= start);
274 
275  extensionmark[stackposition - start] = EXEDGE_FIXED;
276  extensionmark[n - start] = EXEDGE_KILLED;
277  break;
278  }
279 
280  /* can we eliminate old edge? */
281  if( cost_exold >= cost_exnew && extensionmark[stackposition - start] == EXEDGE_FREE )
282  {
283  extensionmark[stackposition - start] = EXEDGE_KILLED;
284  extensionmark[n - start] = EXEDGE_FIXED;
285 
286  assert(head == graph->head[edge_exold]);
287  nodepos[head] = 0;
288 
289  cleanup = TRUE;
290  }
291  }
292  }
293  }
294 
295  if( extensionmark[n - start] != EXEDGE_KILLED )
296  {
297  edgestack[n++] = e;
298 
299  assert(nodepos[vertex_exnew] == 0);
300  nodepos[vertex_exnew] = -n;
301  }
302  }
303 
304  for( int i = start; i < n; i++ )
305  {
306  const int edge = edgestack[i];
307 
308  assert(edge >= 0);
309  assert(nodepos[graph->head[edge]] <= 0);
310  assert(nodepos[graph->head[edge]] < 0 || (extensionmark[i - start] == EXEDGE_KILLED));
311 
312  nodepos[graph->head[edge]] = 0;
313  }
314 
315  if( cleanup )
316  {
317  int newn = start;
318  for( int i = start; i < n; i++ )
319  if( extensionmark[i - start] != EXEDGE_KILLED )
320  edgestack[newn++] = edgestack[i];
321 
322 
323  assert(newn < n);
324  n = newn;
325  }
326 
327 
328  return n;
329 }
330 #endif
331 
332 /** remove latest edge from subtree and returns new tree cost */
333 static
335  SCIP* scip, /**< SCIP */
336  const GRAPH* graph, /**< graph data structure */
337  const SCIP_Real* redcost, /**< reduced costs */
338  const int* treeedges, /**< edges of tree */
339  SCIP_Real treecostold, /**< old tree cost */
340  SCIP_Real treecostoffset, /**< tree cost offset */
341  int dfsdepth, /**< DFS depth */
342  int lastnode, /**< last node */
343  int* nodepos, /**< array to mark node position*/
344  int* ancestormark, /**< ancestor mark array */
345  unsigned int* nstallrounds /**< rounds since latest update */
346 )
347 {
348  const int lastedge = treeedges[dfsdepth - 1];
349 
350  assert(dfsdepth >= 1 && lastnode >= 0);
351  assert(SCIPisGE(scip, treecostold, 0.0) && treecostoffset >= 0.0);
352 
353  nodepos[lastnode] = 0;
354  unmarkAncestors(graph, lastedge, ancestormark);
355 
356  /* recompute tree cost? */
357  if( (*nstallrounds)++ >= 9 )
358  {
359  SCIP_Real treecost = treecostoffset;
360 
361  *nstallrounds = 0;
362 
363  for( int i = 0; i < dfsdepth - 1; i++ )
364  treecost += redcost[treeedges[i]];
365 
366 #if 0 // todo deleteeme
367  if( !SCIPisEQ(scip, treecost, treecostold - redcost[lastedge]) )
368  {
369  printf("%.15f %.15f \n", treecost, treecostold - redcost[lastedge]);
370  assert(0);
371  exit(1);
372  }
373 #endif
374 
375  assert(SCIPisEQ(scip, treecost, treecostold - redcost[lastedge]));
376  }
377 
378  return (treecostold - redcost[lastedge]);
379 }
380 
381 /** extend subtree and return new nadded_edges */
382 static
384  SCIP* scip, /**< SCIP */
385  const GRAPH* graph, /**< graph data structure */
386  const SCIP_Real* redcost, /**< reduced costs */
387  int curredge, /**< added edges */
388  int currhead, /**< latest node */
389  int dfsdepth, /**< DFS depth*/
390  int nadded_edges, /**< added edges */
391  SCIP_Real* treecost, /**< pointer to treecost */
392  SCIP_Real* treecosts, /**< edge costs of tree */
393  int* nodepos, /**< node position in tree */
394  int* treeedges, /**< edges of tree */
395  int* edgestack, /**< stack */
396  int* ancestormark /**< ancestor mark array */
397 )
398 {
399  int n = nadded_edges + 1;
400  nodepos[currhead] = dfsdepth + 1;
401  *treecost += redcost[curredge];
402  treeedges[dfsdepth] = curredge;
403  treecosts[dfsdepth] = graph->cost[curredge];
404 
405  markAncestors(graph, curredge, ancestormark);
406 
407  for( int e = graph->outbeg[currhead]; e != EAT_LAST; e = graph->oeat[e] )
408  if( !nodepos[graph->head[e]] )
409  edgestack[n++] = e;
410 
411  return n;
412 }
413 
414 /** extend subtree and return new nadded_edges */
415 static
417  const GRAPH* graph, /**< graph data structure */
418  const SCIP_Real* redcost, /**< reduced costs */
419  int curredge, /**< added edges */
420  int currtail, /**< latest node */
421  int dfsdepth, /**< DFS depth*/
422  int nadded_edges, /**< added edges */
423  SCIP_Real* treecost, /**< pointer to treecost */
424  SCIP_Real* treecosts, /**< edge costs of tree */
425  int* nodepos, /**< node position in tree */
426  int* treeedges, /**< edges of tree */
427  int* edgestack, /**< stack */
428  int* ancestormark /**< ancestor mark array */
429 )
430 {
431  int n = nadded_edges + 1;
432  nodepos[currtail] = dfsdepth + 1;
433  *treecost += redcost[curredge];
434  treeedges[dfsdepth] = curredge;
435  treecosts[dfsdepth] = graph->cost[curredge];
436 
437  markAncestors(graph, curredge, ancestormark);
438 
439  for( int e = graph->inpbeg[currtail]; e != EAT_LAST; e = graph->ieat[e] )
440  if( !nodepos[graph->tail[e]] )
441  edgestack[n++] = e;
442 
443  return n;
444 }
445 
446 
447 /** small helper */
448 static
450  int edge, /**< the edge */
451  unsigned int* eqstack, /**< stores edges that were used for equality comparison */
452  int* eqstack_size, /**< pointer to size of eqstack */
453  SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
454 )
455 {
456  const unsigned int halfedge = ((unsigned int) edge) / 2;
457 
458  assert(edge >= 0);
459 
460  if( !eqmark[halfedge] )
461  {
462  eqmark[halfedge] = TRUE;
463  eqstack[*eqstack_size] = halfedge;
464 
465  *eqstack_size = *eqstack_size + 1;
466  }
467 }
468 
469 
470 /** compare with tree bottleneck */
471 static
473  SCIP* scip, /**< SCIP */
474  const GRAPH* graph, /**< graph data structure */
475  const SCIP_Real* treecosts, /**< costs of edges of the tree */
476  const SCIP_Real* basebottlenecks, /**< bottleneck costs for innode and basenode */
477  const int* nodepos, /**< node position in tree */
478  const int* treeedges, /**< edges in tree (corresponding to treecosts) */
479  SCIP_Real orgedgecost, /**< cost of original edge */
480  SCIP_Real extedgecost, /**< cost of short-cut edge */
481  int orgedge, /**< original edge */
482  int extedge, /**< short-cut edge */
483  int dfsdepth, /**< dfs depth */
484  SCIP_Bool allow_eq, /**< test for equality? */
485  unsigned int* eqstack, /**< stores edges that were used for equality comparison */
486  int* eqstack_size, /**< pointer to size of eqstack */
487  SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
488 )
489 {
490  int start;
491 
492  assert(!allow_eq ||(eqstack != NULL && eqstack_size != NULL && eqmark != NULL));
493 
494  if( allow_eq )
495  {
496  if( SCIPisLE(scip, extedgecost, orgedgecost) )
497  {
498  if( SCIPisEQ(scip, extedgecost, orgedgecost) )
499  updateEqArrays(orgedge, eqstack, eqstack_size, eqmark);
500 
501  return TRUE;
502  }
503  }
504  else if( SCIPisLT(scip, extedgecost, orgedgecost) )
505  return TRUE;
506 
507  if( nodepos[graph->tail[extedge]] > graph->knots )
508  {
509  start = nodepos[graph->tail[extedge]] - 1 - graph->knots;
510  assert(start >= 0 && start <= 3);
511 
512  if( start <= 2 )
513  {
514  assert(basebottlenecks[start] > 0);
515 
516  if( SCIPisLT(scip, extedgecost, basebottlenecks[start]) )
517  return TRUE;
518  }
519  start = 0;
520  }
521  else
522  {
523  start = nodepos[graph->tail[extedge]]; /* not -1! We save the incoming bottlenecks */
524  assert(start >= 1 && start <= dfsdepth);
525  assert(start < dfsdepth || graph->tail[orgedge] == graph->tail[extedge]);
526 
527  }
528 
529  for( int i = start; i < dfsdepth; i++ )
530  {
531  assert(treecosts[i] >= 0.0);
532 
533  if( allow_eq )
534  {
535  if( SCIPisLE(scip, extedgecost, treecosts[i]) )
536  {
537  if( SCIPisEQ(scip, extedgecost, treecosts[i]) )
538  updateEqArrays(treeedges[i], eqstack, eqstack_size, eqmark);
539 
540  return TRUE;
541  }
542  }
543  else if( SCIPisLT(scip, extedgecost, treecosts[i]) )
544  return TRUE;
545  }
546 
547  return FALSE;
548 }
549 
550 /** can subtree be ruled out? */
551 static
553  SCIP* scip, /**< SCIP */
554  const GRAPH* graph, /**< graph data structure */
555  const SCIP_Real* treecosts, /**< costs of edges of the tree */
556  const SCIP_Real* basebottlenecks, /**< bottleneck costs for innode and basenode */
557  const int* nodepos, /**< node position in tree */
558  const int* treeedges, /**< edges in tree (corresponding to treecosts) */
559  SCIP_Real cutoff, /**< cut-off bound */
560  SCIP_Real extendedcost, /**< cost of subtree */
561  int dfsdepth, /**< dfs depth */
562  int curredge, /**< latest edge */
563  SCIP_Bool allowequality, /**< rule out also in case of equality */
564  const int* ancestormark, /**< ancestor mark array */
565  unsigned int* eqstack, /**< stores edges that were used for equality comparison */
566  int* eqstack_size, /**< pointer to size of eqstack */
567  SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
568 )
569 {
570  if( allowequality ? (extendedcost >= cutoff) : SCIPisGT(scip, extendedcost, cutoff) )
571  {
572  return TRUE;
573  }
574  else
575  {
576  SCIP_Real currcost;
577  int currhead;
578 
579  if( ancestormark != NULL )
580  {
581  int count = 0;
582  for( IDX* curr = graph->ancestors[curredge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
583  {
584  assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
585  if( ancestormark[((unsigned) curr->index) / 2] )
586  return TRUE;
587  }
588  }
589 
590  currcost = graph->cost[curredge];
591  currhead = graph->head[curredge];
592 
593  for( int e = graph->inpbeg[currhead]; e != EAT_LAST; e = graph->ieat[e] )
594  {
595  int tail;
596  SCIP_Real ecost;
597 
598  if( e == curredge )
599  continue;
600 
601  assert(flipedge(e) != curredge);
602 
603  tail = graph->tail[e];
604  ecost = graph->cost[e];
605 
606  if( nodepos[tail] )
607  {
608  const SCIP_Bool allow_eq = (eqmark != NULL && eqmark[((unsigned) e) / 2] == FALSE);
609 
610  if( bottleneckRuleOut(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, currcost,
611  ecost, curredge, e, dfsdepth, allow_eq, eqstack, eqstack_size, eqmark) )
612  return TRUE;
613  }
614  else
615  {
616  const SCIP_Bool is_term = Is_term(graph->term[tail]);
617 
618  for( int e2 = graph->outbeg[tail], count = 0; e2 != EAT_LAST && count < STP_DAEX_MAXGRAD;
619  e2 = graph->oeat[e2], count++ )
620  {
621  int head2;
622 
623  if( e2 == e )
624  continue;
625 
626  assert(flipedge(e2) != e && e2 != curredge && e2 != flipedge(curredge) );
627 
628  head2 = graph->head[e2];
629  assert(head2 != tail);
630 
631  if( nodepos[head2] )
632  {
633  const SCIP_Real extcost = is_term ? MAX(ecost, graph->cost[e2]) : (ecost + graph->cost[e2]);
634  const SCIP_Bool allow_eq = (eqmark != NULL && eqmark[((unsigned) e) / 2] == FALSE && eqmark[((unsigned) e2) / 2] == FALSE);
635 
636  assert(head2 != currhead);
637 
638  if( bottleneckRuleOut(scip, graph, treecosts, basebottlenecks,
639  nodepos, treeedges, currcost, extcost, curredge, flipedge(e2), dfsdepth, allow_eq, eqstack, eqstack_size, eqmark) )
640  {
641  return TRUE;
642  }
643  }
644  }
645  }
646  }
647  }
648 
649  return FALSE;
650 }
651 
652 /** updates node bounds for reduced cost fixings */
653 static
655  SCIP_Real* fixingbounds, /**< fixing bounds */
656  const GRAPH* graph, /**< graph data structure */
657  const SCIP_Real* pathdist, /**< shortest path distances */
658  const PATH* vnoi, /**< Voronoi paths */
659  SCIP_Real lpobjval, /**< LP objective */
660  SCIP_Bool initialize /**< initialize fixing bounds? */
661 )
662 {
663  const int nnodes = graph->knots;
664 
665  assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
666 
667  if( initialize )
668  for( int k = 0; k < nnodes; k++ )
669  fixingbounds[k] = -FARAWAY;
670 
671  for( int k = 0; k < nnodes; k++ )
672  {
673  SCIP_Real fixbnd;
674 
675  if( !graph->mark[k] )
676  continue;
677 
678  assert(!Is_pterm(graph->term[k]));
679 
680  fixbnd = pathdist[k] + vnoi[k].dist + lpobjval;
681 
682  if( fixbnd > fixingbounds[k] )
683  fixingbounds[k] = fixbnd;
684  }
685 }
686 
687 /** updates node bounds for reduced cost replacement */
688 static
690  SCIP* scip, /**< SCIP */
691  SCIP_Real* replacebounds, /**< replacement bounds */
692  const GRAPH* graph, /**< graph data structure */
693  const SCIP_Real* cost, /**< reduced costs */
694  const SCIP_Real* pathdist, /**< shortest path distances */
695  const PATH* vnoi, /**< Voronoi paths */
696  const int* vbase, /**< bases to Voronoi paths */
697  int* nodearrint, /**< for internal computations */
698  SCIP_Real lpobjval, /**< LP objective */
699  SCIP_Real upperbound, /**< upper bound */
700  int root, /**< DA root */
701  SCIP_Bool initialize, /**< initialize fixing bounds? */
702  SCIP_Bool extendedsearch /**< perform extended searching? */
703 )
704 {
705  unsigned int* eqstack = NULL;
706  SCIP_Bool* eqmark = NULL;
707  int outedges[STP_DABD_MAXDEGREE];
708  const int nnodes = graph->knots;
709  const SCIP_Real cutoff = upperbound - lpobjval;
710  const int halfnedges = graph->edges / 2;
711 
712  assert(!SCIPisNegative(scip, cutoff));
713  assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
714 
715  if( initialize )
716  for( int k = 0; k < nnodes; k++ )
717  replacebounds[k] = -FARAWAY;
718 
719  if( extendedsearch )
720  {
721  SCIP_CALL( SCIPallocCleanBufferArray(scip, &eqmark, halfnedges) );
722  SCIP_CALL( SCIPallocBufferArray(scip, &eqstack, halfnedges) );
723  }
724 
725  for( int node = 0; node < nnodes; node++ )
726  {
727  const int degree = graph->grad[node];
728 
729  if( degree >= 3 && !Is_gterm(graph->term[node]) )
730  {
731  SCIP_Real fixbnd;
732 
733  /* bound already good enough? */
734  if( SCIPisLT(scip, upperbound, replacebounds[node]) )
735  continue;
736 
737  fixbnd = pathdist[node] + vnoi[node].dist + vnoi[node + nnodes].dist + lpobjval;
738 
739  assert(!Is_pterm(graph->term[node]));
740 
741  /* Y-test for small degrees */
742  if( degree <= STP_DABD_MAXDEGREE && !SCIPisLT(scip, upperbound, fixbnd) )
743  {
744  int eqstack_size = 0;
745  int edgecount = 0;
746 
747  fixbnd = FARAWAY;
748 
749  for( int e = graph->outbeg[node]; e != EAT_LAST; e = graph->oeat[e] )
750  outedges[edgecount++] = e;
751 
752  assert(edgecount == degree);
753 
754  /* compute cost for each root and save minimum */
755  for( int i = 0; i < degree; i++ )
756  {
757  const int tmproot = graph->head[outedges[i]];
758  const int rootedge = flipedge(outedges[i]);
759 
760  /* loop over all combinations of Y with root tmproot */
761  for( int j = i + 1; j <= i + degree - 2; j++ )
762  {
763  for( int k = j + 1; k <= i + degree - 1; k++ )
764  {
765  const int outedge1 = outedges[j % degree];
766  const int outedge2 = outedges[k % degree];
767  const int leaf1 = graph->head[outedge1];
768  const int leaf2 = graph->head[outedge2];
769 
770  SCIP_Real tmpcostY;
771 
772  assert(leaf1 != leaf2 && tmproot != leaf1 && tmproot != leaf2);
773  assert(vbase[leaf1] >= 0 || vnoi[leaf1].dist >= FARAWAY);
774  assert(vbase[leaf2] >= 0 || vnoi[leaf2].dist >= FARAWAY);
775 
776  if( leaf1 == root || leaf2 == root )
777  continue;
778 
779  tmpcostY = pathdist[tmproot] + cost[rootedge] + cost[outedge1] + cost[outedge2];
780 
781  if( vbase[leaf1] != vbase[leaf2] )
782  {
783  tmpcostY += vnoi[leaf1].dist + vnoi[leaf2].dist;
784  }
785  else
786  {
787  /* also covers the case that leaf is a terminal */
788  const SCIP_Real leaf1far = vnoi[leaf1 + nnodes].dist;
789  const SCIP_Real leaf2far = vnoi[leaf2 + nnodes].dist;
790 
791  assert(vbase[leaf1 + nnodes] >= 0 || leaf1far == FARAWAY);
792  assert(vbase[leaf2 + nnodes] >= 0 || leaf2far == FARAWAY);
793 
794  tmpcostY += MIN(leaf1far + vnoi[leaf2].dist, vnoi[leaf1].dist + leaf2far);
795  }
796 
797  if( tmpcostY < fixbnd )
798  {
799  if( extendedsearch && SCIPisLE(scip, tmpcostY, cutoff) )
800  {
801  int tree3outedges[2];
802  SCIP_Bool ruleout;
803 #ifndef NDEBUG
804  const SCIP_Real tmpcostYorg = tmpcostY;
805 #endif
806  tree3outedges[0] = outedge1;
807  tree3outedges[1] = outedge2;
808 
809  SCIP_CALL( reduce_check3Tree(scip, graph, root, cost, pathdist, vnoi, vbase, cutoff, tree3outedges, rootedge, nodearrint,
810  &tmpcostY, &ruleout, eqstack, &eqstack_size, eqmark) );
811 
812  if( ruleout )
813  tmpcostY = FARAWAY;
814 
815 #ifndef NDEBUG
816  assert(tmpcostY >= tmpcostYorg);
817 #endif
818  }
819 
820  if( tmpcostY < fixbnd )
821  fixbnd = tmpcostY;
822  }
823  }
824  } /* Y loop */
825  } /* root loop */
826 
827  fixbnd += lpobjval;
828 
829  assert(SCIPisGE(scip, fixbnd, pathdist[node] + vnoi[node].dist + vnoi[node + nnodes].dist + lpobjval)
830  || fixbnd >= FARAWAY);
831 
832  if( extendedsearch )
833  {
834  for( int i = 0; i < eqstack_size; i++ )
835  eqmark[eqstack[i]] = FALSE;
836 
837  for( int i = 0; i < halfnedges; i++ )
838  assert(eqmark[i] == FALSE);
839  }
840  }
841 
842  if( fixbnd > replacebounds[node] )
843  replacebounds[node] = fixbnd;
844  }
845  }
846  if( extendedsearch )
847  {
848  assert(eqstack != NULL && eqmark != NULL);
849 
850  SCIPfreeBufferArray(scip, &eqstack);
851  SCIPfreeCleanBufferArray(scip, &eqmark);
852  }
853 
854  return SCIP_OKAY;
855 }
856 
857 /** updates edge fixing bounds for reduced cost fixings */
858 static
860  SCIP_Real* fixingbounds, /**< fixing bounds */
861  const GRAPH* graph, /**< graph data structure */
862  const SCIP_Real* cost, /**< reduced costs */
863  const SCIP_Real* pathdist, /**< shortest path distances */
864  const PATH* vnoi, /**< Voronoi paths */
865  SCIP_Real lpobjval, /**< LP objective */
866  int extnedges, /**< number of edges for extended problem */
867  SCIP_Bool initialize, /**< initialize fixing bounds? */
868  SCIP_Bool undirected /**< only consider undirected edges */
869 )
870 {
871  const int nnodes = graph->knots;
872 
873  assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
874  assert(extnedges > 0);
875 
876  if( initialize )
877  for( int e = 0; e < extnedges; e++ )
878  fixingbounds[e] = -FARAWAY;
879 
880  for( int k = 0; k < nnodes; k++ )
881  {
882  if( !graph->mark[k] )
883  continue;
884 
885  assert(!Is_pterm(graph->term[k]));
886 
887  if( undirected )
888  {
889  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
890  {
891  const int head = graph->head[e];
892 
893  if( graph->mark[head] )
894  {
895  const int erev = flipedge(e);
896  const SCIP_Real fixbnd = pathdist[k] + cost[e] + vnoi[head].dist + lpobjval;
897  const SCIP_Real fixbndrev = pathdist[head] + cost[erev] + vnoi[k].dist + lpobjval;
898  const SCIP_Real minbnd = MIN(fixbnd, fixbndrev);
899 
900  assert(fixingbounds[e] == fixingbounds[erev]);
901 
902  if( minbnd > fixingbounds[e] )
903  {
904  fixingbounds[e] = minbnd;
905  fixingbounds[erev] = minbnd;
906  }
907  }
908  }
909  }
910  else
911  {
912  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
913  {
914  const int head = graph->head[e];
915 
916  if( graph->mark[head] )
917  {
918  const SCIP_Real fixbnd = pathdist[k] + cost[e] + vnoi[head].dist + lpobjval;
919 
920  if( fixbnd > fixingbounds[e] )
921  fixingbounds[e] = fixbnd;
922  }
923  }
924  }
925  }
926 }
927 
928 /** eliminate nodes by using fixing-bounds and reduced costs */
929 static
931  SCIP* scip, /**< SCIP data structure */
932  GRAPH* graph, /**< graph data structure */
933  GRAPH* transgraph, /**< graph data structure or NULL */
934  const SCIP_Real* fixingbounds, /**< fixing bounds */
935  SCIP_Real upperbound /**< best upperbound */
936 )
937 {
938  int nfixed = 0;
939  int nnodes = graph->knots;
940 
941  assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
942 
943  for( int k = 0; k < nnodes; k++ )
944  {
945  if( !graph->mark[k] || Is_term(graph->term[k]) )
946  continue;
947 
948  assert(!Is_pterm(graph->term[k]));
949 
950  if( SCIPisLT(scip, upperbound, fixingbounds[k]) )
951  {
952  SCIPdebugMessage("delete knot %d %f < %f %d\n", k, upperbound, fixingbounds[k], graph->grad[k]);
953 
954  graph_knot_del(scip, graph, k, TRUE);
955  if( transgraph != NULL )
956  graph_knot_del(scip, transgraph, k, FALSE);
957  }
958  }
959  return nfixed;
960 }
961 
962 static
964  SCIP* scip, /**< SCIP data structure */
965  GRAPH* graph, /**< graph data structure */
966  const PATH* vnoi, /**< Voronoi data structure */
967  const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
968  const SCIP_Real* cost, /**< dual ascent costs */
969  const SCIP_Real* replacebounds, /**< replacement bounds */
970  int* nodetouched, /**< node array for internal computations */
971  SCIP_Real lpobjval, /**< lower bound corresponding to pathdist and vnoi */
972  SCIP_Real upperbound /**< best upperbound */
973 )
974 {
975  int adjvert[STP_DABD_MAXDEGREE];
977  SCIP_Real cutoffsrev[STP_DABD_MAXDNEDGES];
978  int nfixed = 0;
979  const int nnodes = graph->knots;
980 
981  for( int k = 0; k < nnodes; k++ )
982  nodetouched[k] = 0;
983 
984  /* main loop */
985  for( int degree = 3; degree <= STP_DABD_MAXDEGREE; degree++ )
986  {
987  for( int k = 0; k < nnodes; k++ )
988  {
989  if( degree != graph->grad[k] || Is_gterm(graph->term[k]) )
990  continue;
991 
992  if( SCIPisLT(scip, upperbound, replacebounds[k]) && nodetouched[k] == 0 )
993  {
994  int edgecount = 0;
995  SCIP_Bool success;
996 
997  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
998  nodetouched[graph->head[e]]++;
999 
1000  /* fill cutoff */
1001 
1002  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
1003  adjvert[edgecount++] = graph->head[e];
1004 
1005  assert(edgecount == degree);
1006 
1007  edgecount = 0;
1008  for( int i = 0; i < degree - 1; i++ )
1009  {
1010  const int vert = adjvert[i];
1011  for( int i2 = i + 1; i2 < degree; i2++ )
1012  {
1013  const int vert2 = adjvert[i2];
1014 
1015  assert(edgecount < STP_DABD_MAXDNEDGES);
1016 
1017  cutoffs[edgecount] = upperbound - lpobjval - (pathdist[vert] + vnoi[vert2].dist);
1018  cutoffsrev[edgecount] = upperbound - lpobjval - (pathdist[vert2] + vnoi[vert].dist);
1019 
1020  edgecount++;
1021  }
1022  }
1023 
1024  assert(edgecount > 0);
1025 
1026  /* try to eliminate */
1027 
1028  SCIP_CALL_ABORT(graph_knot_delPseudo(scip, graph, cost, cutoffs, cutoffsrev, k, &success));
1029 
1030  if( success )
1031  nfixed++;
1032  else
1033  {
1034  assert(graph->grad[k] == degree);
1035 
1036  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
1037  {
1038  nodetouched[graph->head[e]]--;
1039  assert(nodetouched[graph->head[e]] >= 0);
1040  }
1041  }
1042  }
1043  }
1044  }
1045 
1046  return nfixed;
1047 }
1048 /** eliminate edges by using fixing-bounds and reduced costs */
1049 static
1051  SCIP* scip, /**< SCIP data structure */
1052  GRAPH* graph, /**< graph data structure */
1053  GRAPH* transgraph, /**< graph data structure or NULL */
1054  const SCIP_Real* fixingbounds, /**< fixing bounds */
1055  const int* result, /**< solution */
1056  SCIP_Real upperbound /**< best upperbound */
1057 )
1058 {
1059  int nfixed = 0;
1060  int nnodes = graph->knots;
1061  const SCIP_Bool solgiven = (result != NULL);
1062 
1063  assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
1064  assert(!solgiven || graph_sol_valid(scip, graph, result));
1065 
1066  for( int k = 0; k < nnodes; k++ )
1067  {
1068  int e;
1069  if( !graph->mark[k] )
1070  continue;
1071 
1072  assert(!Is_pterm(graph->term[k]));
1073 
1074  e = graph->outbeg[k];
1075  while( e != EAT_LAST )
1076  {
1077  const int enext = graph->oeat[e];
1078 
1079  if( graph->mark[graph->head[e]] )
1080  {
1081  const int erev = flipedge(e);
1082  SCIP_Bool delete;
1083 
1084  if( !solgiven || result[e] == CONNECT || result[erev] == CONNECT )
1085  delete = (SCIPisLT(scip, upperbound, fixingbounds[e]) && SCIPisLT(scip, upperbound, fixingbounds[erev]));
1086  else
1087  delete = (upperbound <= fixingbounds[e] && upperbound <= fixingbounds[erev]);
1088 
1089  if( delete )
1090  {
1091  SCIPdebugMessage("delete edge %d \n", e);
1092 
1093  graph_edge_del(scip, graph, e, TRUE);
1094  if( transgraph != NULL )
1095  graph_edge_del(scip, transgraph, e, FALSE);
1096  nfixed++;
1097  }
1098  }
1099 
1100  e = enext;
1101  }
1102  }
1103 
1104  return nfixed;
1105 }
1106 
1107 
1108 /** find roots for PC and MW during DA reduction */
1109 static
1111  SCIP* scip, /**< SCIP data structure */
1112  GRAPH* graph, /**< graph data structure */
1113  const GRAPH* transgraph, /**< graph data structure */
1114  const SCIP_Real* cost, /**< da reduced cost */
1115  const SCIP_Real* bestcost, /**< best incumbent da reduced cost */
1116  SCIP_Real lpobjval, /**< da lower bound */
1117  SCIP_Real bestlpobjval, /**< best da lower bound */
1118  SCIP_Real upperbound, /**< da upper bound */
1119  SCIP_Real prizesum, /**< sum of positive prizes */
1120  SCIP_Bool rerun, /**< not the first run? */
1121  SCIP_Bool probrooted, /**< is transgraph a rooted RMW or RPC? */
1122  PATH* vnoi, /**< SP array */
1123  int* state, /**< array */
1124  int* pathedge, /**< array */
1125  int* vbase, /**< array */
1126  STP_Bool* nodearrchar, /**< bool array */
1127  int* roots, /**< root array */
1128  int* rootscount /**< number of roots */
1129 )
1130 {
1131  const int root = graph->source;
1132  const int nnodes = graph->knots;
1133  const SCIP_Bool graphextended = graph->extended;
1134  int nroots = *rootscount;
1135 
1136  assert(transgraph->extended);
1137 
1138  if( graphextended )
1139  graph_pc_2org(graph);
1140 
1141  assert(rerun || nroots == 0);
1142 
1143  /*
1144  * get possible roots
1145  */
1146 
1147  BMSclearMemoryArray(nodearrchar, nnodes);
1148 
1149  if( rerun )
1150  for( int i = 0; i < nroots; i++ )
1151  nodearrchar[roots[i]] = TRUE;
1152 
1153  if( probrooted )
1154  {
1155  const int transroot = transgraph->source;
1156 
1157  assert(transgraph->term2edge != NULL);
1158 
1159  for( int e = transgraph->outbeg[transroot]; e != EAT_LAST; e = transgraph->oeat[e] )
1160  {
1161  const int pseudoterm = transgraph->head[e];
1162 
1163  if( Is_term(transgraph->term[pseudoterm]) && transgraph->term2edge[pseudoterm] >= 0 )
1164  {
1165  int realterm;
1166  assert(transgraph->tail[transgraph->term2edge[pseudoterm]] == pseudoterm);
1167 
1168  realterm = transgraph->head[transgraph->term2edge[pseudoterm]];
1169  assert(Is_pterm(transgraph->term[realterm]));
1170 
1171  if( rerun && nodearrchar[realterm] )
1172  continue;
1173 
1174  if( SCIPisGT(scip, cost[e], upperbound - lpobjval) || SCIPisGT(scip, bestcost[e], upperbound - bestlpobjval) )
1175  {
1176  assert(SCIPisPositive(scip, graph->prize[realterm]));
1177  assert(nroots < graph->terms);
1178 
1179  roots[nroots++] = realterm;
1180  nodearrchar[realterm] = TRUE;
1181  }
1182  }
1183  }
1184  }
1185  else
1186  {
1187  for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1188  {
1189  const int pseudoterm = graph->head[e];
1190 
1191  if( Is_pterm(graph->term[pseudoterm]) )
1192  {
1193  int realterm = -1;
1194  int e3 = -1;
1195 
1196  assert(graph->grad[pseudoterm] == 2);
1197 
1198  for( int e2 = transgraph->inpbeg[pseudoterm]; e2 != EAT_LAST; e2 = transgraph->ieat[e2] )
1199  {
1200  if( transgraph->cost[e2] == 0.0 )
1201  realterm = transgraph->tail[e2];
1202  else
1203  e3 = e2;
1204  }
1205 
1206  if( rerun && nodearrchar[realterm] )
1207  continue;
1208 
1209  assert(realterm >= 0);
1210  assert(e3 >= 0);
1211  assert(realterm != root);
1212  assert(SCIPisEQ(scip, graph->cost[e], graph->cost[e3]));
1213  assert(graph->mark[realterm]);
1214 
1215  if( SCIPisGT(scip, cost[e3], upperbound - lpobjval) || SCIPisGT(scip, bestcost[e3], upperbound - bestlpobjval) )
1216  {
1217  assert(SCIPisPositive(scip, graph->prize[realterm]));
1218  assert(nroots < graph->terms);
1219 
1220  roots[nroots++] = realterm;
1221  nodearrchar[realterm] = TRUE;
1222  }
1223  }
1224  }
1225  }
1226 
1227  /* could more roots be found? */
1228  if( nroots > *rootscount && graph->terms > 2 )
1229  {
1230  /*
1231  * try to find additional roots by shortest path computations
1232  */
1233 
1234  SCIP_Real maxprize = 0.0;
1235  int qstart = *rootscount;
1236 
1237  for( int i = 0; i < nnodes; i++ )
1238  {
1239  state[i] = UNKNOWN;
1240  vnoi[i].dist = FARAWAY;
1241  vnoi[i].edge = UNKNOWN;
1242 
1243  if( Is_term(graph->term[i]) && !nodearrchar[i] && graph->mark[i] && maxprize < graph->prize[i] && graph->prize[i] < prizesum )
1244  maxprize = graph->prize[i];
1245  }
1246 
1247  for( int rounds = 0; rounds < 10 && qstart != nroots; rounds++ )
1248  {
1249  const int oldnroots = nroots;
1250 
1251  SCIPdebugMessage("root-finding round %d \n", rounds);
1252 
1253  for( int i = qstart; i < nroots; i++ )
1254  {
1255  int nlabeled;
1256  const int term = roots[i];
1257  assert(nodearrchar[term]);
1258  assert(Is_term(graph->term[term]));
1259  assert(SCIPisPositive(scip, graph->prize[term]));
1260 
1261  /* look for terminals in neighborhood of term */
1262  graph_sdPaths(scip, graph, vnoi, graph->cost, maxprize, pathedge, state, vbase, &nlabeled, term, term, 200);
1263 
1264  /* check labelled nodes */
1265  for( int k = 0; k < nlabeled; k++ )
1266  {
1267  const int l = vbase[k];
1268 
1269  assert(graph->mark[l]);
1270  assert(state[l] != UNKNOWN);
1271 
1272  if( Is_term(graph->term[l]) && !nodearrchar[l] && vnoi[l].dist <= graph->prize[l] )
1273  {
1274  roots[nroots++] = l;
1275  nodearrchar[l] = TRUE;
1276 
1277  SCIPdebugMessage("added new root %d dist: %f, prize: %f \n", l, vnoi[l].dist, graph->prize[l]);
1278  }
1279 
1280  /* restore data structures */
1281  state[l] = UNKNOWN;
1282  vnoi[l].dist = FARAWAY;
1283  vnoi[l].edge = UNKNOWN;
1284  }
1285  }
1286 #ifndef NDEBUG
1287  for( int k = 0; k < nnodes; k++ )
1288  {
1289  assert(state[k] == UNKNOWN);
1290  assert(vnoi[k].dist == FARAWAY);
1291  assert(vnoi[k].edge == UNKNOWN);
1292  }
1293 #endif
1294 
1295  qstart = oldnroots;
1296  }
1297  SCIPdebugMessage("number of terminals found by DA: %d \n", nroots);
1298  }
1299 
1300  if( rerun )
1301  SCIPdebugMessage("new roots in rerun %d \n", nroots - *rootscount);
1302 
1303  *rootscount = nroots;
1304 
1305  if( graphextended )
1306  graph_pc_2trans(graph);
1307 }
1308 
1309 
1310 /** set prize of marked terminals to blocked */
1311 static
1313  SCIP* scip, /**< SCIP data structure */
1314  const int* roots, /**< root array */
1315  int nrootsold, /**< old number of roots */
1316  int nrootsnew, /**< new number of roots */
1317  SCIP_Real prizesum, /**< sum of positive prizes */
1318  GRAPH* graph, /**< graph data structure */
1319  SCIP_Bool* userec, /**< recombination? */
1320  STPSOLPOOL** solpool /**< solution pool */
1321 )
1322 {
1323  const int root = graph->source;
1324  const SCIP_Bool graphextended = graph->extended;
1325 
1326  if( graphextended )
1327  graph_pc_2org(graph);
1328 
1329  if( *userec && *solpool != NULL )
1330  {
1331  *userec = FALSE;
1332  SCIPStpHeurRecFreePool(scip, solpool);
1333  *solpool = NULL;
1334  }
1335 
1336  assert(graph != NULL);
1337  assert(roots != NULL);
1338  assert(!graph->extended);
1339  assert(nrootsnew >= 0 && nrootsold >= 0);
1340 
1341  for( int i = nrootsold; i < nrootsnew; i++ )
1342  {
1343  int e;
1344  const int term = roots[i];
1345  const int pterm = graph->head[graph->term2edge[term]];
1346 
1347  assert(Is_term(graph->term[term]));
1348  assert(SCIPisPositive(scip, graph->prize[term]));
1349  assert(pterm >= 0);
1350  assert(Is_pterm(graph->term[pterm]));
1351 
1352  for( e = graph->inpbeg[pterm]; e != EAT_LAST; e = graph->ieat[e] )
1353  if( root == graph->tail[e] )
1354  break;
1355 
1356  assert(e != EAT_LAST);
1357  assert(SCIPisEQ(scip, graph->prize[term], graph->cost[e]));
1358 
1359  graph->prize[term] = prizesum;
1360  graph->cost[e] = prizesum;
1361  }
1362 
1363  if( graphextended )
1364  graph_pc_2trans(graph);
1365 }
1366 
1367 /** are the dual costs still valid, i.e. are there zero cost paths from the root to all terminals? */
1368 static
1370  SCIP* scip, /**< SCIP data structure */
1371  GRAPH* transgraph, /**< graph data structure */
1372  const SCIP_Real* cost, /**< dual ascent costs */
1373  int* nodearrint, /**< int array */
1374  STP_Bool* nodearrbool /**< bool array */
1375 )
1376 {
1377  int* const queue = nodearrint;
1378  STP_Bool* visited = nodearrbool;
1379  int qsize;
1380  const int root = transgraph->source;
1381  const int nnodes = transgraph->knots;
1382 
1383  /*
1384  * construct new graph corresponding to zero cost paths from the root to all terminals
1385  */
1386 
1387  BMSclearMemoryArray(visited, nnodes);
1388 
1389  qsize = 0;
1390  visited[root] = TRUE;
1391  queue[qsize++] = root;
1392 
1393  /* DFS */
1394  while( qsize )
1395  {
1396  const int k = queue[--qsize];
1397 
1398  /* traverse outgoing arcs */
1399  for( int a = transgraph->outbeg[k]; a != EAT_LAST; a = transgraph->oeat[a] )
1400  {
1401  const int head = transgraph->head[a];
1402 
1403  if( !visited[head] && SCIPisZero(scip, cost[a]) )
1404  {
1405  visited[head] = TRUE;
1406  queue[qsize++] = head;
1407  }
1408  }
1409  assert(qsize <= nnodes);
1410  }
1411 
1412  for( int k = 0; k < nnodes; k++ )
1413  if( Is_term(transgraph->term[k]) && !visited[k] )
1414  {
1415  return FALSE;
1416  }
1417 
1418  return TRUE;
1419 }
1420 
1421 
1422 /** pertubate edge costs for PCMW dual-ascent */
1423 static
1425  SCIP* scip,
1426  const GRAPH* graph,
1427  GRAPH* transgraph,
1428  const int* result,
1429  STP_Bool* nodearrchar,
1430  int randomize
1431 )
1432 {
1433  int e;
1434  const int root = graph->source;
1435  const int newroot = transgraph->source;
1436  const int nnodes = graph->knots;
1437  const int nedges = graph->edges;
1438 
1439  BMSclearMemoryArray(nodearrchar, nnodes);
1440 
1441  /* mark all vertices visited in regular graph */
1442  for( e = 0; e < nedges; e++ )
1443  if( result[e] == CONNECT && graph->tail[e] != root )
1444  nodearrchar[graph->head[e]] = TRUE;
1445  srand((unsigned)graph->terms);
1446 
1447  if( graph->stp_type != STP_MWCSP )
1448  {
1450 
1451  for( int k = 0; k < nnodes; k++ )
1452  {
1453  assert(Is_gterm(graph->term[k]) == Is_gterm(transgraph->term[k]) || transgraph->grad[k] == 0);
1454 
1455  if( randomize > 8 )
1456  pratio = ((SCIP_Real)(rand() % 10)) / (100.0) - 1.0 / 100.0;
1457  else if( randomize > 6 )
1458  pratio = ((SCIP_Real)(rand() % 10)) / (200.0);
1459  else if( randomize > 4 )
1460  pratio = ((SCIP_Real)(rand() % 10)) / (300.0);
1461  else if( randomize > 0 )
1462  pratio = ((SCIP_Real)(rand() % 10)) / 1000.0;
1463  else
1464  pratio = PERTUBATION_RATIO_PC + ((SCIP_Real)(rand() % 10)) / 1000.0;
1465 
1466  assert(SCIPisPositive(scip, 1.0 - pratio));
1467  assert(SCIPisPositive(scip, 1.0 + pratio));
1468 
1469  if( !Is_gterm(graph->term[k]) )
1470  {
1471  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1472  {
1473  assert(transgraph->tail[e] != root);
1474 
1475  if( result[e] == CONNECT || result[flipedge(e)] == CONNECT )
1476  transgraph->cost[e] *= 1.0 - pratio;
1477  else
1478  transgraph->cost[e] *= 1.0 + pratio;
1479  }
1480  }
1481  else if( Is_term(transgraph->term[k]) && k != root && k != newroot )
1482  {
1483  assert(transgraph->grad[k] == 2);
1484 
1485  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1486  if( SCIPisPositive(scip, transgraph->cost[e]) )
1487  {
1488  assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1489  assert(transgraph->tail[e] != root);
1490  assert(result[flipedge(e)] != CONNECT);
1491 
1492  if( result[e] == CONNECT )
1493  transgraph->cost[e] *= 1.0 - pratio;
1494  else
1495  transgraph->cost[e] *= 1.0 + pratio;
1496 
1497  assert(SCIPisPositive(scip, transgraph->cost[e]));
1498  }
1499  }
1500  }
1501 
1502  return;
1503  }
1504 
1505  for( int k = 0; k < nnodes; k++ )
1506  {
1507  SCIP_Real pratio = PERTUBATION_RATIO;
1508 
1509  assert(Is_gterm(graph->term[k]) == Is_gterm(transgraph->term[k]));
1510 
1511  if( randomize > 8 )
1512  pratio = ((SCIP_Real)(rand() % 10)) / (50.0) - 1.0 / 10.0;
1513  else if( randomize > 6 )
1514  pratio = ((SCIP_Real)(rand() % 10)) / (20.0);
1515  else if( randomize > 4 )
1516  pratio = ((SCIP_Real)(rand() % 10)) / (30.0);
1517  else if( randomize > 0 )
1518  pratio = ((SCIP_Real)(rand() % 10)) / 100.0;
1519  else
1520  pratio = PERTUBATION_RATIO + ((SCIP_Real)(rand() % 10)) / 200.0;
1521 
1522  if( !Is_gterm(graph->term[k]) )
1523  {
1524  if( nodearrchar[k] )
1525  {
1526  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1527  transgraph->cost[e] *= 1.0 - pratio;
1528  }
1529  else
1530  {
1531  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1532  transgraph->cost[e] *= 1.0 + pratio;
1533  }
1534  }
1535  else if( Is_term(transgraph->term[k]) && k != root && k != newroot )
1536  {
1537  assert(transgraph->grad[k] == 2);
1538 
1539  if( nodearrchar[k] )
1540  {
1541  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1542  if( SCIPisPositive(scip, transgraph->cost[e]) )
1543  {
1544  assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1545 
1546  transgraph->cost[e] *= 1.0 + pratio;
1547  }
1548  }
1549  else
1550  {
1551  for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1552  if( SCIPisPositive(scip, transgraph->cost[e]) )
1553  {
1554  assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1555  transgraph->cost[e] *= 1.0 - pratio;
1556  }
1557  }
1558  }
1559  }
1560 }
1561 
1562 /** order roots */
1563 static
1565  SCIP* scip, /**< SCIP data structure */
1566  const GRAPH* graph, /**< graph structure */
1567  int* terms, /**< sol int array corresponding to upper bound */
1568  int nterms, /**< number of terminals */
1569  SCIP_Bool randomize, /**< randomize */
1570  SCIP_RANDNUMGEN* randnumgen /**< random number generator */
1571 )
1572 {
1573  int* termdegs;
1574  int maxdeg = 0;
1575 
1576  assert(terms != NULL);
1577  assert(nterms > 0);
1578 
1579  SCIP_CALL( SCIPallocBufferArray(scip, &termdegs, nterms) );
1580 
1581  for( int i = 0; i < nterms; i++ )
1582  {
1583  const int grad = graph->grad[terms[i]];
1584  assert(terms[i] >= 0);
1585 
1586  termdegs[i] = -grad;
1587 
1588  if( grad > maxdeg )
1589  maxdeg = termdegs[i];
1590  }
1591 
1592  if( randomize )
1593  for( int i = 0; i < nterms; i++ )
1594  termdegs[i] -= SCIPrandomGetInt(randnumgen, 0, maxdeg);
1595 
1596  SCIPsortIntInt(termdegs, terms, nterms);
1597 
1598  SCIPfreeBufferArray(scip, &termdegs);
1599 
1600  return SCIP_OKAY;
1601 }
1602 
1603 
1604 /** compute primal solution during dual-ascent routine for PCSTP or MWCSP */
1605 static
1607  SCIP* scip, /**< SCIP data structure */
1608  GRAPH* graph, /**< graph data structure */
1609  STPSOLPOOL* pool, /**< solution pool */
1610  PATH* vnoi, /**< Voronoi data structure */
1611  const SCIP_Real* cost, /**< dual ascent costs */
1612  SCIP_Real* pathdist, /**< distance array from shortest path calculations */
1613  SCIP_Real* upperbound, /**< upperbound pointer */
1614  int* result1, /**< sol int array corresponding to upper bound */
1615  int* result2, /**< sol int array */
1616  int* vbase, /**< int array */
1617  int* pathedge, /**< int array */
1618  STP_Bool* nodearrchar, /**< node array storing solution vertices */
1619  SCIP_Bool* apsol /**< ascend-prune sol? */
1620 )
1621 {
1622  SCIP_Real ub;
1623  const int nedges = graph->edges;
1624  SCIP_Bool success;
1625 
1626  /* compute new solution and store it in result2 */
1627 
1628  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, result2, vbase, -1, nodearrchar, &success, FALSE) );
1629  assert(success);
1630 
1631  SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1632 
1633  assert(graph_sol_valid(scip, graph, result2));
1634 
1635  /* compute objective value */
1636  ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1637  SCIPdebugMessage("DA: first new sol value in computeDaSolPcMw: %f ... old value: %f \n", ub, *upperbound);
1638 
1639  /* try recombination? */
1640  if( pool != NULL )
1641  {
1642  SCIPdebugMessage("ub %f vs best sol %f\n", ub, pool->sols[0]->obj);
1643  SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
1644 
1645 #ifdef SCIP_DEBUG
1646  for( int i = 0; i < pool->size; i++ )
1647  {
1648  printf(" %f ", pool->sols[i]->obj);
1649  }
1650  printf("\n ");
1651 #endif
1652 
1653  if( success && pool->size >= 2 )
1654  {
1655  /* get index of just added solution */
1656  int solindex = pool->maxindex;
1657 
1658  SCIP_Bool solfound;
1659 
1660  SCIPdebugMessage("POOLSIZE %d \n", pool->size);
1661 
1662  SCIP_CALL( SCIPStpHeurRecRun(scip, pool, NULL, NULL, graph, NULL, &solindex, 3, pool->size, FALSE, &solfound) );
1663 
1664  if( solfound )
1665  {
1666  STPSOL* sol = SCIPStpHeurRecSolfromIdx(pool, solindex);
1667 
1668  assert(pool->size >= 2);
1669  assert(sol != NULL);
1670 
1671  SCIPdebugMessage("DA: rec found better solution with obj %f vs %f \n", sol->obj, ub);
1672 
1673  if( SCIPisLE(scip, sol->obj, ub) )
1674  {
1675  BMScopyMemoryArray(result2, sol->soledges, nedges);
1676 
1677  SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1678  ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1679 
1680  if( SCIPisLT(scip, ub, sol->obj) )
1681  SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
1682  }
1683  }
1684  }
1685  }
1686 
1687  if( SCIPisLE(scip, ub, *upperbound) )
1688  {
1689  SCIPdebugMessage("DA: improved incumbent %f vs %f, return \n", ub, *upperbound);
1690 
1691  *apsol = TRUE;
1692  *upperbound = ub;
1693  BMScopyMemoryArray(result1, result2, nedges);
1694  }
1695 
1696  if( graph->stp_type != STP_MWCSP || !(*apsol) )
1697  return SCIP_OKAY;
1698 
1699 #if 1
1700  SCIP_CALL( SCIPStpHeurRecExclude(scip, graph, result1, result2, pathedge, nodearrchar, &success) );
1701 
1702  if( success )
1703  {
1704  BMScopyMemoryArray(result1, result2, nedges);
1705  *upperbound = graph_sol_getObj(graph->cost, result1, 0.0, nedges);
1706  SCIPdebugMessage("DA: afterLastExclusion %f \n", *upperbound);
1707  }
1708 #endif
1709 
1710  return SCIP_OKAY;
1711 }
1712 
1713 
1714 /** try to improve both dual and primal solution */
1715 static
1717  SCIP* scip,
1718  GRAPH* graph,
1719  GRAPH* transgraph,
1720  STPSOLPOOL* pool,
1721  PATH* vnoi,
1722  GNODE** gnodearr,
1723  SCIP_Real* cost,
1724  SCIP_Real* costrev,
1725  SCIP_Real* bestcost,
1726  SCIP_Real* pathdist,
1727  int* state,
1728  int* vbase,
1729  int* pathedge,
1730  int* result,
1731  int* result2,
1732  int* transresult,
1733  STP_Bool* nodearrchar,
1734  SCIP_Real* upperbound,
1735  SCIP_Real* lpobjval,
1736  SCIP_Real* bestlpobjval,
1737  SCIP_Real* minpathcost,
1738  SCIP_Bool* apsol,
1739  SCIP_Real offset,
1740  int extnedges,
1741  int pertubation
1742 )
1743 {
1744  SCIP_Real lb;
1745  int e;
1746  const int root = graph->source;
1747  const int nedges = graph->edges;
1748  const int transnedges = transgraph->edges;
1749 
1750  graph_pc_2transcheck(graph);
1751 
1752  /* pertubate the reduced cost? */
1753  if( graph->stp_type == STP_MWCSP )
1754  {
1755  SCIP_Real* transcost;
1756 
1757  SCIP_CALL( SCIPallocBufferArray(scip, &transcost, transnedges) );
1758  BMScopyMemoryArray(transcost, transgraph->cost, transnedges);
1759 
1760  /* result contains no valid solution?*/
1761  if( !(*apsol) )
1762  {
1763  /* compute new solution */
1764  SCIP_Real bound;
1765  SCIP_Bool success;
1766  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, bestcost, result2, vbase, -1, nodearrchar, &success, FALSE) );
1767  assert(success);
1768 
1769  SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1770 
1771  assert(graph_sol_valid(scip, graph, result2));
1772 
1773  bound = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1774 
1775  if( SCIPisLE(scip, bound, *upperbound) )
1776  {
1777  *upperbound = bound;
1778  *apsol = TRUE;
1779  BMScopyMemoryArray(result, result2, nedges);
1780  }
1781  pertubateEdgeCosts(scip, graph, transgraph, result2, nodearrchar, pertubation);
1782  }
1783  else
1784  {
1785  pertubateEdgeCosts(scip, graph, transgraph, result, nodearrchar, pertubation);
1786  }
1787 
1788  /* todo use result as guiding solution? */
1789  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lb, FALSE, FALSE, gnodearr, NULL, transresult, state, root, TRUE, -1.0, nodearrchar) );
1790 
1791  BMScopyMemoryArray(transgraph->cost, transcost, transnedges);
1792 
1793  SCIPfreeBufferArray(scip, &transcost);
1794  }
1795 
1796  SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, upperbound, result, result2, vbase, pathedge, nodearrchar, apsol) );
1797 
1798  /* does result not contain a valid solution? */
1799  if( !(*apsol) )
1800  BMScopyMemoryArray(result, result2, nedges);
1801 
1802  SCIPdebugMessage("DA: pertubated sol value %f \n", *upperbound);
1803 
1804  /* compute guiding solution */
1805 
1806  BMScopyMemoryArray(transresult, result, nedges);
1807 
1808  for( e = nedges; e < transnedges; e++ )
1809  transresult[e] = UNKNOWN;
1810 
1811  /* non-trivial solution */
1812  for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1813  if( Is_term(graph->term[graph->head[e]]) && result[e] != CONNECT )
1814  break;
1815 
1816  if( e != EAT_LAST)
1817  {
1818  int k;
1819  for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1820  if( !Is_term(graph->term[graph->head[e]]) && result[e] == CONNECT )
1821  break;
1822 
1823  assert(e != EAT_LAST);
1824  k = graph->head[e];
1825 
1826  for( e = transgraph->outbeg[k]; e != EAT_LAST; e = transgraph->oeat[e] )
1827  {
1828  if( transgraph->head[e] == graph->knots )
1829  {
1830  transresult[e] = CONNECT;
1831  break;
1832  }
1833  }
1834  }
1835 
1836  assert(!(*apsol) || graph_sol_valid(scip, graph, result));
1837  assert(graph_valid(transgraph));
1838  assert(root == transgraph->source);
1839 
1840  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lb, FALSE, FALSE, gnodearr, transresult, NULL, state, root, TRUE, -1.0, nodearrchar) );
1841 
1842  lb += offset;
1843  *lpobjval = lb;
1844 
1845  SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, upperbound, result, result2, vbase, pathedge, nodearrchar, apsol) );
1846 
1847  assert(!(*apsol) || graph_sol_valid(scip, graph, result));
1848 
1849  if( SCIPisGE(scip, lb, *bestlpobjval) )
1850  {
1851  *bestlpobjval = lb;
1852  BMScopyMemoryArray(bestcost, cost, extnedges);
1853  }
1854 
1855  /* the required reduced path cost to be surpassed */
1856  *minpathcost = *upperbound - *lpobjval;
1857 
1858  return SCIP_OKAY;
1859 }
1860 
1861 /** compute Voronoi region for dual-ascent elimination for PC/MW */
1862 static
1864  SCIP* scip,
1865  GRAPH* transgraph,
1866  PATH* vnoi,
1867  const SCIP_Real* cost,
1868  SCIP_Real* costrev,
1869  SCIP_Real* pathdist,
1870  int* vbase,
1871  int* pathedge
1872 )
1873 {
1874  const int root = transgraph->source;
1875  const int transnnodes = transgraph->knots;
1876  const int transnedges = transgraph->edges;
1877 
1878  for( int k = 0; k < transnnodes; k++ )
1879  transgraph->mark[k] = (transgraph->grad[k] > 0);
1880 
1881  /* distance from root to all nodes */
1882  graph_path_execX(scip, transgraph, root, cost, pathdist, pathedge);
1883 
1884  for( int k = 0; k < transnnodes; k++ )
1885  if( Is_term(transgraph->term[k]) )
1886  assert(SCIPisEQ(scip, pathdist[k], 0.0));
1887 
1888  for( int e = 0; e < transnedges; e++ )
1889  costrev[e] = cost[flipedge(e)];
1890 
1891  /* no paths should go back to the root */
1892  for( int e = transgraph->outbeg[root]; e != EAT_LAST; e = transgraph->oeat[e] )
1893  costrev[e] = FARAWAY;
1894 
1895  /* build Voronoi diagram wrt incoming arcs */
1896  graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, transgraph->path_state);
1897 }
1898 
1899 
1900 /** reduce SPG graph based on information from dual ascent and given upper bound */
1901 static
1903  SCIP* scip, /**< SCIP data structure */
1904  GRAPH* graph, /**< graph data structure */
1905  SCIP_Real* offset, /**< offset pointer */
1906  STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
1907  STP_Bool* nodearrchar, /**< node array storing solution vertices */
1908  const PATH* vnoi, /**< Voronoi data structure */
1909  const SCIP_Real* cost, /**< dual ascent costs */
1910  const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
1911  const int* result, /**< sol int array */
1912  SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
1913  int root, /**< the root */
1914  SCIP_Bool solgiven /**< is sol given? */
1915 )
1916 {
1917  int nfixed = 0;
1918  const int nedges = graph->edges;
1919  const int nnodes = graph->knots;
1920  const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
1921  const SCIP_Bool keepsol = (solgiven && SCIPisZero(scip, minpathcost));
1922 
1923  if( solgiven )
1924  {
1925  for( int k = 0; k < nnodes; k++ )
1926  nodearrchar[k] = FALSE;
1927  for( int e = 0; e < nedges; e++ )
1928  {
1929  if( result[e] == CONNECT )
1930  {
1931  nodearrchar[graph->tail[e]] = TRUE;
1932  nodearrchar[graph->head[e]] = TRUE;
1933  }
1934  }
1935  }
1936 
1937  /* RPC? If yes, restore original graph */
1938  if( rpc )
1939  graph_pc_2org(graph);
1940 
1941  for( int k = 0; k < nnodes; k++ )
1942  {
1943  SCIP_Real redcost;
1944  int e;
1945 
1946  if( !graph->mark[k] )
1947  continue;
1948 
1949  assert(!Is_pterm(graph->term[k]));
1950 
1951  redcost = pathdist[k] + vnoi[k].dist;
1952 
1953  if( rpc && Is_term(graph->term[k]) && SCIPisGT(scip, redcost, minpathcost) && k != root )
1954  {
1955  (*offset) += graph->prize[k];
1956  nfixed += graph_pc_deleteTerm(scip, graph, k);
1957  continue;
1958  }
1959 
1960  if( !Is_term(graph->term[k]) && !keepsol &&
1961  (SCIPisGT(scip, redcost, minpathcost) || (solgiven && SCIPisEQ(scip, redcost, minpathcost) && !nodearrchar[k])) )
1962  {
1963  nfixed += graph->grad[k];
1964  graph_knot_del(scip, graph, k, TRUE);
1965  continue;
1966  }
1967 
1968  e = graph->outbeg[k];
1969  while( e != EAT_LAST )
1970  {
1971  const int head = graph->head[e];
1972  const int enext = graph->oeat[e];
1973 
1974  /* for rpc no artificial terminal arcs should be deleted */
1975  if( (rpc && !graph->mark[head])
1976  || (keepsol && (result[e] == CONNECT || result[flipedge(e)] == CONNECT)) )
1977  {
1978  e = enext;
1979  continue;
1980  }
1981 
1982  redcost = pathdist[k] + cost[e] + vnoi[head].dist;
1983 
1984  if( SCIPisGT(scip, redcost, minpathcost)
1985  || (solgiven && SCIPisEQ(scip, redcost, minpathcost) && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
1986  {
1987  if( marked[flipedge(e)] )
1988  {
1989  graph_edge_del(scip, graph, e, TRUE);
1990  nfixed++;
1991  }
1992  else
1993  {
1994  marked[e] = TRUE;
1995  }
1996  }
1997  e = enext;
1998  }
1999  }
2000 
2001  for( int k = 0; k < nnodes; k++ )
2002  if( graph->grad[k] == 0 && k != root )
2003  graph->mark[k] = FALSE;
2004 
2005  return nfixed;
2006 }
2007 
2008 /** reduce PCSTP or MWCS graph based on information from dual ascent and given upper bound */
2009 static
2011  SCIP* scip, /**< SCIP data structure */
2012  GRAPH* graph, /**< graph data structure */
2013  GRAPH* transgraph, /**< graph data structure */
2014  const PATH* vnoi, /**< Voronoi data structure */
2015  const SCIP_Real* cost, /**< dual ascent costs */
2016  const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2017  SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
2018  const int* result, /**< sol int array */
2019  STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
2020  STP_Bool* nodearrchar, /**< node array storing solution vertices */
2021  SCIP_Bool solgiven /**< is sol given? */
2022 )
2023 {
2024  int nnodes;
2025  int nfixed;
2026  SCIP_Real tmpcost;
2027  SCIP_Bool keepsol = FALSE;
2028 
2029  assert(scip != NULL);
2030  assert(graph != NULL);
2031  assert(pathdist != NULL);
2032  assert(result != NULL);
2033  assert(cost != NULL);
2034  assert(vnoi != NULL);
2035 
2036  assert(SCIPisGE(scip, minpathcost, 0.0));
2037  assert(!solgiven || graph_sol_valid(scip, graph, result));
2038 
2039  if( minpathcost < 0.0 )
2040  minpathcost = 0.0;
2041 
2042  nfixed = 0;
2043  nnodes = graph->knots;
2044 
2045  graph_pc_2orgcheck(graph);
2046 
2047  if( solgiven )
2048  {
2049  const int nedges = graph->edges;
2050 
2051  for( int k = 0; k < nnodes; k++ )
2052  nodearrchar[k] = FALSE;
2053 
2054  for( int e = 0; e < nedges; e++ )
2055  {
2056  if( result[e] == CONNECT )
2057  {
2058  assert(graph->oeat[e] != EAT_FREE);
2059 
2060  nodearrchar[graph->head[e]] = TRUE;
2061  nodearrchar[graph->tail[e]] = TRUE;
2062  }
2063  }
2064  if( SCIPisZero(scip, minpathcost) )
2065  keepsol = TRUE;
2066  }
2067 
2068  /* try to eliminate vertices and edges */
2069  for( int k = 0; k < nnodes; k++ )
2070  {
2071  if( !graph->mark[k] )
2072  continue;
2073 
2074  assert(!Is_pterm(graph->term[k]));
2075 
2076  if( Is_term(graph->term[k]) )
2077  {
2078  int e = graph->outbeg[k];
2079  while( e != EAT_LAST )
2080  {
2081  const int etmp = graph->oeat[e];
2082  tmpcost = pathdist[k] + cost[e] + vnoi[graph->head[e]].dist;
2083 
2084  if( graph->mark[graph->head[e]] &&
2085  ((SCIPisGT(scip, tmpcost, minpathcost) && (!keepsol || (result[e] != CONNECT && result[flipedge(e)] != CONNECT))) ||
2086  (solgiven && tmpcost >= minpathcost && result[e] != CONNECT && result[flipedge(e)] != CONNECT)) )
2087  {
2088  if( marked[flipedge(e)] )
2089  {
2090  assert(!Is_pterm(graph->term[graph->head[e]]));
2091 
2092  graph_edge_del(scip, graph, e, TRUE);
2093  graph_edge_del(scip, transgraph, e, FALSE);
2094  nfixed++;
2095  }
2096  else
2097  {
2098  marked[e] = TRUE;
2099  }
2100  }
2101  e = etmp;
2102  }
2103  continue;
2104  }
2105 
2106  tmpcost = pathdist[k] + vnoi[k].dist;
2107 
2108  if( (SCIPisGT(scip, tmpcost, minpathcost) && !keepsol) ||
2109  (solgiven && tmpcost >= minpathcost && !nodearrchar[k]))
2110  {
2111  while( transgraph->outbeg[k] != EAT_LAST )
2112  {
2113  const int e = transgraph->outbeg[k];
2114 
2115  graph_edge_del(scip, transgraph, e, FALSE);
2116  graph_edge_del(scip, graph, e, TRUE);
2117  nfixed++;
2118  }
2119  assert(graph->outbeg[k] == EAT_LAST);
2120  }
2121  else
2122  {
2123  int e = transgraph->outbeg[k];
2124  while( e != EAT_LAST )
2125  {
2126  const int etmp = transgraph->oeat[e];
2127  tmpcost = pathdist[k] + cost[e] + vnoi[transgraph->head[e]].dist;
2128 
2129  if( (SCIPisGT(scip, tmpcost, minpathcost) && (!keepsol || (result[e] != CONNECT && result[flipedge(e)] != CONNECT))) ||
2130  (solgiven && tmpcost >= minpathcost && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
2131  {
2132  if( marked[flipedge(e)] )
2133  {
2134  graph_edge_del(scip, graph, e, TRUE);
2135  graph_edge_del(scip, transgraph, e, FALSE);
2136 
2137  nfixed++;
2138  }
2139  else
2140  {
2141  marked[e] = TRUE;
2142  }
2143  }
2144  e = etmp;
2145  }
2146  }
2147  }
2148  SCIPdebugMessage("DA: eliminations %d \n", nfixed);
2149 
2150  for( int k = 0; k < nnodes; k++ )
2151  if( graph->grad[k] == 0 && k != graph->source )
2152  graph->mark[k] = FALSE;
2153 
2154  return nfixed;
2155 }
2156 
2157 /** try to run reduction method for best known reduced costs (if they are valid) */
2158 static
2160  SCIP* scip, /**< SCIP data structure */
2161  GRAPH* graph, /**< graph data structure */
2162  GRAPH* transgraph, /**< graph data structure */
2163  PATH* vnoi, /**< Voronoi data structure */
2164  const SCIP_Real* cost, /**< dual ascent costs */
2165  SCIP_Real* costrev, /**< SCIP_Real array */
2166  SCIP_Real* bestcost, /**< best dual ascent costs */
2167  SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2168  SCIP_Real* upperbound, /**< upper bound */
2169  SCIP_Real* lpobjval, /**< reduced cost value */
2170  SCIP_Real* bestlpobjval, /**< best reduced cost value */
2171  SCIP_Real* minpathcost, /**< the required reduced path cost to be surpassed */
2172  SCIP_Real oldupperbound, /**< old upper bound */
2173  const int* result, /**< sol int array */
2174  int* vbase, /**< array for Voronoi bases */
2175  int* state, /**< int 4 * nnodes array for internal computations */
2176  int* pathedge, /**< array for predecessor edge on a path */
2177  STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
2178  STP_Bool* nodearrchar, /**< node array storing solution vertices */
2179  SCIP_Bool* solgiven, /**< is sol given? */
2180  int extnedges /**< number of edges for transgraph */
2181 )
2182 {
2183  const SCIP_Bool dualisvalid = dualCostIsValid(scip, transgraph, bestcost, state, nodearrchar);
2184 
2185  if( !dualisvalid )
2186  {
2187  *bestlpobjval = *lpobjval;
2188  BMScopyMemoryArray(bestcost, cost, extnedges);
2189  }
2190 
2191  /* has upperbound and changed, and is best reduced cost valid and different from cost? */
2192  if( SCIPisGT(scip, *bestlpobjval, *lpobjval) && SCIPisLT(scip, *upperbound, oldupperbound) )
2193  {
2194  *solgiven = *solgiven && graph_sol_unreduced(scip, graph, result);
2195 
2196  assert(!(*solgiven) || graph_sol_valid(scip, graph, result));
2197 
2198  *minpathcost = *upperbound - *bestlpobjval;
2199  assert(SCIPisGE(scip, *minpathcost, 0.0));
2200 
2201  computeTransVoronoi(scip, transgraph, vnoi, bestcost, costrev, pathdist, vbase, pathedge);
2202 
2203  return reducePcMw(scip, graph, transgraph, vnoi, bestcost, pathdist, *minpathcost, result, marked, nodearrchar, *solgiven);
2204  }
2205  return 0;
2206 }
2207 
2208 
2209 /** check (directed) edge */
2210 static
2212  SCIP* scip, /**< SCIP data structure */
2213  const GRAPH* graph, /**< graph data structure */
2214  int root, /**< graph root from dual ascent */
2215  const SCIP_Real* redcost, /**< reduced costs */
2216  const SCIP_Real* pathdist, /**< shortest path distances */
2217  const PATH* vnoi, /**< Voronoi paths */
2218  SCIP_Real cutoff, /**< cutoff value */
2219  int edge, /**< (directed) edge to be checked */
2220  SCIP_Bool equality, /**< allow equality? */
2221  int* edgestack, /**< array of size nodes for internal computations */
2222  SCIP_Bool* deletable, /**< is edge deletable? */
2223  unsigned int* eqstack, /**< stores edges that were used for equality comparison */
2224  int* eqstack_size, /**< pointer to size of eqstack */
2225  SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
2226 )
2227 {
2228  const int head = graph->head[edge];
2229  const int tail = graph->tail[edge];
2230  const int maxgrad = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINGRAD : STP_DAEX_MAXGRAD;
2231  SCIP_Real edgebound = redcost[edge] + pathdist[tail] + vnoi[head].dist;
2232 
2233  assert(scip != NULL);
2234  assert(graph != NULL);
2235  assert(redcost != NULL);
2236  assert(pathdist != NULL);
2237  assert(vnoi != NULL);
2238  assert(edge >= 0 && edge < graph->edges);
2239 
2240  *deletable = FALSE;
2241 
2242  /* trivial rule-out? */
2243  if( SCIPisGT(scip, edgebound, cutoff) || (equality && SCIPisEQ(scip, edgebound, cutoff)) || head == root )
2244  {
2245  *deletable = TRUE;
2246  }
2247  else if( (graph->grad[tail] <= maxgrad && !Is_term(graph->term[tail]))
2248  || (graph->grad[head] <= maxgrad && !Is_term(graph->term[head])) )
2249  {
2250  SCIP_Real treecosts[STP_DAEX_MAXDFSDEPTH + 1];
2251  int treeedges[STP_DAEX_MAXDFSDEPTH + 1];
2252  SCIP_Real basebottlenecks[3];
2253  int* nodepos;
2254  int* ancestormark;
2255 
2256  const int nnodes = graph->knots;
2257  const int maxdfsdepth = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINDFSDEPTH : STP_DAEX_MAXDFSDEPTH;
2258 
2259  /* allocate clean arrays */
2260  SCIP_CALL( SCIPallocCleanBufferArray(scip, &nodepos, nnodes) );
2261  SCIP_CALL( SCIPallocCleanBufferArray(scip, &ancestormark, (MAX(graph->edges, graph->orgedges) / 2)) );
2262 
2263  basebottlenecks[0] = graph->cost[edge];
2264 #ifndef NEBUG
2265  basebottlenecks[1] = -1.0;
2266  basebottlenecks[2] = -1.0;
2267 
2268  for( int i = 0; i < STP_DAEX_MAXDFSDEPTH + 1; i++ )
2269  {
2270  treecosts[i] = -1.0;
2271  treeedges[i] = -1;
2272  }
2273 #endif
2274 
2275  markAncestors(graph, edge, ancestormark);
2276 
2277  /* check whether to extend from head */
2278  if( !Is_term(graph->term[head]) && graph->grad[head] <= maxgrad )
2279  {
2280  const SCIP_Real treecostoffset = redcost[edge] + pathdist[tail];
2281  SCIP_Real minbound = FARAWAY;
2282  SCIP_Real treecost = treecostoffset;
2283  unsigned int rounds = 0;
2284  int dfsdepth = 0;
2285  int nadded_edges = 0;
2286  SCIP_Bool stopped = FALSE;
2287 
2288  nodepos[tail] = nnodes + 1;
2289  nodepos[head] = nnodes + 4;
2290 
2291  for( int e = graph->outbeg[head]; e != EAT_LAST; e = graph->oeat[e] )
2292  if( graph->head[e] != tail )
2293  edgestack[nadded_edges++] = e;
2294 
2295  /* limited DFS */
2296  while( nadded_edges > 0 )
2297  {
2298  const int curredge = edgestack[--nadded_edges];
2299  const int currhead = graph->head[curredge];
2300 
2301  /* subtree already processed? */
2302  if( nodepos[currhead] )
2303  {
2304  assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2305 
2306  treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2307  currhead, nodepos, ancestormark, &rounds);
2308 
2309  dfsdepth--;
2310  }
2311  else
2312  {
2313  const SCIP_Real extendedcost = treecost + redcost[curredge] + vnoi[currhead].dist;
2314 
2315  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, curredge, equality,
2316  ancestormark, eqstack, eqstack_size, eqmark) )
2317  continue;
2318 
2319  if( truncateSubtree(graph, extendedcost, root, currhead, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2320  {
2321  /* stopped and no further improvement of bound possible? */
2322  if( stopped && minbound <= edgebound )
2323  break;
2324 
2325  continue;
2326  }
2327 
2328  nadded_edges = extendSubtreeHead(scip, graph, redcost, curredge, currhead, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2329  treeedges, edgestack, ancestormark);
2330  dfsdepth++;
2331  }
2332  } /* DFS loop */
2333 
2334  assert(stopped || minbound == FARAWAY);
2335  assert(SCIPisGE(scip, minbound, redcost[edge] + pathdist[tail] + vnoi[head].dist));
2336 
2337  finalizeSubtree(graph, graph->head, treeedges, dfsdepth, stopped, minbound, &edgebound, deletable, nodepos, ancestormark);
2338  } /* extend from head */
2339 
2340  /* check whether to extend from tail */
2341  if( !(*deletable) && !Is_term(graph->term[tail]) && graph->grad[tail] <= maxgrad )
2342  {
2343  const SCIP_Real treecostoffset = edgebound - pathdist[tail];
2344  SCIP_Real minbound = FARAWAY;
2345  SCIP_Real treecost = treecostoffset;
2346  int dfsdepth = 0;
2347  int nadded_edges = 0;
2348  unsigned int rounds = 0;
2349  SCIP_Bool stopped = FALSE;
2350 
2351  nodepos[head] = nnodes + 1;
2352  nodepos[tail] = nnodes + 4;
2353 
2354  for( int e = graph->inpbeg[tail]; e != EAT_LAST; e = graph->ieat[e] )
2355  if( graph->tail[e] != head )
2356  edgestack[nadded_edges++] = e;
2357 
2358  /* limited DFS */
2359  while( nadded_edges > 0 )
2360  {
2361  const int curredge = edgestack[--nadded_edges];
2362  const int currtail = graph->tail[curredge];
2363 
2364  /* subtree already processed? */
2365  if( nodepos[currtail] )
2366  {
2367  assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2368 
2369  treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2370  currtail, nodepos, ancestormark, &rounds);
2371 
2372  dfsdepth--;
2373  }
2374  else
2375  {
2376  const SCIP_Real extendedcost = treecost + redcost[curredge] + pathdist[currtail];
2377 
2378  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, flipedge((unsigned)curredge), equality,
2379  ancestormark, eqstack, eqstack_size, eqmark) )
2380  continue;
2381 
2382  if( truncateSubtree(graph, extendedcost, -1, currtail, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2383  {
2384  if( stopped )
2385  break;
2386  continue;
2387  }
2388 
2389  nadded_edges = extendSubtreeTail(graph, redcost, curredge, currtail, dfsdepth, nadded_edges, &treecost, treecosts, nodepos, treeedges,
2390  edgestack, ancestormark);
2391  dfsdepth++;
2392  }
2393  } /* DFS loop */
2394 
2395  assert(stopped || minbound == FARAWAY);
2396  assert(SCIPisGE(scip, minbound, redcost[edge] + pathdist[tail] + vnoi[head].dist));
2397 
2398  finalizeSubtree(graph, graph->tail, treeedges, dfsdepth, stopped, minbound, &edgebound, deletable, nodepos, ancestormark);
2399  } /* extend from tail */
2400 
2401  /* clean arrays */
2402  nodepos[head] = 0;
2403  nodepos[tail] = 0;
2404  unmarkAncestors(graph, edge, ancestormark);
2405 
2406  /* free memory */
2407  for( int i = 0; i < nnodes; i++ )
2408  assert(nodepos[i] == 0);
2409  for( int i = 0; i < MAX(graph->edges, graph->orgedges) / 2; i++ )
2410  assert(ancestormark[i] == 0);
2411 
2412  SCIPfreeCleanBufferArray(scip, &ancestormark);
2413  SCIPfreeCleanBufferArray(scip, &nodepos);
2414  }
2415 
2416  return SCIP_OKAY;
2417 }
2418 
2419 
2420 /** todo don't use this function, but check for conflicts in misc_stp.c and handle them */
2422  SCIP* scip, /**< SCIP data structure */
2423  GRAPH* g /**< the graph */
2424 )
2425 {
2426  int* edgemark;
2427  const int nedges = g->edges;
2428  const int nancestors = MAX(g->edges, g->orgedges) / 2;
2429 
2430  assert(scip != NULL && g != NULL);
2431  assert(g->ancestors != NULL);
2432 
2433  SCIP_CALL( SCIPallocBufferArray(scip, &edgemark, nancestors) );
2434 
2435  for( int e = 0; e < nancestors; e++ )
2436  edgemark[e] = FALSE;
2437 
2438  for( int e = 0; e < nedges; e += 2 )
2439  {
2440  SCIP_Bool conflict;
2441 
2442  if( g->oeat[e] == EAT_FREE )
2443  continue;
2444 
2445  conflict = markAncestorsConflict(g, e, edgemark);
2446 
2447  unmarkAncestorsConflict(g, e, edgemark);
2448  if( conflict )
2449  {
2450  conflict = markAncestorsConflict(g, e + 1, edgemark);
2451  unmarkAncestorsConflict(g, e + 1, edgemark);
2452 
2453  if( conflict )
2454  graph_edge_del(scip, g, e, TRUE);
2455  }
2456  }
2457 
2458  SCIPfreeBufferArray(scip, &edgemark);
2459 
2460  return SCIP_OKAY;
2461 }
2462 
2463 /** check 3-tree */
2465  SCIP* scip, /**< SCIP data structure */
2466  const GRAPH* graph, /**< graph data structure */
2467  int root, /**< graph root from dual ascent */
2468  const SCIP_Real* redcost, /**< reduced costs */
2469  const SCIP_Real* pathdist, /**< shortest path distances */
2470  const PATH* vnoi, /**< Voronoi paths */
2471  const int* vbase, /**< bases to Voronoi paths */
2472  SCIP_Real cutoff, /**< cutoff value */
2473  const int* outedges, /**< two outgoing edge */
2474  int inedge, /**< incoming edge */
2475  int* edgestack, /**< array of size nodes for internal computations */
2476  SCIP_Real* treebound, /**< to store a lower bound for the tree */
2477  SCIP_Bool* ruleout, /**< could tree be ruled out? */
2478  unsigned int* eqstack, /**< stores edges that were used for equality comparison */
2479  int* eqstack_size, /**< pointer to size of eqstack */
2480  SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
2481 )
2482 {
2483 #ifndef NDEBUG
2484  const SCIP_Real orgtreebound = *treebound;
2485 #endif
2486  assert(scip != NULL);
2487  assert(graph != NULL);
2488  assert(redcost != NULL);
2489  assert(ruleout != NULL);
2490  assert(pathdist != NULL);
2491  assert(vnoi != NULL);
2492  assert(vbase != NULL);
2493  assert(outedges != NULL);
2494  assert(inedge >= 0 && outedges[0] >= 0 && outedges[1] >= 0);
2495  assert(inedge < graph->edges && outedges[0] < graph->edges && outedges[1] < graph->edges);
2496 
2497  /* trivial rule-out? */
2498  if( SCIPisGT(scip, *treebound, cutoff) )
2499  {
2500  *ruleout = TRUE;
2501  }
2502  else
2503  {
2504  SCIP_Real treecosts[STP_DAEX_MAXDFSDEPTH + 1];
2505  int treeedges[STP_DAEX_MAXDFSDEPTH + 1];
2506  SCIP_Real basebottlenecks[3];
2507  const SCIP_Real basecost = redcost[inedge] + redcost[outedges[0]] + redcost[outedges[1]] + pathdist[graph->tail[inedge]];
2508  int* nodepos;
2509  int* ancestormark;
2510  const int nnodes = graph->knots;
2511  const int innode = graph->tail[inedge];
2512  const int basenode = graph->head[inedge];
2513  const int maxdfsdepth = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINDFSDEPTH : STP_DAEX_MAXDFSDEPTH;
2514  const int maxgrad = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINGRAD : STP_DAEX_MAXGRAD;
2515 
2516  assert(basenode == graph->tail[outedges[0]] && basenode == graph->tail[outedges[1]]);
2517 
2518  /* allocate clean array */
2519  SCIP_CALL(SCIPallocCleanBufferArray(scip, &nodepos, nnodes));
2520  SCIP_CALL( SCIPallocCleanBufferArray(scip, &ancestormark, (MAX(graph->edges, graph->orgedges) / 2)) );
2521 
2522  nodepos[basenode] = nnodes + 1; /* basenode */
2523  nodepos[innode] = nnodes + 2; /* innode */
2524 
2525 #ifndef NDEBUG
2526  for( int i = 0; i < STP_DAEX_MAXDFSDEPTH + 1; i++ )
2527  {
2528  treecosts[i] = -1.0;
2529  treeedges[i] = -1;
2530  }
2531 #endif
2532 
2533  *ruleout = FALSE;
2534 
2535  /* can we rule out subgraph beforehand? */
2536  if( markAncestorsConflict(graph, inedge, ancestormark)
2537  || markAncestorsConflict(graph, outedges[0], ancestormark)
2538  || markAncestorsConflict(graph, outedges[1], ancestormark) )
2539  {
2540  *ruleout = TRUE;
2541  }
2542 
2543  /* main loop: extend tree and try to rule it out */
2544  for( int i = 0; i < 2 && !(*ruleout); i++ )
2545  {
2546  SCIP_Real minbound = FARAWAY;
2547  SCIP_Real treecost = basecost;
2548  const int startedge = outedges[i];
2549  const int costartedge = outedges[(i + 1) % 2];
2550  const int startnode = graph->head[startedge];
2551  const int costartnode = graph->head[costartedge];
2552  int dfsdepth = 0;
2553  int nadded_edges = 0;
2554  SCIP_Bool stopped = FALSE;
2555  unsigned int rounds = 0;
2556 
2557  /* try to extend the tree from startnode */
2558 
2559  assert(startnode != root && costartnode != root);
2560 
2561  /* for basenode */
2562  basebottlenecks[0] = graph->cost[startedge];
2563 
2564  /* for innode */
2565  basebottlenecks[1] = MAX(graph->cost[startedge], graph->cost[inedge]);
2566 
2567  /* for costartnode */
2568  basebottlenecks[2] = MAX(graph->cost[startedge], graph->cost[costartedge]);
2569 
2570  nodepos[costartnode] = nnodes + 3;
2571  nodepos[startnode] = nnodes + 4;
2572 
2573  /* can we rule out entire subtree already? */
2574  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, FARAWAY, 0.0, 0, startedge, FALSE,
2575  NULL, NULL, NULL, NULL) )
2576  {
2577  *ruleout = TRUE;
2578  break;
2579  }
2580 
2581  /* cannot extend over terminals or vertices of high degree */
2582  if( Is_term(graph->term[startnode]) || graph->grad[startnode] > maxgrad )
2583  continue;
2584 
2585  assert(nodepos[startnode]);
2586 
2587  for( int e = graph->outbeg[startnode]; e != EAT_LAST; e = graph->oeat[e] )
2588  if( !nodepos[graph->head[e]] )
2589  edgestack[nadded_edges++] = e;
2590 
2591  /* limited DFS starting from startnode */
2592  while ( nadded_edges > 0 )
2593  {
2594  const int curredge = edgestack[--nadded_edges];
2595  const int currhead = graph->head[curredge];
2596 
2597  /* subtree already processed? */
2598  if( nodepos[currhead] )
2599  {
2600  assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2601 
2602  treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, basecost, dfsdepth,
2603  currhead, nodepos, ancestormark, &rounds);
2604 
2605  dfsdepth--;
2606  }
2607  else
2608  {
2609  SCIP_Real extendedcost = treecost + redcost[curredge];
2610 
2611  if( vbase[costartnode] == vbase[currhead] )
2612  {
2613  /* also covers the case that leaf is a terminal */
2614  const SCIP_Real costartfar = vnoi[costartnode + nnodes].dist;
2615  const SCIP_Real currentfar = vnoi[currhead + nnodes].dist;
2616 
2617  assert(vbase[costartnode + nnodes] >= 0 || costartfar == FARAWAY);
2618  assert(vbase[currhead + nnodes] >= 0 || currentfar == FARAWAY);
2619  extendedcost += MIN(costartfar + vnoi[currhead].dist, vnoi[costartnode].dist + currentfar);
2620  }
2621  else
2622  extendedcost += vnoi[costartnode].dist + vnoi[currhead].dist;
2623 
2624  /* can we rule out subtree? */
2625  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, curredge, FALSE,
2626  ancestormark, eqstack, eqstack_size, eqmark) )
2627  continue;
2628 
2629  /* do we need to stop extension? */
2630  if( truncateSubtree(graph, extendedcost, root, currhead, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2631  {
2632  if( stopped && minbound <= (*treebound) )
2633  break;
2634 
2635  continue;
2636  }
2637 
2638  nadded_edges = extendSubtreeHead(scip, graph, redcost, curredge, currhead, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2639  treeedges, edgestack, ancestormark);
2640  dfsdepth++;
2641  }
2642  } /* DFS loop */
2643 
2644  assert(stopped || minbound == FARAWAY);
2645 #ifndef NDEBUG
2646  assert(SCIPisGE(scip, minbound, orgtreebound));
2647 #endif
2648 
2649  finalizeSubtree(graph, graph->head, treeedges, dfsdepth, stopped, minbound, treebound, ruleout, nodepos, ancestormark);
2650  } /* main loop for outgoing edges */
2651 
2652 #ifndef NDEBUG
2653  assert(*treebound >= orgtreebound);
2654 #endif
2655 
2656  if( !(*ruleout) && !Is_term(graph->term[innode]) && graph->grad[innode] <= maxgrad )
2657  {
2658  /* move down the incoming edge */
2659  const SCIP_Real treecostoffset = *treebound - pathdist[innode];
2660  SCIP_Real minbound = FARAWAY;
2661  SCIP_Real treecost = treecostoffset;
2662  int dfsdepth = 0;
2663  int nadded_edges = 0;
2664  SCIP_Bool stopped = FALSE;
2665  unsigned int rounds = 0;
2666 
2667  assert(treecost >= 0.0);
2668  assert(nodepos[innode] == nnodes + 2);
2669  assert(nodepos[basenode] == nnodes + 1);
2670 
2671  nodepos[graph->head[outedges[0]]] = nnodes + 2;
2672  nodepos[graph->head[outedges[1]]] = nnodes + 3;
2673  nodepos[innode] = nnodes + 4;
2674 
2675  /* for basenode */
2676  basebottlenecks[0] = graph->cost[inedge];
2677 
2678  /* for outedges[0] */
2679  basebottlenecks[1] = MAX(graph->cost[outedges[0]], graph->cost[inedge]);
2680 
2681  /* for outedges[1] */
2682  basebottlenecks[2] = MAX(graph->cost[outedges[1]], graph->cost[inedge]);
2683 
2684  /* can we rule out entire subtree already? */
2685  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, FARAWAY, 0.0, 0, flipedge(inedge), FALSE,
2686  NULL, NULL, NULL, NULL) )
2687  {
2688  *ruleout = TRUE;
2689  }
2690  else
2691  {
2692  for( int e = graph->inpbeg[innode]; e != EAT_LAST; e = graph->ieat[e] )
2693  if( !nodepos[graph->tail[e]] )
2694  edgestack[nadded_edges++] = e;
2695  }
2696 
2697  /* limited DFS starting from inedge */
2698  while( nadded_edges > 0 )
2699  {
2700  const int curredge = edgestack[--nadded_edges];
2701  const int currtail = graph->tail[curredge];
2702 
2703  /* subtree already processed? */
2704  if( nodepos[currtail] )
2705  {
2706  assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2707 
2708  treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2709  currtail, nodepos, ancestormark, &rounds);
2710 
2711  dfsdepth--;
2712  }
2713  else
2714  {
2715  const SCIP_Real extendedcost = treecost + redcost[curredge] + pathdist[currtail];
2716 
2717  if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, flipedge((unsigned)curredge), FALSE,
2718  ancestormark, eqstack, eqstack_size, eqmark) )
2719  continue;
2720 
2721  if( truncateSubtree(graph, extendedcost, -1, currtail, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2722  {
2723  if( stopped && minbound <= (*treebound) )
2724  break;
2725 
2726  continue;
2727  }
2728 
2729  nadded_edges = extendSubtreeTail(graph, redcost, curredge, currtail, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2730  treeedges, edgestack, ancestormark);
2731  dfsdepth++;
2732  }
2733  } /* DFS loop */
2734 
2735 #ifndef NDEBUG
2736  assert(stopped || minbound == FARAWAY);
2737  assert(SCIPisGE(scip, minbound, orgtreebound));
2738 #endif
2739 
2740  finalizeSubtree(graph, graph->tail, treeedges, dfsdepth, stopped, minbound, treebound, ruleout, nodepos, ancestormark);
2741  } /* inedge */
2742 
2743  /* clean nodepos array */
2744  nodepos[basenode] = 0;
2745  nodepos[innode] = 0;
2746 
2747  unmarkAncestorsConflict(graph, inedge, ancestormark);
2748 
2749  for( int i = 0; i < 2; i++ )
2750  {
2751  nodepos[graph->head[outedges[i]]] = 0;
2752  unmarkAncestorsConflict(graph, outedges[i], ancestormark);
2753  }
2754 
2755 #ifndef NDEBUG
2756  assert(*treebound >= orgtreebound);
2757 
2758  for( int i = 0; i < nnodes; i++ )
2759  assert(nodepos[i] == 0);
2760 
2761  for( int i = 0; i < MAX(graph->edges, graph->orgedges) / 2; i++ )
2762  assert(ancestormark[i] == 0);
2763 #endif
2764 
2765  SCIPfreeCleanBufferArray(scip, &ancestormark);
2766  SCIPfreeCleanBufferArray(scip, &nodepos);
2767  }
2768 
2769  return SCIP_OKAY;
2770 }
2771 
2772 
2773 /** reduce SPG graph based on reduced cost information and given upper bound */
2775  SCIP* scip, /**< SCIP data structure */
2776  GRAPH* graph, /**< graph data structure */
2777  const PATH* vnoi, /**< Voronoi data structure */
2778  const SCIP_Real* cost, /**< dual ascent costs */
2779  const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2780  const int* result, /**< sol int array */
2781  SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
2782  int root, /**< the root */
2783  int* nodearr, /**< for internal stuff */
2784  STP_Bool* marked /**< edge array to mark which (directed) edge can be removed */
2785 )
2786 {
2787  unsigned int* eqstack;
2788  SCIP_Bool* eqmark;
2789  int nfixed = 0;
2790  const int nedges = graph->edges;
2791  const int nnodes = graph->knots;
2792  const int halfnedges = graph->edges / 2;
2793  const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
2794 
2795  if( rpc )
2796  graph_pc_2orgcheck(graph);
2797 
2798  if( SCIPisZero(scip, minpathcost) )
2799  return 0;
2800 
2801  SCIP_CALL( SCIPallocCleanBufferArray(scip, &eqmark, halfnedges) );
2802  SCIP_CALL( SCIPallocBufferArray(scip, &eqstack, halfnedges) );
2803 
2804  /* main loop */
2805  for( int e = 0; e < nedges; e += 2 )
2806  {
2807  if( graph->oeat[e] != EAT_FREE )
2808  {
2809  const int erev = e + 1;
2810  int eqstack_size = 0;
2811  SCIP_Bool deletable = TRUE;
2812  const SCIP_Bool allowequality = (result != NULL && result[e] != CONNECT && result[erev] != CONNECT);
2813 
2814  assert(graph->oeat[erev] != EAT_FREE);
2815 
2816  if( SCIPisZero(scip, cost[e]) || SCIPisZero(scip, cost[erev]) )
2817  continue;
2818 
2819  if( marked == NULL || !marked[e] )
2820  {
2821  SCIP_CALL_ABORT(reduceCheckEdge(scip, graph, root, cost, pathdist, vnoi, minpathcost, e, allowequality, nodearr,
2822  &deletable, eqstack, &eqstack_size, eqmark));
2823 
2824  if( marked != NULL && deletable )
2825  marked[e] = TRUE;
2826  }
2827 
2828  if( (marked == NULL || !marked[erev]) && deletable )
2829  {
2830  SCIP_CALL_ABORT(reduceCheckEdge(scip, graph, root, cost, pathdist, vnoi, minpathcost, erev, allowequality,
2831  nodearr, &deletable, eqstack, &eqstack_size, eqmark));
2832 
2833  if( marked != NULL && deletable )
2834  marked[erev] = TRUE;
2835  }
2836 
2837  for( int i = 0; i < eqstack_size; i++ )
2838  eqmark[eqstack[i]] = FALSE;
2839  for( int i = 0; i < halfnedges; i++ )
2840  assert(eqmark[i] == FALSE);
2841 
2842  if( deletable )
2843  {
2844  graph_edge_del(scip, graph, e, TRUE);
2845  nfixed++;
2846  }
2847  }
2848  }
2849 
2850  SCIPfreeBufferArray(scip, &eqstack);
2851  SCIPfreeCleanBufferArray(scip, &eqmark);
2852 
2853  for( int k = 0; k < nnodes; k++ )
2854  if( graph->grad[k] == 0 && k != root )
2855  graph->mark[k] = FALSE;
2856 
2857  return nfixed;
2858 }
2859 
2860 /** dual ascent based reductions */
2862  SCIP* scip, /**< SCIP data structure */
2863  GRAPH* graph, /**< graph data structure */
2864  PATH* vnoi, /**< Voronoi data structure */
2865  GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
2866  SCIP_Real* cost, /**< edge costs */
2867  SCIP_Real* costrev, /**< reverse edge costs */
2868  SCIP_Real* pathdist, /**< distance array for shortest path calculations */
2869  SCIP_Real* ub, /**< pointer to provide upper bound and return upper bound found during ascent and prune (if better) */
2870  SCIP_Real* offset, /**< pointer to store offset */
2871  int* edgearrint, /**< int edges array for internal computations or NULL */
2872  int* vbase, /**< array for Voronoi bases */
2873  int* state, /**< int 4 * nnodes array for internal computations */
2874  int* pathedge, /**< array for predecessor edge on a path */
2875  int* nodearrint, /**< int nnodes array for internal computations */
2876  int* heursources, /**< array to store starting points of TM heuristic */
2877  STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
2878  int* nelims, /**< pointer to store number of reduced edges */
2879  int prevrounds, /**< number of reduction rounds that have been performed already */
2880  SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
2881  SCIP_Bool userec, /**< use recombination heuristic? */
2882  SCIP_Bool extended /**< use extended tests? */
2883 )
2884 {
2885  STPSOLPOOL* pool;
2886 
2887  SCIP_Real* edgefixingbounds;
2888  SCIP_Real* nodefixingbounds;
2889  SCIP_Real* nodereplacebounds;
2890  SCIP_Real lpobjval;
2891  SCIP_Real upperbound;
2892  SCIP_Real minpathcost;
2893  SCIP_Real damaxdeviation;
2894  SCIP_Bool success;
2895  const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
2896  const SCIP_Bool directed = (graph->stp_type == STP_SAP || graph->stp_type == STP_NWSPG);
2897  int* terms;
2898  int* result;
2899  int* starts;
2900  int k;
2901  int runs;
2902  int nruns;
2903  int nterms;
2904  const int nedges = graph->edges;
2905  const int nnodes = graph->knots;
2906  int nfixed;
2907  int best_start;
2908  STP_Bool* marked;
2909 
2910  assert(ub != NULL);
2911  assert(scip != NULL);
2912  assert(graph != NULL);
2913  assert(nelims != NULL);
2914  assert(nodearrint != NULL);
2915 
2916  nfixed = 0;
2917  minpathcost = 0.0;
2918 
2919  if( graph->terms <= 1 )
2920  return SCIP_OKAY;
2921 
2922  if( SCIPisGE(scip, *ub, 0.0) )
2923  upperbound = *ub;
2924  else
2925  upperbound = FARAWAY;
2926 
2927  /* allocate memory */
2928  SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
2929  SCIP_CALL( SCIPallocBufferArray(scip, &marked, nedges) );
2930  SCIP_CALL( SCIPallocBufferArray(scip, &edgefixingbounds, nedges) );
2931  SCIP_CALL( SCIPallocBufferArray(scip, &nodefixingbounds, nnodes) );
2932  SCIP_CALL( SCIPallocBufferArray(scip, &nodereplacebounds, nnodes) );
2933 
2934  if( !rpc && !directed )
2935  SCIP_CALL( SCIPallocBufferArray(scip, &terms, graph->terms) );
2936  else
2937  terms = NULL;
2938 
2939  /*
2940  * 1. step: compute upper bound
2941  */
2942  if( userec )
2943  SCIP_CALL( SCIPStpHeurRecInitPool(scip, &pool, nedges, SOLPOOL_SIZE) );
2944  else
2945  pool = NULL;
2946 
2947  if( !rpc && extended )
2948  SCIP_CALL( reduce_deleteConflictEdges(scip, graph) );
2949 
2950  /* initialize */
2951  k = 0;
2952  nterms = 0;
2953  for( int i = 0; i < nnodes; i++ )
2954  {
2955  if( !rpc )
2956  graph->mark[i] = (graph->grad[i] > 0);
2957 
2958  if( graph->mark[i] )
2959  {
2960  k++;
2961  if( Is_term(graph->term[i]) )
2962  {
2963  if( !rpc && !directed )
2964  terms[nterms] = i;
2965  nterms++;
2966  }
2967  }
2968  }
2969 
2970  /* not more than two terminals? */
2971  if( nterms <= 2 )
2972  goto TERMINATE;
2973 
2974  /* number of runs should not exceed number of connected vertices */
2975  if( directed )
2976  runs = MIN(k, DEFAULT_HEURRUNS);
2977  else
2978  runs = MIN(k, DEFAULT_HEURRUNS / 5);
2979 
2980  /* neither PC, MW, RPC nor HC? */
2981  if( !rpc && heursources != NULL )
2982  {
2983  /* choose starting points for TM heuristic */
2984  starts = heursources;
2985  SCIPStpHeurTMCompStarts(graph, starts, &runs);
2986  }
2987  else
2988  {
2989  starts = NULL;
2990  }
2991 
2992  for( int e = 0; e < nedges; e++ )
2993  {
2994  cost[e] = graph->cost[e];
2995  costrev[e] = graph->cost[flipedge(e)];
2996  result[e] = UNKNOWN;
2997  }
2998 
2999  if( directed || rpc )
3000  {
3001  SCIP_Real ubnew;
3002 
3003  SCIP_CALL( SCIPStpHeurTMRun(scip, NULL, graph, starts, &best_start, result, runs, graph->source, cost, costrev, NULL, NULL, 0.0, &success, FALSE) );
3004  assert(success);
3005 
3006  /* calculate objective value of solution */
3007  ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3008 
3009  if( SCIPisLT(scip, ubnew, upperbound) )
3010  upperbound = ubnew;
3011  }
3012 
3013  /*
3014  * 2. step: repeatedly compute lower bound and reduced costs and try to eliminate
3015  */
3016 
3017 
3018  damaxdeviation = -1.0;
3019 
3020  /* if not RPC, select roots for dual ascent */
3021  if( !rpc && !directed )
3022  {
3023  assert(terms != NULL);
3024 
3025  nruns = MIN(graph->terms, DEFAULT_DARUNS);
3026  SCIP_CALL( orderDaRoots(scip, graph, terms, graph->terms, (prevrounds > 0), randnumgen) );
3027 
3028  if( prevrounds > 0 )
3030  }
3031  else
3032  {
3033  nruns = 1;
3034  if( rpc )
3035  graph_pc_2transcheck(graph);
3036  }
3037 
3038  assert(nruns > 0);
3039 
3040  for( int outerrounds = 0; outerrounds < 2; outerrounds++ )
3041  {
3042  /* main reduction loop */
3043  for( int run = 0; run < nruns; run++ )
3044  {
3045  int root = graph->source;
3046  SCIP_Bool apsol = FALSE;
3047  SCIP_Bool usesol = (run > 1);
3048 
3049  if( rpc )
3050  usesol = TRUE;
3051 
3052  /* graph vanished? */
3053  if( graph->grad[graph->source] == 0 )
3054  break;
3055 
3056  if( !rpc && !directed )
3057  {
3058  assert(terms != NULL);
3059  root = terms[run];
3060  usesol = usesol && (SCIPrandomGetInt(randnumgen, 0, 2) < 2) && graph->stp_type != STP_RSMT;
3061  }
3062 
3063  if( usesol )
3064  {
3065  /* solution might not be valid anymore */
3066  if( !graph_sol_valid(scip, graph, result) )
3067  {
3068  SCIPdebugMessage("solution not valid; run normal dual-ascent \n");
3069  SCIP_CALL(
3070  SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3071  }
3072  else
3073  {
3074  if( !rpc )
3075  {
3076  SCIPdebugMessage("reroot solution and run guided dual-ascent \n");
3077  SCIP_CALL(graph_sol_reroot(scip, graph, result, root));
3078 #ifndef NDEBUG
3079  {
3080  const int realroot = graph->source;
3081  graph->source = root;
3082  assert(graph_sol_valid(scip, graph, result));
3083  graph->source = realroot;
3084  }
3085 #endif
3086  }
3087 
3088  SCIP_CALL(
3089  SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, result, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3090  }
3091  }
3092  else
3093  {
3094  SCIPdebugMessage("no rerooting, run normal dual-ascent \n");
3095  SCIP_CALL(
3096  SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3097  }
3098 
3099  /* perform ascent and prune */
3100  if( !directed )
3101  {
3102  SCIP_Real ubnew;
3103  SCIP_Bool soladded = FALSE;
3104 
3105  SCIP_CALL(SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, result, nodearrint, root, nodearrchar, &success, FALSE));
3106  assert(success && graph_sol_valid(scip, graph, result));
3107 
3108  ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3109 
3110  if( userec && !rpc )
3111  {
3112  SCIPdebugMessage("obj before local %f \n", ubnew);
3113 
3114  SCIP_CALL(SCIPStpHeurLocalRun(scip, graph, graph->cost, result));
3115  ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3116 
3117  SCIPdebugMessage("obj after local %f \n", ubnew);
3118  assert(graph_sol_valid(scip, graph, result));
3119 
3120  SCIP_CALL(SCIPStpHeurRecAddToPool(scip, ubnew, result, pool, &soladded));
3121  }
3122 
3123  if( userec && soladded && pool->size >= 2 && ubnew < upperbound )
3124  {
3125  /* get index of just added solution */
3126  int solindex = pool->maxindex;
3127  SCIP_Bool solfound;
3128 
3129  SCIPdebugMessage("POOLSIZE %d \n", pool->size);
3130 
3131  SCIP_CALL(SCIPStpHeurRecRun(scip, pool, NULL, NULL, graph, NULL, &solindex, 1, pool->size, FALSE, &solfound));
3132 
3133  if( solfound )
3134  {
3135  const STPSOL* const sol = SCIPStpHeurRecSolfromIdx(pool, solindex);
3136 
3137  assert(sol != NULL);
3138 
3139  SCIPdebugMessage("DA: rec found better solution with obj %f vs %f \n", sol->obj, ubnew);
3140 
3141  if( sol->obj < ubnew )
3142  {
3143  BMScopyMemoryArray(result, sol->soledges, nedges);
3144 
3145  SCIP_CALL(SCIPStpHeurLocalRun(scip, graph, graph->cost, result));
3146  ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3147 
3148  assert(SCIPisLE(scip, ubnew, sol->obj));
3149 
3150  if( ubnew < sol->obj )
3151  SCIP_CALL(SCIPStpHeurRecAddToPool(scip, ubnew, result, pool, &solfound));
3152  }
3153  }
3154  }
3155 
3156  if( SCIPisLE(scip, ubnew, upperbound) )
3157  {
3158  apsol = TRUE;
3159  upperbound = ubnew;
3160  }
3161  }
3162 
3163  /* the required reduced path cost to be surpassed */
3164  minpathcost = upperbound - lpobjval;
3165  assert(SCIPisGE(scip, minpathcost, 0.0));
3166 
3167  SCIPdebugMessage("upper: %f lower: %f \n", upperbound, lpobjval);
3168 
3169  for( k = 0; k < nnodes; k++ )
3170  graph->mark[k] = (graph->grad[k] > 0);
3171 
3172  /* distance from root to all nodes */
3173  graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
3174 
3175  for( unsigned e = 0; e < (unsigned) nedges; e++ )
3176  {
3177  marked[e] = FALSE;
3178  costrev[e] = cost[flipedge(e)];
3179  }
3180 
3181  /* no paths should go back to the root */
3182  for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
3183  costrev[e] = FARAWAY;
3184 
3185  /* build Voronoi diagram */
3186  if( directed || rpc )
3187  graph_voronoiTerms(scip, graph, costrev, vnoi, vbase, graph->path_heap, state);
3188  else
3189  {
3190  graph_get4nextTerms(scip, graph, costrev, costrev, vnoi, vbase, graph->path_heap, state);
3191 
3192 #ifndef NDEBUG
3193  for( int i = 0; i < nnodes; i++ )
3194  if( !Is_term(graph->term[i]) )
3195  {
3196  assert(vbase[i] != root || vnoi[i].dist >= FARAWAY);
3197  assert(vbase[i + nnodes] != root || vnoi[i + nnodes].dist >= FARAWAY);
3198  }
3199  else
3200  assert(vbase[i] == i);
3201 #endif
3202  updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, (run == 0));
3203  updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, nedges, (run == 0), TRUE);
3204  }
3205 
3206  nfixed += reduceSPG(scip, graph, offset, marked, nodearrchar, vnoi, cost, pathdist, result, minpathcost, root, apsol);
3207 
3208  if( !directed && !rpc )
3209  {
3210  if( !SCIPisZero(scip, minpathcost) )
3211  {
3212  nfixed += reduceWithNodeFixingBounds(scip, graph, NULL, nodefixingbounds, upperbound);
3213  nfixed += reduceWithEdgeFixingBounds(scip, graph, NULL, edgefixingbounds, (apsol ? result : NULL), upperbound);
3214  }
3215 
3216  if( extended )
3217  {
3218  nfixed += reduce_extendedEdge(scip, graph, vnoi, cost, pathdist, (apsol ? result : NULL), minpathcost, root, nodearrint, marked);
3219  }
3220 
3221  if( !SCIPisZero(scip, minpathcost) )
3222  SCIP_CALL(
3223  updateNodeReplaceBounds(scip, nodereplacebounds, graph, cost, pathdist, vnoi, vbase, nodearrint, lpobjval, upperbound, root,
3224  (run == 0), extended));
3225  }
3226 
3227  if( nfixed > 0 )
3228  SCIP_CALL(level0(scip, graph));
3229  assert(graph_valid(graph));
3230 
3231  if( !rpc )
3232  {
3233  for( k = 0; k < nnodes; k++ )
3234  graph->mark[k] = graph->grad[k] > 0;
3235  }
3236  else
3237  {
3238  graph_pc_2trans(graph);
3239  graph_pc_2org(graph);
3240  }
3241  } /* root loop */
3242 
3243  if( !directed && !rpc && !SCIPisZero(scip, minpathcost) )
3244  nfixed += reduceWithNodeReplaceBounds(scip, graph, vnoi, pathdist, cost, nodereplacebounds, nodearrint, lpobjval, upperbound);
3245 
3246  if( nfixed == 0 || !userec )
3247  {
3248  break;
3249  }
3250  else if( userec )
3251  {
3253  }
3254  }
3255 
3256 TERMINATE:
3257  *nelims = nfixed;
3258 
3259  if( SCIPisLT(scip, upperbound, *ub) || SCIPisLT(scip, *ub, 0.0) )
3260  *ub = upperbound;
3261 
3262  if( userec )
3263  SCIPStpHeurRecFreePool(scip, &pool);
3264 
3265  SCIPfreeBufferArrayNull(scip, &terms);
3266  SCIPfreeBufferArray(scip, &nodereplacebounds);
3267  SCIPfreeBufferArray(scip, &nodefixingbounds);
3268  SCIPfreeBufferArray(scip, &edgefixingbounds);
3269  SCIPfreeBufferArray(scip, &marked);
3270  SCIPfreeBufferArray(scip, &result);
3271 
3272  assert(graph_valid(graph));
3273 
3274  return SCIP_OKAY;
3275 }
3276 
3277 
3278 /** dual ascent reduction for slack-and-prune heuristic */
3280  SCIP* scip, /**< SCIP data structure */
3281  SCIP_VAR** vars, /**< problem variables or NULL */
3282  GRAPH* graph, /**< graph data structure */
3283  PATH* vnoi, /**< Voronoi data structure */
3284  GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
3285  SCIP_Real* cost, /**< array to store reduced costs */
3286  SCIP_Real* costrev, /**< reverse edge costs */
3287  SCIP_Real* pathdist, /**< distance array for shortest path calculations */
3288  SCIP_Real* upperbound, /**< pointer to store new upper bound */
3289  int* edgearrint, /**< int edges array to store solution value */
3290  int* edgearrint2, /**< int edges array for internal computations */
3291  int* vbase, /**< array for Voronoi bases */
3292  int* pathedge, /**< array for predecessor edge on a path */
3293  int* state, /**< int 4 * nnodes array for internal computations */
3294  int* solnode, /**< array of nodes of current solution that is not to be destroyed */
3295  STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
3296  STP_Bool* edgearrchar, /**< STP_Bool edge array for internal computations */
3297  int* nelims, /**< pointer to store number of reduced edges */
3298  int minelims, /**< minimum number of edges to eliminate */
3299  SCIP_Bool solgiven /**< solution provided? */
3300  )
3301 {
3302  IDX** ancestors;
3303  IDX** revancestors;
3304  SCIP_Real obj;
3305  SCIP_Real tmpcost;
3306  SCIP_Real lpobjval;
3307  SCIP_Real objprune;
3308  SCIP_Real minpathcost;
3309  SCIP_Real* sd;
3310  SCIP_Real* ecost;
3311  SCIP_Bool rpc;
3312  SCIP_Bool success;
3313  SCIP_Bool eliminate;
3314 
3315  int* grad;
3316  int* adjvert;
3317  int* incedge;
3318  int* reinsert;
3319  int i;
3320  int k;
3321  int e;
3322  int e2;
3323  int e3;
3324  int etmp;
3325  int root;
3326  int head;
3327  int nterms;
3328  int nedges;
3329  int nnodes;
3330  int nfixed;
3331  int redrounds;
3332  STP_Bool* marked;
3333 
3334  assert(scip != NULL);
3335  assert(cost != NULL);
3336  assert(graph != NULL);
3337  assert(nelims != NULL);
3338  assert(solnode != NULL);
3339  assert(costrev != NULL);
3340  assert(pathedge != NULL);
3341  assert(upperbound != NULL);
3342  assert(edgearrint != NULL);
3343  assert(edgearrint2 != NULL);
3344  assert(nodearrchar != NULL);
3345  assert(edgearrchar != NULL);
3346 
3347  /* 1. step: initialize */
3348 
3349  rpc = (graph->stp_type == STP_RPCSPG);
3350  grad = graph->grad;
3351  root = graph->source;
3352  nfixed = 0;
3353  nedges = graph->edges;
3354  nnodes = graph->knots;
3355 
3356  /* graph vanished? */
3357  if( grad[graph->source] == 0 )
3358  return SCIP_OKAY;
3359 
3360  marked = edgearrchar;
3361 
3362  /* allocate length-4 buffer memory */
3363  SCIP_CALL( SCIPallocBufferArray(scip, &sd, 4) );
3364  SCIP_CALL( SCIPallocBufferArray(scip, &ancestors, 4) );
3365  SCIP_CALL( SCIPallocBufferArray(scip, &revancestors, 4) );
3366  SCIP_CALL( SCIPallocBufferArray(scip, &ecost, 4) );
3367  SCIP_CALL( SCIPallocBufferArray(scip, &adjvert, 4) );
3368  SCIP_CALL( SCIPallocBufferArray(scip, &reinsert, 4) );
3369  SCIP_CALL( SCIPallocBufferArray(scip, &incedge, 4) );
3370 
3371  for( i = 0; i < 4; i++ )
3372  {
3373  sd[i] = 0.0;
3374  ancestors[i] = NULL;
3375  revancestors[i] = NULL;
3376  }
3377 
3378  k = 0;
3379  nterms = 0;
3380  for( i = 0; i < nnodes; i++ )
3381  {
3382  if( !rpc )
3383  graph->mark[i] = (grad[i] > 0);
3384  if( graph->mark[i] )
3385  {
3386  k++;
3387  if( Is_term(graph->term[i]) )
3388  nterms++;
3389  }
3390  }
3391 
3392  /* not more than two terminals? */
3393  if( nterms <= 2 )
3394  goto TERMINATE;
3395 
3396  /* 2. step: - if not provided, compute lower bound and reduced costs
3397  * - try to eliminate edges and nodes */
3398 
3399  for( i = 0; i < nnodes; i++ )
3400  if( Is_term(graph->term[i]) )
3401  assert(grad[i] > 0);
3402 
3403  if( rpc )
3404  {
3405  graph_pc_2trans(graph);
3406  }
3407 
3408 
3409  if( !solgiven )
3410  {
3411  /* try to build MST on solnode nodes */
3412  for( i = 0; i < nnodes; i++ )
3413  graph->mark[i] = (solnode[i] == CONNECT);
3414 
3415  for( e = 0; e < nedges; e++ )
3416  edgearrint[e] = UNKNOWN;
3417 
3418  graph_path_exec(scip, graph, MST_MODE, root, graph->cost, vnoi);
3419 
3420  for( i = 0; i < nnodes; i++ )
3421  {
3422  e = vnoi[i].edge;
3423  if( e >= 0 )
3424  {
3425  edgearrint[e] = CONNECT;
3426  }
3427  else if( Is_term(graph->term[i]) && i != root )
3428  {
3429  break;
3430  }
3431  }
3432 
3433  if( i == nnodes )
3434  {
3435  int l;
3436  int count;
3437 
3438  do
3439  {
3440  count = 0;
3441 
3442  for( l = nnodes - 1; l >= 0; --l )
3443  {
3444  if( (solnode[l] != CONNECT) || Is_term(graph->term[l]) )
3445  continue;
3446 
3447  for( e = graph->outbeg[l]; e != EAT_LAST; e = graph->oeat[e] )
3448  if( edgearrint[e] == CONNECT )
3449  break;
3450 
3451  if( e == EAT_LAST )
3452  {
3453  /* there has to be exactly one incoming edge */
3454  for( e = graph->inpbeg[l]; e != EAT_LAST; e = graph->ieat[e] )
3455  {
3456  if( edgearrint[e] == CONNECT )
3457  {
3458  edgearrint[e] = UNKNOWN;
3459  solnode[l] = UNKNOWN;
3460  count++;
3461  break;
3462  }
3463  }
3464  }
3465  }
3466  }
3467  while( count > 0 );
3468  }
3469  }
3470 
3471 
3472  if( solgiven || i == nnodes )
3473  {
3474  obj = graph_sol_getObj(graph->cost, edgearrint, 0.0, nedges);
3475 
3476  SCIP_CALL( SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, edgearrint, edgearrint2, state, root, FALSE, -1.0, nodearrchar) );
3477  }
3478  else
3479  {
3480  obj = FARAWAY;
3481  SCIP_CALL( SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint2, state, root, FALSE, -1.0, nodearrchar) );
3482  }
3483 
3484 #if 0
3485  SCIP_QUEUE* queue;
3486  SCIP_CALL( SCIPqueueCreate(&queue, nnodes, 2.0) );
3487 
3488  GRAPH* prunegraph = graph;
3489  int* mark = graph->mark;
3490  int* pnode;
3491  int a;
3492  for( k = 0; k < nnodes; k++ )
3493  mark[k] = FALSE;
3494 
3495 
3496  /* BFS from root along incoming arcs of zero cost */
3497 
3498  mark[prunegraph->source] = TRUE;
3499 
3500  SCIP_CALL( SCIPqueueInsert(queue, &(prunegraph->source)) );
3501 
3502  while( !SCIPqueueIsEmpty(queue) )
3503  {
3504  pnode = (SCIPqueueRemove(queue));
3505  k = *pnode;
3506 
3507  /* traverse outgoing arcs */
3508  for( a = prunegraph->outbeg[k]; a != EAT_LAST; a = prunegraph->oeat[a] )
3509  {
3510  head = prunegraph->head[a];
3511 
3512  if( SCIPisEQ(scip, cost[a], 0.0) )
3513  {
3514  /* vertex not labeled yet? */
3515  if( !mark[head] )
3516  {
3517  mark[head] = TRUE;
3518  SCIP_CALL( SCIPqueueInsert(queue, &(prunegraph->head[a])) );
3519  }
3520  }
3521  }
3522  }
3523  SCIPqueueFree(&queue);
3524  for( k = 0; k < nnodes; k++ ) // TODO
3525  if( Is_term(prunegraph->term[k]) && !mark[k] )
3526  printf("in bnd FAIL %d not marked, but terminal, \n", k);
3527 #endif
3528 
3529  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, edgearrint2, vbase, root, nodearrchar, &success, FALSE) );
3530 
3531  objprune = graph_sol_getObj(graph->cost, edgearrint2, 0.0, nedges);
3532 
3533  assert(success);
3534 
3535  if( success && SCIPisLT(scip, objprune, obj ) )
3536  {
3537 
3538  for( i = 0; i < nnodes; i++ )
3539  solnode[i] = UNKNOWN;
3540 
3541  for( e = 0; e < nedges; e++ )
3542  {
3543  edgearrint[e] = edgearrint2[e];
3544  if( edgearrint[e] == CONNECT )
3545  {
3546  solnode[graph->tail[e]] = CONNECT;
3547  solnode[graph->head[e]] = CONNECT;
3548  }
3549  }
3550  }
3551 
3552  obj = 0.0;
3553 
3554  for( e = 0; e < nedges; e++ )
3555  {
3556  if( edgearrint[e] == CONNECT )
3557  obj += graph->cost[e];
3558 
3559  marked[e] = FALSE;
3560  costrev[e] = cost[flipedge(e)];
3561  }
3562 
3563  *upperbound = obj;
3564 
3565  for( k = 0; k < nnodes; k++ )
3566  graph->mark[k] = (grad[k] > 0);
3567 
3568  /* distance from root to all nodes */
3569  graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
3570 
3571  /* no paths should go back to the root */
3572  for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
3573  costrev[e] = FARAWAY;
3574 
3575  /* build Voronoi diagram */
3576  graph_get4nextTerms(scip, graph, costrev, costrev, vnoi, vbase, graph->path_heap, state);
3577 
3578  for( k = 0; k < nnodes; k++ )
3579  if( !Is_term(graph->term[k]) )
3580  assert(vbase[k + nnodes] != root );
3581 
3582  /* RPC? If yes, restore original graph */
3583  if( rpc )
3584  {
3585  graph_pc_2org(graph);
3586  graph->mark[root] = FALSE;
3587  }
3588 
3589  for( e = 0; e < nedges; e++ )
3590  costrev[e] = -1.0;
3591 
3592  for( redrounds = 0; redrounds < 3; redrounds++ )
3593  {
3594  if( redrounds == 0 )
3595  {
3596  eliminate = FALSE;
3597  minpathcost = FARAWAY;
3598  }
3599  else if( redrounds == 1 )
3600  {
3601  assert(minelims > 0);
3602  assert(2 * minelims < nedges);
3603 
3604  eliminate = TRUE;
3605  SCIPsortReal(costrev, nedges);
3606 
3607  /* the required reduced path cost to be surpassed */
3608  minpathcost = costrev[nedges - 2 * minelims];
3609 
3610  if( SCIPisLE(scip, minpathcost, 0.0) )
3611  minpathcost = 2 * SCIPepsilon(scip);
3612 
3613  k = nedges - 2 * minelims;
3614 
3615  /* try to first eliminate edges with higher gap */
3616  for( e = nedges - 1; e > k && e >= 2; e = e - 2 )
3617  {
3618  if( SCIPisLE(scip, costrev[e - 2], minpathcost) )
3619  break;
3620  }
3621 
3622  if( SCIPisGT(scip, costrev[e], minpathcost) )
3623  minpathcost = costrev[e];
3624 
3625  if( SCIPisLE(scip, minpathcost, 0.0) )
3626  minpathcost = 2 * SCIPepsilon(scip);
3627 
3628  for( e = 0; e < nedges; e++ )
3629  marked[e] = FALSE;
3630  }
3631  else
3632  {
3633  eliminate = TRUE;
3634 
3635  /* the required reduced path cost to be surpassed */
3636  minpathcost = costrev[nedges - 2 * minelims];
3637 
3638  if( SCIPisLE(scip, minpathcost, 0.0) )
3639  minpathcost = 2 * SCIPepsilon(scip);
3640 
3641  for( e = 0; e < nedges; e++ )
3642  marked[e] = FALSE;
3643  }
3644 
3645  for( k = 0; k < nnodes; k++ )
3646  {
3647  if( grad[k] <= 0 )
3648  continue;
3649 
3650  if( nfixed > minelims )
3651  break;
3652 
3653  if( !Is_term(graph->term[k]) && (!eliminate || pathdist[k] + vnoi[k].dist >= minpathcost) && solnode[k] != CONNECT )
3654  {
3655  if( !eliminate )
3656  {
3657  tmpcost = pathdist[k] + vnoi[k].dist;
3658 
3659  for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
3660  {
3661  if( SCIPisGT(scip, tmpcost, costrev[e]) )
3662  costrev[e] = tmpcost;
3663 
3664  e2 = flipedge(e);
3665 
3666  if( SCIPisGT(scip, tmpcost, costrev[e2]) )
3667  costrev[e2] = tmpcost;
3668  }
3669 
3670  continue;
3671  }
3672  nfixed += grad[k];
3673 
3674  while( graph->outbeg[k] != EAT_LAST )
3675  graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
3676  }
3677  else
3678  {
3679  e = graph->outbeg[k];
3680  while( e != EAT_LAST )
3681  {
3682  etmp = graph->oeat[e];
3683  head = graph->head[e];
3684 
3685  /* for rpc no artificial terminal arcs should be deleted; in general: delete no solution edges */
3686  if( (rpc && !graph->mark[head])
3687  || (edgearrint[e] == CONNECT) || (edgearrint[flipedge(e)] == CONNECT) )
3688  {
3689  e = etmp;
3690  continue;
3691  }
3692 
3693  tmpcost = pathdist[k] + cost[e] + vnoi[head].dist;
3694 
3695  if( (!eliminate) || tmpcost >= minpathcost )
3696  {
3697  if( marked[flipedge(e)] )
3698  {
3699  if( eliminate )
3700  {
3701  graph_edge_del(scip, graph, e, TRUE);
3702  nfixed++;
3703  }
3704  else
3705  {
3706  if( SCIPisGT(scip, tmpcost, costrev[e]) )
3707  costrev[e] = tmpcost;
3708 
3709  if( SCIPisGT(scip, tmpcost, costrev[flipedge(e)]) )
3710  costrev[flipedge(e)] = tmpcost;
3711  }
3712  }
3713  else
3714  {
3715  marked[e] = TRUE;
3716  }
3717  }
3718  e = etmp;
3719  }
3720  }
3721  }
3722 
3723  for( k = 0; k < nnodes; k++ )
3724  {
3725  if( nfixed > minelims )
3726  break;
3727 
3728  if( !graph->mark[k] || Is_term(graph->term[k]) || solnode[k] == CONNECT )
3729  continue;
3730 
3731  if( grad[k] == 3 )
3732  {
3733  tmpcost = pathdist[k] + vnoi[k].dist + vnoi[k + nnodes].dist;
3734 
3735  if( !eliminate || tmpcost >= minpathcost )
3736  {
3737  e = graph->outbeg[k];
3738  assert(graph->oeat[e] != EAT_LAST);
3739  e2 = graph->oeat[e];
3740  assert(graph->oeat[e2] != EAT_LAST);
3741  e3 = graph->oeat[e2];
3742  assert(graph->oeat[e3] == EAT_LAST);
3743 
3744  if( SCIPisLE(scip, cost[e], 0.0) || SCIPisLE(scip, cost[e2], 0.0) || SCIPisLE(scip, cost[e3], 0.0) )
3745  continue;
3746 
3747  if( !eliminate )
3748  {
3749  for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
3750  {
3751  if( SCIPisGT(scip, tmpcost, costrev[e]) )
3752  costrev[e] = tmpcost;
3753 
3754  e2 = flipedge(e);
3755 
3756  if( SCIPisGT(scip, tmpcost, costrev[e2]) )
3757  costrev[e2] = tmpcost;
3758  }
3759 
3760  continue;
3761  }
3762  nfixed += 3;
3763  while( graph->outbeg[k] != EAT_LAST )
3764  graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
3765  }
3766  }
3767  }
3768  }
3769  SCIPdebugMessage("deleted by da: %d \n", nfixed );
3770 
3771  if( rpc )
3772  graph_pc_2trans(graph);
3773  assert(graph->mark[root]);
3774 
3775  TERMINATE:
3776  *nelims = nfixed;
3777 
3778  /* free memory */
3779  for( k = 0; k < 4; k++ )
3780  {
3781  SCIPintListNodeFree(scip, &(ancestors[k]));
3782  SCIPintListNodeFree(scip, &(revancestors[k]));
3783  }
3784 
3785  SCIPfreeBufferArray(scip, &incedge);
3786  SCIPfreeBufferArray(scip, &reinsert);
3787  SCIPfreeBufferArray(scip, &adjvert);
3788  SCIPfreeBufferArray(scip, &ecost);
3789  SCIPfreeBufferArray(scip, &revancestors);
3790  SCIPfreeBufferArray(scip, &ancestors);
3791  SCIPfreeBufferArray(scip, &sd);
3792 
3793  if( edgearrchar == NULL )
3794  SCIPfreeBufferArray(scip, &marked);
3795 
3796  return SCIP_OKAY;
3797 }
3798 
3799 
3800 /** dual ascent based reductions for PCSPG and MWCSP */
3802  SCIP* scip, /**< SCIP data structure */
3803  GRAPH* graph, /**< graph data structure */
3804  PATH* vnoi, /**< Voronoi data structure array */
3805  GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
3806  SCIP_Real* cost, /**< reduced edge costs */
3807  SCIP_Real* costrev, /**< reduced reverse edge costs */
3808  SCIP_Real* pathdist, /**< distance array for shortest path calculations */
3809  int* vbase, /**< Voronoi base array */
3810  int* pathedge, /**< shortest path incoming edge array for shortest path calculations */
3811  int* edgearrint, /**< int edges array for internal computations or NULL */
3812  int* state, /**< int 4 * vertices array */
3813  STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
3814  int* nelims, /**< pointer to store number of reduced edges */
3815  SCIP_Bool solbasedda, /**< rerun Da based on best primal solution */
3816  SCIP_Bool varyroot, /**< vary root for DA if possible */
3817  SCIP_Bool markroots, /**< should terminals proven to be part of an opt. sol. be marked as such? */
3818  SCIP_Bool userec, /**< use recombination heuristic? */
3819  SCIP_Bool fastmode, /**< run heuristic in fast mode? */
3820  SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
3821  SCIP_Real prizesum /**< sum of positive prizes */
3822 )
3823 {
3824  STPSOLPOOL* pool;
3825  GRAPH* transgraph;
3826  SCIP_Real* bestcost;
3827  SCIP_Real* edgefixingbounds;
3828  SCIP_Real* nodefixingbounds;
3829  SCIP_Real ub;
3830  SCIP_Real offset;
3831  SCIP_Real lpobjval;
3832  SCIP_Real bestlpobjval;
3833  SCIP_Real upperbound;
3834  SCIP_Real minpathcost;
3835  SCIP_Real damaxdeviation;
3836  int* roots;
3837  int* result;
3838  int* result2;
3839  int* transresult;
3840  STP_Bool* marked;
3841  int nroots;
3842  int nfixed;
3843  int nusedroots;
3844  const int root = graph->source;
3845  const int nedges = graph->edges;
3846  const int extnedges = nedges + 2 * (graph->terms - 1);
3847  const int nnodes = graph->knots;
3848  SCIP_Bool apsol;
3849  SCIP_Bool success;
3850 
3851  assert(scip != NULL);
3852  assert(graph != NULL);
3853  assert(nelims != NULL);
3854  assert(nodearrchar != NULL);
3855 
3856  /* not more than one terminal? */
3857  if( graph->terms <= 1 )
3858  return SCIP_OKAY;
3859 
3860  nroots = 0;
3861  nfixed = 0;
3862  varyroot = varyroot && userec;
3863 
3864  /* allocate memory */
3865  if( edgearrint == NULL )
3866  SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
3867  else
3868  result = edgearrint;
3869 
3870  SCIP_CALL( SCIPallocBufferArray(scip, &transresult, extnedges) );
3871  SCIP_CALL( SCIPallocBufferArray(scip, &marked, extnedges) );
3872  SCIP_CALL( SCIPallocBufferArray(scip, &result2, nedges) );
3873  SCIP_CALL( SCIPallocBufferArray(scip, &bestcost, extnedges) );
3874  SCIP_CALL( SCIPallocBufferArray(scip, &edgefixingbounds, extnedges) );
3875  SCIP_CALL( SCIPallocBufferArray(scip, &nodefixingbounds, nnodes + 1) );
3876 
3877  if( userec )
3878  SCIP_CALL( SCIPStpHeurRecInitPool(scip, &pool, nedges, SOLPOOL_SIZE) );
3879  else
3880  pool = NULL;
3881 
3882 #ifndef NDEBUG
3883  {
3884  int nterms = 0;
3885  for( int i = 0; i < nnodes; i++ )
3886  if( graph->mark[i] )
3887  if( Is_term(graph->term[i]) )
3888  nterms++;
3889 
3890  assert(nterms == (graph->terms - ((graph->stp_type != STP_RPCSPG)? 1 : 0)));
3891  }
3892 #endif
3893 
3894  /*
3895  * 1. step: compute lower bound and reduced costs
3896  */
3897 
3898  graph_pc_2trans(graph);
3899  offset = 0.0;
3900 
3901  /* transform the problem to a real SAP */
3902  SCIP_CALL( graph_pc_getSap(scip, graph, &transgraph, &offset) );
3903 
3904  /* initialize data structures for shortest paths */
3905  SCIP_CALL( graph_path_init(scip, transgraph) );
3906 
3907  damaxdeviation = fastmode ? DAMAXDEVIATION_FAST : -1.0;
3908 
3909  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
3910  gnodearr, NULL, transresult, state, root, TRUE, damaxdeviation, nodearrchar) );
3911 
3912  lpobjval += offset;
3913  bestlpobjval = lpobjval;
3914  BMScopyMemoryArray(bestcost, cost, extnedges);
3915 
3916  /* compute first primal solution */
3917  upperbound = FARAWAY;
3918  apsol = FALSE;
3919  SCIP_CALL( computeDaSolPcMw(scip, graph, NULL, vnoi, cost, pathdist, &upperbound, result, result2, vbase, pathedge, nodearrchar, &apsol) );
3920 
3921  assert(apsol && upperbound < FARAWAY);
3922  assert(graph_sol_valid(scip, graph, result));
3923 
3924  /* the required reduced path cost to be surpassed */
3925  minpathcost = upperbound - lpobjval;
3926  if( userec)
3927  SCIPdebugMessage("DA first minpathcost %f \n", minpathcost);
3928 
3929  /* initialize data structures for transgraph */
3930  SCIP_CALL( graph_init_history(scip, transgraph) );
3931  computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
3932 
3933  /*
3934  * 2. step: try to eliminate
3935  */
3936 
3937  /* restore original graph */
3938  graph_pc_2org(graph);
3939 
3940  for( int e = 0; e < extnedges; e++ )
3941  marked[e] = FALSE;
3942 
3943  /* try to reduce the graph */
3944  assert(dualCostIsValid(scip, transgraph, cost, state, nodearrchar));
3945 
3946  updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, TRUE, FALSE);
3947  updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, TRUE);
3948 
3949  nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, TRUE);
3950 
3951  /* edges from result might have been deleted! */
3952  apsol = apsol && graph_sol_unreduced(scip, graph, result);
3953  assert(!apsol || graph_sol_valid(scip, graph, result));
3954 
3955  if( userec )
3956  SCIPdebugMessage("DA: 1. NFIXED %d \n", nfixed);
3957 
3958  /* rerun dual ascent? */
3959  if( solbasedda && graph->terms > 2 && SCIPisGT(scip, minpathcost, 0.0) )
3960  {
3961  const SCIP_Real oldupperbound = upperbound;
3962 
3963  /* with recombination? */
3964  if( userec && graph->stp_type != STP_MWCSP )
3965  {
3966  /* compute second solution and add to pool */
3967  SCIP_CALL( SCIPStpHeurTMRun(scip, NULL, graph, NULL, NULL, result2, DEFAULT_HEURRUNS / 5, root, graph->cost, graph->cost, NULL, NULL, 0.0, &success, FALSE) );
3968  assert(success);
3969 
3970  SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
3971  ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
3972 
3973  SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
3974  SCIPdebugMessage("added initial TM sol to pool? %d , ub %f \n", success, ub);
3975  }
3976 
3977  /* try to improve both dual and primal bound */
3978  SCIP_CALL( computePertubedSol(scip, graph, transgraph, pool, vnoi, gnodearr, cost, costrev, bestcost, pathdist, state, vbase, pathedge, result, result2,
3979  transresult, nodearrchar, &upperbound, &lpobjval, &bestlpobjval, &minpathcost, &apsol, offset, extnedges, 0) );
3980 
3981  assert(graph_sol_valid(scip, graph, result));
3982  assert(!apsol || SCIPisEQ(scip, graph_sol_getObj(graph->cost, result, 0.0, nedges), upperbound));
3983 
3984  graph_pc_2orgcheck(graph);
3985 
3986  assert(dualCostIsValid(scip, transgraph, cost, state, nodearrchar));
3987  computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
3988  updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, FALSE);
3989  updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
3990 
3991  nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
3992 
3993  nfixed += reducePcMwTryBest(scip, graph, transgraph, vnoi, cost, costrev, bestcost, pathdist, &upperbound,
3994  &lpobjval, &bestlpobjval, &minpathcost, oldupperbound, result, vbase, state, pathedge, marked, nodearrchar, &apsol, extnedges);
3995 
3996  nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
3997  nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
3998 
3999  if( userec )
4000  SCIPdebugMessage("eliminations after sol based run2 with best dual sol %d bestlb %f newlb %f\n", nfixed, bestlpobjval, lpobjval);
4001  }
4002 
4003  /* pertubation runs for MWCSP */
4004  if( varyroot && graph->stp_type == STP_MWCSP )
4005  {
4006  for( int run = 0; run < DEFAULT_NMAXROOTS && graph->terms > STP_RED_MINBNDTERMS; run++ )
4007  {
4008  SCIP_Real oldupperbound = upperbound;
4009 
4010  graph_pc_2trans(graph);
4011 
4012  apsol = apsol && graph_sol_unreduced(scip, graph, result);
4013  assert(!apsol || graph_sol_valid(scip, graph, result));
4014 
4015  assert(SCIPisEQ(scip, upperbound, graph_sol_getObj(graph->cost, result, 0.0, nedges)));
4016 
4017  /* try to improve both dual and primal bound */
4018  SCIP_CALL( computePertubedSol(scip, graph, transgraph, pool, vnoi, gnodearr, cost, costrev, bestcost, pathdist, state, vbase, pathedge, result, result2,
4019  transresult, nodearrchar, &upperbound, &lpobjval, &bestlpobjval, &minpathcost, &apsol, offset, extnedges, run) );
4020 
4021  SCIPdebugMessage("DA: pertubated run %d ub: %f \n", run, upperbound);
4022  SCIPdebugMessage("DA: pertubated run %d minpathcost: %f \n", run, upperbound - lpobjval);
4023 
4024  computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
4025  updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, FALSE);
4026  updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
4027 
4028  nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
4029 
4030  nfixed += reducePcMwTryBest(scip, graph, transgraph, vnoi, cost, costrev, bestcost, pathdist, &upperbound,
4031  &lpobjval, &bestlpobjval, &minpathcost, oldupperbound, result, vbase, state, pathedge, marked, nodearrchar, &apsol, extnedges);
4032 
4033  nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
4034  nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
4035 
4036  SCIPdebugMessage("DA: pertubated run %d NFIXED %d \n", run, nfixed);
4037  }
4038  }
4039 
4040  if( graph->stp_type == STP_MWCSP && graph->terms < STP_RED_MINBNDTERMS )
4041  varyroot = FALSE;
4042 
4043  /* change or mark roots? */
4044  if( varyroot || markroots )
4045  {
4046  SCIP_CALL( SCIPallocBufferArray(scip, &roots, graph->terms) );
4047 
4048  findDaRoots(scip, graph, transgraph, cost, bestcost, lpobjval, bestlpobjval, upperbound, prizesum, FALSE, FALSE,
4049  vnoi, state, pathedge, vbase, nodearrchar, roots, &nroots);
4050 
4051  /* should prize of terminals be changed? */
4052  if( nroots > 0 && markroots )
4053  markPcMwRoots(scip, roots, 0, nroots, prizesum, graph, &userec, &pool);
4054 
4055  if( nroots > 0 && varyroot )
4056  SCIP_CALL( orderDaRoots(scip, graph, roots, nroots, TRUE, randnumgen) );
4057  }
4058  else
4059  {
4060  roots = NULL;
4061  }
4062 
4063  if( varyroot )
4064  nusedroots = MIN(DEFAULT_NMAXROOTS, nroots);
4065  else
4066  nusedroots = -1;
4067 
4068  graph_path_exit(scip, transgraph);
4069  graph_free(scip, &transgraph, TRUE);
4070 
4071  /* loop and change root for dual ascent run */
4072  for( int run = 0; run < nusedroots; run++ )
4073  {
4074  const int tmproot = roots[run];
4075  int transnnodes;
4076  int transnedges;
4077 
4078  assert(nroots > 0);
4079  assert(roots != NULL);
4080  assert(Is_term(graph->term[tmproot]));
4081 
4082  if( graph->terms <= 2 )
4083  break;
4084 
4085  SCIP_CALL( graph_pc_getRsap(scip, graph, &transgraph, roots, nroots, tmproot) );
4086 
4087  assert(graph_valid(transgraph));
4088 
4089  transnnodes = transgraph->knots;
4090  transnedges = transgraph->edges;
4091 
4092  for( int k = 0; k < transnnodes; k++ )
4093  transgraph->mark[k] = (transgraph->grad[k] > 0);
4094 
4095  /* init data structures for shortest paths and history */
4096  SCIP_CALL( graph_path_init(scip, transgraph) );
4097  SCIP_CALL( graph_init_history(scip, transgraph ) );
4098 
4099  transgraph->stp_type = STP_SAP;
4100 
4101  if( apsol && run > 1 )
4102  {
4103  BMScopyMemoryArray(transresult, result, graph->edges);
4104  SCIP_CALL(graph_sol_reroot(scip, transgraph, transresult, tmproot));
4105  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
4106  gnodearr, transresult, result2, state, tmproot, FALSE, -1.0, nodearrchar));
4107  }
4108  else
4109  {
4110  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
4111  gnodearr, NULL, transresult, state, tmproot, FALSE, -1.0, nodearrchar));
4112  }
4113 
4114  assert(graph_valid(transgraph));
4115 
4116  for( int e = graph->outbeg[tmproot]; e != EAT_LAST; e = graph->oeat[e] )
4117  {
4118  const int k = graph->head[e];
4119  if( Is_term(graph->term[k]) )
4120  {
4121  if( k == root )
4122  cost[flipedge(e)] = 0.0;
4123  else
4124  cost[e] = 0.0;
4125  }
4126  }
4127 
4128  for( int k = 0; k < nnodes; k++ )
4129  {
4130  if( Is_pterm(graph->term[k]) && graph->prize[k] >= prizesum )
4131  {
4132  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
4133  {
4134  const int head = graph->head[e];
4135  if( Is_term(graph->term[head]) && head != root )
4136  cost[e] = 0.0;
4137  }
4138  }
4139  }
4140 
4141  apsol = FALSE;
4142  SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, &upperbound, result, result2, vbase, pathedge, nodearrchar, &apsol) );
4143 
4144  SCIPdebugMessage("ROOTRUNS upperbound %f \n", upperbound);
4145  SCIPdebugMessage("ROOTRUNS best sol in pool %f \n", pool->sols[0]->obj);
4146 
4147  for( int k = 0; k < transnnodes; k++ )
4148  transgraph->mark[k] = (transgraph->grad[k] > 0);
4149 
4150  /* the required reduced path cost to be surpassed */
4151  minpathcost = upperbound - lpobjval;
4152 
4153  if( markroots )
4154  {
4155  const int oldnroots = nroots;
4156  findDaRoots(scip, graph, transgraph, cost, cost, lpobjval, lpobjval, upperbound, prizesum, TRUE, TRUE,
4157  vnoi, state, pathedge, vbase, nodearrchar, roots, &nroots);
4158 
4159  /* should prize of terminals be changed? */
4160  if( nroots > oldnroots )
4161  markPcMwRoots(scip, roots, oldnroots, nroots, prizesum, graph, &userec, &pool);
4162  }
4163 
4164  SCIPdebugMessage("ROOTRUNS: minpathcost %f \n", minpathcost);
4165  SCIPdebugMessage("lb: %f ub: %f \n", lpobjval, upperbound);
4166 
4167  /* distance from root to all nodes */
4168  graph_path_execX(scip, transgraph, tmproot, cost, pathdist, pathedge);
4169 
4170  for( int e = 0; e < transnedges; e++ )
4171  costrev[e] = cost[flipedge(e)];
4172 
4173  /* no paths should go back to the root */
4174  for( int e = transgraph->outbeg[tmproot]; e != EAT_LAST; e = transgraph->oeat[e] )
4175  costrev[e] = FARAWAY;
4176 
4177  /* build Voronoi diagram */
4178  graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, state);
4179  graph_get2next(scip, transgraph, costrev, costrev, vnoi, vbase, transgraph->path_heap, state);
4180 
4181  /* restore original graph */
4182  graph_pc_2org(graph);
4183 
4184  assert(graph->mark[tmproot]);
4185  graph->mark[tmproot] = FALSE;
4186 
4187  /* at first run switch to undirected case */
4188  if( run == 0 )
4189  for( int e = 0; e < extnedges; e++ )
4190  edgefixingbounds[e] = MIN(edgefixingbounds[e], edgefixingbounds[flipedge(e)]);
4191 
4192  updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, TRUE);
4193  updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
4194 
4195  for( int e = 0; e < transnedges; e++ )
4196  marked[e] = FALSE;
4197 
4198  /* try to eliminate vertices and edges */
4199  nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
4200  nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
4201  nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
4202 
4203  apsol = apsol && graph_sol_unreduced(scip, graph, result);
4204  assert(!apsol || graph_sol_valid(scip, graph, result));
4205  SCIPdebugMessage("FIXED with changed root %d \n\n", nfixed);
4206 
4207  graph->mark[tmproot] = TRUE;
4208 
4209  transgraph->stp_type = STP_RPCSPG;
4210  graph_path_exit(scip, transgraph);
4211  graph_free(scip, &transgraph, TRUE);
4212  }
4213 
4214  *nelims = nfixed;
4215 
4216  /* free memory */
4217  if( pool != NULL )
4218  SCIPStpHeurRecFreePool(scip, &pool);
4219  SCIPfreeBufferArrayNull(scip, &roots);
4220  SCIPfreeBufferArray(scip, &nodefixingbounds);
4221  SCIPfreeBufferArray(scip, &edgefixingbounds);
4222  SCIPfreeBufferArray(scip, &bestcost);
4223  SCIPfreeBufferArray(scip, &result2);
4224  SCIPfreeBufferArray(scip, &marked);
4225  SCIPfreeBufferArray(scip, &transresult);
4226 
4227  if( edgearrint == NULL )
4228  SCIPfreeBufferArray(scip, &result);
4229 
4230  assert(graph_valid(graph));
4231 
4232  return SCIP_OKAY;
4233 }
4234 
4235 
4236 /** dual ascent based heuristic reductions for MWCSP */
4238  SCIP* scip, /**< SCIP data structure */
4239  GRAPH* graph, /**< graph data structure */
4240  PATH* vnoi, /**< Voronoi data structure array */
4241  GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
4242  SCIP_Real* cost, /**< edge costs */
4243  SCIP_Real* costrev, /**< reverse edge costs */
4244  SCIP_Real* pathdist, /**< distance array for shortest path calculations */
4245  int* vbase, /**< Voronoi base array */
4246  int* pathedge, /**< shortest path incoming edge array for shortest path calculations */
4247  int* soledge, /**< edge solution array (CONNECT/UNKNOWN) or NULL; needs to contain solution if solgiven == TRUE */
4248  int* state, /**< int 4 * vertices array */
4249  int* solnode, /**< array of nodes of current solution that is not to be destroyed */
4250  STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
4251  int* nelims, /**< pointer to store number of reduced edges */
4252  int minelims, /**< minimum number of edges to eliminate */
4253  SCIP_Bool solgiven /**< solution provided? */
4254  )
4255 {
4256  SCIP_HEURDATA* tmheurdata;
4257  IDX** ancestors;
4258  IDX** revancestors;
4259  GRAPH* transgraph;
4260  SCIP_Real* sd;
4261  SCIP_Real* ecost;
4262  SCIP_Real ub;
4263  SCIP_Real offset;
4264  SCIP_Bool success;
4265  SCIP_Real tmpcost;
4266  SCIP_Real lpobjval;
4267  SCIP_Real hopfactor;
4268  SCIP_Real upperbound;
4269  SCIP_Real minpathcost;
4270  SCIP_Bool eliminate;
4271  int* result;
4272  int* adjvert;
4273  int* incedge;
4274  int* reinsert;
4275  int* transresult;
4276  int i;
4277  int k;
4278  int e;
4279  int e2;
4280  int etmp;
4281  int runs;
4282  int root;
4283  int nterms;
4284  int nedges;
4285  int nnodes;
4286  int nfixed;
4287  int redrounds;
4288  int best_start;
4289  int transnnodes;
4290  int transnedges;
4291  STP_Bool* marked;
4292 
4293  assert(scip != NULL);
4294  assert(graph != NULL);
4295  assert(nelims != NULL);
4296  assert(nodearrchar != NULL);
4297 
4298  root = graph->source;
4299  nfixed = 0;
4300  nnodes = graph->knots;
4301  nedges = graph->edges;
4302  success = FALSE;
4303  hopfactor = DEFAULT_HOPFACTOR;
4304 
4305  /* not more than two terminals? */
4306  if( graph->terms <= 3 )
4307  return SCIP_OKAY;
4308 
4309  /* allocate memory */
4310  if( soledge == NULL )
4311  {
4312  SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
4313  }
4314  else
4315  {
4316  result = soledge;
4317  }
4318 
4319  SCIP_CALL( SCIPallocBufferArray(scip, &transresult, nedges + 2 * (graph->terms - 1)) );
4320  SCIP_CALL( SCIPallocBufferArray(scip, &marked, nedges + 2 * (graph->terms - 1)) );
4321 
4322  /* allocate length-4 buffer memory */
4323  SCIP_CALL( SCIPallocBufferArray(scip, &sd, 4) );
4324  SCIP_CALL( SCIPallocBufferArray(scip, &ancestors, 4) );
4325  SCIP_CALL( SCIPallocBufferArray(scip, &revancestors, 4) );
4326  SCIP_CALL( SCIPallocBufferArray(scip, &ecost, 4) );
4327  SCIP_CALL( SCIPallocBufferArray(scip, &adjvert, 4) );
4328  SCIP_CALL( SCIPallocBufferArray(scip, &reinsert, 4) );
4329  SCIP_CALL( SCIPallocBufferArray(scip, &incedge, 4) );
4330 
4331  for( i = 0; i < 4; i++ )
4332  {
4333  sd[i] = 0.0;
4334  ancestors[i] = NULL;
4335  revancestors[i] = NULL;
4336  }
4337 
4338  /* 1. step: compute upper bound */
4339 
4340  runs = 0;
4341  nterms = 0;
4342  for( i = 0; i < nnodes; i++ )
4343  {
4344  if( graph->grad[i] > 0 )
4345  {
4346  runs++;
4347  if( Is_term(graph->term[i]) )
4348  nterms++;
4349  }
4350  }
4351 
4352  assert(nterms == (graph->terms - ((graph->stp_type != STP_RPCSPG)? 1 : 0)));
4353 
4354  offset = 0.0;
4355 
4356  /* transform the problem to a real SAP */
4357  SCIP_CALL( graph_pc_getSap(scip, graph, &transgraph, &offset) );
4358 
4359  /* initialize data structures for shortest paths and history */
4360  SCIP_CALL( graph_path_init(scip, transgraph) );
4361  SCIP_CALL( graph_init_history(scip, transgraph ) );
4362  transnnodes = transgraph->knots;
4363  transnedges = transgraph->edges;
4364 
4365  if( !solgiven )
4366  {
4367  for( k = 0; k < nnodes; k++ )
4368  nodearrchar[k] = (solnode[k] == CONNECT);
4369 
4370  for( e = 0; e < nedges; e++ )
4371  {
4372  cost[e] = graph->cost[e];
4373  costrev[e] = graph->cost[flipedge(e)];
4374  result[e] = UNKNOWN;
4375  }
4376 
4377  /* build trivial solution */
4378  SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, result, nodearrchar) );
4379  }
4380  else
4381  {
4382  for( e = 0; e < nedges; e++ )
4383  {
4384  cost[e] = graph->cost[e];
4385  costrev[e] = graph->cost[flipedge(e)];
4386  }
4387  }
4388 
4389  upperbound = graph_sol_getObj(graph->cost, result, 0.0, nedges);
4390 
4391  /* number of runs should not exceed number of connected vertices */
4392  runs = MIN(runs, DEFAULT_HEURRUNS);
4393 
4394  /* get TM heuristic data */
4395  assert(SCIPfindHeur(scip, "TM") != NULL);
4396  tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
4397 
4398 
4399  /* compute Steiner tree to obtain upper bound */
4400  SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, graph, NULL, &best_start, transresult, runs, root, cost, costrev, &hopfactor, NULL, 0.0, &success, FALSE) );
4401 
4402  /* feasible solution found? */
4403  if( success )
4404  {
4405  /* calculate objective value of solution */
4406  ub = graph_sol_getObj(graph->cost, transresult, 0.0, nedges);
4407 
4408  SCIPdebugMessage("TM upperbound in reduce_daSlackPruneMw %f \n", ub + SCIPprobdataGetOffset(scip));
4409 
4410  if( SCIPisLE(scip, ub, upperbound) )
4411  {
4412  upperbound = ub;
4413  for( e = 0; e < nedges; e++ )
4414  result[e] = transresult[e];
4415  }
4416  }
4417 
4418  /* compute lower bound and reduced costs todo use SCIPdualAscentStpSol */
4419  SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, transresult, state, root, FALSE, -1.0, nodearrchar) );
4420 
4421  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, transresult, vbase, -1, nodearrchar, &success, FALSE) );
4422 
4423  assert(success);
4424  assert(graph_sol_valid(scip, graph, transresult));
4425 
4426  if( success )
4427  {
4428  ub = graph_sol_getObj(graph->cost, transresult, 0.0, nedges);
4429 
4430  SCIPdebugMessage("AP upperbound in reduce_daSlackPruneMw %f \n", ub + SCIPprobdataGetOffset(scip));
4431 
4432  if( SCIPisLE(scip, ub, upperbound) )
4433  for( e = 0; e < nedges; e++ )
4434  result[e] = transresult[e];
4435  }
4436 
4437  /*
4438  * 2. step: try to eliminate
4439  * */
4440 
4441  for( e = 0; e < transnedges; e++ )
4442  {
4443  costrev[e] = cost[flipedge(e)];
4444  marked[e] = FALSE;
4445  }
4446 
4447  for( k = 0; k < transnnodes; k++ )
4448  transgraph->mark[k] = (transgraph->grad[k] > 0);
4449 
4450  /* distance from root to all nodes */
4451  graph_path_execX(scip, transgraph, root, cost, pathdist, pathedge);
4452 
4453  for( i = 0; i < transnnodes; i++ )
4454  if( Is_term(transgraph->term[i]) )
4455  assert(SCIPisEQ(scip, pathdist[i], 0.0));
4456 
4457  /* no paths should go back to the root */
4458  for( e = transgraph->outbeg[root]; e != EAT_LAST; e = transgraph->oeat[e] )
4459  costrev[e] = FARAWAY;
4460 
4461  /* build Voronoi diagram */
4462  graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, transgraph->path_state);
4463 
4464  /* restore original graph */
4465  graph_pc_2org(graph);
4466 
4467  for( k = 0; k < nnodes; k++ )
4468  solnode[k] = UNKNOWN;
4469 
4470  for( e = 0; e < nedges; e++ )
4471  {
4472  costrev[e] = -1.0;
4473  if( result[e] == CONNECT )
4474  {
4475  solnode[graph->head[e]] = CONNECT;
4476  solnode[graph->tail[e]] = CONNECT;
4477  }
4478  }
4479 
4480  for( redrounds = 0; redrounds < 3; redrounds++ )
4481  {
4482  if( redrounds == 0 )
4483  {
4484  eliminate = FALSE;
4485  minpathcost = FARAWAY;
4486  }
4487  else if( redrounds == 1 )
4488  {
4489  assert(minelims > 0);
4490  assert(2 * minelims < nedges);
4491 
4492  eliminate = TRUE;
4493  SCIPsortReal(costrev, nedges);
4494 
4495  /* the required reduced path cost to be surpassed */
4496  minpathcost = costrev[nedges - 2 * minelims];
4497 
4498  if( SCIPisLE(scip, minpathcost, 0.0) )
4499  minpathcost = 2 * SCIPepsilon(scip);
4500 
4501  k = nedges - 2 * minelims;
4502 
4503  /* try to first eliminate edges with higher gap */
4504  for( e = nedges - 1; e > k && e >= 2; e = e - 2 )
4505  {
4506  if( SCIPisLE(scip, costrev[e - 2], minpathcost) )
4507  break;
4508  }
4509 
4510  if( SCIPisGT(scip, costrev[e], minpathcost) )
4511  minpathcost = costrev[e];
4512 
4513  if( SCIPisLE(scip, minpathcost, 0.0) )
4514  minpathcost = 2 * SCIPepsilon(scip);
4515 
4516 
4517  for( e = 0; e < nedges; e++ )
4518  marked[e] = FALSE;
4519  }
4520  else
4521  {
4522  eliminate = TRUE;
4523 
4524  /* the required reduced path cost to be surpassed */
4525  minpathcost = costrev[nedges - 2 * minelims];
4526 
4527  if( SCIPisLE(scip, minpathcost, 0.0) )
4528  minpathcost = 2 * SCIPepsilon(scip);
4529 
4530  for( e = 0; e < nedges; e++ )
4531  marked[e] = FALSE;
4532  }
4533 
4534  /* try to eliminate vertices and edges */
4535  for( k = 0; k < nnodes; k++ )
4536  {
4537  if( !graph->mark[k] )
4538  continue;
4539 
4540  if( nfixed > minelims )
4541  break;
4542 
4543 
4544  if( Is_term(graph->term[k]) )
4545  continue;
4546 
4547  tmpcost = pathdist[k] + vnoi[k].dist;
4548 
4549  if( SCIPisGT(scip, tmpcost, minpathcost) ||
4550  (SCIPisGE(scip, tmpcost, minpathcost) && solnode[k] != CONNECT) || (!eliminate && solnode[k] != CONNECT) )
4551  {
4552  if( !eliminate )
4553  {
4554  for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
4555  {
4556  if( SCIPisGT(scip, tmpcost, costrev[e]) )
4557  costrev[e] = tmpcost;
4558 
4559  e2 = flipedge(e);
4560 
4561  if( SCIPisGT(scip, tmpcost, costrev[e2]) )
4562  costrev[e2] = tmpcost;
4563  }
4564 
4565  continue;
4566  }
4567 
4568  while( transgraph->outbeg[k] != EAT_LAST )
4569  {
4570  e = transgraph->outbeg[k];
4571  graph_edge_del(scip, transgraph, e, FALSE);
4572  graph_edge_del(scip, graph, e, TRUE);
4573  nfixed++;
4574  }
4575  assert(graph->outbeg[k] == EAT_LAST);
4576  }
4577  else
4578  {
4579  e = transgraph->outbeg[k];
4580  while( e != EAT_LAST )
4581  {
4582  etmp = transgraph->oeat[e];
4583  tmpcost = pathdist[k] + cost[e] + vnoi[transgraph->head[e]].dist;
4584 
4585  if( SCIPisGT(scip, tmpcost, minpathcost) ||
4586  ( (SCIPisGE(scip, tmpcost, minpathcost) || !eliminate) && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
4587  {
4588  if( marked[flipedge(e)] )
4589  {
4590  if( eliminate )
4591  {
4592  graph_edge_del(scip, graph, e, TRUE);
4593  graph_edge_del(scip, transgraph, e, FALSE);
4594  nfixed++;
4595  }
4596  else
4597  {
4598  if( SCIPisGT(scip, tmpcost, costrev[e]) )
4599  costrev[e] = tmpcost;
4600 
4601  if( SCIPisGT(scip, tmpcost, costrev[flipedge(e)]) )
4602  costrev[flipedge(e)] = tmpcost;
4603  }
4604  }
4605  else
4606  {
4607  marked[e] = TRUE;
4608  }
4609  }
4610  e = etmp;
4611  }
4612  }
4613  }
4614  }
4615  assert(root == transgraph->source);
4616 
4617  SCIPdebugMessage("fixed by da: %d \n", nfixed );
4618 
4619  graph_path_exit(scip, transgraph);
4620  graph_free(scip, &transgraph, TRUE);
4621 
4622  /* restore transformed graph */
4623  graph_pc_2trans(graph);
4624 
4625  *nelims = nfixed;
4626 
4627  /* free memory */
4628  for( k = 0; k < 4; k++ )
4629  {
4630  SCIPintListNodeFree(scip, &(ancestors[k]));
4631  SCIPintListNodeFree(scip, &(revancestors[k]));
4632  }
4633 
4634  SCIPfreeBufferArray(scip, &incedge);
4635  SCIPfreeBufferArray(scip, &reinsert);
4636  SCIPfreeBufferArray(scip, &adjvert);
4637  SCIPfreeBufferArray(scip, &ecost);
4638  SCIPfreeBufferArray(scip, &revancestors);
4639  SCIPfreeBufferArray(scip, &ancestors);
4640  SCIPfreeBufferArray(scip, &sd);
4641  SCIPfreeBufferArray(scip, &marked);
4642  SCIPfreeBufferArray(scip, &transresult);
4643 
4644  if( soledge == NULL )
4645  SCIPfreeBufferArray(scip, &result);
4646 
4647  return SCIP_OKAY;
4648 }
4649 
4650 
4651 
4652 /** bound-based reductions for the (R)PCSTP, the MWCSP and the STP */
4654  SCIP* scip, /**< SCIP data structure */
4655  GRAPH* graph, /**< graph data structure */
4656  PATH* vnoi, /**< Voronoi data structure */
4657  SCIP_Real* cost, /**< edge cost array */
4658  SCIP_Real* prize, /**< prize (nodes) array */
4659  SCIP_Real* radius, /**< radius array */
4660  SCIP_Real* costrev, /**< reversed edge cost array */
4661  SCIP_Real* offset, /**< pointer to the offset */
4662  SCIP_Real* upperbound, /**< pointer to an upper bound */
4663  int* heap, /**< heap array */
4664  int* state, /**< array to store state of a node during Voronoi computation*/
4665  int* vbase, /**< Voronoi base to each node */
4666  int* nelims /**< pointer to store number of eliminated edges */
4667  )
4668 {
4669  SCIP_HEURDATA* tmheurdata;
4670  GRAPH* adjgraph;
4671  PATH* mst;
4672  SCIP_Real r;
4673  SCIP_Real obj;
4674  SCIP_Real max;
4675  SCIP_Real radii;
4676  SCIP_Real bound;
4677  SCIP_Real tmpcost;
4678  SCIP_Real mstobj;
4679  SCIP_Real mstobj2;
4680  SCIP_Real maxcost;
4681  SCIP_Real radiim2;
4682  SCIP_Real radiim3;
4683  int* perm;
4684  int* result;
4685  int* starts;
4686  int e;
4687  int k;
4688  int l;
4689  int head;
4690  int tail;
4691  int runs;
4692  int root;
4693  int etemp;
4694  int nterms;
4695  int nnodes;
4696  int nedges;
4697  int best_start;
4698  STP_Bool* stnode;
4699  SCIP_Bool ub;
4700  SCIP_Bool pc;
4701  SCIP_Bool mw;
4702  SCIP_Bool pcmw;
4703  SCIP_Bool success;
4704 
4705  assert(scip != NULL);
4706  assert(graph != NULL);
4707  assert(vnoi != NULL);
4708  assert(cost != NULL);
4709  assert(radius != NULL);
4710  assert(costrev != NULL);
4711  assert(heap != NULL);
4712  assert(state != NULL);
4713  assert(vbase != NULL);
4714  assert(nelims != NULL);
4715  assert(graph->source >= 0);
4716  assert(upperbound != NULL);
4717 
4718  mst = NULL;
4719  obj = DEFAULT_HOPFACTOR;
4720  perm = NULL;
4721  root = graph->source;
4722  nedges = graph->edges;
4723  nnodes = graph->knots;
4724  mstobj = 0.0;
4725  *nelims = 0;
4726  mstobj2 = 0.0;
4727  best_start = 0;
4728  ub = SCIPisGT(scip, *upperbound, 0.0);
4729  mw = (graph->stp_type == STP_MWCSP);
4730  pc = (graph->stp_type == STP_RPCSPG) || (graph->stp_type == STP_PCSPG);
4731  pcmw = pc || mw;
4732  success = TRUE;
4733 
4734  /* no upper bound provided? */
4735  if( !ub )
4736  {
4737  /* allocate memory */
4738  SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
4739  SCIP_CALL( SCIPallocBufferArray(scip, &stnode, nnodes) );
4740  }
4741  else
4742  {
4743  result = NULL;
4744  stnode = NULL;
4745  }
4746 
4747  /* initialize */
4748  e = 0;
4749  nterms = 0;
4750  for( k = 0; k < nnodes; k++ )
4751  {
4752  if( !ub )
4753  {
4754  assert(stnode != NULL);
4755  stnode[k] = FALSE;
4756  }
4757  if( !pcmw )
4758  graph->mark[k] = (graph->grad[k] > 0);
4759  if( graph->mark[k] )
4760  {
4761  e++;
4762  if( Is_term(graph->term[k]) )
4763  nterms++;
4764  }
4765  }
4766 
4767  /* not more than two terminals? */
4768  if( nterms <= 2 || (pcmw && nterms <= 3) )
4769  {
4770  /* free memory and return */
4771  SCIPfreeBufferArrayNull(scip, &stnode);
4772  SCIPfreeBufferArrayNull(scip, &result);
4773  return SCIP_OKAY;
4774  }
4775 
4776  assert(nterms == (graph->terms - ((graph->stp_type == STP_PCSPG || mw)? 1 : 0)));
4777 
4778  runs = MIN(e, DEFAULT_HEURRUNS);
4779 
4780  /* neither PC, MW, RPC nor HC? */
4781  if( !pcmw && graph->stp_type != STP_DHCSTP )
4782  {
4783  /* choose starting points for TM heuristic */
4784 
4785  SCIP_CALL( SCIPallocBufferArray(scip, &starts, nnodes) );
4786 
4787  SCIPStpHeurTMCompStarts(graph, starts, &runs);
4788  }
4789  else
4790  {
4791  starts = NULL;
4792  }
4793 
4794  /* initialise cost and costrev array */
4795  maxcost = 0.0;
4796  for( e = 0; e < nedges; e++ )
4797  {
4798  if( !ub )
4799  {
4800  assert(result != NULL);
4801  result[e] = UNKNOWN;
4802  }
4803  cost[e] = graph->cost[e];
4804  costrev[e] = graph->cost[flipedge(e)];
4805 
4806  assert(SCIPisGE(scip, cost[e], 0.0));
4807 
4808  if( graph->stp_type == STP_DHCSTP && SCIPisLT(scip, graph->cost[e], BLOCKED) && SCIPisGT(scip, graph->cost[e], maxcost) )
4809  maxcost = graph->cost[e];
4810  }
4811 
4812  /* init auxiliary graph */
4813  if( !mw )
4814  {
4815  SCIP_CALL( graph_init(scip, &adjgraph, nterms, MIN(nedges, (nterms - 1) * nterms), 1) );
4816  }
4817  else
4818  {
4819  adjgraph = NULL;
4820  }
4821 
4822  /* build voronoi regions, concomitantly building adjgraph and computing radii values*/
4823  SCIP_CALL( graph_voronoiWithRadius(scip, graph, adjgraph, vnoi, radius, cost, costrev, vbase, heap, state) );
4824 
4825  /* get 2nd next terminals to all non-terminal nodes */
4826  graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4827 
4828  /* get 3th next terminals to all non-terminal nodes */
4829  graph_get3next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4830 
4831  /* for (rooted) prize collecting get 4th next terminals to all non-terminal nodes */
4832  if( pc )
4833  graph_get4next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4834 
4835  /* no MWCS problem? */
4836  if( !mw )
4837  {
4838  assert(adjgraph != NULL);
4839  graph_knot_chg(adjgraph, 0, 0);
4840  adjgraph->source = 0;
4841 
4842  /* compute MST on adjgraph */
4843  SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
4844  SCIP_CALL( graph_path_init(scip, adjgraph) );
4845  graph_path_exec(scip, adjgraph, MST_MODE, 0, adjgraph->cost, mst);
4846 
4847  max = -1.0;
4848  r = -1.0;
4849  for( k = 1; k < nterms; k++ )
4850  {
4851  assert(adjgraph->path_state[k] == CONNECT);
4852  e = mst[k].edge;
4853  assert(e >= 0);
4854  tmpcost = adjgraph->cost[e];
4855  mstobj += tmpcost;
4856  if( SCIPisGT(scip, tmpcost, max) )
4857  max = tmpcost;
4858  else if( SCIPisGT(scip, tmpcost, r) )
4859  r = tmpcost;
4860  }
4861  mstobj -= max;
4862  mstobj2 = mstobj - r;
4863  }
4864 
4865  /* for (rooted) prize collecting an maximum weight problems: correct radius values */
4866  if( graph->stp_type == STP_RPCSPG )
4867  {
4868  assert(graph->mark[graph->source]);
4869  for( k = 0; k < nnodes; k++ )
4870  {
4871  if( !Is_term(graph->term[k]) || !graph->mark[k] )
4872  continue;
4873 
4874  if( SCIPisGT(scip, radius[k], prize[k]) && k != graph->source )
4875  radius[k] = prize[k];
4876  }
4877  }
4878  else if( graph->stp_type == STP_PCSPG )
4879  {
4880  for( k = 0; k < nnodes; k++ )
4881  {
4882  if( !graph->mark[k] )
4883  continue;
4884 
4885  if( Is_term(graph->term[k]) && SCIPisGT(scip, radius[k], prize[k]) )
4886  radius[k] = prize[k];
4887  }
4888  }
4889  else if( graph->stp_type == STP_MWCSP )
4890  {
4891  max = 0.0;
4892  for( k = 0; k < nnodes; k++ )
4893  {
4894  if( !Is_term(graph->term[k]) || !graph->mark[k] )
4895  continue;
4896 
4897  assert(SCIPisGE(scip, prize[k], 0.0));
4898  if( SCIPisGT(scip, prize[k], max) )
4899  max = prize[k];
4900 
4901  if( SCIPisGE(scip, radius[k], 0.0 ) )
4902  radius[k] = 0.0;
4903  else
4904  radius[k] = -radius[k];
4905  }
4906  }
4907 
4908  /* sort radius values */
4909  if( pc )
4910  {
4911  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nnodes) );
4912  for( k = 0; k < nnodes; k++ )
4913  perm[k] = k;
4914  SCIPsortRealInt(radius, perm, nnodes);
4915  }
4916  else
4917  {
4918  SCIPsortReal(radius, nnodes);
4919  }
4920 
4921  radiim2 = 0.0;
4922 
4923  /* sum all but two radius values of highest/lowest value */
4924  if( mw )
4925  {
4926  for( k = 2; k < nterms; k++ )
4927  {
4928  assert( SCIPisGT(scip, FARAWAY, radius[k]) );
4929  radiim2 += radius[k];
4930  }
4931  }
4932  else
4933  {
4934  for( k = 0; k < nterms - 2; k++ )
4935  {
4936  assert( SCIPisGT(scip, FARAWAY, radius[k]) );
4937  radiim2 += radius[k];
4938  }
4939  }
4940  radii = radiim2 + radius[nterms - 2] + radius[nterms - 1];
4941 #if 1
4942  if( nterms >= 3 )
4943  radiim3 = radiim2 - radius[nterms - 3];
4944  else
4945  radiim3 = 0;
4946 #endif
4947  /* get TM heuristic data */
4948  assert(SCIPfindHeur(scip, "TM") != NULL);
4949  tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
4950 
4951  /* no upper bound available? */
4952  if( !ub )
4953  {
4954  assert(stnode != NULL);
4955  assert(result != NULL);
4956 
4957  /* PC or RPC? Then restore transformed graph */
4958  if( pcmw )
4959  graph_pc_2trans(graph);
4960 
4961  SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, graph, starts, &best_start, result, runs, root, cost, costrev, &obj, NULL, maxcost, &success, FALSE) );
4962 
4963  /* PC or RPC? Then restore oringinal graph */
4964  if( pcmw )
4965  graph_pc_2org(graph);
4966 
4967  /* no feasible solution found? */
4968  if( !success )
4969  {
4970  /* free memory and return */
4971  if( !mw )
4972  {
4973  graph_path_exit(scip, adjgraph);
4974  graph_free(scip, &adjgraph, TRUE);
4975  }
4976  SCIPfreeBufferArrayNull(scip, &mst);
4977  SCIPfreeBufferArrayNull(scip, &starts);
4978  SCIPfreeBufferArray(scip, &result);
4979  SCIPfreeBufferArray(scip, &stnode);
4980  return SCIP_OKAY;
4981  }
4982 
4983  /* calculate objective value of solution */
4984  obj = 0.0;
4985  for( e = 0; e < nedges; e++ )
4986  {
4987  if( result[e] == CONNECT )
4988  {
4989  head = graph->head[e];
4990  if( mw )
4991  {
4992  if( graph->mark[head] )
4993  {
4994  assert(stnode[head] == FALSE);
4995  obj += prize[head];
4996  }
4997  }
4998  else
4999  {
5000  obj += graph->cost[e];
5001  stnode[head] = TRUE;
5002  stnode[graph->tail[e]] = TRUE;
5003  }
5004  }
5005  }
5006  *upperbound = obj;
5007  }
5008  else
5009  {
5010  obj = *upperbound;
5011  assert(SCIPisGE(scip, obj, 0.0));
5012  }
5013 
5014  if( SCIPisGT(scip, radiim2, mstobj) )
5015  bound = radiim2;
5016  else
5017  bound = mstobj;
5018 
5019  /* traverse all node, try to eliminate each node or incident edges */
5020  for( k = 0; k < nnodes; k++ )
5021  {
5022  if( (!graph->mark[k] && (pcmw)) || graph->grad[k] == 0 )
5023  continue;
5024 
5025  if( mw )
5026  {
5027  tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + bound - graph->prize[k];
5028 
5029  if( !Is_term(graph->term[k]) && (SCIPisLT(scip, tmpcost, obj)) )
5030  {
5031  while( graph->outbeg[k] != EAT_LAST )
5032  {
5033  (*nelims)++;
5034  graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5035  }
5036  }
5037  else
5038  {
5039  e = graph->outbeg[k];
5040  while( e != EAT_LAST )
5041  {
5042  etemp = graph->oeat[e];
5043  tail = graph->tail[e];
5044  head = graph->head[e];
5045  if( !graph->mark[head] )
5046  {
5047  e = etemp;
5048  continue;
5049  }
5050  tmpcost = bound - graph->prize[k];
5051  if( vbase[tail] != vbase[head] )
5052  {
5053  tmpcost -= vnoi[head].dist + vnoi[tail].dist;
5054  }
5055  else
5056  {
5057  if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5058  tmpcost -= vnoi[tail + nnodes].dist + vnoi[head].dist;
5059  else
5060  tmpcost -= vnoi[tail].dist + vnoi[head + nnodes].dist;
5061  }
5062  if( (SCIPisLT(scip, tmpcost, obj)) )
5063  {
5064  (*nelims)++;
5065  graph_edge_del(scip, graph, e, TRUE);
5066  }
5067  e = etemp;
5068  }
5069  }
5070  continue;
5071  }
5072 
5073  tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + bound;
5074 
5075  /* can node k be deleted? */
5076  if( !Is_term(graph->term[k]) && (SCIPisGT(scip, tmpcost, obj)
5077  || (((stnode != NULL) ? !stnode[k] : 0) && SCIPisGE(scip, tmpcost, obj))) )
5078  {
5079  /* delete all incident edges */
5080  while( graph->outbeg[k] != EAT_LAST )
5081  {
5082  e = graph->outbeg[k];
5083  (*nelims)++;
5084  assert(!pc || graph->tail[e] != root);
5085  assert(!pc || graph->mark[graph->head[e]]);
5086  assert(!Is_pterm(graph->term[graph->head[e]]));
5087  assert(!Is_pterm(graph->term[graph->tail[e]]));
5088 
5089  graph_edge_del(scip, graph, e, TRUE);
5090  }
5091  }
5092  else
5093  {
5094  e = graph->outbeg[k];
5095  while( e != EAT_LAST )
5096  {
5097  etemp = graph->oeat[e];
5098  tail = graph->tail[e];
5099  head = graph->head[e];
5100  tmpcost = graph->cost[e] + bound;
5101 
5102  if( vbase[tail] != vbase[head] )
5103  {
5104  tmpcost += vnoi[head].dist + vnoi[tail].dist;
5105  }
5106  else
5107  {
5108  if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5109  tmpcost += vnoi[tail + nnodes].dist + vnoi[head].dist;
5110  else
5111  tmpcost += vnoi[tail].dist + vnoi[head + nnodes].dist;
5112  assert(SCIPisGE(scip, tmpcost, vnoi[head].dist + vnoi[tail].dist + graph->cost[e] + bound));
5113  }
5114  /* can edge e or arc e be deleted? */
5115  if( (SCIPisGT(scip, tmpcost, obj) || (((result != NULL) ? (result[e] != CONNECT) : 0) && result[flipedge(e)] != CONNECT && SCIPisGE(scip, tmpcost, obj)))
5116  && SCIPisLT(scip, graph->cost[e], FARAWAY) && (!(pc || mw) || graph->mark[head]) )
5117  {
5118  if( graph->stp_type == STP_DHCSTP && SCIPisGT(scip, graph->cost[e], graph->cost[flipedge(e)]) )
5119  {
5120  graph->cost[e] = FARAWAY;
5121  (*nelims)++;
5122  }
5123  else
5124  {
5125  assert(!Is_pterm(graph->term[head]));
5126  assert(!Is_pterm(graph->term[tail]));
5127  graph_edge_del(scip, graph, e, TRUE);
5128  (*nelims)++;
5129  }
5130  }
5131  e = etemp;
5132  }
5133  }
5134  }
5135 #if 1
5136  /* traverse all node, try to eliminate 3 degree nodes */
5137  if( !mw )
5138  {
5139  for( k = 0; k < nnodes; k++ )
5140  {
5141  if( (!graph->mark[k] && pc) || graph->grad[k] == 0 )
5142  continue;
5143 
5144  if( (graph->grad[k] == 3 || graph->grad[k] == 4) && !Is_term(graph->term[k]) )
5145  {
5146  tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5147  if( SCIPisGT(scip, tmpcost, obj) )
5148  {
5149  SCIP_CALL( graph_knot_delPseudo(scip, graph, graph->cost, NULL, NULL, k, &success) );
5150 
5151  assert(graph->grad[k] == 0 || (graph->grad[k] == 4 && !success));
5152  }
5153  }
5154  }
5155  }
5156 #endif
5157 
5158  /* for(R)PC: try to eliminate terminals */
5159  if( pc )
5160  {
5161  SCIP_CALL( graph_get4nextTTerms(scip, graph, cost, vnoi, vbase, heap, state) );
5162 
5163  for( k = 0; k < nnodes; k++ )
5164  {
5165  /* is k a terminal other than the root? */
5166  if( Is_term(graph->term[k]) && graph->mark[k] && graph->grad[k] > 2 && k != graph->source )
5167  {
5168  assert(perm != NULL);
5169  for( l = 0; l < nterms; l++ )
5170  if( perm[l] == k )
5171  break;
5172 
5173  tmpcost = vnoi[k].dist + radii - radius[l];
5174 
5175  if( l == nterms - 1 )
5176  tmpcost -= radius[nterms - 2];
5177  else
5178  tmpcost -= radius[nterms - 1];
5179 
5180 
5181  if( SCIPisGT(scip, tmpcost, obj) )
5182  {
5183  SCIPdebugMessage("alternative bnd elimination!!! \n\n");
5184  (*nelims) += graph_pc_deleteTerm(scip, graph, k);
5185  (*offset) += graph->prize[k];
5186  }
5187  else
5188  {
5189  tmpcost += vnoi[k + nnodes].dist;
5190  if( l >= nterms - 2 )
5191  tmpcost -= radius[nterms - 3];
5192  else
5193  tmpcost -= radius[nterms - 2];
5194  if( SCIPisGT(scip, tmpcost, obj) || SCIPisGT(scip, mstobj2 + vnoi[k].dist + vnoi[k + nnodes].dist, obj) )
5195  {
5196  for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5197  if( graph->mark[graph->head[e]] && SCIPisLT(scip, graph->cost[e], graph->prize[k]) )
5198  break;
5199 
5200  if( e == EAT_LAST )
5201  {
5202  SCIPdebugMessage("second elimination!!! prize: %f \n\n", graph->prize[k]);
5203  (*offset) += graph->prize[k];
5204  (*nelims) += graph_pc_deleteTerm(scip, graph, k);
5205  }
5206  }
5207  }
5208  }
5209  }
5210  }
5211 
5212  SCIPdebugMessage("nelims (edges) in bound reduce: %d,\n", *nelims);
5213 
5214  /* free adjgraph */
5215  if( !mw )
5216  {
5217  graph_path_exit(scip, adjgraph);
5218  graph_free(scip, &adjgraph, TRUE);
5219  }
5220 
5221  /* free memory*/
5222  SCIPfreeBufferArrayNull(scip, &perm);
5223  SCIPfreeBufferArrayNull(scip, &mst);
5224  SCIPfreeBufferArrayNull(scip, &starts);
5225  SCIPfreeBufferArrayNull(scip, &stnode);
5226  SCIPfreeBufferArrayNull(scip, &result);
5227 
5228  assert(graph_valid(graph));
5229 
5230  return SCIP_OKAY;
5231 }
5232 
5233 
5234 
5235 
5236 /** Bound-based reduction method for the MWCSP .
5237  * The reduction method tries to eliminate nodes and vertices
5238  * by checking whether an upper bound for each solution that contains them
5239  * is smaller than the best known solution value.
5240  * The essence of the approach is a decomposition of the graph such that this upper bound
5241  * is minimized.
5242  * */
5244  SCIP* scip, /**< SCIP data structure */
5245  GRAPH* graph, /**< graph data structure */
5246  PATH* vnoi, /**< Voronoi data structure (size 3 * nnodes) */
5247  PATH* path, /**< shortest path data structure (size nnodes) */
5248  SCIP_Real* cost, /**< edge cost array */
5249  SCIP_Real* radius, /**< radius array */
5250  SCIP_Real* costrev, /**< reversed edge cost array */
5251  SCIP_Real* offset, /**< pointer to the offset */
5252  int* heap, /**< heap array */
5253  int* state, /**< array to store state of a node during Voronoi computation*/
5254  int* vbase, /**< Voronoi base to each node */
5255  int* result, /**< solution array or NULL */
5256  int* nelims /**< pointer to store number of eliminated edges */
5257  )
5258 {
5259  PATH* mst;
5260  SCIP_Real* prize;
5261  SCIP_Real obj;
5262  SCIP_Real bound;
5263  SCIP_Real tmpcost;
5264  SCIP_Real radiim2;
5265  int e;
5266  int k;
5267  int head;
5268  int nterms;
5269  int nnodes;
5270  int nedges;
5271 
5272  assert(scip != NULL);
5273  assert(graph != NULL);
5274  assert(vnoi != NULL);
5275  assert(path != NULL);
5276  assert(cost != NULL);
5277  assert(radius != NULL);
5278  assert(costrev != NULL);
5279  assert(heap != NULL);
5280  assert(state != NULL);
5281  assert(vbase != NULL);
5282  assert(nelims != NULL);
5283  assert(graph->source >= 0);
5284 
5285  mst = NULL;
5286  prize = graph->prize;
5287  nedges = graph->edges;
5288  nnodes = graph->knots;
5289  nterms = graph->terms - 1;
5290  *nelims = 0;
5291 
5292  assert(prize != NULL);
5293 
5294  /* not more than two nodes of positive weight? */
5295  if( nterms <= 2 )
5296  /* return */
5297  return SCIP_OKAY;
5298 
5299  /* initialize cost and costrev array */
5300  for( e = 0; e < nedges; e++ )
5301  {
5302  cost[e] = graph->cost[e];
5303  costrev[e] = graph->cost[flipedge(e)];
5304 
5305  assert(SCIPisGE(scip, cost[e], 0.0));
5306  }
5307 
5308  /* compute decomposition of graph and radius values */
5309  graph_voronoiWithRadiusMw(scip, graph, path, cost, radius, vbase, heap, state);
5310 
5311  /* sum all radius values, exclude two radius values of lowest value */
5312  for( k = 0; k < nnodes; k++ )
5313  {
5314  if( !Is_term(graph->term[k]) || !graph->mark[k] )
5315  continue;
5316 
5317  assert(vbase[k] == k);
5318  assert(SCIPisGE(scip, prize[k], 0.0));
5319 
5320  if( SCIPisGE(scip, radius[k], FARAWAY) )
5321  radius[k] = 0.0;
5322  else
5323  {
5324  if( SCIPisGE(scip, radius[k], prize[k] ) )
5325  radius[k] = 0.0;
5326  else
5327  radius[k] = prize[k] - radius[k];
5328  }
5329  }
5330 
5331  for( k = 0; k < nnodes; k++ )
5332  {
5333  if( !graph->mark[k] || graph->grad[k] == 0 || Is_term(graph->term[k]) )
5334  continue;
5335  }
5336 
5337  /* build Voronoi regions */
5338  graph_voronoiMw(scip, graph, costrev, vnoi, vbase, heap, state);
5339 
5340  /* get 2nd next positive node to all non-positive nodes */
5341  graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5342 
5343  for( k = 0; k < nnodes; k++ )
5344  {
5345  if( !Is_term(graph->term[k]) || !graph->mark[k] )
5346  continue;
5347 
5348  assert(vbase[k] == k);
5349  }
5350 
5351  SCIPsortReal(radius, nnodes);
5352 
5353  radiim2 = 0.0;
5354 
5355  for( k = 2; k < nterms; k++ )
5356  {
5357  assert( SCIPisGT(scip, FARAWAY, radius[k]) );
5358  radiim2 += radius[k];
5359  }
5360 
5361  /* solution available? */
5362  if( result != NULL)
5363  {
5364  /* calculate objective value of solution */
5365  obj = 0.0;
5366  for( e = 0; e < nedges; e++ )
5367  {
5368  if( result[e] == CONNECT )
5369  {
5370  head = graph->head[e];
5371 
5372  if( graph->mark[head] )
5373  obj += prize[head];
5374  }
5375  }
5376  }
5377  else
5378  {
5379  obj = 0.0;
5380  }
5381 
5382  bound = radiim2;
5383 
5384  /* traverse all nodes, try to eliminate each non-positive node */
5385  for( k = 0; k < nnodes; k++ )
5386  {
5387  if( !graph->mark[k] || graph->grad[k] == 0 || Is_term(graph->term[k]) )
5388  continue;
5389 
5390  assert(SCIPisLE(scip, graph->prize[k], 0.0));
5391 
5392  tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + bound + graph->prize[k];
5393 
5394  if( (SCIPisLT(scip, tmpcost, obj)) )
5395  {
5396  while( graph->outbeg[k] != EAT_LAST )
5397  {
5398  (*nelims)++;
5399  graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5400  }
5401  }
5402  }
5403 
5404  SCIPdebugMessage("nelims (edges) in MWCSP bound reduce: %d,\n", *nelims);
5405 
5406  /* free memory*/
5407  SCIPfreeBufferArrayNull(scip, &mst);
5408 
5409  return SCIP_OKAY;
5410 }
5411 
5412 
5413 
5414 /** bound-based reductions for the (R)PCSTP, the MWCSP and the STP; used by prune heuristic */
5416  SCIP* scip, /**< SCIP data structure */
5417  GRAPH* graph, /**< graph data structure */
5418  PATH* vnoi, /**< Voronoi data structure */
5419  SCIP_Real* cost, /**< edge cost array */
5420  SCIP_Real* radius, /**< radius array */
5421  SCIP_Real* costrev, /**< reversed edge cost array */
5422  SCIP_Real* offset, /**< pointer to the offset */
5423  int* heap, /**< heap array */
5424  int* state, /**< array to store state of a node during Voronoi computation */
5425  int* vbase, /**< Voronoi base to each node */
5426  const int* solnode, /**< array of nodes of current solution that is not to be destroyed */
5427  const int* soledge, /**< array of edges of current solution that is not to be destroyed */
5428  int* nelims, /**< pointer to store number of eliminated edges */
5429  int minelims /**< minimum number of edges to be eliminated */
5430  )
5431 {
5432  SCIP_Real* const prize = graph->prize;
5433  SCIP_Real obj;
5434  SCIP_Real max;
5435  SCIP_Real bound;
5436  SCIP_Real tmpcost;
5437  SCIP_Real mstobj;
5438  SCIP_Real maxcost;
5439  SCIP_Real radiim2;
5440  SCIP_Real radiim3;
5441  const int root = graph->source;
5442  const int nnodes = graph->knots;
5443  const int nedges = graph->edges;
5444  const SCIP_Bool mw = (graph->stp_type == STP_MWCSP);
5445  const SCIP_Bool pc = (graph->stp_type == STP_RPCSPG) || (graph->stp_type == STP_PCSPG);
5446  const SCIP_Bool pcmw = (pc || mw);
5447  const int nterms = graph->terms - ( (mw || graph->stp_type == STP_PCSPG) ? 1 : 0 );
5448  SCIP_Bool eliminate;
5449 
5450  assert(scip != NULL);
5451  assert(vnoi != NULL);
5452  assert(cost != NULL);
5453  assert(heap != NULL);
5454  assert(graph != NULL);
5455  assert(state != NULL);
5456  assert(vbase != NULL);
5457  assert(nelims != NULL);
5458  assert(radius != NULL);
5459  assert(costrev != NULL);
5460  assert(solnode != NULL);
5461  assert(soledge != NULL);
5462  assert(graph->source >= 0);
5463  assert(graph_valid(graph));
5464 
5465  mstobj = 0.0;
5466  *nelims = 0;
5467 
5468  if( nterms <= 2 )
5469  return SCIP_OKAY;
5470 
5471  assert( graph->stp_type != STP_RPCSPG || Is_term(graph->term[root]) );
5472 
5473  if( !pcmw )
5474  for( int k = 0; k < nnodes; k++ )
5475  graph->mark[k] = (graph->grad[k] > 0);
5476 
5477  /* initialize cost and costrev array */
5478  maxcost = 0.0;
5479  for( int e = 0; e < nedges; e++ )
5480  {
5481  cost[e] = graph->cost[e];
5482  costrev[e] = graph->cost[flipedge(e)];
5483 
5484  assert(SCIPisGE(scip, cost[e], 0.0));
5485 
5486  if( graph->stp_type == STP_DHCSTP && SCIPisLT(scip, graph->cost[e], BLOCKED) && SCIPisGT(scip, graph->cost[e], maxcost) )
5487  maxcost = graph->cost[e];
5488  }
5489 
5490  /* no MWCSP? */
5491  if( !mw )
5492  {
5493  GRAPH* adjgraph;
5494  PATH* mst;
5495 
5496  SCIP_CALL( graph_init(scip, &adjgraph, nterms, MIN(nedges, (nterms - 1) * nterms), 1) );
5497 
5498  /* build Voronoi regions, concomitantly building adjgraph and computing radii values*/
5499  SCIP_CALL( graph_voronoiWithRadius(scip, graph, adjgraph, vnoi, radius, cost, costrev, vbase, heap, state) );
5500 
5501  /* get 2nd next terminals to all non-terminal nodes */
5502  graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5503 
5504  /* get 3th next terminals to all non-terminal nodes */
5505  graph_get3next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5506 
5507  assert(adjgraph != NULL);
5508  graph_knot_chg(adjgraph, 0, 0);
5509  adjgraph->source = 0;
5510  assert(graph_valid(adjgraph));
5511 
5512  /* compute MST on adjgraph */
5513  SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
5514  SCIP_CALL( graph_path_init(scip, adjgraph) );
5515  graph_path_exec(scip, adjgraph, MST_MODE, 0, adjgraph->cost, mst);
5516 
5517  max = -1.0;
5518  for( int k = 1; k < nterms; k++ )
5519  {
5520  const int e = mst[k].edge;
5521  assert(adjgraph->path_state[k] == CONNECT);
5522  assert(e >= 0);
5523  tmpcost = adjgraph->cost[e];
5524  mstobj += tmpcost;
5525  if( SCIPisGT(scip, tmpcost, max) )
5526  max = tmpcost;
5527  }
5528  mstobj -= max;
5529 
5530  /* free memory*/
5531  SCIPfreeBufferArray(scip, &mst);
5532  graph_path_exit(scip, adjgraph);
5533  graph_free(scip, &adjgraph, TRUE);
5534 
5535  }
5536  else
5537  {
5538 
5539  /* build Voronoi regions */
5540  graph_voronoiMw(scip, graph, costrev, vnoi, vbase, heap, state);
5541 
5542  /* get 2nd next positive node to all non-positive nodes */
5543  graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5544  }
5545 
5546 
5547  /* for (rooted) prize collecting an maximum weight problems: correct radius values */
5548  if( graph->stp_type == STP_RPCSPG )
5549  {
5550  assert(graph->mark[graph->source]);
5551  for( int k = 0; k < nnodes; k++ )
5552  {
5553  if( !Is_term(graph->term[k]) || !graph->mark[k] )
5554  continue;
5555 
5556  if( SCIPisGT(scip, radius[k], prize[k]) && k != graph->source )
5557  radius[k] = prize[k];
5558  }
5559  }
5560  else if( graph->stp_type == STP_PCSPG )
5561  {
5562  for( int k = 0; k < nnodes; k++ )
5563  {
5564  if( !graph->mark[k] )
5565  continue;
5566 
5567  if( Is_term(graph->term[k]) && SCIPisGT(scip, radius[k], prize[k]) )
5568  radius[k] = prize[k];
5569  }
5570  }
5571 
5572  /* sort radius values */
5573  if( !mw )
5574  SCIPsortReal(radius, nnodes);
5575 
5576  radiim2 = 0.0;
5577 
5578  if( mw )
5579  {
5580  for( int e = 0; e < nedges; e++ )
5581  costrev[e] = FARAWAY;
5582  bound = 0.0;
5583  radiim3 = 0.0;
5584  }
5585  else
5586  {
5587  for( int e = 0; e < nedges; e++ )
5588  costrev[e] = -1.0;
5589 
5590  /* sum all but two radius values of highest/lowest value */
5591  for( int k = 0; k < nterms - 2; k++ )
5592  {
5593  assert( SCIPisGT(scip, FARAWAY, radius[k]) );
5594  radiim2 += radius[k];
5595  }
5596  if( nterms >= 3 )
5597  radiim3 = radiim2 - radius[nterms - 3];
5598  else
5599  radiim3 = 0.0;
5600 
5601  if( SCIPisGT(scip, radiim2, mstobj) )
5602  bound = radiim2;
5603  else
5604  bound = mstobj;
5605  }
5606 
5607  for( int redrounds = 0; redrounds < 3; redrounds++ )
5608  {
5609  int nrealelims = MIN(2 * minelims, nedges - 1);
5610 
5611  if( redrounds == 0 )
5612  {
5613  eliminate = FALSE;
5614  obj = FARAWAY;
5615  }
5616  else if( redrounds == 1 )
5617  {
5618  assert(minelims > 0);
5619  assert(2 * minelims < nedges);
5620  eliminate = TRUE;
5621  SCIPsortReal(costrev, nedges);
5622 
5623  if( mw )
5624  {
5625  obj = costrev[nrealelims];
5626  }
5627  else
5628  {
5629  obj = costrev[nedges - nrealelims];
5630 
5631  if( SCIPisLT(scip, obj, 0.0) )
5632  obj = 0.0;
5633  }
5634  }
5635  else
5636  {
5637  if( mw )
5638  {
5639  obj = costrev[nrealelims] + 2 * SCIPepsilon(scip);
5640  }
5641  else
5642  {
5643  obj = costrev[nedges - nrealelims] - 2 * SCIPepsilon(scip);
5644 
5645  if( SCIPisLT(scip, obj, 0.0) )
5646  obj = 0.0;
5647  }
5648  eliminate = TRUE;
5649  }
5650 
5651  if( mw )
5652  {
5653  for( int k = 0; k < nnodes; k++ )
5654  {
5655  if( (*nelims) >= minelims )
5656  break;
5657 
5658  if( root == k )
5659  continue;
5660 
5661  if( !graph->mark[k] || graph->grad[k] == 0 || Is_gterm(graph->term[k] ) )
5662  continue;
5663 
5664  tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + graph->prize[k];
5665 
5666  if( (!eliminate || SCIPisLT(scip, tmpcost, obj)) && solnode[k] != CONNECT )
5667  {
5668  if( eliminate )
5669  {
5670  while (graph->outbeg[k] != EAT_LAST)
5671  {
5672  (*nelims)++;
5673  graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5674  }
5675  }
5676  else
5677  {
5678  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5679  {
5680  if( SCIPisLT(scip, tmpcost, costrev[e]) )
5681  {
5682  costrev[e] = tmpcost;
5683  costrev[flipedge(e)] = tmpcost;
5684  }
5685  }
5686  }
5687  }
5688  else if( solnode[k] == CONNECT )
5689  {
5690  int e = graph->outbeg[k];
5691 
5692  while( e != EAT_LAST )
5693  {
5694  const int etemp = graph->oeat[e];
5695  const int tail = graph->tail[e];
5696  const int head = graph->head[e];
5697 
5698  if( tail > head || soledge[e] == CONNECT || soledge[flipedge(e)] == CONNECT )
5699  {
5700  e = etemp;
5701  continue;
5702  }
5703 
5704  tmpcost = graph->prize[head] + graph->prize[tail];
5705 
5706  if( vbase[tail] != vbase[head] )
5707  {
5708  tmpcost -= vnoi[head].dist + vnoi[tail].dist;
5709  }
5710  else
5711  {
5712  if( SCIPisGT(scip, -vnoi[tail].dist -vnoi[head + nnodes].dist, -vnoi[tail + nnodes].dist -vnoi[head].dist) )
5713  tmpcost -= vnoi[tail].dist + vnoi[head + nnodes].dist;
5714  else
5715  tmpcost -= vnoi[tail + nnodes].dist + vnoi[head].dist;
5716  }
5717  /* can edge e or arc e be deleted? */
5718  if( (!eliminate || SCIPisLT(scip, tmpcost, obj))
5719  && SCIPisLT(scip, graph->cost[e], FARAWAY) && (graph->mark[head]) )
5720  {
5721  assert(!Is_pterm(graph->term[head]));
5722  assert(!Is_pterm(graph->term[tail]));
5723 
5724  if( eliminate )
5725  {
5726  graph_edge_del(scip, graph, e, TRUE);
5727  (*nelims)++;
5728  }
5729  else if( SCIPisLT(scip, tmpcost, costrev[e]) )
5730  {
5731  costrev[e] = tmpcost;
5732  costrev[flipedge(e)] = tmpcost;
5733  }
5734 
5735  }
5736  e = etemp;
5737  }
5738  }
5739  }
5740  }
5741  /* no MWCSP */
5742  else
5743  {
5744  /* traverse all nodes, try to eliminate each node or incident edges */
5745  for( int k = 0; k < nnodes; k++ )
5746  {
5747  if( (*nelims) >= minelims )
5748  break;
5749  if( root == k )
5750  continue;
5751 
5752  if( (!graph->mark[k] && (pcmw)) || graph->grad[k] == 0 )
5753  continue;
5754 
5755  tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + bound;
5756 
5757  /* can node k be deleted? */
5758  if( !Is_term(graph->term[k]) && (!eliminate || SCIPisGT(scip, tmpcost, obj)) && solnode[k] != CONNECT )
5759  {
5760  /* delete all incident edges */
5761  if( eliminate )
5762  {
5763  while( graph->outbeg[k] != EAT_LAST )
5764  {
5765  const int e = graph->outbeg[k];
5766  (*nelims)++;
5767 
5768  assert(!pc || graph->tail[e] != graph->source);
5769  assert(!pc || graph->mark[graph->head[e]]);
5770  assert(!Is_pterm(graph->term[graph->head[e]]));
5771  assert(!Is_pterm(graph->term[graph->tail[e]]));
5772 
5773  graph_edge_del(scip, graph, e, TRUE);
5774  }
5775  }
5776  else
5777  {
5778  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5779  {
5780  if( SCIPisGT(scip, tmpcost, costrev[e]) )
5781  {
5782  costrev[e] = tmpcost;
5783  costrev[flipedge(e)] = tmpcost;
5784  }
5785  }
5786  }
5787  }
5788  else
5789  {
5790  int e = graph->outbeg[k];
5791  while( e != EAT_LAST )
5792  {
5793  const int etemp = graph->oeat[e];
5794  const int tail = graph->tail[e];
5795  const int head = graph->head[e];
5796 
5797  if( tail > head )
5798  {
5799  e = etemp;
5800  continue;
5801  }
5802 
5803  if( soledge[e] == CONNECT || soledge[flipedge(e)] == CONNECT )
5804  {
5805  e = etemp;
5806  continue;
5807  }
5808 
5809  tmpcost = graph->cost[e] + bound;
5810 
5811  if( vbase[tail] != vbase[head] )
5812  {
5813  tmpcost += vnoi[head].dist + vnoi[tail].dist;
5814  }
5815  else
5816  {
5817  if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5818  tmpcost += vnoi[tail + nnodes].dist + vnoi[head].dist;
5819  else
5820  tmpcost += vnoi[tail].dist + vnoi[head + nnodes].dist;
5821 
5822  assert(SCIPisGE(scip, tmpcost, vnoi[head].dist + vnoi[tail].dist + graph->cost[e] + bound));
5823  }
5824  /* can edge e or arc e be deleted? */
5825  if( (!eliminate || SCIPisGT(scip, tmpcost, obj))
5826  && SCIPisLT(scip, graph->cost[e], FARAWAY) && (!(pc) || graph->mark[head]) )
5827  {
5828  assert(!Is_pterm(graph->term[head]));
5829  assert(!Is_pterm(graph->term[tail]));
5830 
5831  if( eliminate )
5832  {
5833  graph_edge_del(scip, graph, e, TRUE);
5834  (*nelims)++;
5835  }
5836  else if( SCIPisGT(scip, tmpcost, costrev[e]) )
5837  {
5838  costrev[e] = tmpcost;
5839  costrev[flipedge(e)] = tmpcost;
5840  }
5841  }
5842  e = etemp;
5843  }
5844  }
5845  }
5846 #if 1
5847  /* traverse all nodes, try to eliminate 3,4 degree nodes */
5848  for( int k = 0; k < nnodes; k++ )
5849  {
5850  if( (*nelims) >= minelims )
5851  break;
5852 
5853  if( (!graph->mark[k] && pc) || graph->grad[k] <= 0 )
5854  continue;
5855 
5856  if( solnode[k] == CONNECT )
5857  continue;
5858 
5859  if( !eliminate )
5860  {
5861  if( graph->grad[k] >= 3 && graph->grad[k] <= 4 && !Is_term(graph->term[k]) )
5862  {
5863  tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5864  for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5865  {
5866  if( SCIPisGT(scip, tmpcost, costrev[e]) )
5867  {
5868  costrev[e] = tmpcost;
5869  costrev[flipedge(e)] = tmpcost;
5870  }
5871  }
5872  }
5873  continue;
5874  }
5875 
5876  if( graph->grad[k] >= 3 && graph->grad[k] <= 4 && !Is_term(graph->term[k]) )
5877  {
5878  tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5879  if( SCIPisGT(scip, tmpcost, obj) )
5880  {
5881  SCIP_Bool success;
5882 
5883  SCIP_CALL( graph_knot_delPseudo(scip, graph, graph->cost, NULL, NULL, k, &success) );
5884 
5885  assert(graph->grad[k] == 0 || (graph->grad[k] == 4 && !success));
5886 
5887  if( success )
5888  (*nelims)++;
5889  }
5890  }
5891  }
5892  } /* no MWCS */
5893 #endif
5894  } /* redrounds */
5895