Scippy

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 /* */
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_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 */
58  SCIP_Bool mstcostsReady; /**< are the mstcosts already available? */
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 
190  assert(distgraph->head[ne] == tempnode);
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 
201  /* already finished? */
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 
222  /* already visited? */
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];
424  const int head = distgraph->head[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  {
1053  graph_knot_add(distgraph, STP_TERM_NONE);
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 */
1074  int head, /**< head */
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);
1083  const SCIP_Real profit_head = reduce_sdprofitGetProfit(sdprofit, head, 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;
1087  const SCIP_Real distHeadTail = edgecost + nodes_vdist[head] - profit_head;
1088 
1089  const SCIP_Real dist = miscstp_maxReal((SCIP_Real[])
1090  { distAll, distTailHead, distHeadTail, edgecost,
1091  nodes_vdist[tail], nodes_vdist[head] },
1092  6);
1093  return dist;
1094 }
1095 
1096 
1097 /** gets SD for boundary edge */
1098 static inline
1100  int tail, /**< tail */
1101  int head, /**< head */
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);
1111  const SCIP_Real profit_head = reduce_sdprofitGetProfit(sdprofit, head, 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;
1115  const SCIP_Real distHeadTail = edgecost + dist_head - profit_head;
1116 
1117  const SCIP_Real dist = miscstp_maxReal((SCIP_Real[])
1118  { distAll, distTailHead, distHeadTail, edgecost,
1119  dist_tail, dist_head},
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 */
1131  int head, /**< head */
1132  SCIP_Real edgecost, /**< cost */
1133  const SDPROFIT* sdprofit, /**< profit */
1134  int* base_tail, /**< base of tail */
1135  int* base_head /**< base of head */
1136  )
1137 {
1138  SCIP_Real dists_tail[3];
1139  SCIP_Real dists_head[3];
1140  int terms_tail[3];
1141  int terms_head[3];
1142  int nterms_tail;
1143  int nterms_head;
1144  SCIP_Real dist = FARAWAY;
1145 
1146  assert(g && tpaths && sdprofit);
1147 
1148  *base_tail = -1;
1149  *base_head = -1;
1150 
1151  graph_tpathsGet3CloseTerms(g, tpaths, tail, FARAWAY, terms_tail, NULL, dists_tail, &nterms_tail);
1152  graph_tpathsGet3CloseTerms(g, tpaths, head, FARAWAY, terms_head, NULL, dists_head, &nterms_head);
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 =
1161  distgraphGetBoundaryEdgeDist2(tail, head, terms_tail[i], terms_head[j], edgecost, dists_tail[i], dists_head[j], sdprofit);
1162 
1163  if( distnew < dist )
1164  {
1165  *base_tail = terms_tail[i];
1166  *base_head = terms_head[j];
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);
1212  assert(distgraph->head[ne] == sdnode2);
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  {
1272  const int head = g->head[e];
1273  const int vbase_head = nodes_vbase[head];
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 ?
1281  distgraphGetBoundaryEdgeDist(tail, head, vbase_tail, vbase_head, g->cost[e], nodes_vdist, sdprofit)
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]));
1289  assert(Is_term(g->term[vbase_head]));
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;
1320  g_sd->mstcostsReady = FALSE;
1321  g_sd->edgemarkReady = TRUE;
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  {
1413  const int head = g->head[e];
1414  SCIP_Real dist_head;
1415  int vbase_head;
1416 
1417  graph_tpathsGetClosestTerm(g, tpaths, head, &vbase_head, NULL, &dist_head);
1418 
1419  if( vbase_tail != vbase_head )
1420  {
1421  SCIP_Bool success;
1422  const SCIP_Real distance = useProfit ?
1423  distgraphGetBoundaryEdgeDist2(tail, head, vbase_tail, vbase_head, g->cost[e], dist_tail, dist_head, sdprofit)
1424  : (g->cost[e] + dist_tail + dist_head);
1425 
1426  assert(LE(distance, g->cost[e] + dist_tail + dist_head));
1427  assert(Is_term(g->term[vbase_tail]) && Is_term(g->term[vbase_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;
1468  distgraphAddNodes(g, distnodes_id, 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;
1495  distgraphAddNodes(g, distnodes_id, 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;
1516  STP_Vectype(int) profitnodes_head = 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);
1526  StpVecReserve(scip, profitnodes_head, 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  {
1533  const int head = g->head[e];
1534  int vbase_tail;
1535  int vbase_head;
1536  SCIP_Real distance;
1537 
1538  if( head < tail )
1539  continue;
1540 
1541  distance = distgraphGetBoundaryEdgeDistBest(g, tpaths, tail, head, g->cost[e], sdprofit, &vbase_tail, &vbase_head);
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  }
1556  graph_tpathsGetProfitNodes(scip, g, tpaths, sdprofit, head, vbase_head, profitnodes_head);
1557 
1558  for( int k = 0; k < StpVecGetSize(profitnodes_head); k++ )
1559  {
1560  if( hasharr[profitnodes_head[k]] )
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  }
1579  //SCIPdebugMessage("add biased MST edge %d->%d (%f<%f) \n", vbase_tail, vbase_head, distance, sdgraphGetSd(vbase_tail, vbase_head, g_sd));
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);
1595  StpVecFree(scip, profitnodes_head);
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 */
1857  (*sdgraph)->edgemarkReady = FALSE;
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);
1896  if( sdgraph->edgemarkReady )
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 
1931  return sdgraph->mstcostsReady;
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 
1945  if( !sdgraph->mstcostsReady )
1946  {
1947  sdgraphMstSortCosts(sdgraph);
1948  sdgraph->mstcostsReady = TRUE;
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
int *RESTRICT head
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:182
#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)
void graph_knot_add(GRAPH *, int)
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:149
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)