# SCIP

Solving Constraint Integer Programs

graph_sdpath.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 graph_sdpath.c
17  * @brief Special distance path and walk algorithms for Steiner problems
18  * @author Daniel Rehfeldt
19  *
20  * This file encompasses various special distance path (and walk) algorithms for
21  * Steiner tree and related problems.
22  *
23  */
24
25
26 #include <stdlib.h>
27 #include <stdio.h>
28 #include <stddef.h>
29 #include <assert.h>
30 #include "portab.h"
31 #include "graph.h"
32 #include "graphheaps.h"
33 #include "reduce.h"
34
35
36 /** internal data for clique-paths computations */
38 {
39  SCIP_Real* nodes_distBaseToMax;/**< distance from base to maximum profit */
40  SCIP_Real* nodes_distFromMax; /**< distance from maximum profit to current node */
41  SCIP_Real* nodes_maxprofit; /**< maximum profit on path */
42  int* nodes_pred; /**< per node: predecessor */
43  int* nodes_base; /**< per node: base */
44  int* nodes_baseToClique; /**< per base: clique id
45  NOTE: undefined for non-base nodes! */
46 } CLIQUEPATHS;
47
48
49 /** initializes clique paths */
50 static
52  SCIP* scip, /**< SCIP data structure */
53  const GRAPH* g, /**< graph data structure */
54  CLIQUEPATHS* cliquepaths /**< data */
55 )
56 {
57  const int nnodes = graph_get_nNodes(g);
58
59  assert(scip && cliquepaths);
60
61  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_distBaseToMax), nnodes) );
62  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_maxprofit), nnodes) );
63  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_distFromMax), nnodes) );
64  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_pred), nnodes) );
65  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_base), nnodes) );
66  SCIP_CALL( SCIPallocBufferArray(scip, &(cliquepaths->nodes_baseToClique), nnodes) );
67
68  return SCIP_OKAY;
69 }
70
71
72 /** frees clique paths */
73 static
75  SCIP* scip, /**< SCIP data structure */
76  CLIQUEPATHS* cliquepaths /**< data */
77 )
78 {
79  assert(scip && cliquepaths);
80
81  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_baseToClique));
82  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_base));
83  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_pred));
84  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_distFromMax));
85  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_maxprofit));
86  SCIPfreeBufferArray(scip, &(cliquepaths->nodes_distBaseToMax));
87 }
88
89
90 /** gets distance for extension along edge e */
91 inline static
93  const int* prevedges, /**< previous edges per node */
94  const int* nprevedges, /**< number of previous edges per node */
95  const SCIP_Real* cost, /**< cost */
96  const SCIP_Real* dist, /**< distance */
97  int k, /**< previous node */
98  int e, /**< current outgoing edge */
99  int maxnprevs /**< maximum number of previous */
100 )
101 {
102  const int nprevs = nprevedges[k];
103  SCIP_Real dist_e;
104
105  /* ancestor list not full? */
106  if( nprevs != maxnprevs + 1 )
107  {
108  int i;
109  const int e2 = e / 2;
110  assert(nprevs <= maxnprevs);
111
112  /* check whether m is contained in ancestor list */
113  for( i = 0; i < nprevs; i++ )
114  {
115  const int prevedge = prevedges[maxnprevs * k + i];
116
117  if( e2 == prevedge )
118  break;
119  }
120
121  /* e2 in list? */
122  if( i != nprevs )
123  {
124  assert(e2 == prevedges[maxnprevs * k + i]);
125  dist_e = dist[k];
126  }
127  else
128  dist_e = dist[k] + cost[e];
129  }
130  else
131  {
132  dist_e = dist[k] + cost[e];
133  }
134
135  return dist_e;
136 }
137
138
139 inline static
141  const int* prevNPterms, /**< previous np terminals per node */
142  const int* nprevNPterms, /**< number of previous np terminals per node */
143  const int* termmark, /**< terminal mark */
144  const STP_Bool* visited, /**< visited */
145  const SCIP_Real* prize, /**< prize */
146  int k, /**< current node */
147  int m, /**< next node */
148  SCIP_Real distnew, /**< distance of m */
149  int maxnprevs /**< maximum number of previous */
150 )
151 {
152  SCIP_Real distnewP = distnew;
153
154  assert(termmark[m] == 1 || termmark[m] == 2 );
155
156  if( termmark[m] == 2 || !visited[m] )
157  {
158  distnewP = MAX(0.0, distnewP - prize[m]);
159  }
160  else
161  {
162  const int nprevs = nprevNPterms[k];
163
164  /* ancestor list not full? */
165  if( nprevs != maxnprevs + 1 )
166  {
167  int i;
168  assert(nprevs <= maxnprevs);
169
170  /* check whether m is contained in ancestor list */
171  for( i = 0; i < nprevs; i++ )
172  {
173  const int prevterm = prevNPterms[maxnprevs * k + i];
174
175  if( m == prevterm )
176  break;
177  }
178
179  /* m not in list? */
180  if( i == nprevs )
181  distnewP = MAX(0.0, distnewP - prize[m]);
182  }
183  }
184
185  return distnewP;
186 }
187
188
189
190 inline static
192  const GRAPH* g, /**< graph data structure */
193  int node, /**< the node to be updated */
194  int prednode, /**< the predecessor node */
195  int maxnprevs, /**< maximum number of previous terminals to save */
196  const int* prevterms, /**< previous terminals */
197  const int* nprevterms, /**< number of previous terminals */
198  const SCIP_Bool nodevisited /**< visited node already? */
199  )
200 {
201  const int nprevs = nprevterms[prednode];
202  SCIP_Bool conflict = FALSE;
203
204  assert(Is_term(g->term[node]));
205
206  if( !nodevisited )
207  return FALSE;
208
209  if( nprevs > maxnprevs )
210  {
211  assert(nprevs == maxnprevs + 1);
212  return TRUE;
213  }
214
215  for( int i = 0; i < nprevs; i++ )
216  {
217  const int prevterm = prevterms[maxnprevs * prednode + i];
218  assert(prevterm >= 0);
219
220  if( prevterm == node )
221  {
222  conflict = TRUE;
223  break;
224  }
225  }
226
227  return conflict;
228 }
229
231 inline static
233  const GRAPH* g, /**< graph data structure */
234  int node, /**< the node to be updated */
235  int prednode, /**< the predecessor node */
236  int maxnprevs, /**< maximum number of previous terminals to save */
237  int* prevterms, /**< previous terminals */
238  int* nprevterms /**< number of previous terminals */
239  )
240 {
241  const int predsize = nprevterms[prednode];
242  const SCIP_Bool isterm = Is_term(g->term[node]);
243
244  assert(predsize <= maxnprevs + 1);
245
246  if( predsize == maxnprevs + 1 || (isterm && predsize == maxnprevs) )
247  {
248  nprevterms[node] = maxnprevs + 1;
249  }
250  else
251  {
252 #ifndef NDEBUG
253  for( int j = 0; j < predsize; j++ )
254  assert(prevterms[maxnprevs * prednode + j] != node);
255 #endif
256
257  for( int i = 0; i < predsize; i++ )
258  prevterms[maxnprevs * node + i] = prevterms[maxnprevs * prednode + i];
259
260  nprevterms[node] = predsize;
261
262  if( isterm )
263  {
264  assert(predsize < maxnprevs);
265  prevterms[maxnprevs * node + predsize] = node;
266  nprevterms[node]++;
267  }
268
269  assert(nprevterms[node] <= maxnprevs);
270  }
271 }
272
273 /** copies from of predecessor */
274 inline static
276  int node, /**< the node to be updated */
277  int prednode, /**< the predecessor node */
278  int maxnprevs, /**< maximum number of previous terminals to save */
279  int* prev, /**< previous data elements */
280  int* nprev /**< number of previous data elements */
281  )
282 {
283  const int predsize = nprev[prednode];
284
285  assert(predsize <= maxnprevs);
286
287  /* copy data from predecesseor */
288  for( int i = 0; i < predsize; i++ )
289  prev[maxnprevs * node + i] = prev[maxnprevs * prednode + i];
290
291  nprev[node] = predsize;
292 }
293
294 /** update method for second version of SD walks */
295 static
297  const int* termmark, /**< terminal mark */
298  int node, /**< the node to be updated */
299  int prednode, /**< the predecessor node */
300  int edge, /**< the edge to be updated */
301  int maxnprevs, /**< maximum number of previous terminals to save */
302  SCIP_Bool clear, /**< clear arrays */
303  int* prevterms, /**< previous terminals */
304  int* nprevterms, /**< number of previous terminals */
305  int* prevNPterms, /**< previous non-proper terminals */
306  int* nprevNPterms, /**< number of previous non-proper terminals */
307  int* prevedges, /**< previous edges */
308  int* nprevedges /**< number of previous edges */
309  )
310 {
311  int predsize = nprevterms[prednode];
312
313  /*** 1. proper terminals ***/
314
315  /* not enough space? */
316  if( predsize == maxnprevs + 1 || (termmark[node] == 2 && predsize == maxnprevs) )
317  {
318  nprevterms[node] = maxnprevs + 1;
319  }
320  else
321  {
322 #ifndef NDEBUG
323  for( int j = 0; j < predsize; j++ )
324  assert(prevterms[maxnprevs * prednode + j] != node);
325 #endif
326
327  sdwalkUpdateCopy(node, prednode, maxnprevs, prevterms, nprevterms);
328
329  if( termmark[node] == 2 )
330  {
331  assert(predsize < maxnprevs);
332  prevterms[maxnprevs * node + predsize] = node;
333  nprevterms[node]++;
334  }
335
336  assert(nprevterms[node] <= maxnprevs);
337  }
338
339
340  /*** 2. edges ***/
341
342  if( clear )
343  {
344  nprevNPterms[node] = 0;
345  nprevedges[node] = 0;
346  return;
347  }
348
349  predsize = nprevedges[prednode];
350
351  if( predsize >= maxnprevs )
352  {
353  assert(predsize == maxnprevs || predsize == maxnprevs + 1);
354
355  nprevedges[node] = maxnprevs + 1;
356  nprevNPterms[node] = maxnprevs + 1;
357  return;
358  }
359  assert(predsize < maxnprevs);
360
361  sdwalkUpdateCopy(node, prednode, maxnprevs, prevedges, nprevedges);
362
363  prevedges[maxnprevs * node + predsize] = edge / 2;
364  nprevedges[node]++;
365
366  assert(nprevedges[node] <= maxnprevs);
367
368
369  /*** 3. non-proper terminals ***/
370
371  predsize = nprevNPterms[prednode];
372
373  if( predsize == maxnprevs + 1 || (termmark[node] == 1 && predsize == maxnprevs) )
374  {
375  nprevNPterms[node] = maxnprevs + 1;
376  }
377  else
378  {
379  sdwalkUpdateCopy(node, prednode, maxnprevs, prevNPterms, nprevNPterms);
380
381  if( termmark[node] == 1 )
382  {
383  assert(predsize < maxnprevs);
384  prevNPterms[maxnprevs * node + predsize] = node;
385  nprevNPterms[node]++;
386  }
387
388  assert(nprevNPterms[node] <= maxnprevs);
389  }
390 }
391
392 /** resets temporary data */
393 inline static
395  int nvisits, /**< number of visited nodes */
396  const int* visitlist, /**< stores all visited nodes */
397  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
398  int* state, /**< array to indicate whether a node has been scanned */
399  STP_Bool* visited /**< stores whether a node has been visited */
400 )
401 {
402  for( int k = 0; k < nvisits; k++ )
403  {
404  const int node = visitlist[k];
405  assert(node >= 0);
406
407  visited[node] = FALSE;
408  dist[node] = FARAWAY;
409  state[node] = UNKNOWN;
410  }
411 }
412
413
414 /** corrects heap entry */
415 inline static
417  int* RESTRICT heap,
418  int* RESTRICT state,
419  int* RESTRICT count, /* pointer to store the number of elements on the heap */
420  SCIP_Real* RESTRICT pathdist,
421  int node,
422  SCIP_Real newcost
423  )
424 {
425  int c;
426  int j;
427
428  pathdist[node] = newcost;
429
430  if (state[node] == UNKNOWN)
431  {
432  heap[++(*count)] = node;
433  state[node] = (*count);
434  }
435
436  /* Heap shift up
437  */
438  j = state[node];
439  c = j / 2;
440
441  while( (j > 1) && pathdist[heap[c]] > pathdist[heap[j]] )
442  {
443  const int t = heap[c];
444  heap[c] = heap[j];
445  heap[j] = t;
446  state[heap[j]] = j;
447  state[heap[c]] = c;
448  j = c;
449  c = j / 2;
450  }
451 }
452
453
454 /** gets position in Sd array
455  * todo: bad design, should be somewhere central...probably not too time expensive */
456 inline static
458  int ncliquenodes, /**< number of clique nodes */
459  int newnode, /**< new node */
460  int prednode, /**< predecessor node */
461  const SCIP_Real* nodes_dist, /**< distance per node */
462  const int* nodes_base, /**< base per node */
463  const int* baseToClique /**< mapping to clique ID */
464 )
465 {
466  const int base_new = nodes_base[newnode];
467  const int base_pred = nodes_base[prednode];
468  const int id_new = baseToClique[base_new];
469  const int id_pred = baseToClique[base_pred];
470  const int id_min = MIN(id_new, id_pred);
471  const int id_max = MAX(id_new, id_pred);
472  int pos = 0;
473
474  assert(id_min >= 0);
475  assert(id_max <= ncliquenodes);
476  assert(id_min < id_max);
477
478  for( int k = 1; k <= id_min; k++ )
479  {
480  pos += (ncliquenodes - k);
481  }
482
483  pos += (id_max - id_min - 1);
484  assert(pos < ((ncliquenodes) * (ncliquenodes - 1)) / 2);
485
486  return pos;
487 }
488
489
490 /** updates node data */
491 inline static
493  int newnode, /**< new node */
494  int prednode, /**< predecessor node */
495  const SCIP_Real* nodes_dist, /**< distance per node */
496  SCIP_Real newprofit, /**< new profit */
497  SCIP_Real newedgecost, /**< new edge cost */
498  SCIP_Real* RESTRICT nodes_distBaseToMax,/**< nodes to max */
499  SCIP_Real* RESTRICT nodes_distFromMax, /**< maximum bias */
500  SCIP_Real* RESTRICT nodes_maxpathprofit /**< maximum profit */
501  )
502 {
503  assert(newnode >= 0 && prednode >= 0);
504  assert(newnode != prednode);
505  assert(GE(nodes_distBaseToMax[prednode], 0.0));
506  assert(GE(nodes_maxpathprofit[prednode], 0.0));
507  assert(GE(newprofit, 0.0));
508  assert(GE(newedgecost, 0.0));
509
510  if( newprofit > nodes_maxpathprofit[prednode] )
511  {
512  nodes_distBaseToMax[newnode] = nodes_dist[prednode];
513  nodes_distFromMax[newnode] = newedgecost;
514  nodes_maxpathprofit[newnode] = newprofit;
515  }
516  else
517  {
518  SCIP_Real bias = MIN(newedgecost, newprofit);
519  bias = MIN(nodes_distFromMax[prednode], bias);
520
521  nodes_distBaseToMax[newnode] = nodes_distBaseToMax[prednode];
522  nodes_distFromMax[newnode] = nodes_distFromMax[prednode] + newedgecost - bias;
523  nodes_maxpathprofit[newnode] = nodes_maxpathprofit[prednode];
524  }
525 }
526
527
528 /** updates node data */
529 inline static
531  int newnode, /**< new node */
532  int prednode, /**< predecessor node */
533  SCIP_Real newdist, /**< new distance */
534  DHEAP* RESTRICT dheap, /**< heap */
535  SCIP_Real* RESTRICT nodes_dist, /**< distance per node */
536  int* RESTRICT nodes_base, /**< base per node */
537  int* RESTRICT nodes_pred /**< predecessor per node */
538  )
539 {
540  assert(newnode >= 0 && prednode >= 0);
541  assert(newnode != prednode);
542  assert(GE(newdist, 0.0) && LT(newdist, FARAWAY));
543
544  graph_heap_correct(newnode, newdist, dheap);
545
546  nodes_pred[newnode] = prednode;
547  nodes_dist[newnode] = newdist;
548  nodes_base[newnode] = nodes_base[prednode];
549 }
550
551
552 /** gets biased special distance for two profits */
553 static inline
555  SCIP_Real distBaseToBase,
556  SCIP_Real distBaseToProfit1,
557  SCIP_Real distBaseToProfit2,
558  SCIP_Real profit1,
559  SCIP_Real profit2
560  )
561 {
562  const SCIP_Real distProfitToProfit = distBaseToBase - distBaseToProfit1 - distBaseToProfit2;
563  const SCIP_Real sdAll = distBaseToBase - profit1 - profit2;
564  const SCIP_Real sdProfit1 = distBaseToBase - distBaseToProfit2 - profit1;
565  const SCIP_Real sdProfit2 = distBaseToBase - distBaseToProfit1 - profit2;
566  const SCIP_Real sd_final = miscstp_maxReal((SCIP_Real[])
567  { sdAll,
568  sdProfit1, sdProfit2,
569  distProfitToProfit,
570  distBaseToProfit1, distBaseToProfit2
571  },
572  6);
573
574  assert(GE(distBaseToProfit1, 0.0));
575  assert(GE(distBaseToProfit2, 0.0));
576  assert(GE(distProfitToProfit, 0.0));
577  assert(LE(sd_final, distBaseToBase));
578
579  return sd_final;
580 }
581
582
583 /** gets biased special distance for one profit */
584 static inline
586  SCIP_Real distBaseToBase,
587  SCIP_Real distBaseToProfit,
588  SCIP_Real profit
589  )
590 {
591  const SCIP_Real sdAll = distBaseToBase - profit;
592  const SCIP_Real distBaseToProfit2 = distBaseToBase - distBaseToProfit;
593  SCIP_Real sd_final = MAX(distBaseToProfit, distBaseToProfit2);
594
595  if( sdAll > sd_final )
596  sd_final = sdAll;
597
598  assert(GE(distBaseToProfit, 0.0));
599  assert(GE(distBaseToProfit2, 0.0));
600  assert(GE(sd_final, 0.0));
601  assert(LE(sd_final, distBaseToBase));
602
603  return sd_final;
604 }
605
606 /** gets node biased */
607 static inline
609  const SDPROFIT* sdprofit, /**< profit or NULL */
610  int node, /**< node to get bias for */
611  int nextnode, /**< next node */
612  int prevnode, /**< previous node */
613  SCIP_Real edgecost, /**< edge cost */
614  SCIP_Real dist /**< distance */
615  )
616 {
617  SCIP_Real bias;
618  if( node == prevnode )
619  {
620  bias = 0.0;
621  }
622  else
623  {
624  const SCIP_Real profit = reduce_sdprofitGetProfit(sdprofit, node, nextnode, prevnode);
625  bias = MIN(edgecost, profit);
626
627  if( dist < bias )
628  bias = dist;
629  }
630
631  assert(GE(bias, 0.0));
632
633  return bias;
634 }
635
636
637 #if 0
638
639 /** computes SD of path
640  * USE IN DEBUG MODE ONLY! */
641 static
642 void sdCliqueStarGetPathSd(
643  const SDPROFIT* sdprofit, /**< profit or NULL */
644  const GRAPH* g, /**< the graph */
645  int startnode, /**< start node */
646  const int* nodes_pred, /**< predecessor per node */
647  const int* nodes_base, /**< base per node */
648  SCIP_Real* endToMaxDist,
649  SCIP_Real* startToEndDist,
650  SCIP_Real* maxprofit
651  )
652 {
653  SCIP_Real dist = 0.0;
654  const int endnode = nodes_base[startnode];
655  int prevnode = startnode;
656  *maxprofit = 0.0;
657  *endToMaxDist = 0.0;
658
659  printf("go from %d \n", startnode);
660
661  for( int node = startnode; node != endnode; node = nodes_pred[node] )
662  {
663  int e;
664  SCIP_Real profit;
665  SCIP_Real offset;
666  const int nextnode = nodes_pred[node];
667
668  assert(nextnode != node);
669
670  for( e = g->outbeg[node]; e != EAT_LAST; e = g->oeat[e] )
671  {
672  if( g->head[e] == nextnode )
673  break;
674  }
675
676  assert(e != EAT_LAST);
677
678  profit = reduce_sdprofitGetProfit(sdprofit, node, prevnode, nextnode);
679
680  printf("check node %d, profit=%f, edgecost=%f \n", node, profit, g->cost[e]);
681
682  offset = MIN(profit, g->cost[e]);
683  offset = MIN(offset, dist);
684
685  dist += g->cost[e] - offset;
686
687  if( node != startnode && profit - offset > *maxprofit )
688  {
689  *maxprofit = profit - offset;
690  *endToMaxDist = dist;
691  }
692
693  prevnode = node;
694  }
695
696  assert(LE(*endToMaxDist, dist));
697  *endToMaxDist = dist - *endToMaxDist;
698
699  *startToEndDist = dist;
700 }
701
702
703 /** computes SD one combined path
704  * USE IN DEBUG MODE ONLY! */
705 static
706 SCIP_Real sdCliqueStarGetRecomputedSd(
707  const SDPROFIT* sdprofit, /**< profit or NULL */
708  const GRAPH* g, /**< the graph */
709  const int* nodes_pred, /**< predecessor per node */
710  int edge, /**< edge */
711  const SCIP_Real* nodes_dist, /**< distance per node */
712  const int* nodes_base /**< base per node */
713  )
714 {
715  SCIP_Real distBaseToBase;
716  SCIP_Real distBaseToTail;
718  SCIP_Real maxprofit_tail;
720  SCIP_Real distBaseToMax_tail;
722  const int tail = g->tail[edge];
724  const SCIP_Real edgecost = g->cost[edge];
725  SCIP_Real sd_final;
726
728
729  sdCliqueStarGetPathSd(sdprofit, g, tail, nodes_pred, nodes_base, &distBaseToMax_tail, &distBaseToTail, &maxprofit_tail);
731
732  {
733  SCIP_Real profit_tail = reduce_sdprofitGetProfit(sdprofit, tail, head, nodes_pred[tail]);
735  const SCIP_Real bias_tail = sdCliqueStarGetNodeBias(sdprofit, tail, head, nodes_pred[tail], edgecost, distBaseToTail);
737
738  if( bias_tail > bias_head )
739  {
740  profit_tail -= bias_tail;
741  distBaseToBase = distBaseToTail + distBaseToHead + edgecost - bias_tail;
742  }
743  else
744  {
747  }
748
749  assert(GE(profit_tail, 0.0));
751
753  {
756  }
757
758  if( GE(profit_tail, maxprofit_tail) || Is_term(g->term[tail]) )
759  {
760  maxprofit_tail = profit_tail;
761  distBaseToMax_tail = distBaseToTail;
762  }
763  }
764
766
767  if( GT(maxprofit_head, 0.0) && GT(maxprofit_tail, 0.0) )
768  {
770  }
771  else if( GT(maxprofit_tail, 0.0) )
772  {
773  sd_final = sdGet1ProfitDist(distBaseToBase, distBaseToMax_tail, maxprofit_tail);
774  }
775  else if( GT(maxprofit_head, 0.0) )
776  {
778  }
779  else
780  {
781  sd_final = distBaseToBase;
782  }
783
784  printf("sd_final=%f \n", sd_final);
785
786  return sd_final;
787 }
788 #endif
789
790
791 /** helper */
792 inline static
794  const SDPROFIT* sdprofit, /**< profit or NULL */
795  int centernode, /**< center node */
796  int node, /**< node */
797  int neighbor, /**< neighbor node */
798  const SCIP_Real* nodes_dist, /**< distance */
799  const int* nodes_pred, /**< predecessors */
800  const SCIP_Real* nodes_distBaseToMax,/**< distance */
801  const SCIP_Real* nodes_distFromMax, /**< distance */
802  const SCIP_Real* nodes_maxpathprofit,/**< maximum profit */
803  SCIP_Real* distToMax, /**< pointer */
804  SCIP_Real* distFromMax, /**< pointer */
805  SCIP_Real* maxprofit_node, /**< pointer */
806  SCIP_Bool* nodeHasMaxProfit /**< pointer */
807  )
808 {
809  if( node == nodes_pred[node] || centernode == node )
810  {
811  *maxprofit_node = 0.0;
812  assert(centernode == node || EQ(nodes_distBaseToMax[node], 0.0));
813  assert(centernode == node || EQ(nodes_maxpathprofit[node], 0.0));
814  assert(centernode == node || EQ(nodes_distFromMax[node], 0.0));
815  }
816  else
817  {
818  *maxprofit_node = reduce_sdprofitGetProfit(sdprofit, node, neighbor, nodes_pred[node]);
819  }
820
821  /* is the saved profit better? */
822  if( *maxprofit_node < nodes_maxpathprofit[node] )
823  {
824  *distToMax = nodes_distBaseToMax[node];
825  *maxprofit_node = nodes_maxpathprofit[node];
826  *distFromMax = nodes_distFromMax[node];
827  *nodeHasMaxProfit = FALSE;
828  }
829  else
830  {
831  *distToMax = nodes_dist[node];
832  *distFromMax = 0.0;
833  *nodeHasMaxProfit = TRUE;
834  }
835 }
836
837
838 /** updates SD between nodes */
839 inline static
841  const SDPROFIT* sdprofit, /**< profit or NULL */
842  const int* nodes_pred, /**< predecessors */
843  const SCIP_Real* nodes_baseToMaxDist,/**< distance to max */
844  const SCIP_Real* nodes_distFromMax, /**< distance from max */
845  const SCIP_Real* nodes_maxpathprofit,/**< maximum profit */
846  int ncliquenodes, /**< number of clique nodes */
847  int centernode, /**< center node */
848  int newnode, /**< new node */
849  int prednode, /**< predecessor node */
850  SCIP_Real newdist, /**< the new distance to 'newnode' (along prednode) */
851  SCIP_Real edgecost, /**< the edgecost from 'pred' to 'new' */
852  const SCIP_Real* nodes_dist, /**< distance per node */
853  const int* nodes_base, /**< base per node */
854  const int* baseToClique, /**< mapping to clique ID */
855  SCIP_Real* RESTRICT sds /**< to be filled */
856  )
857 {
858  const int sdposition = sdCliqueStarGetSdPosition(ncliquenodes,
859  newnode, prednode, nodes_dist, nodes_base, baseToClique);
860  /* complete distance along new node: */
861  SCIP_Real sdist = nodes_dist[newnode] + newdist;
862
863  if( sdprofit != NULL )
864  {
865  SCIP_Bool nodeHasMaxProfit_pred;
866  SCIP_Real distToMax_pred;
867  SCIP_Real distFromMax_pred;
868  SCIP_Real maxprofit_pred;
869  SCIP_Bool nodeHasMaxProfit_new;
870  SCIP_Real distToMax_new;
871  SCIP_Real distFromMax_new;
872  SCIP_Real maxprofit_new;
873  SCIP_Real distBaseToBase;
874  SCIP_Real sd_final;
875
876  sdCliqueStarGetFinalProfitData(sdprofit, centernode, newnode, prednode, nodes_dist, nodes_pred, nodes_baseToMaxDist,
877  nodes_distFromMax, nodes_maxpathprofit, &distToMax_new, &distFromMax_new, &maxprofit_new, &nodeHasMaxProfit_new);
878
879  sdCliqueStarGetFinalProfitData(sdprofit, centernode, prednode, newnode, nodes_dist, nodes_pred, nodes_baseToMaxDist,
880  nodes_distFromMax, nodes_maxpathprofit, &distToMax_pred, &distFromMax_pred, &maxprofit_pred, &nodeHasMaxProfit_pred);
881
882  distBaseToBase = distToMax_pred + distToMax_new + distFromMax_pred + distFromMax_new + edgecost;
883
884  if( !nodeHasMaxProfit_pred )
885  {
886  SCIP_Real bias = sdCliqueStarGetNodeBias(sdprofit, prednode, newnode, nodes_pred[prednode], edgecost, distFromMax_pred);
887  distBaseToBase -= bias;
888  }
889  else if( !nodeHasMaxProfit_new )
890  {
891  SCIP_Real bias = sdCliqueStarGetNodeBias(sdprofit, newnode, prednode, nodes_pred[newnode], edgecost, distFromMax_new);
892  distBaseToBase -= bias;
893  }
894
895  if( GT(maxprofit_new, 0.0) && GT(maxprofit_pred, 0.0) )
896  {
897  sd_final = sdGet2ProfitsDist(distBaseToBase, distToMax_pred, distToMax_new, maxprofit_pred, maxprofit_new);
898  }
899  else if( GT(maxprofit_pred, 0.0) )
900  {
901  sd_final = sdGet1ProfitDist(distBaseToBase, distToMax_pred, maxprofit_pred);
902  }
903  else if( GT(maxprofit_new, 0.0) )
904  {
905  sd_final = sdGet1ProfitDist(distBaseToBase, distToMax_new, maxprofit_new);
906  }
907  else
908  {
909  sd_final = distBaseToBase;
910  }
911
912  assert(LE(sd_final, sdist));
913  sdist = sd_final;
914  }
915
916  assert(nodes_base[newnode] != nodes_base[prednode]);
917  assert(GE(newdist, 0.0) && LT(newdist, FARAWAY));
918  assert(GE(sdist, 0.0) && LT(sdist, FARAWAY));
919
920  if( sdist < sds[sdposition] )
921  {
922  sds[sdposition] = sdist;
923  }
924 }
925
926
927 /** initializes */
928 static
930  const GRAPH* g, /**< graph data structure */
931  SDCLIQUE* cliquedata, /**< data */
932  CLIQUEPATHS* cliquepaths /**< paths data */
933  )
934 {
935  DIJK* RESTRICT dijkdata = cliquedata->dijkdata;
936  SCIP_Real* RESTRICT nodes_dist = dijkdata->node_distance;
937  int* RESTRICT visitlist = dijkdata->visitlist;
938  DHEAP* RESTRICT dheap = dijkdata->dheap;
939  const int* const cliquenodes = cliquedata->cliquenodes;
940  const int ncliquenodes = cliquedata->ncliquenodes;
941  SCIP_Real* RESTRICT nodes_distBaseToMax = cliquepaths->nodes_distBaseToMax;
942  SCIP_Real* RESTRICT nodes_distFromMax = cliquepaths->nodes_distFromMax;
943  SCIP_Real* RESTRICT nodes_maxprofit = cliquepaths->nodes_maxprofit;
944  int* RESTRICT nodes_base = cliquepaths->nodes_base;
945  int* RESTRICT nodes_pred = cliquepaths->nodes_pred;
946  int* RESTRICT nodes_baseToClique = cliquepaths->nodes_baseToClique;
947  STP_Bool* RESTRICT visited = dijkdata->node_visited;
948
949  assert(nodes_dist && nodes_base && dheap && cliquenodes);
950  assert(nodes_distBaseToMax && nodes_maxprofit && nodes_distFromMax);
951  assert(graph_heap_isClean(dheap));
952  assert(ncliquenodes >= 2);
953  assert(dijkdata->edgelimit >= 0 && dijkdata->edgelimit < 20000);
954
955 #ifndef NDEBUG
956  for( int i = 0; i < g->knots; i++ )
957  {
958  assert(UNKNOWN == dheap->position[i]);
959  assert(EQ(FARAWAY, nodes_dist[i]));
960  assert(!visited[i]);
961  }
962 #endif
963
965  for( int i = 0; i < ncliquenodes; i++ )
966  {
967  const int node = cliquenodes[i];
968
969  assert(0 <= node && node < g->knots);
970
971  visitlist[i] = node;
972  visited[node] = TRUE;
973  nodes_base[node] = node;
974  nodes_pred[node] = node;
975  nodes_dist[node] = 0.0;
976  nodes_baseToClique[node] = i;
977  nodes_distBaseToMax[node] = 0.0;
978  nodes_maxprofit[node] = 0.0;
979  nodes_distFromMax[node] = 0.0;
980
981  graph_heap_correct(node, 0.0, dheap);
982  }
983
984  dijkdata->nvisits = ncliquenodes;
985  assert(dheap->size == ncliquenodes);
986 }
987
988
989 /** returns distance limit */
990 static inline
992  const SDCLIQUE* cliquedata, /**< data */
993  const SCIP_Real* sds /**< to be filled */
994 )
995 {
996  const int ncliquenodes = cliquedata->ncliquenodes;
997  const int nsds = (ncliquenodes * (ncliquenodes - 1)) / 2;
998  SCIP_Real limit = 0.0;
999
1000  assert(nsds >= 1);
1001
1002  for( int i = 0; i < nsds; i++ )
1003  {
1004  assert(GE(sds[i], 0.0) && LT(sds[i], FARAWAY));
1005
1006  if( sds[i] > limit )
1007  limit = sds[i];
1008  }
1009
1010  assert(GT(limit, 0.0));
1011
1012  return limit;
1013 }
1014
1015
1016 /** computes SDs */
1017 static
1019  const GRAPH* g, /**< graph data structure */
1020  const SDPROFIT* sdprofit, /**< profit or NULL */
1021  SDCLIQUE* cliquedata, /**< data */
1022  CLIQUEPATHS* cliquepaths, /**< paths data */
1023  SCIP_Real* RESTRICT sds /**< to be filled */
1024 )
1025 {
1026  DIJK* RESTRICT dijkdata = cliquedata->dijkdata;
1027  SCIP_Real* RESTRICT nodes_dist = dijkdata->node_distance;
1028  DHEAP* RESTRICT dheap = dijkdata->dheap;
1029  STP_Bool* RESTRICT visited = dijkdata->node_visited;
1030  int* RESTRICT visitlist = dijkdata->visitlist;
1031  int* const state = dheap->position;
1032  const SCIP_Real* const gCost = g->cost;
1033  const int* const gOeat = g->oeat;
1035  const int* const nodes_baseToClique = cliquepaths->nodes_baseToClique;
1036  SCIP_Real* RESTRICT nodes_distBaseToMax = cliquepaths->nodes_distBaseToMax;
1037  SCIP_Real* RESTRICT nodes_distFromMax = cliquepaths->nodes_distFromMax;
1038  SCIP_Real* RESTRICT nodes_maxprofit = cliquepaths->nodes_maxprofit;
1039  int* RESTRICT nodes_base = cliquepaths->nodes_base;
1040  int* RESTRICT nodes_pred = cliquepaths->nodes_pred;
1041  const int ncliquenodes = cliquedata->ncliquenodes;
1042  int nvisits = dijkdata->nvisits;
1043  const SCIP_Real distlimit = sdCliqueStarGetDistLimit(cliquedata, sds);
1044  const SCIP_Bool useProfit = (sdprofit != NULL);
1045  const int centernode = cliquedata->centernode;
1046  int nchecks = 0;
1047  const int limit = dijkdata->edgelimit;
1048
1049  assert(g->knots > 1);
1050  assert(dheap->size > 1);
1051
1052  /* until the heap is empty */
1053  while( dheap->size > 0 && nchecks <= limit )
1054  {
1055  const int k = graph_heap_deleteMinReturnNode(dheap);
1056  const int k_base = nodes_base[k];
1057  const int k_predNode = nodes_pred[k];
1058  const SCIP_Real k_dist = nodes_dist[k];
1059
1060  assert(0 <= k_base && k_base < g->knots);
1061  assert(CONNECT == state[k]);
1062  assert(CONNECT == state[k_base]);
1063
1064  for( int e = g->outbeg[k]; e >= 0; e = gOeat[e] )
1065  {
1066  const int m = gHead[e];
1067  SCIP_Real profit = 0.0;
1068  SCIP_Real bias = 0.0;
1069  SCIP_Real newdist = k_dist + gCost[e];
1070
1071  nchecks++;
1072
1073  /* NOTE: need to make sure that we do not go over the center of the clique!
1074  * todo: Might be an issue if we pseudo-eliminate edges...probably need to block the edges as well */
1075  if( useProfit && k != centernode && k != k_predNode && m != k_predNode && m != centernode )
1076  {
1077  profit = reduce_sdprofitGetProfit(sdprofit, k, k_predNode, m);
1078  bias = MIN(gCost[e], profit);
1079  bias = MIN(k_dist, bias);
1080  newdist -= bias;
1081
1082  assert(k != k_base || EQ(MIN(k_dist, bias), 0.0));
1083  assert(GE(newdist, 0.0));
1084  }
1085
1086
1087  if( GT(newdist, distlimit) )
1088  continue;
1089
1090  /* first time visit of m? */
1091  if( UNKNOWN == state[m] )
1092  {
1093  assert(!visited[m]);
1094
1095  sdCliqueStarUpdateNodeDist(m, k, newdist, dheap, nodes_dist, nodes_base, nodes_pred);
1096  sdCliqueStarUpdateNodeMaxDist(m, k, nodes_dist, profit, gCost[e], nodes_distBaseToMax,
1097  nodes_distFromMax, nodes_maxprofit);
1098
1099  visitlist[nvisits++] = m;
1100  visited[m] = TRUE;
1101
1102  assert(k_base == nodes_base[m]);
1103  continue;
1104  }
1105
1106  assert(visited[m]);
1107
1108  /* can we update the special distances? */
1109  if( nodes_base[m] != k_base )
1110  {
1111  sdCliqueStarUpdateSd(sdprofit, nodes_pred, nodes_distBaseToMax, nodes_distFromMax, nodes_maxprofit,
1112  ncliquenodes, centernode, m, k, newdist, gCost[e], nodes_dist, nodes_base, nodes_baseToClique, sds);
1113  }
1114
1115  if( state[m] != CONNECT )
1116  {
1117  /* check whether the path (to m) including k is shorter than the so far best known */
1118  if( nodes_dist[m] > newdist )
1119  {
1120  sdCliqueStarUpdateNodeDist(m, k, newdist, dheap, nodes_dist, nodes_base, nodes_pred);
1121  sdCliqueStarUpdateNodeMaxDist(m, k, nodes_dist, profit, gCost[e], nodes_distBaseToMax,
1122  nodes_distFromMax, nodes_maxprofit);
1123  }
1124  }
1125  }
1126  }
1127
1128  dijkdata->nvisits = nvisits;
1129 }
1130
1131
1132 /*
1133  * Interface methods
1134  */
1135
1136 /** limited Dijkstra, stopping at terminals */
1138  SCIP* scip, /**< SCIP data structure */
1139  const GRAPH* g, /**< graph data structure */
1140  SCIP_Bool with_zero_edges, /**< telling name */
1141  int star_root, /**< root of the start */
1142  int edgelimit, /**< maximum number of edges to consider during execution */
1143  int* star_base, /**< star base node, must be initially set to SDSTAR_BASE_UNSET */
1144  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
1145  int* visitlist, /**< stores all visited nodes */
1146  int* nvisits, /**< number of visited nodes */
1147  DHEAP* dheap, /**< Dijkstra heap */
1148  STP_Bool* visited, /**< stores whether a node has been visited */
1149  SCIP_Bool* success /**< will be set to TRUE iff at least one edge can be deleted */
1150  )
1151 {
1152  int nchecks;
1153  int nstarhits;
1154  int* const state = dheap->position;
1155  DCSR* const dcsr = g->dcsr_storage;
1156  const RANGE* const RESTRICT range_csr = dcsr->range;
1158  const SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
1159  const int star_degree = range_csr[star_root].end - range_csr[star_root].start;
1160  SCIP_Real distlimit;
1161  /* NOTE: with zero edges case is already covered with state[k] = UNKNOWN if k == star_base[k] */
1162  const SCIP_Real eps = graph_pc_isPcMw(g) ? 0.0 : SCIPepsilon(scip);
1163
1164  assert(dcsr && g && dist && visitlist && nvisits && visited && dheap && success);
1165  assert(graph_pc_isPcMw(g));
1166  assert(!g->extended);
1167  assert(g->mark[star_root] && star_degree >= 1);
1168  assert(dheap->size == 0);
1169  assert(edgelimit >= 1);
1170
1171  *nvisits = 0;
1172  *success = FALSE;
1173
1174 #ifndef NDEBUG
1175  for( int k = 0; k < g->knots; k++ )
1176  {
1177  assert(dist[k] == FARAWAY);
1178  assert(star_base[k] == SDSTAR_BASE_UNSET);
1179  assert(state[k] == UNKNOWN);
1180  }
1181 #endif
1182
1183  distlimit = 0.0;
1184  dist[star_root] = 0.0;
1185  state[star_root] = CONNECT;
1186  visitlist[(*nvisits)++] = star_root;
1187
1188  for( int e = range_csr[star_root].start, end = range_csr[star_root].end; e < end; e++ )
1189  {
1190  const int m = head_csr[e];
1191
1192  assert(g->mark[m]);
1193  assert(!visited[m]);
1194
1195  visitlist[(*nvisits)++] = m;
1196  visited[m] = TRUE;
1197  dist[m] = cost_csr[e];
1198  star_base[m] = m;
1199
1200  /* add epsilon to make sure that m is removed from the heap last in case of equality */
1201  graph_heap_correct(m, cost_csr[e] + eps, dheap);
1202
1203  if( cost_csr[e] > distlimit )
1204  distlimit = cost_csr[e];
1205  }
1206
1207
1208  nchecks = 0;
1209  nstarhits = 0;
1210
1211  while( dheap->size > 0 && nchecks <= edgelimit )
1212  {
1213  /* get nearest labeled node */
1214  const int k = graph_heap_deleteMinReturnNode(dheap);
1215  const int k_start = range_csr[k].start;
1216  const int k_end = range_csr[k].end;
1217
1218  assert(k != star_root);
1219  assert(state[k] == CONNECT);
1220  assert(LE(dist[k], distlimit));
1221
1222  if( with_zero_edges && k == star_base[k] )
1223  state[k] = UNKNOWN;
1224
1225  /* correct incident nodes */
1226  for( int e = k_start; e < k_end; e++ )
1227  {
1228  const int m = head_csr[e];
1229
1230  assert(g->mark[m] && star_base[k] >= 0);
1231
1232  if( state[m] != CONNECT )
1233  {
1234  const SCIP_Real distnew = dist[k] + cost_csr[e];
1235
1236  if( GT(distnew, distlimit) )
1237  continue;
1238
1239  if( LT(distnew, dist[m]) )
1240  {
1241  if( !visited[m] )
1242  {
1243  visitlist[(*nvisits)++] = m;
1244  visited[m] = TRUE;
1245  }
1246
1247  if( star_base[m] == m )
1248  nstarhits++;
1249
1250  dist[m] = distnew;
1251  star_base[m] = star_base[k];
1252  graph_heap_correct(m, distnew, dheap);
1253
1254  assert(star_base[m] != m);
1255  }
1256  else if( EQ(distnew, dist[m]) && star_base[m] == m )
1257  {
1258  if( with_zero_edges && star_base[k] == star_base[m] )
1259  continue;
1260
1261  assert(visited[m]);
1262  nstarhits++;
1263
1264  assert(star_base[m] != star_base[k]);
1265
1266  dist[m] = distnew;
1267  star_base[m] = star_base[k];
1268  graph_heap_correct(m, distnew, dheap);
1269
1270  assert(star_base[m] != m);
1271  }
1272
1273  /* all star nodes hit already? */
1274  if( nstarhits == star_degree )
1275  {
1276  nchecks = edgelimit + 1;
1277  break;
1278  }
1279  }
1280  nchecks++;
1281  }
1282  }
1283
1284  *success = (nstarhits > 0);
1285 }
1286
1287
1288 /** limited Dijkstra with node bias */
1290  SCIP* scip, /**< SCIP data structure */
1291  const GRAPH* g, /**< graph data structure */
1292  const SDPROFIT* sdprofit, /**< SD profit */
1293  int star_root, /**< root of the start */
1294  int* star_base, /**< star base node, must be initially set to SDSTAR_BASE_UNSET */
1295  DIJK* dijkdata, /**< Dijkstra data */
1296  SCIP_Bool* success /**< will be set to TRUE iff at least one edge can be deleted */
1297  )
1298 {
1299  int nchecks;
1300  int nstarhits;
1301  const int nnodes = graph_get_nNodes(g);
1302  SCIP_Real* RESTRICT dist = dijkdata->node_distance;
1303  int* RESTRICT visitlist = dijkdata->visitlist;
1304  STP_Bool* RESTRICT visited = dijkdata->node_visited;
1305  int* node_predNode;
1306  DHEAP* dheap = dijkdata->dheap;
1307  int* const state = dheap->position;
1308  DCSR* const dcsr = g->dcsr_storage;
1309  const RANGE* const RESTRICT range_csr = dcsr->range;
1311  const SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
1312  const int star_degree = range_csr[star_root].end - range_csr[star_root].start;
1313  int nvisits = 0;
1314  SCIP_Real distlimit;
1315  const int edgelimit = dijkdata->edgelimit;
1316  /* NOTE: with zero edges case is already covered with state[k] = UNKNOWN if k == star_base[k] */
1317  const SCIP_Real eps = graph_pc_isPcMw(g) ? 0.0 : 2.0 * SCIPepsilon(scip);
1318
1319  assert(dcsr && dist && visitlist && visited && dheap && success && sdprofit);
1320  assert(!g->extended);
1321  assert(g->mark[star_root] && star_degree >= 1);
1322  assert(dheap->size == 0);
1323  assert(edgelimit >= 1);
1324
1325  nvisits = 0;
1326  *success = FALSE;
1327
1328  SCIP_CALL( SCIPallocBufferArray(scip, &node_predNode, nnodes) );
1329
1330 #ifndef NDEBUG
1331  for( int k = 0; k < nnodes; k++ )
1332  {
1333  assert(dist[k] == FARAWAY);
1334  assert(star_base[k] == SDSTAR_BASE_UNSET);
1335  assert(state[k] == UNKNOWN);
1336  node_predNode[k] = UNKNOWN;
1337  }
1338 #endif
1339
1340  distlimit = 0.0;
1341  dist[star_root] = 0.0;
1342  state[star_root] = CONNECT;
1343  visitlist[(nvisits)++] = star_root;
1344
1345  for( int e = range_csr[star_root].start, end = range_csr[star_root].end; e < end; e++ )
1346  {
1347  const int m = head_csr[e];
1348
1349  assert(g->mark[m]);
1350  assert(!visited[m]);
1351
1352  visitlist[(nvisits)++] = m;
1353  visited[m] = TRUE;
1354  dist[m] = cost_csr[e];
1355  star_base[m] = m;
1356  node_predNode[m] = star_root;
1357
1358  /* add epsilon to make sure that m is removed from the heap last in case of equality */
1359  graph_heap_correct(m, cost_csr[e] + eps, dheap);
1360
1361  if( cost_csr[e] > distlimit )
1362  distlimit = cost_csr[e];
1363  }
1364
1365  nchecks = 0;
1366  nstarhits = 0;
1367
1368  while( dheap->size > 0 && nchecks <= edgelimit )
1369  {
1370  /* get nearest labeled node */
1371  const int k = graph_heap_deleteMinReturnNode(dheap);
1372  const int k_start = range_csr[k].start;
1373  const int k_end = range_csr[k].end;
1374  const int k_pred = node_predNode[k];
1375
1376  assert(k != star_root);
1377  assert(k_pred >= 0 && k_pred < nnodes);
1378  assert(state[k] == CONNECT);
1379  assert(LE(dist[k], distlimit));
1380
1381  if( k == star_base[k] )
1382  state[k] = UNKNOWN;
1383
1384  /* correct incident nodes */
1385  for( int e = k_start; e < k_end; e++ )
1386  {
1387  const int m = head_csr[e];
1388  assert(g->mark[m] && star_base[k] >= 0);
1389
1390  if( state[m] != CONNECT )
1391  {
1392  SCIP_Real distnew = dist[k] + cost_csr[e];
1393  if( m != k_pred )
1394  {
1395  SCIP_Real profitBias = reduce_sdprofitGetProfit(sdprofit, k, m, k_pred);
1396  profitBias = MIN(profitBias, cost_csr[e]);
1397  profitBias = MIN(profitBias, dist[k]);
1398  distnew -= profitBias;
1399  }
1400
1401  if( GT(distnew, distlimit) )
1402  continue;
1403
1404  if( LT(distnew, dist[m]) )
1405  {
1406  if( !visited[m] )
1407  {
1408  visitlist[(nvisits)++] = m;
1409  visited[m] = TRUE;
1410  }
1411
1412  if( star_base[m] == m )
1413  nstarhits++;
1414
1415  node_predNode[m] = k;
1416  dist[m] = distnew;
1417  star_base[m] = star_base[k];
1418  graph_heap_correct(m, distnew, dheap);
1419
1420  assert(star_base[m] != m && m != star_root);
1421  }
1422  else if( EQ(distnew, dist[m]) && star_base[m] == m )
1423  {
1424  if( star_base[k] == star_base[m] )
1425  continue;
1426
1427  assert(visited[m]);
1428  nstarhits++;
1429
1430  assert(star_base[m] != star_base[k]);
1431
1432  node_predNode[m] = k;
1433  dist[m] = distnew;
1434  star_base[m] = star_base[k];
1435  graph_heap_correct(m, distnew, dheap);
1436
1437  assert(star_base[m] != m && m != star_root);
1438  }
1439
1440  /* all star nodes hit already? */
1441  if( nstarhits == star_degree )
1442  {
1443  nchecks = edgelimit + 1;
1444  break;
1445  }
1446  }
1447  nchecks++;
1448  }
1449  }
1450
1451  dijkdata->nvisits = nvisits;
1452  *success = (nstarhits > 0);
1453
1454  SCIPfreeBufferArray(scip, &node_predNode);
1455
1456  return SCIP_OKAY;
1457 }
1458
1459
1460 /** limited Dijkstra with node bias that stops at terminals */
1462  SCIP* scip, /**< SCIP data structure */
1463  const GRAPH* g, /**< graph data structure */
1464  const SDPROFIT* sdprofit, /**< SD profit or NULL */
1465  int sourcenode, /**< node to start from */
1466  DIJK* dijkdata /**< Dijkstra data */
1467  )
1468 {
1469  int nchecks;
1470  const int nnodes = graph_get_nNodes(g);
1471  SCIP_Real* RESTRICT dist = dijkdata->node_distance;
1472  int* RESTRICT visitlist = dijkdata->visitlist;
1473  STP_Bool* RESTRICT visited = dijkdata->node_visited;
1474  int* RESTRICT node_predNode;
1475  DHEAP* dheap = dijkdata->dheap;
1476  int* const state = dheap->position;
1477  DCSR* const dcsr = g->dcsr_storage;
1478  const RANGE* const RESTRICT range_csr = dcsr->range;
1480  const SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
1481  int nvisits = 0;
1482  const int edgelimit = dijkdata->edgelimit;
1483
1484  assert(dcsr && dist && visitlist && visited && dheap);
1485  assert(!g->extended);
1486  assert(dheap->size == 0);
1487  assert(edgelimit >= 1);
1488
1489  nvisits = 0;
1490
1491  SCIP_CALL( SCIPallocBufferArray(scip, &node_predNode, nnodes) );
1492
1493 #ifndef NDEBUG
1494  for( int k = 0; k < nnodes; k++ )
1495  {
1496  assert(EQ(dist[k], FARAWAY));
1497  assert(state[k] == UNKNOWN);
1498  }
1499 #endif
1500
1501  dist[sourcenode] = 0.0;
1502  visitlist[(nvisits)++] = sourcenode;
1503  graph_heap_correct(sourcenode, 0.0, dheap);
1504  node_predNode[sourcenode] = -1;
1505
1506  nchecks = 0;
1507
1508  while( dheap->size > 0 && nchecks <= edgelimit )
1509  {
1510  /* get nearest labeled node */
1511  const int k = graph_heap_deleteMinReturnNode(dheap);
1512  const int k_start = range_csr[k].start;
1513  const int k_end = range_csr[k].end;
1514  const int k_pred = node_predNode[k];
1515
1516  assert(state[k] == CONNECT);
1517
1518  /* correct incident nodes */
1519  for( int e = k_start; e < k_end; e++ )
1520  {
1521  const int m = head_csr[e];
1522
1523  if( state[m] != CONNECT )
1524  {
1525  SCIP_Real distnew = dist[k] + cost_csr[e];
1526  if( sdprofit && m != k_pred && k != sourcenode )
1527  {
1528  SCIP_Real profitBias = reduce_sdprofitGetProfit(sdprofit, k, m, k_pred);
1529  profitBias = MIN(profitBias, cost_csr[e]);
1530  profitBias = MIN(profitBias, dist[k]);
1531  distnew -= profitBias;
1532  }
1533
1534  if( LT(distnew, dist[m]) )
1535  {
1536  if( !visited[m] )
1537  {
1538  visitlist[(nvisits)++] = m;
1539  visited[m] = TRUE;
1540  }
1541
1542  node_predNode[m] = k;
1543  dist[m] = distnew;
1544
1545  if( !Is_term(g->term[m]) )
1546  {
1547  graph_heap_correct(m, distnew, dheap);
1548  }
1549  }
1550
1551  }
1552  nchecks++;
1553  }
1554  }
1555
1556  dijkdata->nvisits = nvisits;
1557  SCIPfreeBufferArray(scip, &node_predNode);
1558
1559  return SCIP_OKAY;
1560 }
1561
1562
1563 /** modified Dijkstra along walks for PcMw, returns special distance between start and end */
1565  SCIP* scip, /**< SCIP data structure */
1566  const GRAPH* g, /**< graph data structure */
1567  const SCIP_Real* cost, /**< edge costs */
1568  const int* termmark, /**< terminal mark (2 for proper terminal, 1 for non-proper terminal, 0 otherwise)*/
1569  SCIP_Real distlimit, /**< distance limit of the search */
1570  int start, /**< start */
1571  int end, /**< end */
1572  int edgelimit, /**< maximum number of edges to consider during execution */
1573  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
1574  int* heap, /**< array representing a heap */
1575  int* state, /**< array to indicate whether a node has been scanned */
1576  int* visitlist, /**< stores all visited nodes */
1577  int* nvisits, /**< number of visited nodes */
1578  STP_Bool* visited /**< stores whether a node has been visited */
1579  )
1580 {
1581  int count;
1582  int nchecks;
1583  SCIP_Bool success = FALSE;
1584  const int edgelimit1 = edgelimit / 2;
1585
1586  assert(g != NULL);
1587  assert(heap != NULL);
1588  assert(dist != NULL);
1589  assert(cost != NULL);
1590  assert(visitlist != NULL);
1591  assert(nvisits != NULL);
1592  assert(visited != NULL);
1593  assert(graph_pc_isPcMw(g));
1594  assert(!g->extended);
1595
1596  *nvisits = 0;
1597
1599  return FALSE;
1600
1601  assert(g->mark[start] && g->mark[end]);
1602
1603  count = 0;
1604  nchecks = 0;
1605  dist[start] = 0.0;
1606  state[start] = CONNECT;
1607  visitlist[(*nvisits)++] = start;
1608
1609  g->mark[start] = FALSE;
1610  g->mark[end] = FALSE;
1611
1612  for( int e = g->outbeg[start]; e != EAT_LAST; e = g->oeat[e] )
1613  {
1614  const int m = g->head[e];
1615
1616  if( g->mark[m] && SCIPisLE(scip, cost[e], distlimit) )
1617  {
1618  assert(!visited[m]);
1619
1620  visitlist[(*nvisits)++] = m;
1621  visited[m] = TRUE;
1622
1623  if( termmark[m] != 0 )
1624  {
1625  const SCIP_Real newcost = MAX(cost[e] - g->prize[m], 0.0);
1626  sdwalkCorrectHeap(heap, state, &count, dist, m, newcost);
1627  }
1628  else
1629  {
1630  sdwalkCorrectHeap(heap, state, &count, dist, m, cost[e]);
1631  }
1632
1633  if( ++nchecks > edgelimit1 )
1634  break;
1635  }
1636  }
1637  g->mark[end] = TRUE;
1638
1639  while( count > 0 && nchecks <= edgelimit )
1640  {
1641  /* get nearest labelled node */
1642  const int k = nearestX(heap, state, &count, dist);
1643  assert(k != end && k != start);
1644  assert(SCIPisLE(scip, dist[k], distlimit));
1645
1646  if( termmark[k] == 2 )
1647  state[k] = CONNECT;
1648  else
1649  state[k] = UNKNOWN;
1650
1651  /* correct incident nodes */
1652  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
1653  {
1654  const int m = g->head[e];
1655
1656  if( (state[m] != CONNECT) && g->mark[m] )
1657  {
1658  SCIP_Real distnew = dist[k] + cost[e];
1659
1660  if( SCIPisGT(scip, distnew, distlimit) )
1661  continue;
1662
1663  if( termmark[m] != 0 )
1664  distnew = MAX(distnew - g->prize[m], 0.0);
1665
1666  if( distnew < dist[m] )
1667  {
1668  if( !visited[m] )
1669  {
1670  visitlist[(*nvisits)++] = m;
1671  visited[m] = TRUE;
1672  }
1673
1675  if( m == end )
1676  {
1677  nchecks = edgelimit + 1;
1678  success = TRUE;
1679  break;
1680  }
1681
1682  sdwalkCorrectHeap(heap, state, &count, dist, m, distnew);
1683  }
1684  }
1685  nchecks++;
1686  }
1687  }
1688
1689  g->mark[start] = TRUE;
1690  return success;
1691 }
1692
1693
1694
1695 /** modified Dijkstra along walks for PcMw, returns special distance between start and end */
1697  SCIP* scip, /**< SCIP data structure */
1698  const GRAPH* g, /**< graph data structure */
1699  const int* termmark, /**< terminal mark (2 for proper terminal, 1 for non-proper terminal, 0 otherwise)*/
1700  SCIP_Real distlimit, /**< distance limit of the search */
1701  int start, /**< start */
1702  int end, /**< end */
1703  int edgelimit, /**< maximum number of edges to consider during execution */
1704  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
1705  int* visitlist, /**< stores all visited nodes */
1706  int* nvisits, /**< number of visited nodes */
1707  DHEAP* dheap, /**< Dijkstra heap */
1708  STP_Bool* visited /**< stores whether a node has been visited */
1709  )
1710 {
1711  int nchecks;
1712  SCIP_Bool success = FALSE;
1713  const int edgelimit1 = edgelimit / 2;
1714  int* const state = dheap->position;
1715  DCSR* const dcsr = g->dcsr_storage;
1716  RANGE* const RESTRICT range_csr = dcsr->range;
1718  SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
1719
1720  assert(dcsr && g && dist && visitlist && nvisits && visited && dheap);
1721  assert(graph_pc_isPcMw(g));
1722  assert(!g->extended);
1724  assert(g->mark[start] && g->mark[end]);
1725  assert(dheap->size == 0);
1726
1727  *nvisits = 0;
1728
1729 #ifndef NDEBUG
1730  for( int k = 0; k < g->knots; k++ )
1731  assert(state[k] == UNKNOWN);
1732 #endif
1733
1734  nchecks = 0;
1735  dist[start] = 0.0;
1736  state[start] = CONNECT;
1737  visitlist[(*nvisits)++] = start;
1738
1739  for( int e = range_csr[start].start; e < range_csr[start].end; e++ )
1740  {
1741  const int m = head_csr[e];
1742  assert(g->mark[m]);
1743
1744  if( SCIPisLE(scip, cost_csr[e], distlimit) && m != end )
1745  {
1746  assert(!visited[m]);
1747
1748  visitlist[(*nvisits)++] = m;
1749  visited[m] = TRUE;
1750
1751  if( termmark[m] != 0 )
1752  {
1753  const SCIP_Real newcost = MAX(cost_csr[e] - g->prize[m], 0.0);
1754
1755  dist[m] = newcost;
1756  graph_heap_correct(m, newcost, dheap);
1757  }
1758  else
1759  {
1760  dist[m] = cost_csr[e];
1761  graph_heap_correct(m, cost_csr[e], dheap);
1762  }
1763
1764  if( ++nchecks > edgelimit1 )
1765  break;
1766  }
1767  }
1768
1769  while( dheap->size > 0 && nchecks <= edgelimit )
1770  {
1771  /* get nearest labelled node */
1772  const int k = graph_heap_deleteMinReturnNode(dheap);
1773  const int k_start = range_csr[k].start;
1774  const int k_end = range_csr[k].end;
1775
1776  assert(k != end && k != start);
1777  assert(SCIPisLE(scip, dist[k], distlimit));
1778
1779  if( termmark[k] == 2 )
1780  state[k] = CONNECT;
1781  else
1782  state[k] = UNKNOWN;
1783
1784  /* correct incident nodes */
1785  for( int e = k_start; e < k_end; e++ )
1786  {
1787  const int m = head_csr[e];
1788
1789  if( state[m] != CONNECT && m != start )
1790  {
1791  SCIP_Real distnew = dist[k] + cost_csr[e];
1792
1793  assert(g->mark[m]);
1794
1795  if( SCIPisGT(scip, distnew, distlimit) )
1796  continue;
1797
1798  if( termmark[m] != 0 )
1799  distnew = MAX(distnew - g->prize[m], 0.0);
1800
1801  if( LT(distnew, dist[m]) )
1802  {
1803  if( !visited[m] )
1804  {
1805  visitlist[(*nvisits)++] = m;
1806  visited[m] = TRUE;
1807  }
1808
1810  if( m == end )
1811  {
1812  nchecks = edgelimit + 1;
1813  success = TRUE;
1814  break;
1815  }
1816
1817  dist[m] = distnew;
1818  graph_heap_correct(m, distnew, dheap);
1819  }
1820  }
1821  nchecks++;
1822  }
1823  }
1824
1825  return success;
1826 }
1827
1828
1829 /** modified Dijkstra along walks for PcMw, returns special distance between start and end */
1831  SCIP* scip, /**< SCIP data structure */
1832  const GRAPH* g, /**< graph data structure */
1833  const int* termmark, /**< terminal mark (2 for proper terminal, 1 for non-proper terminal, 0 otherwise) */
1834  const int* stateprev, /**< state of previous run or NULL */
1835  SCIP_Real distlimit, /**< distance limit of the search */
1836  int start, /**< start */
1837  int end, /**< end */
1838  int edgelimit, /**< maximum number of edges to consider during execution */
1839  SCIP_Real* prizeoffset, /**< array for storing prize offset or NULL */
1840  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
1841  int* visitlist, /**< stores all visited nodes */
1842  int* nvisits, /**< number of visited nodes */
1843  DHEAP* dheap, /**< Dijkstra heap */
1844  STP_Bool* visited /**< stores whether a node has been visited */
1845  )
1846 {
1847  int nchecks;
1848  SCIP_Bool success = FALSE;
1849  const int edgelimit1 = edgelimit / 2;
1850  int* const state = dheap->position;
1851  DCSR* const dcsr = g->dcsr_storage;
1852  RANGE* const RESTRICT range_csr = dcsr->range;
1854  SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
1855
1856  assert(dcsr && g && dist && visitlist && visited && dheap);
1857  assert(graph_pc_isPcMw(g));
1858  assert(!g->extended);
1860  assert(g->mark[start] && g->mark[end]);
1861  assert(dheap->size == 0);
1862
1863  *nvisits = 0;
1864
1865 #ifndef NDEBUG
1866  for( int k = 0; k < g->knots; k++ )
1867  assert(state[k] == UNKNOWN);
1868 #endif
1869
1870  nchecks = 0;
1871  dist[start] = 0.0;
1872  state[start] = CONNECT;
1873  visitlist[(*nvisits)++] = start;
1874
1875  for( int e = range_csr[start].start; e < range_csr[start].end; e++ )
1876  {
1877  const int m = head_csr[e];
1878  assert(g->mark[m]);
1879
1880  if( SCIPisLE(scip, cost_csr[e], distlimit) && m != end )
1881  {
1882  assert(!visited[m]);
1883
1884  if( stateprev && stateprev[m] == CONNECT )
1885  continue;
1886
1887  visitlist[(*nvisits)++] = m;
1888  visited[m] = TRUE;
1889
1890  if( termmark[m] != 0 )
1891  {
1892  const SCIP_Real newcost = MAX(cost_csr[e] - g->prize[m], 0.0);
1893
1894  dist[m] = newcost;
1895  graph_heap_correct(m, newcost, dheap);
1896
1897  if( prizeoffset )
1898  {
1899  if( g->prize[m] > cost_csr[e] )
1900  {
1901  prizeoffset[m] = cost_csr[e];
1902  assert(SCIPisZero(scip, newcost));
1903  }
1904  else
1905  prizeoffset[m] = g->prize[m];
1906  }
1907  }
1908  else
1909  {
1910  dist[m] = cost_csr[e];
1911  graph_heap_correct(m, cost_csr[e], dheap);
1912  }
1913
1914  if( ++nchecks > edgelimit1 )
1915  break;
1916  }
1917  }
1918
1919  while( dheap->size > 0 && nchecks <= edgelimit )
1920  {
1921  /* get nearest labelled node */
1922  const int k = graph_heap_deleteMinReturnNode(dheap);
1923  const int k_start = range_csr[k].start;
1924  const int k_end = range_csr[k].end;
1925
1926  assert(k != end && k != start);
1927  assert(SCIPisLE(scip, dist[k], distlimit));
1928
1929  if( termmark[k] == 2 )
1930  state[k] = CONNECT;
1931  else
1932  state[k] = UNKNOWN;
1933
1934  /* correct incident nodes */
1935  for( int e = k_start; e < k_end; e++ )
1936  {
1937  const int m = head_csr[e];
1938
1939  if( state[m] != CONNECT )
1940  {
1941  SCIP_Real distnew;
1942
1943  assert(m != start);
1944
1945  if( stateprev && stateprev[m] == CONNECT )
1946  continue;
1947
1948  distnew = dist[k] + cost_csr[e];
1949
1950  assert(g->mark[m]);
1951
1952  if( distnew > distlimit )
1953  continue;
1954
1955  if( termmark[m] != 0 )
1956  distnew = MAX(distnew - g->prize[m], 0.0);
1957
1958  if( LT(distnew, dist[m]) )
1959  {
1960  if( prizeoffset && termmark[m] != 0 )
1961  {
1962  const SCIP_Real distnew0 = dist[k] + cost_csr[e];
1963
1964  if( g->prize[m] > distnew0 )
1965  {
1966  prizeoffset[m] = distnew0;
1967  assert(SCIPisZero(scip, distnew));
1968  }
1969  else
1970  prizeoffset[m] = g->prize[m];
1971  }
1972
1973  if( !visited[m] )
1974  {
1975  visitlist[(*nvisits)++] = m;
1976  visited[m] = TRUE;
1977  }
1978
1980  if( m == end )
1981  {
1982  nchecks = edgelimit + 1;
1983  success = TRUE;
1984  break;
1985  }
1986
1987  dist[m] = distnew;
1988  graph_heap_correct(m, distnew, dheap);
1989  }
1990  }
1991  nchecks++;
1992  }
1993  }
1994
1995  return success;
1996 }
1997
1998 /** modified Dijkstra along walks for PcMw, returns special distance between start and end */
2000  SCIP* scip, /**< SCIP data structure */
2001  const GRAPH* g, /**< graph data structure */
2002  const SCIP_Real* cost, /**< edge costs */
2003  SCIP_Real distlimit, /**< distance limit of the search */
2004  int start, /**< start */
2005  int end, /**< end */
2006  int edgelimit, /**< maximum number of edges to consider during execution */
2007  int maxnprevs, /**< maximum number of previous terminals to save */
2008  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
2009  int* prevterms, /**< previous terminals */
2010  int* nprevterms, /**< number of previous terminals */
2011  int* heap, /**< array representing a heap */
2012  int* state, /**< array to indicate whether a node has been scanned */
2013  int* visitlist, /**< stores all visited nodes */
2014  int* nvisits, /**< number of visited nodes */
2015  STP_Bool* visited /**< stores whether a node has been visited */
2016  )
2017 {
2018  int count;
2019  int nchecks;
2020  SCIP_Bool success = FALSE;
2021  const int edgelimit1 = edgelimit / 2;
2022
2023  assert(g != NULL);
2024  assert(heap != NULL);
2025  assert(dist != NULL);
2026  assert(cost != NULL);
2027  assert(visitlist != NULL);
2028  assert(nvisits != NULL);
2029  assert(visited != NULL);
2030  assert(graph_pc_isPcMw(g));
2031  assert(!g->extended);
2032
2033  *nvisits = 0;
2034
2036  return FALSE;
2037
2038  assert(g->mark[start] && g->mark[end]);
2039
2040  count = 0;
2041  nchecks = 0;
2042  dist[start] = 0.0;
2043  state[start] = CONNECT;
2044  visitlist[(*nvisits)++] = start;
2045
2046  g->mark[start] = FALSE;
2047  g->mark[end] = FALSE;
2048
2049  for( int e = g->outbeg[start]; e != EAT_LAST; e = g->oeat[e] )
2050  {
2051  const int m = g->head[e];
2052
2053  if( g->mark[m] && SCIPisLE(scip, cost[e], distlimit) )
2054  {
2055  assert(!visited[m]);
2056
2057  visitlist[(*nvisits)++] = m;
2058  visited[m] = TRUE;
2059  sdwalkUpdate(g, m, start, maxnprevs, prevterms, nprevterms);
2060
2061  if( Is_term(g->term[m]) )
2062  {
2063  const SCIP_Real newcost = MAX(cost[e] - g->prize[m], 0.0);
2064  sdwalkCorrectHeap(heap, state, &count, dist, m, newcost);
2065  }
2066  else
2067  {
2068  sdwalkCorrectHeap(heap, state, &count, dist, m, cost[e]);
2069  }
2070
2071  if( ++nchecks > edgelimit1 )
2072  break;
2073  }
2074  }
2075  assert(nprevterms[start] == 0);
2076
2077  g->mark[end] = TRUE;
2078
2079  while( count > 0 && nchecks <= edgelimit )
2080  {
2081  /* get nearest labelled node */
2082  const int k = nearestX(heap, state, &count, dist);
2083  assert(k != end && k != start);
2084  assert(SCIPisLE(scip, dist[k], distlimit));
2085
2086  state[k] = UNKNOWN;
2087
2088  /* correct incident nodes */
2089  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2090  {
2091  const int m = g->head[e];
2092
2093  if( g->mark[m] )
2094  {
2095  SCIP_Real distnew = dist[k] + cost[e];
2096
2097  assert(state[m] != CONNECT);
2098
2099  if( SCIPisGT(scip, distnew, distlimit) )
2100  continue;
2101
2102  if( Is_term(g->term[m]) )
2103  distnew = MAX(distnew - g->prize[m], 0.0);
2104
2105  if( distnew < dist[m] )
2106  {
2107  const SCIP_Bool mvisited = visited[m];
2108  if( !mvisited )
2109  {
2110  visitlist[(*nvisits)++] = m;
2111  visited[m] = TRUE;
2112  }
2113
2115  if( m == end )
2116  {
2117  nchecks = edgelimit + 1;
2118  success = TRUE;
2119  break;
2120  }
2121
2122  if( Is_term(g->term[m]) && sdwalkHasConflict(g, m, k, maxnprevs, prevterms, nprevterms, mvisited) )
2123  continue;
2124
2125  sdwalkUpdate(g, m, k, maxnprevs, prevterms, nprevterms);
2126  sdwalkCorrectHeap(heap, state, &count, dist, m, distnew);
2127  }
2128  }
2129  nchecks++;
2130  }
2131  }
2132
2133  g->mark[start] = TRUE;
2134  return success;
2135 }
2136
2137
2138
2139 /** modified Dijkstra along walks for PcMw, returns special distance between start and end */
2141  SCIP* scip, /**< SCIP data structure */
2142  const GRAPH* g, /**< graph data structure */
2143  const SCIP_Real* cost, /**< edge costs */
2144  const int* termmark, /**< terminal mark (2 for proper terminal, 1 for non-proper terminal, 0 otherwise) */
2145  SCIP_Real distlimit, /**< distance limit of the search */
2146  int start, /**< start */
2147  int end, /**< end */
2148  int edgelimit, /**< maximum number of edges to consider during execution */
2149  int maxnprevs, /**< maximum number of previous terminals to save */
2150  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
2151  int* prevterms, /**< previous terminals */
2152  int* nprevterms, /**< number of previous terminals */
2153  int* prevNPterms, /**< previous non-proper terminals */
2154  int* nprevNPterms, /**< number of previous non-proper terminals */
2155  int* prevedges, /**< previous edges */
2156  int* nprevedges, /**< number of previous edges */
2157  int* heap, /**< array representing a heap */
2158  int* state, /**< array to indicate whether a node has been scanned */
2159  int* visitlist, /**< stores all visited nodes */
2160  int* nvisits, /**< number of visited nodes */
2161  STP_Bool* visited /**< stores whether a node has been visited */
2162  )
2163 {
2164  int count;
2165  int nchecks;
2166  SCIP_Bool success = FALSE;
2167  const int edgelimit1 = edgelimit / 2;
2168
2169  assert(g != NULL);
2170  assert(heap != NULL);
2171  assert(dist != NULL);
2172  assert(cost != NULL);
2173  assert(visitlist != NULL);
2174  assert(nvisits != NULL);
2175  assert(visited != NULL);
2176  assert(graph_pc_isPcMw(g));
2177  assert(!g->extended);
2178
2179  *nvisits = 0;
2180
2182  return FALSE;
2183
2184  assert(g->mark[start] && g->mark[end]);
2185
2186  count = 0;
2187  nchecks = 0;
2188  dist[start] = 0.0;
2189  state[start] = CONNECT;
2190  visitlist[(*nvisits)++] = start;
2191
2192  g->mark[start] = FALSE;
2193  g->mark[end] = FALSE;
2194
2195  for( int e = g->outbeg[start]; e != EAT_LAST; e = g->oeat[e] )
2196  {
2197  const int m = g->head[e];
2198
2199  if( g->mark[m] && SCIPisLE(scip, cost[e], distlimit) )
2200  {
2201  SCIP_Real distnew = cost[e];
2202
2203  assert(!visited[m]);
2204
2205  visitlist[(*nvisits)++] = m;
2206  visited[m] = TRUE;
2207
2208  if( termmark[m] != 0 )
2209  distnew = MAX(distnew - g->prize[m], 0.0);
2210
2211  sdwalkUpdate2(termmark, m, start, e, maxnprevs, SCIPisZero(scip, distnew),
2212  prevterms, nprevterms, prevNPterms, nprevNPterms, prevedges, nprevedges);
2213  sdwalkCorrectHeap(heap, state, &count, dist, m, distnew);
2214
2215  if( ++nchecks > edgelimit1 )
2216  break;
2217  }
2218  }
2219  assert(nprevterms[start] == 0);
2220
2221  g->mark[end] = TRUE;
2222
2223  while( count > 0 && nchecks <= edgelimit )
2224  {
2225  /* get nearest labelled node */
2226  const int k = nearestX(heap, state, &count, dist);
2227  assert(k != end && k != start);
2228  assert(SCIPisLE(scip, dist[k], distlimit));
2229
2230  state[k] = UNKNOWN;
2231
2232  /* correct incident nodes */
2233  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2234  {
2235  const int m = g->head[e];
2236
2237  if( g->mark[m] )
2238  {
2239  SCIP_Real distnew = sdwalkGetdistnewEdge(prevedges, nprevedges, cost, dist, k, e, maxnprevs);
2240
2241  assert(state[m] != CONNECT);
2242
2243  if( SCIPisGT(scip, distnew, distlimit) )
2244  continue;
2245
2246  if( termmark[m] != 0 )
2247  distnew = sdwalkGetdistnewPrize(prevNPterms, nprevNPterms, termmark, visited, g->prize, k, m, distnew, maxnprevs);
2248
2249  if( SCIPisLT(scip, distnew, dist[m]) )
2250  {
2251  const SCIP_Bool mvisited = visited[m];
2252  if( !mvisited )
2253  {
2254  visitlist[(*nvisits)++] = m;
2255  visited[m] = TRUE;
2256  }
2257
2259  if( m == end )
2260  {
2261  nchecks = edgelimit + 1;
2262  success = TRUE;
2263  break;
2264  }
2265
2266  /* continue if m is proper terminals and is on the walk to k */
2267  if( termmark[m] == 2 && sdwalkHasConflict(g, m, k, maxnprevs, prevterms, nprevterms, mvisited) )
2268  continue;
2269
2270  sdwalkUpdate2(termmark, m, k, e, maxnprevs, SCIPisZero(scip, distnew),
2271  prevterms, nprevterms, prevNPterms, nprevNPterms, prevedges, nprevedges);
2272  sdwalkCorrectHeap(heap, state, &count, dist, m, distnew);
2273  }
2274  }
2275  nchecks++;
2276  }
2277  }
2278
2279  g->mark[start] = TRUE;
2280  return success;
2281 }
2282
2283
2284 /** modified Dijkstra along walks for PcMw */
2286  SCIP* scip, /**< SCIP data structure */
2287  const GRAPH* g, /**< graph data structure */
2288  const int* termmark, /**< terminal mark (2 for proper terminal, 1 for non-proper terminal, 0 otherwise) */
2289  const SCIP_Real* cost, /**< edge costs */
2290  const STP_Bool* endpoint, /**< stores whether search should be ended at vertex */
2291  int start, /**< start vertex */
2292  int edgelimit, /**< maximum number of edges to consider during execution */
2293  SCIP_Real* dist, /**< distances array, initially set to FARAWAY */
2294  int* visitlist, /**< stores all visited nodes */
2295  int* nvisits, /**< number of visited nodes */
2296  STP_Bool* visited, /**< stores whether a node has been visited */
2297  SCIP_Bool resetarrays /**< should arrays be reset? */
2298  )
2299 {
2300  int* const heap = g->path_heap;
2301  int* const state = g->path_state;
2302  int count;
2303  int nchecks;
2304  const SCIP_Real prize = g->prize[start];
2305
2306  assert(heap != NULL);
2307  assert(state != NULL);
2308  assert(dist != NULL);
2309  assert(cost != NULL);
2310  assert(visitlist != NULL);
2311  assert(nvisits != NULL);
2312  assert(visited != NULL);
2313  assert(graph_pc_isPcMw(g));
2314  assert(!g->extended);
2315  assert(Is_term(g->term[start]));
2317  assert(g->mark[start]);
2318
2319 #ifndef NDEBUG
2320  for( int k = 0; k < g->knots; k++ )
2321  {
2322  assert(state[k] == UNKNOWN);
2323  assert(visited[k] == FALSE);
2324  assert(dist[k] == FARAWAY);
2325  }
2326 #endif
2327
2328  *nvisits = 0;
2329  nchecks = 0;
2330  count = 1;
2331  heap[count] = start;
2332  state[start] = count;
2333  dist[start] = 0.0;
2334  visitlist[(*nvisits)++] = start;
2335  g->mark[start] = FALSE;
2336
2337  while( count > 0 && nchecks <= edgelimit )
2338  {
2339  /* get nearest labelled node */
2340  const int k = nearestX(heap, state, &count, dist);
2341  assert(SCIPisLE(scip, dist[k], prize));
2342
2343  if( termmark[k] == 2 )
2344  state[k] = CONNECT;
2345  else
2346  state[k] = UNKNOWN;
2347
2348  /* correct incident nodes */
2349  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2350  {
2351  const int m = g->head[e];
2352
2353  if( (state[m] != CONNECT) && g->mark[m] )
2354  {
2355  SCIP_Real distnew = dist[k] + cost[e];
2356
2357  if( SCIPisGT(scip, distnew, prize) )
2358  continue;
2359
2360  if( termmark[m] != 0 )
2361  distnew -= g->prize[m];
2362
2363  if( distnew < dist[m] )
2364  {
2365  if( !visited[m] )
2366  {
2367  visitlist[(*nvisits)++] = m;
2368  visited[m] = TRUE;
2369  }
2370
2372  if( endpoint != NULL && endpoint[m] )
2373  {
2374  g->mark[start] = TRUE;
2375  if( resetarrays )
2376  sdwalkReset(*nvisits, visitlist, dist, state, visited);
2377
2378  return TRUE;
2379  }
2380
2381  sdwalkCorrectHeap(heap, state, &count, dist, m, distnew);
2382  }
2383  }
2384  nchecks++;
2385  }
2386  }
2387
2388  g->mark[start] = TRUE;
2389
2390  if( resetarrays )
2391  sdwalkReset(*nvisits, visitlist, dist, state, visited);
2392
2393  return FALSE;
2394 }
2395
2396
2397 /** computes (or rather updates) SDs between all */
2399  SCIP* scip, /**< SCIP data structure */
2400  const GRAPH* g, /**< graph data structure */
2401  const SDPROFIT* sdprofit, /**< profit or NULL */
2402  SDCLIQUE* cliquedata /**< data */
2403 )
2404 {
2405  CLIQUEPATHS cliquepaths = { NULL, NULL, NULL, NULL, NULL, NULL };
2406  SCIP_Real* sds = cliquedata->sds;
2407
2408  assert(scip && g && cliquedata);
2409  assert(cliquedata->dijkdata && cliquedata->cliquenodes);
2410  assert(cliquedata->ncliquenodes >= 2);
2411  assert(g->knots >= cliquedata->ncliquenodes);
2412
2413  SCIP_CALL( cliquePathsInitData(scip, g, &cliquepaths) );
2414
2415  sdCliqueStarInit(g, cliquedata, &cliquepaths);
2416  sdCliqueStarComputeSds(g, sdprofit, cliquedata, &cliquepaths, sds);
2417
2418  cliquePathsFreeData(scip, &cliquepaths);
2419
2420  return SCIP_OKAY;
2421 }
SCIP_RETCODE graph_sdStarBiased(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, int star_root, int *star_base, DIJK *dijkdata, SCIP_Bool *success)
int * visitlist
Definition: graphdefs.h:314
Definition: graphdefs.h:224
SCIP_Bool graph_sdWalksExt(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real distlimit, int start, int end, int edgelimit, int maxnprevs, SCIP_Real *dist, int *prevterms, int *nprevterms, int *heap, int *state, int *visitlist, int *nvisits, STP_Bool *visited)
static void sdwalkReset(int nvisits, const int *visitlist, SCIP_Real *dist, int *state, STP_Bool *visited)
Definition: graph_sdpath.c:394
static SCIP_Real sdwalkGetdistnewPrize(const int *prevNPterms, const int *nprevNPterms, const int *termmark, const STP_Bool *visited, const SCIP_Real *prize, int k, int m, SCIP_Real distnew, int maxnprevs)
Definition: graph_sdpath.c:140
#define Is_term(a)
Definition: graphdefs.h:102
SCIP_Bool graph_sdWalksTriangle(SCIP *scip, const GRAPH *g, const int *termmark, const int *stateprev, SCIP_Real distlimit, int start, int end, int edgelimit, SCIP_Real *prizeoffset, SCIP_Real *dist, int *visitlist, int *nvisits, DHEAP *dheap, STP_Bool *visited)
SCIP_Bool graph_sdWalksExt2(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const int *termmark, SCIP_Real distlimit, int start, int end, int edgelimit, int maxnprevs, SCIP_Real *dist, int *prevterms, int *nprevterms, int *prevNPterms, int *nprevNPterms, int *prevedges, int *nprevedges, int *heap, int *state, int *visitlist, int *nvisits, STP_Bool *visited)
struct special_distance_clique_paths CLIQUEPATHS
#define FALSE
Definition: def.h:87
int *RESTRICT path_state
Definition: graphdefs.h:256
#define SDSTAR_BASE_UNSET
Definition: graphdefs.h:76
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
includes various files containing graph methods used for Steiner tree problems
static void sdCliqueStarUpdateSd(const SDPROFIT *sdprofit, const int *nodes_pred, const SCIP_Real *nodes_baseToMaxDist, const SCIP_Real *nodes_distFromMax, const SCIP_Real *nodes_maxpathprofit, int ncliquenodes, int centernode, int newnode, int prednode, SCIP_Real newdist, SCIP_Real edgecost, const SCIP_Real *nodes_dist, const int *nodes_base, const int *baseToClique, SCIP_Real *RESTRICT sds)
Definition: graph_sdpath.c:840
static SCIP_Real sdwalkGetdistnewEdge(const int *prevedges, const int *nprevedges, const SCIP_Real *cost, const SCIP_Real *dist, int k, int e, int maxnprevs)
Definition: graph_sdpath.c:92
SCIP_Real * cost
Definition: graphdefs.h:164
#define EAT_LAST
Definition: graphdefs.h:58
SCIP_Bool graph_heap_isClean(const DHEAP *)
Definition: graph_util.c:962
#define FARAWAY
Definition: graphdefs.h:89
static int nearestX(int *RESTRICT heap, int *RESTRICT state, int *RESTRICT count, const SCIP_Real *pathdist)
Definition: graphheaps.h:41
static SCIP_Real sdCliqueStarGetDistLimit(const SDCLIQUE *cliquedata, const SCIP_Real *sds)
Definition: graph_sdpath.c:991
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
static void sdCliqueStarInit(const GRAPH *g, SDCLIQUE *cliquedata, CLIQUEPATHS *cliquepaths)
Definition: graph_sdpath.c:929
static SCIP_Bool sdwalkHasConflict(const GRAPH *g, int node, int prednode, int maxnprevs, const int *prevterms, const int *nprevterms, const SCIP_Bool nodevisited)
Definition: graph_sdpath.c:191
SCIP_Real * node_distance
Definition: graphdefs.h:316
SCIP_Real miscstp_maxReal(const SCIP_Real *realarr, unsigned nreals)
Definition: misc_stp.c:142
SCIP_Real SCIPepsilon(SCIP *scip)
real eps
static void cliquePathsFreeData(SCIP *scip, CLIQUEPATHS *cliquepaths)
Definition: graph_sdpath.c:74
int *RESTRICT mark
Definition: graphdefs.h:198
static void sdwalkUpdate(const GRAPH *g, int node, int prednode, int maxnprevs, int *prevterms, int *nprevterms)
Definition: graph_sdpath.c:232
static void sdwalkCorrectHeap(int *RESTRICT heap, int *RESTRICT state, int *RESTRICT count, SCIP_Real *RESTRICT pathdist, int node, SCIP_Real newcost)
Definition: graph_sdpath.c:416
int * position
Definition: graphdefs.h:305
int *RESTRICT oeat
Definition: graphdefs.h:231
#define LE(a, b)
Definition: portab.h:83
#define GE(a, b)
Definition: portab.h:85
SCIP_Bool extended
Definition: graphdefs.h:267
SCIP_RETCODE graph_sdCloseNodesBiased(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, int sourcenode, DIJK *dijkdata)
DHEAP * dheap
Definition: graphdefs.h:315
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real * prize
Definition: graphdefs.h:210
Definition: graphdefs.h:201
#define NULL
Definition: lpi_spx1.cpp:155
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
#define LT(a, b)
Definition: portab.h:81
static int sdCliqueStarGetSdPosition(int ncliquenodes, int newnode, int prednode, const SCIP_Real *nodes_dist, const int *nodes_base, const int *baseToClique)
Definition: graph_sdpath.c:457
unsigned char STP_Bool
Definition: portab.h:34
static SCIP_Real reduce_sdprofitGetProfit(const SDPROFIT *sdprofit, int node, int nonsource1, int nonsource2)
Definition: reducedefs.h:178
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define GT(a, b)
Definition: portab.h:84
#define SCIP_Bool
Definition: def.h:84
int *RESTRICT path_heap
Definition: graphdefs.h:255
static void sdCliqueStarUpdateNodeMaxDist(int newnode, int prednode, const SCIP_Real *nodes_dist, SCIP_Real newprofit, SCIP_Real newedgecost, SCIP_Real *RESTRICT nodes_distBaseToMax, SCIP_Real *RESTRICT nodes_distFromMax, SCIP_Real *RESTRICT nodes_maxpathprofit)
Definition: graph_sdpath.c:492
int *RESTRICT tail
Definition: graphdefs.h:223
#define MAX(x, y)
Definition: tclique_def.h:83
static void sdCliqueStarComputeSds(const GRAPH *g, const SDPROFIT *sdprofit, SDCLIQUE *cliquedata, CLIQUEPATHS *cliquepaths, SCIP_Real *RESTRICT sds)
static SCIP_RETCODE cliquePathsInitData(SCIP *scip, const GRAPH *g, CLIQUEPATHS *cliquepaths)
Definition: graph_sdpath.c:51
int *RESTRICT term
Definition: graphdefs.h:196
#define CONNECT
Definition: graphdefs.h:87
Portable definitions.
SCIP_RETCODE graph_sdComputeCliqueStar(SCIP *scip, const GRAPH *g, const SDPROFIT *sdprofit, SDCLIQUE *cliquedata)
static void sdwalkUpdateCopy(int node, int prednode, int maxnprevs, int *prev, int *nprev)
Definition: graph_sdpath.c:275
static SCIP_Real sdCliqueStarGetNodeBias(const SDPROFIT *sdprofit, int node, int nextnode, int prevnode, SCIP_Real edgecost, SCIP_Real dist)
Definition: graph_sdpath.c:608
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static void sdCliqueStarUpdateNodeDist(int newnode, int prednode, SCIP_Real newdist, DHEAP *RESTRICT dheap, SCIP_Real *RESTRICT nodes_dist, int *RESTRICT nodes_base, int *RESTRICT nodes_pred)
Definition: graph_sdpath.c:530
SCIP_Real * cost
Definition: graphdefs.h:221
SCIP_Bool graph_sdWalks_csr(SCIP *scip, const GRAPH *g, const int *termmark, SCIP_Real distlimit, int start, int end, int edgelimit, SCIP_Real *dist, int *visitlist, int *nvisits, DHEAP *dheap, STP_Bool *visited)
SCIP_Bool graph_pc_isPcMw(const GRAPH *)
#define SCIP_Real
Definition: def.h:177
int graph_heap_deleteMinReturnNode(DHEAP *)
Definition: graph_util.c:1076
static SCIP_Real sdGet2ProfitsDist(SCIP_Real distBaseToBase, SCIP_Real distBaseToProfit1, SCIP_Real distBaseToProfit2, SCIP_Real profit1, SCIP_Real profit2)
Definition: graph_sdpath.c:554
static void sdCliqueStarGetFinalProfitData(const SDPROFIT *sdprofit, int centernode, int node, int neighbor, const SCIP_Real *nodes_dist, const int *nodes_pred, const SCIP_Real *nodes_distBaseToMax, const SCIP_Real *nodes_distFromMax, const SCIP_Real *nodes_maxpathprofit, SCIP_Real *distToMax, SCIP_Real *distFromMax, SCIP_Real *maxprofit_node, SCIP_Bool *nodeHasMaxProfit)
Definition: graph_sdpath.c:793
int *RESTRICT outbeg
Definition: graphdefs.h:204
static void sdwalkUpdate2(const int *termmark, int node, int prednode, int edge, int maxnprevs, SCIP_Bool clear, int *prevterms, int *nprevterms, int *prevNPterms, int *nprevNPterms, int *prevedges, int *nprevedges)
Definition: graph_sdpath.c:296
STP_Bool * node_visited
Definition: graphdefs.h:317
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define UNKNOWN
Definition: sepa_mcf.c:4095
SCIP_Bool graph_sdWalksConnected(SCIP *scip, const GRAPH *g, const int *termmark, const SCIP_Real *cost, const STP_Bool *endpoint, int start, int edgelimit, SCIP_Real *dist, int *visitlist, int *nvisits, STP_Bool *visited, SCIP_Bool resetarrays)
#define nnodes
Definition: gastrans.c:65
includes (inline) graph heap methods used for Steiner tree problems
DCSR * dcsr_storage
Definition: graphdefs.h:271
void graph_sdStar(SCIP *scip, const GRAPH *g, SCIP_Bool with_zero_edges, int star_root, int edgelimit, int *star_base, SCIP_Real *dist, int *visitlist, int *nvisits, DHEAP *dheap, STP_Bool *visited, SCIP_Bool *success)
includes various reduction methods for Steiner tree problems
void graph_heap_correct(int, SCIP_Real, DHEAP *)
Definition: graph_util.c:1166
static SCIP_Real sdGet1ProfitDist(SCIP_Real distBaseToBase, SCIP_Real distBaseToProfit, SCIP_Real profit)
Definition: graph_sdpath.c:585
SCIP_Bool graph_sdWalks(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, const int *termmark, SCIP_Real distlimit, int start, int end, int edgelimit, SCIP_Real *dist, int *heap, int *state, int *visitlist, int *nvisits, STP_Bool *visited)