Scippy

SCIP

Solving Constraint Integer Programs

reduce_termsepa.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-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file reduce_termsepa.c
17  * @brief base class for terminal-separator reductions for Steiner tree problems
18  * @author Daniel Rehfeldt
19  *
20  * This file implements base classes for terminal-separator reduction techniques Steiner tree problems.
21  * The actual algorithms can be found in reduce_termsepada.c and reduce_termsepafull.c
22  *
23  * A list of all interface methods can be found in reduce.h.
24  *
25  */
26 
27 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
28 //#define SCIP_DEBUG
29 
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <string.h>
33 #include <assert.h>
34 #include "graph.h"
35 #include "reduce.h"
36 #include "extreduce.h"
37 #include "solstp.h"
38 #include "mincut.h"
39 #include "portab.h"
40 #include "stpvector.h"
41 #include "scip/scip.h"
42 
43 
44 
45 #define MARK_NONACTIVE 0
46 #define MARK_SUBNODE 1
47 #define MARK_SEPARATOR 2
48 
49 /*
50  * Data structures
51  */
52 
53 
54 /** tree bottleneck node */
55 typedef struct tree_bottleneck_node
56 {
57  SCIP_Real edgecost; /**< cost of outgoing edge */
58  SCIP_Real bottleneck; /**< (temporary) bottleneck up to this node; -1.0 for UNSET */
59  int parent; /**< parent node index */
60  int degree; /**< degree of node */
61 } TBNODE;
62 
63 
64 
65 /** terminal separator tree bottleneck structure */
67 {
68  TBNODE* tree; /**< tree for representation; NOTE: oriented towards root */
69  int nnodes; /**< number of nodes of underlying graph */
70 };
71 
72 
73 
74 /*
75  * Local methods
76  */
77 
78 
79 
80 /** initializes */
81 static
83  SCIP* scip, /**< SCIP data structure */
84  const GRAPH* g, /**< graph data structure */
85  const int* soledges, /**< solution tree to be represented */
86  TBOTTLENECK** tbottleneck /**< to initialize */
87  )
88 {
89  TBOTTLENECK* tbneck;
90  TBNODE* tbtree;
91  const int nnodes = graph_get_nNodes(g);
92  const int nedges = graph_get_nEdges(g);
93 
94  assert(scip && soledges);
95  assert(solstp_isValid(scip, g, soledges));
96 
97  SCIP_CALL( SCIPallocMemory(scip, tbottleneck) );
98  tbneck = *tbottleneck;
99 
100  SCIP_CALL( SCIPallocMemoryArray(scip, &(tbneck->tree), nnodes) );
101  tbtree = tbneck->tree;
102  tbneck->nnodes = nnodes;
103 
104  for( int i = 0; i < nnodes; i++ )
105  {
106  tbtree[i].edgecost = -FARAWAY;
107  tbtree[i].parent = UNKNOWN;
108  tbtree[i].bottleneck = -1.0;
109 
110  /* NOTE: slight hack to make sure that terminals are not taken as part of key path */
111  if( Is_term(g->term[i]) )
112  {
113  tbtree[i].degree = nnodes;
114  }
115  else
116  {
117  tbtree[i].degree = 0.0;
118  }
119  }
120 
121 #ifdef SCIP_DEBUG
122  graph_printInfo(g);
123 #endif
124 
125  for( int e = 0; e < nedges; e++ )
126  {
127  if( soledges[e] == CONNECT )
128  {
129  const int tail = g->tail[e];
130  const int head = g->head[e];
131 
132 #ifdef SCIP_DEBUG
133  SCIPdebugMessage("soledge: ");
134  graph_edge_printInfo(g, e);
135 #endif
136 
137  assert(tbtree[head].parent == UNKNOWN);
138 
139  tbtree[head].parent = tail;
140  tbtree[head].edgecost = g->cost[e];
141 
142  tbtree[tail].degree++;
143  tbtree[head].degree++;
144  }
145  }
146 
147  assert(tbtree[g->source].parent == UNKNOWN);
148 
149  return SCIP_OKAY;
150 }
151 
152 // todo really needed? There should be not need for the key-paths to be disjoint!
153 #ifdef SCIP_DISABLED
154 /** helper */
155 static
156 void tbottleneckCut(
157  int startnode, /**< start node */
158  int endnode, /**< end node */
159  SCIP_Real maxlength, /**< length of bottleneck */
160  TBOTTLENECK* tbottleneck /**< tree bottleneck structure */
161  )
162 {
163  TBNODE* tbtree = tbottleneck->tree;
164  int k;
165  int cutstart = startnode;
166  int cutend;
167  SCIP_Real max_local = 0.0;
168 #ifndef NDEBUG
169  SCIP_Real max_debug = 0.0;
170 #endif
171 
172  assert(startnode != endnode);
173  assert(tbottleneck->nnodes >= 2);
174 
175  /* go up to the root, stop at end node */
176  for( k = startnode; k != endnode; k = tbtree[k].parent )
177  {
178  assert(0 <= k && k < tbottleneck->nnodes);
179 
180  if( GE(max_local, maxlength) )
181  {
182  break;
183  }
184 
185  assert(tbtree[k].parent != UNKNOWN);
186 
187  if( tbtree[k].degree > 2 )
188  {
189  cutstart = k;
190  max_local = 0.0;
191  }
192 
193  max_local += tbtree[k].edgecost;
194  }
195 
196  assert(EQ(max_local, maxlength));
197  cutend = k;
198 
199  for( k = cutstart; k != cutend; k = tbtree[k].parent )
200  {
201  assert(k == cutstart || tbtree[k].degree <= 2);
202 #ifndef NDEBUG
203  max_debug += tbtree[k].edgecost;
204 #endif
205 
206  SCIPdebugMessage("removing tree bottleneck node=%d, degree=%d edgecost=%f \n", k, tbtree[k].degree, tbtree[k].edgecost);
207  tbtree[k].edgecost = -FARAWAY;
208  }
209 
210 #ifndef NDEBUG
211  assert(EQ(max_debug, maxlength));
212 #endif
213 
214 }
215 #endif
216 
217 /** gets tree bottleneck cost */
218 static
220  int node1, /**< first node */
221  int node2, /**< second node */
222  TBOTTLENECK* tbottleneck, /**< tree bottleneck structure */
223  SCIP_Real* maxlength /**< IN/OUT! user needs to provide a maximum that is to be surpassed */
224  )
225 {
226  int k;
227  TBNODE* tbtree = tbottleneck->tree;
228  SCIP_Real max_local = 0.0;
229  SCIP_Real max1 = 0.0;
230  SCIP_Real max2 = 0.0;
231  SCIP_Real maxtotal;
232 
233  assert(node1 >= 0 && node2 >= 0);
234  assert(node1 != node2);
235  assert(GT(*maxlength, 0.0));
236 
237  /* go up to the root and mark the tree on the way */
238  for( k = node1; k != UNKNOWN; k = tbtree[k].parent )
239  {
240  assert(0 <= k && k < tbottleneck->nnodes);
241  assert(EQ(tbtree[k].bottleneck, -1.0));
242 
243  tbtree[k].bottleneck = max1;
244 
245  if( tbtree[k].degree > 2 )
246  max_local = 0.0;
247 
248  if( tbtree[k].parent != UNKNOWN )
249  {
250  max_local += tbtree[k].edgecost;
251  }
252 
253  if( max_local > max1 )
254  max1 = max_local;
255  }
256 
257  max_local = 0.0;
258 
259  /* go up from second node */
260  for( k = node2; tbtree[k].bottleneck < -0.5; k = tbtree[k].parent )
261  {
262  assert(0 <= k && k < tbottleneck->nnodes);
263  /* NOTE: we should not reach the root */
264  assert(tbtree[k].parent != UNKNOWN);
265 
266  if( tbtree[k].degree > 2 )
267  max_local = 0.0;
268 
269  max_local += tbtree[k].edgecost;
270 
271  if( max_local > max2 )
272  max2 = max_local;
273  }
274 
275  assert(GE(tbtree[k].bottleneck, 0.0));
276  max1 = tbtree[k].bottleneck;
277 
278  SCIPdebugMessage("joint node: %d \n", k);
279  SCIPdebugMessage("tree bdist from node %d: %f \n", node1, max1);
280  SCIPdebugMessage("tree bdist from node %d: %f \n", node2, max2);
281 
282  maxtotal = MAX(max1, max2);
283 
284  /* not yet removed sub-path and improvement over given distance? */
285  if( GT(maxtotal, *maxlength) )
286  {
287  SCIPdebugMessage("improved tree bottleneck %f -> %f \n", *maxlength, maxtotal);
288 
289  *maxlength = maxtotal;
290 
291 #ifdef SCIP_DISABELD
292  if( GT(max1, max2) )
293  tbottleneckCut(node1, k, max1, tbottleneck);
294  else
295  tbottleneckCut(node2, k, max2, tbottleneck);
296 #endif
297  }
298 
299  /* reset */
300  for( k = node1; k != UNKNOWN; k = tbtree[k].parent )
301  {
302  assert(0 <= k && k < tbottleneck->nnodes);
303  tbtree[k].bottleneck = -1.0;
304  }
305 }
306 
307 
308 /** frees */
309 static
311  SCIP* scip, /**< SCIP data structure */
312  TBOTTLENECK** tbottleneck /**< to initialize */
313  )
314 {
315  TBOTTLENECK* tbneck;
316  tbneck = *tbottleneck;
317 
318  SCIPfreeMemoryArray(scip, &(tbneck->tree));
319  SCIPfreeMemory(scip, tbottleneck);
320 }
321 
322 
323 /** identifies subgraph */
324 static
326  SCIP* scip, /**< SCIP data structure */
327  GRAPH* g, /**< graph data structure */
328  TERMCOMP* termcomp /**< component */
329  )
330 {
331  const COMPBUILDER* const builder = termcomp->builder;
332  const int* const sepaterms = builder->sepaterms;
333  int* const nodemap_orgToSub = termcomp->nodemap_orgToSub;
334  STP_Vectype(int) bfsqueue = NULL;
335  int* RESTRICT gmark = termcomp->nodes_mark;
336  int sub_e = 0;
337  int sub_n = 0;
338  const int nsepaterms = builder->nsepatterms;
339  const int sourceterm = builder->sourceterm;
340 
341  assert(graph_knot_isInRange(g, sourceterm) && Is_term(g->term[sourceterm]));
342  assert(sepaterms && nodemap_orgToSub);
343  assert(nsepaterms >= 2 && nsepaterms < g->terms);
344 
345  BMSclearMemoryArray(gmark, g->knots);
346 
347  /* mark separator */
348  for( int i = 0; i < nsepaterms; i++ )
349  {
350  const int term = sepaterms[i];
351  assert(graph_knot_isInRange(g, term) && Is_term(g->term[term]));
352  assert(!gmark[term]);
353  assert(nodemap_orgToSub[term] == UNKNOWN);
354 
355  SCIPdebugMessage("mapping separator %d to %d \n", term, sub_n);
356 
357  nodemap_orgToSub[term] = sub_n;
358  sub_n++;
359  sub_e += g->grad[term];
360  gmark[term] = MARK_SEPARATOR;
361  }
362 
363  assert(!gmark[sourceterm]);
364  StpVecReserve(scip, bfsqueue, 16);
365  StpVecPushBack(scip, bfsqueue, sourceterm);
366  gmark[sourceterm] = MARK_SUBNODE;
367 
368  /* BFS loop */
369  for( int i = 0; i < StpVecGetSize(bfsqueue); i++ )
370  {
371  const int k = bfsqueue[i];
372 
373  assert(gmark[k]);
374  assert(nodemap_orgToSub[k] == UNKNOWN);
375 
376  nodemap_orgToSub[k] = sub_n;
377 
378  SCIPdebugMessage("mapping component node %d to %d \n", k, sub_n);
379 
380  sub_n++;
381  sub_e += g->grad[k];
382 
383  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
384  {
385  const int head = g->head[e];
386  if( !gmark[head] )
387  {
388  StpVecPushBack(scip, bfsqueue, head);
389  gmark[head] = MARK_SUBNODE;
390  }
391  }
392  }
393 
394  assert(termcomp->bfsqueue == NULL);
395  assert(termcomp->subnnodes == -1);
396  assert(termcomp->subnedges == -1);
397 
398  /* reserve space for the separator terminal clique */
399  sub_e += (nsepaterms) * (nsepaterms - 1);
400 
401  termcomp->bfsqueue = bfsqueue;
402  termcomp->subnnodes = sub_n;
403  termcomp->subnedges = sub_e;
404 }
405 
406 
407 /** creates subgraph */
408 static
410  SCIP* scip, /**< SCIP data structure */
411  const GRAPH* orggraph, /**< graph data structure */
412  SCIP_Bool useSd, /**< use special distances */
413  EXTPERMA* extperma, /**< extension data or NULL */
414  TERMCOMP* termcomp /**< component */
415  )
416 {
417  GRAPH* subgraph;
418  DISTDATA* const distdata = useSd ? extperma->distdata_default : NULL;
419  const COMPBUILDER* const builder = termcomp->builder;
420  const int* const nodemap_orgToSub = termcomp->nodemap_orgToSub;
421  int* nodemap_subToOrg;
422  int* edgemap_subToOrg;
423  const int* const orgmark = termcomp->nodes_mark;
424  const int* const sepaterms = builder->sepaterms;
425  const int nsepaterms = builder->nsepatterms;
426  const int nnodes_sub = termcomp->subnnodes;
427  const int nedges_sub = termcomp->subnedges;
428  STP_Vectype(int) bfsqueue = termcomp->bfsqueue;
429 
430  assert(nsepaterms >= 2 && nsepaterms < orggraph->terms);
431  assert(nnodes_sub > 0 && nedges_sub > 0);
432  assert(nodemap_orgToSub && sepaterms);
433  assert(!termcomp->subgraph);
434  assert((extperma != NULL) == useSd);
435 
436  SCIP_CALL( graph_init(scip, &(termcomp->subgraph), nnodes_sub, nedges_sub, 1) );
437  subgraph = termcomp->subgraph;
438 
439  SCIP_CALL( SCIPallocMemoryArray(scip, &(nodemap_subToOrg), nnodes_sub) );
440  SCIP_CALL( SCIPallocMemoryArray(scip, &(edgemap_subToOrg), nedges_sub) );
441 
442  for( int i = 0; i < nsepaterms; i++ )
443  {
444  assert(graph_knot_isInRange(orggraph, sepaterms[i]) && Is_term(orggraph->term[sepaterms[i]]));
445  assert(nodemap_orgToSub[sepaterms[i]] == subgraph->knots);
446 
447  graph_knot_add(subgraph, STP_TERM);
448  nodemap_subToOrg[i] = sepaterms[i];
449  }
450 
451  for( int i = 0; i < StpVecGetSize(bfsqueue); i++ )
452  {
453  const int orgnode = bfsqueue[i];
454  assert(nodemap_orgToSub[orgnode] == subgraph->knots);
455 
456  nodemap_subToOrg[subgraph->knots] = orgnode;
457  graph_knot_add(subgraph, orggraph->term[orgnode]);
458  }
459 
460 #ifndef NDEBUG
461  assert(nnodes_sub == subgraph->knots);
462  for( int i = 0; i < nnodes_sub; i++ )
463  {
464  assert(i == nodemap_orgToSub[nodemap_subToOrg[i]]);
465  }
466 #endif
467 
468  for( int subsepaterm = 0; subsepaterm < nsepaterms; subsepaterm++ )
469  {
470  const int orgsepaterm = sepaterms[subsepaterm];
471 
472  assert(graph_knot_isInRange(orggraph, orgsepaterm) && Is_term(orggraph->term[orgsepaterm]));
473  assert(orgmark[orgsepaterm]);
474  assert(nodemap_subToOrg[subsepaterm] == orgsepaterm);
475 
476  for( int e = orggraph->outbeg[orgsepaterm]; e != EAT_LAST; e = orggraph->oeat[e] )
477  {
478  const int orghead = orggraph->head[e];
479 
480  /* NOTE: add only small to big edges, to avoid adding an edge twice */
481  if( orghead < orgsepaterm )
482  continue;
483 
484  /* NOTE: ignore separator nodes, to not add edge twice */
485  if( orgmark[orghead] == MARK_SUBNODE )
486  {
487  const int subnode = nodemap_orgToSub[orghead];
488  assert(subnode >= 0);
489  assert(EQ(orggraph->cost[e], orggraph->cost[flipedge(e)]));
490 
491  edgemap_subToOrg[subgraph->edges] = e;
492  edgemap_subToOrg[subgraph->edges + 1] = flipedge(e);
493  graph_edge_addBi(scip, subgraph, subsepaterm, subnode, orggraph->cost[e]);
494 
495 #ifdef SCIP_DEBUG
496  SCIPdebugMessage("adding separator-to-default edge: ");
497  graph_edge_printInfo(subgraph, subgraph->edges - 2);
498 #endif
499  }
500  }
501 
502  for( int subsepaterm2 = subsepaterm + 1; subsepaterm2 < nsepaterms; subsepaterm2++ )
503  {
504  const int orgsepaterm2 = sepaterms[subsepaterm2];
505  const SCIP_Real sd = useSd ?
506  extreduce_distDataGetSdDouble(scip, orggraph, orgsepaterm, orgsepaterm2, distdata)
507  : BLOCKED;
508 
509  assert(GE(sd, 0.0));
510 
511  edgemap_subToOrg[subgraph->edges] = -1;
512  edgemap_subToOrg[subgraph->edges + 1] = -1;
513  graph_edge_addBi(scip, subgraph, subsepaterm, subsepaterm2, sd);
514 
515 #ifdef SCIP_DEBUG
516  SCIPdebugMessage("adding separator-to-separator edge: ");
517  graph_edge_printInfo(subgraph, subgraph->edges - 2);
518 #endif
519  }
520  }
521 
522  for( int i = 0; i < StpVecGetSize(bfsqueue); i++ )
523  {
524  const int orgnode = bfsqueue[i];
525  const int subnode = nodemap_orgToSub[orgnode];
526 
527  for( int e = orggraph->outbeg[orgnode]; e != EAT_LAST; e = orggraph->oeat[e] )
528  {
529  const int orghead = orggraph->head[e];
530 
531  if( orghead < orgnode )
532  continue;
533 
534  if( orgmark[orghead] )
535  {
536  const int subhead = nodemap_orgToSub[orghead];
537  assert(subhead >= 0);
538  assert(EQ(orggraph->cost[e], orggraph->cost[flipedge(e)]));
539 
540  edgemap_subToOrg[subgraph->edges] = e;
541  edgemap_subToOrg[subgraph->edges + 1] = flipedge(e);
542  graph_edge_addBi(scip, subgraph, subnode, subhead, orggraph->cost[e]);
543 
544 #ifdef SCIP_DEBUG
545  SCIPdebugMessage("adding default-to-? edge: ");
546  graph_edge_printInfo(subgraph, subgraph->edges - 2);
547 #endif
548  }
549  }
550  }
551 
552  /* finally, we need to update the edge map on sepa-sepa edges */
553  for( int subsepaterm = 0; subsepaterm < nsepaterms; subsepaterm++ )
554  {
555  const int orgsepaterm = sepaterms[subsepaterm];
556 
557  assert(orgmark[orgsepaterm]);
558  assert(nodemap_subToOrg[subsepaterm] == orgsepaterm);
559 
560  for( int orge = orggraph->outbeg[orgsepaterm]; orge != EAT_LAST; orge = orggraph->oeat[orge] )
561  {
562  const int orghead = orggraph->head[orge];
563  if( orghead < orgsepaterm )
564  continue;
565 
566  if( orgmark[orghead] == MARK_SEPARATOR )
567  {
568  const int subhead = nodemap_orgToSub[orghead];
569  assert(subhead >= 0);
570 
571  /* find and update the corresponding edge in subgraph */
572  for( int e = subgraph->outbeg[subsepaterm]; e != EAT_LAST; e = subgraph->oeat[e] )
573  {
574  if( subgraph->head[e] == subhead )
575  {
576  const SCIP_Real newcost = MIN(orggraph->cost[orge], subgraph->cost[e]);
577  assert(edgemap_subToOrg[e] == -1);
578 
579  edgemap_subToOrg[e] = orge;
580  edgemap_subToOrg[flipedge(e)] = flipedge(orge);
581 
582  subgraph->cost[e] = newcost;
583  subgraph->cost[flipedge(e)] = newcost;
584 
585  break;
586  }
587  }
588  }
589  }
590  }
591 
592  assert(graph_knot_isInRange(orggraph, builder->sourceterm));
593  assert(graph_knot_isInRange(orggraph, nodemap_orgToSub[builder->sourceterm]));
594  assert(orggraph->stp_type != UNKNOWN);
595 
596  subgraph->source = nodemap_orgToSub[builder->sourceterm];
597  subgraph->stp_type = orggraph->stp_type;
598 
599  SCIP_CALL( graph_path_init(scip, subgraph) );
600 
601  assert(!termcomp->nodemap_subToOrg && !termcomp->edgemap_subToOrg);
602 
603  termcomp->nodemap_subToOrg = nodemap_subToOrg;
604  termcomp->edgemap_subToOrg = edgemap_subToOrg;
605 
606  return SCIP_OKAY;
607 }
608 
609 
610 
611 /** gets extended bottleneck distances */
612 static
614  const GRAPH* g, /**< graph data structure */
615  int term1, /**< first terminal */
616  int term2, /**< second terminal */
617  int term1_sub, /**< corresponding terminal in subgraph */
618  int term2_sub, /**< corresponding terminal in subgraph */
619  int mstsource, /**< source node of minimum spanning tree */
620  PATH* mst, /**< minimum spanning tree */
621  TBOTTLENECK* subsolbottleneck, /**< tree bottleneck */
622  SCIP_Real* RESTRICT mstsdist /**< node distance helper */
623  )
624 {
625  SCIP_Real bdist = 0.0;
626  int tempnode = term1;
627 
628  assert(mst && mstsdist);
629  assert(Is_term(g->term[term1]));
630  assert(Is_term(g->term[term2]));
631  assert(term1 != term2);
632 
633  mstsdist[tempnode] = 0.0;
634 
635  while( tempnode != mstsource )
636  {
637  const int ne = mst[tempnode].edge;
638 
639  assert(ne >= 0);
640  assert(g->head[ne] == tempnode);
641  tempnode = g->tail[ne];
642 
643  if( g->cost[ne] > bdist )
644  bdist = g->cost[ne];
645 
646  mstsdist[tempnode] = bdist;
647  if( tempnode == term2 )
648  break;
649  }
650 
651  /* already finished? */
652  if( tempnode == term2 )
653  {
654  tempnode = mstsource;
655  }
656  else
657  {
658  tempnode = term2;
659  bdist = 0.0;
660  }
661 
662  while( tempnode != mstsource )
663  {
664  const int ne = mst[tempnode].edge;
665  assert(ne >= 0);
666 
667  tempnode = g->tail[ne];
668 
669  if( g->cost[ne] > bdist )
670  bdist = g->cost[ne];
671 
672  /* already visited? */
673  if( mstsdist[tempnode] > -0.5 )
674  {
675  if( mstsdist[tempnode] > bdist )
676  bdist = mstsdist[tempnode];
677  break;
678  }
679 
680  assert(EQ(mstsdist[tempnode], -1.0));
681  }
682 
683  /* restore */
684  tempnode = term1;
685  mstsdist[tempnode] = -1.0;
686  while( tempnode != mstsource )
687  {
688  const int ne = mst[tempnode].edge;
689  tempnode = g->tail[ne];
690  mstsdist[tempnode] = -1.0;
691  if( tempnode == term2 )
692  break;
693  }
694 
695 #ifndef NDEBUG
696  for( int i = 0; i < g->knots; i++ )
697  assert(EQ(mstsdist[i], -1.0));
698 #endif
699 
700  /* updates with tree bottleneck of sub-solution if available */
701  if( subsolbottleneck )
702  tbottleneckGetMax(term1_sub, term2_sub, subsolbottleneck, &bdist);
703 
704  return bdist;
705 }
706 
707 
708 
709 /*
710  * Interface methods
711  */
712 
713 
714 
715 /** initializes */
717  SCIP* scip, /**< SCIP data structure */
718  const GRAPH* g, /**< graph data structure */
719  COMPBUILDER** compbuilder /**< to initialize */
720  )
721 {
722  COMPBUILDER* builder;
723  const int nnodes = graph_get_nNodes(g);
724 
725  SCIP_CALL( SCIPallocMemory(scip, compbuilder) );
726  builder = *compbuilder;
727 
728  builder->sepaterms = NULL;
729  builder->sourceterm = -1;
730  builder->nsepatterms = 2;
731  builder->componentnumber = 0;
732  builder->ncomponentnodes = -1;
733  builder->ngraphnodes = -1;
734  builder->rootcompIsProcessed = FALSE;
735  builder->maxncompchecks = -1;
736  builder->maxsepasize = -1;
737 
738  SCIP_CALL( SCIPallocMemoryArray(scip, &(builder->nodes_bdist), nnodes) );
739  graph_get_nVET(g, &(builder->ngraphnodes), NULL, NULL);
740 
741  for( int i = 0; i < nnodes; i++ )
742  {
743  builder->nodes_bdist[i] = -1.0;
744  }
745 
746  return SCIP_OKAY;
747 }
748 
749 
750 /** frees */
752  SCIP* scip, /**< SCIP data structure */
753  COMPBUILDER** compbuilder /**< to free */
754  )
755 {
756  COMPBUILDER* builder;
757  builder = *compbuilder;
758 
759  SCIPfreeMemoryArray(scip, &(builder->nodes_bdist));
760  SCIPfreeMemory(scip, compbuilder);
761 }
762 
763 
764 /** gets nodes ratio of subgraph */
766  const COMPBUILDER* compbuilder /**< builder */
767  )
768 {
769  const SCIP_Real ratio = (SCIP_Real) compbuilder->ncomponentnodes / (SCIP_Real) compbuilder->ngraphnodes;
770 
771  assert(GT(ratio, 0.0));
772  assert(LE(ratio, 1.0));
773 
774  return ratio;
775 }
776 
777 
778 /** prints */
780  const GRAPH* orggraph, /**< graph data structure */
781  const COMPBUILDER* compbuilder /**< builder */
782  )
783 {
784  assert(orggraph && compbuilder);
785 
786  printf("number of separator nodes=%d, nodes: \n", compbuilder->nsepatterms);
787 
788  for( int i = 0; i < compbuilder->nsepatterms; i++ )
789  {
790  const int sepaterm = compbuilder->sepaterms[i];
791  assert(graph_knot_isInRange(orggraph, sepaterm));
792 
793  graph_knot_printInfo(orggraph, sepaterm);
794  }
795 }
796 
797 
798 /** builds extended subgraph with SD weighted edges between terminal-separator nodes */
800  SCIP* scip, /**< SCIP data structure */
801  GRAPH* orggraph, /**< graph data structure */
802  EXTPERMA* extperma, /**< extension data */
803  TERMCOMP* termcomp /**< component */
804  )
805 {
806  const SCIP_Bool useSd = TRUE;
807 
808  assert(scip && orggraph && extperma && termcomp);
809 
810  subgraphIdentify(scip, orggraph, termcomp);
811  SCIP_CALL( subgraphBuild(scip, orggraph, useSd, extperma, termcomp) );
812 
813  return SCIP_OKAY;
814 }
815 
816 
817 /** builds extended subgraph with BLOCKED weighted edges between terminal-separator nodes */
819  SCIP* scip, /**< SCIP data structure */
820  GRAPH* orggraph, /**< graph data structure */
821  TERMCOMP* termcomp /**< component */
822  )
823 {
824  const SCIP_Bool useSd = FALSE;
825 
826  assert(scip && orggraph && termcomp);
827 
828  subgraphIdentify(scip, orggraph, termcomp);
829  SCIP_CALL( subgraphBuild(scip, orggraph, useSd, NULL, termcomp) );
830 
831  return SCIP_OKAY;
832 }
833 
834 
835 /** Changes weights between terminal separator nodes to bottlenecks.
836  * NOTE: also computes the bottleneck distances, so potentially expensive! */
838  SCIP* scip, /**< SCIP data structure */
839  GRAPH* g, /**< graph data structure */
840  TERMCOMP* termcomp, /**< component */
841  SCIP_Bool* success /**< pointer to store whether computation was successful (OUT) */
842  )
843 {
844  GRAPH* subgraph = termcomp->subgraph;
845  PATH* mst;
846  COMPBUILDER* const builder = termcomp->builder;
847  const int* const sepaterms = builder->sepaterms;
848  const int* const nodemap_orgToSub = termcomp->nodemap_orgToSub;
849  const int* const nodemap_subToOrg = termcomp->nodemap_subToOrg;
850  const int* const nodes_mark = termcomp->nodes_mark;
851  const int nsepaterms = builder->nsepatterms;
852  const int nnodes = graph_get_nNodes(g);
853  const int mstroot = sepaterms[0];
854 
855  *success = TRUE;
856 
857  /* mark the anti-component */
858  for( int i = 0; i < nnodes; i++ )
859  g->mark[i] = (nodes_mark[i] != MARK_SUBNODE);
860 
861  SCIPallocBufferArray(scip, &mst, nnodes);
862  graph_path_exec(scip, g, MST_MODE, mstroot, g->cost, mst);
863 
864  for( int i = 0; i < nsepaterms; i++ )
865  {
866  const int term = sepaterms[i];
867  if( term != mstroot && mst[term].edge == -1 )
868  {
869  SCIPfreeBufferArray(scip, &mst);
870  *success = FALSE;
871  graph_mark(g);
872  return SCIP_OKAY;
873  }
874  }
875 
876 #ifndef NDEBUG
877  for( int i = 0; i < g->knots; i++ )
878  {
879  if( Is_term(g->term[i]) && g->mark[i] && i != mstroot )
880  {
881  assert(mst[i].edge != -1);
882  }
883  }
884 #endif
885 
886  for( int i = 0; i < subgraph->knots; i++ )
887  subgraph->mark[i] = MARK_NONACTIVE;
888 
889  for( int i = 0; i < nsepaterms; i++ )
890  {
891  const int term = sepaterms[i];
892  const int term_sub = nodemap_orgToSub[term];
893  assert(graph_knot_isInRange(subgraph, term_sub));
894  assert(subgraph->mark[term_sub] == MARK_NONACTIVE);
895 
896  subgraph->mark[term_sub] = MARK_SEPARATOR;
897  }
898 
899  for( int i = 0; i < nsepaterms; i++ )
900  {
901  const int term = sepaterms[i];
902  const int term_sub = nodemap_orgToSub[term];
903 
904  for( int e = subgraph->outbeg[term_sub]; e != EAT_LAST; e = subgraph->oeat[e] )
905  {
906  const int head_sub = subgraph->head[e];
907 
908  if( head_sub > term_sub && subgraph->mark[head_sub] == MARK_SEPARATOR )
909  {
910  const int head = nodemap_subToOrg[head_sub];
911  const SCIP_Real b = getExtBottleneckDist(g, term, head, term_sub, head_sub, mstroot, mst,
912  termcomp->subsolbottleneck, builder->nodes_bdist);
913 
914  subgraph->cost[e] = b;
915  subgraph->cost[flipedge(e)] = b;
916 
917  SCIPdebugMessage("%d->%d setting b=%f \n", term_sub, head_sub, b);
918  }
919  }
920  }
921 
922  graph_mark(g);
923  SCIPfreeBufferArray(scip, &mst);
924 
925  return SCIP_OKAY;
926 }
927 
928 
929 /** changes weights between terminal separator nodes to original costs (or BLOCKED) */
931  const GRAPH* orggraph, /**< graph data structure */
932  TERMCOMP* termcomp /**< component */
933  )
934 {
935  GRAPH* subgraph = termcomp->subgraph;
936  const COMPBUILDER* const builder = termcomp->builder;
937  const int* const edgemap_subToOrg = termcomp->edgemap_subToOrg;
938  const int nsepaterms = builder->nsepatterms;
939 
940  assert(subgraph && orggraph);
941  assert(nsepaterms > 0);
942 
943  /* NOTE: we exploit the fact that the first nodes of subgraph are the separator terminals */
944 
945  for( int subsepaterm = 0; subsepaterm < nsepaterms; subsepaterm++ )
946  {
947  for( int sube = subgraph->outbeg[subsepaterm]; sube != EAT_LAST; sube = subgraph->oeat[sube] )
948  {
949  const int subhead = subgraph->head[sube];
950  if( subhead < nsepaterms )
951  {
952  const int orge = edgemap_subToOrg[sube];
953 
954  assert(Is_term(subgraph->term[subsepaterm]));
955  assert(Is_term(subgraph->term[subhead]));
956 
957  if( orge >= 0 )
958  {
959  assert(Is_term(orggraph->term[orggraph->head[orge]]));
960  assert(Is_term(orggraph->term[orggraph->tail[orge]]));
961 
962  subgraph->cost[sube] = orggraph->cost[orge];
963  }
964  else
965  {
966  subgraph->cost[sube] = BLOCKED;
967  }
968  }
969  }
970  }
971 
972 #ifndef NDEBUG
973  for( int e = 0; e < subgraph->edges; e += 2 )
974  {
975  assert(EQ(subgraph->cost[e], subgraph->cost[e + 1]));
976  }
977 #endif
978 }
979 
980 
981 /** initializes */
983  SCIP* scip, /**< SCIP data structure */
984  const GRAPH* g, /**< graph data structure */
985  COMPBUILDER* builder, /**< initializer */
986  TERMCOMP** termcomp /**< to initialize */
987  )
988 {
989  TERMCOMP* comp;
990 
991  SCIP_CALL( SCIPallocMemory(scip, termcomp) );
992  comp = *termcomp;
993 
994  comp->builder = builder;
995  comp->subgraph = NULL;
996  comp->subsolution = NULL;
997  comp->edgemap_subToOrg = NULL;
998  comp->nodemap_subToOrg = NULL;
999  comp->bfsqueue = NULL;
1000  comp->subnedges = -1;
1001  comp->subnnodes = -1;
1002  comp->subprimalobj = -FARAWAY;
1003  comp->subsolbottleneck = NULL;
1004 
1005  SCIP_CALL( SCIPallocMemoryArray(scip, &(comp->nodes_mark), g->knots) );
1006  SCIP_CALL( SCIPallocMemoryArray(scip, &(comp->nodemap_orgToSub), g->knots) );
1007 
1008 #ifndef NDEBUG
1009  for( int i = 0; i < g->knots; i++ )
1010  comp->nodemap_orgToSub[i] = UNKNOWN;
1011 #endif
1012 
1013  return SCIP_OKAY;
1014 }
1015 
1016 
1017 
1018 /** initializes tree bottleneck for terminal component */
1020  SCIP* scip, /**< SCIP data structure */
1021  const int* subsoledges, /**< solution tree to be represented */
1022  TERMCOMP* termcomp /**< to initialize for */
1023  )
1024 {
1025  const GRAPH* subgraph;
1026 
1027  assert(scip && subsoledges && termcomp);
1028  assert(!termcomp->subsolbottleneck);
1029 
1030  subgraph = termcomp->subgraph;
1031 
1032  SCIP_CALL( tbottleneckInit(scip, subgraph, subsoledges, &(termcomp->subsolbottleneck)) );
1033 
1034  return SCIP_OKAY;
1035 }
1036 
1037 
1038 /** frees */
1040  SCIP* scip, /**< SCIP data structure */
1041  TERMCOMP** termcomp /**< to initialize */
1042  )
1043 {
1044  TERMCOMP* comp;
1045  comp = *termcomp;
1046 
1047  assert(comp->nodemap_subToOrg && comp->edgemap_subToOrg);
1048 
1049 
1050  if( comp->subsolbottleneck )
1051  tbottleneckFree(scip, &(comp->subsolbottleneck));
1052  StpVecFree(scip, comp->bfsqueue);
1053 
1054  SCIPfreeMemoryArrayNull(scip, &(comp->subsolution));
1055  SCIPfreeMemoryArray(scip, &(comp->nodemap_subToOrg));
1056  SCIPfreeMemoryArray(scip, &(comp->edgemap_subToOrg));
1057 
1058  if( comp->subgraph )
1059  {
1060  graph_free(scip, &(comp->subgraph), TRUE);
1061  }
1062 
1063  SCIPfreeMemoryArray(scip, &(comp->nodemap_orgToSub));
1064  SCIPfreeMemoryArray(scip, &(comp->nodes_mark));
1065 
1066  SCIPfreeMemory(scip, termcomp);
1067 }
1068 
1069 
1070 
1071 /** gets next separator component */
1073  SCIP* scip, /**< SCIP data structure */
1074  const GRAPH* g, /**< graph data structure */
1075  TERMSEPAS* termsepas, /**< terminal separator store */
1076  COMPBUILDER* builder, /**< builder */
1077  SCIP_Bool* compWasFound /**< was new component found? */
1078  )
1079 {
1080  const int* sepaterms = NULL;
1081  int sinkterm;
1082  int nsinknodes;
1083  int nsepaterms = builder->nsepatterms;
1084 
1085  *compWasFound = FALSE;
1086 
1087  assert(nsepaterms >= 2);
1088  assert(builder->componentnumber >= 0);
1089 
1090  if( builder->componentnumber > builder->maxncompchecks )
1091  {
1092  SCIPdebugMessage("maximum number of components reached, stopping sepaDualAscent reductions \n");
1093  return;
1094  }
1095 
1096  for( ; nsepaterms <= builder->maxsepasize; nsepaterms++ )
1097  {
1098  sepaterms = mincut_termsepasGetNext(nsepaterms, termsepas, &sinkterm, &nsinknodes);
1099 
1100  if( sepaterms != NULL )
1101  {
1102  *compWasFound = TRUE;
1103  break;
1104  }
1105  }
1106 
1107  if( *compWasFound )
1108  {
1109  const int nsourcenodes = builder->ngraphnodes - nsinknodes;
1110 
1111  assert(nsinknodes > 0);
1112  assert(nsourcenodes > 0);
1113 
1114  builder->nsepatterms = nsepaterms;
1115  builder->sepaterms = sepaterms;
1116 
1117  /* NOTE: we want to take the smaller component if possible */
1118  if( nsinknodes > nsourcenodes && !builder->rootcompIsProcessed )
1119  {
1120  const int sourceterm = mincut_termsepasGetSource(termsepas);
1121  assert(graph_knot_isInRange(g, sourceterm));
1122 
1123  builder->sourceterm = sourceterm;
1124  builder->ncomponentnodes = nsourcenodes;
1125  /* NOTE: even if the component turns out to not be promising, we only want to check it once */
1126  builder->rootcompIsProcessed = TRUE;
1127  }
1128  else
1129  {
1130  builder->sourceterm = sinkterm;
1131  builder->ncomponentnodes = nsinknodes;
1132  }
1133 
1134  SCIPdebugMessage("selecting component: source=%d, |S|=%d, |V'|=%d \n", builder->sourceterm, nsepaterms, builder->ncomponentnodes);
1135  }
1136  else
1137  {
1138  SCIPdebugMessage("no further components available, stopping sepaDualAscent reductions \n");
1139  }
1140 }
#define STP_Vectype(type)
Definition: stpvector.h:44
int graph_get_nEdges(const GRAPH *)
Definition: graph_stats.c:548
#define StpVecGetSize(vec)
Definition: stpvector.h:139
static void subgraphIdentify(SCIP *scip, GRAPH *g, TERMCOMP *termcomp)
SCIP_RETCODE reduce_termcompInitTbottleneck(SCIP *scip, const int *subsoledges, TERMCOMP *termcomp)
int *RESTRICT head
Definition: graphdefs.h:224
int source
Definition: graphdefs.h:195
#define SCIPfreeMemoryArrayNull(scip, ptr)
Definition: scip_mem.h:72
#define Is_term(a)
Definition: graphdefs.h:102
SCIP_Real extreduce_distDataGetSdDouble(SCIP *, const GRAPH *, int, int, DISTDATA *)
SCIP_RETCODE reduce_compbuilderInit(SCIP *scip, const GRAPH *g, COMPBUILDER **compbuilder)
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
signed int edge
Definition: graphdefs.h:287
void graph_edge_addBi(SCIP *, GRAPH *, int, int, double)
SCIP_RETCODE reduce_termcompInit(SCIP *scip, const GRAPH *g, COMPBUILDER *builder, TERMCOMP **termcomp)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
void graph_path_exec(SCIP *, const GRAPH *, int, int, const SCIP_Real *, PATH *)
Definition: graph_path.c:541
SCIP_RETCODE reduce_termcompBuildSubgraph(SCIP *scip, GRAPH *orggraph, TERMCOMP *termcomp)
void reduce_termsepaGetNextComp(SCIP *scip, const GRAPH *g, TERMSEPAS *termsepas, COMPBUILDER *builder, SCIP_Bool *compWasFound)
void graph_mark(GRAPH *)
Definition: graph_base.c:1130
includes methods for Steiner tree problem solutions
void graph_knot_printInfo(const GRAPH *, int)
Definition: graph_stats.c:151
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: graph_base.c:796
SCIP_RETCODE graph_path_init(SCIP *, GRAPH *)
Definition: graph_path.c:480
#define FALSE
Definition: def.h:87
int mincut_termsepasGetSource(const TERMSEPAS *termsepas)
Definition: mincut.c:2232
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
includes various files containing graph methods used for Steiner tree problems
#define EAT_LAST
Definition: graphdefs.h:58
#define StpVecReserve(scip, vec, _size_)
Definition: stpvector.h:182
#define SCIPdebugMessage
Definition: pub_message.h:87
#define FARAWAY
Definition: graphdefs.h:89
#define MARK_SUBNODE
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
void graph_printInfo(const GRAPH *)
Definition: graph_stats.c:299
header only, simple implementation of an STL like vector
int *RESTRICT mark
Definition: graphdefs.h:198
void reduce_compbuilderFree(SCIP *scip, COMPBUILDER **compbuilder)
int *RESTRICT oeat
Definition: graphdefs.h:231
This file implements extended reduction techniques for several Steiner problems.
#define LE(a, b)
Definition: portab.h:83
#define GE(a, b)
Definition: portab.h:85
int stp_type
Definition: graphdefs.h:264
void graph_get_nVET(const GRAPH *, int *, int *, int *)
Definition: graph_stats.c:571
#define MARK_NONACTIVE
int *RESTRICT grad
Definition: graphdefs.h:201
void graph_knot_add(GRAPH *, int)
Definition: graph_node.c:64
#define NULL
Definition: lpi_spx1.cpp:155
#define MARK_SEPARATOR
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
#define MST_MODE
Definition: graphdefs.h:98
SCIP_RETCODE reduce_termcompBuildSubgraphWithSds(SCIP *scip, GRAPH *orggraph, EXTPERMA *extperma, TERMCOMP *termcomp)
SCIP_Real reduce_compbuilderGetSubNodesRatio(const COMPBUILDER *compbuilder)
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define GT(a, b)
Definition: portab.h:84
#define SCIP_Bool
Definition: def.h:84
static SCIP_RETCODE tbottleneckInit(SCIP *scip, const GRAPH *g, const int *soledges, TBOTTLENECK **tbottleneck)
int *RESTRICT tail
Definition: graphdefs.h:223
#define flipedge(edge)
Definition: graphdefs.h:84
#define MAX(x, y)
Definition: tclique_def.h:83
static SCIP_Real getExtBottleneckDist(const GRAPH *g, int term1, int term2, int term1_sub, int term2_sub, int mstsource, PATH *mst, TBOTTLENECK *subsolbottleneck, SCIP_Real *RESTRICT mstsdist)
#define STP_TERM
Definition: graphdefs.h:61
int *RESTRICT term
Definition: graphdefs.h:196
void graph_edge_printInfo(const GRAPH *, int)
Definition: graph_stats.c:133
#define CONNECT
Definition: graphdefs.h:87
Portable definitions.
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
SCIP_RETCODE reduce_termcompChangeSubgraphToBottleneck(SCIP *scip, GRAPH *g, TERMCOMP *termcomp, SCIP_Bool *success)
SCIP_VAR ** b
Definition: circlepacking.c:56
void reduce_termcompFree(SCIP *scip, TERMCOMP **termcomp)
SCIP_Real * cost
Definition: graphdefs.h:221
static SCIP_RETCODE subgraphBuild(SCIP *scip, const GRAPH *orggraph, SCIP_Bool useSd, EXTPERMA *extperma, TERMCOMP *termcomp)
const int * mincut_termsepasGetNext(int sepasize, TERMSEPAS *termsepas, int *sinkterm, int *nsinknodes)
Definition: mincut.c:2179
void reduce_termcompChangeSubgraphToOrgCosts(const GRAPH *orggraph, TERMCOMP *termcomp)
#define SCIP_Real
Definition: def.h:177
#define BLOCKED
Definition: graphdefs.h:90
#define StpVecFree(scip, vec)
Definition: stpvector.h:149
int *RESTRICT outbeg
Definition: graphdefs.h:204
int edges
Definition: graphdefs.h:219
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
struct tree_bottleneck_node TBNODE
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
#define UNKNOWN
Definition: sepa_mcf.c:4095
static void tbottleneckGetMax(int node1, int node2, TBOTTLENECK *tbottleneck, SCIP_Real *maxlength)
#define nnodes
Definition: gastrans.c:65
TBOTTLENECK * subsolbottleneck
Definition: termsepadefs.h:65
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:163
SCIP_RETCODE graph_init(SCIP *, GRAPH **, int, int, int)
Definition: graph_base.c:607
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:123
SCIP_Bool solstp_isValid(SCIP *scip, const GRAPH *graph, const int *result)
Definition: solstp.c:1650
Minimum cut routines for Steiner problems.
void reduce_compbuilderPrintSeparators(const GRAPH *orggraph, const COMPBUILDER *compbuilder)
includes various reduction methods for Steiner tree problems
SCIP callable library.
static void tbottleneckFree(SCIP *scip, TBOTTLENECK **tbottleneck)