# SCIP

Solving Constraint Integer Programs

reduce_sdgraph.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 /* */
7 /* fuer Informationstechnik Berlin */
8 /* */
10 /* */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file reduce_sdgraph.c
17  * @brief bottleneck distance graph methods for Steiner tree reductions
18  * @author Daniel Rehfeldt
19  *
20  * This file implements methods for Steiner tree problem special distance (bottleneck Steiner distance) graph.
21  * This graph is the (complete) distance graph on the terminal vertex set.
22  * Note that the complete graph is in general not built.
23  *
24  * A list of all interface methods can be found in reduce.h.
25  *
26  */
27
28 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
29
30 //#define SCIP_DEBUG
31 #include "reduce.h"
32 #include "misc_stp.h"
33 #include "portab.h"
34
35 #define STP_SDQUERYFULL_MAXTERMS 200
36
37 //#define SDQUERYFULL_HALFMATRIX
38
39
40 /*
41  * STRUCTS
42  */
43
44
45 /** see reduce.h */
47 {
48  GRAPH* distgraph; /**< (complete) distance graph */
49  PATH* sdmst; /**< MST on sdgraph */
50  SCIP_Real* mstcosts; /**< maximum MST edge costs in descending order */
51  SCIP_Real* mstsdist; /**< helper array; each entry needs to -1.0; of size nnodesorg */
52  int* nodemapOrgToDist; /**< node mapping from original graph to distance graph */
53  STP_Bool* halfedge_isInMst; /**< signifies whether edge of original graph is part of MST
54  NOTE: operates on edges / 2! */
55  SCIP_Real mstmaxcost; /**< maximum edge cost */
56  int nnodesorg; /**< number of nodes of original graph */
57  int nedgesorg; /**< number of edges of original graph */
59  SCIP_Bool edgemarkReady; /**< edge mark available? */
60  SCIP_Bool usingRMQ; /**< is RMQ being used? */
61  /* Full query data: */
62  SCIP_Real* fullq_dists; /**< of size |T| * |T - 1| or NULL */
63  int* fullq_nodeToIdx; /**< maps nodes of original graphs to full index; of size |V| */
64 #ifdef SDQUERYFULL_HALFMATRIX
65  int* fullq_IdxToStart; /**< start for dists */
66 #endif
67  int fullq_size; /**< |T - 1| * |T - 1| or |T - 1| * |T - 1| / 2 */
68  int fullq_dimension; /**< |T - 1| */
69  /* RMQ query data: */
70  int* rmq_nodeToRmqEntry; /**< of size |V| */
71  SCIP_Real* rmq_sparseTable; /**< of size rmq_log * 2 |T| */
72  SCIP_Real* rmq_edgecosts; /**< of size 2|T| - 3 */
73  int rmq_loglength; /**< log2(2|T| - 3) */
74 };
75
76
77
78 /** Simple binary tree node.
79  * Value of UNKNOWN signifies that child does not exist. */
80 typedef struct binary_tree_node
81 {
82  int child1; /**< first child */
83  int child2; /**< second child */
84 } BINARYNODE;
85
86
87 /** data needed for building the RMQ based SD query structures */
89 {
90  BINARYNODE* lcatree; /**< binary tree used for LCA computation; of size |T| - 1 */
91  SCIP_Real* lcatree_costs; /**< cost per node of binary tree used for LCA computation; of size 2|T| - 3 */
92  int* termToLcatreeNode; /**< node mapping from terminals (SD graph nodes) to binary LCA tree nodes */
93  int* lcatreeNodeToInds; /**< mapping from binary LCA tree nodes to RMQ/Full indices; of size |T| - 1 */
94 } LCABUILDER;
95
96
97 /*
98  * local methods
99  */
100
101 static const int deBruijnBits[32] =
102 {
103  0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
104  8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
105 };
106
107 /** returns the floor of log2(number) */
108 static inline
110  uint32_t number /**< number */
111  )
112 {
113  uint32_t v = number;
114  v |= v >> 1;
115  v |= v >> 2;
116  v |= v >> 4;
117  v |= v >> 8;
118  v |= v >> 16;
119
120  assert(deBruijnBits[(uint32_t)(v * 0x07C4ACDDU) >> 27] == (int) log2(number));
121
122  return deBruijnBits[(uint32_t)(v * 0x07C4ACDDU) >> 27];
123 }
124
125 #ifdef SCIP_DEBUG
126 /* prints the tree */
127 static
128 void lcabuilderPrintLcaTree(
129  const SDGRAPH* sdgraph, /**< the SD graph */
130  const LCABUILDER* lcabuilder /**< the builder */
131 )
132 {
133  const BINARYNODE* const lcatree = lcabuilder->lcatree;
134  const SCIP_Real* const lcatree_costs = lcabuilder->lcatree_costs;
135  const int* const termToLcatreeNode = lcabuilder->termToLcatreeNode;
136  const int nterms = graph_get_nNodes(sdgraph->distgraph);
137
138  assert(termToLcatreeNode && lcatree && lcatree_costs && lcabuilder->lcatree && lcabuilder->termToLcatreeNode);
139
140  printf("\n SD QUERY LCA TREE: \n");
141
142  for( int i = 0; i < nterms - 1; i++ )
143  {
144  const SCIP_Real cost = lcatree_costs[i];
145  printf("node %d: \n", i);
146
147  printf("...cost=%f \n", cost);
148  printf("...child1=%d, child2=%d \n", lcatree[i].child1, lcatree[i].child2);
149  }
150
151  printf("\n");
152
153  for( int i = 0; i < nterms; i++ )
154  {
155  printf("term %d attached to tree node %d \n", i, termToLcatreeNode[i]);
156  }
157 }
158 #endif
159
160
161
162 /** Gets special distance (i.e. bottleneck distance) from graph.
163  * Corresponds to bottleneck length of path between term1 and term2 on distance graph */
164 static inline
166  int term1, /**< terminal 1 */
167  int term2, /**< terminal 2 */
168  SDGRAPH* sdgraph /**< the SD graph */
169 )
170 {
171  SCIP_Real* RESTRICT mstsdist = sdgraph->mstsdist;
172  const GRAPH* distgraph = sdgraph->distgraph;
173  const PATH* sdmst = sdgraph->sdmst;
174  const int* nodesid = sdgraph->nodemapOrgToDist;
175  const int sdnode1 = nodesid[term1];
176  const int sdnode2 = nodesid[term2];
177  SCIP_Real sdist = 0.0;
178  int tempnode = sdnode1;
179
180  assert(sdnode1 != sdnode2);
181  assert(distgraph->source == 0);
182
183  mstsdist[tempnode] = 0.0;
184
185  /* not at root? */
186  while( tempnode != 0 )
187  {
188  const int ne = sdmst[tempnode].edge;
189
191  tempnode = distgraph->tail[ne];
192
193  if( distgraph->cost[ne] > sdist )
194  sdist = distgraph->cost[ne];
195
196  mstsdist[tempnode] = sdist;
197  if( tempnode == sdnode2 )
198  break;
199  }
200
202  if( tempnode == sdnode2 )
203  {
204  tempnode = 0;
205  }
206  else
207  {
208  tempnode = sdnode2;
209  sdist = 0.0;
210  }
211
212  while( tempnode != 0 )
213  {
214  const int ne = sdmst[tempnode].edge;
215  tempnode = distgraph->tail[ne];
216
217  if( distgraph->cost[ne] > sdist )
218  sdist = distgraph->cost[ne];
219
220  assert(GE(mstsdist[tempnode], 0.0) == (mstsdist[tempnode] > -0.5));
221
223  if( mstsdist[tempnode] > -0.5 )
224  {
225  if( mstsdist[tempnode] > sdist )
226  sdist = mstsdist[tempnode];
227  break;
228  }
229
230 #ifndef NDEBUG
231  assert(EQ(mstsdist[tempnode], -1.0));
232
233  if( tempnode == 0 )
234  {
235  assert(sdnode1 == 0);
236  }
237 #endif
238  }
239
240  /* restore mstsdist */
241  tempnode = sdnode1;
242  mstsdist[tempnode] = -1.0;
243  while( tempnode != 0 )
244  {
245  const int ne = sdmst[tempnode].edge;
246  tempnode = distgraph->tail[ne];
247  mstsdist[tempnode] = -1.0;
248  if( tempnode == sdnode2 )
249  break;
250  }
251
252  assert(GE(sdist, 0.0));
253
254  return sdist;
255 }
256
257 /** gets length of RMQ array */
258 static inline
260  const SDGRAPH* sdgraph /**< the SD graph */
261  )
262 {
263  const int nterms = graph_get_nNodes(sdgraph->distgraph);
264  assert(nterms >= 2);
265  assert((2 * nterms - 3) >= 1);
266
267  return (2 * nterms - 3);
268 }
269
270 /** initializes LCA builder */
271 static
273  SCIP* scip, /**< SCIP */
274  const SDGRAPH* sdgraph, /**< the SD graph */
275  LCABUILDER** lcabuilder /**< the builder */
276 )
277 {
278  LCABUILDER* rmqb;
279  const int nterms = graph_get_nNodes(sdgraph->distgraph);
280
281  assert(nterms >= 2);
282
283  SCIP_CALL( SCIPallocMemory(scip, lcabuilder) );
284  rmqb = *lcabuilder;
285
286  SCIP_CALL( SCIPallocMemoryArray(scip, &(rmqb->lcatree), nterms - 1) );
287  SCIP_CALL( SCIPallocMemoryArray(scip, &(rmqb->lcatree_costs), nterms - 1) );
288  SCIP_CALL( SCIPallocMemoryArray(scip, &(rmqb->lcatreeNodeToInds), nterms - 1) );
289  SCIP_CALL( SCIPallocMemoryArray(scip, &(rmqb->termToLcatreeNode), nterms) );
290
291  for( int i = 0; i < nterms - 1; i++ )
292  {
293  rmqb->lcatree[i].child1 = UNKNOWN;
294  rmqb->lcatree[i].child2 = UNKNOWN;
295  }
296
297 #ifndef NDEBUG
298  for( int i = 0; i < nterms - 1; i++ )
299  {
300  rmqb->lcatreeNodeToInds[i] = -1;
301  }
302
303  for( int i = 0; i < nterms; i++ )
304  {
305  rmqb->termToLcatreeNode[i] = -1;
306  }
307 #endif
308
309  return SCIP_OKAY;
310 }
311
312
313 /** frees LCA builder */
314 static
316  SCIP* scip, /**< SCIP */
317  LCABUILDER** lcabuilder /**< the builder */
318 )
319 {
320  LCABUILDER* rmqb = *lcabuilder;
321
322  SCIPfreeMemoryArray(scip, &(rmqb->termToLcatreeNode));
323  SCIPfreeMemoryArray(scip, &(rmqb->lcatreeNodeToInds));
324  SCIPfreeMemoryArray(scip, &(rmqb->lcatree_costs));
325  SCIPfreeMemoryArray(scip, &(rmqb->lcatree));
326
327  SCIPfreeMemory(scip, lcabuilder);
328 }
329
330
331 /* attaches node */
332 static inline
334  const GRAPH* distgraph, /**< SD distance graph */
335  int parentposition, /**< position of parent (edge) in binary tree array */
336  int graphnode, /**< graph node to attach */
337  LCABUILDER* lcabuilder, /**< the builder */
338  UF* uf /**< union-find data structure */
339 )
340 {
341  BINARYNODE* RESTRICT lcatree = lcabuilder->lcatree;
342  int* RESTRICT termToLcatreeNode = lcabuilder->termToLcatreeNode;
343  const int childidentifier = SCIPStpunionfindFind(uf, graphnode);
344  const int nterms = graph_get_nNodes(distgraph);
345
346  graph_knot_isInRange(distgraph, graphnode);
347  assert(0 <= parentposition && parentposition < nterms - 1);
348
349  if( childidentifier == graphnode )
350  {
351  assert(graphnode < nterms);
352
353  /* just need to set a pointer in this case */
354  termToLcatreeNode[graphnode] = parentposition;
355  }
356  else
357  {
358  /* child identifier marks an edge position */
359  const int childosition = childidentifier - nterms;
360  assert(0 <= childosition && childosition < nterms - 1);
361
362  if( lcatree[parentposition].child1 == UNKNOWN )
363  {
364  lcatree[parentposition].child1 = childosition;
365  }
366  else
367  {
368  assert(lcatree[parentposition].child2 == UNKNOWN);
369  lcatree[parentposition].child2 = childosition;
370  }
371  }
372
373  SCIPStpunionfindUnion(uf, nterms + parentposition, childidentifier, FALSE);
374 }
375
376
377 /** builds binary tree for LCA computation */
378 static
380  SCIP* scip, /**< SCIP */
381  SDGRAPH* sdgraph, /**< the SD graph */
382  LCABUILDER* lcabuilder /**< the builder */
383 )
384 {
385  const GRAPH* const distgraph = sdgraph->distgraph;
386  const PATH* const sdmst = sdgraph->sdmst;
387  int* RESTRICT mstedges_id;
388  SCIP_Real* RESTRICT lcatree_costs = lcabuilder->lcatree_costs;
389  const int nterms = graph_get_nNodes(sdgraph->distgraph);
390  int edgecount = 0;
391
392  /* Union-find: entries 0,...,|T|-1 are used for terminals,
393  * entries |T|,...,2 |T| - 2 are used for the ordered MST edges */
394  UF uf;
395
396  SCIP_CALL( SCIPStpunionfindInit(scip, &uf, 2 * nterms - 1) );
397  SCIP_CALL( SCIPallocBufferArray(scip, &mstedges_id, nterms - 1) );
398
399  /* fill the MST edge arrays */
400  for( int i = 0; i < nterms; i++ )
401  {
402  const int edge = sdmst[i].edge;
403
404  if( edge >= 0 )
405  {
406  assert(graph_edge_isInRange(sdgraph->distgraph, edge));
407  assert(EQ(sdmst[i].dist, sdgraph->distgraph->cost[edge]));
408
409  mstedges_id[edgecount] = edge;
410  lcatree_costs[edgecount++] = sdmst[i].dist;
411  }
412  }
413  assert(edgecount == nterms - 1);
414  assert(edgecount >= 1);
415
416  /* build the actual tree */
417
418  SCIPsortDownRealInt(lcatree_costs, mstedges_id, edgecount);
419
420  for( int i = edgecount - 1; i >= 0; i-- )
421  {
422  const int edge = mstedges_id[i];
423  const int tail = distgraph->tail[edge];
425
426  sdqueryAttachBinaryTreeNode(distgraph, i, tail, lcabuilder, &uf);
427  sdqueryAttachBinaryTreeNode(distgraph, i, head, lcabuilder, &uf);
428  }
429
430  SCIPfreeBufferArray(scip, &mstedges_id);
431  SCIPStpunionfindFreeMembers(scip, &uf);
432
433 #ifdef SCIP_DEBUG
434  lcabuilderPrintLcaTree(sdgraph, lcabuilder);
435 #endif
436
437  return SCIP_OKAY;
438 }
439
440
441 /** Builds RMQ by DFS.
442  * NOTE: recursive method. */
443 static
445  int root, /**< root node (LCA tree position) */
446  const BINARYNODE* lcatree, /**< tree */
447  const SCIP_Real* lcatree_costs, /**< LCA tree costs */
448  SCIP_Real* rmq_edgecosts, /**< edge costs for RMQ */
449  int* lcatreeNodeToRmq, /**< mapping from binary LCA tree nodes to RMQ indices */
450  int* rmq_count /**< current position */
451 )
452 {
453  const BINARYNODE bnode = lcatree[root];
454
455  assert(lcatreeNodeToRmq[root] == -1);
456  lcatreeNodeToRmq[root] = *rmq_count;
457  rmq_edgecosts[(*rmq_count)++] = lcatree_costs[root];
458
459  SCIPdebugMessage("checking node %d \n", root);
460
461  if( bnode.child1 == UNKNOWN )
462  {
463  assert(bnode.child2 == UNKNOWN );
464  return;
465  }
466
467  sdqueryRmqDfs(bnode.child1, lcatree, lcatree_costs, rmq_edgecosts, lcatreeNodeToRmq, rmq_count);
468  rmq_edgecosts[(*rmq_count)++] = lcatree_costs[root];
469
470  if( bnode.child2 == UNKNOWN )
471  {
472  return;
473  }
474
475  sdqueryRmqDfs(bnode.child2, lcatree, lcatree_costs, rmq_edgecosts, lcatreeNodeToRmq, rmq_count);
476  rmq_edgecosts[(*rmq_count)++] = lcatree_costs[root];
477 }
478
479
480 /** builds RMQ sparse table */
481 static
483  SDGRAPH* sdgraph /**< the SD graph */
484  )
485 {
486  const SCIP_Real* const rmq_edgecosts = sdgraph->rmq_edgecosts;
487  SCIP_Real* const rmq_sparsetable = sdgraph->rmq_sparseTable;
488  const int rmq_length = sdqueryGetRmqLength(sdgraph);
489  const int logt = sdgraph->rmq_loglength;
490  const int rowlength = logt + 1;
491
492  SCIPdebugMessage("\n building sparse table with floor(log2(|T|))=%d \n", logt);
493
494  /* base computation */
495  for( int i = 0; i < rmq_length; i++ )
496  {
497  const int start_i = i * rowlength;
498  rmq_sparsetable[start_i] = rmq_edgecosts[i];
499  SCIPdebugMessage("max[%d, %d)=%f \n", i, i + 1, rmq_sparsetable[start_i]);
500  }
501
502  /* main (dynamic programming) computation */
503  for( int j = 1; j <= logt; j++ )
504  {
505  const int shift = 1 << (unsigned int) (j - 1);
506
507  SCIPdebugMessage("round %d (2^%d), shift=%d \n", j, j, shift);
508
509  for( int i = 0; i < rmq_length; i++ )
510  {
511  const int start_i = i * rowlength;
512
513  if( i + shift < rmq_length )
514  {
515  /* position for interval [i, i + 2^(j-1)) */
516  const int pos_iprev = start_i + j - 1;
517  const int start_shift = (i + shift) * rowlength;
518  /* position for interval [(i+2^(j-1)), (i 2^(j-1)) + 2^(j-1)) */
519  const int pos_shift = start_shift + j - 1 ;
520
521  assert(pos_shift < rmq_length * rowlength);
522
523  rmq_sparsetable[start_i + j] = MAX(rmq_sparsetable[pos_iprev], rmq_sparsetable[pos_shift]);
524  }
525  else
526  {
527  /* NOTE: basically debug marker */
528  rmq_sparsetable[start_i + j] = FARAWAY;
529  }
530
531  SCIPdebugMessage("max[%d, 2^%d)=%f \n", i, j, rmq_sparsetable[start_i + j]);
532  }
533  }
534 }
535
536
537 /** builds RMQ, including sparse table representation */
538 static
540  SDGRAPH* sdgraph, /**< the SD graph */
541  LCABUILDER* lcabuilder /**< the builder */
542 )
543 {
544  const BINARYNODE* const lcatree = lcabuilder->lcatree;
545  const SCIP_Real* const lcatree_costs = lcabuilder->lcatree_costs;
546  SCIP_Real* const rmq_edgecosts = sdgraph->rmq_edgecosts;
547  int* lcatreeNodeToRmq = lcabuilder->lcatreeNodeToInds;
548  int rmq_length = 0;
549
550  /* NOTE: because we have sorted in ascending edge cost order, the root is 0 */
551
552  sdqueryRmqDfs(0, lcatree, lcatree_costs, rmq_edgecosts, lcatreeNodeToRmq, &rmq_length);
553  assert(rmq_length == sdqueryGetRmqLength(sdgraph));
554
555 #ifdef SCIP_DEBUG
556  for( int i = 0; i < rmq_length; i++ )
557  {
558  printf("i=%d ocst=%f \n", i, rmq_edgecosts[i]);
559  }
560 #endif
561
563 }
564
565
566 /** builds RMQ map */
567 static
569  const GRAPH* g, /**< graph */
570  const LCABUILDER* lcabuilder, /**< the builder */
571  SDGRAPH* sdgraph /**< the SD graph */
572 )
573 {
574  const int* const termToLcatreeNode = lcabuilder->termToLcatreeNode;
575  const int* const nodemapOrgToDist = sdgraph->nodemapOrgToDist;
576  const int* const lcatreeNodeToRmq = lcabuilder->lcatreeNodeToInds;
577  int* const rmq_nodeToRmqEntry = sdgraph->rmq_nodeToRmqEntry;
578  const int nnodes = graph_get_nNodes(g);
579
580  assert(termToLcatreeNode && nodemapOrgToDist && lcatreeNodeToRmq && rmq_nodeToRmqEntry);
581
582  for( int i = 0; i < nnodes; i++ )
583  {
584  assert(rmq_nodeToRmqEntry[i] == UNKNOWN);
585
586  if( Is_term(g->term[i]) )
587  {
588  const int node_distgraph = nodemapOrgToDist[i];
589  const int node_lcatree = termToLcatreeNode[node_distgraph];
590  const int node_rmq = lcatreeNodeToRmq[node_lcatree];
591
592  assert(graph_knot_isInRange(sdgraph->distgraph, node_distgraph));
593  rmq_nodeToRmqEntry[i] = node_rmq;
594
595  SCIPdebugMessage("node %d to %d \n", i, node_rmq);
596  }
597  }
598 }
599
600
601
602 /** initializes SD constant time query data */
603 static
605  SCIP* scip, /**< SCIP */
606  const GRAPH* g, /**< graph to initialize from */
607  SDGRAPH* sdgraph /**< the SD graph */
608 )
609 {
610  const int nnodes = graph_get_nNodes(g);
611  const int rmqlength = sdqueryGetRmqLength(sdgraph);
612  const int logrmqlength = (int) log2(rmqlength);
613  LCABUILDER* lcabuilder;
614 #ifndef NDEBUG
615  const int nterms = graph_get_nNodes(sdgraph->distgraph);
616
617  assert(nterms >= 2);
618  assert(nterms == g->terms);
619  assert(logrmqlength >= 0);
620  assert(!sdgraph->rmq_edgecosts);
621  assert(!sdgraph->rmq_sparseTable);
622  assert(sdgraph->usingRMQ);
623 #endif
624
625  sdgraph->rmq_loglength = logrmqlength;
626
627  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->rmq_nodeToRmqEntry), nnodes) );
628  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->rmq_sparseTable), (logrmqlength + 1) * rmqlength) );
629  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->rmq_edgecosts), rmqlength) );
630
631  for( int i = 0; i < nnodes; i++ )
632  sdgraph->rmq_nodeToRmqEntry[i] = UNKNOWN;
633
634  SCIP_CALL( sdqueryLcaBuilderInit(scip, sdgraph, &lcabuilder) );
635  SCIP_CALL( sdqueryBuildBinaryTree(scip, sdgraph, lcabuilder) );
636  sdqueryBuildRmq(sdgraph, lcabuilder);
637  sdqueryBuildNodesToRmqMap(g, lcabuilder, sdgraph);
638
639  sdqueryLcaBuilderFree(scip, &lcabuilder);
640
641  return SCIP_OKAY;
642 }
643
644
645
646
647
648 /** builds map */
649 static
651  const GRAPH* g, /**< graph */
652  const LCABUILDER* lcabuilder, /**< the builder */
653  SDGRAPH* sdgraph /**< the SD graph */
654 )
655 {
656  const int* const termToLcatreeNode = lcabuilder->termToLcatreeNode;
657  const int* const nodemapOrgToDist = sdgraph->nodemapOrgToDist;
658  int* const fullq_nodeToIdx = sdgraph->fullq_nodeToIdx;
659  const int nnodes = graph_get_nNodes(g);
660
661  assert(termToLcatreeNode && nodemapOrgToDist && fullq_nodeToIdx);
662
663  for( int i = 0; i < nnodes; i++ )
664  {
665  if( Is_term(g->term[i]) )
666  {
667  const int node_distgraph = nodemapOrgToDist[i];
668  const int node_lcatree = termToLcatreeNode[node_distgraph];
669
670  assert(graph_knot_isInRange(sdgraph->distgraph, node_distgraph));
671  fullq_nodeToIdx[i] = node_lcatree;
672
673  SCIPdebugMessage("FullMap: node %d to %d \n", i, node_lcatree);
674  }
675  else
676  {
677  fullq_nodeToIdx[i] = UNKNOWN;
678  }
679  }
680 }
681
682
683 /** Gets special distance (i.e. bottleneck distance) from graph.
684  * Corresponds to bottleneck length of path between term1 and term2 on distance graph */
685 static inline
687  int term1, /**< terminal 1 */
688  int term2, /**< terminal 2 */
689  const SDGRAPH* sdgraph /**< the SD graph */
690 )
691 {
692  const int* const rmq_nodeToRmqEntry = sdgraph->rmq_nodeToRmqEntry;
693  SCIP_Real sd;
694  int rmqindex_min;
695  int rmqindex_max;
696
697  assert(term1 != term2);
698
699  SCIPdebugMessage("query %d--%d \n", term1, term2);
700
701  if( rmq_nodeToRmqEntry[term1] <= rmq_nodeToRmqEntry[term2] )
702  {
703  rmqindex_min = rmq_nodeToRmqEntry[term1];
704  rmqindex_max = rmq_nodeToRmqEntry[term2];
705  }
706  else
707  {
708  rmqindex_min = rmq_nodeToRmqEntry[term2];
709  rmqindex_max = rmq_nodeToRmqEntry[term1];
710  }
711
712  assert(rmqindex_min <= rmqindex_max);
713  assert(0 <= rmqindex_min);
714  assert(rmqindex_max < sdqueryGetRmqLength(sdgraph));
715
716  if( rmqindex_min == rmqindex_max )
717  {
718  sd = sdgraph->rmq_edgecosts[rmqindex_min];
719  SCIPdebugMessage("...mapped to same RMQ index \n");
720  }
721  else
722  {
723  const SCIP_Real* const rmq_sparsetable = sdgraph->rmq_sparseTable;
724  const int rowlength = sdgraph->rmq_loglength + 1;
725  const unsigned int msbit = (unsigned int) log2floor(rmqindex_max - rmqindex_min);
726  const int pos_lower = rmqindex_min * rowlength + msbit;
727  const int pos_upper = (rmqindex_max - (int) (1 << msbit)) * rowlength + msbit;
728
729  SCIPdebugMessage("MSB=%u \n", msbit);
730  SCIPdebugMessage("RMQ indices: %d, %d \n", rmqindex_min, rmqindex_max);
731  SCIPdebugMessage("start indices: %d, %d \n", rmqindex_min, (rmqindex_max - (int) (1 << msbit)));
732  SCIPdebugMessage("max{%f, %f} \n", rmq_sparsetable[pos_lower], rmq_sparsetable[pos_upper]);
733
734  sd = MAX(rmq_sparsetable[pos_lower], rmq_sparsetable[pos_upper]);
735  }
736
737  assert(LT(sd, FARAWAY));
738  //printf("%f vs %f \n", sd, sdgraphGetSd(term1, term2, (SDGRAPH*) sdgraph));
739
740  assert(EQ(sd, sdgraphGetSd(term1, term2, (SDGRAPH*) sdgraph)));
741
742  return sd;
743 }
744
745
746 /** Gets special distance (i.e. bottleneck distance) from graph.
747  * Corresponds to bottleneck length of path between term1 and term2 on distance graph */
748 static inline
750  int term1, /**< terminal 1 */
751  int term2, /**< terminal 2 */
752  const SDGRAPH* sdgraph /**< the SD graph */
753 )
754 {
755  int pos;
756  const int index1 = sdgraph->fullq_nodeToIdx[term1];
757  const int index2 = sdgraph->fullq_nodeToIdx[term2];
758
759  SCIPdebugMessage("query %d--%d \n", term1, term2);
760  assert(term1 != term2);
761  assert(sdgraph->fullq_dimension >= 1);
762  assert(0 <= index1 && index1 < sdgraph->fullq_dimension);
763  assert(0 <= index2 && index2 < sdgraph->fullq_dimension);
764
765 #ifdef SDQUERYFULL_HALFMATRIX
766  if( index1 <= index2 )
767  {
768  pos = sdgraph->fullq_IdxToStart[index2] + index1;
769  //pos = (index2 * (index2 + 1)) / 2 + index1;
770  }
771  else
772  {
773  pos = sdgraph->fullq_IdxToStart[index1] + index2;
774  // pos = (index1 * (index1 + 1)) / 2 + index2;
775  }
776 #else
777  pos = index1 * sdgraph->fullq_dimension + index2;
778 #endif
779
780  assert(pos < sdgraph->fullq_size);
781  return sdgraph->fullq_dists[pos];
782 }
783
784
785 /** initializes full SD query data structure by DFS */
786 static
788  int root, /**< root node (LCA tree position) */
789  const BINARYNODE* lcatree, /**< tree */
790  int nlcanodes, /**< number of nodes */
791  const SCIP_Real* lcatree_costs, /**< LCA tree costs */
792  SCIP_Bool* RESTRICT nodes_isVisited, /**< per node: is visited? */
793  UF* uf, /**< union-find data structure */
794  SCIP_Real* RESTRICT fullq_dists /**< distances to be computed */
795 )
796 {
797  const BINARYNODE bnode = lcatree[root];
798 #ifdef SDQUERYFULL_HALFMATRIX
799  const int offset_root = (root * (root + 1)) / 2;
800 #endif
801
802  assert(!nodes_isVisited[root]);
803  SCIPdebugMessage("checking node %d (cost=%f)\n", root, lcatree_costs[root]);
804
805  nodes_isVisited[root] = TRUE;
806
807  if( bnode.child1 != UNKNOWN )
808  {
809  if( !nodes_isVisited[bnode.child1] )
810  {
811  sdqueryFullDfs(bnode.child1, lcatree, nlcanodes, lcatree_costs, nodes_isVisited, uf, fullq_dists);
812  /* NOTE: root is also made the identifier of the set */
813  SCIPStpunionfindUnion(uf, root, bnode.child1, FALSE);
814  }
815  }
816
817  if( bnode.child2 != UNKNOWN )
818  {
819  if( !nodes_isVisited[bnode.child2] )
820  {
821  sdqueryFullDfs(bnode.child2, lcatree, nlcanodes, lcatree_costs, nodes_isVisited, uf, fullq_dists);
822  SCIPStpunionfindUnion(uf, root, bnode.child2, FALSE);
823  }
824  }
825
826  /* 2. fill distances, also from root to root! */
827
828 #ifdef SDQUERYFULL_HALFMATRIX
829  for( int i = 0; i <= root; i++ )
830  {
831  if( nodes_isVisited[i] )
832  {
833  const int ancestor = SCIPStpunionfindFind(uf, i);
834  const int pos = offset_root + i;
835
836  assert(0 <= ancestor && ancestor < nlcanodes);
837
838  SCIPdebugMessage("add %d->%d=%f \n", root, i, lcatree_costs[ancestor]);
839  fullq_dists[pos] = lcatree_costs[ancestor];
840  }
841  }
842
843  for( int i = root + 1; i < nlcanodes; i++ )
844  {
845  if( nodes_isVisited[i] )
846  {
847  const int ancestor = SCIPStpunionfindFind(uf, i);
848  int pos = (i * (i + 1)) / 2 + root;
849
850  assert(0 <= ancestor && ancestor < nlcanodes);
851
852  SCIPdebugMessage("add %d->%d=%f \n", i, root, lcatree_costs[ancestor]);
853  fullq_dists[pos] = lcatree_costs[ancestor];
854  }
855  }
856 #else
857  for( int i = 0; i < nlcanodes; i++ )
858  {
859  if( nodes_isVisited[i] )
860  {
861  const int ancestor = SCIPStpunionfindFind(uf, i);
862
863  assert(0 <= ancestor && ancestor < nlcanodes);
864
865  SCIPdebugMessage("double add %d->%d=%f \n", i, root, lcatree_costs[ancestor]);
866  fullq_dists[root * nlcanodes + i] = lcatree_costs[ancestor];
867  fullq_dists[i * nlcanodes + root] = lcatree_costs[ancestor];
868  }
869  }
870 #endif
871 }
872
873
874 /** builds full representation */
875 static
877  SCIP* scip, /**< SCIP */
878  SDGRAPH* sdgraph, /**< the SD graph */
879  LCABUILDER* lcabuilder /**< the builder */
880 )
881 {
882  UF uf;
883  const BINARYNODE* const lcatree = lcabuilder->lcatree;
884  const SCIP_Real* const lcatree_costs = lcabuilder->lcatree_costs;
885  SCIP_Real* const fullq_dists = sdgraph->fullq_dists;
886  SCIP_Bool* nodes_isVisited;
887  const int nlcanodes = sdgraph->distgraph->knots - 1;
888
889  assert(fullq_dists);
890  assert(nlcanodes >= 1);
891
892  SCIP_CALL( SCIPallocBufferArray(scip, &(nodes_isVisited), nlcanodes) );
893
894  SCIP_CALL( SCIPStpunionfindInit(scip, &uf, nlcanodes) );
895
896  for( int i = 0; i < nlcanodes; i++ )
897  nodes_isVisited[i] = FALSE;
898
899 #ifdef SDQUERYFULL_HALFMATRIX
900  for( int i = 0; i < nlcanodes; i++ )
901  {
902  const int pos = (i * (i + 1)) / 2;
903  sdgraph->fullq_IdxToStart[i] = pos;
904  }
905 #endif
906
907  /* NOTE: because we have sorted in ascending edge cost order, the root is 0 */
908  sdqueryFullDfs(0, lcatree, nlcanodes, lcatree_costs, nodes_isVisited, &uf, fullq_dists);
909
910  SCIPStpunionfindFreeMembers(scip, &uf);
911  SCIPfreeBufferArray(scip, &nodes_isVisited);
912
913  return SCIP_OKAY;
914 }
915
916
917 /** initializes full SD constant time query data */
918 static
920  SCIP* scip, /**< SCIP */
921  const GRAPH* g, /**< graph to initialize from */
922  SDGRAPH* sdgraph /**< the SD graph */
923 )
924 {
925  const int nterms = graph_get_nNodes(sdgraph->distgraph);
926  const int fullq_dimension = (nterms - 1);
927 #ifdef SDQUERYFULL_HALFMATRIX
928  const int fullq_size = (fullq_dimension * fullq_dimension) / 2 + fullq_dimension;
929 #else
930  const int fullq_size = (fullq_dimension * fullq_dimension);
931 #endif
932  LCABUILDER* lcabuilder;
933
934  assert(!sdgraph->fullq_dists);
935  assert(!sdgraph->usingRMQ);
936  assert(nterms == g->terms);
937  assert(g->terms >= 2 && g->terms <= 1000);
938  SCIPdebugMessage("fullq_dimension=%d, fullq_size=%d \n", fullq_dimension, fullq_size);
939
940  sdgraph->fullq_dimension = fullq_dimension;
941  sdgraph->fullq_size = fullq_size;
942
943  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->fullq_dists), fullq_size) );
944  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->fullq_nodeToIdx), g->knots) );
945 #ifdef SDQUERYFULL_HALFMATRIX
946  SCIP_CALL( SCIPallocMemoryArray(scip, &(sdgraph->fullq_IdxToStart), g->knots) );
947 #endif
948
949  SCIP_CALL( sdqueryLcaBuilderInit(scip, sdgraph, &lcabuilder) );
950  SCIP_CALL( sdqueryBuildBinaryTree(scip, sdgraph, lcabuilder) );
951  SCIP_CALL( sdqueryFullBuild(scip, sdgraph, lcabuilder) );
952  sdqueryBuildNodesToFullMap(g, lcabuilder, sdgraph);
953
954  sdqueryLcaBuilderFree(scip, &lcabuilder);
955
956  return SCIP_OKAY;
957 }
958
959
960 /** initializes SD constant time query data */
961 static
963  SCIP* scip, /**< SCIP */
964  const GRAPH* g, /**< graph to initialize from */
965  SDGRAPH* sdgraph /**< the SD graph */
966 )
967 {
969  {
970  sdgraph->usingRMQ = TRUE;
971  SCIP_CALL( sdqueryRmqInit(scip, g, sdgraph) );
972  }
973  else
974  {
975  sdgraph->usingRMQ = FALSE;
976  SCIP_CALL( sdqueryFullInit(scip, g, sdgraph) );
977  }
978
979  return SCIP_OKAY;
980 }
981
982
983 /** frees SD constant time query data */
984 static
986  SCIP* scip, /**< SCIP */
987  SDGRAPH* sdgraph /**< the SD graph */
988 )
989 {
990  if( sdgraph->usingRMQ )
991  {
992  assert(sdgraph->rmq_edgecosts && sdgraph->rmq_sparseTable && sdgraph->rmq_nodeToRmqEntry);
993  assert(!sdgraph->fullq_dists && !sdgraph->fullq_nodeToIdx);
994
995  SCIPfreeMemoryArray(scip, &(sdgraph->rmq_edgecosts));
996  SCIPfreeMemoryArray(scip, &(sdgraph->rmq_sparseTable));
997  SCIPfreeMemoryArray(scip, &(sdgraph->rmq_nodeToRmqEntry));
998  }
999  else
1000  {
1001  assert(sdgraph->fullq_dists);
1002  assert(!sdgraph->rmq_edgecosts && !sdgraph->rmq_sparseTable && !sdgraph->rmq_nodeToRmqEntry);
1003 #ifdef SDQUERYFULL_HALFMATRIX
1004  SCIPfreeMemoryArray(scip, &(sdgraph->fullq_IdxToStart));
1005 #endif
1006  SCIPfreeMemoryArray(scip, &(sdgraph->fullq_nodeToIdx));
1007  SCIPfreeMemoryArray(scip, &(sdgraph->fullq_dists));
1008  }
1009 }
1010
1011
1012 /** gets number of edges */
1013 static
1015  const GRAPH* g /**< graph to initialize from */
1016  )
1017 {
1018  int maxnedges;
1019  const int nedges = graph_get_nEdges(g);
1020  const int nterms = g->terms;
1021  const SCIP_Longint terms2 = (SCIP_Longint) (nterms - 1) * nterms;
1022
1023  if( nedges >= terms2 )
1024  {
1025  assert(terms2 <= INT_MAX);
1026  maxnedges = terms2;
1027  }
1028  else
1029  {
1030  maxnedges = nedges;
1031  }
1032
1033  return maxnedges;
1034 }
1035
1036
1037 /** adds nodes to distance graph */
1038 static
1040  const GRAPH* g, /**< graph to initialize from */
1041  int* RESTRICT distnodes_id, /**< IDs of nodes */
1042  GRAPH* distgraph /**< distance graph */
1043 )
1044 {
1045  const int nnodes = graph_get_nNodes(g);
1046  int nnodes_distgraph = 0;
1047
1048  /* add the nodes */
1049  for( int k = 0; k < nnodes; k++ )
1050  {
1051  if( Is_term(g->term[k]) )
1052  {
1054  distnodes_id[k] = nnodes_distgraph++;
1055  }
1056  else
1057  {
1058  distnodes_id[k] = UNKNOWN;
1059  }
1060  }
1061
1062  assert(distgraph->knots == nnodes_distgraph);
1063  assert(distgraph->knots == g->terms);
1064
1065  graph_knot_chg(distgraph, 0, STP_TERM);
1066  distgraph->source = 0;
1067 }
1068
1069
1070 /** gets SD for boundary edge */
1071 static inline
1073  int tail, /**< tail */
1075  int vbase_tail, /**< base */
1076  int vbase_head, /**< base */
1077  SCIP_Real edgecost, /**< cost */
1078  const SCIP_Real* nodes_vdist, /**< distance */
1079  const SDPROFIT* sdprofit /**< profit */
1080  )
1081 {
1082  const SCIP_Real profit_tail = reduce_sdprofitGetProfit(sdprofit, tail, vbase_head, vbase_tail);
1084  SCIP_Real distSimple = edgecost + nodes_vdist[tail] + nodes_vdist[head];
1085  const SCIP_Real distAll = distSimple - profit_tail - profit_head;
1086  const SCIP_Real distTailHead = edgecost + nodes_vdist[tail] - profit_tail;
1088
1089  const SCIP_Real dist = miscstp_maxReal((SCIP_Real[])
1092  6);
1093  return dist;
1094 }
1095
1096
1097 /** gets SD for boundary edge */
1098 static inline
1100  int tail, /**< tail */
1102  int vbase_tail, /**< base */
1103  int vbase_head, /**< base */
1104  SCIP_Real edgecost, /**< cost */
1105  SCIP_Real dist_tail, /**< distance */
1106  SCIP_Real dist_head, /**< distance */
1107  const SDPROFIT* sdprofit /**< profit */
1108  )
1109 {
1110  const SCIP_Real profit_tail = reduce_sdprofitGetProfit(sdprofit, tail, vbase_head, vbase_tail);
1112  SCIP_Real distSimple = edgecost + dist_tail + dist_head;
1113  const SCIP_Real distAll = distSimple - profit_tail - profit_head;
1114  const SCIP_Real distTailHead = edgecost + dist_tail - profit_tail;
1116
1117  const SCIP_Real dist = miscstp_maxReal((SCIP_Real[])
1120  6);
1121  return dist;
1122 }
1123
1124
1125 /** gets SD for boundary edge ... choose among nearest terminals (w.r.t. implied SD) */
1126 static inline
1128  const GRAPH* g, /**< graph */
1129  const TPATHS* tpaths, /**< terminal paths */
1130  int tail, /**< tail */
1132  SCIP_Real edgecost, /**< cost */
1133  const SDPROFIT* sdprofit, /**< profit */
1134  int* base_tail, /**< base of tail */
1136  )
1137 {
1138  SCIP_Real dists_tail[3];
1140  int terms_tail[3];
1142  int nterms_tail;
1144  SCIP_Real dist = FARAWAY;
1145
1146  assert(g && tpaths && sdprofit);
1147
1148  *base_tail = -1;
1150
1151  graph_tpathsGet3CloseTerms(g, tpaths, tail, FARAWAY, terms_tail, NULL, dists_tail, &nterms_tail);
1153
1154  for( int i = 0; i < nterms_tail; ++i )
1155  {
1156  for( int j = 0; j < nterms_head; ++j )
1157  {
1158  if( terms_tail[i] != terms_head[j] )
1159  {
1160  const SCIP_Real distnew =
1162
1163  if( distnew < dist )
1164  {
1165  *base_tail = terms_tail[i];
1167  dist = distnew;
1168  }
1169  }
1170  }
1171  }
1172
1173  return dist;
1174 }
1175
1176
1177 /** inserts new edge */
1178 static inline
1180  SCIP* scip, /**< SCIP data structure */
1181  int sdnode1, /**< end node 1 */
1182  int sdnode2, /**< end node 2 */
1183  SCIP_Real edgecost, /**< cost */
1184  int edgeid, /**< ID or -1 */
1185  int* RESTRICT edgeorg, /**< IDs of edges or NULL */
1186  GRAPH* distgraph, /**< the SD graph */
1187  SCIP_Bool* success /**< could the edge be added? */
1188 )
1189 {
1190  int ne;
1191  *success = TRUE;
1192
1193 #ifndef NDEBUG
1194  assert(sdnode1 != sdnode2);
1195  assert(0 <= sdnode1 && sdnode1 < distgraph->knots);
1196  assert(0 <= sdnode2 && sdnode2 < distgraph->knots);
1197  assert(GE(edgecost, 0.0));
1198  assert(edgeorg != NULL || edgeid == -1);
1199 #endif
1200
1201  /* find the corresponding edge in the distance graph */
1202  for( ne = distgraph->outbeg[sdnode1]; ne != EAT_LAST; ne = distgraph->oeat[ne] )
1203  {
1204  if( distgraph->head[ne] == sdnode2 )
1205  break;
1206  }
1207
1208  /* edge exists already? */
1209  if( ne != EAT_LAST )
1210  {
1211  assert(ne >= 0);
1213  assert(distgraph->tail[ne] == sdnode1);
1214
1215  if( distgraph->cost[ne] > edgecost )
1216  {
1217  distgraph->cost[ne] = edgecost;
1218  distgraph->cost[Edge_anti(ne)] = edgecost;
1219
1220  if( edgeorg != NULL )
1221  edgeorg[ne / 2] = edgeid;
1222  }
1223  }
1224  else
1225  {
1226  assert(distgraph->edges <= distgraph->esize);
1227
1228  if( distgraph->edges == distgraph->esize )
1229  {
1230  *success = FALSE;
1231  }
1232  else
1233  {
1234  if( edgeorg != NULL )
1235  edgeorg[distgraph->edges / 2] = edgeid;
1236
1237  graph_edge_add(scip, distgraph, sdnode1, sdnode2, edgecost, edgecost);
1238  }
1239  }
1240 }
1241
1242
1243
1244 /** adds edges to distance graph, given a Voronoi diagram */
1245 static
1247  SCIP* scip, /**< SCIP data structure */
1248  const GRAPH* g, /**< graph to initialize from */
1249  const int* distnodes_id, /**< IDs of nodes */
1250  const VNOI* vnoi, /**< Voronoi */
1251  const SDPROFIT* sdprofit, /**< profit or NULL */
1252  int* RESTRICT edgeorg, /**< IDs of edges */
1253  GRAPH* distgraph /**< distance graph */
1254 )
1255 {
1256  const int nnodes = graph_get_nNodes(g);
1257  const int nedges = graph_get_nEdges(g);
1258  const SCIP_Real* nodes_vdist = vnoi->nodes_dist;
1259  const int* RESTRICT nodes_vbase = vnoi->nodes_base;
1260  const SCIP_Bool useProfit = (sdprofit != NULL);
1261
1262  for( int e = 0; e < nedges / 2; e++ )
1263  edgeorg[e] = UNKNOWN;
1264
1265  /* add the edges */
1266  for( int tail = 0; tail < nnodes; tail++ )
1267  {
1268  const int vbase_tail = nodes_vbase[tail];
1269
1270  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
1271  {
1274
1275  assert(tail == g->tail[e]);
1276
1277  if( vbase_tail != vbase_head )
1278  {
1279  SCIP_Bool success;
1280  const SCIP_Real distance = useProfit ?
1282  : (g->cost[e] + nodes_vdist[tail] + nodes_vdist[head]);
1283
1284  // if( LT(distance, g->cost[e] + nodes_vdist[tail] + nodes_vdist[head]) )
1285  // printf("distance: %f < %f \n", distance, g->cost[e] + nodes_vdist[tail] + nodes_vdist[head]);
1286
1287  assert(LE(distance, g->cost[e] + nodes_vdist[tail] + nodes_vdist[head]));
1288  assert(Is_term(g->term[vbase_tail]));
1290
1291  distgraphInsertEdge(scip, distnodes_id[vbase_tail], distnodes_id[vbase_head], distance, e,
1292  edgeorg, distgraph, &success);
1293
1294  assert(success);
1295  }
1296  }
1297  }
1298 }
1299
1300
1301
1302 /** helper */
1303 static
1305  const GRAPH* g, /**< graph to initialize from */
1306  SDGRAPH* g_sd /**< the SD graph */
1307 )
1308 {
1309  const int nnodes = graph_get_nNodes(g);
1310  SCIP_Real* mstsdist = g_sd->mstsdist;
1311  const int nedges = graph_get_nEdges(g);
1312
1313  for( int i = 0; i < nnodes; i++ )
1314  {
1315  mstsdist[i] = -1.0;
1316  }
1317
1318  g_sd->nnodesorg = nnodes;
1319  g_sd->nedgesorg = nedges;
1322  g_sd->usingRMQ = TRUE;
1323  g_sd->fullq_dists = NULL;
1324  g_sd->fullq_nodeToIdx = NULL;
1325  g_sd->fullq_size = -1;
1326  g_sd->fullq_dimension = -1;
1327  g_sd->rmq_edgecosts = NULL;
1328  g_sd->rmq_sparseTable = NULL;
1329  g_sd->rmq_nodeToRmqEntry = NULL;
1330  g_sd->rmq_loglength = -1;
1331 }
1332
1333
1334 /** allocates memory */
1335 static
1337  SCIP* scip, /**< SCIP */
1338  const GRAPH* g, /**< graph to initialize from */
1339  SDGRAPH** sdgraph /**< the SD graph */
1340 )
1341 {
1342  const int nnodes = graph_get_nNodes(g);
1343  const int nedges = graph_get_nEdges(g);
1344  const int nterms = g->terms;
1345  SDGRAPH* g_sd;
1346
1347  SCIP_CALL( SCIPallocMemory(scip, sdgraph) );
1348  g_sd = *sdgraph;
1349
1350  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->sdmst), nterms) );
1351  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->mstcosts), nterms) );
1352  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->mstsdist), nnodes) );
1353  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->nodemapOrgToDist), nnodes) );
1354  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->halfedge_isInMst), nedges / 2) );
1355
1356  sdgraphSetDefaults(g, g_sd);
1357
1358  return SCIP_OKAY;
1359 }
1360
1361
1362
1363 /** allocates memory */
1364 static
1366  SCIP* scip, /**< SCIP */
1367  const GRAPH* g, /**< graph to initialize from */
1368  SDGRAPH** sdgraph /**< the SD graph */
1369 )
1370 {
1371  const int nnodes = graph_get_nNodes(g);
1372  const int nterms = g->terms;
1373  SDGRAPH* g_sd;
1374
1375  SCIP_CALL( SCIPallocMemory(scip, sdgraph) );
1376  g_sd = *sdgraph;
1377
1378  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->sdmst), nterms) );
1379  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->mstcosts), nterms) );
1380  SCIP_CALL( SCIPallocMemoryArray(scip, &(g_sd->mstsdist), nnodes) );
1381
1382  g_sd->nodemapOrgToDist = NULL;
1383  g_sd->halfedge_isInMst = NULL;
1384
1385  sdgraphSetDefaults(g, g_sd);
1386
1387  return SCIP_OKAY;
1388 }
1389
1390 /** adds edges to distance graph, given terminal paths */
1391 static
1393  SCIP* scip, /**< SCIP data structure */
1394  const GRAPH* g, /**< graph to initialize from */
1395  const int* distnodes_id, /**< IDs of nodes */
1396  const SDPROFIT* sdprofit, /**< profit or NULL */
1397  const TPATHS* tpaths, /**< terminal paths */
1398  GRAPH* distgraph /**< distance graph */
1399 )
1400 {
1401  const int nnodes = graph_get_nNodes(g);
1402  const SCIP_Bool useProfit = (sdprofit != NULL);
1403
1404  /* add the edges */
1405  for( int tail = 0; tail < nnodes; tail++ )
1406  {
1407  SCIP_Real dist_tail;
1408  int vbase_tail;
1409  graph_tpathsGetClosestTerm(g, tpaths, tail, &vbase_tail, NULL, &dist_tail);
1410
1411  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
1412  {
1416
1418
1419  if( vbase_tail != vbase_head )
1420  {
1421  SCIP_Bool success;
1422  const SCIP_Real distance = useProfit ?
1424  : (g->cost[e] + dist_tail + dist_head);
1425
1426  assert(LE(distance, g->cost[e] + dist_tail + dist_head));
1428
1429  distgraphInsertEdge(scip, distnodes_id[vbase_tail], distnodes_id[vbase_head], distance, -1, NULL, distgraph, &success);
1430  assert(success);
1431  }
1432  }
1433  }
1434 }
1435
1436
1437 /** builds distance graph */
1438 static
1440  SCIP* scip, /**< SCIP */
1441  const GRAPH* g, /**< graph to initialize from */
1442  const SDPROFIT* sdprofit, /**< profit or NULL */
1443  SDGRAPH* g_sd, /**< the SD graph */
1444  VNOI** vnoi, /**< Voronoi */
1445  int** distedge2org /**< array of size nedges / 2 */
1446 )
1447 {
1448  GRAPH* distgraph;
1449  int* RESTRICT distnodes_id = g_sd->nodemapOrgToDist;
1450  const int nedges = graph_get_nEdges(g);
1451  const int nedges_distgraph = distgraphGetNedges(g);
1452
1453  assert(g->terms >= 1);
1454
1455  SCIP_CALL( SCIPallocBufferArray(scip, distedge2org, nedges / 2) );
1456
1457  /* build biased Voronoi diagram */
1458  SCIP_CALL( graph_vnoiInit(scip, g, TRUE, vnoi) );
1459
1460  if( sdprofit )
1461  graph_vnoiComputeImplied(scip, g, sdprofit, *vnoi);
1462  else
1463  graph_vnoiCompute(scip, g, *vnoi);
1464
1465  /* build distance graph from Voronoi diagram */
1466  SCIP_CALL( graph_init(scip, &(g_sd->distgraph), g->terms, nedges_distgraph, 1) );
1467  distgraph = g_sd->distgraph;
1469  distgraphAddEdges(scip, g, distnodes_id, *vnoi, sdprofit, *distedge2org, distgraph);
1470
1471  assert(graph_valid(scip, distgraph));
1472
1473  return SCIP_OKAY;
1474 }
1475
1476
1477 /** builds distance graph */
1478 static
1480  SCIP* scip, /**< SCIP */
1481  const GRAPH* g, /**< graph to initialize from */
1482  const SDPROFIT* sdprofit, /**< profit or NULL */
1483  const TPATHS* tpaths, /**< terminal paths */
1484  SDGRAPH* g_sd /**< the SD graph */
1485 )
1486 {
1487  GRAPH* distgraph;
1488  int* RESTRICT distnodes_id = g_sd->nodemapOrgToDist;
1489  const int nedges_distgraph = distgraphGetNedges(g);
1490
1491  assert(g->terms >= 1);
1492
1493  SCIP_CALL( graph_init(scip, &(g_sd->distgraph), g->terms, nedges_distgraph, 1) );
1494  distgraph = g_sd->distgraph;
1496  distgraphAddEdgesFromTpaths(scip, g, distnodes_id, sdprofit, tpaths, distgraph);
1497
1498  assert(graph_valid(scip, distgraph));
1499
1500  return SCIP_OKAY;
1501 }
1502
1503
1504 /** updates distance graph */
1505 static
1507  SCIP* scip, /**< SCIP */
1508  const GRAPH* g, /**< graph to initialize from */
1509  const SDPROFIT* sdprofit, /**< profit or NULL */
1510  const TPATHS* tpaths, /**< terminal paths */
1511  SDGRAPH* g_sd /**< the SD graph */
1512 )
1513 {
1514  int* hasharr;
1515  STP_Vectype(int) profitnodes_tail = NULL;
1517
1518  GRAPH* RESTRICT distgraph = g_sd->distgraph;
1519  const int* distnodes_id = g_sd->nodemapOrgToDist;
1520  const int nnodes = graph_get_nNodes(g);
1521  const int edgelimit = MIN(2 * distgraph->edges, distgraph->esize);
1522
1523  assert(LE(distgraph->edges, distgraph->esize));
1524
1525  StpVecReserve(scip, profitnodes_tail, nnodes);
1527  SCIP_CALL( SCIPallocCleanBufferArray(scip, &hasharr, nnodes) );
1528
1529  for( int tail = 0; tail < nnodes && distgraph->edges < edgelimit; tail++ )
1530  {
1531  for( int e = g->outbeg[tail]; e != EAT_LAST; e = g->oeat[e] )
1532  {
1534  int vbase_tail;
1536  SCIP_Real distance;
1537
1538  if( head < tail )
1539  continue;
1540
1542
1543  /* NOTE: we should not take the fast query method here, because sdgraphGetSd at least partly reflects
1544  * the change of the distance graph */
1545  if( LT(distance, FARAWAY) && LT(distance, sdgraphGetSd(vbase_tail, vbase_head, g_sd)) )
1546  {
1547  SCIP_Bool success = TRUE;
1548  assert(vbase_tail >= 0 && vbase_head >= 0);
1549
1550  graph_tpathsGetProfitNodes(scip, g, tpaths, sdprofit, tail, vbase_tail, profitnodes_tail);
1551  for( int k = 0; k < StpVecGetSize(profitnodes_tail); k++ )
1552  {
1553  assert(!hasharr[profitnodes_tail[k]]);
1554  hasharr[profitnodes_tail[k]] = 1;
1555  }
1557
1558  for( int k = 0; k < StpVecGetSize(profitnodes_head); k++ )
1559  {
1561  {
1562  success = FALSE;
1563  // printf("shared node: %d \n", profitnodes_head[k]);
1564  break;
1565  }
1566  }
1567
1568  for( int k = 0; k < StpVecGetSize(profitnodes_tail); k++ )
1569  {
1570  assert(1 == hasharr[profitnodes_tail[k]]);
1571  hasharr[profitnodes_tail[k]] = 0;
1572  }
1573
1574  if( !success )
1575  {
1576  // printf("fail for %d %d \n", vbase_tail, vbase_head);
1577  continue;
1578  }
1580
1581  distgraphInsertEdge(scip, distnodes_id[vbase_tail], distnodes_id[vbase_head], distance, -1, NULL, distgraph, &success);
1582  assert(success);
1583
1584  if( distgraph->edges >= edgelimit )
1585  break;
1586  }
1587  }
1588  }
1589 #ifndef NDEBUG
1590  for( int i = 0; i < nnodes; i++ )
1591  assert(hasharr[i] == 0);
1592 #endif
1593
1594  SCIPfreeCleanBufferArray(scip, &hasharr);
1596  StpVecFree(scip, profitnodes_tail);
1597
1598  return SCIP_OKAY;
1599 }
1600
1601
1602 /** builds MST costs (ordered) for distance graph */
1603 static
1605  SDGRAPH* g_sd /**< the SD graph */
1606 )
1607 {
1608  const GRAPH* const distgraph = g_sd->distgraph;
1609  const PATH* const mst = g_sd->sdmst;
1610  SCIP_Real* RESTRICT mstcosts = g_sd->mstcosts;
1611  const int nnodes_distgraph = graph_get_nNodes(distgraph);
1612
1613  assert(mst[0].edge == UNKNOWN);
1614  assert(nnodes_distgraph >= 1);
1615
1616  for( int k = 1; k < nnodes_distgraph; k++ )
1617  {
1618 #ifndef NDEBUG
1619  const int e = mst[k].edge;
1620  assert(e >= 0);
1621  assert(GE(distgraph->cost[e], 0.0));
1622  assert(EQ(mst[k].dist, distgraph->cost[e]));
1623 #endif
1624
1625  mstcosts[k - 1] = mst[k].dist;
1626  }
1627
1628  SCIPsortDownReal(mstcosts, nnodes_distgraph - 1);
1629
1630  /* debug sentinel */
1631  if( nnodes_distgraph > 1 )
1632  {
1633  mstcosts[nnodes_distgraph - 1] = -FARAWAY;
1634  }
1635  else
1636  {
1637  mstcosts[0] = 0.0;
1638  }
1639 }
1640
1641
1642 /** marks original edges corresponding to MST */
1643 static
1645  const GRAPH* g, /**< graph to initialize from */
1646  const VNOI* vnoi, /**< Voronoi */
1647  const int* distedge2org, /**< array of size nedges / 2 */
1648  SDGRAPH* g_sd /**< the SD graph */
1649 )
1650 {
1651  GRAPH* distgraph = g_sd->distgraph;
1652  PATH* RESTRICT mst = g_sd->sdmst;
1653  const int nnodes_distgraph = graph_get_nNodes(distgraph);
1654  const int* const nodes_vbase = vnoi->nodes_base;
1655  const int* const nodes_vpred = vnoi->nodes_predEdge;
1656  STP_Bool* RESTRICT orgedges_isInMst = g_sd->halfedge_isInMst;
1657
1658 #ifndef NDEBUG
1659  const int nedges = graph_get_nEdges(g);
1660
1661  for( int e = 0; e < nedges / 2; e++ )
1662  assert(!orgedges_isInMst[e]);
1663 #endif
1664
1665  for( int k = 1; k < nnodes_distgraph; k++ )
1666  {
1667  const int e = mst[k].edge;
1668  const int ne = distedge2org[e / 2];
1669
1670  assert(e >= 0);
1671  assert(e < nedges);
1672
1673  orgedges_isInMst[ne / 2] = TRUE;
1674
1675  for( int v = g->head[ne]; v != nodes_vbase[v]; v = g->tail[nodes_vpred[v]] )
1676  {
1677  orgedges_isInMst[nodes_vpred[v] / 2] = TRUE;
1678  }
1679
1680  for( int v = g->tail[ne]; v != nodes_vbase[v]; v = g->tail[nodes_vpred[v]] )
1681  {
1682  orgedges_isInMst[nodes_vpred[v]/ 2] = TRUE;
1683  }
1684
1685  assert(e != EAT_LAST);
1686  }
1687 }
1688
1689 /** builds MST on distance graph */
1690 static
1692  SCIP* scip, /**< SCIP */
1693  const GRAPH* g, /**< graph to initialize from */
1694  SDGRAPH* g_sd /**< the SD graph */
1695 )
1696 {
1697  PATH* RESTRICT mst = g_sd->sdmst;
1698  GRAPH* distgraph = g_sd->distgraph;
1699  SCIP_Real maxcost;
1700  const int nnodes_distgraph = graph_get_nNodes(distgraph);
1701  STP_Bool* RESTRICT orgedges_isInMst = g_sd->halfedge_isInMst;
1702  const int nedges = graph_get_nEdges(g);
1703  const SCIP_Bool distgraphIsInit = (distgraph->path_heap != NULL);
1704
1705  if( orgedges_isInMst)
1706  {
1707  for( int e = 0; e < nedges / 2; e++ )
1708  orgedges_isInMst[e] = FALSE;
1709  }
1710
1711  for( int k = 0; k < nnodes_distgraph; k++ )
1712  distgraph->mark[k] = TRUE;
1713
1714  if( !distgraphIsInit )
1715  {
1716  SCIP_CALL( graph_path_init(scip, distgraph) );
1717  }
1718  graph_path_exec(scip, distgraph, MST_MODE, distgraph->source, distgraph->cost, mst);
1719
1720  if( !distgraphIsInit )
1721  {
1722  graph_path_exit(scip, distgraph);
1723  }
1724
1725  assert(mst[0].edge == -1);
1726
1727  maxcost = 0.0;
1728  for( int k = 1; k < nnodes_distgraph; k++ )
1729  {
1730  const int e = mst[k].edge;
1731  const SCIP_Real cost = distgraph->cost[e];
1732
1733  assert(graph_edge_isInRange(distgraph, e));
1734
1735  if( cost > maxcost )
1736  maxcost = cost;
1737  }
1738
1739  g_sd->mstmaxcost = maxcost;
1740
1741  return SCIP_OKAY;
1742 }
1743
1744
1745 /** finalizes distance graph */
1746 static
1748  SCIP* scip, /**< SCIP data structure */
1749  VNOI** vnoi, /**< Voronoi data structure */
1750  int** edgeorg /**< array of size nedges / 2 */
1751 )
1752 {
1753  graph_vnoiFree(scip, vnoi);
1754  SCIPfreeBufferArray(scip, edgeorg);
1755 }
1756
1757
1758 /*
1759  * Interface methods
1760  */
1761
1762
1763 /** initializes SD graph */
1765  SCIP* scip, /**< SCIP */
1766  const GRAPH* g, /**< graph to initialize from */
1767  SDGRAPH** sdgraph /**< the SD graph */
1768 )
1769 {
1770  VNOI* vnoi;
1771  int* edgeorg;
1772  assert(scip && g && sdgraph);
1773
1774  SCIP_CALL( sdgraphAlloc(scip, g, sdgraph) );
1775  SCIP_CALL( sdgraphBuildDistgraph(scip, g, NULL, *sdgraph, &vnoi, &edgeorg) );
1776  SCIP_CALL( sdgraphMstBuild(scip, g, *sdgraph) );
1777  sdgraphMstMarkOrgEdges(g, vnoi, edgeorg, *sdgraph);
1778
1779  sdgraphFinalize(scip, &vnoi, &edgeorg);
1780
1781  SCIP_CALL( sdqueryInit(scip, g, *sdgraph) );
1782
1783  return SCIP_OKAY;
1784 }
1785
1786
1787 /** initializes SD graph */
1789  SCIP* scip, /**< SCIP */
1790  const GRAPH* g, /**< graph to initialize from */
1791  GRAPH* distgraph, /**< distance graph */
1792  int* node2dist, /**< map */
1793  SDGRAPH** sdgraph /**< the SD graph */
1794 )
1795 {
1796  assert(scip && g && sdgraph && distgraph && node2dist);
1797
1798  SCIP_CALL( sdgraphAllocRestricted(scip, g, sdgraph) );
1799  (*sdgraph)->distgraph = distgraph;
1800  (*sdgraph)->nodemapOrgToDist = node2dist;
1801
1802  SCIP_CALL( sdgraphMstBuild(scip, g, *sdgraph) );
1803
1804  SCIP_CALL( sdqueryInit(scip, g, *sdgraph) );
1805
1806  return SCIP_OKAY;
1807 }
1808
1809
1810 /** initializes biased SD graph */
1812  SCIP* scip, /**< SCIP */
1813  const GRAPH* g, /**< graph to initialize from */
1814  const SDPROFIT* sdprofit, /**< SD profit */
1815  SDGRAPH** sdgraph /**< the SD graph */
1816 )
1817 {
1818  VNOI* vnoi;
1819  int* edgeorg;
1820  assert(scip && g && sdgraph);
1821
1822  SCIP_CALL( sdgraphAlloc(scip, g, sdgraph) );
1823  SCIP_CALL( sdgraphBuildDistgraph(scip, g, sdprofit, *sdgraph, &vnoi, &edgeorg) );
1824  SCIP_CALL( sdgraphMstBuild(scip, g, *sdgraph) );
1825  sdgraphMstMarkOrgEdges(g, vnoi, edgeorg, *sdgraph);
1826
1827  sdgraphFinalize(scip, &vnoi, &edgeorg);
1828
1829  SCIP_CALL( sdqueryInit(scip, g, *sdgraph) );
1830
1831  return SCIP_OKAY;
1832 }
1833
1834
1835 /** initializes biased SD graph from given terminal paths */
1837  SCIP* scip, /**< SCIP */
1838  GRAPH* g, /**< graph to initialize from */
1839  const SDPROFIT* sdprofit, /**< SD profit */
1840  const TPATHS* tpaths, /**< terminal paths */
1841  SDGRAPH** sdgraph /**< the SD graph */
1842 )
1843 {
1844  assert(scip && g && sdgraph && tpaths);
1845
1846  SCIP_CALL( sdgraphAlloc(scip, g, sdgraph) );
1847  SCIP_CALL( sdgraphBuildDistgraphFromTpaths(scip, g, sdprofit, tpaths, *sdgraph) );
1848  SCIP_CALL( sdgraphMstBuild(scip, g, *sdgraph) );
1849
1850  if( sdprofit )
1851  {
1852  SCIP_CALL( sdgraphUpdateDistgraphFromTpaths(scip, g, sdprofit, tpaths, *sdgraph) );
1853  SCIP_CALL( sdgraphMstBuild(scip, g, *sdgraph) );
1854  }
1855
1856  /* NOTE probably we never need that...for extending reductions we anyway should only take biased paths */
1858  SCIPfreeMemoryArray(scip, &(*sdgraph)->halfedge_isInMst);
1859
1860  SCIP_CALL( sdqueryInit(scip, g, *sdgraph) );
1861
1862  return SCIP_OKAY;
1863 }
1864
1865
1866 /** return maximum MST edge cost */
1868  const SDGRAPH* sdgraph /**< the SD graph */
1869 )
1870 {
1871  assert(sdgraph);
1872  assert(GE(sdgraph->mstmaxcost, 0.0));
1873
1874  return sdgraph->mstmaxcost;
1875 }
1876
1877
1878 /** returns edge mark */
1880  const SDGRAPH* sdgraph /**< the SD graph */
1881 )
1882 {
1883  assert(sdgraph);
1884  assert(reduce_sdgraphHasMstHalfMark(sdgraph));
1885
1886  return sdgraph->halfedge_isInMst;
1887 }
1888
1889
1890 /** has edge mark? */
1892  const SDGRAPH* sdgraph /**< the SD graph */
1893 )
1894 {
1895  assert(sdgraph);
1897  {
1898  assert(sdgraph->halfedge_isInMst);
1899  return TRUE;
1900  }
1901
1902  assert(!sdgraph->halfedge_isInMst);
1903
1904  return FALSE;
1905 }
1906
1907
1908 /** is edge in current SD MST? */
1910  const SDGRAPH* sdgraph, /**< the SD graph */
1911  int edge /**< edge */
1912 )
1913 {
1914  assert(sdgraph);
1915  assert(reduce_sdgraphHasMstHalfMark(sdgraph));
1916  assert(sdgraph->halfedge_isInMst);
1917  assert(edge >= 0);
1918
1919  return sdgraph->halfedge_isInMst[edge / 2];
1920 }
1921
1922
1923 /** MST costs in descending order available? */
1925  const SDGRAPH* sdgraph /**< the SD graph */
1926 )
1927 {
1928  assert(sdgraph);
1929  assert(sdgraph->mstcosts);
1930
1932 }
1933
1934
1935
1936
1937 /** initializes all MST costs in descending order */
1939  SDGRAPH* sdgraph /**< the SD graph */
1940 )
1941 {
1942  assert(sdgraph);
1943  assert(sdgraph->mstcosts);
1944
1946  {
1947  sdgraphMstSortCosts(sdgraph);
1949  }
1950
1951  assert(reduce_sdgraphHasOrderedMstCosts(sdgraph));
1952  assert(sdgraph->mstcosts[0] == sdgraph->mstmaxcost);
1953 }
1954
1955
1956 /** Gets special distance (e.g. bottleneck distance) from distance graph.
1957  * Only works if both nodes are terminals! */
1959  int term1, /**< node 1 */
1960  int term2, /**< node 2 */
1961  SDGRAPH* sdgraph /**< the SD graph */
1962 )
1963 {
1964 #ifndef NDEBUG
1965  assert(sdgraph);
1966  assert(term1 != term2);
1967  assert(sdgraph->mstcosts);
1968  assert(sdgraph->mstsdist);
1969  assert(0 <= term1 && term1 < sdgraph->nnodesorg);
1970  assert(0 <= term2 && term2 < sdgraph->nnodesorg);
1971  assert(sdgraph->nodemapOrgToDist[term1] != UNKNOWN);
1972  assert(sdgraph->nodemapOrgToDist[term2] != UNKNOWN);
1973
1974  for( int i = 0; i < sdgraph->nnodesorg; i++ )
1975  assert(EQ(sdgraph->mstsdist[i], -1.0));
1976 #endif
1977
1978  if( sdgraph->usingRMQ )
1979  {
1980  return sdqueryGetSd(term1, term2, sdgraph);
1981  }
1982  else
1983  {
1984  return sdqueryFullGetSd(term1, term2, sdgraph);
1985  }
1986 }
1987
1988 /** returns all MST costs in descending order */
1990  const SDGRAPH* sdgraph /**< the SD graph */
1991 )
1992 {
1993  assert(sdgraph);
1994  assert(reduce_sdgraphHasOrderedMstCosts(sdgraph));
1995  assert(sdgraph->mstcosts[0] == sdgraph->mstmaxcost);
1996
1997  return sdgraph->mstcosts;
1998 }
1999
2000 /** returns mapping from original nodes to distance graph nodes */
2002  const SDGRAPH* sdgraph /**< the SD graph */
2003 )
2004 {
2005  assert(sdgraph);
2006
2007  return sdgraph->nodemapOrgToDist;
2008 }
2009
2010
2011
2012
2013 /** Inserts new edge.
2014  * NOTE: just a wrapper, should only be used by other reduce_sd* methods */
2016  SCIP* scip, /**< SCIP data structure */
2017  int sdnode1, /**< end node 1 */
2018  int sdnode2, /**< end node 2 */
2019  SCIP_Real edgecost, /**< cost */
2020  int edgeid, /**< ID or -1 */
2021  int* RESTRICT edgeorg, /**< IDs of edges or NULL */
2022  SDGRAPH* sdgraph, /**< the SD graph */
2023  SCIP_Bool* success /**< could the edge be added? */
2024 )
2025 {
2026  distgraphInsertEdge(scip, sdnode1, sdnode2, edgecost, edgeid, edgeorg, sdgraph->distgraph, success);
2027 }
2028
2029
2030 /** Builds MST on distance graph.
2031  * NOTE: just a wrapper, should only be used by other reduce_sd* methods */
2033  SCIP* scip, /**< SCIP */
2034  const GRAPH* g, /**< graph to initialize from */
2035  SDGRAPH* g_sd /**< the SD graph */
2036 )
2037 {
2038  SCIP_CALL( sdgraphMstBuild(scip, g, g_sd) );
2039
2040  return SCIP_OKAY;
2041 }
2042
2043
2044 /** Builds MST costs (ordered) for distance graph
2045 * NOTE: just a wrapper, should only be used by other reduce_sd* methods */
2047  SDGRAPH* g_sd /**< the SD graph */
2048 )
2049 {
2050  sdgraphMstSortCosts(g_sd);
2051 }
2052
2053
2054 /** frees SD graph */
2056  SCIP* scip, /**< SCIP */
2057  SDGRAPH** sdgraph /**< the SD graph */
2058 )
2059 {
2060  SDGRAPH* g_sd;
2061  assert(scip && sdgraph);
2062
2063  g_sd = *sdgraph;
2064  assert(g_sd);
2065
2066  SCIPfreeMemoryArrayNull(scip, &(g_sd->halfedge_isInMst));
2067  SCIPfreeMemoryArray(scip, &(g_sd->nodemapOrgToDist));
2068  SCIPfreeMemoryArray(scip, &(g_sd->mstsdist));
2069  SCIPfreeMemoryArray(scip, &(g_sd->mstcosts));
2070  SCIPfreeMemoryArray(scip, &(g_sd->sdmst));
2071
2072  graph_free(scip, &(g_sd->distgraph), TRUE);
2073  sdqueryFree(scip, g_sd);
2074
2075  SCIPfreeMemory(scip, sdgraph);
2076 }
2077
2078
2079 /** frees SD graph, but does not free actual graph and node-map (assumed to be non-owned) */
2081  SCIP* scip, /**< SCIP */
2082  SDGRAPH** sdgraph /**< the SD graph */
2083 )
2084 {
2085  SDGRAPH* g_sd;
2086  assert(scip && sdgraph);
2087
2088  g_sd = *sdgraph;
2089  assert(g_sd);
2090  assert(!(g_sd->halfedge_isInMst));
2091
2092  SCIPfreeMemoryArray(scip, &(g_sd->mstsdist));
2093  SCIPfreeMemoryArray(scip, &(g_sd->mstcosts));
2094  SCIPfreeMemoryArray(scip, &(g_sd->sdmst));
2095
2096  sdqueryFree(scip, g_sd);
2097
2098  SCIPfreeMemory(scip, sdgraph);
2099 }
#define STP_Vectype(type)
Definition: stpvector.h:44
int graph_get_nEdges(const GRAPH *)
Definition: graph_stats.c:548
static SCIP_Real sdqueryFullGetSd(int term1, int term2, const SDGRAPH *sdgraph)
void graph_knot_chg(GRAPH *, int, int)
Definition: graph_node.c:86
static volatile int nterms
Definition: interrupt.c:38
static void sdqueryBuildRmq(SDGRAPH *sdgraph, LCABUILDER *lcabuilder)
#define StpVecGetSize(vec)
Definition: stpvector.h:139
SCIP_RETCODE SCIPStpunionfindInit(SCIP *scip, UF *uf, int length)
Definition: misc_stp.c:817
Definition: graphdefs.h:224
int source
Definition: graphdefs.h:195
static SCIP_RETCODE sdqueryRmqInit(SCIP *scip, const GRAPH *g, SDGRAPH *sdgraph)
void reduce_sdgraphInitOrderedMstCosts(SDGRAPH *sdgraph)
static void sdqueryFree(SCIP *scip, SDGRAPH *sdgraph)
#define SCIPfreeMemoryArrayNull(scip, ptr)
Definition: scip_mem.h:72
#define STP_SDQUERYFULL_MAXTERMS
#define Is_term(a)
Definition: graphdefs.h:102
static SCIP_RETCODE sdqueryLcaBuilderInit(SCIP *scip, const SDGRAPH *sdgraph, LCABUILDER **lcabuilder)
struct lowest_common_ancestor_tree_builder LCABUILDER
void graph_tpathsGetClosestTerm(const GRAPH *, const TPATHS *, int, int *RESTRICT, int *RESTRICT, SCIP_Real *RESTRICT)
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
signed int edge
Definition: graphdefs.h:287
SCIP_Bool graph_valid(SCIP *, const GRAPH *)
Definition: graph_base.c:1480
int terms
Definition: graphdefs.h:192
int * nodes_predEdge
Definition: graphdefs.h:329
SCIP_RETCODE reduce_sdgraphMstBuild(SCIP *scip, const GRAPH *g, SDGRAPH *g_sd)
#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_Real reduce_sdgraphGetMaxCost(const SDGRAPH *sdgraph)
static SCIP_RETCODE sdgraphUpdateDistgraphFromTpaths(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, const TPATHS *tpaths, SDGRAPH *g_sd)
SCIP_Bool reduce_sdgraphHasOrderedMstCosts(const SDGRAPH *sdgraph)
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
const int * reduce_sdgraphGetOrgnodesToSdMap(const SDGRAPH *sdgraph)
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_Real * rmq_sparseTable
void graph_tpathsGet3CloseTerms(const GRAPH *, const TPATHS *, int, SCIP_Real, int *RESTRICT, int *RESTRICT, SCIP_Real *RESTRICT, int *RESTRICT)
static void sdgraphMstMarkOrgEdges(const GRAPH *g, const VNOI *vnoi, const int *distedge2org, SDGRAPH *g_sd)
static SCIP_Real distgraphGetBoundaryEdgeDistBest(const GRAPH *g, const TPATHS *tpaths, int tail, int head, SCIP_Real edgecost, const SDPROFIT *sdprofit, int *base_tail, int *base_head)
#define EAT_LAST
Definition: graphdefs.h:58
#define StpVecReserve(scip, vec, _size_)
Definition: stpvector.h:186
#define SCIPdebugMessage
Definition: pub_message.h:87
#define FARAWAY
Definition: graphdefs.h:89
static void sdqueryRmqDfs(int root, const BINARYNODE *lcatree, const SCIP_Real *lcatree_costs, SCIP_Real *rmq_edgecosts, int *lcatreeNodeToRmq, int *rmq_count)
void SCIPsortDownRealInt(SCIP_Real *realarray, int *intarray, int len)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
static void sdgraphSetDefaults(const GRAPH *g, SDGRAPH *g_sd)
void graph_vnoiFree(SCIP *, VNOI **)
Definition: graph_vnoi.c:1348
SCIP_Real miscstp_maxReal(const SCIP_Real *realarr, unsigned nreals)
Definition: misc_stp.c:142
static int sdqueryGetRmqLength(const SDGRAPH *sdgraph)
static void sdgraphMstSortCosts(SDGRAPH *g_sd)
#define STP_TERM_NONE
Definition: graphdefs.h:62
int *RESTRICT mark
Definition: graphdefs.h:198
SCIP_RETCODE graph_vnoiCompute(SCIP *, const GRAPH *, VNOI *)
Definition: graph_vnoi.c:1302
struct binary_tree_node BINARYNODE
#define SCIPallocCleanBufferArray(scip, ptr, num)
Definition: scip_mem.h:133
void reduce_sdgraphInsertEdge(SCIP *scip, int sdnode1, int sdnode2, SCIP_Real edgecost, int edgeid, int *RESTRICT edgeorg, SDGRAPH *sdgraph, SCIP_Bool *success)
int *RESTRICT oeat
Definition: graphdefs.h:231
static void sdgraphFinalize(SCIP *scip, VNOI **vnoi, int **edgeorg)
#define LE(a, b)
Definition: portab.h:83
#define GE(a, b)
Definition: portab.h:85
miscellaneous methods used for solving Steiner problems
void graph_tpathsGetProfitNodes(SCIP *, const GRAPH *, const TPATHS *, const SDPROFIT *, int, int, STP_Vectype(int))
Definition: graph_tpath.c:1842
static SCIP_Real sdgraphGetSd(int term1, int term2, SDGRAPH *sdgraph)
static int distgraphGetNedges(const GRAPH *g)
static SCIP_RETCODE sdqueryInit(SCIP *scip, const GRAPH *g, SDGRAPH *sdgraph)
void reduce_sdgraphFreeFromDistGraph(SCIP *scip, SDGRAPH **sdgraph)
static void sdqueryAttachBinaryTreeNode(const GRAPH *distgraph, int parentposition, int graphnode, LCABUILDER *lcabuilder, UF *uf)
SCIP_RETCODE reduce_sdgraphInitFromDistGraph(SCIP *scip, const GRAPH *g, GRAPH *distgraph, int *node2dist, SDGRAPH **sdgraph)
static SCIP_Real distgraphGetBoundaryEdgeDist2(int tail, int head, int vbase_tail, int vbase_head, SCIP_Real edgecost, SCIP_Real dist_tail, SCIP_Real dist_head, const SDPROFIT *sdprofit)
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
SCIP_Real dist
Definition: graphdefs.h:286
static SCIP_RETCODE sdgraphBuildDistgraph(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, SDGRAPH *g_sd, VNOI **vnoi, int **distedge2org)
Definition: graph_node.c:64
#define NULL
Definition: lpi_spx1.cpp:155
static SCIP_RETCODE sdgraphMstBuild(SCIP *scip, const GRAPH *g, SDGRAPH *g_sd)
int * nodes_base
Definition: graphdefs.h:330
#define Edge_anti(a)
Definition: graphdefs.h:106
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
void SCIPStpunionfindUnion(UF *uf, int p, int q, SCIP_Bool compress)
Definition: misc_stp.c:912
SCIP_RETCODE reduce_sdgraphInit(SCIP *scip, const GRAPH *g, SDGRAPH **sdgraph)
#define MST_MODE
Definition: graphdefs.h:98
int SCIPStpunionfindFind(UF *uf, int element)
Definition: misc_stp.c:885
#define LT(a, b)
Definition: portab.h:81
static void distgraphAddNodes(const GRAPH *g, int *RESTRICT distnodes_id, GRAPH *distgraph)
unsigned char STP_Bool
Definition: portab.h:34
static void sdqueryLcaBuilderFree(SCIP *scip, LCABUILDER **lcabuilder)
static SCIP_Real reduce_sdprofitGetProfit(const SDPROFIT *sdprofit, int node, int nonsource1, int nonsource2)
Definition: reducedefs.h:178
static SCIP_RETCODE sdqueryBuildBinaryTree(SCIP *scip, SDGRAPH *sdgraph, LCABUILDER *lcabuilder)
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
static void sdqueryBuildRmqSparseTable(SDGRAPH *sdgraph)
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT path_heap
Definition: graphdefs.h:255
static SCIP_RETCODE sdgraphAlloc(SCIP *scip, const GRAPH *g, SDGRAPH **sdgraph)
static SCIP_RETCODE sdqueryFullBuild(SCIP *scip, SDGRAPH *sdgraph, LCABUILDER *lcabuilder)
static SCIP_RETCODE sdgraphAllocRestricted(SCIP *scip, const GRAPH *g, SDGRAPH **sdgraph)
int *RESTRICT tail
Definition: graphdefs.h:223
#define MAX(x, y)
Definition: tclique_def.h:83
#define STP_TERM
Definition: graphdefs.h:61
static void distgraphInsertEdge(SCIP *scip, int sdnode1, int sdnode2, SCIP_Real edgecost, int edgeid, int *RESTRICT edgeorg, GRAPH *distgraph, SCIP_Bool *success)
SCIP_Real reduce_sdgraphGetSd(int term1, int term2, SDGRAPH *sdgraph)
int *RESTRICT term
Definition: graphdefs.h:196
const SCIP_Real * reduce_sdgraphGetOrderedMstCosts(const SDGRAPH *sdgraph)
SCIP_RETCODE graph_vnoiInit(SCIP *, const GRAPH *, SCIP_Bool, VNOI **)
Definition: graph_vnoi.c:1265
static void sdqueryBuildNodesToRmqMap(const GRAPH *g, const LCABUILDER *lcabuilder, SDGRAPH *sdgraph)
static long * number
void graph_path_exit(SCIP *, GRAPH *)
Definition: graph_path.c:515
static SCIP_RETCODE sdgraphBuildDistgraphFromTpaths(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, const TPATHS *tpaths, SDGRAPH *g_sd)
static int log2floor(uint32_t number)
Portable definitions.
const STP_Bool * reduce_sdgraphGetMstHalfMark(const SDGRAPH *sdgraph)
SCIP_Bool reduce_sdgraphHasMstHalfMark(const SDGRAPH *sdgraph)
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
static void sdqueryBuildNodesToFullMap(const GRAPH *g, const LCABUILDER *lcabuilder, SDGRAPH *sdgraph)
void reduce_sdgraphFree(SCIP *scip, SDGRAPH **sdgraph)
SCIP_Bool reduce_sdgraphEdgeIsInMst(const SDGRAPH *sdgraph, int edge)
static SCIP_RETCODE sdqueryFullInit(SCIP *scip, const GRAPH *g, SDGRAPH *sdgraph)
SCIP_Real * cost
Definition: graphdefs.h:221
SCIP_Bool graph_edge_isInRange(const GRAPH *, int)
Definition: graph_stats.c:110
void reduce_sdgraphMstSortCosts(SDGRAPH *g_sd)
static const int deBruijnBits[32]
#define SCIP_Real
Definition: def.h:177
int esize
Definition: graphdefs.h:218
#define SCIPfreeCleanBufferArray(scip, ptr)
Definition: scip_mem.h:137
static void sdqueryFullDfs(int root, const BINARYNODE *lcatree, int nlcanodes, const SCIP_Real *lcatree_costs, SCIP_Bool *RESTRICT nodes_isVisited, UF *uf, SCIP_Real *RESTRICT fullq_dists)
#define StpVecFree(scip, vec)
Definition: stpvector.h:153
int *RESTRICT outbeg
Definition: graphdefs.h:204
#define SCIP_Longint
Definition: def.h:162
SCIP_RETCODE reduce_sdgraphInitBiasedFromTpaths(SCIP *scip, GRAPH *g, const SDPROFIT *sdprofit, const TPATHS *tpaths, SDGRAPH **sdgraph)
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
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
#define UNKNOWN
Definition: sepa_mcf.c:4095
static void distgraphAddEdgesFromTpaths(SCIP *scip, const GRAPH *g, const int *distnodes_id, const SDPROFIT *sdprofit, const TPATHS *tpaths, GRAPH *distgraph)
SCIP_Real * nodes_dist
Definition: graphdefs.h:328
#define nnodes
Definition: gastrans.c:65
SCIP_RETCODE graph_vnoiComputeImplied(SCIP *, const GRAPH *, const SDPROFIT *, VNOI *)
Definition: graph_vnoi.c:1323
SCIP_RETCODE graph_init(SCIP *, GRAPH **, int, int, int)
Definition: graph_base.c:607
static SCIP_Real distgraphGetBoundaryEdgeDist(int tail, int head, int vbase_tail, int vbase_head, SCIP_Real edgecost, const SCIP_Real *nodes_vdist, const SDPROFIT *sdprofit)
static SCIP_Real sdqueryGetSd(int term1, int term2, const SDGRAPH *sdgraph)
void SCIPStpunionfindFreeMembers(SCIP *scip, UF *uf)
Definition: misc_stp.c:955
includes various reduction methods for Steiner tree problems
void SCIPsortDownReal(SCIP_Real *realarray, int len)
static void distgraphAddEdges(SCIP *scip, const GRAPH *g, const int *distnodes_id, const VNOI *vnoi, const SDPROFIT *sdprofit, int *RESTRICT edgeorg, GRAPH *distgraph)
SCIP_RETCODE reduce_sdgraphInitBiased(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, SDGRAPH **sdgraph)