Scippy

SCIP

Solving Constraint Integer Programs

dpterms_core.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 dpterms_core.c
17  * @brief Core of dynamic programming solver for Steiner tree (sub-) problems with small number of terminals
18  * @author Daniel Rehfeldt
19  *
20  * Contains core methods.
21  *
22  */
23 
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25 //#define SCIP_DEBUG
26 #include "scip/scipdefplugins.h"
27 #include "scip/rbtree.h"
28 #include "dpterms.h"
29 #include "dptermsinterns.h"
30 #include "stpbitset.h"
31 #include "stpvector.h"
32 #include "stpprioqueue.h"
33 
34 
35 /*
36  * Data structures
37  */
38 
39 /** helper data */
40 typedef struct trace_triplet
41 {
42  SCIP_Real bdist_local; /**< temporary tree bottleneck distance */
43  SCIP_Real bdist_global; /**< tree bottleneck distance */
44  int index; /**< index */
45 } TTRIPLET;
46 
47 
48 /** saves some data updated in every iteration */
50 {
52  STP_Vectype(int) stack; /**< general purpose stack */
53  STP_Vectype(TTRIPLET) tripletstack; /**< special stack */
54  STP_Vectype(SOLTRACE) sol_traces; /**< traces of current sub-solution */
55  STP_Bitset sol_termbits; /**< marks terminals of sub-solution */
56  STP_Vectype(SOLTRACE) valid_traces; /**< traces of valid extension */
57  STP_Bitset valid_bitset; /**< marks valid roots */
58  SCIP_Real* nodes_dist; /**< weight of sub-ST rooted at node */
59 #ifdef STP_DPTERM_USEDA
60  SCIP_Real* nodes_reddist; /**< reduced cost weight of sub-ST rooted at node */
61 #endif
62  SCIP_Real* nodes_ub; /**< upper bounds for rule-out */
63  int* nodes_previdx0; /**< predecessor NOTE: with shift! */
64  int* nodes_previdx1; /**< predecessor */
65  SCIP_Bool* nodes_isValidRoot; /**< is node a valid root? */
66  int nnodes; /**< number of nodes */
67  int sol_nterms; /**< popcount */
68  int extterm; /**< extension terminal */
69 } DPITER;
70 
71 
72 /*
73  * Local methods
74  */
75 
76 
77 /** prints separator nodes in SCIP_DEBUG mode */
78 static
80  SCIP_Real maxvaliddist, /**< maximum value distance */
81  const DPITER* dpiterator /**< iterator */
82 )
83 {
84  assert(dpiterator);
85 
86 #ifdef SCIP_DEBUG
87  {
88  const SCIP_Real* const nodes_ub = dpiterator->nodes_ub;
89  const SCIP_Real* const nodes_dist = dpiterator->nodes_dist;
90  const int nnodes = dpiterator->nnodes;
91 
92  SCIPdebugMessage("separator nodes: \n");
93 
94  for( int i = 0; i < nnodes; i++ )
95  {
96  if( GT(nodes_ub[i], 0.0) )
97  {
98  printf("%d (nodes_ub=%f) \n", i, nodes_ub[i]);
99 
100  }
101  else if( LE(nodes_dist[i], maxvaliddist) )
102  {
103  printf("%d (dist=%f) \n", i, nodes_dist[i]);
104  }
105  }
106  }
107 #endif
108 
109 }
110 
111 
112 
113 /** assembles final solution nodes */
114 static
115 STP_Vectype(int) getSolnodesFinal(
116  SCIP* scip, /**< SCIP data structure */
117  DPMISC* dpmisc, /**< misc */
118  DPITER* dpiterator /**< iterator */
119 )
120 {
121  STP_Vectype(int) solnodes = NULL;
122  STP_Vectype(int) stack = dpiterator->stack;
123  STP_Vectype(SOLTRACE) traces_all = dpmisc->global_traces;
124 
125  assert(dpmisc->opt_root);
126  StpVecClear(stack);
127 
128  SCIPdebugMessage("add solution node %d \n", traces_all[dpmisc->opt_root].root);
129  StpVecPushBack(scip, solnodes, traces_all[dpmisc->opt_root].root);
130 
131  StpVecPushBack(scip, stack, dpmisc->opt_root);
132 
133  if( dpmisc->opt_prev[0] != -1 )
134  StpVecPushBack(scip, stack, dpmisc->opt_prev[0]);
135 
136  if( dpmisc->opt_prev[1] != -1 )
137  StpVecPushBack(scip, stack, dpmisc->opt_prev[1]);
138 
139  while( StpVecGetSize(stack) > 0 )
140  {
141  const int i = stack[StpVecGetSize(stack) - 1];
142  const int prev0 = traces_all[i].prevs[0];
143 
144  StpVecPopBack(stack);
145 
146  if( prev0 != -1 )
147  {
148  StpVecPushBack(scip, stack, prev0);
149 
150  if( traces_all[i].prevs[1] != -1 )
151  {
152  StpVecPushBack(scip, stack, traces_all[i].prevs[1]);
153  }
154  else
155  {
156  SCIPdebugMessage("add solution node %d \n", traces_all[prev0].root);
157  StpVecPushBack(scip, solnodes, traces_all[prev0].root);
158  }
159  }
160  }
161 
162  dpiterator->stack = stack;
163 
164  return solnodes;
165 }
166 
167 
168 /** helper */
169 static inline
171  STP_Bitset sol_bitset, /**< bitset */
172  const int* nodes_termId, /**< ID */
173  int node /**< node to check */
174 )
175 {
176  return( nodes_termId[node] != -1 && !stpbitset_bitIsTrue(sol_bitset, nodes_termId[node]) );
177 }
178 
179 
180 /** helper */
181 static
183  const GRAPH* graph, /**< graph */
184  const DPSOLVER* dpsolver, /**< solver */
185  const DPITER* dpiterator /**< iterator */
186 )
187 {
188  const SCIP_Real* const nodes_ub = dpiterator->nodes_ub;
189  STP_Bitset sol_termbits = dpiterator->sol_termbits;
190  const int* const terminals = dpsolver->dpgraph->terminals;
191  const int nterms = graph->terms;
192 
193  for( int i = 0; i < nterms; i++ )
194  {
195  if( !stpbitset_bitIsTrue(sol_termbits, i)
196  && GT(nodes_ub[terminals[i]], 0.0) )
197  {
198  SCIPdebugMessage("terminal %d has positive UB, all extension are invalid! \n", terminals[i]);
199  return TRUE;
200  }
201  }
202 
203  return FALSE;
204 }
205 
206 
207 /** gets ordered root indices according to the solution costs */
208 static
210  SCIP* scip, /**< SCIP data structure */
211  DPITER* dpiterator, /**< iterator */
212  int* roots_indices /**< to initialize */
213 )
214 {
215  const int nnodes = dpiterator->nnodes;
216  SCIP_Real* roots_cost;
217  STP_Vectype(SOLTRACE) soltraces = dpiterator->sol_traces;
218  const int nsolcands = StpVecGetSize(dpiterator->sol_traces);
219 
220  SCIP_CALL( SCIPallocBufferArray(scip, &roots_cost, nnodes) );
221 
222  for( int i = 0; i < nsolcands; i++ )
223  {
224  roots_indices[i] = i;
225  roots_cost[i] = -soltraces[i].cost;
226  }
227 
228  SCIPsortDownRealInt(roots_cost, roots_indices, nsolcands);
229 
230  SCIPfreeBufferArray(scip, &roots_cost);
231 
232 #ifndef NDEBUG
233  for( int i = 1; i < nsolcands; i++ )
234  {
235  const int ind = roots_indices[i];
236  const int ind_prev = roots_indices[i - 1];
237  assert(LE(soltraces[ind_prev].cost, soltraces[ind].cost));
238  }
239 #endif
240 
241  return SCIP_OKAY;
242 }
243 
244 
245 /** initializes */
246 static
248  SCIP* scip, /**< SCIP data structure */
249  const GRAPH* graph, /**< original graph */
250  DPITER** dpiterator /**< to initialize */
251 )
252 {
253  DPITER* iter;
254  const int nnodes = graph_get_nNodes(graph);
255 
256  SCIP_CALL( SCIPallocMemory(scip, dpiterator) );
257  iter = *dpiterator;
258 
259  iter->stack = NULL;
260  iter->tripletstack = NULL;
261  iter->sol_traces = NULL;
262  iter->sol_termbits = NULL;
263  iter->valid_traces = NULL;
264  iter->valid_bitset = NULL;
265  iter->nnodes = nnodes;
266  iter->sol_nterms = -1;
267  iter->extterm = -1;
268 
269  /* NOTE: could be any positive number */
270  StpVecReserve(scip, iter->stack, 2);
271  StpVecReserve(scip, iter->valid_traces, 2);
272 
273 #ifdef STP_DPTERM_USEDA
274  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_reddist), nnodes) );
275 #endif
276  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_dist), nnodes) );
277  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_ub), nnodes) );
278  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_previdx0), nnodes) );
279  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_previdx1), nnodes) );
280  SCIP_CALL( SCIPallocMemoryArray(scip, &(iter->nodes_isValidRoot), nnodes) );
281 
282  return SCIP_OKAY;
283 }
284 
285 
286 /** frees */
287 static
289  SCIP* scip, /**< SCIP data structure */
290  DPITER** dpiterator /**< to free */
291 )
292 {
293  DPITER* iter = *dpiterator;
294 
295  assert(!iter->sol_traces && !iter->sol_termbits);
296  StpVecFree(scip, iter->valid_traces);
297  StpVecFree(scip, iter->tripletstack);
298  StpVecFree(scip, iter->stack);
299 
300  SCIPfreeMemoryArray(scip, &(iter->nodes_isValidRoot));
301  SCIPfreeMemoryArray(scip, &(iter->nodes_previdx1));
302  SCIPfreeMemoryArray(scip, &(iter->nodes_previdx0));
303  SCIPfreeMemoryArray(scip, &(iter->nodes_ub));
304  SCIPfreeMemoryArray(scip, &(iter->nodes_dist));
305 #ifdef STP_DPTERM_USEDA
306  SCIPfreeMemoryArray(scip, &(iter->nodes_reddist));
307 #endif
308 
309  SCIPfreeMemory(scip, dpiterator);
310 }
311 
312 
313 /** sets arrays to default values */
314 static
316  SCIP* scip, /**< SCIP data structure */
317  DPITER* dpiterator /**< to update */
318 )
319 {
320  SCIP_Real* RESTRICT nodes_dist = dpiterator->nodes_dist;
321  SCIP_Real* RESTRICT nodes_ub = dpiterator->nodes_ub;
322  int* RESTRICT nodes_pred1 = dpiterator->nodes_previdx0;
323  int* RESTRICT nodes_pred2 = dpiterator->nodes_previdx1;
324  SCIP_Bool* RESTRICT nodes_isValidRoot = dpiterator->nodes_isValidRoot;
325  const int nnodes = dpiterator->nnodes;
326 
327 #ifdef STP_DPTERM_USEDA
328  SCIP_Real* RESTRICT nodes_reddist = dpiterator->nodes_reddist;
329 
330  for( int i = 0; i < nnodes; i++ )
331  nodes_reddist[i] = 0.0;
332 #endif
333 
334  dpiterator->extterm = -1;
335 
336 
337  for( int i = 0; i < nnodes; i++ )
338  nodes_dist[i] = FARAWAY;
339 
340  for( int i = 0; i < nnodes; i++ )
341  nodes_ub[i] = FARAWAY;
342 
343  for( int i = 0; i < nnodes; i++ )
344  nodes_pred1[i] = -1;
345 
346  for( int i = 0; i < nnodes; i++ )
347  nodes_pred2[i] = -1;
348 
349  BMSclearMemoryArray(nodes_isValidRoot, nnodes);
350 }
351 
352 
353 /** gets solution from tree */
354 static
356  SCIP* scip, /**< SCIP data structure */
357  DPSOLVER* dpsolver, /**< solver */
358  DPITER* dpiterator /**< to update */
359 )
360 {
361  DPSUBSOL* subsol;
362 
363  stpprioqueue_deleteMin((void**) &(dpiterator->sol_termbits), &(dpiterator->sol_nterms), dpsolver->solpqueue);
364 
365  if( findSubsol(dpsolver->soltree_root, dpiterator->sol_termbits, &subsol) == 0 )
366  {
367  dpiterator->sol_traces = subsol->traces;
368  SCIPdebugMessage("number of traces: %d \n", StpVecGetSize(dpiterator->sol_traces));
369 
370  SCIPrbtreeDelete(&(dpsolver->soltree_root), subsol);
371 
372  assert(stpbitset_areEqual(dpiterator->sol_termbits, subsol->bitkey));
373  stpbitset_free(scip, &(subsol->bitkey));
374  }
375  else
376  {
377  assert(0 && "should never happen");
378  }
379 
380  assert(stpbitset_getPopcount(dpiterator->sol_termbits) == dpiterator->sol_nterms);
381  dpiterator->dpsubsol = subsol;
382 }
383 
384 
385 /** updates */
386 static
388  SCIP* scip, /**< SCIP data structure */
389  DPSOLVER* dpsolver, /**< solver */
390  DPITER* dpiterator /**< to update */
391 )
392 {
393  assert(!dpiterator->sol_termbits);
394  assert(!dpiterator->valid_bitset); /* NOTE: should be moved */
395 
396  dpiterSetDefault(scip, dpiterator);
397  dpiterPopSol(scip, dpsolver, dpiterator);
398 
399 #ifdef SCIP_DEBUG
400  SCIPdebugMessage("processing solution: \n");
401  stpbitset_print(dpiterator->sol_termbits);
402  SCIPdebugMessage("terminals: \n");
403  for( int i = 0; i < dpsolver->dpgraph->nterms; i++ )
404  {
405  if( stpbitset_bitIsTrue(dpiterator->sol_termbits, i) )
406  printf("%d ", dpsolver->dpgraph->terminals[i]);
407  }
408  printf(" \n");
409 #endif
410 }
411 
412 
413 /** (Implicitly) constructs sub-Steiner trees */
414 static
416  SCIP* scip, /**< SCIP data structure */
417  DPMISC* dpmisc, /**< misc */
418  DPITER* dpiterator /**< to update */
419 )
420 {
421  SCIP_Real* RESTRICT nodes_dist = dpiterator->nodes_dist;
422  SCIP_Real* RESTRICT nodes_ub = dpiterator->nodes_ub;
423  int* RESTRICT nodes_previdx0 = dpiterator->nodes_previdx0;
424  int* RESTRICT nodes_previdx1 = dpiterator->nodes_previdx1;
425  SCIP_Bool* RESTRICT nodes_isValidRoot = dpiterator->nodes_isValidRoot;
426  STP_Vectype(TTRIPLET) tripletstack = dpiterator->tripletstack;
427  STP_Vectype(SOLTRACE) soltraces = dpiterator->sol_traces;
428  STP_Vectype(SOLTRACE) global_traces = dpmisc->global_traces;
429  int* nodes_nvisits;
430  int* roots_indices;
431  const int nsolcands = StpVecGetSize(soltraces);
432  const int nnodes = dpiterator->nnodes;
433  int nvalidroots = 0;
434 
435  SCIP_CALL( SCIPallocClearBufferArray(scip, &nodes_nvisits, nnodes) );
436  SCIP_CALL( SCIPallocBufferArray(scip, &roots_indices, nnodes) );
437  SCIP_CALL( getOrderedRootIndices(scip, dpiterator, roots_indices) );
438 
439  SCIPdebugMessage("building Steiner trees from %d candidates \n", nsolcands);
440 
441  for( int i = 0; i < nsolcands; i++ )
442  {
443  const SOLTRACE sol_trace = soltraces[roots_indices[i]];
444  const int sol_root = sol_trace.root;
445 
446  // todo GE enough?
447  if( GT(sol_trace.cost, nodes_dist[sol_root]) )
448  continue;
449 
450  SCIPdebugMessage("building solution with root %d and prev0=%d prev1=%d \n", sol_root,
451  sol_trace.prevs[0], sol_trace.prevs[1]);
452 
453 #ifdef STP_DPTERM_USEDA
454  dpiterator->nodes_reddist[sol_root] = sol_trace.redcost;
455 #endif
456  nodes_isValidRoot[sol_root] = TRUE;
457  nodes_dist[sol_root] = sol_trace.cost;
458  nodes_ub[sol_root] = 0.0;
459  nodes_nvisits[sol_root]++;
460  nvalidroots++;
461  nodes_previdx0[sol_root] = sol_trace.prevs[0];
462  nodes_previdx1[sol_root] = sol_trace.prevs[1];
463 
464  assert(sol_trace.prevs[0] != -1 || sol_trace.prevs[1] == -1);
465 
466  if( sol_trace.prevs[0] == -1 )
467  continue;
468 
469  assert(StpVecGetSize(tripletstack) == 0);
470  StpVecPushBack(scip, tripletstack, ((TTRIPLET) {0.0, 0.0, sol_trace.prevs[0]}));
471  StpVecPushBack(scip, tripletstack, ((TTRIPLET) {0.0, 0.0, sol_trace.prevs[1]}));
472 
473  while( StpVecGetSize(tripletstack) > 0 )
474  {
475  TTRIPLET triplet = tripletstack[StpVecGetSize(tripletstack) - 1];
476  SOLTRACE parent_trace = global_traces[triplet.index];
477  const int previdx0 = parent_trace.prevs[0];
478 
479  StpVecPopBack(tripletstack);
480 
481  if( previdx0 != -1 )
482  {
483  const SCIP_Real bdist_local = triplet.bdist_local;
484  const SCIP_Real bdist_global = triplet.bdist_global;
485  const int previdx1 = parent_trace.prevs[1];
486 
487  assert(previdx0 >= 0);
488 
489  /* merged solution? */
490  if( previdx1 != -1 )
491  {
492  assert(previdx0 != -1);
493 
494  StpVecPushBack(scip, tripletstack, ((TTRIPLET) {0.0, bdist_global, previdx0}));
495  StpVecPushBack(scip, tripletstack, ((TTRIPLET) {0.0, bdist_global, previdx1}));
496 
497  SCIPdebugMessage("...merged solution from prev0=%d prev1=%d", previdx0, previdx1);
498  }
499  else
500  {
501  const int curr_root = global_traces[previdx0].root;
502  const SCIP_Real curr_bdist_local = bdist_local + parent_trace.cost - global_traces[previdx0].cost;
503  const SCIP_Real curr_bdist_global = MAX(bdist_global, curr_bdist_local);
504 
505  SCIPdebugMessage("at pred. solution with index=%d root=%d \n", previdx0, global_traces[previdx0].root);
506 
507  assert(GE(parent_trace.cost - global_traces[previdx0].cost, 0.0));
508 
509  if( LT(sol_trace.cost, nodes_dist[curr_root]) )
510  {
511  nodes_dist[curr_root] = sol_trace.cost;
512 #ifdef STP_DPTERM_USEDA
513  dpiterator->nodes_reddist[curr_root] = sol_trace.redcost;
514 #endif
515  }
516 
517  nodes_nvisits[curr_root]++;
518  nodes_ub[curr_root] = MIN(nodes_ub[curr_root], curr_bdist_global);
519 
520  StpVecPushBack(scip, tripletstack, ((TTRIPLET) {curr_bdist_local, curr_bdist_global, previdx0}));
521  }
522  }
523  }
524  }
525 
526  for( int i = 0; i < nnodes; i++ )
527  {
528  if( nodes_nvisits[i] < nvalidroots )
529  nodes_ub[i] = 0.0;
530  }
531 
532  SCIPfreeBufferArray(scip, &roots_indices);
533  SCIPfreeBufferArray(scip, &nodes_nvisits);
534  dpiterator->tripletstack = tripletstack;
535 
536  return SCIP_OKAY;
537 }
538 
539 
540 /** extends sub-Steiner trees */
541 static
543  SCIP* scip, /**< SCIP data structure */
544  const GRAPH* graph, /**< graph */
545  DPSOLVER* dpsolver, /**< solver */
546  DPITER* dpiterator /**< iterator */
547 )
548 {
549  const CSR* const csr = graph->csr_storage;
550  DHEAP* const dheap = dpsolver->dheap;
551  const SCIP_Real* const cost_csr = csr->cost;
552  const int* const head_csr = csr->head;
553  const int* const start_csr = csr->start;
554  int* RESTRICT nodes_previdx0 = dpiterator->nodes_previdx0;
555  int* RESTRICT nodes_previdx1 = dpiterator->nodes_previdx1;
556  SCIP_Real* RESTRICT nodes_dist = dpiterator->nodes_dist;
557  SCIP_Bool* RESTRICT nodes_isValidRoot = dpiterator->nodes_isValidRoot;
558  const int nnodes = dpiterator->nnodes;
559  STP_Bitset sol_termbits = dpiterator->sol_termbits;
560  const int* const nodes_termId = dpsolver->dpgraph->nodes_termId;
561  int* terms_adjCount;
562  const SCIP_Bool breakEarly = (graph->terms - dpiterator->sol_nterms > 1);
563  const int* const grad = graph->grad;
564 #ifdef STP_DPTERM_USEDA
565  SCIP_Real* RESTRICT nodes_reddist = dpiterator->nodes_reddist;
566  const SCIP_Real* const csr_redcost = dpsolver->dpredcosts->csr_redcosts;
567 #endif
568 
569  assert(nnodes == graph->knots);
570  assert(dpiterator->extterm == -1);
571  assert(dpiterator->sol_nterms > 0 && dpiterator->sol_nterms <= graph->terms);
572 
573  graph_heap_clean(TRUE, dheap);
574 
575  SCIP_CALL( SCIPallocClearBufferArray(scip, &terms_adjCount, graph->terms) );
576 
577  for( int i = 0; i < nnodes; i++ )
578  {
579  if( LT(nodes_dist[i], FARAWAY) )
580  graph_heap_correct(i, nodes_dist[i], dheap);
581  }
582 
583  SCIPdebugMessage("starting Dijkstra with %d roots \n", dheap->size);
584 
585  /* run Dijkstra */
586  while( dheap->size > 0 )
587  {
588  int k;
589  SCIP_Real k_dist;
590 
591  graph_heap_deleteMin(&k, &k_dist, dheap);
592  assert(EQ(nodes_dist[k], k_dist));
593 
594  SCIPdebugMessage("updated node %d: dist=%f \n", k, k_dist);
595 
596  if( nodeIsNonSolTerm(sol_termbits, nodes_termId, k) )
597  {
598  assert(dpiterator->extterm == -1);
599  dpiterator->extterm = k;
600  break;
601  }
602  else
603  {
604  const int k_start = start_csr[k];
605  const int k_end = start_csr[k + 1];
606  const SCIP_Bool k_isValid = nodes_isValidRoot[k];
607 
608  for( int e = k_start; e != k_end; e++ )
609  {
610  const int m = head_csr[e];
611  const SCIP_Real distnew = k_dist + cost_csr[e];
612 
613  if( LT(distnew, nodes_dist[m]) )
614  {
615  nodes_isValidRoot[m] = k_isValid;
616  nodes_dist[m] = distnew;
617  nodes_previdx0[m] = k;
618  nodes_previdx1[m] = -1;
619 #ifdef STP_DPTERM_USEDA
620  nodes_reddist[m] = nodes_reddist[k] + csr_redcost[e];
621 #endif
622 
623  graph_heap_correct(m, distnew, dheap);
624  }
625  else if( EQ(distnew, nodes_dist[m]) && !k_isValid )
626  {
627  nodes_isValidRoot[m] = FALSE;
628  }
629 
630  if( nodeIsNonSolTerm(sol_termbits, nodes_termId, m) && breakEarly )
631  {
632  assert(nodes_termId[m] >= 0);
633  terms_adjCount[nodes_termId[m]]++;
634 
635  /* all neighbors hit? */
636  if( terms_adjCount[nodes_termId[m]] == grad[m] )
637  {
638  graph_heap_clean(FALSE, dheap);
639  assert(dheap->size == 0);
640  assert(dpiterator->extterm == -1);
641  dpiterator->extterm = m;
642  break;
643  }
644  }
645  }
646  }
647  }
648 
649  SCIPdebugMessage("...ending extension with root %d \n", dpiterator->extterm);
650 
651  SCIPfreeBufferArray(scip, &terms_adjCount);
652 
653  return SCIP_OKAY;
654 }
655 
656 
657 /** helper */
658 static
660  SCIP* scip, /**< SCIP data structure */
661  DPMISC* dpmisc, /**< DP misc data structure */
662  DPITER* dpiterator, /**< iterator */
663  SCIP_Bool* hasExtension /**< extensions existing? */
664 )
665 {
666  STP_Vectype(SOLTRACE) valid_traces = dpiterator->valid_traces;
667  const SCIP_Bool* const nodes_isValidRoot = dpiterator->nodes_isValidRoot;
668  const SCIP_Real* const nodes_dist = dpiterator->nodes_dist;
669  const int* const nodes_previdx0 = dpiterator->nodes_previdx0;
670  const int* const nodes_previdx1 = dpiterator->nodes_previdx1;
671  const int nnodes = dpiterator->nnodes;
672  STP_Bitset valid_bits = stpbitset_new(scip, nnodes);
673  int* rootmap;
674  const int nall = dpmisc->global_size;
675 
676  StpVecClear(valid_traces);
677 
678  for( int i = 0; i < nnodes; i++ )
679  {
680  if( nodes_isValidRoot[i] )
681  {
682  SCIPdebugMessage("valid trace node %d \n", i);
683 
684  assert(GE(nodes_dist[i], 0.0) && LT(nodes_dist[i], FARAWAY));
685 
686 #ifdef STP_DPTERM_USEDA
687  StpVecPushBack(scip, valid_traces,
688  ((SOLTRACE) { {nodes_previdx0[i], nodes_previdx1[i]}, nodes_dist[i], dpiterator->nodes_reddist[i], i }) );
689 #else
690  StpVecPushBack(scip, valid_traces,
691  ((SOLTRACE) { {nodes_previdx0[i], nodes_previdx1[i]}, nodes_dist[i], i }) );
692 #endif
693 
694  stpbitset_setBitTrue(valid_bits, i);
695  }
696  }
697 
698  /* no valid extensions? */
699  if( 0 == StpVecGetSize(valid_traces) )
700  {
701  SCIPdebugMessage("no valid extensions! \n");
702  *hasExtension = FALSE;
703  stpbitset_free(scip, &valid_bits);
704  return SCIP_OKAY;
705  }
706 
707  *hasExtension = TRUE;
708 
709  SCIP_CALL( SCIPallocBufferArray(scip, &rootmap, nnodes) );
710 #ifndef NDEBUG
711  for( int i = 0; i < nnodes; i++ )
712  rootmap[i] = -1;
713 #endif
714  for( int i = 0; i < StpVecGetSize(valid_traces); i++ )
715  {
716  assert(valid_traces[i].root >= 0);
717  rootmap[valid_traces[i].root] = nall + i;
718  }
719 
720  for( int i = 0; i < StpVecGetSize(valid_traces); i++ )
721  {
722  if( valid_traces[i].prevs[0] != -1 && valid_traces[i].prevs[1] == -1 )
723  {
724  assert(rootmap[valid_traces[i].prevs[0]] >= 0);
725  valid_traces[i].prevs[0] = rootmap[valid_traces[i].prevs[0]];
726  }
727  }
728 
729  assert(!dpiterator->valid_bitset);
730  dpiterator->valid_bitset = valid_bits;
731  dpiterator->valid_traces = valid_traces;
732  SCIPfreeBuffer(scip, &rootmap);
733 
734  return SCIP_OKAY;
735 }
736 
737 
738 
739 /** helper */
740 static
741 STP_Vectype(SOLTRACE) combineTraces(
742  SCIP* scip, /**< SCIP data structure */
743  STP_Vectype(SOLTRACE) traces1, /**< ordered traces */
744  int prevoffset1,
745  SOLTRACE* traces2, /**< ordered traces */
746  int ntraces2, /**< number of traces */
747  int prevoffset2
748 )
749 {
750  STP_Vectype(SOLTRACE) combined = NULL;
751  int pos1 = 0;
752  int pos2 = 0;
753  const int ntraces1 = StpVecGetSize(traces1);
754  assert(ntraces1 > 0 && ntraces2 > 0);
755 
756  while( pos1 < ntraces1 && pos2 < ntraces2 )
757  {
758  assert(pos1 == ntraces1 - 1 || traces1[pos1].root < traces1[pos1 + 1].root);
759  assert(pos2 == ntraces2 - 1 || traces2[pos2].root < traces2[pos2 + 1].root);
760 
761  if( traces1[pos1].root == traces2[pos2].root )
762  {
763 #ifdef STP_DPTERM_USEDA
764  StpVecPushBack(scip, combined, ((SOLTRACE)
765  { {prevoffset1 + pos1, prevoffset2 + pos2},
766  traces1[pos1].cost + traces2[pos2].cost,
767  traces1[pos1].redcost + traces2[pos2].redcost,
768  traces1[pos1].root
769  } ) );
770 #else
771  StpVecPushBack(scip, combined, ((SOLTRACE)
772  { {prevoffset1 + pos1, prevoffset2 + pos2},
773  traces1[pos1].cost + traces2[pos2].cost,
774  traces1[pos1].root
775  } ) );
776 #endif
777 
778  pos1++;
779  pos2++;
780  }
781  else if( traces1[pos1].root < traces2[pos2].root )
782  {
783  pos1++;
784  }
785  else
786  {
787  assert(traces1[pos1].root > traces2[pos2].root);
788  pos2++;
789  }
790  }
791 
792  assert(combined);
793 
794 #ifndef NDEBUG
795  for( int i = 1; i < StpVecGetSize(combined); i++ )
796  assert(combined[i - 1].root <= combined[i].root);
797 #endif
798 
799  return combined;
800 }
801 
802 
803 /** updates existing subsol */
804 static
806  SCIP* scip, /**< SCIP data structure */
807  STP_Vectype(SOLTRACE) combined_traces, /**< for updating */
808  DPITER* dpiterator, /**< iterator */
809  DPSUBSOL* dpsubsol /**< to be updated */
810 )
811 {
812  STP_Vectype(SOLTRACE) subsol_traces = dpsubsol->traces;
813  STP_Vectype(SOLTRACE) updated_traces = NULL;
814  int pos1 = 0;
815  int pos2 = 0;
816  const int nsubsol_traces = StpVecGetSize(subsol_traces);
817  const int ncombined_traces = StpVecGetSize(combined_traces);
818  const int nnodes = dpiterator->nnodes;
819 
820  assert(subsol_traces && combined_traces);
821 
822  while( pos1 < nsubsol_traces || pos2 < ncombined_traces )
823  {
824  const int root1 = (pos1 == nsubsol_traces) ? nnodes : subsol_traces[pos1].root;
825  const int root2 = (pos2 == ncombined_traces) ? nnodes : combined_traces[pos2].root;
826  assert(root1 < nnodes || root2 < nnodes);
827 
828  if( root1 == root2 )
829  {
830  if( LE(subsol_traces[pos1].cost, combined_traces[pos2].cost) )
831  StpVecPushBack(scip, updated_traces, subsol_traces[pos1]);
832  else
833  StpVecPushBack(scip, updated_traces, combined_traces[pos2]);
834  pos1++;
835  pos2++;
836  }
837  else if( root1 < root2 )
838  {
839  StpVecPushBack(scip, updated_traces, subsol_traces[pos1]);
840  pos1++;
841  }
842  else
843  {
844  StpVecPushBack(scip, updated_traces, combined_traces[pos2]);
845  pos2++;
846  }
847 
848  }
849  StpVecFree(scip, subsol_traces);
850  dpsubsol->traces = updated_traces;
851 }
852 
853 
854 /** combines new sub-Steiner trees */
855 static
857  SCIP* scip, /**< SCIP data structure */
858  const GRAPH* graph, /**< graph */
859  int index, /**< index of intersecting */
860  DPSOLVER* dpsolver, /**< solver */
861  DPITER* dpiterator /**< iterator */
862 )
863 {
864  DPSUBSOL* subsol;
865  DPMISC* dpmisc = dpsolver->dpmisc;
866  STP_Vectype(int) offsets = dpmisc->global_starts;
867  STP_Bitset composite_termbits = dpmisc->global_termbits[index];
868  SOLTRACE* composite_traces = &(dpmisc->global_traces[offsets[index]]);
869  const int composite_ntraces = offsets[index + 1] - offsets[index];
870  STP_Bitset combined_termbits = stpbitset_newOr(scip, composite_termbits, dpiterator->sol_termbits);
871  int subsol_pos;
872 
873  STP_Vectype(SOLTRACE) combined_traces = combineTraces(scip,
874  dpiterator->valid_traces, dpmisc->global_size,
875  composite_traces, composite_ntraces, offsets[index]);
876 
877  if( (subsol_pos = findSubsol(dpsolver->soltree_root, combined_termbits, &subsol)) == 0 )
878  {
879  assert(subsol);
880 
881  subsolUpdate(scip, combined_traces, dpiterator, subsol);
882  StpVecFree(scip, combined_traces);
883  stpbitset_free(scip, &combined_termbits);
884  }
885  else
886  {
887  DPSUBSOL* newsol;
888  const int ncombinedterms = stpbitset_getPopcount(combined_termbits);
889  if( 2 * ncombinedterms <= graph->terms )
890  {
891  SCIP_CALL( stpprioqueue_insert(scip, stpbitset_newCopy(scip, combined_termbits),
892  ncombinedterms, dpsolver->solpqueue) );
893  }
894 
895  SCIP_CALL( dpterms_dpsubsolInit(scip, &newsol) );
896  newsol->bitkey = combined_termbits;
897  newsol->traces = combined_traces;
898 
899  SCIPrbtreeInsert(&(dpsolver->soltree_root), subsol, subsol_pos, newsol);
900  }
901 
902  return SCIP_OKAY;
903 }
904 
905 
906 
907 /** updates best (global) solution */
908 static
910  SCIP* scip, /**< SCIP data structure */
911  DPSOLVER* dpsolver, /**< solver */
912  DPITER* dpiterator /**< iterator */
913 )
914 {
915  DPSUBSOL* dpsubsol;
916  DPMISC* dpmisc = dpsolver->dpmisc;
917  STP_Bitset sol_termstoggled = stpbitset_newNot(scip, dpiterator->sol_termbits,
918  dpsolver->dpgraph->nterms);
919 
920  if( findSubsol(dpsolver->soltree_root, sol_termstoggled, &dpsubsol) == 0 )
921  {
922  STP_Vectype(SOLTRACE) valid_traces = dpiterator->valid_traces;
923  STP_Vectype(SOLTRACE) toggled_traces = dpsubsol->traces;
924  int pos1 = 0;
925  int pos2 = 0;
926  const int ntraces1 = StpVecGetSize(valid_traces);
927  const int ntraces2 = StpVecGetSize(toggled_traces);
928  assert(ntraces1 > 0 && ntraces2 > 0);
929  assert(valid_traces && toggled_traces);
930 
931 #ifdef SCIP_DEBUG
932  SCIPdebugMessage("found toggled terminal bits: \n");
933  stpbitset_print(dpiterator->sol_termbits);
934 #endif
935 
936  while( pos1 < ntraces1 && pos2 < ntraces2 )
937  {
938  assert(pos1 == ntraces1 - 1 || valid_traces[pos1].root < valid_traces[pos1 + 1].root);
939  assert(pos2 == ntraces2 - 1 || toggled_traces[pos2].root < toggled_traces[pos2 + 1].root);
940 
941  if( valid_traces[pos1].root == toggled_traces[pos2].root )
942  {
943  const SCIP_Real newcost = valid_traces[pos1].cost + toggled_traces[pos2].cost;
944  if( LT(newcost, dpmisc->opt_obj) )
945  {
946  SCIPdebugMessage("updating incumbent obj %f->%f \n", dpmisc->opt_obj, newcost);
947  dpmisc->opt_obj = newcost;
948  dpmisc->opt_root = dpmisc->global_size + pos1;
949  dpmisc->opt_prev[0] = toggled_traces[pos2].prevs[0];
950  dpmisc->opt_prev[1] = toggled_traces[pos2].prevs[1];
951  }
952  pos1++;
953  pos2++;
954  }
955  else if( valid_traces[pos1].root < toggled_traces[pos2].root )
956  {
957  pos1++;
958  }
959  else
960  {
961  assert(valid_traces[pos1].root > toggled_traces[pos2].root);
962  pos2++;
963  }
964  }
965  }
966  stpbitset_free(scip, &sol_termstoggled);
967 }
968 
969 
970 /** finalizes */
971 static
973  SCIP* scip, /**< SCIP data structure */
974  const GRAPH* graph, /**< graph */
975  DPSOLVER* dpsolver, /**< solver */
976  DPITER* dpiterator /**< iterator */
977 )
978 {
979  DPMISC* dpmisc = dpsolver->dpmisc;
980  const int nextensions = StpVecGetSize(dpiterator->valid_traces);
981 
982  assert(nextensions > 0);
983  assert(dpiterator->sol_termbits);
984  assert(dpiterator->valid_bitset);
985 
986  if( 3 * dpiterator->sol_nterms <= graph->terms )
987  {
988  const int nsubsets = StpVecGetSize(dpmisc->global_termbits);
989 
990  SCIP_CALL( dpterms_streeInsert(scip, stpbitset_newCopy(scip, dpiterator->sol_termbits),
991  dpiterator->valid_bitset, nsubsets, dpsolver->dpstree) );
992  }
993  else
994  {
995  stpbitset_free(scip, &(dpiterator->valid_bitset));
996  }
997 
998  for( int i = 0; i < nextensions; i++ )
999  {
1000  StpVecPushBack(scip, dpmisc->global_traces, dpiterator->valid_traces[i]);
1001  }
1002  StpVecPushBack(scip, dpmisc->global_termbits, dpiterator->sol_termbits);
1003  StpVecPushBack(scip, dpmisc->global_termbitscount, dpiterator->sol_nterms);
1004  dpmisc->global_size += nextensions;
1005  StpVecPushBack(scip, dpmisc->global_starts, dpmisc->global_size);
1006 
1007  dpiterator->sol_termbits = NULL;
1008  dpiterator->valid_bitset = NULL;
1009 
1010  return SCIP_OKAY;
1011 }
1012 
1013 /** adds new sub-Steiner trees */
1014 static
1016  SCIP* scip, /**< SCIP data structure */
1017  const GRAPH* graph, /**< graph */
1018  DPSOLVER* dpsolver, /**< solver */
1019  DPITER* dpiterator /**< iterator */
1020 )
1021 {
1022  DPMISC* dpmisc = dpsolver->dpmisc;
1023  STP_Vectype(int) intersections;
1024  const int nterms = graph->terms;
1025  SCIP_Bool hasExtensions;
1026 
1027  SCIP_CALL( dpiterAddNewPrepare(scip, dpmisc, dpiterator, &hasExtensions) );
1028 
1029  if( !hasExtensions )
1030  {
1031  return SCIP_OKAY;
1032  }
1033 
1034  intersections = dpterms_streeCollectIntersects(scip, dpiterator->sol_termbits,
1035  dpiterator->valid_bitset, dpsolver->dpstree);
1036 
1037  for( int i = 0; i < StpVecGetSize(intersections); i++ )
1038  {
1039  const int pos = intersections[i];
1040  assert(0 <= pos && pos < dpmisc->global_size);
1041  assert( stpbitset_getPopcount(dpmisc->global_termbits[pos]) == dpmisc->global_termbitscount[pos]);
1042 
1043  if( 2 * dpiterator->sol_nterms + dpmisc->global_termbitscount[pos] > nterms )
1044  continue;
1045 
1046 #ifdef SCIP_DEBUG
1047  SCIPdebugMessage("intersecting bitset: \n");
1048  stpbitset_print(dpmisc->global_termbits[pos]);
1049 #endif
1050 
1051  SCIP_CALL( combineWithIntersecting(scip, graph, pos, dpsolver, dpiterator) );
1052  }
1053 
1054  updateIncumbent(scip, dpsolver, dpiterator);
1055  SCIP_CALL( subtreesAddNewFinalize(scip, graph, dpsolver, dpiterator) );
1056 
1057  StpVecFree(scip, intersections);
1058  return SCIP_OKAY;
1059 }
1060 
1061 
1062 /** remove non-valid sub-Steiner trees */
1063 static
1065  const GRAPH* graph, /**< graph */
1066  DPSOLVER* dpsolver, /**< solver */
1067  DPITER* dpiterator /**< iterator */
1068 )
1069 {
1070  const CSR* const csr = graph->csr_storage;
1071  DHEAP* const dheap = dpsolver->dheap;
1072  const SCIP_Real* const cost_csr = csr->cost;
1073  const int* const head_csr = csr->head;
1074  const int* const start_csr = csr->start;
1075  SCIP_Real* RESTRICT nodes_ub = dpiterator->nodes_ub;
1076  const int nnodes = dpiterator->nnodes;
1077 
1078  assert(nnodes == graph->knots);
1079 
1080  graph_heap_clean(TRUE, dheap);
1081 
1082  for( int i = 0; i < nnodes; i++ )
1083  {
1084  if( GT(nodes_ub[i], 0.0) )
1085  graph_heap_correct(i, -nodes_ub[i], dheap);
1086  }
1087 
1088  /* propagate UBs */
1089  while( dheap->size > 0 )
1090  {
1091  int k;
1092  SCIP_Real k_ub;
1093 
1094  graph_heap_deleteMin(&k, &k_ub, dheap);
1095  k_ub *= -1.0;
1096  assert(EQ(nodes_ub[k], k_ub));
1097 
1098  {
1099  const int k_start = start_csr[k];
1100  const int k_end = start_csr[k + 1];
1101 
1102  for( int e = k_start; e != k_end; e++ )
1103  {
1104  const int m = head_csr[e];
1105  const SCIP_Real ubnew = k_ub - cost_csr[e];
1106 
1107  if( GT(ubnew, nodes_ub[m]) )
1108  {
1109  nodes_ub[m] = ubnew;
1110  graph_heap_correct(m, -ubnew, dheap);
1111  }
1112  }
1113  }
1114  }
1115 }
1116 
1117 
1118 /** remove non-valid sub-Steiner trees */
1119 static
1121  SCIP* scip, /**< SCIP data structure */
1122  const GRAPH* graph, /**< graph */
1123  DPSOLVER* dpsolver, /**< solver */
1124  DPITER* dpiterator /**< iterator */
1125 )
1126 {
1127  const CSR* const csr = graph->csr_storage;
1128  DHEAP* const dheap = dpsolver->dheap;
1129  const int* const head_csr = csr->head;
1130  const int* const start_csr = csr->start;
1131  const SCIP_Real* const nodes_dist = dpiterator->nodes_dist;
1132  const SCIP_Real* const nodes_ub = dpiterator->nodes_ub;
1133  SCIP_Bool* RESTRICT nodes_isValidRoot = dpiterator->nodes_isValidRoot;
1134  STP_Bitset sol_bitset = dpiterator->sol_termbits;
1135  const int* const nodes_termId = dpsolver->dpgraph->nodes_termId;
1136  const int nnodes = dpiterator->nnodes;
1137  SCIP_Real* nodes_mindist;
1138  SCIP_Real maxvaliddist = -1.0;
1139  STP_Vectype(int) stack = dpiterator->stack;
1140  const int extterm = dpiterator->extterm;
1141  int termcount = graph->terms - dpiterator->sol_nterms;
1142 #ifdef STP_DPTERM_USEDA
1143  SCIP_Real* RESTRICT nodes_reddist = dpiterator->nodes_reddist;
1144  const SCIP_Real cutoffbound = dpsolver->dpredcosts->cutoffbound;
1145 #endif
1146 
1147  assert(nnodes == graph->knots);
1148  assert(termcount >= 1);
1149  assert(graph_knot_isInRange(graph, extterm));
1150 
1151  if( allExtensionsAreInvalid(graph, dpsolver, dpiterator) )
1152  {
1153  BMSclearMemoryArray(nodes_isValidRoot, nnodes);
1154  return SCIP_OKAY;
1155  }
1156 
1157  graph_heap_clean(TRUE, dheap);
1158  SCIP_CALL( SCIPallocBufferArray(scip, &nodes_mindist, nnodes) );
1159  StpVecClear(stack);
1160 
1161  for( int i = 0; i < nnodes; i++ )
1162  nodes_mindist[i] = -1.0;
1163 
1164  nodes_mindist[extterm] = nodes_dist[extterm];
1165  graph_heap_correct(extterm, -nodes_mindist[extterm], dheap);
1166  SCIPdebugMessage("starting to exclude vertices with separators \n");
1167 
1168  /* propagate UBs */
1169  while( dheap->size > 0 )
1170  {
1171  int heapnode;
1172  SCIP_Real heapnode_dist;
1173 
1174  graph_heap_deleteMin(&heapnode, &heapnode_dist, dheap);
1175  heapnode_dist *= -1.0;
1176  assert(EQ(nodes_mindist[heapnode], heapnode_dist));
1177  assert(GE(heapnode_dist, 0.0));
1178  assert(StpVecGetSize(stack) == 0);
1179 
1180  StpVecPushBack(scip, stack, heapnode);
1181  SCIPdebugMessage("heapnode=%d dist=%f \n", heapnode, heapnode_dist);
1182 
1183  while( StpVecGetSize(stack) > 0 )
1184  {
1185  const int m = stack[StpVecGetSize(stack) - 1];
1186  const int m_start = start_csr[m];
1187  const int m_end = start_csr[m + 1];
1188 
1189  StpVecPopBack(stack);
1190 
1191  if( nodeIsNonSolTerm(sol_bitset, nodes_termId, m) )
1192  {
1193  if( --termcount == 0 )
1194  {
1195  maxvaliddist = heapnode_dist;
1196  graph_heap_clean(FALSE, dheap);
1197  assert(dheap->size == 0);
1198  break;
1199  }
1200  }
1201 
1202  for( int e = m_start; e != m_end; e++ )
1203  {
1204  SCIP_Real dist_new;
1205  const int q = head_csr[e];
1206 
1207  if( GT(nodes_ub[q], 0.0) )
1208  {
1209  continue;
1210  }
1211 
1212  dist_new = MIN(heapnode_dist, nodes_dist[q]);
1213  assert(GE(dist_new, 0.0));
1214 
1215  if( GT(dist_new, nodes_mindist[q]) )
1216  {
1217  if( LE(heapnode_dist, nodes_dist[q]) )
1218  {
1219  SCIPdebugMessage("update and push-back node %d (dist=%f) \n", q, nodes_dist[q]);
1220  StpVecPushBack(scip, stack, q);
1221  }
1222  else
1223  {
1224  assert(EQ(dist_new, nodes_dist[q]));
1225  graph_heap_correct(q, -dist_new, dheap);
1226  }
1227 
1228  nodes_mindist[q] = dist_new;
1229  }
1230  }
1231  }
1232  }
1233 
1234 #ifdef STP_DPTERM_USEDA
1235  maxvaliddist = MIN(maxvaliddist, dpsolver->dpredcosts->upperbound);
1236 #endif
1237 
1238  for( int i = 0; i < nnodes; i++ )
1239  {
1240  if( GT(nodes_dist[i], maxvaliddist) )
1241  nodes_isValidRoot[i] = FALSE;
1242 #ifdef STP_DPTERM_USEDA
1243  // todo also check dpsolver->dpredcosts->nodes_rootdist[i]
1244  else if( GT(nodes_reddist[i], cutoffbound) )
1245  {
1246  //printf("kill %f > %f \n", nodes_reddist[i], cutoffbound);
1247  nodes_isValidRoot[i] = FALSE;
1248  }
1249 #endif
1250  }
1251 
1252  debugPrintSeparator(maxvaliddist, dpiterator);
1253 
1254  SCIPdebugMessage("maxvaliddist=%f \n", maxvaliddist);
1255 
1256  dpiterator->stack = stack;
1257  SCIPfreeBufferArray(scip, &nodes_mindist);
1258 
1259  return SCIP_OKAY;
1260 }
1261 
1262 
1263 /** frees etc. */
1264 static
1266  SCIP* scip, /**< SCIP data structure */
1267  DPITER* dpiterator /**< to update */
1268 )
1269 {
1270  assert(dpiterator->sol_traces == dpiterator->dpsubsol->traces);
1271 
1272  if( dpiterator->sol_termbits )
1273  {
1274  stpbitset_free(scip, &(dpiterator->sol_termbits));
1275  }
1276 
1277  dpterms_dpsubsolFree(scip, &(dpiterator->dpsubsol));
1278  dpiterator->sol_traces = NULL;
1279 }
1280 
1281 
1282 /*
1283  * Interface methods
1284  */
1285 
1286 
1287 /** solves problem */
1289  SCIP* scip, /**< SCIP data structure */
1290  GRAPH* graph, /**< graph */
1291  DPSOLVER* dpsolver, /**< solver */
1292  SCIP_Bool* wasSolved /**< was problem solved to optimality? */
1293 )
1294 {
1295  DPITER* dpiterator;
1296  STP_PQ* const solpqueue = dpsolver->solpqueue;
1297  DPMISC* const dpmisc = dpsolver->dpmisc;
1298 
1299  assert(scip && graph && dpsolver && wasSolved);
1300  assert(solpqueue && dpsolver->dpgraph && dpsolver->soltree_root && dpsolver->dpstree && dpmisc);
1301  assert(!stpprioqueue_isClean(solpqueue));
1302 
1303  *wasSolved = TRUE;
1304 
1305  SCIP_CALL( dpiterInit(scip, graph, &dpiterator) );
1306 
1307  /* DP loop */
1308  while( !stpprioqueue_isClean(solpqueue) )
1309  {
1310  dpiterGetNextSol(scip, dpsolver, dpiterator);
1311 
1312  /* construct implicit sub-Steiner trees by setting traces */
1313  SCIP_CALL( subtreesBuild(scip, dpmisc, dpiterator) );
1314 
1315  SCIP_CALL( subtreesExtend(scip, graph, dpsolver, dpiterator) );
1316 
1317  if( graph->terms - dpiterator->sol_nterms > 1 )
1318  {
1319  propagateUBs(graph, dpsolver, dpiterator);
1320  SCIP_CALL( subtreesRemoveNonValids(scip, graph, dpsolver, dpiterator) );
1321  }
1322 
1323  /* add all valid extensions (including merges) */
1324  SCIP_CALL( subtreesAddNew(scip, graph, dpsolver, dpiterator) );
1325 
1326  dpiterFinalizeSol(scip, dpiterator);
1327 
1328  if( dpmisc->global_size % 16 == 0 && SCIPisStopped(scip) )
1329  {
1330  *wasSolved = FALSE;
1331  break;
1332  }
1333  }
1334 
1335  assert(!dpsolver->solnodes);
1336 
1337  if( *wasSolved )
1338  {
1339  SCIPdebugMessage("OBJ=%f \n", dpmisc->opt_obj);
1340  dpsolver->solnodes = getSolnodesFinal(scip, dpmisc, dpiterator);
1341  }
1342 
1343  dpiterFree(scip, &dpiterator);
1344 
1345  return SCIP_OKAY;
1346 }
static SCIP_Bool allExtensionsAreInvalid(const GRAPH *graph, const DPSOLVER *dpsolver, const DPITER *dpiterator)
Definition: dpterms_core.c:182
static SCIP_RETCODE dpiterInit(SCIP *scip, const GRAPH *graph, DPITER **dpiterator)
Definition: dpterms_core.c:247
static volatile int nterms
Definition: interrupt.c:38
static SCIP_RETCODE dpiterAddNewPrepare(SCIP *scip, DPMISC *dpmisc, DPITER *dpiterator, SCIP_Bool *hasExtension)
Definition: dpterms_core.c:659
#define StpVecGetSize(vec)
Definition: stpvector.h:139
#define StpVecClear(vec)
Definition: stpvector.h:134
static void dpiterPopSol(SCIP *scip, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:355
static void dpiterGetNextSol(SCIP *scip, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:387
int * head
Definition: graphdefs.h:141
#define SCIPallocClearBufferArray(scip, ptr, num)
Definition: scip_mem.h:117
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
int terms
Definition: graphdefs.h:192
struct dynamic_programming_iterator DPITER
static void debugPrintSeparator(SCIP_Real maxvaliddist, const DPITER *dpiterator)
Definition: dpterms_core.c:79
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
#define FALSE
Definition: def.h:87
SCIP_Real cost
SCIP_Real * cost
Definition: graphdefs.h:142
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_RETCODE dpterms_coreSolve(SCIP *scip, GRAPH *graph, DPSOLVER *dpsolver, SCIP_Bool *wasSolved)
void graph_heap_clean(SCIP_Bool, DHEAP *)
Definition: graph_util.c:938
Dynamic programming solver for Steiner tree (sub-) problems with small number of terminals.
static void stpbitset_print(STP_Bitset bitset)
Definition: stpbitset.h:452
#define StpVecPopBack(vec)
Definition: stpvector.h:182
static void updateIncumbent(SCIP *scip, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:909
void graph_heap_deleteMin(int *, SCIP_Real *, DHEAP *)
Definition: graph_util.c:1053
#define StpVecReserve(scip, vec, _size_)
Definition: stpvector.h:186
#define SCIPdebugMessage
Definition: pub_message.h:87
#define SCIPrbtreeDelete(root, node)
Definition: rbtree.h:60
#define FARAWAY
Definition: graphdefs.h:89
void SCIPsortDownRealInt(SCIP_Real *realarray, int *intarray, int len)
static void propagateUBs(const GRAPH *graph, DPSOLVER *dpsolver, DPITER *dpiterator)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
SCIP_RETCODE stpprioqueue_insert(SCIP *scip, void *data, int newkey, STP_PQ *prioqueue)
Definition: stpprioqueue.c:236
int * start
Definition: graphdefs.h:140
header only, simple implementation of an STL like vector
CSR * csr_storage
Definition: graphdefs.h:270
static void dpiterSetDefault(SCIP *scip, DPITER *dpiterator)
Definition: dpterms_core.c:315
SCIP_Real bdist_global
Definition: dpterms_core.c:43
void stpprioqueue_deleteMin(void **data, int *key, STP_PQ *prioqueue)
Definition: stpprioqueue.c:157
static SCIP_RETCODE subtreesAddNew(SCIP *scip, const GRAPH *graph, DPSOLVER *dpsolver, DPITER *dpiterator)
#define LE(a, b)
Definition: portab.h:83
STP_Bitset
#define GE(a, b)
Definition: portab.h:85
static SCIP_Bool stpbitset_areEqual(STP_Bitset bitset1, STP_Bitset bitset2)
Definition: stpbitset.h:200
static void dpiterFree(SCIP *scip, DPITER **dpiterator)
Definition: dpterms_core.c:288
priority queue with integer keys
static SCIP_RETCODE getOrderedRootIndices(SCIP *scip, DPITER *dpiterator, int *roots_indices)
Definition: dpterms_core.c:209
SCIP_Real bdist_local
Definition: dpterms_core.c:42
int *RESTRICT grad
Definition: graphdefs.h:201
#define NULL
Definition: lpi_spx1.cpp:155
#define EQ(a, b)
Definition: portab.h:79
static SCIP_RETCODE combineWithIntersecting(SCIP *scip, const GRAPH *graph, int index, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:856
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
static STP_Bitset stpbitset_newNot(SCIP *scip, STP_Bitset bitset1, int size)
Definition: stpbitset.h:390
header only, simple implementation of a bitset
#define LT(a, b)
Definition: portab.h:81
struct trace_triplet TTRIPLET
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define GT(a, b)
Definition: portab.h:84
#define SCIP_Bool
Definition: def.h:84
static STP_Bitset stpbitset_newCopy(SCIP *scip, STP_Bitset bitsetOrg)
Definition: stpbitset.h:267
static void dpterms_dpsubsolFree(SCIP *scip, DPSUBSOL **subsol)
static STP_Vectype(int)
Definition: dpterms_core.c:115
static STP_Bitset stpbitset_new(SCIP *scip, int maxnbits)
Definition: stpbitset.h:75
#define MAX(x, y)
Definition: tclique_def.h:83
static SCIP_RETCODE subtreesRemoveNonValids(SCIP *scip, const GRAPH *graph, DPSOLVER *dpsolver, DPITER *dpiterator)
static void subsolUpdate(SCIP *scip, STP_Vectype(SOLTRACE) combined_traces, DPITER *dpiterator, DPSUBSOL *dpsubsol)
Definition: dpterms_core.c:805
static void stpbitset_free(SCIP *scip, STP_Bitset *bitset)
Definition: stpbitset.h:100
static SCIP_Bool stpbitset_bitIsTrue(STP_Bitset bitset, int index)
Definition: stpbitset.h:223
Dynamic programming internals for Steiner tree (sub-) problems with small number of terminals...
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:69
static void dpiterFinalizeSol(SCIP *scip, DPITER *dpiterator)
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:125
SCIP_Bool stpprioqueue_isClean(const STP_PQ *prioqueue)
Definition: stpprioqueue.c:78
static STP_Bitset stpbitset_newOr(SCIP *scip, STP_Bitset bitset1, STP_Bitset bitset2)
Definition: stpbitset.h:348
#define SCIPrbtreeInsert(r, p, c, n)
Definition: rbtree.h:61
SCIP_RETCODE dpterms_streeInsert(SCIP *scip, STP_Bitset termsmark, STP_Bitset rootsmark, int64_t nsubsets, DPSTREE *dpstree)
Definition: dpterms_util.c:449
#define SCIP_Real
Definition: def.h:177
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:694
#define StpVecFree(scip, vec)
Definition: stpvector.h:153
static SCIP_RETCODE subtreesBuild(SCIP *scip, DPMISC *dpmisc, DPITER *dpiterator)
Definition: dpterms_core.c:415
static int stpbitset_getPopcount(STP_Bitset bitset)
Definition: stpbitset.h:247
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
static DPSUBSOL ** subsol
intrusive red black tree datastructure
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
#define nnodes
Definition: gastrans.c:65
#define StpVecPushBack(scip, vec, value)
Definition: stpvector.h:167
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:123
static void stpbitset_setBitTrue(STP_Bitset bitset, int index)
Definition: stpbitset.h:129
static SCIP_RETCODE subtreesExtend(SCIP *scip, const GRAPH *graph, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:542
static SCIP_Bool nodeIsNonSolTerm(STP_Bitset sol_bitset, const int *nodes_termId, int node)
Definition: dpterms_core.c:170
default SCIP plugins
void graph_heap_correct(int, SCIP_Real, DHEAP *)
Definition: graph_util.c:1166
static SCIP_RETCODE subtreesAddNewFinalize(SCIP *scip, const GRAPH *graph, DPSOLVER *dpsolver, DPITER *dpiterator)
Definition: dpterms_core.c:972