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