Scippy

GCG

Branch-and-Price & Column Generation for Everyone

heur_setcover.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program */
4 /* GCG --- Generic Column Generation */
5 /* a Dantzig-Wolfe decomposition based extension */
6 /* of the branch-cut-and-price framework */
7 /* SCIP --- Solving Constraint Integer Programs */
8 /* */
9 /* Copyright (C) 2010-2021 Operations Research, RWTH Aachen University */
10 /* Zuse Institute Berlin (ZIB) */
11 /* */
12 /* This program is free software; you can redistribute it and/or */
13 /* modify it under the terms of the GNU Lesser General Public License */
14 /* as published by the Free Software Foundation; either version 3 */
15 /* of the License, or (at your option) any later version. */
16 /* */
17 /* This program is distributed in the hope that it will be useful, */
18 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
19 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
20 /* GNU Lesser General Public License for more details. */
21 /* */
22 /* You should have received a copy of the GNU Lesser General Public License */
23 /* along with this program; if not, write to the Free Software */
24 /* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
25 /* */
26 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
27 
28 /**@file heur_setcover.c
29  * @brief set covering primal heuristic according to Caprara, Fischetti, and Toth (1999)
30  * @author Tobias Oelschlaegel
31  * @author Christian Puchert
32  */
33 
34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35 
36 #include <assert.h>
37 #include <string.h>
38 #include "gcg.h"
39 #include "pricer_gcg.h"
40 #include "scip/clock.h"
41 #include "scip/cons_linear.h"
42 #include "heur_setcover.h"
43 
44 
45 #define HEUR_NAME "setcover"
46 #define HEUR_DESC "set covering primal heuristic according to Caprara, Fischetti, and Toth (1999)"
47 #define HEUR_DISPCHAR 'S'
48 #define HEUR_PRIORITY 0
49 #define HEUR_FREQ 0
50 #define HEUR_FREQOFS 0
51 #define HEUR_MAXDEPTH -1
52 #define HEUR_TIMING SCIP_HEURTIMING_AFTERNODE
53 #define HEUR_USESSUBSCIP FALSE /**< does the heuristic use a secondary SCIP instance? */
54 
55 #define DEF_CORE_TENT_SIZE 10 /**< number of columns covering each row that are added to the tentative core at the beginning */
56 #define DEF_LAMBDA_ADJUSTMENTS TRUE /**< adjust step size during the subgradient phase */
57 #define DEF_LAMBDA_P 50 /**< number of iterations after which lambda is adjusted */
58 #define DEF_LAMBDA 0.1 /**< initial step size for the subgradient phase */
59 #define DEF_STOP_CRIT_ITER 300 /**< number of iterations of the subgradient phase after which the stopping criterion is tested again */
60 #define DEF_STOP_CRIT_DIFF 1.0 /**< stop if absolute difference between best lower and upper bound is less than SCP_STOP_CRIT_DIFF, and */
61 #define DEF_STOP_CRIT_RATIO 0.99 /**< the relative ratio between best lower and upper bound is less than DEF_STOP_CRIT_RATIO */
62 #define DEF_PI_MIN 0.3 /**< percentage of rows to be removed after fixing columns */
63 #define DEF_PI_ALPHA 1.1 /**< increase of pi when no improvement was made, i.e. more columns will be fixed */
64 #define DEF_BETA 1.025 /**< allowed gap between lowerbound and upper bound during the subgradient phase */
65 #define DEF_MAX_ITER 2 /**< maximum number of iterations of three-phase */
66 #define DEF_MAX_ITER_NO_IMP 1 /**< stop of no improvements during the last X iterations of three-phase */
67 #define DEF_THREEPHASE_MAX_ITER 10 /**< maximum number of iterations inside three-phase */
68 #define DEF_GREEDY_MAX_ITER 100 /**< number of multipliers that are used for computing greedy solutions during each iteration */
69 #define DEF_MIN_PROB_SIZE 1000 /**< minimum number of variables the problem has to contain for the heuristic to start */
70 #define DEFAULT_RANDSEED 19 /**< initial random seed */
71 
72 /*
73  * Data structures
74  */
75 
76 /** primal heuristic data */
77 
78 typedef struct
79 {
80  int size; /**< actual number of elements stored in the queue */
81  int reserved; /**< number of elements for which memory is allocated */
82  SCIP_Real* keys; /**< array of keys */
83  int* data; /**< array of values */
84  int** positions; /**< array of pointers to save positions of elements */
85 } PQUEUE;
86 
87 typedef struct
88 {
89  SCIP_Bool* varsfixed; /**< boolean array that indicates which variables are fixed */
90  int nvarsfixed; /**< number of fixed variables */
91  SCIP_Bool* rowscovered; /**< boolean array that indicates which rows are covered by the fixed variables */
92  int nrowscovered; /**< number of rows that are covered by the fixed variables */
93  SCIP_Real costsfixed; /**< total costs of variables that are fixed */
94 } SCP_INSTANCE;
95 
96 typedef struct
97 {
98  SCIP_Bool* varincore; /**< boolean array that indicates whether a variable is in the core */
99  int* listcorevariables; /**< array of indices of core variables */
100  int ncorevariables; /**< number of core variables */
101  SCIP_HASHMAP* mapvariables; /**< maps variables-indices to [0, nvariables) in array 'variables' */
102  SCIP_VAR** variables; /**< all variables of the problem */
103  int* nvarconss; /**< array that contains for each variable the number of constraints it is part of */
104  int nvariables; /**< total number of variables */
105  SCIP_Bool columnsavailable; /**< if set then 'columns' contains the columns for all core variables */
106  int** columns; /**< columns of core variables, NULL if not a core variable */
107  SCIP_Bool rowsavailable; /**< if set then 'rows' contains all rows reduced to core variables */
108  int** rows; /**< rows that only contain core variables */
109  int* nrowvars; /**< number of core variables a row contains */
110  int nconss; /**< total number of constraints (including inactive ones) */
111  int nactiveconss; /**< total number of active constraints for which the variables can be retrieved */
112  int maxconstraintvariables; /**< greatest number of variables some constraint contains */
113  SCIP_CONS** conss; /**< all constraints of the problem */
114  int* constraintid; /**< unique id for each constraint */
115  SCIP_Real* delta; /**< delta values of variables */
116  int* delta_pos;
117  int* solgreedy;
119 } SCP_CORE;
120 
121 typedef struct
122 {
123  SCIP_Bool* xgreedylocal; /**< contains variables that are part of a greedy solution, this is not necessarily a global solution */
124  SCIP_Real* u; /**< lagrange multipliers for the rows */
125  SCIP_Real* subgradient;
126  SCIP_Real* lagrangiancostslocal; /**< lagrangian costs (for a certain instance) when only uncovered rows are considered */
127  SCIP_Real* lagrangiancostsglobal; /**< lagrangian costs for the whole instance when all rows and columns are considered */
128  SCIP_Real ubgreedylocal; /**< bound computed by the greedy set cover algorithm for the restricted instance */
129  SCIP_Real lblagrangelocal; /**< lower bound by lagrange relaxation for the restricted instance */
130  SCIP_Real lblagrangeglobal; /**< lower bound by lagrange relaxation for the unrestricted instance */
132 
133 struct SCIP_HeurData
134 {
135  int param_core_tent_size; /**< number of columns covering each row that are added to the tentative core at the beginning */
136  SCIP_Bool param_lambda_adjustments; /**< adjust step size during the subgradient phase */
137  int param_lambda_p; /**< number of iterations after which lambda is adjusted */
138  SCIP_Real param_lambda; /**< initial step size for the subgradient phase */
139  int param_stop_crit_iter; /**< number of iterations of the subgradient phase after which the stopping criterion is tested again */
140  SCIP_Real param_stop_crit_diff; /**< stop if absolute difference between best lower and upper bound is less than SCP_STOP_CRIT_DIFF, and */
141  SCIP_Real param_stop_crit_ratio; /**< the relative gap between best lower and upper bound is less than (1 - SCP_STOP_CRIT_PER */
142  SCIP_Real param_pi_min; /**< percentage of rows to be removed after fixing columns */
143  SCIP_Real param_pi_alpha; /**< increase of pi when no improvement was made, i.e. more columns will be fixed */
144  SCIP_Real param_beta; /**< allowed a gap between lowerbound and upper bound during the subgradient phase */
145  int param_max_iter; /**< maximum number of iterations of three-phase */
146  int param_max_iter_no_imp; /**< stop of no improvements during the last X iterations of three-phase */
147  int param_threephase_max_iter; /**< maximum number of iterations inside three-phase */
148  int param_greedy_max_iter; /**< number of multipliers that are used for computing greedy solutions during each iteration */
149  int param_min_prob_size; /**< minimum number of variables the master problem needs to contain before the heuristic starts at all */
150 
151  SCP_CORE core; /**< core (subcollection of columns) of the problem covering all rows */
152  SCP_INSTANCE inst; /**< reduced instance where some variables may be fixed and some rows be covered */
153  SCP_INSTANCE subinst; /**< reduced instance of 'inst', used during the three-phase */
154 
155  SCP_Lagrange_Sol multbestlbtotal; /**< lagrange multiplier that gives the best lower bound for the complete problem */
156  SCP_Lagrange_Sol multbestlbinst; /**< best multiplier for instance 'inst' */
157  SCP_Lagrange_Sol multbestlbsubinst; /**< best multiplier for instance 'subinst' */
158 
159  SCIP_Real bestub; /**< best upper bound that could be obtained so far */
160  SCIP_Bool* bestubsol; /**< actual solution that gives the best upper bound */
161  SCIP_Real bestubinst; /**< best upper bound for the reduced instance 'inst' (including fixed costs) */
162  SCIP_Bool* bestubinst_sol; /**< actual solution for the instance 'inst' (including fixed variables) */
163  SCIP_Real bestubsubinst; /**< best upper bound for the reduced instance 'subinst' (including fixed costs) */
164  SCIP_Bool* bestubsubinstsol; /**< actual solution for the instance 'subinst' (including fixed variables) */
165 
166  /* local data that is used by (most) local procedures */
167  SCIP_VAR** vars; /**< used to iterate through the variables of a constraint */
168 
169  /* memory that is used locally by 'threephase' */
170  SCP_Lagrange_Sol tpmultlbsubinst; /**< lagrange multiplier */
171  SCIP_Bool useinitialmultiplier; /**< compute an own initial lagrange multiplier instead of using the best known */
172 
173  /* memory that is used locally by 'greedysetcover' */
174  PQUEUE greedyqueue; /**< priority queue used by the greedy algorithm */
175  int* greedycolpos; /**< stores position of a variable within the priority queue */
176  int* greedycolmu; /**< value mu for each variable */
177  SCIP_Real* greedycolgamma; /**< value gamma for each variable */
178  SCIP_Real* greedycolscore; /**< score of each variable */
179  SCP_INSTANCE greedyinst; /**< instance that is used to hold covered rows */
180 
181  /* memory that is used locally by redefineCore */
182  int* rccols; /**< stores columns covering some row sorted by increasing costs */
183  SCIP_Real* rccoldelta; /**< delta values of the columns stored in 'rccols' */
184 
185  /* memory that is used locally by subgradientOptimization */
186  SCIP_Real* sglastlb; /**< stores lower bounds of the last iterations */
187  SCIP_RANDNUMGEN* randnumgen; /**< random number generator */
188 };
189 
190 /*
191  * Local methods
192  */
193 
194 /** initializes a priority queue */
195 static
196 SCIP_RETCODE pqueue_init(
197  SCIP* scip, /**< master SCIP data structure */
198  PQUEUE* queue /**< priority queue instance */
199  )
200 {
201  queue->size = 0;
202  queue->reserved = 32;
203 
204  SCIP_CALL( SCIPallocBufferArray(scip, &queue->keys, queue->reserved) );
205  SCIP_CALL( SCIPallocBufferArray(scip, &queue->data, queue->reserved) );
206  SCIP_CALL( SCIPallocBufferArray(scip, &queue->positions, queue->reserved) );
207 
208  return SCIP_OKAY;
209 }
210 
211 /** removes all elements from a priority queue, but does not release any memory */
212 static
214  SCIP* scip, /**< master SCIP data structure */
215  PQUEUE* queue /**< priority queue instance */
216  )
217 {
218  queue->size = 0;
219 }
220 
221 /** releases all memory that is used by the priority queue */
222 static
224  SCIP* scip, /**< master SCIP data structure */
225  PQUEUE* queue /**< priority queue instance */
226  )
227 {
228  SCIPfreeBufferArrayNull(scip, &queue->positions);
229  SCIPfreeBufferArrayNull(scip, &queue->data);
230  SCIPfreeBufferArrayNull(scip, &queue->keys);
231 }
232 
233 /** inserts an element with key 'key' and value 'elem' into the queue.
234  * If 'position' is not NULL, the referenced memory location will always contain the internal position of the element
235  * */
236 static
237 SCIP_RETCODE pqueue_insert(
238  SCIP* scip, /**< master SCIP data structure */
239  PQUEUE* queue, /**< priority queue instance */
240  SCIP_Real key, /**< key of the new item */
241  int elem, /**< value of the new item */
242  int* position /**< address where the items position within the queue is stored */
243  )
244 {
245  int pos;
246  int parent = 0;
247 
248  if( queue->size == queue->reserved )
249  {
250  queue->reserved += 128;
251  SCIP_CALL( SCIPreallocBufferArray(scip, &queue->keys, queue->reserved) );
252  SCIP_CALL( SCIPreallocBufferArray(scip, &queue->data, queue->reserved) );
253  SCIP_CALL( SCIPreallocBufferArray(scip, &queue->positions, queue->reserved) );
254  }
255 
256  pos = queue->size++;
257  queue->keys[pos] = key;
258  queue->data[pos] = elem;
259 
260  if( pos > 0 )
261  parent = (pos - 1) / 2;
262 
263  while( (pos > 0) && (key < queue->keys[parent]) )
264  {
265  /* swap element with parent */
266  queue->keys[pos] = queue->keys[parent];
267  queue->data[pos] = queue->data[parent];
268  queue->positions[pos] = queue->positions[parent];
269  /*queue->data[parent] = elem;*/
270 
271  if( queue->positions[pos] != NULL )
272  *queue->positions[pos] = pos;
273 
274  pos = (pos - 1) / 2;
275  parent = (pos - 1) / 2;
276  }
277 
278  queue->keys[pos] = key;
279  queue->data[pos] = elem;
280  queue->positions[pos] = position;
281 
282  if( position != NULL )
283  *position = pos;
284 
285  return SCIP_OKAY;
286 }
287 
288 /** decreases the key to 'key' of the element that is currently at position 'pos' */
289 static
291  SCIP* scip, /**< master SCIP data structure */
292  PQUEUE* queue, /**< priority queue instance */
293  int pos, /**< position of the item within the queue that is to be changed */
294  SCIP_Real key /**< new key of the item, needs to be greater than the old key */
295  )
296 {
297  int parent;
298  int* positionptr;
299  int elem;
300 
301  if( pos >= queue->size )
302  return;
303 
304  parent = (pos - 1) / 2;
305  positionptr = queue->positions[pos];
306  elem = queue->data[pos];
307 
308  while( (pos > 0) && (key < queue->keys[parent]) )
309  {
310  /* swap element with parent */
311  queue->keys[pos] = queue->keys[parent];
312  queue->data[pos] = queue->data[parent];
313  queue->positions[pos] = queue->positions[parent];
314  queue->data[parent] = elem;
315 
316  if( queue->positions[pos] != NULL )
317  *queue->positions[pos] = pos;
318 
319  pos = (pos - 1) / 2;
320  parent = (pos - 1) / 2;
321  }
322 
323  queue->keys[pos] = key;
324  queue->data[pos] = elem;
325  queue->positions[pos] = positionptr;
326 
327  if( positionptr != NULL )
328  *positionptr = pos;
329 }
330 
331 /** increases the key to 'key' of the element that is currently at position 'pos' */
332 static
334  SCIP* scip, /**< master SCIP data structure */
335  PQUEUE* queue, /**< priority queue instance */
336  int pos, /**< position of the item within the queue that is to be changed */
337  SCIP_Real key /**< new key of the item, needs to be smaller than the old one */
338  )
339 {
340  int elem;
341  int* positionptr;
342  int left;
343  int right;
344 
345  left = 2 * pos + 1;
346  right = 2 * pos + 2;
347 
348  if( pos >= queue->size )
349  return;
350 
351  elem = queue->data[pos];
352  positionptr = queue->positions[pos];
353 
354  while( left < queue->size )
355  {
356  if( right < queue->size )
357  {
358  /* both children exist, so swap with smallest child */
359  if( key <= queue->keys[left] )
360  {
361  /* left child is not smaller than element */
362  if( key <= queue->keys[right] )
363  break;
364  else
365  {
366  /* swap with right child */
367  queue->keys[pos] = queue->keys[right];
368  queue->data[pos] = queue->data[right];
369  queue->positions[pos] = queue->positions[right];
370 
371  if( queue->positions[pos] != NULL )
372  *(queue->positions[pos]) = pos;
373 
374  pos = right;
375  }
376  }
377  else
378  {
379  /* left child is smaller than the element */
380  if( queue->keys[left] <= queue->keys[right] )
381  {
382  /* left child is smaller than right child and smaller than the element */
383  queue->keys[pos] = queue->keys[left];
384  queue->data[pos] = queue->data[left];
385  queue->positions[pos] = queue->positions[left];
386 
387  if( queue->positions[pos] != NULL )
388  *(queue->positions[pos]) = pos;
389 
390  pos = left;
391  }
392  else
393  {
394  /* right child is smallest element */
395  queue->keys[pos] = queue->keys[right];
396  queue->data[pos] = queue->data[right];
397  queue->positions[pos] = queue->positions[right];
398 
399  if( queue->positions[pos] != NULL )
400  *(queue->positions[pos]) = pos;
401 
402  pos = right;
403  }
404  }
405  }
406  else
407  {
408  /* only left child exists */
409  if( key > queue->keys[left] )
410  {
411  /* left child is smaller than element, so swap them */
412  queue->keys[pos] = queue->keys[left];
413  queue->data[pos] = queue->data[left];
414  queue->positions[pos] = queue->positions[left];
415 
416  if( queue->positions[pos] != NULL )
417  *(queue->positions[pos]) = pos;
418 
419  pos = left;
420  }
421  else
422  break;
423  }
424 
425  left = 2 * pos + 1;
426  right = 2 * pos + 2;
427  }
428 
429  queue->keys[pos] = key;
430  queue->data[pos] = elem;
431  queue->positions[pos] = positionptr;
432 
433  if( positionptr != NULL )
434  *positionptr = pos;
435 }
436 
437 /** returns the value of a minimum element in 'elem'. Sets 'elem' to -1 if the queue is empty */
438 static
440  SCIP* scip, /**< master SCIP data structure */
441  PQUEUE* queue, /**< priority queue instance */
442  int* elem /**< address where the value of a minimum key item will be stored */
443  )
444 {
445  *elem = -1;
446 
447  if( queue->size > 0 )
448  {
449  *elem = queue->data[0];
450 
451  queue->size--;
452 
453  /* copy last element to front and repair heap structure */
454  if( queue->size > 0 )
455  {
456  queue->keys[0] = queue->keys[queue->size];
457  queue->data[0] = queue->data[queue->size];
458  queue->positions[0] = queue->positions[queue->size];
459 
460  if( queue->positions[0] != NULL )
461  *(queue->positions[0]) = 0;
462 
463  pqueue_increase_key(scip, queue, 0, queue->keys[0]);
464  }
465  }
466 }
467 
468 /** allocates memory for a lagrange multiplier and a set covering solution */
469 static
471  SCIP* scip, /**< master SCIP data structure */
472  SCP_CORE* core, /**< core data structure */
473  SCP_Lagrange_Sol* mult /**< uninitialized solution data structure */
474  )
475 {
476  int i;
477 
478  SCIP_CALL( SCIPallocBufferArray(scip, &mult->u, core->nconss) );
479  SCIP_CALL( SCIPallocBufferArray(scip, &mult->subgradient, core->nconss) );
480  SCIP_CALL( SCIPallocBufferArray(scip, &mult->lagrangiancostslocal, core->nvariables) );
481  SCIP_CALL( SCIPallocBufferArray(scip, &mult->lagrangiancostsglobal, core->nvariables) );
482  SCIP_CALL( SCIPallocClearBufferArray(scip, &mult->xgreedylocal, SCIPgetNVars(scip)) );
483 
484  for( i = 0; i < core->nvariables; i++ )
485  {
486  mult->lagrangiancostslocal[i] = 0.0;
487  mult->lagrangiancostsglobal[i] = 0.0;
488  }
489 
490  for( i = 0; i < core->nconss; i++ )
491  {
492  mult->u[i] = 0.0;
493  mult->subgradient[i] = 0.0;
494  }
495 
496  mult->lblagrangeglobal = SCIP_REAL_MIN;
497  mult->lblagrangelocal = SCIP_REAL_MIN;
498  mult->ubgreedylocal = SCIP_REAL_MAX;
499 
500  return SCIP_OKAY;
501 }
502 
503 /** checks if 'var' is a core variable */
504 static
505 SCIP_Bool isCoreVariable(
506  SCP_CORE* core, /**< initialized core data structure */
507  SCIP_VAR* var /**< SCIP variable */
508  )
509 {
510  int varidx;
511 
512  assert(core != NULL);
513  assert(var != NULL);
514 
515  varidx = SCIPvarGetProbindex(var);
516  assert(varidx != -1);
517 
518  return core->varincore[varidx];
519 }
520 
521 /** checks if 'var' is fixed within instance 'inst' */
522 static
523 SCIP_Bool isFixedVariable(
524  SCP_INSTANCE* inst, /**< instance of the set covering problem */
525  SCIP_VAR* var /**< SCIP variable */
526  )
527 {
528  int varidx;
529 
530  assert(inst != NULL);
531  assert(var != NULL);
532 
533  varidx = SCIPvarGetProbindex(var);
534  assert(varidx != -1);
535 
536  return inst->varsfixed[varidx];
537 }
538 
539 /** checks if 'variable' is part of 'solution' */
540 static
541 SCIP_Bool isVarInSolution(
542  SCIP_Bool* solution, /**< SCIP hashtable that contains SCIP variables */
543  SCIP_VAR* var /**< SCIP variable */
544  )
545 {
546  int varidx;
547 
548  assert(solution != NULL);
549  assert(var != NULL);
550 
551  varidx = SCIPvarGetProbindex(var);
552  assert(varidx != -1);
553 
554  return solution[varidx];
555 }
556 
557 /** adds a variable to the core */
558 static
560  SCP_CORE* core, /**< initialized core data structure */
561  SCIP_VAR* var /**< SCIP variable */
562  )
563 {
564  int varidx;
565 
566  assert(core != NULL);
567  assert(var != NULL);
568 
569  varidx = SCIPvarGetProbindex(var);
570  assert(varidx != -1);
571 
572  if( !core->varincore[varidx] )
573  ++core->ncorevariables;
574  core->varincore[varidx] = TRUE;
575 }
576 
577 /** fixes 'var' within instance 'inst' */
578 static
580  SCP_INSTANCE* inst, /**< SCP instance where 'variable' is to be fixed */
581  SCIP_VAR* var /**< SCIP variable */
582  )
583 {
584  int varidx;
585 
586  assert(inst != NULL);
587  assert(var != NULL);
588 
589  varidx = SCIPvarGetProbindex(var);
590  assert(varidx != -1);
591 
592  if( !inst->varsfixed[varidx] )
593  ++inst->nvarsfixed;
594  inst->varsfixed[varidx] = TRUE;
595 }
596 
597 /** adds variable 'var' to 'solution' */
598 static
600  SCIP_Bool* solution, /**< solution represented as a boolean array */
601  SCIP_VAR* var /**< SCIP variable */
602  )
603 {
604  int varidx;
605 
606  assert(solution != NULL);
607  assert(var != NULL);
608 
609  varidx = SCIPvarGetProbindex(var);
610  assert(varidx != -1);
611 
612  /* no test whether the variable is already part of the solution */
613  solution[varidx] = TRUE;
614 }
615 
616 /** removes all variables from the core */
617 static
619  SCIP* scip, /**< SCIP data structure */
620  SCP_CORE* core /**< initialized core data structure */
621  )
622 {
623  int i;
624 
625  assert(core != NULL);
626 
627  for( i = 0; i < SCIPgetNVars(scip); ++i )
628  core->varincore[i] = FALSE;
629  core->ncorevariables = 0;
630 }
631 
632 /** removes a 'var' from 'solution' */
633 static
635  SCIP_Bool* solution, /**< solution represented as a boolean array */
636  SCIP_VAR* var /**< SCIP variable */
637  )
638 {
639  int varidx;
640 
641  assert(solution != NULL);
642  assert(var != NULL);
643 
644  varidx = SCIPvarGetProbindex(var);
645  assert(varidx != -1);
646 
647  /* no test whether the variable is already part of the solution */
648  solution[varidx] = FALSE;
649 }
650 
651 /** clears a solution */
652 static
654  SCIP* scip, /**< SCIP data structure */
655  SCIP_Bool* solution /**< solution represented as a boolean array */
656  )
657 {
658  int i;
659 
660  assert(scip != NULL);
661  assert(solution != NULL);
662 
663  for( i = 0; i < SCIPgetNVars(scip); ++i )
664  solution[i] = FALSE;
665 }
666 
667 /** checks if the row at position 'rowpos' is covered by fixed variables of 'inst' */
668 static
669 SCIP_Bool isRowCovered(
670  SCP_INSTANCE* inst, /**< SCP instance */
671  int rowpos /**< position of the row within 'core' */
672  )
673 {
674  assert(inst != NULL);
675 
676  return inst->rowscovered[rowpos];
677 }
678 
679 /** marks the row at position 'rowpos' as covered within instance 'inst' */
680 static
682  SCP_INSTANCE* inst, /**< SCP instance */
683  int rowpos /**< position of the row within 'core' */
684  )
685 {
686  assert(inst != NULL);
687 
688  if( !inst->rowscovered[rowpos] )
689  ++inst->nrowscovered;
690  inst->rowscovered[rowpos] = TRUE;
691 }
692 
693 /** returns the position of 'variable' within the array core->variables */
694 static
695 SCIP_RETCODE getVarIndex(
696  SCIP* scip, /**< master SCIP data structure */
697  SCP_CORE* core, /**< SCP core data structure */
698  SCIP_VAR* variable, /**< SCIP variable */
699  int* pos /**< address where the position is to be stored */
700  )
701 {
702  int varidx;
703 
704  assert(scip != NULL);
705  assert(core != NULL);
706  assert(variable != NULL);
707  assert(pos != NULL);
708 
709  varidx = SCIPvarGetIndex(variable);
710  *pos = (int) (size_t) SCIPhashmapGetImage(core->mapvariables, (void *) ((size_t) (varidx + 1))); /*lint !e507*/
711 
712  return SCIP_OKAY;
713 }
714 
715 /** get all variables that are part of the constraint at position 'pos' (within SCIP) and saves them into 'vars' */
716 static
717 SCIP_RETCODE getConsVars(
718  SCIP* scip, /**< master SCIP data structure */
719  SCP_CORE* core, /**< SCP core data structure */
720  int pos, /**< position of the constraint within SCIP */
721  SCIP_VAR** vars, /**< allocated memory that can hold all variables */
722  int* nvars, /**< address where the number of variables will be stored */
723  SCIP_Bool* success /**< address where the result of the operation will be stored */
724  )
725 {
726  *success = FALSE;
727 
728  if( SCIPconsIsActive(core->conss[pos]) == FALSE )
729  return SCIP_OKAY;
730 
731  SCIP_CALL( SCIPgetConsNVars(scip, core->conss[pos], nvars, success) );
732 
733  if( *success == FALSE )
734  return SCIP_OKAY;
735 
736  SCIP_CALL( SCIPgetConsVars(scip, core->conss[pos], vars, core->maxconstraintvariables, success) );
737 
738  return SCIP_OKAY;
739 }
740 
741 /** releases all memory of a lagrange multiplier */
742 static
744  SCIP* scip, /**< master SCIP data structure */
745  SCP_Lagrange_Sol* mult /**< initialized SCP lagrange multiplier that is to be freed */
746  )
747 {
748  SCIPfreeBufferArray(scip, &mult->xgreedylocal);
749  SCIPfreeBufferArray(scip, &mult->lagrangiancostsglobal);
750  SCIPfreeBufferArray(scip, &mult->lagrangiancostslocal);
751  SCIPfreeBufferArray(scip, &mult->subgradient);
752  SCIPfreeBufferArray(scip, &mult->u);
753 }
754 
755 /** creates a set covering solution. Adds all fixed variables of 'inst' and all variables of 'source' to 'dest'.
756  * 'destcosts' contains the total costs of the solution.
757  */
758 static
760  SCIP* scip, /**< master SCIP data structure */
761  SCP_CORE* core, /**< SCP core data structure */
762  SCP_INSTANCE* inst, /**< SCP instance */
763  SCIP_Bool* dest, /**< array where variables belonging to the solution are inserted */
764  SCIP_Real* destcosts, /**< address where the total costs of the solution will be stored */
765  SCIP_Bool* source /**< SCP solution to the (reduced) instance 'inst' */
766  )
767 {
768  int i;
769  SCIP_Real costs = 0.0;
770 
771  /* remove all entries from dest */
772  clearSolution(scip, dest);
773 
774  /* iterate through all variables of the problem */
775  for( i = 0; i < core->nvariables; i++ )
776  {
777  if( isVarInSolution(source, core->variables[i]) || isFixedVariable(inst, core->variables[i]) )
778  {
779  /* add variable 'i' if it is in the solution 'source' or fixed in instance 'inst' */
780  costs += SCIPvarGetObj(core->variables[i]);
781  addVarToSolution(dest, core->variables[i]);
782  }
783  }
784 
785  /* save total costs if required */
786  if( destcosts != NULL )
787  *destcosts = costs;
788 }
789 
790 /** copies all data of the lagrange multiplier 'source' to the lagrange multiplier 'dest' */
791 static
793  SCIP* scip, /**< master SCIP data structure */
794  SCP_CORE* core, /**< SCP core data structure */
795  SCP_Lagrange_Sol* dest, /**< SCP lagrange multiplier where data will be copied to */
796  SCP_Lagrange_Sol* source /**< SCP lagrange multiplier where data will be copied from */
797  )
798 {
799  int i;
800 
801  clearSolution(scip, dest->xgreedylocal);
802 
803  memcpy(dest->lagrangiancostslocal, source->lagrangiancostslocal, sizeof(*dest->lagrangiancostslocal) * core->nvariables);
804  memcpy(dest->lagrangiancostsglobal, source->lagrangiancostsglobal, sizeof(*dest->lagrangiancostsglobal) * core->nvariables);
805  memcpy(dest->u, source->u, sizeof(*dest->u) * core->nconss);
806  memcpy(dest->subgradient, source->subgradient, sizeof(*dest->subgradient) * core->nconss);
807 
808  for( i = 0; i < core->nvariables; i++ )
809  {
810  if( isVarInSolution(source->xgreedylocal, core->variables[i]) )
811  addVarToSolution(dest->xgreedylocal, core->variables[i]);
812  }
813 
814  dest->lblagrangeglobal = source->lblagrangeglobal;
815  dest->lblagrangelocal = source->lblagrangelocal;
816  dest->ubgreedylocal = source->ubgreedylocal;
817 }
818 
819 /** initializes an instance where no variables are fixed and no rows are covered */
820 static
821 SCIP_RETCODE initInstance(
822  SCIP* scip, /**< master SCIP data structure */
823  SCP_INSTANCE* inst /**< unitialized SCP instance */
824  )
825 {
826  SCIP_CALL( SCIPallocClearBufferArray(scip, &inst->varsfixed, SCIPgetNVars(scip)) );
827  SCIP_CALL( SCIPallocClearBufferArray(scip, &inst->rowscovered, SCIPgetNConss(scip)) );
828  inst->nvarsfixed = 0;
829  inst->nrowscovered = 0;
830  inst->costsfixed = 0.0;
831 
832  return SCIP_OKAY;
833 }
834 
835 /** copies the fixed variables from 'source' to 'dest', but not the set of covered rows. this must be done separately. */
836 static
838  SCIP* scip, /**< master SCIP data structure */
839  SCP_CORE* core, /**< SCP core data structure */
840  SCP_INSTANCE* dest, /**< initialized SCP instance */
841  SCP_INSTANCE* source /**< initialized SCP instance that is to be copied */
842  )
843 {
844  int i;
845 
846  for( i = 0; i < SCIPgetNVars(scip); ++i )
847  dest->varsfixed[i] = FALSE;
848  dest->nvarsfixed = 0;
849  for( i = 0; i < SCIPgetNConss(scip); ++i )
850  dest->rowscovered[i] = FALSE;
851  dest->nrowscovered = 0;
852  dest->costsfixed = 0.0;
853 
854  for( i = 0; i < core->nvariables; i++ )
855  {
856  if( isFixedVariable(source, core->variables[i]) == TRUE )
857  {
858  fixVariable(dest, core->variables[i]);
859  dest->costsfixed += SCIPvarGetObj(core->variables[i]);
860  }
861  }
862 }
863 
864 /** releases all memory used by an instance */
865 static
867  SCIP* scip, /**< master SCIP data structure */
868  SCP_INSTANCE* inst /**< initialized SCP instance that is to be freed */
869  )
870 {
871  SCIPfreeBufferArray(scip, &inst->rowscovered);
872  SCIPfreeBufferArray(scip, &inst->varsfixed);
873 }
874 
875 /** initializes a tentative core: for each row the first few columns covering this row are added to the core */
876 static
877 SCIP_RETCODE initTentativeCore(
878  SCIP* scip, /**< master SCIP data structure */
879  SCP_CORE* core, /**< SCP core data structure */
880  SCIP_HEURDATA* heurdata /**< pointer to SCIP's heurdata */
881  )
882 {
883  SCIP_VAR** vars;
884  SCIP_Bool success;
885  int nvars;
886  int i;
887  int j;
888 
889  assert(scip != NULL);
890  assert(core != NULL);
891 
892  /* mapvariables is a mapping: SCIPvarGetIndex(core->variables[i]) -> i */
893  SCIP_CALL( SCIPhashmapCreate(&core->mapvariables, SCIPblkmem(scip), SCIPgetNVars(scip)) );
894  SCIP_CALL( SCIPallocClearBufferArray(scip, &core->varincore, SCIPgetNVars(scip)) );
895 
896  core->rowsavailable = FALSE;
897  core->rows = NULL;
898  core->nrowvars = NULL;
899  core->columnsavailable = FALSE;
900  core->columns = NULL;
901  core->nvariables = SCIPgetNVars(scip);
902  core->variables = SCIPgetVars(scip);
903 
904  SCIP_CALL( SCIPallocBufferArray(scip, &core->delta, core->nvariables) );
905  SCIP_CALL( SCIPallocBufferArray(scip, &core->delta_pos, core->nvariables) );
906  SCIP_CALL( SCIPallocBufferArray(scip, &core->solgreedy, core->nvariables) );
907  core->nsolgreedy = 0;
908 
909  SCIP_CALL( SCIPallocBufferArray(scip, &core->nvarconss, core->nvariables) );
910  for( i = 0; i< core->nvariables; i++ )
911  core->nvarconss[i] = 0;
912 
913  /* construct mapping of variable-indices to array 'variables' */
914  for( i = 0; i < core->nvariables; i++ )
915  {
916  int varidx = SCIPvarGetIndex(core->variables[i]);
917  SCIP_CALL( SCIPhashmapInsert(core->mapvariables, (void *) ((size_t) (varidx + 1)), (void *) (size_t) i) );
918  }
919 
920  core->nactiveconss = 0;
921  core->maxconstraintvariables = 0;
922  core->nconss = SCIPgetNConss(scip);
923  core->conss = SCIPgetConss(scip);
924  core->ncorevariables = 0;
925  core->listcorevariables = NULL;
926  SCIP_CALL( SCIPallocBufferArray(scip, &core->constraintid, core->nconss) );
927 
928  SCIP_CALL( SCIPallocBufferArray(scip, &vars, core->nvariables) );
929 
930  for( i = 0; i < core->nconss; i++ )
931  {
932  core->constraintid[i] = i;
933  if( SCIPconsIsActive(core->conss[i]) == FALSE )
934  {
935  SCIPdebugMessage("constraint %i (%s) is inactive\n", i, SCIPconsGetName(core->conss[i]));
936  continue;
937  }
938 
939  /* get all variables that are part of this constraint */
940  SCIP_CALL( SCIPgetConsNVars(scip, core->conss[i], &nvars, &success) );
941  if( success == FALSE )
942  {
943  SCIPdebugMessage("constraint %i (%s): can't get number of variables\n", i, SCIPconsGetName(core->conss[i]));
944  continue;
945  }
946 
947  SCIP_CALL( SCIPgetConsVars(scip, core->conss[i], vars, core->nvariables, &success) );
948  if( success == FALSE )
949  {
950  SCIPdebugMessage("constraint %i (%s): can't get variables\n", i, SCIPconsGetName(core->conss[i]));
951  continue;
952  }
953 
954  if( nvars > core->maxconstraintvariables )
955  core->maxconstraintvariables = nvars;
956 
957  for( j = 0; j < nvars; j++ )
958  {
959  int varpos;
960 
961  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
962  if( j < heurdata->param_core_tent_size )
963  {
964  /* add this variable to the core if it's not already in there */
965  if( !isCoreVariable(core, core->variables[varpos]) )
966  addVarToCore(core, core->variables[varpos]);
967  }
968 
969  /* increase the number of constraints this variable is part of */
970  core->nvarconss[varpos]++;
971  }
972 
973  core->nactiveconss++;
974  }
975 
976  SCIPfreeBufferArray(scip, &vars);
977 
978  /* create list of core variables, so it is easy to traverse them */
979  j = 0;
980  SCIP_CALL( SCIPallocBufferArray(scip, &core->listcorevariables, core->ncorevariables) );
981  for( i = 0; i < core->nvariables; i++ )
982  {
983  if( isCoreVariable(core, core->variables[i]) )
984  core->listcorevariables[j++] = i;
985  }
986 
987  SCIPdebugMessage("%d variables in the tentative core\n", core->ncorevariables);
988 
989  return SCIP_OKAY;
990 }
991 
992 /** adds all fixed variables of 'inst' to a set covering solution 'solution' */
993 static
995  SCIP* scip, /**< master SCIP data structure */
996  SCP_CORE* core, /**< SCP core data structure */
997  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
998  SCIP_Bool* solution /**< SCP solution (a hashtable that contains SCIP variables) */
999  )
1000 {
1001  int i;
1002 
1003  for( i = 0; i < core->nvariables; i++ )
1004  {
1005  if( !isFixedVariable(inst, core->variables[i]) )
1006  continue;
1007 
1008  /* add fixed variable if it is not already part of the solution */
1009  if( !isVarInSolution(solution, core->variables[i]) )
1010  addVarToSolution(solution, core->variables[i]);
1011  }
1012 }
1013 
1014 /** constructs rows of all constraints, but only includes core variables */
1015 static
1016 SCIP_RETCODE computeCoreRows(
1017  SCIP* scip, /**< master SCIP data structure */
1018  SCP_CORE* core, /**< SCP core data structure */
1019  SCIP_HEURDATA* heurdata /**< pointer to the heuristics's data */
1020  )
1021 {
1022  SCIP_Bool success;
1023  SCIP_VAR** vars;
1024  int i;
1025  int j;
1026  int ncorevars;
1027  int nvars;
1028 
1029  assert(scip != NULL);
1030  assert(core != NULL);
1031  assert(heurdata != NULL);
1032 
1033  /* don't compute again if already computed */
1034  if( core->rowsavailable )
1035  return SCIP_OKAY;
1036 
1037  assert(core->rows == NULL);
1038 
1039  SCIP_CALL( SCIPallocBufferArray(scip, &core->rows, core->nconss) );
1040  SCIP_CALL( SCIPallocBufferArray(scip, &core->nrowvars, core->nconss) );
1041  SCIP_CALL( SCIPallocBufferArray(scip, &vars, core->maxconstraintvariables) );
1042 
1043  /* iterate through list of constraints */
1044  for( i = 0; i < core->nconss; i++ )
1045  {
1046  core->rows[i] = NULL;
1047  core->nrowvars[i] = 0;
1048 
1049  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1050  if( success == FALSE )
1051  continue;
1052 
1053  /* count number of core variables this constraint contains */
1054  ncorevars = 0;
1055  for( j = 0; j < nvars; j++ )
1056  {
1057  int varpos;
1058 
1059  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1060 
1061  if( isCoreVariable(core, core->variables[varpos]) )
1062  ncorevars++;
1063  }
1064 
1065  if( ncorevars > 0 )
1066  {
1067  /* create array of core variables of this constraint */
1068  core->nrowvars[i] = ncorevars;
1069  SCIP_CALL( SCIPallocBufferArray(scip, &(core->rows[i]), core->nrowvars[i]) ); /*lint !e866*/
1070  for( j = 0; j < nvars; j++ )
1071  {
1072  int varpos;
1073 
1074  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1075 
1076  if( isCoreVariable(core, core->variables[varpos]) )
1077  {
1078  ncorevars--;
1079  core->rows[i][ncorevars] = varpos;
1080  }
1081  }
1082  }
1083  }
1084 
1085  SCIPfreeBufferArray(scip, &vars);
1086  core->rowsavailable = TRUE;
1087  return SCIP_OKAY;
1088 }
1089 
1090 /** constructs columns of core variables to provide faster access */
1091 static
1092 SCIP_RETCODE computeCoreColumns(
1093  SCIP* scip, /**< master SCIP data structure */
1094  SCP_CORE* core /**< SCP core data structure */
1095  )
1096 {
1097  SCIP_Bool success;
1098  SCIP_VAR** vars;
1099  int i;
1100  int j;
1101  int k;
1102  int nvars;
1103 
1104  assert(scip != NULL);
1105  assert(core != NULL);
1106 
1107  /* don't compute columns again if already computed */
1108  if( core->columnsavailable )
1109  return SCIP_OKAY;
1110 
1111  SCIP_CALL( SCIPallocBufferArray(scip, &core->columns, core->nvariables) );
1112  for( i = 0; i < core->nvariables; i++ )
1113  {
1114  /* columns[i] = NULL for all non-core variables */
1115  core->columns[i] = NULL;
1116 
1117  /* check if variable is a core variable */
1118  if( !isCoreVariable(core, core->variables[i]) )
1119  continue;
1120 
1121  /* only allocate memory of it is part of any constraints at all (this should always be the case for core variables) */
1122  if( core->nvarconss[i] == 0 )
1123  continue;
1124 
1125  SCIP_CALL( SCIPallocBufferArray(scip, &(core->columns[i]), core->nvarconss[i]) ); /*lint !e866*/
1126 
1127  for( j = 0; j < core->nvarconss[i]; j++ )
1128  core->columns[i][j] = -1;
1129  }
1130 
1131  SCIP_CALL( SCIPallocBufferArray(scip, &vars, core->maxconstraintvariables) );
1132 
1133  for( i = 0; i < core->nconss; i++ )
1134  {
1135  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1136  if( success == FALSE )
1137  continue;
1138 
1139  for( j = 0; j < nvars; j++ )
1140  {
1141  int varpos;
1142 
1143  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1144 
1145  if( !isCoreVariable(core, core->variables[varpos]) )
1146  continue;
1147 
1148  /* add this constraint to the column of the variable */
1149  for( k = 0; k < core->nvarconss[varpos]; k++ )
1150  {
1151  /* find first position that is not used yet in column 'varpos' */
1152  if( core->columns[varpos][k] == -1 )
1153  {
1154  core->columns[varpos][k] = i;
1155  break;
1156  }
1157  }
1158  }
1159  }
1160 
1161  SCIPfreeBufferArray(scip, &vars);
1162  core->columnsavailable = TRUE;
1163  return SCIP_OKAY;
1164 }
1165 
1166 /** releases memory that is used by the core, including rows and columns, if they were computed */
1167 static
1169  SCIP* scip, /**< master SCIP data structure */
1170  SCP_CORE* core /**< SCP core data structure that is to be freed */
1171  )
1172 {
1173  assert(core != NULL);
1174  assert(core->varincore != NULL);
1175  assert(core->mapvariables != NULL);
1176  assert(core->nvarconss != NULL);
1177 
1178  if( core->rowsavailable )
1179  {
1180  int i;
1181 
1182  for( i = 0; i < core->nconss; i++ )
1183  if( core->rows[i] )
1184  SCIPfreeBufferArray(scip, &core->rows[i]);
1185 
1186  SCIPfreeBufferArray(scip, &core->nrowvars);
1187  SCIPfreeBufferArray(scip, &core->rows);
1188  }
1189 
1190  if( core->columnsavailable )
1191  {
1192  int i;
1193 
1194  for( i = 0; i < core->nvariables; i++ )
1195  if( core->columns[i] )
1196  SCIPfreeBufferArray(scip, &core->columns[i]);
1197 
1198  SCIPfreeBufferArray(scip, &core->columns);
1199  }
1200 
1201  SCIPfreeBufferArrayNull(scip, &core->listcorevariables);
1202  SCIPfreeBufferArray(scip, &core->constraintid);
1203  SCIPfreeBufferArray(scip, &core->nvarconss);
1204  SCIPfreeBufferArray(scip, &core->solgreedy);
1205  SCIPfreeBufferArray(scip, &core->delta_pos);
1206  SCIPfreeBufferArray(scip, &core->delta);
1207  SCIPfreeBufferArray(scip, &core->varincore);
1208  SCIPhashmapFree(&core->mapvariables);
1209 }
1210 
1211 /** computes a new core based on the delta values of variables, see formula (9) in the paper */
1212 static
1213 SCIP_RETCODE redefineCore(
1214  SCIP* scip, /**< master SCIP data structure */
1215  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data; this contains the SCP core */
1216  )
1217 {
1218  SCIP_Bool recomputecolumns = FALSE;
1219  SCIP_Bool recomputerows = FALSE;
1220  SCP_CORE* core;
1221  int* deltaperm;
1222  int nvars;
1223  SCIP_VAR** vars;
1224 
1225  int i;
1226  int j;
1227 
1228  /* assumption: delta values were already computed and are sorted in increasing order, this happens in computeDelta */
1229 
1230  core = &heurdata->core;
1231  vars = heurdata->vars;
1232 
1233  SCIP_CALL( SCIPallocBufferArray(scip, &deltaperm, core->nvariables) );
1234 
1235  /* construct mapping deltaperm: variables -> delta_pos */
1236  for( i = 0; i < core->nvariables; i++ )
1237  {
1238  deltaperm[core->delta_pos[i]] = i;
1239  }
1240 
1241  /* release memory that is used for the core rows and columns. it will be allocated again later on */
1242  removeVarsFromCore(scip, core);
1243  if( core->columnsavailable )
1244  {
1245  recomputecolumns = TRUE;
1246 
1247  for( i = 0; i < core->nvariables; i++ )
1248  {
1249  if( core->columns[i] != NULL )
1250  SCIPfreeBufferArray(scip, &core->columns[i]);
1251  }
1252 
1253  SCIPfreeBufferArray(scip, &core->columns);
1254  core->columns = NULL;
1255  core->columnsavailable = FALSE;
1256  }
1257 
1258  if( core->rowsavailable )
1259  {
1260  recomputerows = TRUE;
1261 
1262  for( i = 0; i < core->nconss; i++ )
1263  {
1264  if( core->rows[i] )
1265  SCIPfreeBufferArray(scip, &core->rows[i]);
1266  }
1267 
1268  SCIPfreeBufferArray(scip, &core->rows);
1269  SCIPfreeBufferArray(scip, &core->nrowvars);
1270 
1271  core->rows = NULL;
1272  core->nrowvars = NULL;
1273  core->rowsavailable = FALSE;
1274  }
1275 
1276  if( core->listcorevariables != NULL )
1277  {
1278  SCIPfreeBufferArray(scip, &core->listcorevariables);
1279  core->listcorevariables = NULL;
1280  }
1281 
1282  /* pick the first 'SCP_CORE_TENT_SIZE'*m columns with lowest delta values to be in the core */
1283  for( i = 0; i < core->nvariables; i++ )
1284  {
1285  if( i >= heurdata->param_core_tent_size * core->nactiveconss )
1286  break;
1287 
1288  addVarToCore(core, core->variables[core->delta_pos[i]]);
1289  }
1290 
1291  /* then add the first 'SCP_CORE_TENT_SIZE' columns covering each row in increasing order of their delta values */
1292  for( i = 0; i < core->nconss; i++ )
1293  {
1294  SCIP_Bool success = FALSE;
1295  int *cols = heurdata->rccols; /*[SCP_CORE_TENT_SIZE];*/
1296  SCIP_Real *coldelta = heurdata->rccoldelta; /*[SCP_CORE_TENT_SIZE];*/
1297 
1298  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1299  if( success == FALSE )
1300  continue;
1301 
1302  for( j = 0; j < heurdata->param_core_tent_size; j++ )
1303  {
1304  cols[j] = 0;
1305  coldelta[j] = SCIP_REAL_MAX;
1306  }
1307 
1308  /* iterate through list of variables that are part of this constraint, and find the 'param_core_tent_size' variables with lowest delta values */
1309  for( j = 0; j < nvars; j++ )
1310  {
1311  int varpos;
1312  int k;
1313  SCIP_Real value;
1314 
1315  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1316 
1317  value = core->delta[deltaperm[varpos]];
1318 
1319  if( j < heurdata->param_core_tent_size )
1320  k = j;
1321  else
1322  k = heurdata->param_core_tent_size - 1;
1323 
1324  if( (j < heurdata->param_core_tent_size) || (coldelta[k] > value) )
1325  {
1326  cols[k] = varpos;
1327  coldelta[k] = value;
1328 
1329  /* insertion sort */
1330  while( (k > 0) && (coldelta[k] < coldelta[k - 1]) )
1331  {
1332  int tmp1 = cols[k - 1];
1333  SCIP_Real tmp2 = coldelta[k - 1];
1334 
1335  cols[k - 1] = cols[k];
1336  coldelta[k - 1] = coldelta[k];
1337  cols[k] = tmp1;
1338  coldelta[k] = tmp2;
1339  k--;
1340  }
1341  }
1342  }
1343 
1344  /* add the variables with lowest delta values to the core */
1345  for( j = 0; j < nvars; j++ )
1346  {
1347  if( j >= heurdata->param_core_tent_size )
1348  break;
1349 
1350  if( !isCoreVariable(core, core->variables[cols[j]]) )
1351  addVarToCore(core, core->variables[cols[j]]);
1352  }
1353  }
1354 
1355  /* construct list of core variables */
1356  SCIP_CALL( SCIPallocBufferArray(scip, &core->listcorevariables, core->ncorevariables) );
1357  j = 0;
1358  for( i = 0; i < core->nvariables; i++ )
1359  {
1360  if( isCoreVariable(core, core->variables[i]) )
1361  core->listcorevariables[j++] = i;
1362  }
1363  assert(j == core->ncorevariables);
1364 
1365  /* recompute core rows and columns if they were available before */
1366  if( recomputecolumns )
1367  SCIP_CALL( computeCoreColumns(scip, core) );
1368 
1369  if( recomputerows )
1370  SCIP_CALL( computeCoreRows(scip, core, heurdata) );
1371 
1372  SCIPfreeBufferArray(scip, &deltaperm);
1373 
1374  SCIPdebugMessage("%d variables are in the refined core\n", core->ncorevariables);
1375  return SCIP_OKAY;
1376 }
1377 
1378 /** for all rows that are covered by the variables in inst->varsfixed, this function adds their indices to inst->rowscovered */
1379 static
1381  SCIP* scip, /**< master SCIP data structure */
1382  SCP_CORE* core, /**< SCP core data structure */
1383  SCP_INSTANCE* inst, /**< initialized SCP instance where some variables may be fixed */
1384  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
1385  )
1386 {
1387  int i;
1388  int j;
1389 
1390  for( i = 0; i < SCIPgetNConss(scip); ++i )
1391  inst->rowscovered[i] = FALSE;
1392  inst->nrowscovered = 0;
1393 
1394  if( core->columnsavailable == TRUE )
1395  {
1396  /* iterate through list of core variables */
1397  for( i = 0; i < core->ncorevariables; i++ )
1398  {
1399  int corevar = core->listcorevariables[i];
1400 
1401  if( !isFixedVariable(inst, core->variables[corevar]) )
1402  continue;
1403 
1404  /* this variable is fixed, so mark its rows as covered */
1405  for( j = 0; j < core->nvarconss[corevar]; j++ )
1406  {
1407  int rowidx = core->columns[corevar][j];
1408 
1409  if( !isRowCovered(inst, rowidx) )
1410  markRowAsCovered(inst, rowidx);
1411  }
1412  }
1413  }
1414  else
1415  {
1416  SCIP_VAR **vars;
1417  int nvars;
1418 
1419  vars = heurdata->vars;
1420 
1421  /* if the core columns are not available, iterate through list of constraints and test which conss are covered */
1422  for( i = 0; i < core->nconss; i++ )
1423  {
1424  SCIP_Bool success = FALSE;
1425 
1426  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1427  if( success == FALSE )
1428  continue;
1429 
1430  for( j = 0; j < nvars; j++ )
1431  {
1432  if( isFixedVariable(inst, vars[j]) )
1433  {
1434  if( !isRowCovered(inst, i) )
1435  markRowAsCovered(inst, i);
1436 
1437  break;
1438  }
1439  }
1440  }
1441  }
1442 
1443  return SCIP_OKAY;
1444 }
1445 
1446 /** checks whether the given 'solution' is actually a set cover */
1447 static
1448 SCIP_RETCODE checkSetCover(
1449  SCIP* scip, /**< master SCIP data structure */
1450  SCP_CORE* core, /**< SCP core data structure */
1451  SCP_INSTANCE* inst, /**< SCP instance */
1452  SCIP_Bool* solution, /**< possible SCP solution */
1453  SCIP_Bool* issetcover, /**< address where TRUE will be stored if 'solution' is a set cover */
1454  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
1455  )
1456 {
1457  SCIP_Bool success;
1458  SCIP_VAR** vars;
1459  int nvars;
1460  int i;
1461  int j;
1462  int varpos;
1463 
1464  assert(scip != NULL);
1465  assert(core != NULL);
1466  assert(inst != NULL);
1467  assert(solution != NULL);
1468  assert(issetcover != NULL);
1469  assert(heurdata != NULL);
1470 
1471  *issetcover = TRUE;
1472  vars = heurdata->vars;
1473 
1474  /* iterate through all constraints and check whether each of them contains a variable that is part of the cover */
1475  for( i = 0; i < core->nconss; i++ )
1476  {
1477  SCIP_Bool rowcovered = FALSE;
1478 
1479  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1480  if( success == FALSE )
1481  continue;
1482 
1483  for( j = 0; j < nvars; j++ )
1484  {
1485  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1486 
1487  if( isVarInSolution(solution, core->variables[varpos]) )
1488  {
1489  rowcovered = TRUE;
1490  break;
1491  }
1492  }
1493 
1494  if( rowcovered == FALSE )
1495  {
1496  SCIPdebugMessage("check set cover: row %i is not covered by any column\n", i);
1497  *issetcover = FALSE;
1498  break;
1499  }
1500  }
1501 
1502  return SCIP_OKAY;
1503 }
1504 
1505 /** computes delta values for variables and fixes new variables within inst, see formula (9) of the paper */
1506 static
1507 SCIP_RETCODE computeDelta(
1508  SCIP* scip, /**< master SCIP data structure */
1509  SCP_CORE* core, /**< SCP core data structure */
1510  SCP_INSTANCE* inst, /**< SCP instance */
1511  SCIP_Real* lagrangiancosts, /**< array with lagrangian costs of the rows */
1512  SCIP_Bool* solution, /**< valid solution to the global SCP instance */
1513  SCIP_Real pi, /**< percentage of rows that need be covered by new fixed variables */
1514  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
1515  )
1516 {
1517  SCIP_Bool success;
1518  SCIP_Real delta_sum;
1519  SCIP_Real delta_max;
1520  int* nvarcovering;
1521  SCIP_VAR** vars;
1522  int nvars;
1523  int i;
1524  int j;
1525  int* delta_pos;
1526  SCIP_Real* delta;
1527 
1528  delta = core->delta;
1529  delta_pos = core->delta_pos;
1530 
1531  SCIP_CALL( SCIPallocBufferArray(scip, &nvarcovering, core->nconss) );
1532  vars = heurdata->vars;
1533 
1534  /* compute nvarcovering[i] = 'number of columns covering row i' */
1535  if( core->rowsavailable )
1536  {
1537  for( i = 0; i < core->nconss; i++ )
1538  {
1539  nvarcovering[i] = 0;
1540 
1541  if( !SCIPconsIsActive(core->conss[i]) )
1542  continue;
1543 
1544  for( j = 0; j < core->nrowvars[i]; j++ )
1545  {
1546  int varpos = core->rows[i][j];
1547  if( isVarInSolution(solution, core->variables[varpos]) )
1548  nvarcovering[i]++;
1549  }
1550  }
1551  }
1552  else
1553  {
1554  for( i = 0; i < core->nconss; i++ )
1555  {
1556  nvarcovering[i] = 0;
1557 
1558  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1559  if( success == FALSE )
1560  continue;
1561 
1562  for( j = 0; j < nvars; j++ )
1563  {
1564  int varpos;
1565 
1566  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
1567 
1568  if( isVarInSolution(solution, core->variables[varpos]) )
1569  nvarcovering[i]++;
1570  }
1571  }
1572  }
1573 
1574  /* compute delta values of variables that are part of 'solution', see formula (9) */
1575  for( i = 0; i < core->nvariables; i++ )
1576  {
1577  delta[i] = SCIP_REAL_MAX;
1578  delta_pos[i] = i;
1579 
1580  if( !isVarInSolution(solution, core->variables[i]) )
1581  continue;
1582 
1583  delta[i] = lagrangiancosts[i];
1584  if( delta[i] < 0.0 )
1585  delta[i] = 0.0;
1586 
1587  for( j = 0; j < core->nvarconss[i]; j++ )
1588  {
1589  int rowpos = core->columns[i][j];
1590 
1591  if( isRowCovered(inst, rowpos) )
1592  continue;
1593 
1594  delta[i] += lagrangiancosts[rowpos] * (nvarcovering[rowpos] - 1) / ((SCIP_Real) nvarcovering[rowpos]);
1595  }
1596  }
1597 
1598  /* sort the variables in increasing order of their delta values, delta_pos is used to track the variables within the array */
1599  SCIPsortRealInt(delta, delta_pos, core->nvariables);
1600 
1601  /* fix variables as long as the sum of their delta values is not greater than delta_max.
1602  * We assume that variables with small delta values are likely part of an optimal solution. */
1603  delta_sum = 0.0;
1604  delta_max = core->nactiveconss * pi;
1605 
1606  for( i = 0; i < SCIPgetNVars(scip); ++i )
1607  inst->varsfixed[i] = FALSE;
1608  inst->nvarsfixed = 0;
1609 
1610  inst->costsfixed = 0.0;
1611 
1612  for( i = 0; i < core->nvariables; i++ )
1613  {
1614  int varpos = delta_pos[i];
1615 
1616  if( !isVarInSolution(solution, core->variables[varpos]) )
1617  break;
1618 
1619  inst->costsfixed += SCIPvarGetObj(core->variables[varpos]);
1620  delta_sum += delta[i];
1621 
1622  /* fix variable delta_pos[i] */
1623  fixVariable(inst, core->variables[varpos]);
1624 
1625  if( delta_sum >= delta_max )
1626  break;
1627  }
1628 
1629  SCIPfreeBufferArray(scip, &nvarcovering);
1630  return SCIP_OKAY;
1631 }
1632 
1633 /* removes redundant columns from 'solution' by removing them one by one (in the order of decreasing costs),
1634  * and checking whether all rows are still covered */
1635 static
1637  SCIP* scip, /**< master SCIP data structure */
1638  SCIP_HEURDATA* heurdata, /**< pointer to the heuristic's data */
1639  SCIP_Bool* solution, /**< solution to the global SCP instance */
1640  SCIP_Real* solcosts /**< pointer to original costs of 'solution'; contains updated costs */
1641  )
1642 {
1643  SCP_CORE* core;
1644  SCIP_Real* costs;
1645  int* varpos;
1646  int* nvarcovering;
1647  SCIP_VAR** vars;
1648  int nvars;
1649  SCIP_Bool success;
1650  int solsize;
1651 
1652  int i;
1653  int j;
1654 
1655  core = &heurdata->core;
1656 
1657  if( core->columnsavailable == FALSE )
1658  {
1659  SCIPdebugMessage("can only remove redundant columns if all core columns are available\n");
1660  return SCIP_OKAY;
1661  }
1662 
1663  SCIP_CALL( SCIPallocBufferArray(scip, &costs, core->nvariables) );
1664  SCIP_CALL( SCIPallocBufferArray(scip, &varpos, core->nvariables) );
1665  SCIP_CALL( SCIPallocBufferArray(scip, &nvarcovering, core->nconss) );
1666 
1667  vars = heurdata->vars;
1668 
1669  /* compute nvarcovering[i] = 'number of columns covering row i' */
1670  if( core->rowsavailable )
1671  {
1672  for( i = 0; i < core->nconss; i++ )
1673  {
1674  nvarcovering[i] = 0;
1675 
1676  if( SCIPconsIsActive(core->conss[i]) == FALSE )
1677  continue;
1678 
1679  if( core->nrowvars[i] == 0 )
1680  continue;
1681 
1682  for( j = 0; j < core->nrowvars[i]; j++ )
1683  {
1684  int vpos = core->rows[i][j];
1685 
1686  if( isVarInSolution(solution, core->variables[vpos]) )
1687  nvarcovering[i]++;
1688  }
1689  }
1690  }
1691  else
1692  {
1693  /* iterate through the constraints if the core rows are not available */
1694  for( i = 0; i < core->nconss; i++ )
1695  {
1696  nvarcovering[i] = 0;
1697 
1698  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1699  if( success == FALSE )
1700  continue;
1701 
1702  for( j = 0; j < nvars; j++ )
1703  {
1704  int vpos;
1705 
1706  SCIP_CALL( getVarIndex(scip, core, vars[j], &vpos) );
1707 
1708  if( isVarInSolution(solution, core->variables[vpos]) )
1709  nvarcovering[i]++;
1710  }
1711  }
1712  }
1713 
1714  /* construct array of costs of variables that are part of the solution.
1715  * If v_1, ..., v_solsize are the variables of the solution, then
1716  * costs[i] = -costs(v_i), and varpos[i] = i.
1717  * We use negative costs because we only sort them in increasing order and want to remove the ones with highest costs first.
1718  * The 'varpos' array is used to track which variable is where within the sorted array (it is permutated in the same way as 'costs') */
1719  solsize = 0;
1720  for( i = 0; i < core->nvariables; i++ )
1721  {
1722  if( isVarInSolution(solution, core->variables[i]) )
1723  {
1724  costs[solsize] = -SCIPvarGetObj(core->variables[i]);
1725  varpos[solsize] = i;
1726  solsize++;
1727  }
1728  }
1729 
1730  /* sort variables by decreasing costs */
1731  SCIPsortRealInt(costs, varpos, solsize);
1732 
1733  /* for each variable of the solution: try to remove it and check whether it is still a set cover */
1734  for( i = 0; i < solsize; i++ )
1735  {
1736  SCIP_Bool canremove = TRUE;
1737  int vpos = varpos[i];
1738 
1739  /* we can only remove this variable if every column it covers is covered by another variable */
1740  for( j = 0; j < core->nvarconss[vpos]; j++ )
1741  {
1742  if( nvarcovering[core->columns[vpos][j]] == 1 )
1743  {
1744  canremove = FALSE;
1745  break;
1746  }
1747  }
1748 
1749  if( canremove == TRUE )
1750  {
1751  removeVarFromSolution(solution, core->variables[vpos]);
1752  *solcosts -= SCIPvarGetObj(core->variables[vpos]);
1753 
1754  for( j = 0; j < core->nvarconss[vpos]; j++ )
1755  {
1756  nvarcovering[core->columns[vpos][j]]--;
1757  }
1758  }
1759  }
1760 
1761  SCIPfreeBufferArray(scip, &nvarcovering);
1762  SCIPfreeBufferArray(scip, &varpos);
1763  SCIPfreeBufferArray(scip, &costs);
1764 
1765  return SCIP_OKAY;
1766 }
1767 
1768 /** greedy set cover algorithm that uses lagrangian costs instead of original costs. */
1769 static
1770 SCIP_RETCODE greedySetCover(
1771  SCIP* scip, /**< master SCIP data structure */
1772  SCP_CORE* core, /**< SCP core data structure */
1773  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
1774  SCP_Lagrange_Sol* mult, /**< lagrange multiplier that is used to compute lagrangian costs */
1775  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
1776  )
1777 {
1778  SCP_INSTANCE* greedyinst;
1779  int nrowsuncovered = 0;
1780  SCIP_Bool success = FALSE;
1781  int nvars;
1782  PQUEUE* prioqueue;
1783  int* colpos;
1784  int* colmu;
1785  SCIP_Real* colgamma;
1786  SCIP_Real* colscore;
1787  SCIP_VAR** vars;
1788 
1789  int i;
1790  int j;
1791  int k;
1792 
1793  core->nsolgreedy = 0;
1794  mult->ubgreedylocal = 0.0;
1795  clearSolution(scip, mult->xgreedylocal);
1796 
1797  greedyinst = &heurdata->greedyinst;
1798  for( i = 0; i < SCIPgetNConss(scip); ++i )
1799  greedyinst->rowscovered[i] = FALSE;
1800  greedyinst->nrowscovered = 0;
1801 
1802  /* count number of uncovered rows and add covered rows to 'greedyinst' */
1803  for( i = 0; i < core->nconss; i++ )
1804  {
1805  if( SCIPconsIsActive(core->conss[i]) == FALSE )
1806  continue;
1807 
1808  /* this is actually necessary because there exist constraints where this fails, and we simply need to ignore them */
1809  SCIP_CALL( SCIPgetConsNVars(scip, core->conss[i], &nvars, &success) );
1810  if( success == FALSE )
1811  continue;
1812 
1813  if( isRowCovered(inst, i) )
1814  markRowAsCovered(greedyinst, i);
1815  else
1816  nrowsuncovered++;
1817  }
1818 
1819  if( core->columnsavailable == FALSE )
1820  {
1821  SCIPerrorMessage("greedy algorithm requires core columns to be available\n");
1822  SCIPABORT();
1823  }
1824 
1825  /* compute scores of all (unfixed) variables and add them to the priority queue */
1826  colpos = heurdata->greedycolpos; /* this array contains the position of the element within the queue and is needed to update the costs */
1827  colmu = heurdata->greedycolmu; /* for each column this contains the number of uncovered rows that it covers */
1828  colgamma = heurdata->greedycolgamma; /* for each column this contains the lagrangian costs of uncovered rows that it covers */
1829  colscore = heurdata->greedycolscore; /* for each column this contains its score */
1830  vars = heurdata->vars;
1831 
1832  prioqueue = &heurdata->greedyqueue;
1833  pqueue_clear(scip, prioqueue);
1834 
1835  /* compute initial scores for all unfixed columns, see "3. Heuristic Phase" of the paper, and add them to the priority queue */
1836  for( i = 0; i < core->nvariables; i++ )
1837  {
1838  int varpos = i;
1839  int mu = 0;
1840  SCIP_Real gamma = SCIPvarGetObj(core->variables[varpos]);
1841 
1842  colmu[i] = 0;
1843  colgamma[i] = 0;
1844  colpos[i] = 0;
1845  colscore[i] = 0.0;
1846 
1847  if( !isCoreVariable(core, core->variables[varpos]) )
1848  continue;
1849  else if( isFixedVariable(inst, core->variables[varpos]) )
1850  continue;
1851 
1852  for( j = 0; j < core->nvarconss[varpos]; j++ )
1853  {
1854  if( !isRowCovered(greedyinst, core->columns[varpos][j]) )
1855  {
1856  gamma -= mult->u[core->columns[varpos][j]];
1857  mu++;
1858  }
1859  }
1860 
1861  /* skip columns that do not cover anything */
1862  if( mu > 0 )
1863  {
1864  SCIP_Real score = (gamma > 0) ? (gamma / ((SCIP_Real) mu)) : (gamma * mu);
1865  colmu[i] = mu;
1866  colgamma[i] = gamma;
1867  colscore[i] = score;
1868 
1869  SCIP_CALL( pqueue_insert(scip, prioqueue, colscore[i], i, &(colpos[i])) );
1870  }
1871  }
1872 
1873  /* main loop: add variables to the solution as long as there are uncovered rows */
1874  while( nrowsuncovered > 0 )
1875  {
1876  int mincolumn;
1877 
1878  /* get variable with minimum score */
1879  pqueue_get_min(scip, prioqueue, &mincolumn);
1880 
1881  /* add variable 'variables[mincolumn]' to the set cover */
1882  addVarToSolution(mult->xgreedylocal, core->variables[mincolumn]);
1883  mult->ubgreedylocal += SCIPvarGetObj(core->variables[mincolumn]);
1884  core->solgreedy[core->nsolgreedy++] = mincolumn;
1885 
1886  /* mark this variable, so its score won't be updated in the future */
1887  colmu[mincolumn] = 0;
1888 
1889  /* mark all rows covered by this variable as covered and update scores of variables that are affected by this */
1890  for( j = 0; j < core->nvarconss[mincolumn]; j++ )
1891  {
1892  int columnpos = core->columns[mincolumn][j];
1893 
1894  if( !isRowCovered(greedyinst, columnpos) )
1895  {
1896  markRowAsCovered(greedyinst, columnpos);
1897  nrowsuncovered--;
1898 
1899  /* update scores of columns covering this row */
1900  SCIP_CALL( getConsVars(scip, core, columnpos, vars, &nvars, &success) );
1901  if( success == FALSE )
1902  continue;
1903 
1904  /* for each variable: subtract u[i] from the variable's costs */
1905  for( k = 0; k < nvars; k++ )
1906  {
1907  SCIP_Real score;
1908  int varpos;
1909 
1910  SCIP_CALL( getVarIndex(scip, core, vars[k], &varpos) );
1911 
1912  /* skip non-core variables */
1913  if( !isCoreVariable(core, core->variables[varpos]) )
1914  continue;
1915 
1916  if( colmu[varpos] == 0 )
1917  continue;
1918 
1919  score = colscore[varpos];
1920 
1921  /* compute a new score for this variable */
1922  colscore[varpos] = SCIP_REAL_MAX;
1923  colmu[varpos]--;
1924  colgamma[varpos] += mult->u[columnpos];
1925 
1926  if( colmu[varpos] > 0 )
1927  colscore[varpos] = (colgamma[varpos] > 0) ? (colgamma[varpos] / ((SCIP_Real) colmu[varpos])) : (colgamma[varpos] * colmu[varpos]);
1928 
1929  /* either decrease or increase the key, depending on whether it decreased or increased */
1930  if( score > colscore[varpos] )
1931  pqueue_decrease_key(scip, prioqueue, colpos[varpos], colscore[varpos]);
1932  else
1933  pqueue_increase_key(scip, prioqueue, colpos[varpos], colscore[varpos]);
1934  }
1935  }
1936  }
1937  }
1938 
1939  return SCIP_OKAY;
1940 }
1941 
1942 /** computes lagrangian costs for all columns, only considering rows that are uncovered by fixed variables in 'inst' */
1943 static
1945  SCIP* scip, /**< master SCIP data structure */
1946  SCP_CORE* core, /**< SCP core data structure */
1947  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
1948  SCP_Lagrange_Sol* mult, /**< lagrange multiplier that is used to compute the lagrangian costs*/
1949  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
1950  )
1951 {
1952  SCIP_VAR** vars;
1953  int nvars;
1954 
1955  int i;
1956  int j;
1957 
1958  /* set all lagrangian costs to objective values */
1959  for( i = 0; i < core->nvariables; i++ )
1960  mult->lagrangiancostslocal[i] = SCIPvarGetObj(core->variables[i]);
1961 
1962  vars = heurdata->vars;
1963 
1964  if( core->rowsavailable )
1965  {
1966  for( i = 0; i < core->nconss; i++ )
1967  {
1968  if( SCIPconsIsActive(core->conss[i]) == FALSE )
1969  continue;
1970 
1971  if( isRowCovered(inst, i) )
1972  continue;
1973 
1974  for( j = 0; j < core->nrowvars[i]; j++ )
1975  {
1976  int varpos = core->rows[i][j];
1977  mult->lagrangiancostslocal[varpos] -= mult->u[i];
1978  }
1979  }
1980  }
1981  else
1982  {
1983  for( i = 0; i < core->nconss; i++ )
1984  {
1985  SCIP_Bool success;
1986 
1987  /* skip rows that are not part of the reduced instance */
1988  if( isRowCovered(inst, i) )
1989  continue;
1990 
1991  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
1992  if( success == FALSE )
1993  continue;
1994 
1995  /* for each core variable: subtract u[i] from the variable's costs */
1996  for( j = 0; j < nvars; j++ )
1997  {
1998  int varpos;
1999 
2000  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2001 
2002  mult->lagrangiancostslocal[varpos] -= mult->u[i];
2003  }
2004  }
2005  }
2006 
2007  return SCIP_OKAY;
2008 }
2009 
2010 /** computes lagrangian costs for all columns of the unrestricted instance and a global lower bound. */
2011 static
2013  SCIP* scip, /**< master SCIP data structure */
2014  SCP_CORE* core, /**< SCP core data structure */
2015  SCP_Lagrange_Sol* mult, /**< lagrange multiplier that is used to compute the lagrangian costs*/
2016  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
2017  )
2018 {
2019  SCIP_VAR** vars;
2020  int nvars;
2021  int i;
2022  int j;
2023 
2024  /* set all lagrangian costs to objective values */
2025  for( i = 0; i < core->nvariables; i++ )
2026  mult->lagrangiancostsglobal[i] = SCIPvarGetObj(core->variables[i]);
2027 
2028  vars = heurdata->vars;
2029 
2030  for( i = 0; i < core->nconss; i++ )
2031  {
2032  SCIP_Bool success = FALSE;
2033 
2034  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
2035  if( success == FALSE )
2036  continue;
2037 
2038  /* for each core variable: subtract u[i] from the variable's costs */
2039  for( j = 0; j < nvars; j++ )
2040  {
2041  int varpos;
2042 
2043  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2044 
2045  mult->lagrangiancostsglobal[varpos] -= mult->u[i];
2046  }
2047 
2048  mult->lblagrangeglobal += mult->u[i];
2049  }
2050 
2051  /* compute global lower bound */
2052  for( i = 0; i < core->nvariables; i++ )
2053  {
2054  if( mult->lagrangiancostsglobal[i] < 0.0 )
2055  mult->lblagrangeglobal += mult->lagrangiancostsglobal[i];
2056  }
2057 
2058  return SCIP_OKAY;
2059 }
2060 
2061 /** computes an optimal solution to the lagrangian relaxation, see formulae (4), (5) in the paper */
2062 static
2064  SCIP* scip, /**< master SCIP data structure */
2065  SCP_CORE* core, /**< SCP core data structure */
2066  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
2067  SCP_Lagrange_Sol* mult, /**< lagrange multiplier */
2068  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
2069  )
2070 {
2071  SCIP_VAR** vars;
2072  SCIP_Bool success;
2073  int nvars;
2074 
2075  int i;
2076  int j;
2077 
2078  mult->lblagrangelocal = 0.0;
2079  mult->lblagrangeglobal = 0.0;
2080 
2081  /* compute lagrangian costs for the reduced instance */
2082  SCIP_CALL( computeLocalLagrangianCosts(scip, core, inst, mult, heurdata) );
2083 
2084  /* also, compute costs for the unrestricted instance. This even computes mult->lblagrangeglobal */
2085  SCIP_CALL( computeGlobalLagrangianCosts(scip, core, mult, heurdata) );
2086 
2087  /* compute an optimal solution for 'inst', but only save its costs in lblagrangelocal */
2088  for( i = 0; i < core->ncorevariables; i++ )
2089  {
2090  int varpos = core->listcorevariables[i];
2091 
2092  if( mult->lagrangiancostslocal[varpos] < 0.0 )
2093  {
2094  if( !isFixedVariable(inst, core->variables[varpos]) )
2095  mult->lblagrangelocal = mult->lblagrangelocal + mult->lagrangiancostslocal[varpos];
2096  }
2097  }
2098 
2099  /* the subgradient can be computed faster when the core rows are available */
2100  if( core->rowsavailable )
2101  {
2102  /* compute the subgradient */
2103  for( i = 0; i < core->nconss; i++ )
2104  {
2105  mult->subgradient[i] = 0.0;
2106 
2107  if( SCIPconsIsActive(core->conss[i]) == FALSE )
2108  continue;
2109 
2110  if( isRowCovered(inst, i) )
2111  continue;
2112 
2113  if( core->nrowvars[i] == 0 )
2114  continue;
2115 
2116  mult->subgradient[i] = 1.0;
2117 
2118  for( j = 0; j < core->nrowvars[i]; j++ )
2119  {
2120  int varpos = core->rows[i][j];
2121 
2122  if( mult->lagrangiancostslocal[varpos] < 0.0 )
2123  mult->subgradient[i] -= 1.0;
2124  }
2125 
2126  /* to get a lower bound from the multiplier we still need to subtract sum_i u[i] from lblagrangelocal */
2127  mult->lblagrangelocal = mult->lblagrangelocal + mult->u[i];
2128  }
2129  }
2130  else
2131  {
2132  /* the core rows are NOT available, so we need to iterate through the constraints */
2133  vars = heurdata->vars;
2134 
2135  for( i = 0; i < core->nconss; i++ )
2136  {
2137  mult->subgradient[i] = 0.0;
2138 
2139  if( isRowCovered(inst, i) )
2140  continue;
2141 
2142  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
2143  if( success == FALSE )
2144  continue;
2145 
2146  mult->subgradient[i] = 1.0;
2147 
2148  for( j = 0; j < nvars; j++ )
2149  {
2150  int varpos;
2151 
2152  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2153 
2154  if( !isCoreVariable(core, core->variables[varpos]) )
2155  continue;
2156 
2157  if( mult->lagrangiancostslocal[varpos] < 0.0 )
2158  mult->subgradient[i] -= 1.0;
2159  }
2160 
2161  mult->lblagrangelocal = mult->lblagrangelocal + mult->u[i];
2162  }
2163  }
2164 
2165  return SCIP_OKAY;
2166 }
2167 
2168 /** subgradient method with Held-Karp update of subgradients */
2169 static
2171  SCIP* scip, /**< master SCIP data structure */
2172  SCP_CORE* core, /**< SCP core data structure */
2173  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
2174  SCP_Lagrange_Sol* best_mult_lb, /**< lagrange multiplier that gives the best lower bound for 'inst' */
2175  SCIP_Real bestub, /**< costs of the best known solution for 'inst' */
2176  SCIP_HEURDATA* heurdata /**< pointer the heuristic's data */
2177  )
2178 {
2179  int max_iter;
2180  SCP_Lagrange_Sol last_mult;
2181  SCP_Lagrange_Sol next_mult;
2182  SCP_Lagrange_Sol tmp_mult;
2183  SCIP_Real norm;
2184  SCIP_Real lambda = heurdata->param_lambda;
2185  SCIP_Bool lambda_adjustments = heurdata->param_lambda_adjustments;
2186  SCIP_Real* last_lb = heurdata->sglastlb;
2187  SCIP_Real stop_crit_lb = 0.0;
2188  int iter;
2189  int last_data_pos = 0;
2190  int i;
2191 
2192  SCIP_CALL( allocateMemoryForSolution(scip, core, &next_mult) );
2193  SCIP_CALL( allocateMemoryForSolution(scip, core, &last_mult) );
2194 
2195  /* save data from best lower bound multiplier in last_mult */
2196  copySolution(scip, core, &last_mult, best_mult_lb);
2197 
2198  /* permutate best u by multiplying each entry with a uniformly random value in the range [0.9, 1.1] */
2199  for( i = 0; i < core->nconss; i++ )
2200  {
2201  if( !isRowCovered(inst, i) )
2202  last_mult.u[i] = SCIPrandomGetReal(heurdata->randnumgen, 0.9, 1.1) * last_mult.u[i];
2203  else
2204  last_mult.u[i] = 0.0;
2205  }
2206 
2207  /* bestub is an upper bound for the restricted instance, without the fixed costs */
2208  bestub -= inst->costsfixed;
2209 
2210  /* maximum number of subgradient iterations */
2211  max_iter = 10 * core->nconss;
2212 
2213  for( iter = 0; iter < max_iter; iter++ )
2214  {
2215  /* compute norm of the subgradient */
2216  norm = 0.0;
2217  for( i = 0; i < core->nconss; i++ )
2218  norm = norm + last_mult.subgradient[i] * last_mult.subgradient[i]; /* we have subgradient[i] = 0.0 if row i is not to be considered */
2219 
2220  /* the norm is zero if an optimal solution (to the restricted instance) was found */
2221  if( SCIPisZero(scip, norm) )
2222  break;
2223 
2224  /* Held-Karp update */
2225  for( i = 0; i < core->nconss; i++ )
2226  {
2227  SCIP_Real hk = last_mult.u[i] + lambda * (bestub - last_mult.lblagrangelocal) * last_mult.subgradient[i] / norm; /*lint !e795*/
2228  next_mult.u[i] = 0.0;
2229 
2230  if( SCIPisPositive(scip, 0.0) )
2231  next_mult.u[i] = hk;
2232  }
2233 
2234  /* for each multiplier compute an optimal solution to get a lower bound */
2235  SCIP_CALL( computeOptimalSolution(scip, core, inst, &next_mult, heurdata) );
2236 
2237  /* save the best incumbent */
2238  if( SCIPisGT(scip, next_mult.lblagrangelocal, best_mult_lb->lblagrangelocal) )
2239  copySolution(scip, core, best_mult_lb, &next_mult);
2240 
2241  if( SCIPisGT(scip, next_mult.lblagrangeglobal, heurdata->multbestlbtotal.lblagrangeglobal) )
2242  {
2243  copySolution(scip, core, &heurdata->multbestlbtotal, &next_mult);
2244  }
2245 
2246  /* adjust lambda, the step size, if necessary */
2247  if( lambda_adjustments )
2248  {
2249  /* save last 'p' lower and upper bounds */
2250  last_lb[last_data_pos++] = next_mult.lblagrangelocal;
2251 
2252  /* lambda is adjusted after param_lambda_p iterations */
2253  if( last_data_pos >= heurdata->param_lambda_p )
2254  {
2255  SCIP_Real max_lb = last_lb[0];
2256  SCIP_Real min_lb = last_lb[0];
2257 
2258  /* find best and worst bounds that occurred within the last lambda_p iterations */
2259  for( i = 1; i < heurdata->param_lambda_p; i++ )
2260  {
2261  if( SCIPisGT(scip, last_lb[i], max_lb) )
2262  max_lb = last_lb[i];
2263  if( SCIPisLT(scip, last_lb[i], min_lb) )
2264  min_lb = last_lb[i];
2265  }
2266 
2267  /* if min_lb and max_lb differ by more than 1%, lambda is halved */
2268  if( SCIPisGT(scip, max_lb - min_lb, 0.01) )
2269  lambda = lambda / 2;
2270 
2271  /* if min_lb and max_lb differ by less than 0.1%, lambda is multiplied by 1.5 */
2272  if( SCIPisLT(scip, max_lb - min_lb, 0.001) )
2273  lambda = lambda * 1.5;
2274  last_data_pos = 0;
2275  }
2276  }
2277 
2278  /* swap next_mult and last_mult */
2279  tmp_mult = last_mult;
2280  last_mult = next_mult;
2281  next_mult = tmp_mult;
2282 
2283  /* check stopping criteria after stop_crit_iter iterations */
2284  if( iter % heurdata->param_stop_crit_iter == 0 )
2285  {
2286  /* we test how the current best lower bound compares to the best lower bound that was known 'stop_crit_iter' iterations ago:
2287  * if the absolute difference is smaller than param_stop_crit_diff and
2288  * the ratio lb_now / lb_old is at least stop_crit_ration
2289  * then we stop */
2290  if( (iter > 0) && (best_mult_lb->lblagrangelocal - stop_crit_lb <= heurdata->param_stop_crit_diff )
2291  && (stop_crit_lb / best_mult_lb->lblagrangelocal >= heurdata->param_stop_crit_ratio) )
2292  {
2293  break;
2294  }
2295 
2296  stop_crit_lb = best_mult_lb->lblagrangelocal;
2297  }
2298  }
2299 
2300  freeMemoryForSolution(scip, &last_mult);
2301  freeMemoryForSolution(scip, &next_mult);
2302 
2303  return SCIP_OKAY;
2304 }
2305 
2306 /** computes an initial lagrange multiplier, see formula (8) in the paper */
2307 static
2309  SCIP* scip, /**< master SCIP data structure */
2310  SCP_CORE* core, /**< SCP core data structure */
2311  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
2312  SCP_Lagrange_Sol* mult, /**< initialized lagrange multiplier were the new one will be stored */
2313  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
2314  )
2315 {
2316  int i;
2317  int j;
2318 
2319  if( core->columnsavailable == TRUE )
2320  {
2321  for( i = 0; i < core->nconss; i++ )
2322  mult->u[i] = SCIP_REAL_MAX;
2323 
2324  for( i = 0; i < core->nvariables; i++ )
2325  {
2326  int nuncovered = 0;
2327 
2328  if( !isCoreVariable(core, core->variables[i]) )
2329  continue;
2330 
2331  if( isFixedVariable(inst, core->variables[i]) )
2332  continue;
2333 
2334  /* count how many uncovered, active rows this column covers */
2335  for( j = 0; j < core->nvarconss[i]; j++ )
2336  {
2337  int colpos = core->columns[i][j];
2338 
2339  /* skip inactive constraints */
2340  if( SCIPconsIsActive(core->conss[colpos]) == FALSE )
2341  continue;
2342 
2343  /* skip covered rows */
2344  if( isRowCovered(inst, colpos) )
2345  continue;
2346 
2347  nuncovered++;
2348  }
2349 
2350  /* if this column covers any rows, update their cost if necessary */
2351  if( nuncovered > 0 )
2352  {
2353  SCIP_Real costs = SCIPvarGetObj(core->variables[i]) / ((SCIP_Real) nuncovered);
2354  int nvars;
2355 
2356  for( j = 0; j < core->nvarconss[i]; j++ )
2357  {
2358  int colpos = core->columns[i][j];
2359  SCIP_Bool success = FALSE;
2360 
2361  if( isRowCovered(inst, colpos) )
2362  continue;
2363 
2364  if( SCIPconsIsActive(core->conss[colpos]) == FALSE )
2365  continue;
2366 
2367  SCIP_CALL( SCIPgetConsNVars(scip, core->conss[colpos], &nvars, &success) );
2368  if( success == FALSE )
2369  continue;
2370 
2371  if( costs < mult->u[colpos] )
2372  mult->u[colpos] = costs;
2373  }
2374  }
2375  }
2376  }
2377  else
2378  {
2379  SCIP_VAR **vars;
2380  int nvars;
2381  int *nuncoveredactive;
2382 
2383  vars = heurdata->vars;
2384  SCIP_CALL( SCIPallocBufferArray(scip, &nuncoveredactive, core->nvariables) );
2385 
2386  for( i = 0; i < core->nvariables; i++ )
2387  nuncoveredactive[i] = 0;
2388 
2389  for( i = 0; i < core->nconss; i++ )
2390  {
2391  SCIP_Bool success = FALSE;
2392 
2393  if( isRowCovered(inst, i) )
2394  continue;
2395 
2396  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
2397  if( success == FALSE )
2398  continue;
2399 
2400  for( j = 0; j < nvars; j++ )
2401  {
2402  int varpos;
2403 
2404  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2405 
2406  if( !isCoreVariable(core, vars[j]) )
2407  continue;
2408  else
2409  nuncoveredactive[varpos]++;
2410  }
2411  }
2412 
2413  for( i = 0; i < core->nconss; i++ )
2414  {
2415  SCIP_Bool found = FALSE;
2416  SCIP_Bool success = FALSE;
2417 
2418  if( isRowCovered(inst, i) )
2419  {
2420  mult->u[i] = 0.0;
2421  continue;
2422  }
2423 
2424  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
2425  if( success == FALSE )
2426  continue;
2427 
2428  for( j = 0; j < nvars; j++ )
2429  {
2430  if( !isCoreVariable(core, vars[j]) )
2431  continue;
2432  else
2433  {
2434  int varpos;
2435  SCIP_Real costs = SCIP_REAL_MAX;
2436 
2437  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2438 
2439  if( nuncoveredactive[varpos] > 0 )
2440  costs = SCIPvarGetObj(core->variables[varpos]) / ((SCIP_Real) nuncoveredactive[varpos]);
2441 
2442  if( (found == FALSE) || (costs < mult->u[i]) )
2443  {
2444  found = TRUE;
2445  mult->u[i] = costs;
2446  }
2447  }
2448  }
2449  }
2450 
2451  SCIPfreeBufferArray(scip, &nuncoveredactive);
2452  }
2453 
2454  return SCIP_OKAY;
2455 }
2456 
2457 /* Explores a neighborhood of lagrange multipliers around 'startmult',
2458  * and calls the greedy algorithm for each multiplier. The neighborhood is
2459  * simply computed by doing the Held-Karp subgradient updates. The best solutions
2460  * are extended to solutions of the unrestricted instance.
2461  */
2462 static
2463 SCIP_RETCODE exploreNeighborhood(
2464  SCIP* scip, /**< master SCIP data structure */
2465  SCP_Lagrange_Sol* startmult, /**< lagrange multiplier where the search is started */
2466  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
2467  )
2468 {
2469  SCP_CORE* core;
2470  SCP_INSTANCE* subinst;
2471  SCP_Lagrange_Sol mult;
2472  SCIP_VAR** vars;
2473  SCIP_Bool success;
2474  SCIP_Real bestub;
2475  SCIP_Real norm;
2476  SCIP_Real costs;
2477  int nvars;
2478  int iter;
2479  int i;
2480  int j;
2481 
2482  core = &heurdata->core;
2483  subinst = &heurdata->subinst;
2484  vars = heurdata->vars;
2485  bestub = heurdata->bestub - subinst->costsfixed;
2486 
2487  SCIP_CALL( allocateMemoryForSolution(scip, &heurdata->core, &mult) );
2488  copySolution(scip, core, &mult, startmult);
2489 
2490  /* perform 'param_greedy_max_iter' subgradient iterations */
2491  for( iter = 0; iter < heurdata->param_greedy_max_iter; iter++ )
2492  {
2493  /* compute subgradient for 'mult' */
2494  if( core->rowsavailable )
2495  {
2496  for( i = 0; i < core->nconss; i++ )
2497  {
2498  mult.subgradient[i] = 0.0;
2499 
2500  if( SCIPconsIsActive(core->conss[i]) == FALSE )
2501  continue;
2502 
2503  if( isRowCovered(subinst, i) )
2504  continue;
2505 
2506  if( core->nrowvars[i] == 0 )
2507  continue;
2508 
2509  mult.subgradient[i] = 1.0;
2510 
2511  for( j = 0; j < core->nrowvars[i]; j++ )
2512  {
2513  int varpos = core->rows[i][j];
2514 
2515  if( mult.lagrangiancostslocal[varpos] < 0.0 )
2516  mult.subgradient[i] -= 1.0;
2517  }
2518  }
2519  }
2520  else
2521  {
2522  for( i = 0; i < core->nconss; i++ )
2523  {
2524  mult.subgradient[i] = 0.0;
2525 
2526  if( isRowCovered(subinst, i) )
2527  continue;
2528 
2529  SCIP_CALL( getConsVars(scip, core, i, vars, &nvars, &success) );
2530  if( success == FALSE )
2531  continue;
2532 
2533  mult.subgradient[i] = 1.0;
2534 
2535  for( j = 0; j < nvars; j++ )
2536  {
2537  int varpos;
2538 
2539  SCIP_CALL( getVarIndex(scip, core, vars[j], &varpos) );
2540 
2541  if( !isCoreVariable(core, core->variables[varpos]) )
2542  continue;
2543 
2544  if( mult.lagrangiancostslocal[varpos] < 0.0 )
2545  mult.subgradient[i] -= 1.0;
2546  }
2547  }
2548  }
2549 
2550  /* compute norm of subgradient */
2551  norm = 0.0;
2552  for( i = 0; i < core->nconss; i++ )
2553  norm = norm + mult.subgradient[i] * mult.subgradient[i];
2554 
2555  /* the norm is zero if an optimal solution (to the restricted instance) was found */
2556  if( SCIPisZero(scip, norm) )
2557  break;
2558 
2559  /* Held-Karp update */
2560  for( i = 0; i < core->nconss; i++ )
2561  {
2562  SCIP_Real hk = mult.u[i] + heurdata->param_lambda * (bestub - mult.lblagrangelocal) * mult.subgradient[i] / norm; /*lint !e795*/
2563  mult.u[i] = 0.0;
2564 
2565  if( hk > 0.0 )
2566  mult.u[i] = hk;
2567  }
2568 
2569  /* first, compute lagrangian costs */
2570  SCIP_CALL( computeOptimalSolution(scip, core, subinst, &mult, heurdata) );
2571 
2572  /* save this multiplier if it gives a better local lower bound than the one known */
2573  if( mult.lblagrangelocal > heurdata->multbestlbsubinst.lblagrangelocal )
2574  copySolution(scip, core, &heurdata->multbestlbsubinst, &mult);
2575 
2576  /* save this multiplier if it gives a better global lower bound */
2577  if( mult.lblagrangeglobal > heurdata->multbestlbtotal.lblagrangeglobal )
2578  copySolution(scip, core, &heurdata->multbestlbtotal, &mult);
2579 
2580  /* apply the greedy algorithm and extend it to a solution of the unrestricted instance by adding fixed variables */
2581  SCIP_CALL( greedySetCover(scip, core, subinst, &mult, heurdata) );
2582  extendSolution(scip, core, subinst, mult.xgreedylocal);
2583 
2584  costs = subinst->costsfixed + mult.ubgreedylocal;
2585  SCIP_CALL( removeRedundantColumns(scip, heurdata, mult.xgreedylocal, &costs) );
2586 
2587  /* save solution if it is better */
2588  if( costs < heurdata->bestubsubinst )
2589  copySetCoverSolution(scip, core, subinst, heurdata->bestubsubinstsol, &heurdata->bestubsubinst, mult.xgreedylocal);
2590  }
2591 
2592  freeMemoryForSolution(scip, &mult);
2593 
2594  return SCIP_OKAY;
2595 }
2596 
2597 /* reports a solution to SCIP */
2598 static
2599 SCIP_RETCODE reportSolution(
2600  SCIP* scip, /**< master SCIP data structure */
2601  SCP_INSTANCE* inst, /**< (reduced) SCP instance */
2602  SCIP_Bool* solution, /**< solution for 'inst'; will be extended to global solution */
2603  SCIP_HEUR* heur /**< pointer to the heuristic's data */
2604  )
2605 {
2606  SCIP_VAR** solvars;
2607  SCIP_Real* solvals;
2608  int nsolvars;
2609  int i;
2610  SCIP_SOL* newsol;
2611  SCIP_Bool success = FALSE;
2612  SCIP_Bool foundsol = FALSE;
2613  SCIP_CONS** conss;
2614  int nconss;
2615 
2616  SCIP_CALL( SCIPgetVarsData(scip, &solvars, &nsolvars, NULL, NULL, NULL, NULL) );
2617  SCIP_CALL( SCIPcreateSol(scip, &newsol, heur) );
2618  SCIP_CALL( SCIPallocBufferArray(scip, &solvals, nsolvars) );
2619 
2620  for( i = 0; i < nsolvars; i++ )
2621  {
2622  if( isVarInSolution(solution, solvars[i]) )
2623  solvals[i] = 1.0;
2624  else if( (isFixedVariable(inst, solvars[i]) == TRUE) && (SCIPisZero(scip, SCIPvarGetObj(solvars[i])) == FALSE) )
2625  solvals[i] = 1.0;
2626  else
2627  solvals[i] = 0.0;
2628  }
2629 
2630  SCIP_CALL( SCIPsetSolVals(scip, newsol, nsolvars, solvars, solvals) );
2631 
2632  /* test all constraints and check if the activity is correct, adjust some free variable if necessary */
2633  conss = SCIPgetConss(scip);
2634  nconss = SCIPgetNConss(scip);
2635  for( i = 0; (i < nconss) && (foundsol == FALSE); i++ )
2636  {
2637  SCIP_Real lhs = 0.0, activity = 0.0;
2638  SCIP_Real *vals = NULL;
2639  SCIP_Bool valuesallones = FALSE;
2640 
2641  if( !strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[i])), "linear") )
2642  {
2643  lhs = SCIPgetLhsLinear(scip, conss[i]);
2644  activity = SCIPgetActivityLinear(scip, conss[i], newsol);
2645  vals = SCIPgetValsLinear(scip, conss[i]);
2646  }
2647  else if( !strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[i])), "logicor") )
2648  {
2649  valuesallones = TRUE;
2650  }
2651  else if( !strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[i])), "masterbranch") )
2652  {
2653  /* do nothing */
2654  SCIPdebugMessage("constraint %i is a masterbranch\n", SCIPconsGetPos(conss[i]));
2655  continue;
2656  }
2657  else
2658  {
2659  SCIPerrorMessage("constraint is '%s', can't handle this\n", SCIPconshdlrGetName(SCIPconsGetHdlr(conss[i])));
2660  SCIPABORT();
2661  }
2662 
2663  if( lhs > activity )
2664  {
2665  int j;
2666  int nvars;
2667  SCIP_Bool changed = FALSE;
2668  SCIP_VAR **vars;
2669 
2670 
2671  /*SCIPdebugMessage("constraint %i: left hand side is violated by %f\n", i, lhs - activity);*/
2672  SCIP_CALL( SCIPgetConsNVars(scip, conss[i], &nvars, &success) );
2673  SCIP_CALL( SCIPallocBufferArray(scip, &vars, nvars) );
2674  SCIP_CALL( SCIPgetConsVars(scip, conss[i], vars, nvars, &success) );
2675 
2676  for( j = 0; (changed == FALSE) && (j < nvars); j++ )
2677  {
2678  SCIP_Bool costszero = SCIPisZero(scip, SCIPvarGetObj(vars[j]));
2679  if( costszero == TRUE )
2680  {
2681  assert(valuesallones || vals != NULL);
2682 
2683  /* constraint is "logicor" */
2684  if( valuesallones )
2685  SCIP_CALL( SCIPincSolVal(scip, newsol, vars[j], lhs - activity) );
2686  /* constraint is "linear" */
2687  else
2688  {
2689  assert(vals != NULL);
2690  if( !SCIPisZero(scip, vals[j]) )
2691  SCIP_CALL( SCIPincSolVal(scip, newsol, vars[j], (lhs - activity) / vals[j]) );
2692  else
2693  SCIPdebugMessage("could not adjust activity\n");
2694  }
2695 
2696  SCIP_CALL( SCIPtrySol(scip, newsol, FALSE, FALSE, TRUE, TRUE, TRUE, &foundsol) );
2697  changed = TRUE;
2698  }
2699  }
2700 
2701  if( changed == FALSE )
2702  SCIPdebugMessage("could not find variable with zero costs\n");
2703 
2704  SCIPfreeBufferArray(scip, &vars);
2705  }
2706 
2707  if( !strcmp(SCIPconshdlrGetName(SCIPconsGetHdlr(conss[i])), "linear") )
2708  {
2709  if( SCIPgetLhsLinear(scip, conss[i]) > SCIPgetActivityLinear(scip, conss[i], newsol) )
2710  SCIPdebugMessage("activity is still smaller than lhs\n");
2711  }
2712  /* take care of the case rhs < activity. Question: can this occur?
2713  * YES it can, but only the convexity constraint can be violated (is this really true?).
2714  * We simply ignore this issue. SCIP will not accept the solution if it is not good enough.
2715  */
2716  }
2717 
2718  SCIP_CALL( SCIPtrySolFree(scip, &newsol, FALSE, FALSE, TRUE, TRUE, TRUE, &success) );
2719  SCIPfreeBufferArray(scip, &solvals);
2720 
2721  if( success == TRUE )
2722  SCIPdebugMessage("new solution found by set covering heuristic\n");
2723 
2724  return SCIP_OKAY;
2725 }
2726 
2727 /** the three-phase procedure from the paper. It tries to find near-optimal lagrange multipliers and
2728  * reduces the size of an instance by fixing variables that are likely to be in an optimal solution. */
2729 static
2730 SCIP_RETCODE threePhase(
2731  SCIP* scip, /**< master SCIP data structure */
2732  SCIP_HEUR* heur, /**< master SCIP heur data structure */
2733  SCIP_HEURDATA* heurdata /**< pointer to the heuristic's data */
2734  )
2735 {
2736  SCP_Lagrange_Sol* mult_lb_subinst;
2737  int i;
2738  SCIP_Bool success;
2739  SCP_CORE* core;
2740  SCP_INSTANCE* inst;
2741  SCP_INSTANCE* subinst;
2742  int niter;
2743 
2744  core = &heurdata->core;
2745  inst = &heurdata->inst;
2746  subinst = &heurdata->subinst;
2747  mult_lb_subinst = &heurdata->tpmultlbsubinst;
2748 
2749  /* we first create our own copy of the instance, as we need to mark variables as fixed until all variables are fixed */
2750  copyInstance(scip, core, subinst, inst);
2751  SCIP_CALL( markRowsCoveredByFixedVariables(scip, core, subinst, heurdata) );
2752 
2753  /* if this function is called for the first time, no multiplier is available and we need to compute one to start with */
2754  if( heurdata->useinitialmultiplier )
2755  {
2756  /* next, compute initial lagrange multipliers and find a first lower bound */
2757  SCIP_CALL( computeInitialLagrangeMultiplier(scip, core, subinst, mult_lb_subinst, heurdata) );
2758  heurdata->useinitialmultiplier = FALSE;
2759  }
2760  else
2761  {
2762  /* take the best lagrange multiplier of the unrestricted instance as starting point */
2763  copySolution(scip, core, mult_lb_subinst, &heurdata->multbestlbtotal);
2764  }
2765 
2766  /* computeOptimalSolution computes an optimal solution to the lagrange relaxation. Also, it computes the subgradient */
2767  SCIP_CALL( computeOptimalSolution(scip, core, subinst, mult_lb_subinst, heurdata) );
2768  SCIP_CALL( greedySetCover(scip, core, subinst, mult_lb_subinst, heurdata) );
2769 
2770  /* we now have a lower and upper bound in mult_lb_local for the instance 'inst' and take these as our starting values */
2771  copySetCoverSolution(scip, core, inst, heurdata->bestubinst_sol, &heurdata->bestubinst, mult_lb_subinst->xgreedylocal);
2772  copySetCoverSolution(scip, core, subinst, heurdata->bestubsubinstsol, &heurdata->bestubsubinst, mult_lb_subinst->xgreedylocal);
2773  copySolution(scip, core, &heurdata->multbestlbinst, mult_lb_subinst);
2774 
2775  /* check whether 'bestubinst_sol' is a solution of the reduced instance 'subinst' */
2776  SCIP_CALL( checkSetCover(scip, core, subinst, heurdata->bestubinst_sol, &success, heurdata) );
2777  if( !success )
2778  {
2779  SCIPerrorMessage("three-phase: initial solution is not a valid set cover\n");
2780  SCIPABORT();
2781  }
2782 
2783  /* in the first call of this proceduce there is no best upper bound for the unrestricted instance */
2784  if( heurdata->bestubinst < heurdata->bestub )
2785  {
2786  copySetCoverSolution(scip, core, inst, heurdata->bestubsol, &heurdata->bestub, heurdata->bestubinst_sol);
2787  SCIPdebugMessage("new upper bound: %f\n", heurdata->bestub);
2788  SCIP_CALL( reportSolution(scip, inst, heurdata->bestubsol, heur) );
2789  }
2790 
2791  if( core->nactiveconss <= subinst->nrowscovered )
2792  {
2793  SCIPdebugMessage("threephase: all rows are already covered\n");
2794  }
2795 
2796  /* The next loop works as follows:
2797  * 1. Derive the subinstance 'subinst' by marking all rows as covered that are covered by fixed variables of 'subinst'. If all rows are covered by fixed variables, stop.
2798  * 2. Perform subgradient optimization to find a lagrange multiplier that gives a better lower bound for 'subinst'.
2799  * 3. Explore the neighborhood of this best multiplier by doing additional subgradient steps and calling the greedy algorithm with each multiplier as argument.
2800  * 4. Check whether any of these local solutions can be extended to better global solutions.
2801  * 5. Fix several variables to reduce the size of 'subinst'.
2802  */
2803 
2804  /* stop if all rows are covered by fixed variables or a maximum number of iterations has been reached */
2805  for( niter = 0; niter < heurdata->param_threephase_max_iter && core->nactiveconss > subinst->nrowscovered; ++niter )
2806  {
2807  /* we mark all rows covered by fixed variables, in addition to the ones that were already covered */
2808  SCIP_CALL( markRowsCoveredByFixedVariables(scip, core, subinst, heurdata) );
2809  SCIP_CALL( subgradientOptimization(scip, core, subinst, mult_lb_subinst, heurdata->bestubsubinst, heurdata) );
2810 
2811  /* explore neighborhood of the best lower bound multiplier by applying the held-karp update again.
2812  * this updates all upper bounds (and sometimes also lower bounds) in 'heurdata' */
2813  SCIP_CALL( exploreNeighborhood(scip, mult_lb_subinst, heurdata) );
2814 
2815  /* extend a solution of 'subinst' to one of 'inst' if it gives a better upper bound */
2816  if( heurdata->bestubsubinst < heurdata->bestubinst )
2817  {
2818  copySetCoverSolution(scip, core, subinst, heurdata->bestubinst_sol, &heurdata->bestubinst, heurdata->bestubsubinstsol);
2819 
2820  /* extend a solution of 'inst' to a solution of the whole instance and check whether this gives a better global upper bound */
2821  if( heurdata->bestubinst < heurdata->bestub )
2822  {
2823  copySetCoverSolution(scip, core, inst, heurdata->bestubsol, &heurdata->bestub, heurdata->bestubinst_sol);
2824 
2825  SCIP_CALL( removeRedundantColumns(scip, heurdata, heurdata->bestubsol, &heurdata->bestub) );
2826  SCIPdebugMessage("new upper bound: %f\n", heurdata->bestub);
2827  SCIP_CALL( reportSolution(scip, inst, heurdata->bestubsol, heur) );
2828  }
2829  }
2830 
2831  /* stop if any solution can not be better than the best known solution */
2832  if( subinst->costsfixed + mult_lb_subinst->lblagrangelocal >= heurdata->bestub )
2833  break;
2834 
2835  /* At each step we fix new variables: fix the first max(1, nconstraints / 200) variables
2836  * that were picked by the greedy solution (in the order they were picked) */
2837  for( i = 0; i < core->nsolgreedy && i <= core->nconss / 200; i++ )
2838  {
2839  fixVariable(subinst, core->variables[core->solgreedy[i]]);
2840  subinst->costsfixed += SCIPvarGetObj(core->variables[core->solgreedy[i]]);
2841  }
2842  }
2843 
2844  SCIP_CALL( checkSetCover(scip, core, inst, heurdata->bestubinst_sol, &success, heurdata) );
2845  if( !success )
2846  {
2847  SCIPdebugMessage("three-phase: final solution is not a valid set cover\n");
2848  }
2849 
2850  return SCIP_OKAY;
2851 }
2852 
2853 /** driver for the three-phase procedure */
2854 static
2856  SCIP* scip, /**< master SCIP data structure */
2857  SCIP_HEUR* heur /**< original SCIP heuristic data structure */
2858  )
2859 {
2860  SCIP_HEURDATA* heurdata;
2861  SCIP_Bool stopcft = FALSE;
2862  int niter = 0; /* total number of iterations of the three-phase loop */
2863  int nitercore = 0; /* total number of iterations on the current core */
2864  int niternoimp = 0; /* number of iterations for how long no improvement on the upper bound was made */
2865  int coret = 10; /* compute new core after this number of iterations */
2866  SCIP_Real pi; /* percentage of rows that need to be covered in each fixing phase */
2867  SCIP_Real corelb = 0.0;
2868  SCP_CORE* core;
2869  SCP_INSTANCE* inst;
2870  SCIP_Bool success;
2871 
2872  int i;
2873 
2874  heurdata = SCIPheurGetData(heur);
2875  assert(heurdata != NULL);
2876 
2877  /* allocate memory that is used locally by redefineCore */
2878  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->rccols, heurdata->param_core_tent_size) );
2879  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->rccoldelta, heurdata->param_core_tent_size) );
2880 
2881  /* allocate memory that is used locally by subgradientOptimization */
2882  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->sglastlb, heurdata->param_lambda_p) );
2883 
2884  /* we always consider core variables (columns) only to cover rows, so we need to define the first core.
2885  * This is done by choosing for each row the first five columns covering it and adding them to the core. */
2886  SCIP_CALL( initTentativeCore(scip, &heurdata->core, heurdata) );
2887 
2888  /* allocate memory that is used to get all variables that are part of some constraint.
2889  * 'maxconstraintvariables' is the maximum number of variables any constraint of the problem contains. */
2890  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->vars, heurdata->core.maxconstraintvariables) );
2891 
2892  /* after defining the core, we compute all core columns and core rows */
2893  SCIP_CALL( computeCoreColumns(scip, &heurdata->core) );
2894  SCIP_CALL( computeCoreRows(scip, &heurdata->core, heurdata) );
2895 
2896  /* set up basic instance. so far no variables are fixed */
2897  SCIP_CALL( initInstance(scip, &heurdata->inst) );
2898  SCIP_CALL( initInstance(scip, &heurdata->subinst) );
2899 
2900  core = &heurdata->core;
2901  inst = &heurdata->inst;
2902  pi = heurdata->param_pi_min;
2903 
2904  /* allocate memory that is used by the greedy algorithm locally */
2905  SCIP_CALL( pqueue_init(scip, &heurdata->greedyqueue) );
2906  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->greedycolpos, core->nvariables) );
2907  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->greedycolmu, core->nvariables) );
2908  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->greedycolgamma, core->nvariables) );
2909  SCIP_CALL( SCIPallocBufferArray(scip, &heurdata->greedycolscore, core->nvariables) );
2910  SCIP_CALL( initInstance(scip, &heurdata->greedyinst) );
2911 
2912  SCIP_CALL( allocateMemoryForSolution(scip, &heurdata->core, &heurdata->multbestlbinst) );
2913  SCIP_CALL( allocateMemoryForSolution(scip, &heurdata->core, &heurdata->multbestlbsubinst) );
2914  SCIP_CALL( allocateMemoryForSolution(scip, &heurdata->core, &heurdata->multbestlbtotal) );
2915  SCIP_CALL( allocateMemoryForSolution(scip, &heurdata->core, &heurdata->tpmultlbsubinst) );
2916 
2917  SCIP_CALL( SCIPallocClearBufferArray(scip, &heurdata->bestubinst_sol, SCIPgetNVars(scip)) );
2918  SCIP_CALL( SCIPallocClearBufferArray(scip, &heurdata->bestubsubinstsol, SCIPgetNVars(scip)) );
2919  SCIP_CALL( SCIPallocClearBufferArray(scip, &heurdata->bestubsol, SCIPgetNVars(scip)) );
2920 
2921  heurdata->multbestlbtotal.lblagrangeglobal = 0;
2922  heurdata->bestub = SCIP_REAL_MAX;
2923  heurdata->useinitialmultiplier = TRUE;
2924 
2925  while( stopcft == FALSE )
2926  {
2927  SCIP_Bool redefine = FALSE;
2928  SCIP_Real currentub = heurdata->bestub;
2929 
2930 
2931  /* derive the reduced subinstance by marking rows covered by fixed variables */
2932  for( i = 0; i < SCIPgetNConss(scip); ++i )
2933  inst->rowscovered[i] = FALSE;
2934  inst->nrowscovered = 0;
2935  SCIP_CALL( markRowsCoveredByFixedVariables(scip, core, inst, heurdata) );
2936 
2937  /* call three-phase as long as 'inst' contains uncovered rows */
2938  if( core->nactiveconss > inst->nrowscovered )
2939  {
2940  SCIPdebugMessage("%d variables are fixed, %d rows are covered\n", inst->nvarsfixed, inst->nrowscovered);
2941 
2942  /* apply procedure three-phase to find an near optimal lagrange multiplier */
2943  SCIP_CALL( threePhase(scip, heur, heurdata) );
2944  }
2945  else
2946  {
2947  /* compute new core when all of its rows are covered by fixed variables */
2948  redefine = TRUE;
2949  }
2950 
2951  /* stop if maximum number of iterations is reached */
2952  if( niter++ >= heurdata->param_max_iter )
2953  stopcft = TRUE;
2954 
2955  /* set pi to PI_MIN if the current best solution was found in this iteration */
2956  if( heurdata->bestub < currentub )
2957  {
2958  niternoimp = 0;
2959  pi = heurdata->param_pi_min;
2960  }
2961  else
2962  {
2963  /* increase pi if no better solution was found, i.e. fix more variables in order to cover more rows */
2964  pi = pi * heurdata->param_pi_alpha;
2965  niternoimp++;
2966  }
2967 
2968  /* stop if there was no improvement during the last 'SCP_MAX_ITER_NO_IMP' iterations */
2969  if( niternoimp >= heurdata->param_max_iter_no_imp )
2970  stopcft = TRUE;
2971 
2972  /* stop if UB <= beta * LB */
2973  if( heurdata->multbestlbtotal.lblagrangeglobal * heurdata->param_beta >= heurdata->bestub )
2974  stopcft = TRUE;
2975 
2976  /* redefine core if the current core was worked on for 'coret' iterations */
2977  if( nitercore++ == coret )
2978  redefine = TRUE;
2979 
2980  if( stopcft )
2981  break;
2982 
2983  /* compute delta fixes variables of 'bestubsol' within the core that are likely to be part of an optimal solution */
2984  SCIP_CALL( markRowsCoveredByFixedVariables(scip, core, inst, heurdata) );
2985  SCIP_CALL( computeDelta(scip, core, inst, heurdata->multbestlbtotal.lagrangiancostsglobal, heurdata->bestubsol, pi, heurdata) );
2986 
2987  if( redefine == TRUE )
2988  {
2989  /* stop if the last core did not lead to any improvements */
2990  if( SCIPisEQ(scip, corelb, heurdata->multbestlbtotal.lblagrangeglobal) == TRUE )
2991  stopcft = TRUE;
2992  else
2993  {
2994  SCIP_CALL( redefineCore(scip, heurdata) );
2995  for( i = 0; i < SCIPgetNVars(scip); ++i )
2996  inst->varsfixed[i] = FALSE;
2997  inst->nvarsfixed = 0;
2998  inst->costsfixed = 0.0;
2999  pi = heurdata->param_pi_min;
3000  nitercore = 0;
3001  corelb = heurdata->multbestlbtotal.lblagrangeglobal;
3002  }
3003  }
3004 
3005  SCIPdebugMessage("iteration %i: best lower bound: %f, best upper bound: %f\n", niter, heurdata->multbestlbtotal.lblagrangeglobal, heurdata->bestub);
3006  }
3007 
3008  for( i = 0; i < SCIPgetNVars(scip); ++i )
3009  inst->varsfixed[i] = FALSE;
3010  inst->nvarsfixed = 0;
3011 
3012  SCIP_CALL( checkSetCover(scip, core, inst, heurdata->bestubsol, &success, heurdata) );
3013  if( success == FALSE )
3014  SCIPdebugMessage("final solution is not a valid set cover\n");
3015  else
3016  SCIPdebugMessage("final solution has costs %f\n", heurdata->bestub);
3017 
3018  /* report the final solution to SCIP */
3019  SCIP_CALL( reportSolution(scip, inst, heurdata->bestubsol, heur) );
3020 
3021  /* release all memory that was allocated before */
3022  SCIPfreeBufferArray(scip, &heurdata->bestubsol);
3023  SCIPfreeBufferArray(scip, &heurdata->bestubsubinstsol);
3024  SCIPfreeBufferArray(scip, &heurdata->bestubinst_sol);
3025 
3026  freeMemoryForSolution(scip, &heurdata->tpmultlbsubinst);
3027  freeMemoryForSolution(scip, &heurdata->multbestlbtotal);
3028  freeMemoryForSolution(scip, &heurdata->multbestlbsubinst);
3029  freeMemoryForSolution(scip, &heurdata->multbestlbinst);
3030 
3031  freeInstance(scip, &heurdata->greedyinst);
3032  SCIPfreeBufferArray(scip, &heurdata->greedycolscore);
3033  SCIPfreeBufferArray(scip, &heurdata->greedycolgamma);
3034  SCIPfreeBufferArray(scip, &heurdata->greedycolmu);
3035  SCIPfreeBufferArray(scip, &heurdata->greedycolpos);
3036 
3037  pqueue_destroy(scip, &heurdata->greedyqueue);
3038 
3039  freeInstance(scip, &heurdata->subinst);
3040  freeInstance(scip, &heurdata->inst);
3041 
3042  SCIPfreeBufferArray(scip, &heurdata->vars);
3043 
3044  freeCore(scip, &heurdata->core);
3045 
3046  SCIPfreeBufferArray(scip, &heurdata->sglastlb);
3047 
3048  SCIPfreeBufferArray(scip, &heurdata->rccoldelta);
3049  SCIPfreeBufferArray(scip, &heurdata->rccols);
3050 
3051  return SCIP_OKAY;
3052 }
3053 
3054 /*
3055  * Callback methods of primal heuristic
3056  */
3057 
3058 /** destructor of primal heuristic to free user data (called when SCIP is exiting) */
3059 static
3060 SCIP_DECL_HEURFREE(heurFreeSetcover)
3061 { /*lint --e{715}*/
3062  SCIP_HEURDATA* heurdata;
3063 
3064  assert(heur != NULL);
3065  assert(scip != NULL);
3066 
3067  /* get heuristic data */
3068  heurdata = SCIPheurGetData(heur);
3069  assert(heurdata != NULL);
3070 
3071  /* free heuristic data */
3072  SCIPfreeMemory(scip, &heurdata);
3073  SCIPheurSetData(heur, NULL);
3074 
3075  return SCIP_OKAY;
3076 }
3077 
3078 /** initialization method of primal heuristic (called after problem was transformed) */
3079 static
3080 SCIP_DECL_HEURINIT(heurInitSetcover)
3081 { /*lint --e{715}*/
3082  SCIP_HEURDATA* heurdata;
3083 
3084  assert(heur != NULL);
3085  assert(scip != NULL);
3086 
3087  /* get heuristic data */
3088  heurdata = SCIPheurGetData(heur);
3089  assert(heurdata != NULL);
3090 
3091  /* create random number generator */
3092  SCIP_CALL( SCIPcreateRandom(scip, &heurdata->randnumgen,
3093  SCIPinitializeRandomSeed(scip, DEFAULT_RANDSEED), TRUE) );
3094 
3095  return SCIP_OKAY;
3096 }
3097 
3098 /** deinitialization method of primal heuristic (called before transformed problem is freed) */
3099 static
3100 SCIP_DECL_HEUREXIT(heurExitSetcover)
3101 { /*lint --e{715}*/
3102  SCIP_HEURDATA* heurdata;
3103 
3104  assert(heur != NULL);
3105  assert(scip != NULL);
3106 
3107  /* get heuristic data */
3108  heurdata = SCIPheurGetData(heur);
3109  assert(heurdata != NULL);
3110 
3111  /* free random number generator */
3112  SCIPfreeRandom(scip, &heurdata->randnumgen);
3113 
3114  return SCIP_OKAY;
3115 }
3116 
3117 
3118 /** execution method of primal heuristic */
3119 static
3120 SCIP_DECL_HEUREXEC(heurExecSetcover)
3121 { /*lint --e{715}*/
3122  SCIP* origprob;
3123  SCIP_HEURDATA* heurdata;
3124 
3125  assert(heur != NULL);
3126  assert(scip != NULL);
3127  assert(result != NULL);
3128 
3129  heurdata = SCIPheurGetData(heur);
3130  assert(heurdata != NULL);
3131 
3132  origprob = GCGmasterGetOrigprob(scip);
3133  assert(origprob != NULL);
3134 
3135  *result = SCIP_DIDNOTRUN;
3136 
3137  /* the heuristic is only executed on problems with many variables */
3138  if( SCIPgetNVars(scip) < heurdata->param_min_prob_size )
3139  {
3140  SCIPdebugMessage("not running set covering heuristic because instance is too small (only %i variables)\n", SCIPgetNVars(scip));
3141  return SCIP_OKAY;
3142  }
3143 
3144  /* only run on set covering problems */
3145  if( GCGisMasterSetCovering(origprob) == FALSE )
3146  return SCIP_OKAY;
3147 
3148  if( SCIPgetNVars(scip) == 0 )
3149  return SCIP_OKAY;
3150 
3151  SCIP_CALL( setCoveringHeuristic(scip, heur) );
3152 
3153  *result = SCIP_FOUNDSOL;
3154  return SCIP_OKAY;
3155 }
3156 
3157 
3158 /*
3159  * primal heuristic specific interface methods
3160  */
3161 
3162 /** creates the setcover primal heuristic and includes it in SCIP */
3164  SCIP* scip /**< SCIP data structure */
3165  )
3166 {
3167  SCIP_HEURDATA* heurdata = NULL;
3168  SCIP_HEUR* heur;
3169 
3170  heur = NULL;
3171  SCIP_CALL( SCIPallocMemory(scip, &heurdata) );
3172 
3173  /* include primal heuristic */
3174 
3175  SCIP_CALL( SCIPincludeHeurBasic(scip, &heur,
3177  HEUR_MAXDEPTH, HEUR_TIMING, HEUR_USESSUBSCIP, heurExecSetcover, heurdata) );
3178 
3179  assert(heur != NULL);
3180 
3181  /* set non fundamental callbacks via setter functions */
3182  SCIP_CALL( SCIPsetHeurFree(scip, heur, heurFreeSetcover) );
3183  SCIP_CALL( SCIPsetHeurInit(scip, heur, heurInitSetcover) );
3184  SCIP_CALL( SCIPsetHeurExit(scip, heur, heurExitSetcover) );
3185 
3186  /* add setcover primal heuristic parameters */
3187 
3188  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/tentcoresize",
3189  "how many columns are added to the tentative core for each row",
3190  &heurdata->param_core_tent_size, FALSE, DEF_CORE_TENT_SIZE, 1, INT_MAX, NULL, NULL) );
3191 
3192  SCIP_CALL( SCIPaddBoolParam(scip, "heuristics/"HEUR_NAME"/lambdaadjustments",
3193  "should the step size during the subgradient phase be adjusted?",
3194  &heurdata->param_lambda_adjustments, FALSE, DEF_LAMBDA_ADJUSTMENTS, NULL, NULL) );
3195 
3196  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/lambda",
3197  "default step size",
3198  &heurdata->param_lambda, FALSE, DEF_LAMBDA, 0.001, SCIP_REAL_MAX, NULL, NULL) );
3199 
3200  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/lambdap",
3201  "number of iterations after which lambda is adjusted",
3202  &heurdata->param_lambda_p, FALSE, DEF_LAMBDA_P, 1, INT_MAX, NULL, NULL) );
3203 
3204  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/beta",
3205  "greatest gap between lower and upper bounds that is allowed",
3206  &heurdata->param_beta, FALSE, DEF_BETA, 1.0, SCIP_REAL_MAX, NULL, NULL) );
3207 
3208  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/sciter",
3209  "number of subgradient iterations after which the stopping criterion is tested again",
3210  &heurdata->param_stop_crit_iter, FALSE, DEF_STOP_CRIT_ITER, 1, INT_MAX, NULL, NULL) );
3211 
3212  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/scdiff",
3213  "maximum absolute difference between lower and upper bound that is allowed for the stopping criterion",
3214  &heurdata->param_stop_crit_diff, FALSE, DEF_STOP_CRIT_DIFF, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3215 
3216  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/scratio",
3217  "minimum ratio between lower and upper bound that is allowed that is allowed for the stopping criterion",
3218  &heurdata->param_stop_crit_ratio, FALSE, DEF_STOP_CRIT_RATIO, 0.0, 1.0, NULL, NULL) );
3219 
3220  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/pimin",
3221  "percentage of rows to be removed after column fixing",
3222  &heurdata->param_pi_min, FALSE, DEF_PI_MIN, 0.0, 1.0, NULL, NULL) );
3223 
3224  SCIP_CALL( SCIPaddRealParam(scip, "heuristics/"HEUR_NAME"/pialpha",
3225  "increase of pi when no improvement was made",
3226  &heurdata->param_pi_alpha, FALSE, DEF_PI_ALPHA, 0.0, SCIP_REAL_MAX, NULL, NULL) );
3227 
3228  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/tpmaxiter",
3229  "maximum iterations of the three-phase procedure",
3230  &heurdata->param_max_iter, FALSE, DEF_MAX_ITER, 1, INT_MAX, NULL, NULL) );
3231 
3232  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/tpmaxiternoimp",
3233  "maximum allowed number of iterations of the three-phase procedure without improvements",
3234  &heurdata->param_max_iter_no_imp, FALSE, DEF_MAX_ITER_NO_IMP, 1, INT_MAX, NULL, NULL) );
3235 
3236  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/threephasemaxiter",
3237  "maximum number of iterations inside three-phase",
3238  &heurdata->param_threephase_max_iter, FALSE, DEF_THREEPHASE_MAX_ITER, 1, INT_MAX, NULL, NULL) );
3239 
3240  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/greedymaxiter",
3241  "maximum number of iterations of the greedy algorithm",
3242  &heurdata->param_greedy_max_iter, FALSE, DEF_GREEDY_MAX_ITER, 1, INT_MAX, NULL, NULL) );
3243 
3244  SCIP_CALL( SCIPaddIntParam(scip, "heuristics/"HEUR_NAME"/minprobsize",
3245  "minimum number of variables in the problem for the heuristic to start",
3246  &heurdata->param_min_prob_size, FALSE, DEF_MIN_PROB_SIZE, 1, INT_MAX, NULL, NULL) );
3247 
3248  return SCIP_OKAY;
3249 }
#define HEUR_NAME
Definition: heur_setcover.c:45
static void clearSolution(SCIP *scip, SCIP_Bool *solution)
SCIP_HASHMAP * mapvariables
static SCIP_RETCODE checkSetCover(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCIP_Bool *solution, SCIP_Bool *issetcover, SCIP_HEURDATA *heurdata)
SCIP_Real * rccoldelta
static SCIP_Bool isFixedVariable(SCP_INSTANCE *inst, SCIP_VAR *var)
static SCIP_Bool isVarInSolution(SCIP_Bool *solution, SCIP_VAR *var)
#define DEF_CORE_TENT_SIZE
Definition: heur_setcover.c:55
GCG interface methods.
SCIP_Real * greedycolscore
SCIP_Real lblagrangelocal
static SCIP_DECL_HEUREXEC(heurExecSetcover)
#define DEF_MIN_PROB_SIZE
Definition: heur_setcover.c:69
int ncorevariables
SCIP_Real * delta
SCIP_Real * subgradient
static void copySolution(SCIP *scip, SCP_CORE *core, SCP_Lagrange_Sol *dest, SCP_Lagrange_Sol *source)
#define DEF_LAMBDA
Definition: heur_setcover.c:58
static void copySetCoverSolution(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCIP_Bool *dest, SCIP_Real *destcosts, SCIP_Bool *source)
static SCIP_RETCODE computeCoreRows(SCIP *scip, SCP_CORE *core, SCIP_HEURDATA *heurdata)
SCP_Lagrange_Sol multbestlbinst
static void freeMemoryForSolution(SCIP *scip, SCP_Lagrange_Sol *mult)
static void fixVariable(SCP_INSTANCE *inst, SCIP_VAR *var)
SCIP_VAR ** vars
static SCIP_RETCODE greedySetCover(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCP_Lagrange_Sol *mult, SCIP_HEURDATA *heurdata)
int * data
Definition: heur_setcover.c:83
#define HEUR_TIMING
Definition: heur_setcover.c:52
static SCIP_RETCODE pqueue_init(SCIP *scip, PQUEUE *queue)
static SCIP_Bool isRowCovered(SCP_INSTANCE *inst, int rowpos)
static void addVarToCore(SCP_CORE *core, SCIP_VAR *var)
SCIP_Bool useinitialmultiplier
static void pqueue_destroy(SCIP *scip, PQUEUE *queue)
static SCIP_RETCODE markRowsCoveredByFixedVariables(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCIP_HEURDATA *heurdata)
SCIP_Real * sglastlb
#define HEUR_MAXDEPTH
Definition: heur_setcover.c:51
SCIP_Real bestubinst
#define DEF_MAX_ITER_NO_IMP
Definition: heur_setcover.c:66
SCIP_Bool param_lambda_adjustments
#define DEF_STOP_CRIT_ITER
Definition: heur_setcover.c:59
#define DEF_LAMBDA_ADJUSTMENTS
Definition: heur_setcover.c:56
static SCIP_RETCODE computeDelta(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCIP_Real *lagrangiancosts, SCIP_Bool *solution, SCIP_Real pi, SCIP_HEURDATA *heurdata)
SCIP_Real * keys
Definition: heur_setcover.c:82
static SCIP_DECL_HEURFREE(heurFreeSetcover)
static SCIP_RETCODE setCoveringHeuristic(SCIP *scip, SCIP_HEUR *heur)
SCIP_RETCODE SCIPincludeHeurSetcover(SCIP *scip)
#define DEF_THREEPHASE_MAX_ITER
Definition: heur_setcover.c:67
static SCIP_RETCODE computeOptimalSolution(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCP_Lagrange_Sol *mult, SCIP_HEURDATA *heurdata)
static void removeVarsFromCore(SCIP *scip, SCP_CORE *core)
static SCIP_RETCODE threePhase(SCIP *scip, SCIP_HEUR *heur, SCIP_HEURDATA *heurdata)
SCP_INSTANCE subinst
#define HEUR_FREQ
Definition: heur_setcover.c:49
SCIP_Real param_beta
SCIP_RANDNUMGEN * randnumgen
GCG variable pricer.
SCIP_Real costsfixed
Definition: heur_setcover.c:93
int size
Definition: heur_setcover.c:80
static SCIP_RETCODE initInstance(SCIP *scip, SCP_INSTANCE *inst)
static void pqueue_get_min(SCIP *scip, PQUEUE *queue, int *elem)
#define HEUR_PRIORITY
Definition: heur_setcover.c:48
int ** positions
Definition: heur_setcover.c:84
static SCIP_Bool isCoreVariable(SCP_CORE *core, SCIP_VAR *var)
#define DEF_STOP_CRIT_DIFF
Definition: heur_setcover.c:60
int * constraintid
static void copyInstance(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *dest, SCP_INSTANCE *source)
#define HEUR_DISPCHAR
Definition: heur_setcover.c:47
#define DEF_LAMBDA_P
Definition: heur_setcover.c:57
int * nvarconss
SCP_INSTANCE greedyinst
static void freeInstance(SCIP *scip, SCP_INSTANCE *inst)
static SCIP_RETCODE computeInitialLagrangeMultiplier(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCP_Lagrange_Sol *mult, SCIP_HEURDATA *heurdata)
SCIP_Bool * bestubinst_sol
static SCIP_RETCODE computeCoreColumns(SCIP *scip, SCP_CORE *core)
SCP_Lagrange_Sol multbestlbtotal
SCIP_Bool * bestubsubinstsol
static SCIP_RETCODE getConsVars(SCIP *scip, SCP_CORE *core, int pos, SCIP_VAR **vars, int *nvars, SCIP_Bool *success)
SCIP_Real ubgreedylocal
int nsolgreedy
SCIP_Real bestub
static void markRowAsCovered(SCP_INSTANCE *inst, int rowpos)
#define HEUR_DESC
Definition: heur_setcover.c:46
int * nrowvars
static void pqueue_decrease_key(SCIP *scip, PQUEUE *queue, int pos, SCIP_Real key)
#define DEF_MAX_ITER
Definition: heur_setcover.c:65
SCIP_Bool * rowscovered
Definition: heur_setcover.c:91
static SCIP_RETCODE exploreNeighborhood(SCIP *scip, SCP_Lagrange_Sol *startmult, SCIP_HEURDATA *heurdata)
SCIP * GCGmasterGetOrigprob(SCIP *scip)
static SCIP_DECL_HEUREXIT(heurExitSetcover)
SCIP_Real lblagrangeglobal
SCIP_Real * lagrangiancostsglobal
#define DEF_PI_ALPHA
Definition: heur_setcover.c:63
int * solgreedy
int param_stop_crit_iter
static SCIP_RETCODE computeLocalLagrangianCosts(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCP_Lagrange_Sol *mult, SCIP_HEURDATA *heurdata)
static SCIP_RETCODE subgradientOptimization(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCP_Lagrange_Sol *best_mult_lb, SCIP_Real bestub, SCIP_HEURDATA *heurdata)
static SCIP_RETCODE computeGlobalLagrangianCosts(SCIP *scip, SCP_CORE *core, SCP_Lagrange_Sol *mult, SCIP_HEURDATA *heurdata)
static SCIP_RETCODE getVarIndex(SCIP *scip, SCP_CORE *core, SCIP_VAR *variable, int *pos)
static SCIP_DECL_HEURINIT(heurInitSetcover)
#define DEF_GREEDY_MAX_ITER
Definition: heur_setcover.c:68
SCP_INSTANCE inst
SCIP_Real param_pi_alpha
SCIP_Real param_lambda
static void freeCore(SCIP *scip, SCP_CORE *core)
#define DEF_BETA
Definition: heur_setcover.c:64
SCIP_Real param_stop_crit_ratio
int ** rows
static void addVarToSolution(SCIP_Bool *solution, SCIP_VAR *var)
SCP_Lagrange_Sol tpmultlbsubinst
set covering primal heuristic according to Caprara, Fischetti, and Toth (1999)
static SCIP_RETCODE initTentativeCore(SCIP *scip, SCP_CORE *core, SCIP_HEURDATA *heurdata)
int reserved
Definition: heur_setcover.c:81
SCIP_Real bestubsubinst
SCIP_VAR ** variables
#define DEFAULT_RANDSEED
Definition: heur_setcover.c:70
SCIP_Real param_stop_crit_diff
int * delta_pos
int nactiveconss
int * listcorevariables
Definition: heur_setcover.c:99
int param_max_iter_no_imp
PQUEUE greedyqueue
#define HEUR_USESSUBSCIP
Definition: heur_setcover.c:53
SCIP_Bool * xgreedylocal
#define DEF_STOP_CRIT_RATIO
Definition: heur_setcover.c:61
static void extendSolution(SCIP *scip, SCP_CORE *core, SCP_INSTANCE *inst, SCIP_Bool *solution)
SCIP_Bool * varincore
Definition: heur_setcover.c:98
SCIP_Real * greedycolgamma
SCIP_Bool * varsfixed
Definition: heur_setcover.c:89
static SCIP_RETCODE removeRedundantColumns(SCIP *scip, SCIP_HEURDATA *heurdata, SCIP_Bool *solution, SCIP_Real *solcosts)
int param_core_tent_size
#define HEUR_FREQOFS
Definition: heur_setcover.c:50
static void pqueue_increase_key(SCIP *scip, PQUEUE *queue, int pos, SCIP_Real key)
static SCIP_RETCODE redefineCore(SCIP *scip, SCIP_HEURDATA *heurdata)
SCIP_CONS ** conss
SCIP_Bool * bestubsol
SCIP_Bool rowsavailable
SCIP_Bool GCGisMasterSetCovering(SCIP *scip)
Definition: relax_gcg.c:4221
static void pqueue_clear(SCIP *scip, PQUEUE *queue)
SCP_Lagrange_Sol multbestlbsubinst
SCIP_Real param_pi_min
SCIP_Bool columnsavailable
SCIP_Real * lagrangiancostslocal
static SCIP_RETCODE pqueue_insert(SCIP *scip, PQUEUE *queue, SCIP_Real key, int elem, int *position)
int param_greedy_max_iter
int param_threephase_max_iter
int maxconstraintvariables
static SCIP_RETCODE allocateMemoryForSolution(SCIP *scip, SCP_CORE *core, SCP_Lagrange_Sol *mult)
static void removeVarFromSolution(SCIP_Bool *solution, SCIP_VAR *var)
int nvariables
static SCIP_RETCODE reportSolution(SCIP *scip, SCP_INSTANCE *inst, SCIP_Bool *solution, SCIP_HEUR *heur)
int ** columns
#define DEF_PI_MIN
Definition: heur_setcover.c:62