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