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