Scippy

SCIP

Solving Constraint Integer Programs

heur_rec.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-2018 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_rec.c
17  * @brief Primal recombination heuristic for Steiner problems
18  * @author Daniel Rehfeldt
19  *
20  * This file implements a recombination heuristic for Steiner problems, see
21  * "SCIP-Jack - A solver for STP and variants with parallelization extensions" (2017) by
22  * Gamrath, Koch, Maher, Rehfeldt and Shinano
23  *
24  * A list of all interface methods can be found in heur_rec.h
25  *
26  */
27 
28 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
29 
30 #include <assert.h>
31 #include <string.h>
32 #include <stdio.h>
33 #include "scip/scip.h"
34 #include "scip/scipdefplugins.h"
35 #include "scip/cons_linear.h"
36 #include "heur_rec.h"
37 #include "heur_prune.h"
38 #include "heur_slackprune.h"
39 #include "heur_local.h"
40 #include "grph.h"
41 #include "heur_tm.h"
42 #include "cons_stp.h"
43 #include "scip/pub_misc.h"
44 #include "probdata_stp.h"
45 #include "math.h"
46 
47 #define HEUR_NAME "rec"
48 #define HEUR_DESC "recombination heuristic for Steiner problems"
49 #define HEUR_DISPCHAR 'R'
50 #define HEUR_PRIORITY 100
51 #define HEUR_FREQ 1
52 #define HEUR_FREQOFS 0
53 #define HEUR_MAXDEPTH -1
54 #define HEUR_TIMING (SCIP_HEURTIMING_DURINGLPLOOP | SCIP_HEURTIMING_AFTERLPLOOP | SCIP_HEURTIMING_AFTERNODE)
55 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
56 
57 #define DEFAULT_MAXFREQREC FALSE /**< executions of the rec heuristic at maximum frequency? */
58 #define DEFAULT_MAXNSOLS 15 /**< default maximum number of (good) solutions be regarded in the subproblem */
59 #define DEFAULT_NUSEDSOLS 4 /**< number of solutions that will be taken into account */
60 #define DEFAULT_RANDSEED 1984 /**< random seed */
61 #define DEFAULT_NTMRUNS 100 /**< number of runs in TM heuristic */
62 #define DEFAULT_NWAITINGSOLS 4 /**< max number of new solutions to be available before executing the heuristic again */
63 
64 #define BOUND_MAXNTERMINALS 1000
65 #define BOUND_MAXNEDGES 20000
66 #define RUNS_RESTRICTED 3
67 #define RUNS_NORMAL 10
68 #define CYCLES_PER_RUN 3
69 #define REC_MAX_FAILS 4
70 #define REC_MIN_NSOLS 4
71 #define REC_ADDEDGE_FACTOR 1.0
72 
73 #define COST_MAX_POLY_x0 500
74 #define COST_MIN_POLY_x0 100
75 
76 #ifdef WITH_UG
77 extern
78 int getUgRank(void);
79 #endif
80 
81 
82 /*
83  * Data structures
84  */
85 
86 /** primal heuristic data */
87 struct SCIP_HeurData
88 {
89  int lastsolindex;
90  int bestsolindex; /**< best solution during the previous run */
91  int maxnsols; /**< maximum number of (good) solutions be regarded in the subproblem */
92  SCIP_Longint ncalls; /**< number of calls */
93  SCIP_Longint nlastsols; /**< number of solutions during the last run */
94  int ntmruns; /**< number of runs in TM heuristic */
95  int nusedsols; /**< number of solutions that will be taken into account */
96  int nselectedsols; /**< number of solutions actually selected */
97  int nwaitingsols; /**< number of new solutions before executing the heuristic again */
98  int nfailures; /**< number of failures since last successful call */
99  unsigned int randseed; /**< seed value for random number generator */
100  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
101  SCIP_Bool maxfreq; /**< should the heuristic be called at maximum frequency? */
102 };
103 
104 
105 /*
106  * Local methods
107  */
108 
109 /** information method for a parameter change of random seed */
110 static
111 SCIP_DECL_PARAMCHGD(paramChgdRandomseed)
112 { /*lint --e{715}*/
113  SCIP_HEURDATA* heurdata;
114  int newrandseed;
115 
116  newrandseed = SCIPparamGetInt(param);
117 
118  heurdata = (SCIP_HEURDATA*)SCIPparamGetData(param);
119  assert(heurdata != NULL);
120 
121  heurdata->randseed = (int) newrandseed;
122 
123  return SCIP_OKAY;
124 }
125 
126 
127 /** edge cost multiplier */
128 static
130  SCIP* scip, /**< SCIP data structure */
131  SCIP_HEURDATA* heurdata, /**< SCIP data structure */
132  SCIP_Real avg /**< number of solutions containing this edge */
133  )
134 {
135  SCIP_Real normed;
136  int maxpolyx0;
137  int minpolyx0;
138 
139  int upper;
140  int lower;
141  int factor;
142 
143  /* if STP, then avg >= 2.0 */
144  assert(SCIPisGE(scip, avg, 1.0));
145  assert(heurdata->nusedsols >= 2);
146 
147  maxpolyx0 = COST_MAX_POLY_x0 * (heurdata->nusedsols - 1);
148  minpolyx0 = COST_MIN_POLY_x0 * (heurdata->nusedsols - 1);
149 
150  /* norm to [0,1] (STP) or [-1,1] and evaluate polynomials */
151  normed = (avg - 2.0) / ((double) heurdata->nusedsols - 1.0);
152  upper = (int) (maxpolyx0 - 2 * maxpolyx0 * normed + maxpolyx0 * (normed * normed));
153  lower = (int) (minpolyx0 - 2 * minpolyx0 * normed + minpolyx0 * (normed * normed));
154 
155  assert(SCIPisGE(scip, normed, -1.0) && SCIPisLE(scip, normed, 1.0));
156 
157  lower = MAX(0, lower);
158  upper = MAX(0, upper);
159 
160  factor = SCIPrandomGetInt(heurdata->randnumgen, lower, upper);
161  factor++;
162 
163  assert(factor <= COST_MAX_POLY_x0 * 3 + 1 || avg <= 2.0);
164  assert(factor >= 1);
165 
166  return factor;
167 }
168 
169 /** select solutions to be merged */
170 static
172  SCIP* scip, /**< SCIP data structure */
173  const STPSOLPOOL* pool, /**< solution pool or NULL */
174  const GRAPH* graph, /**< graph data structure */
175  SCIP_HEURDATA* heurdata, /**< SCIP data structure */
176  SCIP_VAR** vars, /**< problem variables */
177  const int* incumbentedges, /**< edges of solution to be used as recombination root */
178  int* incumbentindex, /**< index of ancestor of incumbent solution */
179  int* selection, /**< selected solutions */
180  SCIP_Bool* success /**< could at least two solutions be selected? */
181  )
182 {
183  SCIP_SOL** sols = NULL;
184  STPSOL** poolsols = NULL;
185  int* perm;
186  int* soledgestmp;
187  STP_Bool* solselected;
188  STP_Bool* soledges;
189  int pos;
190  int nsols;
191  int maxnsols = heurdata->maxnsols;
192  int nselectedsols = 0;
193  const int nedges = graph->edges;
194  const int nusedsols = heurdata->nusedsols;
195  const SCIP_Bool usestppool = (pool != NULL);
196 
197  assert(selection != NULL);
198  assert(graph != NULL);
199 
200  if( !usestppool )
201  {
202  assert(vars != NULL);
203 
204  /* get solution data */
205  sols = SCIPgetSols(scip);
206  nsols = SCIPgetNSols(scip);
207  }
208  else
209  {
210  poolsols = pool->sols;
211  nsols = pool->size;
212  assert(nsols > 1);
213  }
214 
215  assert(nusedsols > 1);
216  assert(nsols >= nusedsols);
217 
218  /* allocate memory */
219  SCIP_CALL( SCIPallocBufferArray(scip, &solselected, nsols) );
220  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nsols) );
221  SCIP_CALL( SCIPallocBufferArray(scip, &soledges, nedges / 2) );
222  SCIP_CALL( SCIPallocBufferArray(scip, &soledgestmp, nedges / 2) );
223 
224  for( int i = 0; i < nsols; i++ )
225  {
226  perm[i] = i;
227  solselected[i] = FALSE;
228  }
229 
230  if( usestppool )
231  {
232  for( pos = 0; pos < nsols; pos++ )
233  if( *incumbentindex == poolsols[pos]->index )
234  break;
235  }
236  else
237  {
238  for( pos = 0; pos < nsols; pos++ )
239  if( *incumbentindex == SCIPsolGetIndex(sols[pos]) )
240  break;
241  }
242  assert(pos < nsols);
243  solselected[pos] = TRUE;
244  selection[nselectedsols++] = pos;
245 
246  for( int e = 0; e < nedges; e += 2 )
247  soledges[e / 2] = (incumbentedges[e] == CONNECT
248  || incumbentedges[e + 1] == CONNECT);
249 
250  maxnsols = MIN(nsols, maxnsols);
251 
252  if( !usestppool )
253  {
254  int sqrtnallsols = (int) sqrt((double) nsols);
255 
256  if( sqrtnallsols >= REC_MIN_NSOLS && sqrtnallsols < maxnsols )
257  maxnsols = sqrtnallsols;
258  }
259 
260  SCIPrandomPermuteIntArray(heurdata->randnumgen, perm, 0, maxnsols);
261 
262  for( int i = 0; i < nsols; i++ )
263  {
264  if( solselected[perm[i]] == FALSE )
265  {
266  int eqnedges = 0;
267  int diffnedges = 0;
268  const int k = perm[i];
269  SCIP_SOL* solk = NULL;
270  const int* solKedges = NULL;
271 
272  if( usestppool )
273  solKedges = poolsols[k]->soledges;
274  else
275  solk = sols[k];
276 
277  for( int e = 0; e < nedges; e += 2 )
278  {
279  SCIP_Bool hit;
280 
281  if( usestppool )
282  hit = (solKedges[e] == CONNECT || solKedges[e + 1] == CONNECT);
283  else
284  hit = SCIPisEQ(scip, SCIPgetSolVal(scip, solk, vars[e]), 1.0)
285  || SCIPisEQ(scip, SCIPgetSolVal(scip, solk, vars[e + 1]), 1.0);
286 
287  if( hit )
288  {
289  if( soledges[e / 2] == FALSE )
290  soledgestmp[diffnedges++] = e / 2;
291  else
292  {
293  const int tail = graph->tail[e];
294  const int head = graph->head[e];
295 
296  /* no dummy edge? */
297  if( !(Is_term(graph->term[tail]) && Is_term(graph->term[head])) )
298  eqnedges++;
299  }
300  }
301  }
302 
303  /* enough similarities and differences with new solution? */
304  if( diffnedges > 5 && eqnedges > 0 )
305  {
306  SCIPdebugMessage("REC: different edges: %d same edges: %d \n", diffnedges, eqnedges);
307  selection[nselectedsols++] = k;
308  solselected[k] = TRUE;
309  *success = TRUE;
310 
311  for( int j = 0; j < diffnedges; j++ )
312  soledges[soledgestmp[j]] = TRUE;
313 
314  if( nselectedsols >= nusedsols )
315  break;
316  }
317  }
318  }
319 
320  assert(nselectedsols <= nusedsols);
321 
322  heurdata->nselectedsols = nselectedsols;
323 
324  SCIPdebugMessage("REC: selected %d sols \n", nselectedsols);
325 
326  /* free memory */
327  SCIPfreeBufferArray(scip, &soledgestmp);
328  SCIPfreeBufferArray(scip, &soledges);
329  SCIPfreeBufferArray(scip, &perm);
330  SCIPfreeBufferArray(scip, &solselected);
331 
332  return SCIP_OKAY;
333 }
334 
335 /** select solutions to be merged */
336 static
338  SCIP* scip, /**< SCIP data structure */
339  const STPSOLPOOL* pool, /**< solution pool or NULL */
340  SCIP_HEURDATA* heurdata, /**< SCIP data structure */
341  int* incumbentindex, /**< index of ancestor of incumbent solution */
342  int* selection, /**< selected solutions */
343  SCIP_Bool randomize /**< select solutions randomly? */
344  )
345 {
346  SCIP_SOL** sols = NULL; /* array of all solutions found so far */
347  STPSOL** poolsols = NULL;
348  int* perm;
349  STP_Bool* solselected;
350  int i;
351  int nsols; /* number of all available solutions */
352  int maxnsols;
353  int nusedsols; /* number of solutions to use in rec */
354  int nselectedsols;
355 
356  const SCIP_Bool usestppool = (pool != NULL);
357  assert(selection != NULL);
358 
359  if( usestppool )
360  {
361  poolsols = pool->sols;
362  nsols = pool->size;
363  assert(nsols > 1);
364  }
365  else
366  {
367  /* get solution data */
368  sols = SCIPgetSols(scip);
369  nsols = SCIPgetNSols(scip);
370  }
371 
372  maxnsols = heurdata->maxnsols;
373  nusedsols = heurdata->nusedsols;
374  assert(nusedsols <= nsols);
375  nselectedsols = 0;
376 
377  assert(nusedsols > 1);
378  assert(nsols >= nusedsols);
379 
380  /* allocate memory */
381  SCIP_CALL( SCIPallocBufferArray(scip, &solselected, nsols) );
382  SCIP_CALL( SCIPallocBufferArray(scip, &perm, nsols) );
383 
384  for( i = 0; i < nsols; i++ )
385  {
386  perm[i] = i;
387  solselected[i] = FALSE;
388  }
389 
390  if( usestppool )
391  {
392  for( i = 0; i < nsols; i++ )
393  if( *incumbentindex == poolsols[i]->index )
394  break;
395  }
396  else
397  {
398  for( i = 0; i < nsols; i++ )
399  if( *incumbentindex == SCIPsolGetIndex(sols[i]) )
400  break;
401  }
402 
403  assert(i < nsols);
404  solselected[i] = TRUE;
405  selection[nselectedsols++] = i;
406 
407  if( !randomize )
408  {
409  const int end = SCIPrandomGetInt(heurdata->randnumgen, 1, nusedsols - 1);
410  int shift = SCIPrandomGetInt(heurdata->randnumgen, end, 2 * nusedsols - 1);
411 
412  if( shift > nsols )
413  shift = nsols;
414 
415  SCIPrandomPermuteIntArray(heurdata->randnumgen, perm, 0, shift);
416 
417  for( i = 0; i < end; i++ )
418  {
419  if( solselected[perm[i]] == FALSE )
420  {
421  selection[nselectedsols++] = perm[i];
422  solselected[perm[i]] = TRUE;
423  }
424  }
425  }
426 
427  maxnsols = MIN(nsols, maxnsols);
428 
429  if( !usestppool )
430  {
431  int sqrtnallsols = (int) sqrt(nsols);
432 
433  if( sqrtnallsols >= REC_MIN_NSOLS && sqrtnallsols < maxnsols )
434  maxnsols = sqrtnallsols;
435  }
436 
437  SCIPdebugMessage("maxnsols in REC %d \n", maxnsols);
438 
439  if( nselectedsols < nusedsols )
440  {
441  SCIPrandomPermuteIntArray(heurdata->randnumgen, perm, 0, maxnsols);
442  for( i = 0; i < maxnsols; i++ )
443  {
444  if( solselected[perm[i]] == FALSE )
445  {
446  selection[nselectedsols++] = perm[i];
447  if( nselectedsols >= nusedsols )
448  break;
449  }
450  }
451  }
452 
453  assert(nselectedsols <= nusedsols);
454 
455  SCIPdebugMessage("REC: selected %d sols \n", nselectedsols);
456 
457  heurdata->nselectedsols = nselectedsols;
458  SCIPfreeBufferArray(scip, &perm);
459  SCIPfreeBufferArray(scip, &solselected);
460 
461  return SCIP_OKAY;
462 }
463 
464 /** merge selected solutions to a new graph */
465 static
467  SCIP* scip, /**< SCIP data structure */
468  const STPSOLPOOL* pool, /**< solution pool or NULL */
469  SCIP_HEURDATA* heurdata, /**< SCIP data structure */
470  const GRAPH* graph, /**< graph structure */
471  GRAPH** solgraph, /**< pointer to store new graph */
472  const int* incumbentedges, /**< edges of solution to be used as recombination root */
473  int* incumbentindex, /**< index of ancestor of incumbent solution */
474  int** edgeancestor, /**< ancestor to edge edge */
475  int** edgeweight, /**< for each edge: number of solution that contain this edge */
476  SCIP_Bool* success, /**< new graph constructed? */
477  SCIP_Bool randomize, /**< select solution randomly? */
478  SCIP_Bool addedges /**< add additional edges between solution vertices? */
479  )
480 {
481  GRAPH* newgraph = NULL;
482  SCIP_SOL** sols = NULL;
483  STPSOL** poolsols = NULL;
484  SCIP_VAR** vars = NULL;
485  int* solselection;
486  const SCIP_Bool pcmw = (graph->stp_type == STP_PCSPG || graph->stp_type == STP_MWCSP
487  || graph->stp_type == STP_RPCSPG || graph->stp_type == STP_RMWCSP );
488  const SCIP_Bool usestppool = (pool != NULL);
489 
490  assert(scip != NULL);
491  assert(graph != NULL);
492 
493  if( !usestppool )
494  {
495  sols = SCIPgetSols(scip);
496  vars = SCIPprobdataGetEdgeVars(scip);
497  assert(vars != NULL);
498  assert(sols != NULL);
499  }
500  else
501  {
502  poolsols = pool->sols;
503  }
504 
505  *success = TRUE;
506  *edgeweight = NULL;
507  *edgeancestor = NULL;
508 
509  /* allocate memory */
510  SCIP_CALL( SCIPallocBufferArray(scip, &solselection, heurdata->nusedsols) );
511 
512  /* select solutions to be merged */
513  if( pcmw || graph->stp_type == STP_DCSTP )
514  SCIP_CALL( selectdiffsols(scip, pool, graph, heurdata, vars, incumbentedges, incumbentindex, solselection, success) );
515  else
516  SCIP_CALL( selectsols(scip, pool, heurdata, incumbentindex, solselection, randomize) );
517 
518  if( *success )
519  {
520  int* dnodemap;
521  STP_Bool* solnode; /* marks nodes contained in at least one solution */
522  STP_Bool* soledge; /* marks edges contained in at least one solution */
523  int j;
524  int nsoledges = 0;
525  int nsolnodes = 0;
526  const int nedges = graph->edges;
527  const int nnodes = graph->knots;
528  const int selectedsols = heurdata->nselectedsols;
529 
530  SCIPdebugMessage("REC: successfully selected new solutions \n");
531 
532  assert(selectedsols > 0);
533 
534  SCIP_CALL( SCIPallocBufferArray(scip, &solnode, nnodes) );
535  SCIP_CALL( SCIPallocBufferArray(scip, &dnodemap, nnodes) );
536  SCIP_CALL( SCIPallocBufferArray(scip, &soledge, nedges / 2) );
537 
538  for( int i = 0; i < nnodes; i++ )
539  {
540  solnode[i] = FALSE;
541  dnodemap[i] = UNKNOWN;
542  }
543 
544  /* count and mark selected nodes and edges */
545  for( int i = 0; i < nedges; i += 2 )
546  {
547  const int ihalf = i / 2;
548  soledge[ihalf] = FALSE;
549 
550  if( graph->oeat[i] == EAT_FREE )
551  continue;
552 
553  for( j = 0; j < selectedsols; j++ )
554  {
555  SCIP_Bool hit;
556 
557  if( j == 0 )
558  {
559  hit = (incumbentedges[i] == CONNECT
560  || incumbentedges[i + 1] == CONNECT);
561  }
562  else if( usestppool )
563  hit = (poolsols[solselection[j]]->soledges[i] == CONNECT
564  || poolsols[solselection[j]]->soledges[i + 1] == CONNECT);
565  else
566  hit = (SCIPisEQ(scip, SCIPgetSolVal(scip, sols[solselection[j]], vars[i]), 1.0)
567  || SCIPisEQ(scip, SCIPgetSolVal(scip, sols[solselection[j]], vars[i + 1]), 1.0));
568 
569  if( hit )
570  {
571  nsoledges++;
572  soledge[ihalf] = TRUE;
573  if( !solnode[graph->tail[i]] )
574  {
575  solnode[graph->tail[i]] = TRUE;
576  nsolnodes++;
577  }
578  if( !solnode[graph->head[i]] )
579  {
580  solnode[graph->head[i]] = TRUE;
581  nsolnodes++;
582  }
583  break;
584  }
585  }
586  }
587  if( pcmw ) /* todo this probably won't work for RMWCSP */
588  {
589  const int oldroot = graph->source;
590  for( int i = graph->outbeg[oldroot]; i != EAT_LAST; i = graph->oeat[i] )
591  {
592  if( Is_gterm(graph->term[graph->head[i]]) )
593  {
594  const int ihalf = i / 2;
595  const int head = graph->head[i];
596  if( soledge[ihalf] == FALSE )
597  {
598  nsoledges++;
599  soledge[ihalf] = TRUE;
600  if( !solnode[head] && SCIPisEQ(scip, graph->cost[flipedge(i)], FARAWAY) )
601  {
602  solnode[head] = TRUE;
603  nsolnodes++;
604  }
605  assert(solnode[graph->head[i]]);
606  }
607 
608  if( Is_pterm(graph->term[head]) )
609  {
610  int e2;
611  for( e2 = graph->outbeg[head]; e2 != EAT_LAST; e2 = graph->oeat[e2] )
612  if( Is_term(graph->term[graph->head[e2]]) && graph->head[e2] != oldroot )
613  break;
614 
615  assert(e2 != EAT_LAST);
616 
617  if( soledge[e2 / 2] == FALSE )
618  {
619  nsoledges++;
620  soledge[e2 / 2] = TRUE;
621  }
622  }
623  else
624  {
625  int e2;
626  assert(Is_term(graph->term[head]));
627  for( e2 = graph->outbeg[head]; e2 != EAT_LAST; e2 = graph->oeat[e2] )
628  if( Is_pterm(graph->term[graph->head[e2]]) && graph->head[e2] != oldroot )
629  break;
630 
631  assert(e2 != EAT_LAST);
632 
633  if( soledge[e2 / 2] == FALSE )
634  {
635  nsoledges++;
636  soledge[e2 / 2] = TRUE;
637  }
638  }
639  }
640  }
641  }
642 
643  /* add additional edges? */
644  if( addedges )
645  {
646  int naddedges = 0;
647 
648  for( int e = 0; e < nedges && naddedges <= (int)(REC_ADDEDGE_FACTOR * nsoledges); e += 2 )
649  {
650  int tail;
651  int head;
652 
653  if( soledge[e / 2] )
654  continue;
655  // todo if fixed to zero, continue?
656  if( graph->oeat[e] == EAT_FREE )
657  continue;
658 
659  tail = graph->tail[e];
660  head = graph->head[e];
661 
662  if( solnode[tail] && solnode[head] )
663  {
664  soledge[e / 2] = TRUE;
665  naddedges++;
666  }
667  }
668  SCIPdebugMessage("additional edges %d (orig: %d )\n", naddedges, nsoledges);
669 
670  nsoledges += naddedges;
671  }
672 
673  if( graph->stp_type == STP_GSTP )
674  {
675  for( int k = 0; k < nnodes; k++ )
676  {
677  if( Is_term(graph->term[k]) )
678  {
679  assert(solnode[k]);
680  for( int i = graph->outbeg[k]; i != EAT_LAST; i = graph->oeat[i] )
681  {
682  if( solnode[graph->head[i]] && !soledge[i / 2] )
683  {
684  soledge[i / 2] = TRUE;
685  nsoledges++;
686  }
687  }
688  }
689  }
690  }
691 
692  /* initialize new graph */
693  SCIP_CALL( graph_init(scip, &newgraph, nsolnodes, 2 * nsoledges, 1) );
694 
695  if( graph->stp_type == STP_RSMT || graph->stp_type == STP_OARSMT || graph->stp_type == STP_GSTP )
696  newgraph->stp_type = STP_SPG;
697  else
698  newgraph->stp_type = graph->stp_type;
699 
700  if( pcmw )
701  {
702  SCIP_CALL( graph_pc_init(scip, newgraph, nsolnodes, nsolnodes) );
703  }
704 
705  newgraph->hoplimit = graph->hoplimit;
706  j = 0;
707  for( int i = 0; i < nnodes; i++ )
708  {
709  if( solnode[i] )
710  {
711  if( pcmw )
712  {
713  if( (!Is_term(graph->term[i])) )
714  newgraph->prize[j] = graph->prize[i];
715  else
716  newgraph->prize[j] = 0.0;
717  }
718 
719  graph_knot_add(newgraph, graph->term[i]);
720  dnodemap[i] = j++;
721  }
722  }
723 
724  if( pcmw )
725  {
726  newgraph->extended = TRUE;
727  newgraph->norgmodelknots = newgraph->knots - newgraph->terms;
728  }
729 
730  SCIPdebugMessage("REC: sol graph with nodes: %d, edges: %d, terminals: %d \n", nsolnodes, 2 * nsoledges, newgraph->terms);
731 
732  /* set root */
733  newgraph->source = dnodemap[graph->source];
734  if( newgraph->stp_type == STP_RPCSPG || newgraph->stp_type == STP_RMWCSP )
735  newgraph->prize[newgraph->source] = FARAWAY;
736 
737  assert(newgraph->source >= 0);
738 
739  /* copy max degrees*/
740  if( graph->stp_type == STP_DCSTP )
741  {
742  SCIP_CALL( SCIPallocMemoryArray(scip, &(newgraph->maxdeg), nsolnodes) );
743  for( int i = 0; i < nnodes; i++ )
744  if( solnode[i] )
745  newgraph->maxdeg[dnodemap[i]] = graph->maxdeg[i];
746  }
747  /* allocate memory */
748  SCIP_CALL( SCIPallocMemoryArray(scip, edgeancestor, 2 * nsoledges) );
749  SCIP_CALL( SCIPallocMemoryArray(scip, edgeweight, 2 * nsoledges) );
750 
751  for( int i = 0; i < 2 * nsoledges; i++ )
752  (*edgeweight)[i] = 1;
753 
754  /* store original ID of each new edge (i.e. edge in the merged graph) */
755  j = 0;
756  assert(selectedsols == heurdata->nselectedsols);
757  for( int i = 0; i < nedges; i += 2 )
758  {
759  if( soledge[i / 2] )
760  {
761  const int orgtail = graph->tail[i];
762  const int orghead = graph->head[i];
763 
764  (*edgeancestor)[j++] = i;
765  (*edgeancestor)[j++] = i + 1;
766 
767  if( pcmw )
768  {
769  assert(newgraph->term2edge != NULL);
770  graph_pc_updateTerm2edge(newgraph, graph, dnodemap[orgtail], dnodemap[orghead], orgtail, orghead);
771  }
772 
773  graph_edge_add(scip, newgraph, dnodemap[orgtail], dnodemap[orghead], graph->cost[i], graph->cost[i + 1]);
774 
775  /* (*edgeweight)[e]: number of solutions containing edge e */
776  for( int k = 0; k < selectedsols; k++ )
777  {
778  SCIP_Bool hit;
779 
780  if( k == 0 )
781  {
782  hit = (incumbentedges[i] == CONNECT
783  || incumbentedges[i + 1] == CONNECT);
784  }
785  else if( usestppool )
786  hit = (poolsols[solselection[k]]->soledges[i] == CONNECT
787  || poolsols[solselection[k]]->soledges[i + 1] == CONNECT);
788  else
789  hit = (SCIPisEQ(scip, SCIPgetSolVal(scip, sols[solselection[k]], vars[i]), 1.0)
790  || SCIPisEQ(scip, SCIPgetSolVal(scip, sols[solselection[k]], vars[i + 1]), 1.0));
791 
792  /* is edge i or reversed edge in current solution? */
793  if( hit )
794  {
795  (*edgeweight)[j - 2]++;
796  (*edgeweight)[j - 1]++;
797  }
798  }
799  }
800  }
801 
802  assert(j == 2 * nsoledges);
803  SCIPfreeBufferArray(scip, &soledge);
804  SCIPfreeBufferArray(scip, &dnodemap);
805  SCIPfreeBufferArray(scip, &solnode);
806  }
807 
808  SCIPfreeBufferArray(scip, &solselection);
809  assert(graph_valid(newgraph));
810  *solgraph = newgraph;
811 
812  return SCIP_OKAY;
813 }
814 
815 static
817  GRAPH* g,
818  IDX* curr,
819  int* unodemap,
820  STP_Bool* stvertex
821  )
822 {
823  while( curr != NULL )
824  {
825  int i = curr->index;
826 
827  stvertex[unodemap[ g->orghead[i] ]] = TRUE;
828  stvertex[unodemap[ g->orgtail[i] ]] = TRUE;
829 
830  curr = curr->parent;
831  }
832 }
833 
834 static
836  const int* soledges, /**< edge array of solution to be checked */
837  const STPSOLPOOL* pool /**< the pool */
838 )
839 {
840  STPSOL** poolsols = pool->sols;
841  const int poolsize = pool->size;
842  const int nedges = pool->nedges;
843 
844  for( int i = 0; i < poolsize; i++ )
845  {
846  int j;
847  const int* pooledges = poolsols[i]->soledges;
848  assert(pooledges != NULL);
849 
850  for( j = 0; j < nedges; j++ )
851  if( pooledges[j] != soledges[j] )
852  break;
853 
854  /* pooledges == soledges? */
855  if( j == nedges )
856  {
857  SCIPdebugMessage("Pool: solution is already contained \n");
858  return TRUE;
859  }
860  }
861  return FALSE;
862 }
863 
864 /*
865  * primal heuristic specific interface methods
866  */
867 
868 
869 /** get solution from index */
871  STPSOLPOOL* pool, /**< the pool */
872  const int soindex /**< the index */
873  )
874 {
875  int i;
876  int size;
877 
878  assert(pool != NULL);
879  assert(soindex >= 0 && soindex <= pool->maxindex);
880 
881  size = pool->size;
882 
883  for( i = 0; i < size; i++ )
884  if( pool->sols[i]->index == soindex )
885  break;
886 
887  if( i == size )
888  return NULL;
889  else
890  return pool->sols[i];
891 }
892 
893 /** initializes STPSOL pool */
895  SCIP* scip, /**< SCIP data structure */
896  STPSOLPOOL** pool, /**< the pool */
897  const int nedges, /**< number of edges of solutions to be stored in the pool */
898  const int maxsize /**< capacity of pool */
899  )
900 {
901  STPSOLPOOL* dpool;
902 
903  assert(pool != NULL);
904  assert(nedges > 0);
905  assert(maxsize > 0);
906 
907  SCIP_CALL( SCIPallocBlockMemory(scip, pool) );
908 
909  dpool = *pool;
910  SCIP_CALL( SCIPallocMemoryArray(scip, &(dpool->sols), maxsize) );
911 
912  for( int i = 0; i < maxsize; i++ )
913  dpool->sols[i] = NULL;
914 
915  dpool->size = 0;
916  dpool->nedges = nedges;
917  dpool->maxsize = maxsize;
918  dpool->maxindex = -1;
919 
920  return SCIP_OKAY;
921 }
922 
923 
924 /** frees STPSOL pool */
926  SCIP* scip, /**< SCIP data structure */
927  STPSOLPOOL** pool /**< the pool */
928  )
929 {
930  STPSOLPOOL* dpool = *pool;
931  const int poolsize = dpool->size;
932 
933  assert(pool != NULL);
934  assert(dpool != NULL);
935  assert(poolsize == dpool->maxsize || dpool->sols[poolsize] == NULL);
936 
937  for( int i = poolsize - 1; i >= 0; i-- )
938  {
939  STPSOL* sol = dpool->sols[i];
940 
941  assert(sol != NULL);
942 
943  SCIPfreeMemoryArray(scip, &(sol->soledges));
944  SCIPfreeBlockMemory(scip, &sol);
945  }
946 
947  SCIPfreeMemoryArray(scip, &(dpool->sols));
948  SCIPfreeBlockMemory(scip, pool);
949 }
950 
951 
952 /** tries to add STPSOL to pool */
954  SCIP* scip, /**< SCIP data structure */
955  const SCIP_Real obj, /**< objective of solution to be added */
956  const int* soledges, /**< edge array of solution to be added */
957  STPSOLPOOL* pool, /**< the pool */
958  SCIP_Bool* success /**< has solution been added? */
959  )
960 {
961  STPSOL** poolsols = pool->sols;
962  STPSOL* sol;
963  int i;
964  int poolsize = pool->size;
965  const int nedges = pool->nedges;
966  const int poolmaxsize = pool->maxsize;
967 
968  assert(scip != NULL);
969  assert(pool != NULL);
970  assert(poolsols != NULL);
971  assert(poolsize >= 0);
972  assert(poolmaxsize >= 0);
973  assert(poolsize <= poolmaxsize);
974 
975  *success = FALSE;
976 
977  /* is solution in pool? */
978  if( isInPool(soledges, pool) )
979  return SCIP_OKAY;
980 
981  SCIPdebugMessage("Pool: add to pool (current size: %d, max: %d) \n", poolsize, poolmaxsize);
982 
983  /* enlarge pool if possible */
984  if( poolsize < poolmaxsize )
985  {
986  SCIPdebugMessage("Pool: alloc memory at position %d \n", poolsize);
987 
988  SCIP_CALL( SCIPallocBlockMemory(scip, &(poolsols[poolsize])) );
989  SCIP_CALL( SCIPallocMemoryArray(scip, &(poolsols[poolsize]->soledges), nedges) );
990 
991  poolsize++;
992  pool->size++;
993  }
994  /* pool is full; new solution worse than worst solution in pool? */
995  else if( SCIPisGT(scip, obj, poolsols[poolsize - 1]->obj) )
996  {
997  return SCIP_OKAY;
998  }
999 
1000  /* overwrite last element of pool (either empty or inferior to current solution) */
1001  sol = poolsols[poolsize - 1];
1002  assert(sol != NULL);
1003  sol->obj = obj;
1004  sol->index = ++(pool->maxindex);
1005  BMScopyMemoryArray(sol->soledges, soledges, nedges);
1006 
1007  /* shift solution up */
1008  for( i = poolsize - 1; i >= 1; i-- )
1009  {
1010  if( SCIPisGT(scip, obj, poolsols[i - 1]->obj) )
1011  break;
1012 
1013  poolsols[i] = poolsols[i - 1];
1014  }
1015 
1016  poolsols[i] = sol;
1017  SCIPdebugMessage("Pool: put new solution to position %d \n", i);
1018  *success = TRUE;
1019 
1020  return SCIP_OKAY;
1021 }
1022 
1023 
1024 
1025 /** runs STP recombination heuristic */
1027  SCIP* scip, /**< SCIP data structure */
1028  STPSOLPOOL* pool, /**< STP solution pool or NULL */
1029  SCIP_HEUR* heur, /**< heuristic or NULL */
1030  SCIP_HEURDATA* heurdata, /**< heuristic data or NULL */
1031  const GRAPH* graph, /**< graph data */
1032  SCIP_VAR** vars, /**< variables or NULL */
1033  int* newsolindex, /**< index of new solution */
1034  int runs, /**< number of runs */
1035  int nsols, /**< number of solutions in pool (SCIP or STP) */
1036  SCIP_Bool restrictheur, /**< use restricted version of heur? */
1037  SCIP_Bool* solfound /**< new solution found? */
1038 )
1039 {
1040  int* newsoledges;
1041  int* incumbentedges;
1042  SCIP_Real incumentobj;
1043  SCIP_Real hopfactor = 0.1;
1044 
1045  const STP_Bool usestppool = (pool != NULL);
1046  const int nnodes = graph->knots;
1047  const int nedges = graph->edges;
1048  const int probtype = graph->stp_type;
1049  const SCIP_Bool pcmw = (probtype == STP_PCSPG || probtype == STP_MWCSP || probtype == STP_RPCSPG || probtype == STP_RMWCSP );
1050  STP_Bool* stnodes;
1051 
1052  assert(runs >= 0);
1053  assert(graph != NULL);
1054  assert(scip != NULL);
1055  assert(*newsolindex >= 0);
1056 
1057  SCIPdebugMessage("REC: start \n");
1058 
1059  *solfound = FALSE;
1060 
1061  SCIP_CALL( SCIPallocBufferArray(scip, &newsoledges, nedges) );
1062  SCIP_CALL( SCIPallocBufferArray(scip, &incumbentedges, nedges) );
1063  SCIP_CALL( SCIPallocBufferArray(scip, &stnodes, nnodes) );
1064 
1065  /*
1066  * initialize incumbent solution
1067  */
1068 
1069  if( usestppool )
1070  {
1071  int i;
1072  STPSOL** poolsols = pool->sols;
1073 
1074  assert(pool->nedges == nedges);
1075 
1076  for( i = 0; i < nsols; i++ )
1077  if( *newsolindex == poolsols[i]->index )
1078  break;
1079  assert(i < nsols);
1080 
1081  BMScopyMemoryArray(incumbentedges, poolsols[i]->soledges, nedges);
1082  }
1083  else
1084  {
1085  SCIP_SOL* newsol;
1086  SCIP_SOL** sols = SCIPgetSols(scip);
1087  int i;
1088 
1089  assert(sols != NULL);
1090 
1091  for( i = 0; i < nsols; i++ )
1092  if( *newsolindex == SCIPsolGetIndex(sols[i]) )
1093  break;
1094  assert(i < nsols);
1095 
1096  newsol = sols[i];
1097 
1098  for( int e = 0; e < nedges; e++ )
1099  if( SCIPisEQ(scip, SCIPgetSolVal(scip, newsol, vars[e]), 1.0) )
1100  incumbentedges[e] = CONNECT;
1101  else
1102  incumbentedges[e] = UNKNOWN;
1103  }
1104 
1105  /* get objective value of incumbent */
1106  incumentobj = graph_sol_getObj(graph->cost, incumbentedges, 0.0, nedges);
1107 
1108  SCIPdebugMessage("REC: incumbent obj: %f \n", incumentobj);
1109 
1110  if( heur == NULL )
1111  {
1112  heur = SCIPfindHeur(scip, "rec");
1113  assert(heur != NULL);
1114  heurdata = SCIPheurGetData(heur);
1115  assert(heurdata != NULL);
1116  }
1117 
1118  /* main loop (for recombination) */
1119  for( int v = 0, failcount = 0; v < CYCLES_PER_RUN * runs && !SCIPisStopped(scip); v++ )
1120  {
1121  GRAPH* solgraph;
1122  IDX* curr;
1123  int* edgeweight;
1124  int* edgeancestor;
1125  SCIP_Real pobj;
1126  SCIP_Bool success;
1127  SCIP_Bool randomize;
1128  int randn;
1129 
1130  if( SCIPrandomGetInt(heurdata->randnumgen, 0, 1) == 1 )
1131  randomize = TRUE;
1132  else
1133  randomize = FALSE;
1134 
1135  if( restrictheur )
1136  randn = SCIPrandomGetInt(heurdata->randnumgen, 0, 3);
1137  else
1138  randn = SCIPrandomGetInt(heurdata->randnumgen, 0, 6);
1139 
1140  if( (randn <= 2) || (nsols < 3) )
1141  heurdata->nusedsols = 2;
1142  else if( (randn <= 4) || (nsols < 4) )
1143  heurdata->nusedsols = 3;
1144  else
1145  heurdata->nusedsols = 4;
1146 
1147  SCIPdebugMessage("REC: merge %d solutions \n", heurdata->nusedsols);
1148 
1149  /* build a new graph, consisting of several solutions */
1150  SCIP_CALL( buildsolgraph(scip, pool, heurdata, graph, &solgraph, incumbentedges, newsolindex,
1151  &edgeancestor, &edgeweight, &success, randomize, TRUE) );
1152 
1153  /* valid graph built? */
1154  if( success )
1155  {
1156  IDX** ancestors;
1157  GRAPH* psolgraph;
1158  STP_Bool* solnodes = NULL;
1159  int* soledges = NULL;
1160  SCIP_HEURDATA* tmheurdata;
1161  int nsoledges;
1162 
1163  SCIPdebugMessage("REC: solution successfully built \n");
1164 
1165  assert(SCIPfindHeur(scip, "TM") != NULL);
1166  tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
1167 
1168  assert(graph_valid(solgraph));
1169 
1170  /* reduce new graph */
1171  if( probtype == STP_RPCSPG || probtype == STP_DHCSTP || probtype == STP_DCSTP
1172  || probtype == STP_NWSPG || probtype == STP_SAP || probtype == STP_RMWCSP )
1173  SCIP_CALL( reduce(scip, &solgraph, &pobj, 0, 5, FALSE) );
1174  else
1175  SCIP_CALL( reduce(scip, &solgraph, &pobj, 2, 5, FALSE) );
1176 
1177  SCIP_CALL( graph_pack(scip, solgraph, &psolgraph, FALSE) );
1178 
1179  solgraph = psolgraph;
1180  ancestors = solgraph->ancestors;
1181  nsoledges = solgraph->edges;
1182 
1183  /* if graph reduction solved the whole problem, solgraph has only one node */
1184  if( solgraph->terms > 1 )
1185  {
1186  SCIP_Real* cost;
1187  SCIP_Real* costrev;
1188  SCIP_Real* orgprize = NULL;
1189  SCIP_Real* nodepriority;
1190  SCIP_Real maxcost;
1191  int best_start;
1192  const int nsolnodes = solgraph->knots;
1193 
1194  SCIPdebugMessage("REC: graph not completely reduced, nodes: %d, edges: %d, terminals: %d \n", solgraph->knots, nsoledges, solgraph->terms);
1195 
1196  /* allocate memory */
1197  SCIP_CALL( SCIPallocBufferArray(scip, &soledges, nsoledges) );
1198  SCIP_CALL( SCIPallocBufferArray(scip, &cost, nsoledges) );
1199  SCIP_CALL( SCIPallocBufferArray(scip, &costrev, nsoledges) );
1200  SCIP_CALL( SCIPallocBufferArray(scip, &nodepriority, nsolnodes) );
1201 
1202  if( pcmw )
1203  SCIP_CALL( SCIPallocBufferArray(scip, &orgprize, nsolnodes) );
1204 
1205  for( int i = 0; i < nsolnodes; i++ )
1206  nodepriority[i] = 0.0;
1207 
1208  /*
1209  * 1. modify edge costs
1210  */
1211 
1212  /* copy edge costs */
1213  BMScopyMemoryArray(cost, solgraph->cost, nsoledges);
1214 
1215  maxcost = 0.0;
1216  for( int e = 0; e < nsoledges; e++ )
1217  {
1218  SCIP_Real avg = 0.0;
1219  int i = 0;
1220  SCIP_Bool fixed = FALSE;
1221 
1222  curr = ancestors[e];
1223 
1224  if( curr != NULL )
1225  {
1226  while( curr != NULL )
1227  {
1228  i++;
1229  avg += edgeweight[curr->index];
1230  if( !usestppool && SCIPvarGetUbGlobal(vars[edgeancestor[curr->index]] ) < 0.5 )
1231  fixed = TRUE;
1232 
1233  curr = curr->parent;
1234  }
1235  avg = avg / (double) i;
1236  assert(avg >= 1);
1237  }
1238  /* is an ancestor edge fixed? */
1239  if( fixed )
1240  {
1241  cost[e] = BLOCKED;
1242 
1243  nodepriority[solgraph->head[e]] /= 2.0;
1244  nodepriority[solgraph->tail[e]] /= 2.0;
1245  }
1246 
1247  nodepriority[solgraph->head[e]] += avg - 1.0;
1248  nodepriority[solgraph->tail[e]] += avg - 1.0;
1249 
1250  cost[e] *= edgecostmultiplier(scip, heurdata, avg);
1251 
1252  if( probtype == STP_DHCSTP && SCIPisLT(scip, cost[e], BLOCKED) && SCIPisGT(scip, cost[e], maxcost) )
1253  maxcost = cost[e];
1254  }
1255 
1256  /* adapted prizes */
1257  if( pcmw )
1258  {
1259  const int solgraphroot = solgraph->source;
1260  SCIP_Real* const prize = solgraph->prize;
1261 
1262  assert(prize != NULL);
1263  assert(orgprize != NULL);
1264  assert(solgraph->extended);
1265 
1266 #ifndef NDEBUG
1267  graph_pc_2org(solgraph);
1268  assert(graph_pc_term2edgeConsistent(solgraph));
1269  graph_pc_2trans(solgraph);
1270 #endif
1271 
1272  for( int k = 0; k < nsolnodes; k++ )
1273  {
1274  if( Is_pterm(solgraph->term[k]) && k != solgraphroot )
1275  {
1276  int e;
1277  const int term = solgraph->head[solgraph->term2edge[k]];
1278  orgprize[k] = prize[k];
1279 
1280  assert(term >= 0);
1281  assert(Is_term(solgraph->term[term]));
1282 
1283  for( e = solgraph->inpbeg[term]; e != EAT_LAST; e = solgraph->ieat[e] )
1284  if( solgraph->tail[e] == solgraphroot )
1285  break;
1286 
1287  assert(e != EAT_LAST);
1288 
1289  prize[k] = cost[e];
1290  assert(solgraph->cost[e] > 0);
1291  }
1292  }
1293  }
1294 
1295  for( int e = 0; e < nsoledges; e++ )
1296  {
1297  costrev[e] = cost[flipedge(e)];
1298  soledges[e] = UNKNOWN;
1299  }
1300 
1301  /* initialize shortest path algorithm */
1302  SCIP_CALL( graph_path_init(scip, solgraph) );
1303 
1304  /*
1305  * 2. compute solution
1306  */
1307 
1308  // todo: run prune heuristic with changed weights!
1309 
1310  /* run TM heuristic */
1311  SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, solgraph, NULL, &best_start, soledges, heurdata->ntmruns,
1312  solgraph->source, cost, costrev, &hopfactor, nodepriority, maxcost, &success, FALSE) );
1313 
1314  assert(SCIPisStopped(scip) || success);
1315  assert(SCIPisStopped(scip) || graph_sol_valid(scip, solgraph, soledges));
1316 
1317  /* reset vertex weights */
1318  if( pcmw )
1319  {
1320  SCIP_Real* const prize = solgraph->prize;
1321 
1322  assert(orgprize != NULL);
1323 
1324  for( int k = 0; k < nsolnodes; k++ )
1325  if( Is_pterm(solgraph->term[k]) && k != solgraph->source )
1326  prize[k] = orgprize[k];
1327 
1328  SCIPfreeBufferArray(scip, &orgprize);
1329  }
1330 
1331  assert(graph_valid(solgraph));
1332 
1333  /* run local heuristic (with original costs) */
1334  if( !SCIPisStopped(scip) && probtype != STP_DHCSTP && probtype != STP_DCSTP
1335  && probtype != STP_SAP && probtype != STP_NWSPG && probtype != STP_RMWCSP )
1336  {
1337  SCIP_CALL( SCIPStpHeurLocalRun(scip, solgraph, solgraph->cost, soledges) );
1338  assert(graph_sol_valid(scip, solgraph, soledges));
1339  }
1340 
1341  graph_path_exit(scip, solgraph);
1342 
1343  SCIPfreeBufferArray(scip, &nodepriority);
1344  SCIPfreeBufferArray(scip, &costrev);
1345  SCIPfreeBufferArray(scip, &cost);
1346  } /* graph->terms > 1 */
1347 
1348  for( int i = 0; i < nedges; i++ )
1349  newsoledges[i] = UNKNOWN;
1350 
1351  for( int i = 0; i < nnodes; i++ )
1352  stnodes[i] = FALSE;
1353 
1354  if( pcmw )
1355  {
1356  SCIP_CALL( SCIPallocBufferArray(scip, &solnodes, solgraph->orgknots) );
1357 
1358  for( int i = 0; i < solgraph->orgknots; i++ )
1359  solnodes[i] = FALSE;
1360  }
1361 
1362  /*
1363  * retransform solution found by heuristic
1364  */
1365 
1366  if( solgraph->terms > 1 )
1367  {
1368  assert(soledges != NULL);
1369  for( int e = 0; e < nsoledges; e++ )
1370  {
1371  if( soledges[e] == CONNECT )
1372  {
1373  /* iterate through list of ancestors */
1374  if( probtype != STP_DCSTP )
1375  {
1376  curr = ancestors[e];
1377 
1378  while( curr != NULL )
1379  {
1380  int i = edgeancestor[curr->index];
1381 
1382  stnodes[graph->head[i]] = TRUE;
1383  stnodes[graph->tail[i]] = TRUE;
1384 
1385  if( pcmw )
1386  {
1387  assert(solgraph->orghead[curr->index] < solgraph->orgknots);
1388  assert(solgraph->orgtail[curr->index] < solgraph->orgknots);
1389 
1390  solnodes[solgraph->orghead[curr->index]] = TRUE;
1391  solnodes[solgraph->orgtail[curr->index]] = TRUE;
1392  }
1393 
1394  curr = curr->parent;
1395  }
1396  }
1397  else
1398  {
1399  curr = ancestors[e];
1400  while( curr != NULL )
1401  {
1402  int i = edgeancestor[curr->index];
1403  newsoledges[i] = CONNECT;
1404 
1405  curr = curr->parent;
1406  }
1407  }
1408  }
1409  }
1410  } /* if( solgraph->terms > 1 ) */
1411 
1412  /* retransform edges fixed during graph reduction */
1413  if( probtype != STP_DCSTP )
1414  {
1415  curr = solgraph->fixedges;
1416 
1417  while( curr != NULL )
1418  {
1419  const int i = edgeancestor[curr->index];
1420 
1421  stnodes[graph->head[i]] = TRUE;
1422  stnodes[graph->tail[i]] = TRUE;
1423  if( pcmw )
1424  {
1425  assert(solgraph->orghead[curr->index] < solgraph->orgknots);
1426  assert(solgraph->orgtail[curr->index] < solgraph->orgknots);
1427 
1428  solnodes[solgraph->orghead[curr->index]] = TRUE;
1429  solnodes[solgraph->orgtail[curr->index]] = TRUE;
1430  }
1431 
1432  curr = curr->parent;
1433  }
1434  }
1435  else
1436  {
1437  curr = solgraph->fixedges;
1438  while( curr != NULL )
1439  {
1440  const int i = edgeancestor[curr->index];
1441  newsoledges[i] = CONNECT;
1442  }
1443  }
1444 
1445  if( pcmw )
1446  {
1447  STP_Bool* pcancestoredges;
1448  SCIP_CALL( SCIPallocBufferArray(scip, &pcancestoredges, solgraph->orgedges) );
1449 
1450  for( int k = 0; k < solgraph->orgedges; k++ )
1451  pcancestoredges[k] = FALSE;
1452 
1453  SCIP_CALL( graph_sol_markPcancestors(scip, solgraph->pcancestors, solgraph->orgtail, solgraph->orghead, solgraph->orgknots,
1454  solnodes, pcancestoredges, NULL, NULL, NULL ) );
1455 
1456  for( int k = 0; k < solgraph->orgedges; k++ )
1457  if( pcancestoredges[k] )
1458  {
1459  const int i = edgeancestor[k];
1460  stnodes[graph->tail[i]] = TRUE;
1461  stnodes[graph->head[i]] = TRUE;
1462  }
1463 
1464  SCIPfreeBufferArray(scip, &pcancestoredges);
1465  SCIPfreeBufferArray(scip, &solnodes);
1466  }
1467 
1468  SCIPfreeMemoryArray(scip, &edgeancestor);
1469 
1470  graph_free(scip, &solgraph, TRUE);
1471 
1472  /* prune solution (in the original graph) */
1473  if( pcmw )
1474  SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, newsoledges, stnodes) );
1475  else if( probtype == STP_DCSTP )
1476  SCIP_CALL( SCIPStpHeurTMBuildTreeDc(scip, graph, newsoledges, stnodes) );
1477  else
1478  SCIP_CALL( SCIPStpHeurTMPrune(scip, graph, graph->cost, 0, newsoledges, stnodes) );
1479 
1480  assert(graph_sol_valid(scip, graph, newsoledges));
1481  pobj = graph_sol_getObj(graph->cost, newsoledges, 0.0, nedges);
1482 
1483  SCIPdebugMessage("REC: new obj: %f \n", pobj);
1484 
1485  /* improved incumbent solution? */
1486  if( !SCIPisStopped(scip) && SCIPisLT(scip, pobj, incumentobj) )
1487  {
1488  SCIPdebugMessage("...is better! \n");
1489  incumentobj = pobj;
1490  *solfound = TRUE;
1491  BMScopyMemoryArray(incumbentedges, newsoledges, nedges);
1492  }
1493  else
1494  {
1495  failcount++;
1496  }
1497 
1498  SCIPfreeBufferArrayNull(scip, &soledges);
1499  }
1500  else
1501  {
1502  failcount++;
1503  }
1504 
1505  SCIPfreeMemoryArray(scip, &edgeweight);
1506 
1507  /* too many fails? */
1508  if( failcount >= REC_MAX_FAILS )
1509  {
1510  SCIPdebugMessage("REC: too many fails, exit! \n");
1511  break;
1512  }
1513  }
1514 
1515  SCIPfreeBufferArray(scip, &stnodes);
1516 
1517  /* improving solution found? */
1518  if( *solfound )
1519  {
1520  /*
1521  * store best solution in pool
1522  */
1523 
1524  SCIPdebugMessage("REC: better solution found ... ");
1525 
1526  if( usestppool )
1527  {
1528  const int maxindex = pool->maxindex;
1529  SCIP_CALL( SCIPStpHeurRecAddToPool(scip, incumentobj, incumbentedges, pool, solfound) );
1530 
1531  if( *solfound )
1532  {
1533  assert(pool->maxindex == maxindex + 1);
1534  *newsolindex = maxindex + 1;
1535  SCIPdebugMessage("and added! \n");
1536  }
1537  else
1538  {
1539  *newsolindex = -1;
1540  }
1541  }
1542  else
1543  {
1544  SCIP_SOL* sol = NULL;
1545  SCIP_Real* nval = NULL;
1546 
1547  SCIP_CALL( SCIPallocBufferArray(scip, &nval, nedges) );
1548 
1549  for( int e = 0; e < nedges; e++ )
1550  {
1551  if( newsoledges[e] == CONNECT )
1552  nval[e] = 1.0;
1553  else
1554  nval[e] = 0.0;
1555  }
1556 
1557  SCIP_CALL( SCIPprobdataAddNewSol(scip, nval, sol, heur, solfound) );
1558 
1559  if( *solfound )
1560  {
1561  SCIP_SOL** sols = SCIPgetSols(scip);
1562  int idx = -1;
1563 
1564  SCIPdebugMessage("and added! \n");
1565 
1566  nsols = SCIPgetNSols(scip);
1567 
1568  assert(nsols > 0);
1569 
1570  for( int i = 1; i < nsols; i++ )
1571  if( SCIPsolGetIndex(sols[i]) > idx )
1572  idx = SCIPsolGetIndex(sols[i]);
1573 
1574  assert(idx >= 0);
1575 
1576  *newsolindex = idx;
1577  }
1578  SCIPfreeBufferArray(scip, &nval);
1579 
1580  }
1581  }
1582  else
1583  {
1584  *newsolindex = -1;
1585  }
1586 
1587  SCIPfreeBufferArray(scip, &incumbentedges);
1588  SCIPfreeBufferArray(scip, &newsoledges);
1589 
1590 
1591  return SCIP_OKAY;
1592 }
1593 
1594 /** heuristic to exclude vertices or edges from a given solution (and inserting other edges) to improve objective */
1596  SCIP* scip, /**< SCIP data structure */
1597  const GRAPH* graph, /**< graph structure */
1598  const int* result, /**< edge solution array (UNKNOWN/CONNECT) */
1599  int* newresult, /**< new edge solution array (UNKNOWN/CONNECT) */
1600  int* dnodemap, /**< node array for internal use */
1601  STP_Bool* stvertex, /**< node array for internally marking solution vertices */
1602  SCIP_Bool* success /**< solution improved? */
1603  )
1604 {
1605  SCIP_HEURDATA* tmheurdata;
1606  GRAPH* newgraph;
1607  int* unodemap;
1608  STP_Bool* solnodes;
1609  SCIP_Real dummy;
1610  int j;
1611  int nsolterms;
1612  int nsoledges;
1613  int nsolnodes;
1614  int best_start;
1615  const int root = graph->source;
1616  const int nedges = graph->edges;
1617  const int nnodes = graph->knots;
1618 
1619  assert(scip != NULL);
1620  assert(graph != NULL);
1621  assert(result != NULL);
1622  assert(stvertex != NULL);
1623  assert(success != NULL);
1624  assert(dnodemap != NULL);
1625  assert(newresult != NULL);
1626 
1627  *success = FALSE;
1628  assert(graph->stp_type == STP_MWCSP);
1629  assert(graph_valid(graph));
1630 
1631  /* killed solution edge? */
1632  for( int e = 0; e < nedges; e++ )
1633  if( result[e] == CONNECT && graph->oeat[e] == EAT_FREE )
1634  return SCIP_OKAY;
1635 
1636  /*** 1. step: for solution S and original graph (V,E) initialize new graph (V[S], (V[S] X V[S]) \cup E) ***/
1637 
1638  BMSclearMemoryArray(stvertex, nnodes);
1639 
1640  nsolnodes = 1;
1641  nsolterms = 0;
1642  stvertex[root] = 1;
1643 
1644  /* mark nodes in solution */
1645  for( int e = 0; e < nedges; e++ )
1646  {
1647  if( result[e] == CONNECT )
1648  {
1649  int tail = graph->tail[e];
1650  int head = graph->head[e];
1651 
1652  if( tail == root )
1653  {
1654  /* there might be only one node */
1655  if( Is_pterm(graph->term[head]) )
1656  {
1657  stvertex[head] = 1;
1658  nsolterms++;
1659  nsolnodes++;
1660  }
1661 
1662  continue;
1663  }
1664 
1665  assert(!stvertex[head] );
1666  stvertex[head] = 1;
1667 
1668  if( Is_pterm(graph->term[head]) )
1669  nsolterms++;
1670  nsolnodes++;
1671  }
1672  }
1673 
1674  assert(nsolterms > 0);
1675 
1676  nsoledges = 0;
1677 
1678  /* count edges of new graph */
1679  for( int i = 0; i < nedges; i += 2 )
1680  if( stvertex[graph->tail[i]] && stvertex[graph->head[i]] && graph->oeat[i] != EAT_FREE )
1681  nsoledges++;
1682 
1683  nsoledges *= 2;
1684 
1685  /* create new graph */
1686  SCIP_CALL( graph_init(scip, &newgraph, nsolnodes, nsoledges, 1) );
1687 
1688  newgraph->stp_type = graph->stp_type;
1689  newgraph->extended = TRUE;
1690 
1691  SCIP_CALL( SCIPallocBufferArray(scip, &unodemap, nsolnodes) );
1692  SCIP_CALL( graph_pc_init(scip, newgraph, nsolnodes, nsolnodes) );
1693 
1694  j = 0;
1695  for( int i = 0; i < nnodes; i++ )
1696  {
1697  if( stvertex[i] )
1698  {
1699  if( (!Is_term(graph->term[i])) )
1700  newgraph->prize[j] = graph->prize[i];
1701  else
1702  newgraph->prize[j] = 0.0;
1703 
1704  graph_knot_add(newgraph, graph->term[i]);
1705  unodemap[j] = i;
1706  dnodemap[i] = j++;
1707  }
1708  else
1709  {
1710  dnodemap[i] = -1;
1711  }
1712  }
1713 
1714  assert(j == nsolnodes);
1715 
1716  /* set root */
1717  newgraph->source = dnodemap[root];
1718  assert(newgraph->source >= 0 && newgraph->source < nsolnodes);
1719 
1720  /* add edges */
1721  for( int i = 0; i < nedges; i += 2 )
1722  if( stvertex[graph->tail[i]] && stvertex[graph->head[i]] && graph->oeat[i] != EAT_FREE )
1723  {
1724  const int tail = graph->tail[i];
1725  const int head = graph->head[i];
1726  const int dtail = dnodemap[tail];
1727  const int dhead = dnodemap[head];
1728 
1729  graph_pc_updateTerm2edge(newgraph, graph, dtail, dhead, tail, head);
1730  graph_edge_add(scip, newgraph, dtail, dhead, graph->cost[i], graph->cost[i + 1]);
1731  }
1732 
1733  assert(newgraph->edges == nsoledges);
1734 
1735 
1736  /*** step 2: presolve ***/
1737 
1738  newgraph->norgmodelknots = nsolnodes;
1739 
1740  dummy = 0.0;
1741  SCIP_CALL( reduce(scip, &newgraph, &dummy, 1, 5, FALSE) );
1742 
1743 
1744  /*** step 3: compute solution on new graph ***/
1745 
1746 
1747  /* get TM heuristic data */
1748  assert(SCIPfindHeur(scip, "TM") != NULL);
1749  tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
1750 
1751  SCIP_CALL( graph_path_init(scip, newgraph) );
1752 
1753  /* compute Steiner tree to obtain upper bound */
1754  best_start = newgraph->source;
1755  SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, newgraph, NULL, &best_start, newresult, MIN(50, nsolterms), newgraph->source, newgraph->cost,
1756  newgraph->cost, &dummy, NULL, 0.0, success, FALSE) );
1757 
1758  graph_path_exit(scip, newgraph);
1759 
1760  assert(*success);
1761  assert(graph_sol_valid(scip, newgraph, newresult));
1762 
1763 
1764  /*** step 4: retransform solution to original graph ***/
1765 
1766 
1767  BMSclearMemoryArray(stvertex, nnodes);
1768 
1769  for( int e = 0; e < nsoledges; e++ )
1770  if( newresult[e] == CONNECT )
1771  marksolverts(newgraph, newgraph->ancestors[e], unodemap, stvertex);
1772 
1773  /* retransform edges fixed during graph reduction */
1774  marksolverts(newgraph, newgraph->fixedges, unodemap, stvertex);
1775 
1776  for( int k = 0; k < nsolnodes; k++ )
1777  if( stvertex[unodemap[k]] )
1778  marksolverts(newgraph, newgraph->pcancestors[k], unodemap, stvertex);
1779 
1780  SCIP_CALL(SCIPallocBufferArray(scip, &solnodes, nsolnodes));
1781 
1782  for( int k = 0; k < nsolnodes; k++ )
1783  solnodes[k] = FALSE;
1784 
1785  for( int k = 0; k < nnodes; k++ )
1786  if( stvertex[k] )
1787  {
1788  assert(dnodemap[k] < nsolnodes);
1789  solnodes[dnodemap[k]] = TRUE;
1790  }
1791 
1792  SCIP_CALL( graph_sol_markPcancestors(scip, newgraph->pcancestors, newgraph->orgtail, newgraph->orghead, nsolnodes, solnodes, NULL, NULL, NULL, NULL ) );
1793 
1794  for( int k = 0; k < nsolnodes; k++ )
1795  if( solnodes[k] )
1796  stvertex[unodemap[k]] = TRUE;
1797 
1798  SCIPfreeBufferArray(scip, &solnodes);
1799 
1800  for( int e = 0; e < nedges; e++ )
1801  newresult[e] = UNKNOWN;
1802 
1803  SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, newresult, stvertex) );
1804 
1805  /* solution better than original one? */
1806 
1807  if( SCIPisLT(scip, graph_sol_getObj(graph->cost, newresult, 0.0, nedges),
1808  graph_sol_getObj(graph->cost, result, 0.0, nedges)) )
1809  {
1810  *success = TRUE;
1811  SCIPdebugMessage("success %f < %f \n", graph_sol_getObj(graph->cost, newresult, 0.0, nedges), graph_sol_getObj(graph->cost, result, 0.0, nedges));
1812  }
1813  else
1814  {
1815  SCIPdebugMessage("no improvements %f >= %f \n", graph_sol_getObj(graph->cost, newresult, 0.0, nedges), graph_sol_getObj(graph->cost, result, 0.0, nedges));
1816  *success = FALSE;
1817  }
1818 
1819  assert(graph_sol_valid(scip, graph, newresult));
1820 
1821  if( !graph_sol_valid(scip, graph, newresult) )
1822  *success = FALSE;
1823 
1824  SCIPfreeBufferArray(scip, &unodemap);
1825  graph_free(scip, &newgraph, TRUE);
1826 
1827  return SCIP_OKAY;
1828 }
1829 
1830 
1831 
1832 
1833 /*
1834  * Callback methods of primal heuristic
1835  */
1836 
1837 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
1838 static
1840 {
1841  SCIP_HEURDATA* heurdata;
1842 
1843  heurdata = SCIPheurGetData(heur);
1844  assert(heurdata != NULL);
1845 
1846  SCIPheurSetData(heur, heurdata);
1847 
1848  return SCIP_OKAY;
1849 }
1850 
1851 /** copy method for primal heuristic plugins (called when SCIP copies plugins) */
1852 static
1854 { /*lint --e{715}*/
1855  assert(scip != NULL);
1856  assert(heur != NULL);
1857  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1858 
1859  /* call inclusion method of primal heuristic */
1861 
1862  return SCIP_OKAY;
1863 }
1864 
1865 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
1866 static
1868 { /*lint --e{715}*/
1869  SCIP_HEURDATA* heurdata;
1870 
1871  assert(heur != NULL);
1872  assert(scip != NULL);
1873 
1874  /* get heuristic data */
1875  heurdata = SCIPheurGetData(heur);
1876  assert(heurdata != NULL);
1877 
1878  /* free random number generator */
1879  SCIPfreeRandom(scip, &heurdata->randnumgen);
1880 
1881  /* free heuristic data */
1882  SCIPfreeMemory(scip, &heurdata);
1883  SCIPheurSetData(heur, NULL);
1884 
1885  return SCIP_OKAY;
1886 }
1887 
1888 /** initialization method of primal heuristic (called after problem was transformed) */
1889 static
1891 { /*lint --e{715}*/
1892 
1893  assert(heur != NULL);
1894  assert(scip != NULL);
1895 
1896  return SCIP_OKAY;
1897 }
1898 
1899 /** solving process initialization method of primal heuristic (called when branch and bound process is about to begin) */
1900 static
1901 SCIP_DECL_HEURINITSOL(heurInitsolRec)
1902 { /*lint --e{715}*/
1903  SCIP_HEURDATA* heurdata;
1904 
1905  assert(heur != NULL);
1906  assert(scip != NULL);
1907 
1908  heurdata = SCIPheurGetData(heur);
1909  assert(heurdata != NULL);
1910 
1911  /* initialize data */
1912  heurdata->nselectedsols = 0;
1913  heurdata->nusedsols = DEFAULT_NUSEDSOLS;
1914  heurdata->randseed = DEFAULT_RANDSEED;
1915  heurdata->nselectedsols = 0;
1916  heurdata->ncalls = 0;
1917  heurdata->nlastsols = 0;
1918  heurdata->lastsolindex = -1;
1919  heurdata->bestsolindex = -1;
1920  heurdata->nfailures = 0;
1921 
1922 
1923  return SCIP_OKAY;
1924 }
1925 
1926 /** solving process deinitialization method of primal heuristic (called before branch and bound process data is freed) */
1927 static
1928 SCIP_DECL_HEUREXITSOL(heurExitsolRec)
1929 { /*lint --e{715}*/
1930 
1931  return SCIP_OKAY;
1932 }
1933 
1934 /** execution method of primal heuristic */
1935 static
1937 { /*lint --e{715}*/
1938  SCIP_HEURDATA* heurdata;
1939  SCIP_PROBDATA* probdata;
1940  SCIP_VAR** vars;
1941  SCIP_SOL** sols;
1942  GRAPH* graph;
1943  SCIP_Longint nallsols;
1944  SCIP_Bool pcmw;
1945  SCIP_Bool solfound;
1946  SCIP_Bool restrictheur;
1947 
1948  int runs;
1949  int nsols;
1950  int solindex;
1951  int probtype;
1952  int newsolindex;
1953  int nreadysols;
1954  int bestsolindex;
1955 
1956  assert(heur != NULL);
1957  assert(scip != NULL);
1958  assert(strcmp(SCIPheurGetName(heur), HEUR_NAME) == 0);
1959  assert(result != NULL);
1960 
1961  /* get heuristic data */
1962  heurdata = SCIPheurGetData(heur);
1963  assert(heurdata != NULL);
1964 
1965  /* get problem data */
1966  probdata = SCIPgetProbData(scip);
1967  assert(probdata != NULL);
1968 
1969  /* get graph */
1970  graph = SCIPprobdataGetGraph(probdata);
1971  assert(graph != NULL);
1972 
1973  probtype = graph->stp_type;
1974  *result = SCIP_DIDNOTRUN;
1975 
1976  if( probtype == STP_RMWCSP )
1977  return SCIP_OKAY;
1978 
1979  SCIPdebugMessage("REC: checking ... \n");
1980 
1981  pcmw = (probtype == STP_PCSPG || probtype == STP_MWCSP || probtype == STP_RPCSPG || probtype == STP_RMWCSP);
1982  nallsols = SCIPgetNSolsFound(scip);
1983  nreadysols = SCIPgetNSols(scip);
1984 
1985  /* only call heuristic if sufficiently many solutions are available */
1986  if( nreadysols < DEFAULT_NUSEDSOLS )
1987  return SCIP_OKAY;
1988 
1989  if( probtype == STP_DCSTP )
1990  return SCIP_OKAY;
1991 
1992  /* suspend heuristic? */
1993  if( pcmw || probtype == STP_DHCSTP || probtype == STP_DCSTP )
1994  {
1995  int i;
1996  if( heurdata->ncalls == 0 )
1997  i = 0;
1998  else if( heurdata->maxfreq )
1999  i = 1;
2000  else if( probtype == STP_RPCSPG || probtype == STP_DCSTP )
2001  i = MIN(2 * heurdata->nwaitingsols, 2 * heurdata->nfailures);
2002  else
2003  i = MIN(heurdata->nwaitingsols, heurdata->nfailures);
2004 
2005  if( nallsols <= heurdata->nlastsols + i )
2006  return SCIP_OKAY;
2007  }
2008  else
2009  {
2010  int i;
2011  if( heurdata->maxfreq )
2012  i = 1;
2013  else
2014  i = MIN(heurdata->nwaitingsols, heurdata->nfailures);
2015 
2016  if( nallsols <= heurdata->nlastsols + i && heurdata->bestsolindex == SCIPsolGetIndex(SCIPgetBestSol(scip)) )
2017  return SCIP_OKAY;
2018  }
2019 
2020  /* get edge variables */
2021  vars = SCIPprobdataGetVars(scip);
2022  assert(vars != NULL);
2023  assert(vars[0] != NULL);
2024 
2025  heurdata->ncalls++;
2026 
2027  restrictheur = (graph->terms > BOUND_MAXNTERMINALS && graph->edges > BOUND_MAXNEDGES);
2028 
2029  if( restrictheur )
2030  runs = RUNS_RESTRICTED;
2031  else
2032  runs = RUNS_NORMAL;
2033 
2034  if( runs > nreadysols )
2035  runs = nreadysols;
2036 
2037  assert(runs > 0);
2038 
2039  sols = SCIPgetSols(scip);
2040  assert(sols != NULL);
2041 
2042  bestsolindex = SCIPsolGetIndex(SCIPgetBestSol(scip));
2043 
2044  /* first run or exotic variant? */
2045  if( heurdata->lastsolindex == -1 || probtype == STP_DHCSTP || probtype == STP_DCSTP )
2046  newsolindex = bestsolindex;
2047  else
2048  {
2049  /* get best new solution */
2050 
2051  SCIP_Real minobj = FARAWAY;
2052  const int lastindex = heurdata->lastsolindex;
2053 
2054  newsolindex = -1;
2055  assert(lastindex >= 0);
2056 
2057  for( int i = 0; i < nreadysols; i++ )
2058  {
2059  const int currindex = SCIPsolGetIndex(sols[i]);
2060  if( currindex > lastindex && SCIPisLT(scip, SCIPgetSolOrigObj(scip, sols[i]), minobj) )
2061  {
2062  newsolindex = currindex;
2063  minobj = SCIPgetSolOrigObj(scip, sols[i]);
2064  SCIPdebugMessage("REC: better start solution, obj: %f \n", minobj);
2065  }
2066  }
2067  if( newsolindex == -1 )
2068  {
2069  SCIPdebugMessage("REC: random start solution\n");
2070  newsolindex = SCIPsolGetIndex(sols[SCIPrandomGetInt(heurdata->randnumgen, 0, nreadysols)]);
2071  }
2072  }
2073 
2074  /* run the actual heuristic */
2075  SCIP_CALL( SCIPStpHeurRecRun(scip, NULL, heur, heurdata, graph, vars, &newsolindex, runs, nreadysols, restrictheur, &solfound) );
2076 
2077  /* save latest solution index */
2078  solindex = 0;
2079  nsols = SCIPgetNSols(scip);
2080  assert(nsols > 0);
2081  sols = SCIPgetSols(scip);
2082 
2083  for( int i = 1; i < nsols; i++ )
2084  if( SCIPsolGetIndex(sols[i]) > SCIPsolGetIndex(sols[solindex]) )
2085  solindex = i;
2086 
2087  if( solfound )
2088  *result = SCIP_FOUNDSOL;
2089 
2090  if( SCIPsolGetIndex(SCIPgetBestSol(scip)) == bestsolindex )
2091  heurdata->nfailures++;
2092  else
2093  heurdata->nfailures = 0;
2094 
2095  heurdata->lastsolindex = SCIPsolGetIndex(sols[solindex]);
2096  heurdata->bestsolindex = SCIPsolGetIndex(SCIPgetBestSol(scip));
2097  heurdata->nlastsols = SCIPgetNSolsFound(scip);
2098 
2099  return SCIP_OKAY;
2100 }
2101 
2102 
2103 /** creates the rec primal heuristic and includes it in SCIP */
2105  SCIP* scip /**< SCIP data structure */
2106  )
2107 {
2108  SCIP_HEURDATA* heurdata;
2109  SCIP_HEUR* heur;
2110 
2111  /* create Rec primal heuristic data */
2112  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
2113 
2114  /* include primal heuristic */
2115  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
2117  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecRec, heurdata) );
2118 
2119  assert(heur != NULL);
2120 
2121  /* set non-NULL pointers to callback methods */
2122  SCIP_CALL( SCIPsetHeurCopy(scip, heur, heurCopyRec) );
2123  SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeRec) );
2124  SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitRec) );
2125  SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitRec) );
2126  SCIP_CALL( SCIPsetHeurInitsol(scip, heur, heurInitsolRec) );
2127  SCIP_CALL( SCIPsetHeurExitsol(scip, heur, heurExitsolRec) );
2128 
2129  /* add rec primal heuristic parameters */
2130  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/nwaitingsols",
2131  "number of solution findings to be in abeyance",
2132  &heurdata->nwaitingsols, FALSE, DEFAULT_NWAITINGSOLS, 1, INT_MAX, NULL, NULL) );
2133 
2134  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/randseed",
2135  "random seed for heuristic",
2136  NULL, FALSE, DEFAULT_RANDSEED, 1, INT_MAX, paramChgdRandomseed, (SCIP_PARAMDATA*)heurdata) );
2137 
2138  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/maxnsols",
2139  "max size of solution pool for heuristic",
2140  &heurdata->maxnsols, FALSE, DEFAULT_MAXNSOLS, 5, INT_MAX, NULL, NULL) );
2141 
2142  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/ntmruns",
2143  "number of runs in TM",
2144  &heurdata->ntmruns, FALSE, DEFAULT_NTMRUNS, 1, INT_MAX, NULL, NULL) );
2145 
2146  SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/"HEUR_NAME"/maxfreq",
2147  "should the heuristic be executed at maximum frequeny?",
2148  &heurdata->maxfreq, FALSE, DEFAULT_MAXFREQREC, NULL, NULL) );
2149 
2150  heurdata->nusedsols = DEFAULT_NUSEDSOLS;
2151  heurdata->randseed = DEFAULT_RANDSEED;
2152 
2153 #ifdef WITH_UG
2154  heurdata->randseed += getUgRank();
2155 #endif
2156 
2157  /* create random number generator */
2158  SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen, heurdata->randseed, TRUE) );
2159 
2160  return SCIP_OKAY;
2161 }
SCIP_RETCODE graph_path_init(SCIP *, GRAPH *)
Definition: grphpath.c:444
SCIP_RETCODE SCIPsetHeurExitsol(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEUREXITSOL((*heurexitsol)))
Definition: scip_heur.c:312
#define BOUND_MAXNEDGES
Definition: heur_rec.c:65
void SCIPfreeRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen)
#define REC_MAX_FAILS
Definition: heur_rec.c:69
#define NULL
Definition: def.h:239
SCIP_RETCODE graph_pc_init(SCIP *, GRAPH *, int, int)
Definition: grphbase.c:766
#define HEUR_FREQ
Definition: heur_rec.c:51
int *RESTRICT head
Definition: grph.h:96
int *RESTRICT orgtail
Definition: grph.h:97
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
#define COST_MIN_POLY_x0
Definition: heur_rec.c:74
static SCIP_Bool isInPool(const int *soledges, const STPSOLPOOL *pool)
Definition: heur_rec.c:835
#define STP_GSTP
Definition: grph.h:49
static int edgecostmultiplier(SCIP *scip, SCIP_HEURDATA *heurdata, SCIP_Real avg)
Definition: heur_rec.c:129
Constraint handler for Steiner problems.
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:88
SCIP_PARAMDATA * SCIPparamGetData(SCIP_PARAM *param)
Definition: paramset.c:661
SCIP_Longint SCIPgetNSolsFound(SCIP *scip)
SCIP_Bool graph_valid(const GRAPH *)
Definition: grphbase.c:4324
int terms
Definition: grph.h:64
SCIP_RETCODE SCIPsetHeurExit(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEUREXIT((*heurexit)))
Definition: scip_heur.c:280
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
#define DEFAULT_NTMRUNS
Definition: heur_rec.c:61
struct SCIP_ParamData SCIP_PARAMDATA
Definition: type_paramset.h:76
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:72
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define REC_MIN_NSOLS
Definition: heur_rec.c:70
int *RESTRICT maxdeg
Definition: grph.h:78
#define EAT_LAST
Definition: grph.h:31
dual-ascent and reduction based primal heuristic for Steiner problems
#define BLOCKED
Definition: grph.h:157
SCIP_SOL ** SCIPgetSols(SCIP *scip)
Definition: scip_sol.c:2329
#define FALSE
Definition: def.h:65
#define HEUR_PRIORITY
Definition: heur_rec.c:50
int *RESTRICT inpbeg
Definition: grph.h:74
static SCIP_DECL_PARAMCHGD(paramChgdRandomseed)
Definition: heur_rec.c:111
static SCIP_DECL_HEURCOPY(heurCopyRec)
Definition: heur_rec.c:1853
#define STP_RMWCSP
Definition: grph.h:50
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: grphbase.c:3674
Problem data for stp problem.
#define TRUE
Definition: def.h:64
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
void graph_path_exit(SCIP *, GRAPH *)
Definition: grphpath.c:466
#define STP_DHCSTP
Definition: grph.h:48
SCIP_RETCODE SCIPStpIncludeHeurRec(SCIP *scip)
Definition: heur_rec.c:2104
struct SCIP_HeurData SCIP_HEURDATA
Definition: type_heur.h:51
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:114
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
#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:9372
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:142
void SCIPheurSetData(SCIP_HEUR *heur, SCIP_HEURDATA *heurdata)
Definition: heur.c:1175
int *RESTRICT orghead
Definition: grph.h:98
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:97
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
static SCIP_DECL_HEUREXEC(heurExecRec)
Definition: heur_rec.c:1936
void graph_pc_updateTerm2edge(GRAPH *, const GRAPH *, int, int, int, int)
Definition: grphbase.c:928
GRAPH * SCIPprobdataGetGraph(SCIP_PROBDATA *probdata)
#define DEFAULT_MAXNSOLS
Definition: heur_rec.c:58
SCIP_Bool graph_pc_term2edgeConsistent(const GRAPH *)
Definition: grphbase.c:853
IDX * fixedges
Definition: grph.h:85
SCIP_RETCODE graph_pack(SCIP *, GRAPH *, GRAPH **, SCIP_Bool)
Definition: grphbase.c:3984
#define RUNS_NORMAL
Definition: heur_rec.c:67
STPSOL ** sols
Definition: heur_rec.h:48
SCIP_Real SCIPvarGetUbGlobal(SCIP_VAR *var)
Definition: var.c:17353
int *RESTRICT oeat
Definition: grph.h:103
void graph_pc_2org(GRAPH *)
Definition: grphbase.c:964
#define CONNECT
Definition: grph.h:154
static SCIP_RETCODE buildsolgraph(SCIP *scip, const STPSOLPOOL *pool, SCIP_HEURDATA *heurdata, const GRAPH *graph, GRAPH **solgraph, const int *incumbentedges, int *incumbentindex, int **edgeancestor, int **edgeweight, SCIP_Bool *success, SCIP_Bool randomize, SCIP_Bool addedges)
Definition: heur_rec.c:466
SCIP_RETCODE SCIPsetHeurInitsol(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURINITSOL((*heurinitsol)))
Definition: scip_heur.c:296
const char * SCIPheurGetName(SCIP_HEUR *heur)
Definition: heur.c:1254
SCIP_Bool extended
Definition: grph.h:128
#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
IDX ** ancestors
Definition: grph.h:86
int orgedges
Definition: grph.h:93
#define CYCLES_PER_RUN
Definition: heur_rec.c:68
void graph_pc_2trans(GRAPH *)
Definition: grphbase.c:1002
static void marksolverts(GRAPH *g, IDX *curr, int *unodemap, STP_Bool *stvertex)
Definition: heur_rec.c:816
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
#define HEUR_FREQOFS
Definition: heur_rec.c:52
SCIP_RETCODE SCIPprobdataAddNewSol(SCIP *scip, SCIP_Real *nval, SCIP_SOL *sol, SCIP_HEUR *heur, SCIP_Bool *success)
#define STP_DCSTP
Definition: grph.h:43
SCIPInterval sqrt(const SCIPInterval &x)
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:143
SCIP_Real * prize
Definition: grph.h:82
#define RUNS_RESTRICTED
Definition: heur_rec.c:66
SCIP_Bool graph_sol_valid(SCIP *, const GRAPH *, const int *)
Definition: grphbase.c:3066
SCIP_RETCODE SCIPStpHeurRecExclude(SCIP *scip, const GRAPH *graph, const int *result, int *newresult, int *dnodemap, STP_Bool *stvertex, SCIP_Bool *success)
Definition: heur_rec.c:1595
int knots
Definition: grph.h:62
#define SCIP_CALL(x)
Definition: def.h:351
int * term2edge
Definition: grph.h:80
IDX ** pcancestors
Definition: grph.h:87
#define STP_NWSPG
Definition: grph.h:42
#define DEFAULT_NUSEDSOLS
Definition: heur_rec.c:59
int orgknots
Definition: grph.h:63
Improvement heuristic for Steiner problems.
#define Is_gterm(a)
Definition: grph.h:170
reduction-based primal heuristic for Steiner problems
#define FARAWAY
Definition: grph.h:156
SCIP_RETCODE SCIPcreateRandom(SCIP *scip, SCIP_RANDNUMGEN **randnumgen, unsigned int initialseed, SCIP_Bool useglobalseed)
#define STP_SPG
Definition: grph.h:38
static SCIP_DECL_HEURFREE(heurFreeRec)
Definition: heur_rec.c:1867
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:130
static SCIP_DECL_HEURINIT(heurInitRec)
Definition: heur_rec.c:1890
public data structures and miscellaneous methods
#define DEFAULT_NWAITINGSOLS
Definition: heur_rec.c:62
#define SCIP_Bool
Definition: def.h:62
int *RESTRICT ieat
Definition: grph.h:102
SCIP_Real obj
Definition: heur_rec.h:40
#define STP_MWCSP
Definition: grph.h:47
int *RESTRICT tail
Definition: grph.h:95
static SCIP_RETCODE selectsols(SCIP *scip, const STPSOLPOOL *pool, SCIP_HEURDATA *heurdata, int *incumbentindex, int *selection, SCIP_Bool randomize)
Definition: heur_rec.c:337
static SCIP_DECL_HEUREXITSOL(heurExitsolRec)
Definition: heur_rec.c:1928
void SCIPrandomPermuteIntArray(SCIP_RANDNUMGEN *randnumgen, int *array, int begin, int end)
Definition: misc.c:9413
SCIP_VAR ** SCIPprobdataGetEdgeVars(SCIP *scip)
#define HEUR_USESSUBSCIP
Definition: heur_rec.c:55
#define HEUR_DESC
Definition: heur_rec.c:48
SCIP_RETCODE SCIPStpHeurRecInitPool(SCIP *scip, STPSOLPOOL **pool, const int nedges, const int maxsize)
Definition: heur_rec.c:894
#define MIN(x, y)
Definition: def.h:209
#define DEFAULT_MAXFREQREC
Definition: heur_rec.c:57
SCIP_RETCODE graph_sol_markPcancestors(SCIP *, IDX **, const int *, const int *, int, STP_Bool *, STP_Bool *, int *, int *, int *)
Definition: grphbase.c:3288
static const unsigned int randseed
Definition: circle.c:46
int SCIPgetNSols(SCIP *scip)
Definition: scip_sol.c:2280
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:116
#define HEUR_DISPCHAR
Definition: heur_rec.c:49
SCIP_Real SCIPgetSolOrigObj(SCIP *scip, SCIP_SOL *sol)
Definition: scip_sol.c:1493
Constraint handler for linear constraints in their most general form, .
#define DEFAULT_RANDSEED
Definition: heur_rec.c:60
static SCIP_RETCODE selectdiffsols(SCIP *scip, const STPSOLPOOL *pool, const GRAPH *graph, SCIP_HEURDATA *heurdata, SCIP_VAR **vars, const int *incumbentedges, int *incumbentindex, int *selection, SCIP_Bool *success)
Definition: heur_rec.c:171
includes various files containing graph methods used for Steiner tree problems
int SCIPparamGetInt(SCIP_PARAM *param)
Definition: paramset.c:716
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:86
int * soledges
Definition: heur_rec.h:41
#define REC_ADDEDGE_FACTOR
Definition: heur_rec.c:71
#define Is_term(a)
Definition: grph.h:168
#define MAX(x, y)
Definition: def.h:208
#define HEUR_MAXDEPTH
Definition: heur_rec.c:53
SCIP_RETCODE SCIPStpHeurRecRun(SCIP *scip, STPSOLPOOL *pool, SCIP_HEUR *heur, SCIP_HEURDATA *heurdata, const GRAPH *graph, SCIP_VAR **vars, int *newsolindex, int runs, int nsols, SCIP_Bool restrictheur, SCIP_Bool *solfound)
Definition: heur_rec.c:1026
#define EAT_FREE
Definition: grph.h:30
SCIP_SOL * SCIPgetBestSol(SCIP *scip)
Definition: scip_sol.c:2379
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
struct SCIP_ProbData SCIP_PROBDATA
Definition: type_prob.h:44
SCIP_PROBDATA * SCIPgetProbData(SCIP *scip)
Definition: scip_prob.c:1020
SCIP_Real * cost
Definition: grph.h:94
#define COST_MAX_POLY_x0
Definition: heur_rec.c:73
SCIP_RETCODE SCIPsetHeurInit(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURINIT((*heurinit)))
Definition: scip_heur.c:264
STPSOL * SCIPStpHeurRecSolfromIdx(STPSOLPOOL *pool, const int soindex)
Definition: heur_rec.c:870
static SCIP_DECL_HEUREXIT(heurExitRec)
Definition: heur_rec.c:1839
#define SCIP_Real
Definition: def.h:150
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:739
Primal recombination heuristic for Steiner problems.
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
SCIP_RETCODE SCIPStpHeurLocalRun(SCIP *scip, GRAPH *graph, const SCIP_Real *cost, int *best_result)
Definition: heur_local.c:208
shortest paths based primal heuristics for Steiner problems
SCIP_RETCODE SCIPStpHeurRecAddToPool(SCIP *scip, const SCIP_Real obj, const int *soledges, STPSOLPOOL *pool, SCIP_Bool *success)
Definition: heur_rec.c:953
#define SCIP_Longint
Definition: def.h:135
int edges
Definition: grph.h:92
#define flipedge(edge)
Definition: grph.h:150
void graph_knot_add(GRAPH *, int)
Definition: grphbase.c:2196
SCIP_RETCODE reduce(SCIP *, GRAPH **, SCIP_Real *, int, int, SCIP_Bool)
Definition: reduce.c:1673
#define BOUND_MAXNTERMINALS
Definition: heur_rec.c:64
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:70
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define UNKNOWN
Definition: sepa_mcf.c:4081
#define STP_RSMT
Definition: grph.h:45
SCIP_RETCODE SCIPsetHeurCopy(SCIP *scip, SCIP_HEUR *heur, SCIP_DECL_HEURCOPY((*heurcopy)))
Definition: scip_heur.c:232
#define nnodes
Definition: gastrans.c:65
void SCIPStpHeurRecFreePool(SCIP *scip, STPSOLPOOL **pool)
Definition: heur_rec.c:925
#define STP_OARSMT
Definition: grph.h:46
#define BMSclearMemoryArray(ptr, num)
Definition: memory.h:112
struct Int_List_Node * parent
Definition: misc_stp.h:77
SCIP_HEURDATA * SCIPheurGetData(SCIP_HEUR *heur)
Definition: heur.c:1165
int hoplimit
Definition: grph.h:89
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1410
SCIP_RETCODE SCIPStpHeurTMPrune(SCIP *scip, const GRAPH *g, const SCIP_Real *cost, int layer, int *result, STP_Bool *connected)
Definition: heur_tm.c:556
void graph_edge_add(SCIP *, GRAPH *, int, int, double, double)
default SCIP plugins
#define STP_RPCSPG
Definition: grph.h:41
int SCIPsolGetIndex(SCIP_SOL *sol)
Definition: sol.c:2584
#define HEUR_TIMING
Definition: heur_rec.c:54
SCIP callable library.
SCIP_RETCODE graph_init(SCIP *, GRAPH **, int, int, int)
Definition: grphbase.c:3491
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:129
int norgmodelknots
Definition: grph.h:60
static SCIP_DECL_HEURINITSOL(heurInitsolRec)
Definition: heur_rec.c:1901
#define HEUR_NAME
Definition: heur_rec.c:47