Scippy

SCIP

Solving Constraint Integer Programs

heur_tm.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-2019 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 heur_tm.c
17  * @brief Shortest paths based primal heuristics for Steiner problems
18  * @author Gerald Gamrath
19  * @author Thorsten Koch
20  * @author Daniel Rehfeldt
21  * @author Michael Winkler
22  *
23  * This file implements several shortest paths based primal heuristics for Steiner problems, see
24  * "SCIP-Jack - A solver for STP and variants with parallelization extensions" by
25  * Gamrath, Koch, Maher, Rehfeldt and Shinano
26  *
27  * A list of all interface methods can be found in heur_tm.h
28  *
29  */
30 
31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32 #include <assert.h>
33 #include <string.h>
34 #include "heur_tm.h"
35 #include "probdata_stp.h"
36 #include "portab.h"
37 #include "scip/misc.h"
38 #include <math.h>
39 #define HEUR_NAME "TM"
40 #define HEUR_DESC "shortest path based primal heuristics for Steiner trees"
41 #define HEUR_DISPCHAR '+'
42 #define HEUR_PRIORITY 10000000
43 #define HEUR_FREQ 1
44 #define HEUR_FREQOFS 0
45 #define HEUR_MAXDEPTH -1
46 #define HEUR_TIMING (SCIP_HEURTIMING_BEFORENODE | SCIP_HEURTIMING_DURINGLPLOOP | SCIP_HEURTIMING_AFTERLPLOOP | SCIP_HEURTIMING_AFTERNODE)
47 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
48 
49 #define DEFAULT_EVALRUNS 25 /**< number of runs */
50 #define DEFAULT_INITRUNS 100 /**< number of initial runs */
51 #define DEFAULT_LEAFRUNS 15 /**< number of runs at leafs */
52 #define DEFAULT_ROOTRUNS 50 /**< number of runs at the root */
53 #define DEFAULT_DURINGLPFREQ 5 /**< frequency during LP solving */
54 #define DEFAULT_TYPE 0 /**< heuristic to execute */
55 #define DEFAULT_RANDSEED 5 /**< seed for pseudo-random functions */
56 
57 #define AUTO 0
58 #define TM_SP 1
59 #define TM_VORONOI 2
60 #define TM_DIJKSTRA 3
61 
62 #ifdef WITH_UG
63 int getUgRank(void);
64 #endif
65 
66 /*
67  * Data structures
68  */
69 
70 
71 /** primal heuristic data */
72 struct SCIP_HeurData
73 {
74  SCIP_Longint nlpiterations; /**< number of total LP iterations*/
75  SCIP_Longint ncalls; /**< number of total calls (of TM) */
76  SCIP_Longint nexecs; /**< number of total executions (of TM) */
77  SCIP_Real hopfactor; /**< edge multiplication factor for hop constrained problems */
78  int stp_type; /**< problem type */
79  int evalruns; /**< number of runs */
80  int initruns; /**< number of initial runs */
81  int leafruns; /**< number of runs at leafs */
82  int rootruns; /**< number of runs at the root */
83  int duringlpfreq; /**< frequency during LP solving */
84  int type; /**< Heuristic type: 0 automatic, 1 TM_SP, 2 TM_VORONOI, 3 TM_DIJKSTRA */
85  int beststartnode; /**< start node of the so far best found solution */
86  unsigned int randseed; /**< seed value for random number generator */
87  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
88  unsigned int timing; /**< timing for timing mask */
89 };
90 
91 /*
92  * Static methods
93  */
94 
95 /** information method for a parameter change of random seed */
96 static
97 SCIP_DECL_PARAMCHGD(paramChgdRandomseed)
98 { /*lint --e{715}*/
99  SCIP_HEURDATA* heurdata;
100  int newrandseed;
101 
102  newrandseed = SCIPparamGetInt(param);
103 
104  heurdata = (SCIP_HEURDATA*)SCIPparamGetData(param);
105  assert(heurdata != NULL);
106 
107  heurdata->randseed = (unsigned int)newrandseed;
108 
109  return SCIP_OKAY;
110 }
111 
112 
113 /** compute starting points among marked (w.r.t. g->mark) vertices for constructive heuristics */
115  GRAPH* graph, /**< graph data structure */
116  int* starts, /**< starting points array */
117  int* runs /**< pointer to number of runs */
118  )
119 {
120  int r;
121  int k;
122  int l;
123  int root;
124  int nruns;
125  int nnodes;
126  int nterms;
127  int randval;
128 
129  assert(runs != NULL);
130  assert(graph != NULL);
131  assert(starts != NULL);
132 
133  nruns = *runs;
134  root = graph->source;
135  nnodes = graph->knots;
136  nterms = graph->terms;
137 
138  r = 0;
139  if( graph->mark[root] && nruns > 0 )
140  starts[r++] = root;
141 
142  randval = nnodes - nterms;
143 
144  /* use non-isolated terminals as starting points for TM heuristic */
145  for( k = 0; k < nnodes; k++ )
146  {
147  if( r >= nruns || r >= nterms )
148  break;
149 
150  l = (k + randval) % nnodes;
151  if( Is_term(graph->term[l]) && graph->mark[l] && l != root )
152  starts[r++] = l;
153  }
154 
155  /* fill empty slots randomly */
156  for( k = 0; k < nnodes && r < nruns; k++ )
157  {
158  l = (k + randval) % nnodes;
159  if( !Is_term(graph->term[l]) && graph->mark[l] )
160  starts[r++] = l;
161  }
162 
163  *runs = r;
164 }
165 
166 /* prune the (rooted) prize collecting Steiner tree in such a way that all leaves are terminals */
168  SCIP* scip, /**< SCIP data structure */
169  const GRAPH* g, /**< graph structure */
170  const SCIP_Real* cost, /**< edge costs */
171  int* result, /**< ST edges (need to be set to UNKNOWN) */
172  STP_Bool* connected /**< ST nodes */
173  )
174 {
175  PATH* mst;
176  int i;
177  int j;
178  int e1;
179  int e2;
180  int k1;
181  int k2;
182  int root;
183  int count;
184  int nnodes;
185  SCIP_Bool rmw = (g->stp_type == STP_RMWCSP);
186  SCIP_Bool rpc = (g->stp_type == STP_RPCSPG);
187  nnodes = g->knots;
188  root = g->source;
189 
190  if( rmw )
191  {
192  for( i = 0; i < nnodes; i++ )
193  {
194  if( connected[i] && g->mark[i] )
195  g->mark[i] = TRUE;
196  else
197  g->mark[i] = FALSE;
198  }
199  if( !g->mark[root] )
200  {
201  int nedges = g->edges;
202  for( i = 0; i < nedges; i++ )
203  result[i] = CONNECT;
204  return SCIP_OKAY;
205  }
206  }
207  else
208  {
209  /* compute the MST, exclude all terminals */
210  for( i = 0; i < nnodes; i++ )
211  {
212  if( connected[i] && !Is_term(g->term[i]) )
213  g->mark[i] = TRUE;
214  else
215  g->mark[i] = FALSE;
216  }
217  }
218 
219  if( rpc )
220  {
221  g->mark[root] = TRUE;
222  }
223  else if( !rmw )
224  {
225  int a;
226  for( a = g->outbeg[root]; a != EAT_LAST; a = g->oeat[a] )
227  {
228  i = g->head[a];
229  if( !Is_term(g->term[i]) && connected[i] )
230  break;
231  }
232 
233  /* trivial solution? */
234  if( a == EAT_LAST )
235  {
236  printf("trivial solution in pruning \n");
237  for( a = g->outbeg[g->source]; a != EAT_LAST; a = g->oeat[a] )
238  {
239  i = g->head[a];
240  if( Is_term(g->term[i]) )
241  {
242  assert(connected[i]);
243  result[a] = CONNECT;
244  }
245  }
246  return SCIP_OKAY;
247  }
248 
249  assert(g->mark[i]);
250  root = i;
251  }
252  assert(root >= 0);
253  assert(root < nnodes);
254 
255  SCIP_CALL( SCIPallocBufferArray(scip, &mst, nnodes) );
256  graph_path_exec(scip, g, MST_MODE, root, cost, mst);
257 
258  for( i = 0; i < nnodes; i++ )
259  {
260  if( g->mark[i] && (mst[i].edge != -1) )
261  {
262  assert(g->path_state[i] == CONNECT);
263  assert(g->head[mst[i].edge] == i);
264  assert(result[mst[i].edge] == -1);
265  result[mst[i].edge] = CONNECT;
266  }
267  }
268 
269  /* connect all terminals */
270  for( i = 0; i < nnodes; i++ )
271  {
272  if( Is_term(g->term[i]) && i != g->source )
273  {
274  if( rmw )
275  if( g->mark[i] )
276  continue;
277 
278  e1 = g->inpbeg[i];
279  assert(e1 >= 0);
280  e2 = g->ieat[e1];
281 
282  if( e2 == EAT_LAST )
283  {
284  result[e1] = CONNECT;
285  }
286  else
287  {
288  assert(e2 >= 0);
289 
290  assert(g->ieat[e2] == EAT_LAST);
291  k1 = g->tail[e1];
292  k2 = g->tail[e2];
293  assert(k1 == g->source || k2 == g->source);
294  if( k1 != g->source && g->path_state[k1] == CONNECT )
295  {
296  result[e1] = CONNECT;
297  }
298  else if( k2 != g->source && g->path_state[k2] == CONNECT )
299  {
300  result[e2] = CONNECT;
301  }
302  else if( k1 == g->source )
303  {
304  result[e1] = CONNECT;
305  }
306  else if( k2 == g->source )
307  {
308  result[e2] = CONNECT;
309  }
310  }
311  }
312  else if( i == root && !rpc && !rmw )
313  {
314  for( e1 = g->inpbeg[i]; e1 != EAT_LAST; e1 = g->ieat[e1] )
315  if( g->tail[e1] == g->source )
316  break;
317  assert(e1 != EAT_LAST);
318  result[e1] = CONNECT;
319  }
320  }
321 
322  /* prune */
323  do
324  {
325  count = 0;
326 
327  for( i = nnodes - 1; i >= 0; --i )
328  {
329  if( !g->mark[i] || g->path_state[i] != CONNECT || Is_term(g->term[i]) )
330  continue;
331 
332  for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
333  if( result[j] == CONNECT )
334  break;
335 
336  if( j == EAT_LAST )
337  {
338  /* there has to be exactly one incoming edge
339  */
340  for( j = g->inpbeg[i]; j != EAT_LAST; j = g->ieat[j] )
341  {
342  if( result[j] == CONNECT )
343  {
344  result[j] = -1;
345  g->mark[i] = FALSE;
346  connected[i] = FALSE;
347  count++;
348  break;
349  }
350  }
351  assert(j != EAT_LAST);
352  }
353  }
354  }
355  while( count > 0 );
356 
357 #ifndef NDEBUG
358  /* make sure there is no unconnected vertex */
359  for( i = 0; i < nnodes; i++ )
360  {
361  if( connected[i] && i != g->source )
362  {
363  for( j = g->inpbeg[i]; j != EAT_LAST; j = g->ieat[j] )
364  if( result[j] == CONNECT )
365  break;
366 
367  assert(j != EAT_LAST);
368 
369  }
370  }
371 #endif
372 
373  assert(graph_sol_valid(scip, g, result));
374  SCIPfreeBufferArray(scip, &mst);
375 
376  return SCIP_OKAY;
377 }
378 
379 
380 
381 /** build (rooted) prize collecting Steiner tree in such a way that all leaves are terminals */
383  SCIP* scip, /**< SCIP data structure */
384  const GRAPH* g, /**< graph structure */
385  PATH* mst, /**< path data structure array */
386  const SCIP_Real* cost, /**< edge costs */
387  SCIP_Real* objresult, /**< pointer to store objective value of result */
388  int* connected /**< CONNECT/UNKNOWN */
389  )
390 {
391  SCIP_Real obj;
392  int i;
393  int j;
394  int e1;
395  int e2;
396  int k1;
397  int k2;
398  int root;
399  int count;
400  int nnodes;
401  int orgroot;
402 
403  assert(g != NULL);
404  assert(mst != NULL);
405  assert(scip != NULL);
406  assert(cost != NULL);
407  assert(connected != NULL);
408 
409  obj = 0.0;
410  nnodes = g->knots;
411  orgroot = g->source;
412 
413  /* compute the MST, exclude all terminals */
414  for( i = nnodes - 1; i >= 0; --i )
415  {
416  if( connected[i] == CONNECT && !Is_term(g->term[i]) )
417  g->mark[i] = TRUE;
418  else
419  g->mark[i] = FALSE;
420  }
421 
422  if( g->stp_type == STP_RPCSPG )
423  {
424  root = orgroot;
425  g->mark[root] = TRUE;
426  }
427  else
428  {
429  int a;
430  for( a = g->outbeg[orgroot]; a != EAT_LAST; a = g->oeat[a] )
431  {
432  i = g->head[a];
433  if( !Is_term(g->term[i]) && connected[i] == CONNECT )
434  break;
435  }
436 
437  /* trivial solution? */
438  if( a == EAT_LAST )
439  {
440  for( i = 0; i < nnodes; i++ )
441  mst[i].edge = UNKNOWN;
442 
443  printf("trivial solution in buildPcMwTree \n");
444  for( a = g->outbeg[orgroot]; a != EAT_LAST; a = g->oeat[a] )
445  {
446  i = g->head[a];
447  if( Is_term(g->term[i]) )
448  {
449  obj += cost[a];
450  mst[i].edge = a;
451  }
452  }
453  (*objresult) = obj;
454  return SCIP_OKAY;
455  }
456 
457  assert(g->mark[i]);
458  root = i;
459  }
460  assert(root >= 0);
461  assert(root < nnodes);
462 
463  graph_path_exec(scip, g, MST_MODE, root, cost, mst);
464 
465  /* connect all terminals */
466  for( i = nnodes - 1; i >= 0; --i )
467  {
468  if( Is_term(g->term[i]) && i != orgroot )
469  {
470  e1 = g->inpbeg[i];
471  assert(e1 >= 0);
472  e2 = g->ieat[e1];
473 
474  if( e2 == EAT_LAST )
475  {
476  mst[i].edge = e1;
477  }
478  else
479  {
480  assert(e2 >= 0);
481 
482  assert(g->ieat[e2] == EAT_LAST);
483  k1 = g->tail[e1];
484  k2 = g->tail[e2];
485  assert(k1 == orgroot || k2 == orgroot);
486 
487  if( k1 != orgroot && g->path_state[k1] == CONNECT )
488  {
489  mst[i].edge = e1;
490  }
491  else if( k2 != orgroot && g->path_state[k2] == CONNECT )
492  {
493  mst[i].edge = e2;
494  }
495  else if( k1 == orgroot )
496  {
497  mst[i].edge = e1;
498  }
499  else if( k2 == orgroot )
500  {
501  mst[i].edge = e2;
502  }
503  }
504  }
505  else if( i == root && g->stp_type != STP_RPCSPG )
506  {
507  for( e1 = g->inpbeg[i]; e1 != EAT_LAST; e1 = g->ieat[e1] )
508  if( g->tail[e1] == orgroot )
509  break;
510  assert(e1 != EAT_LAST);
511  mst[i].edge = e1;
512  }
513  }
514 
515  /* prune */
516  do
517  {
518  count = 0;
519 
520  for( i = nnodes - 1; i >= 0; --i )
521  {
522  if( !g->mark[i] || g->path_state[i] != CONNECT || Is_term(g->term[i]) )
523  continue;
524 
525  for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
526  {
527  e1 = mst[g->head[j]].edge;
528  if( e1 == j )
529  break;
530  }
531 
532  if( j == EAT_LAST )
533  {
534  mst[i].edge = UNKNOWN;
535  g->mark[i] = FALSE;
536  connected[i] = UNKNOWN;
537  count++;
538  break;
539  }
540  }
541  }
542  while( count > 0 );
543 
544  for( i = nnodes - 1; i >= 0; --i )
545  if( mst[i].edge >= 0 )
546  obj += cost[mst[i].edge];
547 
548  (*objresult) = obj;
549 
550  return SCIP_OKAY;
551 }
552 
553 
554 /** prune a Steiner tree in such a way that all leaves are terminals */
556  SCIP* scip, /**< SCIP data structure */
557  const GRAPH* g, /**< graph structure */
558  const SCIP_Real* cost, /**< edge costs */
559  int layer, /**< layer (usually 0) */
560  int* result, /**< ST edges, which need to be set to UNKNOWN */
561  STP_Bool* connected /**< ST nodes */
562  )
563 {
564  PATH* mst;
565  int i;
566  int j;
567  int count;
568  int nnodes;
569 
570  assert(g != NULL);
571  assert(scip != NULL);
572  assert(cost != NULL);
573  assert(layer == 0);
574  assert(result != NULL);
575  assert(connected != NULL);
576 
577  j = 0;
578  nnodes = g->knots;
579 
580  SCIP_CALL( SCIPallocBufferArray(scip, &mst, nnodes) );
581 
582  /* compute the MST */
583  for( i = nnodes - 1; i >= 0; --i )
584  {
585  if( connected[i] )
586  j++;
587  g->mark[i] = connected[i];
588  }
589 
590  assert(g->source >= 0);
591  assert(g->source < nnodes);
592  assert(j >= g->terms);
593 
594  graph_path_exec(scip, g, MST_MODE, g->source, cost, mst);
595 
596  for( i = nnodes - 1; i >= 0; --i )
597  {
598  if( connected[i] && (mst[i].edge != -1) )
599  {
600  assert(g->head[mst[i].edge] == i);
601  assert(result[mst[i].edge] == UNKNOWN);
602 
603  result[mst[i].edge] = layer;
604  }
605  }
606 
607  /* prune */
608  do
609  {
610  SCIPdebug(fputc('C', stdout));
611  SCIPdebug(fflush(stdout));
612 
613  count = 0;
614 
615  for( i = nnodes - 1; i >= 0; --i )
616  {
617  if( !g->mark[i] )
618  continue;
619 
620  if( g->term[i] == layer )
621  continue;
622 
623  for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
624  if( result[j] == layer )
625  break;
626 
627  if( j == EAT_LAST )
628  {
629  /* there has to be exactly one incoming edge
630  */
631  for( j = g->inpbeg[i]; j != EAT_LAST; j = g->ieat[j] )
632  {
633  if( result[j] == layer )
634  {
635  result[j] = -1;
636  g->mark[i] = FALSE;
637  connected[i] = FALSE;
638  count++;
639  break;
640  }
641  }
642  }
643  }
644  }
645  while( count > 0 );
646 
647  SCIPfreeBufferArray(scip, &mst);
648 
649  return SCIP_OKAY;
650 }
651 
652 /** build Steiner tree in such a way that all leaves are terminals */
654  SCIP* scip, /**< SCIP data structure */
655  const GRAPH* g, /**< graph structure */
656  PATH* mst, /**< path data structure array */
657  const SCIP_Real* cost, /**< edge costs */
658  SCIP_Real* objresult, /**< pointer to store objective value of result */
659  int* connected /**< CONNECT/UNKNOWN */
660  )
661 {
662  SCIP_Real obj;
663  int i;
664  int j;
665  int e1;
666  int root;
667  int count;
668  int nnodes;
669 
670  assert(g != NULL);
671  assert(mst != NULL);
672  assert(scip != NULL);
673  assert(cost != NULL);
674  assert(connected != NULL);
675 
676  obj = 0.0;
677  root = g->source;
678  nnodes = g->knots;
679 
680  /* compute the MST */
681  for( i = nnodes - 1; i >= 0; --i )
682  g->mark[i] = (connected[i] == CONNECT);
683 
684  graph_path_exec(scip, g, MST_MODE, root, cost, mst);
685 
686  /* prune */
687  do
688  {
689  count = 0;
690 
691  for( i = nnodes - 1; i >= 0; --i )
692  {
693  if( !g->mark[i] || Is_term(g->term[i]) )
694  continue;
695 
696  for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
697  {
698  e1 = mst[g->head[j]].edge;
699  if( e1 == j )
700  break;
701  }
702 
703  if( j == EAT_LAST )
704  {
705  mst[i].edge = UNKNOWN;
706  g->mark[i] = FALSE;
707  connected[i] = UNKNOWN;
708  count++;
709  break;
710  }
711  }
712  }
713  while( count > 0 );
714 
715  for( i = nnodes - 1; i >= 0; --i )
716  {
717  if( mst[i].edge >= 0 )
718  obj += cost[mst[i].edge];
719  else if( Is_term(g->term[i]) && i != root )
720  {
721  obj = FARAWAY;
722  break;
723  }
724  }
725 
726  *objresult = obj;
727 
728  return SCIP_OKAY;
729 }
730 
731 /** prune a degree constrained Steiner tree in such a way that all leaves are terminals */
733  SCIP* scip, /**< SCIP data structure */
734  const GRAPH* g, /**< graph structure */
735  int* result, /**< ST edges */
736  STP_Bool* connected /**< ST nodes (to be set) */
737  )
738 {
739  int* queue;
740  int i;
741  int j;
742  int qsize;
743  int count;
744  int nnodes;
745  assert(scip != NULL);
746  assert(g != NULL);
747  assert(result != NULL);
748  assert(connected != NULL);
749 
750  nnodes = g->knots;
751 
752  /*
753  * DFS until all terminals are reached
754  */
755 
756  for( i = 0; i < nnodes; i++ )
757  connected[i] = FALSE;
758 
759  SCIP_CALL( SCIPallocBufferArray(scip, &queue, nnodes) );
760 
761  qsize = 0;
762  queue[qsize++] = g->source;
763  connected[g->source] = TRUE;
764 
765  while( qsize )
766  {
767  const int node = queue[--qsize];
768  for( int e = g->outbeg[node]; e != EAT_LAST; e = g->oeat[e] )
769  {
770  if( (result[e] == CONNECT || result[flipedge(e)] == CONNECT) && !(connected[g->head[e]]) )
771  {
772  i = g->head[e];
773  result[e] = CONNECT;
774  result[flipedge(e)] = UNKNOWN;
775  connected[i] = TRUE;
776  queue[qsize++] = i;
777  }
778  }
779  }
780 
781  SCIPfreeBufferArray(scip, &queue);
782 
783  /* prune */
784  do
785  {
786  SCIPdebug(fputc('C', stdout));
787  SCIPdebug(fflush(stdout));
788 
789  count = 0;
790 
791  for( i = 0; i < nnodes; i++ )
792  {
793  if( !connected[i] )
794  continue;
795 
796  if( Is_term(g->term[i]) )
797  continue;
798 
799  for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
800  if( result[j] == CONNECT )
801  break;
802 
803  if( j == EAT_LAST )
804  {
805  /* there has to be exactly one incoming edge
806  */
807  for( j = g->inpbeg[i]; j != EAT_LAST; j = g->ieat[j] )
808  {
809  if( result[j] == CONNECT )
810  {
811  result[j] = UNKNOWN;
812  connected[i] = FALSE;
813  count++;
814  break;
815  }
816  }
817  assert(j != EAT_LAST);
818  }
819  }
820  }
821  while( count > 0 );
822 
823  return SCIP_OKAY;
824 }
825 
826 /*
827  * local functions
828  */
829 
830 static
832  SCIP* scip, /**< SCIP data structure */
833  const GRAPH* g, /**< graph structure */
834  const SCIP_Real* cost, /**< edge costs for DHCSTP */
835  int* result, /**< ST edges */
836  STP_Bool* connected /**< ST nodes */
837  )
838 {
839  const int nedges = g->edges;
840 
841  if( g->stp_type != STP_DHCSTP )
842  for( int e = 0; e < nedges; e++ )
843  result[e] = UNKNOWN;
844 
845  if( g->stp_type == STP_MWCSP || g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP )
846  SCIP_CALL( SCIPStpHeurTMPrunePc(scip, g, g->cost, result, connected) );
847  else
848  SCIP_CALL( SCIPStpHeurTMPrune(scip, g, (g->stp_type != STP_DHCSTP) ? g->cost : cost, 0, result, connected) );
849 
850  return SCIP_OKAY;
851 }
852 
853 /** Dijkstra based shortest paths heuristic */
854 static
856  SCIP* scip, /**< SCIP data structure */
857  const GRAPH* g, /**< graph structure */
858  SCIP_Real* cost, /**< edge costs */
859  SCIP_Real* dijkdist, /**< distance array */
860  int* result, /**< solution array (on edges) */
861  int* dijkedge, /**< predecessor edge array */
862  int start, /**< start vertex*/
863  SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
864  STP_Bool* connected /**< array marking all solution vertices*/
865  )
866 {
867  int k;
868  int nnodes;
869 
870  assert(g != NULL);
871 
872  nnodes = g->knots;
873 
874  for( k = 0; k < nnodes; k++ )
875  g->mark[k] = (g->grad[k] > 0);
876 
877  graph_path_st(scip, g, cost, dijkdist, dijkedge, start, randnumgen, connected);
878 
879  SCIP_CALL(prune(scip, g, cost, result, connected));
880 
881  return SCIP_OKAY;
882 }
883 
884 /** Dijkstra based shortest paths heuristic for PCSTP and MWCSP */
885 static
887  SCIP* scip, /**< SCIP data structure */
888  const GRAPH* g, /**< graph structure */
889  SCIP_Real* cost, /**< edge costs */
890  SCIP_Real* dijkdist, /**< distance array */
891  int* result, /**< solution array (on edges) */
892  int* dijkedge, /**< predecessor edge array */
893  int start, /**< start vertex */
894  STP_Bool* connected /**< array marking all solution vertices*/
895  )
896 {
897  if( g->stp_type == STP_RMWCSP )
898  graph_path_st_rmw(scip, g, cost, dijkdist, dijkedge, start, connected);
899  else if( g->stp_type == STP_RPCSPG )
900  graph_path_st_rpc(scip, g, cost, dijkdist, dijkedge, start, connected);
901  else
902  graph_path_st_pcmw(scip, g, cost, dijkdist, dijkedge, start, connected);
903 
904  SCIP_CALL(prune(scip, g, cost, result, connected));
905 
906  return SCIP_OKAY;
907 }
908 
909 /** Dijkstra based shortest paths heuristic for PCSTP and MWCSP that computes tree spanning all positive
910  * vertex weights and subsequently prunes this tree */
911 static
913  SCIP* scip, /**< SCIP data structure */
914  const GRAPH* g, /**< graph structure */
915  const SCIP_Real* cost, /**< edge costs */
916  SCIP_Real* dijkdist, /**< distance array */
917  int* result, /**< solution array (on edges) */
918  int* dijkedge, /**< predecessor edge array */
919  int start, /**< start vertex */
920  STP_Bool* connected /**< array marking all solution vertices*/
921  )
922 {
923  graph_path_st_pcmw_full(scip, g, cost, dijkdist, dijkedge, start, connected);
924 
925  SCIP_CALL(prune(scip, g, cost, result, connected));
926 
927  return SCIP_OKAY;
928 }
929 
930 
931 /** shortest paths based heuristic */
932 static
934  SCIP* scip, /**< SCIP data structure */
935  const GRAPH* g, /**< graph structure */
936  SCIP_Real* cost, /**< edge costs */
937  SCIP_Real* costrev, /**< reversed edge costs */
938  SCIP_Real** pathdist, /**< distance array */
939  int start, /**< start vertex */
940  int* perm, /**< permutation array (on nodes) */
941  int* result, /**< solution array (on edges) */
942  int* cluster, /**< array used to store current vertices of each subtree during construction */
943  int** pathedge, /**< predecessor edge array */
944  STP_Bool* connected, /**< array marking all solution vertices */
945  SCIP_RANDNUMGEN* randnumgen /**< random number generator */
946  )
947 {
948  SCIP_Real min;
949  SCIP_Bool directed = (g->stp_type == STP_SAP || g->stp_type == STP_DHCSTP);
950  int k;
951  int e;
952  int i;
953  int j;
954  int l;
955  int z;
956  int old;
957  int root;
958  int csize;
959  int newval;
960  int nnodes;
961 
962  assert(scip != NULL);
963  assert(g != NULL);
964  assert(result != NULL);
965  assert(connected != NULL);
966  assert(cost != NULL);
967  assert(costrev != NULL);
968  assert(pathdist != NULL);
969  assert(pathedge != NULL);
970  assert(cluster != NULL);
971  assert(perm != NULL);
972 
973  root = g->source;
974  csize = 0;
975  nnodes = g->knots;
976 
977  assert(0 <= start && start < nnodes);
978 
979  SCIPdebugMessage("Heuristic: Start=%5d \n", start);
980 
981  cluster[csize++] = start;
982 
983  /* initialize arrays */
984  for( i = 0; i < nnodes; i++ )
985  {
986  g->mark[i] = (g->grad[i] > 0);
987  connected[i] = FALSE;
988  perm[i] = i;
989  }
990 
991  connected[start] = TRUE;
992 
993  SCIPrandomPermuteIntArray(randnumgen, perm, 0, nnodes - 1);
994 
995  assert(graph_valid(g));
996 
997  /* CONSTCOND */
998  for( ;; )
999  {
1000  /* Find a terminal with minimal distance to current ST */
1001  min = FARAWAY;
1002  old = -1;
1003  newval = -1;
1004 
1005  /* time limit exceeded? */
1006  if( SCIPisStopped(scip) )
1007  return SCIP_OKAY;
1008 
1009  /* find shortest path from current tree to unconnected terminal */
1010  for( l = nnodes - 1; l >= 0; --l )
1011  {
1012  i = perm[l];
1013  if( !Is_term(g->term[i]) || connected[i] || !g->mark[i] || (directed && !connected[root] && i != root) )
1014  continue;
1015 
1016  z = SCIPrandomGetInt(randnumgen, 0, nnodes - 1);
1017 
1018  for( k = csize - 1; k >= 0; k-- )
1019  {
1020  j = cluster[(k + z) % csize];
1021  assert(i != j);
1022  assert(connected[j]);
1023 
1024  if( SCIPisLT(scip, pathdist[i][j], min) )
1025  {
1026  min = pathdist[i][j];
1027  newval = i;
1028  old = j;
1029  }
1030  }
1031  }
1032 
1033  /* no new path found? */
1034  if( newval == -1 )
1035  break;
1036 
1037  /* mark the new path */
1038  assert(old > -1);
1039  assert(newval > -1);
1040  assert(pathdist[newval] != NULL);
1041  assert(pathdist[newval][old] < FARAWAY);
1042  assert(g->term[newval] == 0);
1043  assert(!connected[newval]);
1044  assert(connected[old]);
1045 
1046  SCIPdebug(fputc('R', stdout));
1047  SCIPdebug(fflush(stdout));
1048 
1049  /* start from current tree */
1050  k = old;
1051 
1052  while( k != newval )
1053  {
1054  e = pathedge[newval][k];
1055  k = g->tail[e];
1056  if (!connected[k])
1057  {
1058  connected[k] = TRUE;
1059  cluster[csize++] = k;
1060  }
1061  }
1062  }
1063 
1064  /* prune the tree */
1065  SCIP_CALL( prune(scip, g, cost, result, connected) );
1066 
1067  return SCIP_OKAY;
1068 }
1069 
1070 
1071 /** heuristic for degree constrained STPs */
1072 static
1074  SCIP* scip, /**< SCIP data structure */
1075  const GRAPH* g, /**< graph data structure */
1076  SCIP_Real* cost, /**< edge costs */
1077  SCIP_Real* costrev, /**< reversed edge costs */
1078  SCIP_Real** pathdist, /**< distances from each terminal to all nodes */
1079  int start, /**< start vertex */
1080  int* perm, /**< permutation array */
1081  int* result, /**< array to indicate whether an edge is in the solution */
1082  int* cluster, /**< array for internal computations */
1083  int** pathedge, /**< ancestor edges for shortest paths from each terminal to all nodes */
1084  SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
1085  STP_Bool* connected, /**< array to indicate whether a vertex is in the solution */
1086  STP_Bool* solfound /**< pointer to store whether a solution has been found */
1087  )
1088 {
1089  SCIP_Real min;
1090  int csize = 0;
1091  int k;
1092  int e;
1093  int l;
1094  int i;
1095  int j;
1096  int t;
1097  int u;
1098  int z;
1099  int n;
1100  int old;
1101  int tldegcount;
1102  int degcount;
1103  int mindegsum;
1104  int degmax;
1105  int newval;
1106  int nnodes;
1107  int nterms;
1108  int termcount;
1109  int* degs;
1110  int* maxdegs;
1111  assert(scip != NULL);
1112  assert(g != NULL);
1113  assert(g->maxdeg != NULL);
1114  assert(result != NULL);
1115  assert(connected != NULL);
1116  assert(cost != NULL);
1117  assert(costrev != NULL);
1118  assert(pathdist != NULL);
1119  assert(pathedge != NULL);
1120  assert(cluster != NULL);
1121  assert(perm != NULL);
1122 
1123  z = 0;
1124  nterms = g->terms;
1125  nnodes = g->knots;
1126  mindegsum = 2;
1127  maxdegs = g->maxdeg;
1128 
1129  SCIPdebugMessage("Heuristic: Start=%5d ", start);
1130 
1131  SCIP_CALL( SCIPallocBufferArray(scip, &degs, nnodes) );
1132 
1133  cluster[csize++] = start;
1134 
1135  for( i = 0; i < nnodes; i++ )
1136  {
1137  degs[i] = 0;
1138  g->mark[i] = (g->grad[i] > 0);
1139  connected[i] = FALSE;
1140  perm[i] = i;
1141  }
1142 
1143  for( e = 0; e < g->edges; e++ )
1144  {
1145  assert(SCIPisGT(scip, cost[e], 0.0));
1146  assert(result[e] == UNKNOWN);
1147  }
1148  connected[start] = TRUE;
1149  tldegcount = MIN(g->grad[start], maxdegs[start]);
1150  SCIPrandomPermuteIntArray(randnumgen, perm, 0, nnodes - 1);
1151 
1152  if( Is_term(g->term[start]) )
1153  termcount = 1;
1154  else
1155  termcount = 0;
1156 
1157  for( n = 0; n < nnodes; n++ )
1158  {
1159  /* Find a terminal with minimal distance to the current ST */
1160  min = FARAWAY - 1;
1161  /* is the free degree sum at most one and are we not connecting the last terminal? */
1162  if( tldegcount <= 1 && termcount < nterms - 1)
1163  degmax = 1;
1164  else
1165  degmax = 0;
1166  old = -1;
1167  newval = -1;
1168  for( t = 0; t < nnodes; t++ )
1169  {
1170  i = perm[t];
1171  if( !Is_term(g->term[i]) || connected[i] || g->grad[i] == 0 )
1172  continue;
1173 
1174  z = SCIPrandomGetInt(randnumgen, 0, nnodes - 1);
1175 
1176  for( k = 0; k < csize; k++ )
1177  {
1178  j = cluster[(k + z) % csize];
1179  assert(i != j);
1180  assert(connected[j]);
1181 
1182  if( SCIPisLE(scip, pathdist[i][j], min) && degs[j] < maxdegs[j])
1183  {
1184  u = j;
1185  degcount = 1;
1186  while( u != i )
1187  {
1188  u = g->tail[pathedge[i][u]];
1189  if( !connected[u] )
1190  {
1191  if( (MIN(g->grad[u], maxdegs[u]) < 2 || Is_term(g->term[u])) && u != i )
1192  {
1193  degcount = -2;
1194  break;
1195  }
1196  degcount += MIN(g->grad[u] - 2, maxdegs[u] - 2);
1197  }
1198  else
1199  {
1200  assert(u != i);
1201  l = g->tail[pathedge[i][u]];
1202  if( !connected[l] && degs[u] >= maxdegs[u] )
1203  {
1204  degcount = -2;
1205  break;
1206  }
1207  }
1208  }
1209  if( degcount >= degmax || (degcount >= mindegsum && SCIPisLT(scip, pathdist[i][j], min)) )
1210  {
1211  degmax = degcount;
1212  min = pathdist[i][j];
1213  newval = i;
1214  old = j;
1215  }
1216  }
1217  }
1218  }
1219 
1220  if( newval == -1 )
1221  {
1222  j = UNKNOWN;
1223  for( k = 0; k < csize; k++ )
1224  {
1225  j = cluster[(k + z) % csize];
1226  if( degs[j] < maxdegs[j] )
1227  break;
1228  }
1229  if( j != UNKNOWN )
1230  {
1231  assert(k != csize);
1232 
1233  min = FARAWAY + 1;
1234  newval = UNKNOWN;
1235  for( e = g->outbeg[j]; e != EAT_LAST; e = g->oeat[e] )
1236  {
1237  u = g->head[e];
1238  if( !Is_term(g->term[u]) && !connected[u] && SCIPisGE(scip, min, cost[e]) )
1239  {
1240  min = cost[e];
1241  k = e;
1242  newval = u;
1243  }
1244  }
1245  if( newval != UNKNOWN )
1246  {
1247  result[flipedge(k)] = CONNECT;
1248  degs[newval]++;
1249  degs[j]++;
1250  connected[newval] = TRUE;
1251  cluster[csize++] = newval;
1252  tldegcount += MIN(maxdegs[newval], g->grad[newval]) - 2;
1253  continue;
1254  }
1255 
1256  }
1257  }
1258  tldegcount += degmax - 1;
1259  /* break? */
1260  if( newval == -1 )
1261  break;
1262  /* Weg setzten
1263  */
1264  assert(old > -1);
1265  assert(newval > -1);
1266  assert(pathdist[newval] != NULL);
1267  assert(g->term[newval] == 0);
1268  assert(!connected[newval]);
1269  assert(connected[old]);
1270 
1271  SCIPdebug(fputc('R', stdout));
1272  SCIPdebug(fflush(stdout));
1273 
1274  /* mark new tree nodes/edges */
1275  k = old;
1276  if( Is_term(g->term[newval]) )
1277  termcount++;
1278 
1279  while(k != newval)
1280  {
1281  e = pathedge[newval][k];
1282  u = k;
1283  k = g->tail[e];
1284 
1285  if( !connected[k])
1286  {
1287  result[flipedge(e)] = CONNECT;
1288  degs[u]++;
1289  degs[k]++;
1290  connected[k] = TRUE;
1291  cluster[csize++] = k;
1292  if( k != newval )
1293  assert(!Is_term(g->term[k]));
1294  }
1295  }
1296  if( termcount == nterms )
1297  break;
1298  assert(degs[newval] == 1);
1299  }
1300 
1301  *solfound = TRUE;
1302 
1303  for( i = 0; i < nnodes; i++ )
1304  {
1305  if( Is_term(g->term[i]) && !connected[i] )
1306  {
1307  *solfound = FALSE;
1308  break;
1309  }
1310  }
1311 
1312  if( *solfound )
1313  {
1314  /* prune the solution */
1315  SCIP_CALL( SCIPStpHeurTMBuildTreeDc(scip, g, result, connected) );
1316 
1317  for( t = 0; t < nnodes; t++ )
1318  if( degs[t] > maxdegs[t] )
1319  *solfound = FALSE;
1320  }
1321 
1322  SCIPfreeBufferArray(scip, &degs);
1323  return SCIP_OKAY;
1324 }
1325 
1326 /** Voronoi based shortest path heuristic */
1327 static
1329  SCIP* scip, /**< SCIP data structure */
1330  const GRAPH* g, /**< graph structure */
1331  SCIP_PQUEUE* pqueue, /**< priority queue */
1332  GNODE** gnodearr, /**< internal array */
1333  SCIP_Real* cost, /**< edge costs */
1334  SCIP_Real* costrev, /**< reversed edge costs */
1335  SCIP_Real** node_dist, /**< internal array */
1336  int start, /**< start vertex */
1337  int* result, /**< array to indicate whether an edge is in the solution */
1338  int* vcount, /**< internal array */
1339  int* nodenterms, /**< internal array */
1340  int** node_base, /**< internal array */
1341  int** node_edge, /**< internal array */
1342  STP_Bool firstrun, /**< method called for the first time? (during one heuristic round) */
1343  STP_Bool* connected /**< array to indicate whether a vertex is in the solution */
1344  )
1345 {
1346  int k;
1347  int i;
1348  int j;
1349  int best;
1350  int term;
1351  int count;
1352  int nnodes;
1353  int nterms;
1354 
1355  assert(scip != NULL);
1356  assert(g != NULL);
1357  assert(result != NULL);
1358  assert(connected != NULL);
1359  assert(cost != NULL);
1360  assert(costrev != NULL);
1361  nnodes = g->knots;
1362  nterms = g->terms;
1363 
1364  SCIPdebugMessage("TM_Polzin Heuristic: Start=%5d ", start);
1365 
1366  /* if the heuristic is called for the first time several data structures have to be set up */
1367  if( firstrun )
1368  {
1369  PATH* vnoi;
1370  SCIP_Real* vcost;
1371  int old;
1372  int oedge;
1373  int root = g->source;
1374  int ntovisit;
1375  int nneighbnodes;
1376  int nneighbterms;
1377  int nreachednodes;
1378  int* state;
1379  int* vbase;
1380  int* terms;
1381  int* tovisit;
1382  int* reachednodes;
1383  STP_Bool* termsmark;
1384  STP_Bool* visited;
1385  int e;
1386  /* PHASE I: */
1387  for( i = 0; i < nnodes; i++ )
1388  g->mark[i] = (g->grad[i] > 0);
1389 
1390  /* allocate memory needed in PHASE I */
1391  SCIP_CALL( SCIPallocBufferArray(scip, &terms, nterms) );
1392  SCIP_CALL( SCIPallocBufferArray(scip, &termsmark, nnodes) );
1393  SCIP_CALL( SCIPallocBufferArray(scip, &vnoi, nnodes) );
1394  SCIP_CALL( SCIPallocBufferArray(scip, &visited, nnodes) );
1395  SCIP_CALL( SCIPallocBufferArray(scip, &reachednodes, nnodes) );
1396  SCIP_CALL( SCIPallocBufferArray(scip, &vbase, nnodes) );
1397  SCIP_CALL( SCIPallocBufferArray(scip, &tovisit, nnodes) );
1398  SCIP_CALL( SCIPallocBufferArray(scip, &vcost, nnodes) );
1399 
1400  j = 0;
1401  for( i = 0; i < nnodes; i++ )
1402  {
1403  visited[i] = FALSE;
1404  if( Is_term(g->term[i]) )
1405  {
1406  termsmark[i] = TRUE;
1407  terms[j++] = i;
1408  }
1409  else
1410  {
1411  termsmark[i] = FALSE;
1412  }
1413  }
1414 
1415  for( e = 0; e < g->edges; e++)
1416  {
1417  assert(SCIPisGE(scip, cost[e], 0.0));
1418  assert(SCIPisGE(scip, costrev[e], 0.0));
1419  }
1420 
1421  assert(j == nterms);
1422  graph_voronoi(scip, g, cost, costrev, termsmark, vbase, vnoi);
1423  state = g->path_state;
1424 
1425  for( i = 0; i < nnodes; i++ )
1426  if( Is_term(g->term[i]) )
1427  assert(vbase[i] == i);
1428 
1429  for( k = 0; k < nnodes; k++ )
1430  {
1431  connected[k] = FALSE;
1432  vcount[k] = 0;
1433  gnodearr[k]->number = k;
1434  if( !Is_term(g->term[k]) )
1435  {
1436  node_dist[k][0] = vnoi[k].dist;
1437  node_edge[k][0] = vnoi[k].edge;
1438  node_base[k][0] = vbase[k];
1439  nodenterms[k] = 1;
1440  }
1441  else
1442  {
1443  nodenterms[k] = 0;
1444  node_edge[k][0] = UNKNOWN;
1445  termsmark[k] = FALSE;
1446  }
1447  state[k] = UNKNOWN;
1448  vcost[k] = vnoi[k].dist;
1449  vnoi[k].dist = FARAWAY;
1450  }
1451 
1452  /* for each terminal: extend the voronoi regions until all neighbouring terminals have been visited */
1453  for( i = 0; i < nterms; i++ )
1454  {
1455  term = terms[i];
1456  nneighbterms = 0;
1457  nneighbnodes = 0;
1458  nreachednodes = 0;
1459  for( k = 0; k < nnodes; k++ )
1460  assert(termsmark[k] == FALSE);
1461  /* DFS (starting from terminal i) until the entire voronoi region has been visited */
1462  tovisit[0] = term;
1463  ntovisit = 1;
1464  visited[term] = TRUE;
1465  state[term] = CONNECT;
1466  while( ntovisit > 0 )
1467  {
1468  /* iterate all incident edges */
1469  old = tovisit[--ntovisit];
1470 
1471  for( oedge = g->outbeg[old]; oedge != EAT_LAST; oedge = g->oeat[oedge] )
1472  {
1473  k = g->head[oedge];
1474 
1475  /* is node k in the voronoi region of the i-th terminal ? */
1476  if( vbase[k] == term )
1477  {
1478  if( !visited[k] )
1479  {
1480  state[k] = CONNECT;
1481  assert(nnodes - (nneighbnodes + 1) > ntovisit);
1482  tovisit[ntovisit++] = k;
1483  visited[k] = TRUE;
1484  reachednodes[nreachednodes++] = k;
1485  }
1486  }
1487  else
1488  {
1489  if( !visited[k] )
1490  {
1491  visited[k] = TRUE;
1492  vnoi[k].dist = vcost[old] + ((vbase[k] == root)? cost[oedge] : costrev[oedge]);
1493  vnoi[k].edge = oedge;
1494 
1495  if( termsmark[vbase[k]] == FALSE )
1496  {
1497  termsmark[vbase[k]] = TRUE;
1498  nneighbterms++;
1499  }
1500  assert(nnodes - (nneighbnodes + 1) > ntovisit - 1);
1501  tovisit[nnodes - (++nneighbnodes)] = k;
1502  }
1503  else
1504  {
1505  /* if edge 'oedge' allows a shorter connection of node k, update */
1506  if( SCIPisGT(scip, vnoi[k].dist, vcost[old] + ((vbase[k] == root)? cost[oedge] : costrev[oedge])) )
1507  {
1508  vnoi[k].dist = vcost[old] + ((vbase[k] == root)? cost[oedge] : costrev[oedge]);
1509  vnoi[k].edge = oedge;
1510  }
1511  }
1512  }
1513  }
1514  }
1515 
1516  count = 0;
1517  for( j = 0; j < nneighbnodes; j++ )
1518  {
1519  assert(termsmark[vbase[tovisit[nnodes - j - 1]]]);
1520  heap_add(g->path_heap, state, &count, tovisit[nnodes - j - 1], vnoi);
1521  }
1522  SCIP_CALL( graph_voronoiExtend(scip, g, ((term == root)? cost : costrev), vnoi, node_dist, node_base, node_edge, termsmark, reachednodes, &nreachednodes, nodenterms,
1523  nneighbterms, term, nneighbnodes) );
1524 
1525  reachednodes[nreachednodes++] = term;
1526 
1527  for( j = 0; j < nreachednodes; j++ )
1528  {
1529  vnoi[reachednodes[j]].dist = FARAWAY;
1530  state[reachednodes[j]] = UNKNOWN;
1531  visited[reachednodes[j]] = FALSE;
1532  }
1533 
1534  for( j = 0; j < nneighbnodes; j++ )
1535  {
1536  vnoi[tovisit[nnodes - j - 1]].dist = FARAWAY;
1537  state[tovisit[nnodes - j - 1]] = UNKNOWN;
1538  visited[tovisit[nnodes - j - 1]] = FALSE;
1539  }
1540  }
1541 
1542  /* for each node v: sort the terminal arrays according to their distance to v */
1543  for( i = 0; i < nnodes && !SCIPisStopped(scip); i++ )
1544  SCIPsortRealIntInt(node_dist[i], node_base[i], node_edge[i], nodenterms[i]);
1545 
1546  /* free memory */
1547  SCIPfreeBufferArray(scip, &vcost);
1548  SCIPfreeBufferArray(scip, &tovisit);
1549  SCIPfreeBufferArray(scip, &vbase);
1550  SCIPfreeBufferArray(scip, &reachednodes);
1551  SCIPfreeBufferArray(scip, &visited);
1552  SCIPfreeBufferArray(scip, &vnoi);
1553  SCIPfreeBufferArray(scip, &termsmark);
1554  SCIPfreeBufferArray(scip, &terms);
1555  }
1556 
1557  /* PHASE II */
1558  else
1559  {
1560  for( k = 0; k < nnodes; k++ )
1561  {
1562  connected[k] = FALSE;
1563  vcount[k] = 0;
1564  }
1565  }
1566 
1567  connected[start] = TRUE;
1568  gnodearr[start]->dist = node_dist[start][0];
1569  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[start]) );
1570 
1571  while( SCIPpqueueNElems(pqueue) > 0 && !SCIPisStopped(scip) )
1572  {
1573  best = ((GNODE*) SCIPpqueueRemove(pqueue))->number;
1574 
1575  term = node_base[best][vcount[best]];
1576  assert( Is_term(g->term[term]) );
1577  /* has the terminal already been connected? */
1578  if( !connected[term] )
1579  {
1580  /* connect the terminal */
1581  k = g->tail[node_edge[best][vcount[best]]];
1582  while( k != term )
1583  {
1584  j = 0;
1585 
1586  while( node_base[k][vcount[k] + j] != term )
1587  j++;
1588 
1589  assert(vcount[k] + j < nodenterms[k]);
1590 
1591  if( !connected[k] )
1592  {
1593  assert(vcount[k] == 0);
1594 
1595  connected[k] = TRUE;
1596  while( vcount[k] < nodenterms[k] && connected[node_base[k][vcount[k]]] )
1597  {
1598  vcount[k]++;
1599  j--;
1600  }
1601 
1602  if( vcount[k] < nodenterms[k] )
1603  {
1604  gnodearr[k]->dist = node_dist[k][vcount[k]];
1605  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[k]) );
1606  }
1607  }
1608 
1609  assert( vcount[k] + j < nodenterms[k] );
1610  assert(node_base[k][vcount[k] + j] == term);
1611  k = g->tail[node_edge[k][vcount[k] + j]];
1612  }
1613 
1614  /* finally, connected the terminal */
1615  assert( k == term );
1616  assert( !connected[k] );
1617  connected[k] = TRUE;
1618 
1619  assert( vcount[k] == 0 );
1620  while( vcount[k] < nodenterms[k] && connected[node_base[k][vcount[k]]] )
1621  {
1622  vcount[k]++;
1623  }
1624  if( vcount[k] < nodenterms[k] )
1625  {
1626  gnodearr[k]->dist = node_dist[k][vcount[k]];
1627  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[k]) );
1628  }
1629  }
1630 
1631  while( vcount[best] + 1 < nodenterms[best] )
1632  {
1633  if( !connected[node_base[best][++vcount[best]]] )
1634  {
1635  gnodearr[best]->dist = node_dist[best][vcount[best]];
1636  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[best]) );
1637  break;
1638  }
1639  }
1640  }
1641 
1642  /* prune the ST, so that all leaves are terminals */
1643  SCIP_CALL( prune(scip, g, cost, result, connected) );
1644 
1645  return SCIP_OKAY;
1646 }
1647 
1648 /** execute shortest paths heuristic to obtain a Steiner tree */
1650  SCIP* scip, /**< SCIP data structure */
1651  SCIP_HEURDATA* heurdata, /**< SCIP data structure */
1652  GRAPH* graph, /**< graph data structure */
1653  int* starts, /**< array containing start vertices (NULL to not provide any) */
1654  int* bestnewstart, /**< pointer to the start vertex resulting in the best solution */
1655  int* best_result, /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
1656  int runs, /**< number of runs */
1657  int bestincstart, /**< best incumbent start vertex */
1658  SCIP_Real* cost, /**< arc costs */
1659  SCIP_Real* costrev, /**< reversed arc costs */
1660  SCIP_Real* hopfactor, /**< edge cost multiplicator for HC problems */
1661  SCIP_Real* nodepriority, /**< vertex priorities for vertices to be starting points (NULL for no priorities) */
1662  SCIP_Real maxcost, /**< maximal edge cost (only for HC) */
1663  SCIP_Bool* success, /**< pointer to store whether a solution could be found */
1664  SCIP_Bool pcmwfull /**< use full computation of tree (i.e. connect all terminals and prune), only for prize-collecting variants */
1665  )
1666 {
1667  SCIP_PQUEUE* pqueue = NULL;
1668  SCIP_Longint nexecs;
1669  SCIP_Real obj;
1670  SCIP_Real min = FARAWAY;
1671  SCIP_Real* dijkdist = NULL;
1672  SCIP_Real** pathdist = NULL;
1673  SCIP_Real** node_dist = NULL;
1674  SCIP_Bool lsuccess;
1675  GNODE** gnodearr = NULL;
1676  int best;
1677  int k;
1678  int r;
1679  int e;
1680  int root;
1681  int mode;
1682  int nedges;
1683  int nnodes;
1684  int nterms;
1685  int* perm = NULL;
1686  int* start;
1687  int* vcount = NULL;
1688  int* RESTRICT result;
1689  int* cluster = NULL;
1690  int* dijkedge = NULL;
1691  int* nodenterms = NULL;
1692  int** pathedge = NULL;
1693  int** node_base = NULL;
1694  int** node_edge = NULL;
1695  SCIP_Bool startsgiven;
1696  STP_Bool* connected;
1697 
1698  assert(scip != NULL);
1699  assert(cost != NULL);
1700  assert(graph != NULL);
1701  assert(costrev != NULL);
1702  assert(best_result != NULL);
1703 
1704  if( heurdata == NULL )
1705  {
1706  assert(SCIPfindHeur(scip, "TM") != NULL);
1707  heurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
1708  }
1709 
1710  best = bestincstart;
1711  root = graph->source;
1712  nnodes = graph->knots;
1713  nedges = graph->edges;
1714  nterms = graph->terms;
1715  nexecs = heurdata->nexecs;
1716  (*success) = FALSE;
1717 
1718  if( runs < 1 )
1719  {
1720  startsgiven = FALSE;
1721  runs = 1;
1722  }
1723  else
1724  startsgiven = (starts != NULL);
1725 
1726  for( e = 0; e < nedges; e++)
1727  {
1728  assert(SCIPisGE(scip, cost[e], 0.0));
1729  assert(SCIPisGE(scip, costrev[e], 0.0));
1730  }
1731 
1732  if( SCIPisStopped(scip) )
1733  return SCIP_OKAY;
1734 
1735  assert(root >= 0);
1736  assert(nedges > 0);
1737  assert(nnodes > 0);
1738  assert(nterms > 0);
1739  assert(nexecs >= 0);
1740 
1741  /* allocate memory */
1742  SCIP_CALL( SCIPallocBufferArray(scip, &connected, nnodes) );
1743  SCIP_CALL( SCIPallocBufferArray(scip, &start, MIN(runs, nnodes)) );
1744  SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
1745 
1746  /* get user parameter */
1747  mode = heurdata->type;
1748  assert(mode == AUTO || mode == TM_SP || mode == TM_VORONOI || mode == TM_DIJKSTRA);
1749 
1750  if( graph->stp_type == STP_RPCSPG || graph->stp_type == STP_PCSPG || graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP )
1751  {
1752  mode = TM_DIJKSTRA;
1753  graph_pc_2transcheck(graph);
1754  }
1755  else if( graph->stp_type == STP_DHCSTP )
1756  {
1757  mode = TM_SP;
1758  }
1759  else if( graph->stp_type == STP_DCSTP )
1760  {
1761  mode = TM_SP;
1762  }
1763  else if( graph->stp_type == STP_SAP )
1764  {
1765  mode = TM_SP;
1766  }
1767  else
1768  {
1769  if( mode == AUTO )
1770  mode = TM_DIJKSTRA;
1771  }
1772 
1773  /* allocate memory according to selected mode */
1774  if( mode == TM_DIJKSTRA || graph->stp_type == STP_DHCSTP )
1775  {
1776  SCIP_CALL( SCIPallocBufferArray(scip, &dijkdist, nnodes) );
1777  SCIP_CALL( SCIPallocBufferArray(scip, &dijkedge, nnodes) );
1778  }
1779  if( mode == TM_VORONOI )
1780  {
1781  SCIP_CALL( SCIPallocBufferArray(scip, &nodenterms, nnodes) );
1782  SCIP_CALL( SCIPallocBufferArray(scip, &gnodearr, nnodes) );
1783  SCIP_CALL( SCIPallocBufferArray(scip, &node_base, nnodes) );
1784  SCIP_CALL( SCIPallocBufferArray(scip, &node_dist, nnodes) );
1785  SCIP_CALL( SCIPallocBufferArray(scip, &node_edge, nnodes) );
1786 
1787  for( k = 0; k < nnodes; k++ )
1788  {
1789  SCIP_CALL( SCIPallocBuffer(scip, &gnodearr[k]) ); /*lint !e866*/
1790  SCIP_CALL( SCIPallocBufferArray(scip, &node_base[k], nterms) ); /*lint !e866*/
1791  SCIP_CALL( SCIPallocBufferArray(scip, &node_dist[k], nterms) ); /*lint !e866*/
1792  SCIP_CALL( SCIPallocBufferArray(scip, &node_edge[k], nterms) ); /*lint !e866*/
1793  }
1794 
1795  SCIP_CALL( SCIPallocBufferArray(scip, &vcount, nnodes) );
1796  SCIP_CALL( SCIPpqueueCreate( &pqueue, nnodes, 2.0, GNODECmpByDist) );
1797  }
1798  if( mode == TM_SP )
1799  {
1800  SCIP_CALL( SCIPallocBufferArray(scip, &pathdist, nnodes) );
1801  SCIP_CALL( SCIPallocBufferArray(scip, &pathedge, nnodes) );
1802  BMSclearMemoryArray(pathdist, nnodes);
1803  BMSclearMemoryArray(pathedge, nnodes);
1804 
1805  for( k = 0; k < nnodes; k++ )
1806  {
1807  graph->mark[k] = (graph->grad[k] > 0);
1808 
1809  if( Is_term(graph->term[k]) )
1810  {
1811  SCIP_CALL( SCIPallocBufferArray(scip, &(pathdist[k]), nnodes) ); /*lint !e866*/
1812  SCIP_CALL( SCIPallocBufferArray(scip, &(pathedge[k]), nnodes) ); /*lint !e866*/
1813  }
1814  }
1815 
1816  SCIP_CALL( SCIPallocBufferArray(scip, &cluster, nnodes) );
1817  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nnodes) );
1818  }
1819 
1820  /* compute number of iterations and starting points for SPH */
1821  if( graph->stp_type == STP_PCSPG || graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP || graph->stp_type == STP_RPCSPG )
1822  {
1823  if( runs > (nterms - 1) && graph->stp_type != STP_RMWCSP )
1824  runs = nterms - 1;
1825  else if( runs > (nterms) && graph->stp_type == STP_RMWCSP )
1826  runs = nterms;
1827  }
1828  else if( startsgiven )
1829  {
1830  for( k = 0; k < MIN(runs, nnodes); k++ )
1831  start[k] = starts[k];
1832  }
1833  else if( runs < nnodes || graph->stp_type == STP_DHCSTP )
1834  {
1835  if( best < 0 )
1836  {
1837  best = root;
1838  }
1839  else
1840  {
1841  int randint = SCIPrandomGetInt(heurdata->randnumgen, 0, 2);
1842  if( randint == 0 )
1843  best = -1;
1844  else if( randint == 1 )
1845  best = root;
1846  }
1847  r = 0;
1848  if( graph->stp_type == STP_DHCSTP )
1849  graph_path_execX(scip, graph, root, cost, dijkdist, dijkedge);
1850 
1851  /* allocate memory for permutation array */
1852  if( perm == NULL )
1853  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nnodes) );
1854 
1855  /* are there no nodes are to be priorized a priori? */
1856  if( nodepriority == NULL )
1857  {
1858  for( k = 0; k < nnodes; k++ )
1859  perm[k] = k;
1860  SCIPrandomPermuteIntArray(heurdata->randnumgen, perm, 0, nnodes);
1861 
1862  /* use terminals (randomly permutated) as starting points for TM heuristic */
1863  for( k = 0; k < nnodes; k++ )
1864  {
1865  if( r >= runs || r >= nterms )
1866  break;
1867 
1868  if( Is_term(graph->term[perm[k]]) )
1869  start[r++] = perm[k];
1870  }
1871 
1872  /* fill empty slots */
1873  for( k = 0; k < nnodes && r < runs; k++ )
1874  {
1875  if( graph->stp_type == STP_DHCSTP )
1876  {
1877  assert(dijkdist != NULL);
1878  if( SCIPisGE(scip, dijkdist[perm[k]], BLOCKED) )
1879  continue;
1880  }
1881  if( !Is_term(graph->term[perm[k]]) && graph->mark[k] )
1882  start[r++] = perm[k];
1883  }
1884  }
1885  else
1886  {
1887  SCIP_Real max = 0.0;
1888  int bbound = runs - runs / 3;
1889 
1890  for( k = 0; k < nnodes; k++ )
1891  {
1892  perm[k] = k;
1893  if( SCIPisLT(scip, max, nodepriority[k]) && Is_term(graph->term[k]) )
1894  max = nodepriority[k];
1895  }
1896  for( k = 0; k < nnodes; k++ )
1897  {
1898  if( Is_term(graph->term[k]) )
1899  {
1900  nodepriority[k] += SCIPrandomGetReal(heurdata->randnumgen, 0.0, max);
1901  }
1902  else if( SCIPisLE(scip, 1.0, nodepriority[k]) )
1903  {
1904  nodepriority[k] = nodepriority[k] * SCIPrandomGetReal(heurdata->randnumgen, 1.0, 2.0);
1905  }
1906  }
1907 
1908  SCIPsortRealInt(nodepriority, perm, nnodes);
1909 
1910  for( k = nnodes - 1; k >= 0; k-- )
1911  {
1912  if( r >= nterms || r >= bbound )
1913  break;
1914 
1915  if( Is_term(graph->term[perm[k]]) )
1916  {
1917  start[r++] = perm[k];
1918  perm[k] = -1;
1919  }
1920  }
1921 
1922  /* fill empty slots */
1923  for( k = nnodes - 1; k >= 0 && r < runs; k-- )
1924  {
1925  if( perm[k] == -1 )
1926  continue;
1927  if( graph->stp_type == STP_DHCSTP )
1928  {
1929  assert(dijkdist != NULL);
1930  if( SCIPisGE(scip, dijkdist[perm[k]], BLOCKED) )
1931  continue;
1932  }
1933  if( graph->mark[k] )
1934  start[r++] = perm[k];
1935  }
1936  }
1937  /* not all slots filled? */
1938  if( r < runs )
1939  runs = r;
1940 
1941  if( best >= 0 )
1942  {
1943  /* check whether we have a already selected the best starting node */
1944  for( r = 0; r < runs; r++ )
1945  if( start[r] == best )
1946  break;
1947 
1948  /* no, we still have to */
1949  if( r == runs && runs > 0 )
1950  start[nexecs % runs] = best;
1951  }
1952  } /* STP_DHCSTP */
1953  else
1954  {
1955  runs = nnodes;
1956  for( k = 0; k < nnodes; k++ )
1957  start[k] = k;
1958  }
1959 
1960  /* perform SPH computations, differentiate between STP variants */
1961  if( graph->stp_type == STP_PCSPG || graph->stp_type == STP_RPCSPG || graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP )
1962  {
1963  SCIP_Real* terminalprio;
1964  int* terminalperm;
1965  int t;
1966  const SCIP_Bool rmw = (graph->stp_type == STP_RMWCSP);
1967  const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
1968 
1969  /* allocate memory */
1970  SCIP_CALL( SCIPallocBufferArray(scip, &terminalperm, nterms) );
1971  SCIP_CALL( SCIPallocBufferArray(scip, &terminalprio, nterms) );
1972 
1973  /* initialize arrays */
1974 
1975  t = 0;
1976  if( nodepriority == NULL )
1977  {
1978  for( k = nnodes - 1; k >= 0; --k )
1979  if( Is_pterm(graph->term[k]) && graph->grad[k] > 0 )
1980  {
1981  assert(SCIPisGT(scip, graph->prize[k], 0.0));
1982  terminalperm[t] = k;
1983  terminalprio[t++] = SCIPrandomGetReal(heurdata->randnumgen, 0.0, graph->prize[k]);
1984  }
1985  }
1986  else
1987  {
1988  for( k = nnodes - 1; k >= 0; --k )
1989  if( Is_pterm(graph->term[k]) && graph->grad[k] > 0 )
1990  {
1991  assert(SCIPisGT(scip, graph->prize[k], 0.0));
1992  terminalperm[t] = k;
1993  terminalprio[t++] = SCIPrandomGetReal(heurdata->randnumgen, nodepriority[k] / 2.0, nodepriority[k]);
1994  }
1995  }
1996 
1997  if( rmw )
1998  {
1999  SCIP_Real max = 0.0;
2000  int head;
2001  for( k = t - 1; k >= 0; --k )
2002  {
2003  if( SCIPisGT(scip, terminalprio[k], max) )
2004  max = terminalprio[k];
2005  }
2006 
2007  for( k = nnodes - 1; k >= 0; --k )
2008  graph->mark[k] = (graph->grad[k] > 0);
2009 
2010  for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
2011  {
2012  if( SCIPisGT(scip, graph->cost[e], 0.0) && Is_term(graph->term[graph->head[e]]) )
2013  {
2014  head = graph->head[e];
2015  graph->mark[head] = FALSE;
2016  assert(graph->grad[head] == 2);
2017  }
2018  }
2019 
2020  for( k = nnodes - 1; k >= 0; --k )
2021  {
2022  if( Is_term(graph->term[k]) && graph->mark[k] )
2023  {
2024  assert(SCIPisGT(scip, graph->prize[k], 0.0));
2025  terminalperm[t] = k;
2026  terminalprio[t++] = SCIPrandomGetReal(heurdata->randnumgen, max / 2.0, 1.5 * max);
2027  }
2028  }
2029  assert(nterms == t);
2030  SCIPsortRealInt(terminalprio, terminalperm, nterms);
2031  }
2032  else
2033  {
2034  if( rpc )
2035  {
2036  terminalperm[t] = root;
2037  terminalprio[t++] = FARAWAY;
2038  runs++;
2039 
2040  SCIPsortRealInt(terminalprio, terminalperm, nterms);
2041  }
2042  else
2043  {
2044  SCIPsortRealInt(terminalprio, terminalperm, nterms - 1);
2045  }
2046  }
2047 
2048  if( best >= 0 && best < nnodes && Is_pterm(graph->term[best]) && SCIPrandomGetInt(heurdata->randnumgen, 0, 2) == 1 )
2049  terminalperm[nterms - 1] = best;
2050 
2051  /* local main loop */
2052  for( r = runs - 1; r >= 0; --r )
2053  {
2054  k = terminalperm[r];
2055 
2056  if( pcmwfull )
2057  {
2058  SCIP_CALL( computeSteinerTreeDijkPcMwFull(scip, graph, cost, dijkdist, result, dijkedge, k, connected) );
2059  }
2060  else
2061  {
2062  SCIP_CALL( computeSteinerTreeDijkPcMw(scip, graph, cost, dijkdist, result, dijkedge, k, connected) );
2063  }
2064 
2065  if( SCIPisStopped(scip) )
2066  break;
2067 
2068  /* compute objective value (wrt original costs) */
2069  obj = 0.0;
2070  for( e = nedges - 1; e >= 0; e-- ) /* todo: save result array as char put into computeSteinerTreeDijkPcMw */
2071  if( result[e] == CONNECT )
2072  obj += graph->cost[e];
2073 
2074  if( SCIPisLT(scip, obj, min) )
2075  {
2076  if( bestnewstart != NULL )
2077  *bestnewstart = k;
2078  min = obj;
2079 
2080  SCIPdebugMessage("\n Obj(run: %d, ncall: %d)=%.12e\n\n", r, (int) nexecs, obj);
2081 
2082  for( e = 0; e < nedges; e++ )
2083  best_result[e] = result[e];
2084  (*success) = TRUE;
2085  }
2086  }
2087 
2088  /* free memory */
2089  SCIPfreeBufferArray(scip, &terminalprio);
2090  SCIPfreeBufferArray(scip, &terminalperm);
2091  } /* STP_PCSPG or STP_(R)MWCSP */
2092  else
2093  {
2094  SCIP_Real* orgcost = NULL;
2095  int edgecount;
2096  STP_Bool solfound = FALSE;
2097 
2098  if( SCIPisLE(scip, maxcost, 0.0) )
2099  maxcost = 1.0;
2100 
2101  if( graph->stp_type == STP_DHCSTP )
2102  {
2103  SCIP_Real bestfactor = -1;
2104 
2105  assert(hopfactor != NULL);
2106  assert(SCIPisGT(scip, (*hopfactor), 0.0));
2107 
2108  SCIP_CALL( SCIPallocBufferArray(scip, &orgcost, nedges) );
2109 
2110  BMScopyMemoryArray(orgcost, cost, nedges);
2111 
2112  /* do a warm-up run */
2113  for( r = 0; r < 10; r++ )
2114  {
2115  for( e = 0; e < nedges; e++ )
2116  {
2117  if( (SCIPisLT(scip, cost[e], BLOCKED )) )
2118  cost[e] = 1.0 + orgcost[e] / ((*hopfactor) * maxcost);
2119  result[e] = UNKNOWN;
2120  }
2121 
2122  SCIP_CALL( computeSteinerTreeDijk(scip, graph, cost, dijkdist, result, dijkedge, root, heurdata->randnumgen, connected) );
2123 
2124  obj = 0.0;
2125  edgecount = 0;
2126 
2127  for( e = 0; e < nedges; e++)
2128  {
2129  if( result[e] == CONNECT )
2130  {
2131  obj += graph->cost[e];
2132  edgecount++;
2133  }
2134  }
2135  lsuccess = FALSE;
2136  if( SCIPisLT(scip, obj, min) && edgecount <= graph->hoplimit )
2137  {
2138  min = obj;
2139  if( bestnewstart != NULL )
2140  *bestnewstart = root;
2141 
2142  for( e = 0; e < nedges; e++ )
2143  best_result[e] = result[e];
2144  (*success) = TRUE;
2145  lsuccess = TRUE;
2146  bestfactor = (*hopfactor);
2147  }
2148 
2149  if( !lsuccess || SCIPisGT(scip, fabs((double) edgecount - graph->hoplimit) / (double) graph->hoplimit, 0.05) )
2150  {
2151  if( !lsuccess )
2152  {
2153  if( (*success) )
2154  {
2155  (*hopfactor) = (*hopfactor) * (1.0 + fabs((double) edgecount - graph->hoplimit) / (double) graph->hoplimit);
2156  }
2157  else
2158  {
2159  (*hopfactor) = (*hopfactor) * (1.0 + 3 * fabs((double) edgecount - graph->hoplimit) / (double) graph->hoplimit);
2160  bestfactor = (*hopfactor);
2161  }
2162  }
2163  else
2164  {
2165  (*hopfactor) = (*hopfactor) / (1.0 + fabs((double) edgecount - graph->hoplimit) / (double) graph->hoplimit);
2166  }
2167 
2168  assert(SCIPisGT(scip, (*hopfactor), 0.0));
2169  }
2170  else
2171  {
2172  break;
2173  }
2174  }
2175  (*hopfactor) = bestfactor;
2176 
2177  for( e = 0; e < nedges; e++ )
2178  if( (SCIPisLT(scip, cost[e], BLOCKED )) )
2179  cost[e] = 1.0 + orgcost[e] / ((*hopfactor) * maxcost);
2180  for( e = 0; e < nedges; e++)
2181  costrev[e] = cost[flipedge(e)];
2182  }
2183  for( r = 0; r < runs; r++ )
2184  {
2185  assert(start[r] >= 0);
2186  assert(start[r] < nnodes);
2187 
2188  if( mode == TM_DIJKSTRA && graph->stp_type != STP_DCSTP )
2189  {
2190  SCIP_CALL( computeSteinerTreeDijk(scip, graph, cost, dijkdist, result, dijkedge, start[r], heurdata->randnumgen, connected) );
2191  }
2192  else if( graph->stp_type == STP_DCSTP )
2193  {
2194  /* first run? */
2195  if( r == 0 )
2196  {
2197  assert(pathdist != NULL);
2198  assert(pathedge != NULL);
2199 
2200  for( k = 0; k < nnodes; k++ )
2201  graph->mark[k] = (graph->grad[k] > 0);
2202 
2203  /* initialize shortest paths from all terminals */
2204  for( k = 0; k < nnodes; k++ )
2205  {
2206  if( Is_term(graph->term[k]) )
2207  {
2208  if( root == k )
2209  graph_path_execX(scip, graph, k, cost, pathdist[k], pathedge[k]);
2210  else
2211  graph_path_execX(scip, graph, k, costrev, pathdist[k], pathedge[k]);
2212  }
2213  }
2214  }
2215  for( e = 0; e < nedges; e++ )
2216  result[e] = UNKNOWN;
2217 
2218  SCIP_CALL( computeDegConsTree(scip, graph, cost, costrev, pathdist, start[r], perm, result, cluster, pathedge, heurdata->randnumgen, connected, &solfound) );
2219  }
2220  else if( mode == TM_SP )
2221  {
2222  for( e = nedges - 1; e >= 0; --e )
2223  result[e] = UNKNOWN;
2224 
2225  if( r == 0 )
2226  {
2227  int i;
2228  assert(pathdist != NULL);
2229  assert(pathedge != NULL);
2230 
2231  for( i = 0; i < nnodes; i++ )
2232  graph->mark[i] = (graph->grad[i] > 0);
2233 
2234  /* initialize shortest paths from all terminals */
2235  for( k = 0; k < nnodes; k++ )
2236  {
2237  if( Is_term(graph->term[k]) )
2238  {
2239  if( root == k )
2240  graph_path_execX(scip, graph, k, cost, pathdist[k], pathedge[k]);
2241  else
2242  graph_path_execX(scip, graph, k, costrev, pathdist[k], pathedge[k]);
2243  }
2244  }
2245  }
2246 
2247  SCIP_CALL( computeSteinerTree(scip, graph, cost, costrev, pathdist, start[r], perm, result, cluster, pathedge, connected, heurdata->randnumgen) );
2248  }
2249  else
2250  {
2251  SCIP_CALL( computeSteinerTreeVnoi(scip, graph, pqueue, gnodearr, cost, costrev, node_dist, start[r], result, vcount,
2252  nodenterms, node_base, node_edge, (r == 0), connected) );
2253  }
2254  obj = 0.0;
2255  edgecount = 0;
2256 
2257  /* here another measure than in the do_(...) heuristics is being used */
2258  for( e = 0; e < nedges; e++)
2259  {
2260  if( result[e] >= 0 )
2261  {
2262  obj += graph->cost[e];
2263  edgecount++;
2264  }
2265  }
2266 
2267  SCIPdebugMessage(" Obj=%.12e\n", obj);
2268 
2269  if( SCIPisLT(scip, obj, min) && (graph->stp_type != STP_DCSTP || solfound) && !SCIPisStopped(scip) && r < runs )
2270  {
2271  if( graph->stp_type != STP_DHCSTP || edgecount <= graph->hoplimit )
2272  {
2273  min = obj;
2274  if( bestnewstart != NULL )
2275  *bestnewstart = start[r];
2276 
2277  for( e = 0; e < nedges; e++ )
2278  best_result[e] = result[e];
2279  (*success) = TRUE;
2280  }
2281  }
2282 
2283  /* time limit exceeded?*/
2284  if( SCIPisStopped(scip) )
2285  break;
2286  }
2287 
2288  if( graph->stp_type == STP_DHCSTP )
2289  {
2290  assert(orgcost != NULL);
2291  for( e = 0; e < nedges; e++ )
2292  {
2293  cost[e] = orgcost[e];
2294  costrev[e] = orgcost[flipedge(e)];
2295  }
2296  SCIPfreeBufferArray(scip, &orgcost);
2297  }
2298  }
2299 
2300  /* free allocated memory */
2301  SCIPfreeBufferArrayNull(scip, &perm);
2302  if( mode == TM_SP )
2303  {
2304  assert(pathedge != NULL);
2305  assert(pathdist != NULL);
2306  SCIPfreeBufferArray(scip, &cluster);
2307  for( k = nnodes - 1; k >= 0; k-- )
2308  {
2309  SCIPfreeBufferArrayNull(scip, &(pathedge[k]));
2310  SCIPfreeBufferArrayNull(scip, &(pathdist[k]));
2311  }
2312 
2313  SCIPfreeBufferArray(scip, &pathedge);
2314  SCIPfreeBufferArray(scip, &pathdist);
2315  }
2316  else if( mode == TM_VORONOI )
2317  {
2318  SCIPpqueueFree(&pqueue);
2319 
2320  SCIPfreeBufferArray(scip, &vcount);
2321  assert(node_edge != NULL);
2322  assert(node_dist != NULL);
2323  assert(node_base != NULL);
2324  assert(gnodearr != NULL);
2325  for( k = nnodes - 1; k >= 0; k-- )
2326  {
2327  SCIPfreeBufferArray(scip, &node_edge[k]);
2328  SCIPfreeBufferArray(scip, &node_dist[k]);
2329  SCIPfreeBufferArray(scip, &node_base[k]);
2330  SCIPfreeBuffer(scip, &gnodearr[k]);
2331  }
2332  SCIPfreeBufferArray(scip, &node_edge);
2333  SCIPfreeBufferArray(scip, &node_dist);
2334  SCIPfreeBufferArray(scip, &node_base);
2335  SCIPfreeBufferArray(scip, &gnodearr);
2336  SCIPfreeBufferArray(scip, &nodenterms);
2337  }
2338  if( mode == TM_DIJKSTRA || graph->stp_type == STP_DHCSTP )
2339  {
2340  SCIPfreeBufferArray(scip, &dijkedge);
2341  SCIPfreeBufferArray(scip, &dijkdist);
2342  }
2343 
2344  SCIPfreeBufferArray(scip, &result);
2345  SCIPfreeBufferArray(scip, &start);
2346  SCIPfreeBufferArray(scip, &connected);
2347 
2348  return SCIP_OKAY;
2349 }
2350 
2351 /** run shortest path heuristic, but bias edge costs towards best current LP solution */
2353  SCIP* scip, /**< SCIP data structure */
2354  GRAPH* graph, /**< graph data structure */
2355  SCIP_HEUR* heur, /**< heuristic or NULL */
2356  int* result, /**< array indicating whether an arc is part of the solution (CONNECTED/UNKNOWN) */
2357  int runs, /**< number of runs */
2358  SCIP_Real* cost, /**< arc costs (uninitialized) */
2359  SCIP_Real* costrev, /**< reversed arc costs (uninitialized) */
2360  SCIP_Bool* success /**< pointer to store whether a solution could be found */
2361  )
2362 {
2363  SCIP_VAR** vars;
2364  SCIP_HEURDATA* heurdata;
2365  SCIP_Real* xval;
2366  SCIP_Real* nodepriority = NULL;
2367  SCIP_Real maxcost = 0.0;
2368  int beststart;
2369  const int nedges = graph->edges;
2370  const int nnodes = graph->knots;
2371 
2372  assert(scip != NULL);
2373  assert(graph != NULL);
2374  assert(result != NULL);
2375  assert(cost != NULL);
2376  assert(costrev != NULL);
2377  assert(success != NULL);
2378 
2379  assert(SCIPfindHeur(scip, "TM") != NULL);
2380  heurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
2381 
2382  vars = SCIPprobdataGetVars(scip);
2383  assert(vars != NULL);
2384 
2385  /* LP was not solved */
2387  {
2388  xval = NULL;
2389  }
2390  else
2391  {
2392  SCIP_SOL* sol = NULL;
2393  SCIP_CALL(SCIPcreateSol(scip, &sol, heur));
2394 
2395  /* copy the current LP solution to the working solution */
2396  SCIP_CALL(SCIPlinkLPSol(scip, sol));
2397 
2398  xval = SCIPprobdataGetXval(scip, sol);
2399 
2400  SCIP_CALL(SCIPfreeSol(scip, &sol));
2401  }
2402 
2403  /* set (edge) result array to default */
2404  for( int e = 0; e < nedges; e++ )
2405  result[e] = UNKNOWN;
2406 
2407  {
2408  const SCIP_Real randupper = SCIPrandomGetReal(heurdata->randnumgen, 1.1, 2.5);
2409  const SCIP_Real randlower = SCIPrandomGetReal(heurdata->randnumgen, 1.1, randupper);
2410 
2411  if( xval == NULL )
2412  {
2413  BMScopyMemoryArray(cost, graph->cost, nedges);
2414 
2415  /* hop constraint problem? */
2416  if( graph->stp_type == STP_DHCSTP )
2417  {
2418  for( int e = 0; e < nedges; e++ )
2419  {
2420  if( SCIPvarGetUbGlobal(vars[e]) < 0.5 )
2421  cost[e] = BLOCKED;
2422  else if( SCIPisGT(scip, cost[e], maxcost) && SCIPisLT(scip, cost[e], FARAWAY) )
2423  maxcost = cost[e];
2424  }
2425  for( int e = 0; e < nedges; e++ )
2426  costrev[e] = cost[flipedge(e)];
2427  }
2428  else
2429  {
2430  if( graph->stp_type != STP_PCSPG && graph->stp_type != STP_MWCSP )
2431  {
2432  SCIP_CALL(SCIPallocBufferArray(scip, &nodepriority, nnodes));
2433  for( int k = 0; k < nnodes; k++ )
2434  {
2435  if( Is_term(graph->term[k]) )
2436  nodepriority[k] = (SCIP_Real) graph->grad[k];
2437  else
2438  nodepriority[k] = SCIPrandomGetReal(heurdata->randnumgen, 0.0, 1.0);
2439  }
2440  }
2441 
2442  for( int e = 0; e < nedges; e += 2 )
2443  {
2444  if( SCIPvarGetUbGlobal(vars[e + 1]) < 0.5 )
2445  {
2446  costrev[e] = BLOCKED;
2447  cost[e + 1] = BLOCKED;
2448  }
2449  else
2450  {
2451  costrev[e] = cost[e + 1];
2452  }
2453 
2454  if( SCIPvarGetUbGlobal(vars[e]) < 0.5 )
2455  {
2456  costrev[e + 1] = BLOCKED;
2457  cost[e] = BLOCKED;
2458  }
2459  else
2460  {
2461  costrev[e + 1] = cost[e];
2462  }
2463  }
2464  }
2465  }
2466  else
2467  {
2468  SCIP_Bool partrand = FALSE;
2469  SCIP_Bool totalrand = FALSE;
2470 
2471  if( (heurdata->nlpiterations == SCIPgetNLPIterations(scip) && SCIPrandomGetInt(heurdata->randnumgen, 0, 5) != 1)
2472  || SCIPrandomGetInt(heurdata->randnumgen, 0, 15) == 5 )
2473  partrand = TRUE;
2474 
2475  if( !partrand && (heurdata->nlpiterations == SCIPgetNLPIterations(scip) ) )
2476  totalrand = TRUE;
2477  else if( graph->stp_type == STP_DCSTP && heurdata->ncalls != 1 && SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 1
2478  && (graph->maxdeg[graph->source] == 1 || SCIPrandomGetInt(heurdata->randnumgen, 0, 5) == 5) )
2479  {
2480  totalrand = TRUE;
2481  partrand = FALSE;
2482  }
2483 
2484  assert(nodepriority == NULL);
2485 
2486  SCIP_CALL(SCIPallocBufferArray(scip, &nodepriority, nnodes));
2487 
2488  if( graph->stp_type != STP_MWCSP && graph->stp_type != STP_RMWCSP )
2489  {
2490  for( int k = 0; k < nnodes; k++ )
2491  nodepriority[k] = 0.0;
2492 
2493  for( int e = 0; e < nedges; e++ )
2494  {
2495  nodepriority[graph->head[e]] += xval[e];
2496  nodepriority[graph->tail[e]] += xval[e];
2497  }
2498  }
2499 
2500  if( graph->stp_type == STP_DHCSTP )
2501  {
2502  for( int e = 0; e < nedges; e++ )
2503  {
2504  if( SCIPvarGetUbGlobal(vars[e]) < 0.5 )
2505  {
2506  cost[e] = BLOCKED;
2507  }
2508  else
2509  {
2510  if( totalrand )
2511  {
2512  const SCIP_Real randval = SCIPrandomGetReal(heurdata->randnumgen, randlower, randupper);
2513  cost[e] = graph->cost[e] * randval;
2514  }
2515  else
2516  {
2517  cost[e] = ((1.0 - xval[e]) * graph->cost[e]);
2518  }
2519  }
2520  if( partrand )
2521  {
2522  const SCIP_Real randval = SCIPrandomGetReal(heurdata->randnumgen, randlower, randupper);
2523  cost[e] = cost[e] * randval;
2524  }
2525  if( SCIPisLT(scip, cost[e], BLOCKED) && SCIPisGT(scip, cost[e], maxcost) )
2526  maxcost = cost[e];
2527  assert(SCIPisGE(scip, cost[e], 0.0));
2528  }
2529  for( int e = 0; e < nedges; e++ )
2530  costrev[e] = cost[flipedge(e)];
2531  }
2532  else
2533  {
2534  /* swap costs; set a high cost if the variable is fixed to 0 */
2535  if( graph->stp_type == STP_MWCSP || graph->stp_type == STP_RMWCSP )
2536  {
2537  for( int e = 0; e < nnodes; e++ )
2538  nodepriority[e] = 0.0;
2539  for( int e = 0; e < nedges; e++ )
2540  nodepriority[graph->head[e]] += xval[e];
2541 
2542  for( int e = 0; e < nedges; e++ )
2543  {
2544  if( graph->cost[e] >= FARAWAY )
2545  cost[e] = graph->cost[e];
2546 
2547  if( SCIPvarGetUbLocal(vars[e]) < 0.5 )
2548  cost[e] = BLOCKED;
2549  else
2550  cost[e] = graph->cost[e] * (1.0 - MIN(1.0, nodepriority[graph->head[e]]));
2551  }
2552 
2553  for( int e = 0; e < nedges; e++ )
2554  {
2555  nodepriority[graph->tail[e]] += xval[e];
2556  costrev[flipedge(e)] = cost[e];
2557  }
2558  }
2559  else
2560  {
2561  for( int e = 0; e < nedges; e += 2 )
2562  {
2563  const SCIP_Real randval = SCIPrandomGetReal(heurdata->randnumgen, randlower, randupper);
2564 
2565  if( SCIPvarGetUbLocal(vars[e + 1]) < 0.5 )
2566  {
2567  costrev[e] = BLOCKED;
2568  cost[e + 1] = BLOCKED;
2569  }
2570  else
2571  {
2572  if( totalrand )
2573  costrev[e] = graph->cost[e + 1] * randval;
2574  else
2575  costrev[e] = ((1.0 - xval[e + 1]) * graph->cost[e + 1]);
2576 
2577  if( partrand )
2578  costrev[e] = costrev[e] * randval;
2579 
2580  cost[e + 1] = costrev[e];
2581  }
2582 
2583  if( SCIPvarGetUbLocal(vars[e]) < 0.5 )
2584  {
2585  costrev[e + 1] = BLOCKED;
2586  cost[e] = BLOCKED;
2587  }
2588  else
2589  {
2590  if( totalrand )
2591  costrev[e + 1] = graph->cost[e] * randval;
2592  else
2593  costrev[e + 1] = ((1.0 - xval[e]) * graph->cost[e]);
2594 
2595  if( partrand )
2596  costrev[e + 1] = costrev[e + 1] * randval;
2597  cost[e] = costrev[e + 1];
2598  }
2599  assert(SCIPisGE(scip, cost[e], 0.0));
2600  assert(SCIPisGE(scip, costrev[e], 0.0));
2601  }
2602  }
2603  }
2604  }
2605 
2606  for( int e = 0; e < nedges; e++ )
2607  {
2608  if( SCIPisZero(scip, cost[e]) )
2609  {
2610  cost[e] = SCIPepsilon(scip) * 2.0;
2611  assert(!SCIPisZero(scip, cost[e]));
2612  assert(SCIPisZero(scip, costrev[flipedge(e)]));
2613  costrev[flipedge(e)] = cost[e];
2614  }
2615  }
2616 
2617  /* can we connect the network */
2618  SCIP_CALL( SCIPStpHeurTMRun(scip, heurdata, graph, NULL, &beststart, result, runs, heurdata->beststartnode,
2619  cost, costrev, &(heurdata->hopfactor), nodepriority, maxcost, success, FALSE) );
2620  }
2621 
2622  SCIPfreeBufferArrayNull(scip, &nodepriority);
2623 
2624  return SCIP_OKAY;
2625 }
2626 
2627 
2628 /*
2629  * Callback methods of primal heuristic
2630  */
2631 
2632 
2633 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
2634 static
2636 { /*lint --e{715}*/
2637  assert(scip != NULL);
2638  assert(heur != NULL);
2639  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2640 
2641  /* call inclusion method of primal heuristic */
2643 
2644  return SCIP_OKAY;
2645 }
2646 
2647 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
2648 static
2650 { /*lint --e{715}*/
2651  SCIP_HEURDATA* heurdata;
2652 
2653  assert(heur != NULL);
2654  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2655  assert(scip != NULL);
2656 
2657  heurdata = SCIPheurGetData(heur);
2658  assert(heurdata != NULL);
2659 
2660  /* free random number generator */
2661  SCIPfreeRandom(scip, &heurdata->randnumgen);
2662 
2663  /* free heuristic data */
2664  SCIPfreeMemory(scip, &heurdata);
2665  SCIPheurSetData(heur, NULL);
2666 
2667  return SCIP_OKAY;
2668 }
2669 
2670 
2671 /** initialization method of primal heuristic (called after problem was transformed) */
2672 static
2674 { /*lint --e{715}*/
2675  SCIP_HEURDATA* heurdata;
2676  SCIP_PROBDATA* probdata;
2677  GRAPH* graph;
2678 
2679  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2680 
2681  heurdata = SCIPheurGetData(heur);
2682  assert(heurdata != NULL);
2683 
2684  probdata = SCIPgetProbData(scip);
2685  assert(probdata != NULL);
2686 
2687  graph = SCIPprobdataGetGraph(probdata);
2688  if( graph == NULL )
2689  {
2690  heurdata->stp_type = STP_SPG;
2691  return SCIP_OKAY;
2692  }
2693  heurdata->stp_type = graph->stp_type;
2694  heurdata->beststartnode = -1;
2695  heurdata->ncalls = 0;
2696  heurdata->nlpiterations = -1;
2697  heurdata->nexecs = 0;
2698 
2699 #ifdef WITH_UG
2700  heurdata->randseed += getUgRank();
2701 #endif
2702 
2703  return SCIP_OKAY;
2704 }
2705 
2706 /** execution method of primal heuristic */
2707 static
2709 { /*lint --e{715}*/
2710  SCIP_VAR** vars;
2711  SCIP_PROBDATA* probdata;
2712  SCIP_HEURDATA* heurdata;
2713  GRAPH* graph;
2714  SCIP_Real* cost;
2715  SCIP_Real* costrev;
2716  int* soledges;
2717  int runs;
2718  int nedges;
2719  SCIP_Bool success = FALSE;
2720 
2721  assert(scip != NULL);
2722  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
2723  assert(scip != NULL);
2724  assert(result != NULL);
2725 
2726  *result = SCIP_DIDNOTRUN;
2727 
2728  /* get heuristic data */
2729  heurdata = SCIPheurGetData(heur);
2730  assert(heurdata != NULL);
2731 
2732  probdata = SCIPgetProbData(scip);
2733  assert(probdata != NULL);
2734 
2735  graph = SCIPprobdataGetGraph(probdata);
2736  assert(graph != NULL);
2737 
2738  runs = 0;
2739 
2740  /* set the runs, i.e. number of different starting points for the heuristic */
2741  if( heurtiming & SCIP_HEURTIMING_BEFORENODE )
2742  {
2743  if( SCIPgetDepth(scip) > 0 )
2744  return SCIP_OKAY;
2745 
2746  runs = heurdata->initruns;
2747  }
2748  else if( ((heurtiming & SCIP_HEURTIMING_DURINGLPLOOP) && (heurdata->ncalls % heurdata->duringlpfreq == 0)) || (heurtiming & SCIP_HEURTIMING_AFTERLPLOOP) )
2749  {
2750  if( graph->stp_type == STP_PCSPG || graph->stp_type == STP_MWCSP )
2751  runs = 2 * heurdata->evalruns;
2752  else
2753  runs = heurdata->evalruns;
2754  }
2755  else if( heurtiming & SCIP_HEURTIMING_AFTERNODE )
2756  {
2757  if( SCIPgetDepth(scip) == 0 )
2758  runs = heurdata->rootruns;
2759  else
2760  runs = heurdata->leafruns;
2761  }
2762 
2763  /* increase counter for number of (TM) calls */
2764  heurdata->ncalls++;
2765 
2766  if( runs == 0 )
2767  return SCIP_OKAY;
2768 
2769  heurdata->nexecs++;
2770 
2771  SCIPdebugMessage("Heuristic Start\n");
2772 
2773  /* get all variables (corresponding to the edges) */
2774  vars = SCIPprobdataGetVars(scip);
2775  if( vars == NULL )
2776  return SCIP_OKAY;
2777 
2778  assert(vars[0] != NULL);
2779 
2780  nedges = graph->edges;
2781 
2782  /* allocate memory */
2783  SCIP_CALL(SCIPallocBufferArray(scip, &cost, nedges));
2784  SCIP_CALL(SCIPallocBufferArray(scip, &costrev, nedges));
2785  SCIP_CALL(SCIPallocBufferArray(scip, &soledges, nedges));
2786 
2787  *result = SCIP_DIDNOTFIND;
2788 
2789  /* call the actual heuristic */
2790  SCIP_CALL( SCIPStpHeurTMRunLP(scip, graph, heur, soledges, runs, cost, costrev, &success) );
2791 
2792  if( success )
2793  {
2794  SCIP_Real* nval;
2795  const int nvars = SCIPprobdataGetNVars(scip);
2796 
2797  SCIP_CALL(SCIPallocBufferArray(scip, &nval, nvars));
2798 
2799  for( int v = 0; v < nvars; v++ )
2800  nval[v] = (soledges[v % nedges] == (v / nedges)) ? 1.0 : 0.0;
2801 
2802  SCIP_CALL( SCIPStpValidateSol(scip, graph, nval, &success) );
2803  if( success )
2804  {
2805  SCIP_SOL* sol = NULL;
2806  SCIP_CALL( SCIPprobdataAddNewSol(scip, nval, sol, heur, &success) );
2807 
2808  if( success )
2809  {
2810  SCIPdebugMessage("TM solution added, value %f \n",
2811  graph_sol_getObj(graph->cost, soledges, SCIPprobdataGetOffset(scip), nedges));
2812 
2813  *result = SCIP_FOUNDSOL;
2814  }
2815  }
2816  SCIPfreeBufferArray(scip, &nval);
2817  }
2818 
2819  heurdata->nlpiterations = SCIPgetNLPIterations(scip);
2820  SCIPfreeBufferArray(scip, &soledges);
2821  SCIPfreeBufferArray(scip, &costrev);
2822  SCIPfreeBufferArray(scip, &cost);
2823 
2824  return SCIP_OKAY;
2825 }
2826 
2827 /*
2828  * primal heuristic specific interface methods
2829  */
2830 
2831 /** creates the TM primal heuristic and includes it in SCIP */
2833  SCIP* scip /**< SCIP data structure */
2834  )
2835 {
2836  SCIP_HEURDATA* heurdata;
2837  SCIP_HEUR* heur;
2838  char paramdesc[SCIP_MAXSTRLEN];
2839 
2840  /* create TM primal heuristic data */
2841  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
2842  heur = NULL;
2843 
2844  /* include primal heuristic */
2845  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
2847  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecTM, heurdata) );
2848 
2849  assert(heur != NULL);
2850 
2851  /* set non fundamental callbacks via setter functions */
2852  SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyTM) );
2853  SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeTM) );
2854  SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitTM) );
2855 
2856  heurdata->ncalls = 0;
2857  heurdata->nlpiterations = -1;
2858  heurdata->nexecs = 0;
2859  heurdata->randseed = DEFAULT_RANDSEED;
2860 
2861  /* add TM primal heuristic parameters */
2862  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/evalruns",
2863  "number of runs for eval",
2864  &heurdata->evalruns, FALSE, DEFAULT_EVALRUNS, -1, INT_MAX, NULL, NULL) );
2865  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/randseed",
2866  "random seed for heuristic",
2867  NULL, FALSE, DEFAULT_RANDSEED, 1, INT_MAX, paramChgdRandomseed, (SCIP_PARAMDATA*)heurdata) );
2868  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/initruns",
2869  "number of runs for init",
2870  &heurdata->initruns, FALSE, DEFAULT_INITRUNS, -1, INT_MAX, NULL, NULL) );
2871  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/leafruns",
2872  "number of runs for leaf",
2873  &heurdata->leafruns, FALSE, DEFAULT_LEAFRUNS, -1, INT_MAX, NULL, NULL) );
2874  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/rootruns",
2875  "number of runs for root",
2876  &heurdata->rootruns, FALSE, DEFAULT_ROOTRUNS, -1, INT_MAX, NULL, NULL) );
2877  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/duringlpfreq",
2878  "frequency for calling heuristic during LP loop",
2879  &heurdata->duringlpfreq, FALSE, DEFAULT_DURINGLPFREQ, 1, INT_MAX, NULL, NULL) );
2880  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/type",
2881  "Heuristic: 0 automatic, 1 TM_SP, 2 TM_VORONOI, 3 TM_DIJKSTRA",
2882  &heurdata->type, FALSE, DEFAULT_TYPE, 0, 3, NULL, NULL) );
2883  heurdata->hopfactor = DEFAULT_HOPFACTOR;
2884 
2885 #ifdef WITH_UG
2886  heurdata->randseed += getUgRank();
2887 #endif
2888 
2889  /* create random number generator */
2890  SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, heurdata->randseed, TRUE) );
2891 
2892  (void) SCIPsnprintf(paramdesc, SCIP_MAXSTRLEN, "timing when heuristic should be called (%u:BEFORENODE, %u:DURINGLPLOOP, %u:AFTERLPLOOP, %u:AFTERNODE)", SCIP_HEURTIMING_BEFORENODE, SCIP_HEURTIMING_DURINGLPLOOP, SCIP_HEURTIMING_AFTERLPLOOP, SCIP_HEURTIMING_AFTERNODE);
2893  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/timing", paramdesc,
2894  (int*) &heurdata->timing, TRUE, (int) HEUR_TIMING, (int) SCIP_HEURTIMING_BEFORENODE, 2 * (int) SCIP_HEURTIMING_AFTERPSEUDONODE - 1, NULL, NULL) ); /*lint !e713*/
2895 
2896  return SCIP_OKAY;
2897 }
#define HEUR_TIMING
Definition: heur_tm.c:46
#define HEUR_FREQ
Definition: heur_tm.c:43
static SCIP_DECL_HEURFREE(heurFreeTM)
Definition: heur_tm.c:2649
void * SCIPpqueueRemove(SCIP_PQUEUE *pqueue)
Definition: misc.c:1307
#define SCIP_HEURTIMING_DURINGLPLOOP
Definition: type_timing.h:71
static volatile int nterms
Definition: interrupt.c:37
SCIP_RETCODE SCIPfreeSol(SCIP *scip, SCIP_SOL **sol)
Definition: scip_sol.c:976
SCIP_Real SCIPepsilon(SCIP *scip)
static SCIP_RETCODE computeDegConsTree(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, SCIP_Real **pathdist, int start, int *perm, int *result, int *cluster, int **pathedge, SCIP_RANDNUMGEN *randnumgen, STP_Bool *connected, STP_Bool *solfound)
Definition: heur_tm.c:1073
#define NULL
Definition: def.h:253
void graph_voronoi(SCIP *scip, const GRAPH *, SCIP_Real *, SCIP_Real *, STP_Bool *, int *, PATH *)
Definition: grphpath.c:1751
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:686
void SCIPpqueueFree(SCIP_PQUEUE **pqueue)
Definition: misc.c:1259
void graph_path_st_rpc(SCIP *, const GRAPH *, const SCIP_Real *, SCIP_Real *, int *, int, STP_Bool *)
Definition: grphpath.c:1033
void graph_path_st(SCIP *, const GRAPH *, SCIP_Real *, SCIP_Real *, int *, int, SCIP_RANDNUMGEN *, STP_Bool *)
Definition: grphpath.c:926
const char * SCIPheurGetName(SCIP_HEUR *heur)
Definition: heur.c:1254
int *RESTRICT head
Definition: grph.h:96
Definition: grph.h:57
SCIP_RETCODE SCIPStpHeurTMBuildTreeDc(SCIP *scip, const GRAPH *g, int *result, STP_Bool *connected)
Definition: heur_tm.c:732
int source
Definition: grph.h:67
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_HEURDATA * SCIPheurGetData(SCIP_HEUR *heur)
Definition: heur.c:1165
SCIP_PARAMDATA * SCIPparamGetData(SCIP_PARAM *param)
Definition: paramset.c:661
signed int edge
Definition: grph.h:146
#define DEFAULT_HOPFACTOR
Definition: heur_tm.h:37
#define MST_MODE
Definition: grph.h:162
SCIP_RETCODE SCIPStpIncludeHeurTM(SCIP *scip)
Definition: heur_tm.c:2832
#define SCIP_MAXSTRLEN
Definition: def.h:274
SCIP_Bool graph_valid(const GRAPH *)
Definition: grphbase.c:4324
int terms
Definition: grph.h:64
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
static SCIP_DECL_HEURINIT(heurInitTM)
Definition: heur_tm.c:2673
struct SCIP_ParamData SCIP_PARAMDATA
Definition: type_paramset.h:76
#define SCIP_HEURTIMING_BEFORENODE
Definition: type_timing.h:70
#define HEUR_DISPCHAR
Definition: heur_tm.c:41
int *RESTRICT maxdeg
Definition: grph.h:78
#define EAT_LAST
Definition: grph.h:31
#define RESTRICT
Definition: def.h:265
void graph_path_st_rmw(SCIP *, const GRAPH *, const SCIP_Real *, SCIP_Real *, int *, int, STP_Bool *)
Definition: grphpath.c:1536
static SCIP_DECL_HEUREXEC(heurExecTM)
Definition: heur_tm.c:2708
#define BLOCKED
Definition: grph.h:157
#define DEFAULT_LEAFRUNS
Definition: heur_tm.c:51
#define FALSE
Definition: def.h:73
int *RESTRICT inpbeg
Definition: grph.h:74
int *RESTRICT path_state
Definition: grph.h:119
SCIP_Real SCIPrandomGetReal(SCIP_RANDNUMGEN *randnumgen, SCIP_Real minrandval, SCIP_Real maxrandval)
Definition: misc.c:9640
#define STP_RMWCSP
Definition: grph.h:50
Problem data for stp problem.
#define TRUE
Definition: def.h:72
#define SCIPdebug(x)
Definition: pub_message.h:74
SCIP_RETCODE SCIPStpHeurTMRun(SCIP *scip, SCIP_HEURDATA *heurdata, GRAPH *graph, int *starts, int *bestnewstart, int *best_result, int runs, int bestincstart, SCIP_Real *cost, SCIP_Real *costrev, SCIP_Real *hopfactor, SCIP_Real *nodepriority, SCIP_Real maxcost, SCIP_Bool *success, SCIP_Bool pcmwfull)
Definition: heur_tm.c:1649
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
SCIP_RETCODE SCIPlinkLPSol(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1017
#define SCIP_HEURTIMING_AFTERNODE
Definition: type_timing.h:92
#define STP_DHCSTP
Definition: grph.h:48
int SCIPpqueueNElems(SCIP_PQUEUE *pqueue)
Definition: misc.c:1362
SCIP_PROBDATA * SCIPgetProbData(SCIP *scip)
Definition: scip_prob.c:963
int SCIPprobdataGetNVars(SCIP *scip)
#define SCIP_HEURTIMING_AFTERLPLOOP
Definition: type_timing.h:72
struct SCIP_HeurData SCIP_HEURDATA
Definition: type_heur.h:51
static SCIP_RETCODE computeSteinerTreeVnoi(SCIP *scip, const GRAPH *g, SCIP_PQUEUE *pqueue, GNODE **gnodearr, SCIP_Real *cost, SCIP_Real *costrev, SCIP_Real **node_dist, int start, int *result, int *vcount, int *nodenterms, int **node_base, int **node_edge, STP_Bool firstrun, STP_Bool *connected)
Definition: heur_tm.c:1328
#define STP_PCSPG
Definition: grph.h:40
#define SCIPdebugMessage
Definition: pub_message.h:77
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
int number
Definition: misc_stp.h:43
static SCIP_RETCODE computeSteinerTreeDijkPcMw(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *dijkdist, int *result, int *dijkedge, int start, STP_Bool *connected)
Definition: heur_tm.c:886
SCIP_RETCODE SCIPsetHeurInit(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURINIT((*heurinit)))
Definition: scip_heur.c:184
GRAPH * SCIPprobdataGetGraph(SCIP_PROBDATA *probdata)
SCIP_LPSOLSTAT SCIPgetLPSolstat(SCIP *scip)
Definition: scip_lp.c:158
SCIP_Bool SCIPhasCurrentNodeLP(SCIP *scip)
Definition: scip_lp.c:73
int *RESTRICT mark
Definition: grph.h:70
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
SCIP_HEUR * SCIPfindHeur(SCIP *scip, const char *name)
Definition: scip_heur.c:248
SCIP_RETCODE SCIPcreateSol(SCIP *scip, SCIP_SOL **sol, SCIP_HEUR *heur)
Definition: scip_sol.c:319
#define HEUR_DESC
Definition: heur_tm.c:40
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPincludeHeurBasic(SCIP *scip, SCIP_HEUR **heur, const char *name, const char *desc, char dispchar, int priority, int freq, int freqofs, int maxdepth, SCIP_HEURTIMING timingmask, SCIP_Bool usessubscip, SCIP_DECL_HEUREXEC((*heurexec)), SCIP_HEURDATA *heurdata)
Definition: scip_heur.c:107
static SCIP_DECL_PARAMCHGD(paramChgdRandomseed)
Definition: heur_tm.c:97
void graph_path_execX(SCIP *, const GRAPH *, int, const SCIP_Real *, SCIP_Real *, int *)
Definition: grphpath.c:781
void graph_path_st_pcmw(SCIP *, const GRAPH *, const SCIP_Real *, SCIP_Real *, int *, int, STP_Bool *)
Definition: grphpath.c:1154
#define HEUR_FREQOFS
Definition: heur_tm.c:44
int *RESTRICT oeat
Definition: grph.h:103
#define CONNECT
Definition: grph.h:154
static SCIP_RETCODE prune(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, int *result, STP_Bool *connected)
Definition: heur_tm.c:831
#define STP_SAP
Definition: grph.h:39
int stp_type
Definition: grph.h:127
void heap_add(int *, int *, int *, int, PATH *)
Definition: grphpath.c:244
#define Is_pterm(a)
Definition: grph.h:169
unsigned char STP_Bool
Definition: grph.h:52
static SCIP_RETCODE computeSteinerTreeDijkPcMwFull(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, SCIP_Real *dijkdist, int *result, int *dijkedge, int start, STP_Bool *connected)
Definition: heur_tm.c:912
#define DEFAULT_DURINGLPFREQ
Definition: heur_tm.c:53
SCIP_RETCODE SCIPprobdataAddNewSol(SCIP *scip, SCIP_Real *nval, SCIP_SOL *sol, SCIP_HEUR *heur, SCIP_Bool *success)
SCIP_Real SCIPprobdataGetOffset(SCIP *scip)
#define SCIPallocBuffer(scip, ptr)
Definition: scip_mem.h:109
void SCIPStpHeurTMCompStarts(GRAPH *graph, int *starts, int *runs)
Definition: heur_tm.c:114
#define STP_DCSTP
Definition: grph.h:43
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:124
SCIP_Real dist
Definition: misc_stp.h:44
SCIP_Real * prize
Definition: grph.h:82
SCIP_Real dist
Definition: grph.h:145
int *RESTRICT grad
Definition: grph.h:73
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool graph_sol_valid(SCIP *, const GRAPH *, const int *)
Definition: grphbase.c:3066
void SCIPheurSetData(SCIP_HEUR *heur, SCIP_HEURDATA *heurdata)
Definition: heur.c:1175
internal miscellaneous methods
SCIPInterval fabs(const SCIPInterval &x)
void graph_path_exec(SCIP *, const GRAPH *, const int, int, const SCIP_Real *, PATH *)
Definition: grphpath.c:491
static SCIP_DECL_HEURCOPY(heurCopyTM)
Definition: heur_tm.c:2635
void graph_path_st_pcmw_full(SCIP *, const GRAPH *, const SCIP_Real *, SCIP_Real *, int *, int, STP_Bool *)
Definition: grphpath.c:1316
int knots
Definition: grph.h:62
#define SCIP_CALL(x)
Definition: def.h:365
#define DEFAULT_TYPE
Definition: heur_tm.c:54
#define AUTO
Definition: heur_tm.c:57
#define FARAWAY
Definition: grph.h:156
SCIP_RETCODE SCIPpqueueCreate(SCIP_PQUEUE **pqueue, int initsize, SCIP_Real sizefac, SCIP_DECL_SORTPTRCOMP((*ptrcomp)))
Definition: misc.c:1234
int GNODECmpByDist(void *first_arg, void *second_arg)
#define STP_SPG
Definition: grph.h:38
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
int SCIPrandomGetInt(SCIP_RANDNUMGEN *randnumgen, int minrandval, int maxrandval)
Definition: misc.c:9618
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:637
SCIP_EXPORT void SCIPsortRealInt(SCIP_Real *realarray, int *intarray, int len)
#define SCIP_Bool
Definition: def.h:70
int *RESTRICT ieat
Definition: grph.h:102
int *RESTRICT path_heap
Definition: grph.h:118
#define STP_MWCSP
Definition: grph.h:47
SCIP_EXPORT SCIP_Real SCIPvarGetUbGlobal(SCIP_VAR *var)
Definition: var.c:17362
int *RESTRICT tail
Definition: grph.h:95
#define MIN(x, y)
Definition: def.h:223
static SCIP_RETCODE computeSteinerTree(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *costrev, SCIP_Real **pathdist, int start, int *perm, int *result, int *cluster, int **pathedge, STP_Bool *connected, SCIP_RANDNUMGEN *randnumgen)
Definition: heur_tm.c:933
static const unsigned int randseed
Definition: circle.c:46
void SCIPrandomPermuteIntArray(SCIP_RANDNUMGEN *randnumgen, int *array, int begin, int end)
Definition: misc.c:9659
SCIP_RETCODE SCIPStpValidateSol(SCIP *, const GRAPH *, const double *, SCIP_Bool *)
Definition: validate.c:136
int *RESTRICT term
Definition: grph.h:68
SCIP_Real graph_sol_getObj(const SCIP_Real *, const int *, SCIP_Real, int)
Definition: grphbase.c:3196
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:124
static long * number
Portable defintions.
SCIP_RETCODE graph_voronoiExtend(SCIP *, const GRAPH *, SCIP_Real *, PATH *, SCIP_Real **, int **, int **, STP_Bool *, int *, int *, int *, int, int, int)
Definition: grphpath.c:1671
int SCIPparamGetInt(SCIP_PARAM *param)
Definition: paramset.c:716
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:67
SCIP_Real * r
Definition: circlepacking.c:50
SCIP_EXPORT SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17418
SCIP_EXPORT void SCIPsortRealIntInt(SCIP_Real *realarray, int *intarray1, int *intarray2, int len)
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:121
#define Is_term(a)
Definition: grph.h:168
struct SCIP_ProbData SCIP_PROBDATA
Definition: type_prob.h:44
SCIP_Longint SCIPgetNLPIterations(SCIP *scip)
static SCIP_RETCODE computeSteinerTreeDijk(SCIP *scip, const GRAPH *g, SCIP_Real *cost, SCIP_Real *dijkdist, int *result, int *dijkedge, int start, SCIP_RANDNUMGEN *randnumgen, STP_Bool *connected)
Definition: heur_tm.c:855
SCIP_Real * cost
Definition: grph.h:94
SCIP_RETCODE SCIPsetHeurCopy(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURCOPY((*heurcopy)))
Definition: scip_heur.c:152
SCIP_VAR * a
Definition: circlepacking.c:57
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10263
#define TM_DIJKSTRA
Definition: heur_tm.c:60
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:73
#define SCIP_HEURTIMING_AFTERPSEUDONODE
Definition: type_timing.h:76
#define SCIP_Real
Definition: def.h:164
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real * SCIPprobdataGetXval(SCIP *scip, SCIP_SOL *sol)
int *RESTRICT outbeg
Definition: grph.h:76
SCIP_RETCODE SCIPStpHeurTMPrunePc(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, int *result, STP_Bool *connected)
Definition: heur_tm.c:167
#define HEUR_PRIORITY
Definition: heur_tm.c:42
shortest paths based primal heuristics for Steiner problems
#define SCIP_Longint
Definition: def.h:149
int edges
Definition: grph.h:92
SCIP_RETCODE SCIPStpHeurTMBuildTree(SCIP *scip, const GRAPH *g, PATH *mst, const SCIP_Real *cost, SCIP_Real *objresult, int *connected)
Definition: heur_tm.c:653
#define flipedge(edge)
Definition: grph.h:150
SCIP_RETCODE SCIPStpHeurTMRunLP(SCIP *scip, GRAPH *graph, SCIP_HEUR *heur, int *result, int runs, SCIP_Real *cost, SCIP_Real *costrev, SCIP_Bool *success)
Definition: heur_tm.c:2352
SCIP_RETCODE SCIPsetHeurFree(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURFREE((*heurfree)))
Definition: scip_heur.c:168
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
#define HEUR_USESSUBSCIP
Definition: heur_tm.c:47
#define UNKNOWN
Definition: sepa_mcf.c:4081
#define nnodes
Definition: gastrans.c:65
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:120
#define TM_VORONOI
Definition: heur_tm.c:59
#define DEFAULT_EVALRUNS
Definition: heur_tm.c:49
#define DEFAULT_INITRUNS
Definition: heur_tm.c:50
int hoplimit
Definition: grph.h:89
#define HEUR_NAME
Definition: heur_tm.c:39
SCIP_RETCODE SCIPStpHeurTMPrune(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, int layer, int *result, STP_Bool *connected)
Definition: heur_tm.c:555
#define STP_RPCSPG
Definition: grph.h:41
#define TM_SP
Definition: heur_tm.c:58
#define HEUR_MAXDEPTH
Definition: heur_tm.c:45
#define DEFAULT_ROOTRUNS
Definition: heur_tm.c:52
void graph_pc_2transcheck(GRAPH *)
Definition: grphbase.c:1041
#define DEFAULT_RANDSEED
Definition: heur_tm.c:55
SCIP_RETCODE SCIPStpHeurTMBuildTreePcMw(SCIP *scip, const GRAPH *g, PATH *mst, const SCIP_Real *cost, SCIP_Real *objresult, int *connected)
Definition: heur_tm.c:382
SCIP_RETCODE SCIPpqueueInsert(SCIP_PQUEUE *pqueue, void *elem)
Definition: misc.c:1280