Scippy

SCIP

Solving Constraint Integer Programs

sepa_cgmip.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 /* #define SCIP_WRITEPROB */
16 /* #define SCIP_OUTPUT */
17 /**@file sepa_cgmip.c
18  * @brief Chvatal-Gomory cuts computed via a sub-MIP
19  * @author Marc Pfetsch
20  *
21  * Separate Chvátal-Gomory cuts using a sub-MIP. The approach is based on the following papers.
22  *
23  * M. Fischetti and A. Lodi@n
24  * Optimizing over the first Chvátal closure,@n
25  * in: M. Jünger and V. Kaibel (eds.) Integer Programming and Combinatorial Optimization IPCO 2005,@n
26  * LNCS 3509, pp. 12-22. Springer, Berlin Heidelberg New York (2005)
27  *
28  * M. Fischetti and A. Lodi@n
29  * Optimizing over the first Chvátal closure,@n
30  * Mathematical Programming 110, 3-20 (2007)
31  *
32  * P. Bonami, G. Cornuéjols, S. Dash, M. Fischetti, and A. Lodi@n
33  * Projected Chvátal-Gomory cuts for mixed integer linear programs,@n
34  * Mathematical Programming 113, No. 2 (2008)
35  *
36  *
37  * There are several versions to generate the final cut:
38  *
39  * - The CMIR-routines of SCIP can be used (if @p usecmir is true). One can determine which bound is
40  * used in the rounding operation (if cmirownbounds is true) or let SCIP choose the best. This
41  * version is generally numerically the most stable.
42  * - If @p usestrongcg is true, we try to generate Strong-CG cuts (as done in sepa_strongcg.c).
43  * - One can directly generate the CG-cut as computed (if @p usecmir and @p usestrongcg are
44  * false). The cut is not take from the solution of the MIP, but is recomputed, and some care (but
45  * not as much as in the first version) has been taken to create a valid cut.
46  *
47  * The computation time of the separation MIP is limited as follows:
48  * - There is a node limit (parameters @a minnodelimit and @a maxnodelimit).
49  * - There is a time limit (parameter @a timelimit).
50  * - If paramter @a earlyterm is true, the separation is run until the first cut that is violated is
51  * found. (Note that these cuts are not necessarily added to the LP, because here also the norm of
52  * the cuts are taken into account - which cannot easily be included into the separation subscip.)
53  * Then the solution is continued for a certain number of nodes.
54  *
55  * @todo Check whether one can weaken the conditions on the continuous variables.
56  * @todo Use pointers to originating separators to sort out cuts that should not be used.
57  *
58  * @warning This separator should be used carefully - it may require a long separation time.
59  */
60 
61 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
62 
63 #include <assert.h>
64 #include <string.h>
65 
66 #include "scip/sepa_cgmip.h"
67 #include "scip/scipdefplugins.h"
68 #include "scip/cons_linear.h"
69 #include "scip/pub_misc.h"
70 #include "scip/pub_lp.h"
71 
72 
73 #define SEPA_NAME "cgmip"
74 #define SEPA_DESC "Chvatal-Gomory cuts via MIPs separator"
75 #define SEPA_PRIORITY -1000
76 #define SEPA_FREQ -1
77 #define SEPA_MAXBOUNDDIST 0.0
78 #define SEPA_USESSUBSCIP TRUE /**< does the separator use a secondary SCIP instance? */
79 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
80 
81 #define DEFAULT_MAXROUNDS 5 /**< maximal number of separation rounds per node (-1: unlimited) */
82 #define DEFAULT_MAXROUNDSROOT 50 /**< maximal number of separation rounds in the root node (-1: unlimited) */
83 #define DEFAULT_MAXDEPTH -1 /**< maximal depth at which the separator is applied */
84 #define DEFAULT_DECISIONTREE FALSE /**< Use decision tree to turn separation on/off? */
85 #define DEFAULT_TIMELIMIT 1e20 /**< time limit for sub-MIP (set to infinity in order to be deterministic) */
86 #define DEFAULT_MEMORYLIMIT 1e20 /**< memory limit for sub-MIP */
87 #define DEFAULT_CUTCOEFBND 1000.0 /**< bounds on the values of the coefficients in the CG-cut */
88 #define DEFAULT_MINNODELIMIT 500LL /**< minimum number of nodes considered for sub-MIP (-1: unlimited) */
89 #define DEFAULT_MAXNODELIMIT 5000LL /**< maximum number of nodes considered for sub-MIP (-1: unlimited) */
90 #define DEFAULT_ONLYACTIVEROWS FALSE /**< Use only active rows to generate cuts? */
91 #define DEFAULT_MAXROWAGE -1 /**< maximal age of rows to consider if onlyactiverows is false */
92 #define DEFAULT_ONLYRANKONE FALSE /**< Separate rank 1 inequalities w.r.t. CG-MIP separator? */
93 #define DEFAULT_ONLYINTVARS FALSE /**< Generate cuts for problems with only integer variables? */
94 #define DEFAULT_ALLOWLOCAL FALSE /**< Allow to generate local cuts? */
95 #define DEFAULT_CONTCONVERT FALSE /**< Convert some integral variables to be continuous to reduce the size of the sub-MIP? */
96 #define DEFAULT_CONTCONVFRAC 0.1 /**< fraction of integral variables converted to be continuous (if contconvert) */
97 #define DEFAULT_CONTCONVMIN 100 /**< minimum number of integral variables before some are converted to be continuous */
98 #define DEFAULT_INTCONVERT FALSE /**< Convert some integral variables attaining fractional values to have integral value? */
99 #define DEFAULT_INTCONVFRAC 0.1 /**< fraction of fractional integral variables converted to have integral value (if intconvert) */
100 #define DEFAULT_INTCONVMIN 100 /**< minimum number of integral variables before some are converted to have integral value */
101 #define DEFAULT_SKIPMULTBOUNDS TRUE /**< Skip the upper bounds on the multipliers in the sub-MIP? */
102 #define DEFAULT_OBJLONE FALSE /**< Should the objective of the sub-MIP only minimize the l1-norm of the multipliers? */
103 #define DEFAULT_OBJWEIGHT 1e-03 /**< objective weight for artificial variables */
104 #define DEFAULT_OBJWEIGHTSIZE TRUE /**< Weight each row by its size? */
105 #define DEFAULT_DYNAMICCUTS TRUE /**< Should generated cuts be removed from the LP if they are no longer tight? */
106 #define DEFAULT_USECMIR TRUE /**< Use CMIR-generator (otherwise add cut directly)? */
107 #define DEFAULT_USESTRONGCG FALSE /**< Use strong CG-function to strengthen cut? */
108 #define DEFAULT_CMIROWNBOUNDS FALSE /**< Tell CMIR-generator which bounds to used in rounding? */
109 #define DEFAULT_USECUTPOOL TRUE /**< Use cutpool to store CG-cuts even if the are not efficient? */
110 #define DEFAULT_PRIMALSEPARATION TRUE /**< Only separate cuts that are tight for the best feasible solution? */
111 #define DEFAULT_EARLYTERM TRUE /**< Terminate separation if a violated (but possibly sub-optimal) cut has been found? */
112 #define DEFAULT_ADDVIOLATIONCONS FALSE /**< Add constraint to subscip that only allows violated cuts (otherwise add obj. limit)?*/
113 #define DEFAULT_ADDVIOLCONSHDLR FALSE /**< Add constraint handler to filter out violated cuts? */
114 #define DEFAULT_CONSHDLRUSENORM TRUE /**< Should the violation constraint handler use the norm of a cut to check for feasibility? */
115 #define DEFAULT_USEOBJUB FALSE /**< Use upper bound on objective function (via primal solution)? */
116 #define DEFAULT_USEOBJLB FALSE /**< Use lower bound on objective function (via lower bound)? */
117 #define DEFAULT_SUBSCIPFAST TRUE /**< Should the settings for the sub-MIP be optimized for speed? */
118 #define DEFAULT_OUTPUT FALSE /**< Should information about the sub-MIP and cuts be displayed? */
119 
120 #define NROWSTOOSMALL 5 /**< only separate if the number of rows is larger than this number */
121 #define NCOLSTOOSMALL 5 /**< only separate if the number of columns is larger than this number */
122 
123 #define EPSILONVALUE 1e-03 /**< epsilon value needed to model strict-inequalities */
124 #define BETAEPSILONVALUE 1e-02 /**< epsilon value for fracbeta - is larger than EPSILONVALUE for numerical stability */
125 #define STALLNODELIMIT 1000LL /**< number of stalling nodes if earlyterm is true */
126 #define CONSHDLRFULLNORM FALSE /**< compute real cut and compute norm for this (if addviolconshdlr and conshdlrusenorm are true) */
127 #define MINEFFICACY 0.05 /**< minimum efficacy of a cut - compare set.c */
128 #define MAXNSOLS 1000 /**< maximal number of solutions stored in sub-SCIP */
129 #define OBJWEIGHTRANGE 0.01 /**< maximal range of scaling of objective w.r.t. size of rows */
130 
131 /* parameters used for CMIR-generation (taken from sepa_gomory) */
132 #define BOUNDSWITCH 0.9999
133 #define USEVBDS TRUE
134 #define MINFRAC 0.0009 /* to allow a deviation of the same size as EPSILONVALUE */
135 #define MAXFRAC 0.9991 /* to allow a deviation of the same size as EPSILONVALUE */
136 #define FIXINTEGRALRHS FALSE
137 #define MAKECONTINTEGRAL FALSE
138 #define MAXWEIGHTRANGE 1e+05 /**< maximal valid range max(|weights|)/min(|weights|) of row weights */
139 
140 #if 0
141 #define MAXAGGRLEN(nvars) (0.1*(nvars)+1000) /**< maximal length of base inequality */
142 #endif
143 #define MAXAGGRLEN(nvars) nvars /**< currently very large to allow any generation */
144 
145 
146 /** separator data */
147 struct SCIP_SepaData
148 {
149  int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */
150  int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */
151  int maxdepth; /**< maximal depth at which the separator is applied */
152  SCIP_Bool decisiontree; /**< Use decision tree to turn separation on/off? */
153  SCIP_Real timelimit; /**< time limit for subscip */
154  SCIP_Real memorylimit; /**< memory limit for subscip */
155  SCIP_Longint minnodelimit; /**< minimum number of nodes considered for sub-MIP (-1: unlimited) */
156  SCIP_Longint maxnodelimit; /**< maximum number of nodes considered for sub-MIP (-1: unlimited) */
157  SCIP_Real cutcoefbnd; /**< bounds on the values of the coefficients in the CG-cut */
158  SCIP_Bool onlyactiverows; /**< Use only active rows to generate cuts? */
159  int maxrowage; /**< maximal age of rows to consider if onlyactiverows is false */
160  SCIP_Bool onlyrankone; /**< Separate only rank 1 inequalities w.r.t. CG-MIP separator? */
161  SCIP_Bool onlyintvars; /**< Generate cuts for problems with only integer variables? */
162  SCIP_Bool allowlocal; /**< Allow local cuts? */
163  SCIP_Bool contconvert; /**< Convert some integral variables to be continuous to reduce the size of the sub-MIP? */
164  SCIP_Real contconvfrac; /**< fraction of integral variables converted to be continuous (if contconvert) */
165  int contconvmin; /**< minimum number of integral variables before some are converted to be continuous */
166  SCIP_Bool intconvert; /**< Convert some integral variables attaining fractional values to have integral value? */
167  SCIP_Real intconvfrac; /**< fraction of frac. integral variables converted to have integral value (if intconvert) */
168  int intconvmin; /**< minimum number of integral variables before some are converted to have integral value */
169  SCIP_Bool skipmultbounds; /**< Skip the upper bounds on the multipliers in the sub-MIP? */
170  SCIP_Bool objlone; /**< Should the objective of the sub-MIP only minimize the l1-norm of the multipliers? */
171  SCIP_Real objweight; /**< objective weight for artificial variables */
172  SCIP_Bool objweightsize; /**< Weight each row by its size? */
173  SCIP_Bool dynamiccuts; /**< Should generated cuts be removed from the LP if they are no longer tight? */
174  SCIP_Bool usecmir; /**< Use CMIR-generator (otherwise add cut directly)? */
175  SCIP_Bool usestrongcg; /**< Use strong CG-function to strengthen cut? */
176  SCIP_Bool cmirownbounds; /**< Tell CMIR-generator which bounds to used in rounding? */
177  SCIP_Bool usecutpool; /**< Use cutpool to store CG-cuts even if the are not efficient? */
178  SCIP_Bool primalseparation; /**< Only separate cuts that are tight for the best feasible solution? */
179  SCIP_Bool earlyterm; /**< Terminate separation if a violated (but possibly sub-optimal) cut has been found? */
180  SCIP_Bool addviolationcons; /**< Add constraint to subscip that only allows violated cuts? */
181  SCIP_Bool addviolconshdlr; /**< Add constraint handler to filter out violated cuts? */
182  SCIP_Bool conshdlrusenorm; /**< Should the violation constraint handler use the cut-norm to check for feasibility? */
183  SCIP_Bool useobjub; /**< Use upper bound on objective function (via primal solution)? */
184  SCIP_Bool useobjlb; /**< Use lower bound on objective function (via lower bound)? */
185  SCIP_Bool subscipfast; /**< Should the settings for the sub-MIP be optimized for speed? */
186  SCIP_Bool output; /**< Should information about the sub-MIP and cuts be displayed? */
187 };
188 
189 
190 /** what happens for columns in the LP */
192 {
193  colPresent = 0, /**< column is present in the separating MIP */
194  colContinuous = 1, /**< column corresponds to a continuous variable */
195  colConverted = 2, /**< column is converted to be continuous */
196  colAtUb = 3, /**< variable corresponding to column was at it's upper bound and was complemented */
197  colAtLb = 4 /**< variable corresponding to column was at it's lower bound (possibly complemented) */
198 };
200 
201 
202 /** data for the sub-MIP */
203 struct CGMIP_MIPData
204 {
205  SCIP* subscip; /**< pointer to (sub)SCIP data structure containing the auxiliary IP */
206  unsigned int m; /**< number of constraints of subscip */
207  unsigned int n; /**< number of variables of subscip */
208  unsigned int nrows; /**< number of rows of original LP */
209  unsigned int ncols; /**< number of columns of original LP */
210  unsigned int ntotalrows; /**< number of total rows used (possibly including objective rows) */
211 
212  SCIP_VAR** alpha; /**< cut coefficient variable (NULL if not in separating MIP) */
213  SCIP_VAR* beta; /**< rhs of cut */
214  SCIP_VAR** fracalpha; /**< fractional part of lhs of cut (NULL if not present) */
215  SCIP_VAR* fracbeta; /**< fractional part of rhs of cut */
216  CGMIP_COLTYPE* coltype; /**< type for the columns */
217  SCIP_Bool* iscomplemented; /**< whether the variable was complemented */
218  SCIP_Bool* isshifted; /**< whether the variable was shifted to have 0 lower bound */
219 
220  SCIP_VAR** ylhs; /**< auxiliary row variables for lhs (NULL if not present) */
221  SCIP_VAR** yrhs; /**< auxiliary row variables for rhs (NULL if not present) */
222 
223  SCIP_VAR** z; /**< auxiliary variables for upper bounds (NULL if not present) */
224 
225  char normtype; /**< type of norm to use for efficacy norm calculation */
226 
227  /* additional redundant data */
228  SCIP_Bool conshdlrusenorm; /**< copy from sepadata */
229  SCIP_Bool conshdlrfullnorm; /**< compute real cut and compute norm for this (if addviolconshdlr and conshdlrusenorm are true) */
230  SCIP* scip; /**< original SCIP */
231  SCIP_SEPA* sepa; /**< CG-cut separator */
232  SCIP_SEPADATA* sepadata; /**< CG-cut separator data */
233 };
234 typedef struct CGMIP_MIPData CGMIP_MIPDATA;
235 
236 
237 /*
238  * constraint handler to filter out violated cuts
239  */
240 
241 /* constraint handler properties */
242 #define CONSHDLR_NAME "violatedCuts"
243 #define CONSHDLR_DESC "only allow solutions corresponding to violated cuts"
244 
245 /** constraint handler data */
246 struct SCIP_ConshdlrData
247 {
248  CGMIP_MIPDATA* mipdata; /**< data of separating sub-MIP */
249 };
250 
251 /* temporary forward declaration */
252 static
254  SCIP* scip, /**< original scip */
255  SCIP_SEPA* sepa, /**< separator */
256  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
257  SCIP_SEPADATA* sepadata, /**< separator data */
258  SCIP_SOL* sol, /**< current solution for sub-MIP */
259  SCIP_Real* cutcoefs, /**< coefficients of the cut */
260  SCIP_Real* cutrhs, /**< rhs of the cut */
261  SCIP_Bool* localrowsused, /**< pointer to store whether local rows were used in summation */
262  SCIP_Bool* localboundsused, /**< pointer to store whether local bounds were used in summation */
263  int * cutrank, /**< pointer to store the cut rank */
264  SCIP_Bool* success /**< whether we produced a valid cut */
265  );
266 
267 
268 /** check whether cut corresponding to solution is violated */
269 static
271  SCIP* scip, /**< SCIP data structure */
272  CGMIP_MIPDATA* mipdata, /**< data of separating sub-MIP */
273  SCIP_SOL* sol, /**< solution to be checked */
274  SCIP_Bool* violated /**< pointer to store if the cut is violated */
275  )
276 {
277  SCIP_Real cutsqrnorm = 0.0;
278  SCIP* subscip;
279  SCIP_Real act;
280  SCIP_Real norm;
281  SCIP_Real val;
282  SCIP_VAR* var;
283  SCIP_Real rhs;
284  unsigned int j;
285  int len = 0;
286 
287  assert( mipdata != NULL );
288  subscip = mipdata->subscip;
289  assert( subscip != NULL );
290  assert( violated != NULL );
291 
292  /* initialize activity and norm */
293  act = 0.0;
294  norm = 1.0;
295  *violated = FALSE;
296 
297  /* compute activity and norm */
298  if ( mipdata->conshdlrusenorm )
299  {
300  /* check whether we should compute the full cut and then compute the norm */
301  if ( mipdata->conshdlrfullnorm )
302  {
303  SCIP_Real* cutcoefs;
304  SCIP_Bool localrowsused;
305  SCIP_Bool localboundsused;
306  SCIP_Bool success;
307  SCIP_VAR** vars;
308  int cutrank = 0;
309  int nvars;
310 
311  /* get data */
312  SCIP_CALL( SCIPgetVarsData(mipdata->scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
313  assert(nvars >= 0);
314  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
315 
316  /* compute coefficients */
317  SCIP_CALL( computeCut(mipdata->scip, mipdata->sepa, mipdata, mipdata->sepadata, sol, cutcoefs, &rhs, &localrowsused, &localboundsused, &cutrank, &success) );
318 
319 #if 0
320  for (j = 0; j < (unsigned int) nvars; ++j)
321  {
322  if ( ! SCIPisZero(scip, cutcoefs[j]) )
323  SCIPinfoMessage(scip, NULL, "+ %f x%d", cutcoefs[j], j);
324  }
325  SCIPinfoMessage(scip, NULL, "\n");
326 #endif
327 
328  /* ignore solution if cut was not valid */
329  if ( ! success )
330  return SCIP_OKAY;
331 
332  /* compute activity and Euclidean norm (todo: use arbitrary norm) */
333  cutsqrnorm = 0.0;
334  for (j = 0; j < (unsigned int) nvars; ++j)
335  {
336  if ( ! SCIPisZero(scip, cutcoefs[j]) )
337  {
338  act += cutcoefs[j] * SCIPvarGetLPSol(vars[j]);
339  cutsqrnorm += SQR(cutcoefs[j]);
340  }
341  }
342  norm = SQRT(cutsqrnorm);
343 
344  SCIPfreeBufferArray(scip, &cutcoefs);
345  }
346  else
347  {
348  switch ( mipdata->normtype )
349  {
350  case 'e':
351  cutsqrnorm = 0.0;
352  for (j = 0; j < mipdata->ncols; ++j)
353  {
354  var = mipdata->alpha[j];
355  if ( var == NULL )
356  continue;
357 
358  val = SCIPgetSolVal(subscip, sol, var);
359  if ( !SCIPisZero(scip, val) )
360  {
361  act += val * SCIPvarGetObj(var);
362  cutsqrnorm += SQR(val);
363  }
364  }
365  norm = SQRT(cutsqrnorm);
366  break;
367  case 'm':
368  for (j = 0; j < mipdata->ncols; ++j)
369  {
370  var = mipdata->alpha[j];
371  if ( var == NULL )
372  continue;
373 
374  val = SCIPgetSolVal(subscip, sol, var);
375  if ( !SCIPisZero(scip, val) )
376  {
377  act += val * SCIPvarGetObj(var);
378  if ( REALABS(val) > norm )
379  norm = REALABS(val);
380  }
381  }
382  break;
383  case 's':
384  for (j = 0; j < mipdata->ncols; ++j)
385  {
386  var = mipdata->alpha[j];
387  if ( var == NULL )
388  continue;
389 
390  val = SCIPgetSolVal(subscip, sol, var);
391  if ( !SCIPisZero(scip, val) )
392  {
393  act += val * SCIPvarGetObj(var);
394  norm += REALABS(val);
395  }
396  }
397  break;
398  case 'd':
399  for (j = 0; j < mipdata->ncols; ++j)
400  {
401  var = mipdata->alpha[j];
402  if ( var == NULL )
403  continue;
404 
405  val = SCIPgetSolVal(subscip, sol, var);
406  if ( !SCIPisZero(scip, val) )
407  {
408  act += val * SCIPvarGetObj(var);
409  ++len;
410  }
411  }
412  if ( len > 0 )
413  norm = 1.0;
414  break;
415  default:
416  SCIPerrorMessage("invalid efficacy norm parameter '%c'\n", mipdata->normtype);
417  return SCIP_INVALIDDATA;
418  }
419  /* get rhs */
420  rhs = SCIPgetSolVal(subscip, sol, mipdata->beta);
421  }
422 
423  /* if norm is 0, the cut is trivial */
424  if ( SCIPisZero(subscip, norm) )
425  return SCIP_OKAY;
426  }
427  else
428  {
429  for (j = 0; j < mipdata->ncols; ++j)
430  {
431  var = mipdata->alpha[j];
432  if ( var == NULL )
433  continue;
434 
435  val = SCIPgetSolVal(subscip, sol, var);
436  if ( !SCIPisZero(subscip, val) )
437  act += SCIPvarGetObj(var) * val;
438  }
439 
440  /* get rhs */
441  rhs = SCIPgetSolVal(subscip, sol, mipdata->beta);
442  }
443 
444 #ifdef SCIP_DEBUG
445  if ( SCIPisEfficacious(subscip, (act - rhs)/norm) )
446  {
447  SCIPdebugMessage("Violated cut from solution - act: %f, rhs: %f, norm: %f, eff.: %f\n", act, rhs, norm, (act-rhs)/norm);
448  }
449  else
450  {
451  SCIPdebugMessage("Rejected cut from solution - act: %f, rhs: %f, norm: %f, eff.: %f\n", act, rhs, norm, (act-rhs)/norm);
452  }
453 #endif
454 
455  *violated = SCIPisEfficacious(subscip, (act - rhs)/norm);
456 
457  return SCIP_OKAY;
458 }
459 
460 
461 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
462 static
463 SCIP_DECL_CONSFREE(consFreeViolatedCuts)
464 { /*lint --e{715}*/
465  SCIP_CONSHDLRDATA* conshdlrdata;
466 
467  assert( scip != NULL );
468  assert( conshdlr != NULL );
469  conshdlrdata = SCIPconshdlrGetData(conshdlr);
470  assert( conshdlrdata != NULL );
471 
472  SCIPfreeMemory(scip, &conshdlrdata);
473 
474  return SCIP_OKAY;
475 }
476 
477 
478 /** constraint enforcing method of constraint handler for LP solutions */
479 static
480 SCIP_DECL_CONSENFOLP(consEnfolpViolatedCuts)
481 { /*lint --e{715}*/
482  SCIP_CONSHDLRDATA* conshdlrdata;
483  SCIP_Bool violated;
484 
485  assert( scip != NULL );
486  assert( conshdlr != NULL );
487  assert( result != NULL );
488 
489  assert( SCIPgetNLPBranchCands(scip) == 0 );
490 
491  conshdlrdata = SCIPconshdlrGetData(conshdlr);
492  assert( conshdlrdata != NULL );
493 
494  SCIP_CALL( solCutIsViolated(scip, conshdlrdata->mipdata, NULL, &violated) );
495 
496  if ( violated )
497  *result = SCIP_FEASIBLE;
498  else
499  *result = SCIP_CUTOFF; /* cutoff, since all integer variables are integer, but the solution is not feasible */
500 
501  return SCIP_OKAY;
502 }
503 
504 
505 /** constraint enforcing method of constraint handler for pseudo solutions */
506 static
507 SCIP_DECL_CONSENFOPS(consEnfopsViolatedCuts)
508 { /*lint --e{715}*/
509  assert( result != NULL );
510 
511  /* this function should better not be called, since we need an LP solution for the sub-MIP to
512  * make sense, because of the multiplier variables. We therefore return SCIP_FEASIBLE. */
513  *result = SCIP_FEASIBLE;
514 
515  return SCIP_OKAY;
516 }
517 
518 
519 /** feasibility check method of constraint handler for integral solutions */
520 static
521 SCIP_DECL_CONSCHECK(consCheckViolatedCuts)
522 { /*lint --e{715}*/
523  SCIP_CONSHDLRDATA* conshdlrdata;
524  SCIP_Bool violated;
525 
526  assert( scip != NULL );
527  assert( conshdlr != NULL );
528  assert( sol != NULL );
529  assert( result != NULL );
530 
531  conshdlrdata = SCIPconshdlrGetData(conshdlr);
532  assert( conshdlrdata != NULL );
533 
534  SCIP_CALL( solCutIsViolated(scip, conshdlrdata->mipdata, sol, &violated) );
535 
536  if ( violated )
537  *result = SCIP_FEASIBLE;
538  else
539  *result = SCIP_INFEASIBLE;
540 
541  return SCIP_OKAY;
542 }
543 
544 
545 /** variable rounding lock method of constraint handler */
546 static
547 SCIP_DECL_CONSLOCK(consLockViolatedCuts)
548 { /*lint --e{715}*/
549  /* do not lock variables */
550  return SCIP_OKAY;
551 }
552 
553 
554 /** creates the violated CG-cut constraint handler and includes it in SCIP */
555 static
557  SCIP* scip, /**< SCIP data structure */
558  CGMIP_MIPDATA* mipdata /**< data of separating sub-MIP */
559  )
560 {
561  SCIP_CONSHDLRDATA* conshdlrdata;
562  SCIP_CONSHDLR* conshdlr;
563 
564  SCIP_CALL( SCIPallocMemory(scip, &conshdlrdata) );
565  conshdlrdata->mipdata = mipdata;
566 
567  /* include constraint handler */
569  -1000000, -1000000, 100, FALSE,
570  consEnfolpViolatedCuts, consEnfopsViolatedCuts, consCheckViolatedCuts, consLockViolatedCuts,
571  conshdlrdata) );
572 
573  assert(conshdlr != NULL);
574 
575  /* set non-fundamental callbacks via specific setter functions */
576  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeViolatedCuts) );
577 
578  return SCIP_OKAY;
579 }
580 
581 
582 /*
583  * local methods
584  */
585 
586 
587 /** stores nonzero elements of dense coefficient vector as sparse vector and calculates activity and norm
588  *
589  * copied from sepa_gomory.c
590  */
591 static
593  SCIP* scip, /**< SCIP data structure */
594  int nvars, /**< number of problem variables */
595  SCIP_VAR** vars, /**< problem variables */
596  SCIP_Real* cutcoefs, /**< dense coefficient vector */
597  SCIP_Real* varsolvals, /**< dense variable LP solution vector */
598  char normtype, /**< type of norm to use for efficacy norm calculation */
599  SCIP_VAR** cutvars, /**< array to store variables of sparse cut vector */
600  SCIP_Real* cutvals, /**< array to store coefficients of sparse cut vector */
601  int* cutlen, /**< pointer to store number of nonzero entries in cut */
602  SCIP_Real* cutact, /**< pointer to store activity of cut */
603  SCIP_Real* cutnorm /**< pointer to store norm of cut vector */
604  )
605 {
606  SCIP_Real val;
607  SCIP_Real cutsqrnorm;
608  SCIP_Real act;
609  SCIP_Real norm;
610  int len;
611  int v;
612 
613  assert( nvars == 0 || cutcoefs != NULL );
614  assert( nvars == 0 || varsolvals != NULL );
615  assert( cutvars != NULL );
616  assert( cutvals != NULL );
617  assert( cutlen != NULL );
618  assert( cutact != NULL );
619  assert( cutnorm != NULL );
620 
621  len = 0;
622  act = 0.0;
623  norm = 0.0;
624  switch ( normtype )
625  {
626  case 'e':
627  cutsqrnorm = 0.0;
628  for (v = 0; v < nvars; ++v)
629  {
630  val = cutcoefs[v];
631  if ( !SCIPisZero(scip, val) )
632  {
633  act += val * varsolvals[v];
634  cutsqrnorm += SQR(val);
635  cutvars[len] = vars[v];
636  cutvals[len++] = val;
637  }
638  }
639  norm = SQRT(cutsqrnorm);
640  break;
641  case 'm':
642  for (v = 0; v < nvars; ++v)
643  {
644  val = cutcoefs[v];
645  if ( !SCIPisZero(scip, val) )
646  {
647  act += val * varsolvals[v];
648  if ( REALABS(val) > norm )
649  norm = REALABS(val);
650  cutvars[len] = vars[v];
651  cutvals[len++] = val;
652  }
653  }
654  break;
655  case 's':
656  for (v = 0; v < nvars; ++v)
657  {
658  val = cutcoefs[v];
659  if ( !SCIPisZero(scip, val) )
660  {
661  act += val * varsolvals[v];
662  norm += REALABS(val);
663  cutvars[len] = vars[v];
664  cutvals[len++] = val;
665  }
666  }
667  break;
668  case 'd':
669  for (v = 0; v < nvars; ++v)
670  {
671  val = cutcoefs[v];
672  if ( !SCIPisZero(scip, val) )
673  {
674  act += val * varsolvals[v];
675  cutvars[len] = vars[v];
676  cutvals[len++] = val;
677  }
678  }
679  if ( len > 0 )
680  norm = 1.0;
681  break;
682  default:
683  SCIPerrorMessage("invalid efficacy norm parameter '%c'\n", normtype);
684  return SCIP_INVALIDDATA;
685  }
686 
687  *cutlen = len;
688  *cutact = act;
689  *cutnorm = norm;
690 
691  return SCIP_OKAY;
692 }
693 
694 
695 /** Compute lhs/rhs for transformed column
696  *
697  * Consider a variable \f$x_j\f$ and some row of the original system:
698  * \f[
699  * \gamma \leq a^T x \leq \delta, \quad \ell_j \leq x_j \leq u_j.
700  * \f]
701  * We perform the transformation
702  * \f[
703  * x_i' = \left\{
704  * \begin{array}{ll}
705  * s + \frac{1}{\sigma}\, x_j & \mbox{if }i = j\\
706  * x_i & \mbox{otherwise},
707  * \end{array}
708  * \right.
709  * \f]
710  * where \f$s\f$ is the offset value and \f$\sigma\f$ is a scaling factor. The new system is
711  * \f[
712  * \gamma + \sigma\, a_j\,s \leq \sum_{i \neq j} a_i\, x_i' + \sigma a_j\, x_j' \leq \delta + \sigma\, a_j\, s
713  * \f]
714  * with bounds
715  * \f[
716  * \frac{1}{\sigma} \ell_j + s \leq x_j' \leq \frac{1}{\sigma} u_j + s, \qquad \mbox{ if }\sigma > 0
717  * \f]
718  * and
719  * \f[
720  * \frac{1}{\sigma} u_j + s \leq x_j' \leq \frac{1}{\sigma} \ell_j + s, \qquad \mbox{ if }\sigma < 0.
721  * \f]
722  *
723  * This can be used as follows:
724  *
725  * - If \f$x_j \geq \ell_j\f$ has a (nonzero) lower bound, one can use \f$s = -\ell_j\f$, \f$\sigma = 1\f$,
726  * and obtain \f$\gamma - a_j\,\ell_j \leq a^T x' \leq \delta - a_j\,\ell_j\f$, \f$0 \leq x_j' \leq u_j - \ell_j\f$.
727  *
728  * - If \f$x_j \leq u_j\f$ has a (nonzero) upper bound, one can use \f$s = u_j\f$, \f$\sigma = -1\f$,
729  * and obtain \f$\gamma - a_j\,u_j \leq \sum_{i \neq j} a_i\, x_i' - a_j\, x_j' \leq \delta - a_j\, u_j\f$,
730  * \f$0 \leq x_j' \leq u_j - \ell_j\f$.
731  */
732 static
734  SCIP* scip, /**< SCIP data structure */
735  SCIP_SEPADATA* sepadata, /**< separator data */
736  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
737  SCIP_COL* col, /**< column that should be complemented */
738  SCIP_Real offset, /**< offset by which column should be shifted */
739  SCIP_Real sigma, /**< scaling factor */
740  SCIP_Real* lhs, /**< array of lhs of rows */
741  SCIP_Real* rhs, /**< array rhs of rows */
742  SCIP_Real* lb, /**< pointer to lb of column */
743  SCIP_Real* ub, /**< pointer to ub of column */
744  SCIP_Real* primsol /**< pointer to solution value */
745  )
746 {
747  SCIP_ROW** colrows;
748  SCIP_Real* colvals;
749  int pos, i;
750 
751  assert( scip != NULL );
752  assert( lhs != NULL );
753  assert( rhs != NULL );
754  assert( col != NULL );
755 
756  colrows = SCIPcolGetRows(col);
757  colvals = SCIPcolGetVals(col);
758  assert( SCIPcolGetNLPNonz(col) == 0 || colrows != NULL );
759  assert( SCIPcolGetNLPNonz(col) == 0 || colvals != NULL );
760  assert( ! SCIPisZero(scip, sigma) );
761 
762  /* loop through rows that contain column */
763  for (i = 0; i < SCIPcolGetNLPNonz(col); ++i)
764  {
765  SCIP_ROW* row;
766 
767  row = colrows[i];
768  assert( row != NULL );
769 
770  /* skip modifiable rows and local rows, unless allowed */
771  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
772  continue;
773 
774  pos = SCIProwGetLPPos(row);
775  assert( 0 <= pos && pos < (int) mipdata->nrows );
776 
777  assert( ! SCIPisInfinity(scip, lhs[pos]) );
778  if ( ! SCIPisInfinity(scip, -lhs[pos]) )
779  lhs[pos] += sigma * colvals[i] * offset;
780 
781  assert( ! SCIPisInfinity(scip, -rhs[pos]) );
782  if ( ! SCIPisInfinity(scip, rhs[pos]) )
783  rhs[pos] += sigma * colvals[i] * offset;
784  }
785 
786  /* check objective function */
787  if ( sepadata->useobjub || sepadata->useobjlb )
788  {
789  assert( SCIPisEQ(scip, SCIPcolGetObj(col), SCIPvarGetObj(SCIPcolGetVar(col))) );
790  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
791 
792  if ( ! SCIPisInfinity(scip, -lhs[mipdata->nrows]) )
793  lhs[mipdata->nrows] += sigma * SCIPcolGetObj(col) * offset;
794 
795  if ( ! SCIPisInfinity(scip, rhs[mipdata->nrows]) )
796  rhs[mipdata->nrows] += sigma * SCIPcolGetObj(col) * offset;
797  }
798 
799  /* correct lower and upper bounds and solution */
800  if ( SCIPisNegative(scip, sigma) )
801  {
802  SCIP_Real l;
803 
804  assert( ! SCIPisInfinity(scip, -*ub) );
805  if ( ! SCIPisInfinity(scip, *ub) )
806  l = *ub/sigma + offset;
807  else
808  l = -SCIPinfinity(scip);
809 
810  assert( ! SCIPisInfinity(scip, *lb) );
811  if ( ! SCIPisInfinity(scip, -*lb) )
812  *ub = *lb/sigma + offset;
813  else
814  *ub = SCIPinfinity(scip);
815  *lb = l;
816  }
817  else
818  {
819  assert( ! SCIPisInfinity(scip, *lb) );
820  if ( ! SCIPisInfinity(scip, -*lb) )
821  *lb = *lb/sigma + offset;
822  assert( ! SCIPisInfinity(scip, -*ub) );
823  if ( ! SCIPisInfinity(scip, *ub) )
824  *ub = *ub/sigma + offset;
825  }
826  *primsol = *primsol/sigma + offset;
827 
828  return SCIP_OKAY;
829 }
830 
831 
832 /** compute objective coefficient for rows that are weighted by size
833  *
834  * The objective is computed by multiplying a default value by
835  * \f[
836  * 1 - (r_{\mbox{max}} - r) \frac{1 - a}{r_{\mbox{max}} - r_{\mbox{min}}},
837  * \f]
838  * where \f$r\f$ is the size of the current row, \f$a \in [0,1]\f$ is a parameter, and \f$r_{\mbox{max}}\f$ and
839  * \f$r_{\mbox{min}}\f$ are the maximal and minimal size of a row, respectively.
840  *
841  * Thus, if \f$r = r_{\mbox{max}}\f$, we get 1 and if \f$r = r_{\mbox{min}}\f$, we get \f$a\f$.
842  */
843 static
845  int rowsize, /**< size of current row */
846  int minrowsize, /**< maximal size of rows */
847  int maxrowsize /**< minimal size of rows */
848  )
849 {
850  SCIP_Real a;
851 
852  assert( maxrowsize > 0 );
853  assert( minrowsize < INT_MAX );
854  assert( minrowsize <= maxrowsize );
855  assert( minrowsize <= rowsize && rowsize <= maxrowsize );
856 
857  if ( minrowsize == maxrowsize )
858  return 1.0;
859 
860  a = (1.0 - OBJWEIGHTRANGE)/((SCIP_Real) (maxrowsize - minrowsize));
861 
862  return 1.0 - a * ((SCIP_Real) (maxrowsize - rowsize));
863 }
864 
865 
866 /** Creates a subscip representing the separating MIP.
867  *
868  * Let the constraints of the original MIP be of the following form:
869  * \f[
870  * \begin{array}{l@{\;}ll}
871  * a \leq A x + & C r & \leq b\\
872  * \ell \leq x & & \leq u\\
873  * c \leq & r & \leq d\\
874  * x \in Z^n.
875  * \end{array}
876  * \f]
877  * Here, some of the bounds may have value \f$\infty\f$ or \f$-\infty\f$. Written in
878  * \f$\leq\f$-form this becomes:
879  * \f[
880  * \begin{array}{r@{\;}l}
881  * \tilde{A} x + \tilde{C} r & \leq \tilde{b}\\
882  * -x & \leq -\ell\\
883  * x & \leq u\\
884  * -r & \leq -c\\
885  * r & \leq d\\
886  * x \in Z^n,
887  * \end{array}
888  * \f]
889  * where we use
890  * \f[
891  * \tilde{A} =
892  * \left[
893  * \begin{array}{r}
894  * -A \\
895  * A
896  * \end{array}
897  * \right],
898  * \quad
899  * \tilde{C} =
900  * \left[
901  * \begin{array}{r}
902  * - C\\
903  * C
904  * \end{array}
905  * \right]
906  * \qquad\mbox{ and }\qquad
907  * \tilde{b} =
908  * \left[
909  * \begin{array}{r}
910  * -a\\
911  * b
912  * \end{array}
913  * \right].
914  * \f]
915  * For the moment we assume that \f$c = 0\f$, i.e., the lower bounds on the continuous variables
916  * are 0. To obtain a Chv&aacute;tal-Gomory cut we have to find nonnegative multipliers \f$y\f$,
917  * \f$\underline{z}\f$, and \f$\overline{z}\f$ such that
918  * \f[
919  * y^T \tilde{A} - \underline{z}^T + \overline{z}^T \in Z \qquad\mbox{ and }\qquad
920  * y^T \tilde{C} \geq 0.
921  * \f]
922  * Note that we use zero multipliers for the bounds on the continuous variables \f$r\f$. Moreover,
923  * if some bounds are infinity, the corresponding multipliers are assumed to be 0. From these
924  * conditions, we obtain
925  * \f[
926  * (y^T \tilde{A} - \underline{z}^T + \overline{z}^T)\, x +
927  * y^T \tilde{C} \, r \leq
928  * y^T \tilde{b} - \underline{z}^T \ell + \overline{z}^T u.
929  * \f]
930  * Because \f$r \geq 0\f$, we can ignore the term \f$y^T \tilde{C} \, r \geq 0\f$ and obtain the
931  * following cut:
932  * \f[
933  * (y^T \tilde{A} - \underline{z}^T + \overline{z}^T )\, x \leq
934  * \lfloor y^T \tilde{b} - \underline{z}^T \ell + \overline{z}^T u \rfloor.
935  * \f]
936  * Assume that \f$\ell = 0\f$ for the meantime. Then the cut can be written as:
937  * \f[
938  * \lfloor y^T \tilde{A} + \overline{z}^T \rfloor \, x \leq
939  * \lfloor y^T \tilde{b} + \overline{z}^T u \rfloor.
940  * \f]
941  *
942  * Following Fischetti and Lodi [2005], let \f$(x^*,r^*)\f$ be a fractional solution of the above
943  * original system. The separating MIP created below is
944  * \f[
945  * \begin{array}{rlr@{\;}l}
946  * \max & \multicolumn{2}{@{}l}{(x^*)^T \alpha - \beta - w^T y} &\\
947  * & f = & \tilde{A}^T y + \overline{z} - \alpha & \\
948  * & \tilde{f} = & \tilde{b}^T y + u^T \overline{z} - \beta &\\
949  * & & \tilde{C}^T y & \geq 0\\
950  * & & 0 \leq f & \leq 1 - \epsilon \\
951  * & & 0 \leq \tilde{f} & \leq 1 - \epsilon\\
952  * & & 0 \leq y, \overline{z} & \leq 1 - \epsilon.\\
953  * & & \alpha \in Z^m, \beta & \in Z.
954  * \end{array}
955  * \f]
956  * Here, \f$w\f$ is a weight vector; it's idea is to make the sum over all components of \f$y\f$ as
957  * small as possible, in order to generate sparse cuts.
958  *
959  * We perform the following additional computations:
960  *
961  * - If the lower bounds on \f$x_i\f$ or \f$r_j\f$ are finite, we shift the variable to have a zero
962  * lower bound, i.e., we replace it by \f$x_i - \ell_i\f$ (or \f$r_j - u_j\f$). This is helpful in
963  * several ways: As seen above, the resulting inequalities/formulations simplify. Moreover, it
964  * allows to drop a variable if \f$x^*_i = 0\f$, see the next comment. If the lower bounds are not
965  * finite, but the upper bounds are finite, we can complement the variable. If the variables are
966  * free, the above formulation changes as follows: For free continuous variables, we require
967  * \f$\tilde{C}^T y = 0\f$. For a free integer variable \f$x_j\f$ (which rarely occurs in
968  * practice), we require \f$f_j = 0\f$, i.e., we force that \f$(\tilde{A}^T y + \overline{z})_j =
969  * \alpha_j\f$.
970  *
971  * - If \f$x^*_j = 0 = \ell_j\f$ (after the above preprocessing), we drop variable \f$\alpha_j\f$
972  * from the formulation. Let \f$(\alpha^*, \beta^*, y^*, \overline{z}^*)\f$ be an
973  * optimal solution to the separating MIP. Then we can compute \f$\alpha_j =
974  * \lfloor(\tilde{A}_j^T y^* + \overline{z}^*)\rfloor\f$.
975  *
976  * - If \f$x^*_i = u_i\f$, we complement the variable and drop it from the formulation, since the
977  * lower bound is 0 afterwards.
978  *
979  * - If a variable has been shifted or complemented, we have to recompute \f$\beta\f$ with the
980  * original lhs/rhs.
981  *
982  * - If a continuous variable \f$r_j\f$ is free, we have to force equality for the corresponding components in
983  * \f$y^T \tilde{C} \, r \geq 0\f$.
984  *
985  * - If an integer variable \f$x_i\f$ is free, we are not allowed to round the cut down. In this
986  * case, the combintation of rows and bounds has to be integral. We force this by requiring that
987  * \f$f_i = 0\f$.
988  *
989  * - If @p contconvert is true some integral variables are randomly treated as if they were
990  * continuous. This has the effect that in the resulting cut the corresponding coefficient has
991  * value 0. This makes the cuts more sparse. Moreover, the separation problems should become
992  * easier.
993  *
994  * - If required, i.e., parameter @p primalseparation is true, we force a primal separation step. For
995  * this we require that the cut is tight at the currently best solution. To get reliable solutions
996  * we relax equality by EPSILONVALUE.
997  *
998  * - If required (via parameters @p useobjub or @p useobjlb), we add a row corresponding to the objective function with
999  * respect to the current lower and upper bounds.
1000  */
1001 static
1003  SCIP* scip, /**< SCIP data structure */
1004  SCIP_SEPA* sepa, /**< separator */
1005  SCIP_SEPADATA* sepadata, /**< separator data */
1006  CGMIP_MIPDATA* mipdata /**< data for sub-MIP */
1007  )
1008 {
1009  SCIP* subscip;
1010  SCIP_COL** cols;
1011  SCIP_ROW** rows;
1012  SCIP_Real* lhs;
1013  SCIP_Real* rhs;
1014  SCIP_Real* lb;
1015  SCIP_Real* ub;
1016  SCIP_Real* primsol;
1017  SCIP_Real multvarub;
1018 
1019  unsigned int cnt;
1020  unsigned int ucnt;
1021  unsigned int nshifted;
1022  unsigned int ncomplemented;
1023  unsigned int ncontconverted;
1024  unsigned int nintconverted;
1025  unsigned int nlbounds;
1026  unsigned int nubounds;
1027 
1028  SCIP_VAR** consvars;
1029  SCIP_Real* consvals;
1030  SCIP_CONS* cons;
1031  int nconsvars;
1032  char name[SCIP_MAXSTRLEN];
1033 
1034  int ncols;
1035  int nrows;
1036  int ntotalrows;
1037  int maxrowsize = 0;
1038  int minrowsize = INT_MAX;
1039  int i, j;
1040 
1041  assert( scip != NULL );
1042  assert( sepadata != NULL );
1043 
1044  assert( mipdata->subscip == NULL );
1045 
1046  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
1047  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
1048  assert( ncols > 0 && nrows > 0 );
1049 
1050  mipdata->m = 0;
1051  mipdata->n = 0;
1052  mipdata->nrows = (unsigned int) nrows;
1053  mipdata->ncols = (unsigned int) ncols;
1054  mipdata->ntotalrows = mipdata->nrows;
1055 
1056  if ( sepadata->useobjub || sepadata->useobjlb )
1057  mipdata->ntotalrows = mipdata->nrows + 1;
1058 
1059  assert(mipdata->ntotalrows <= INT_MAX);
1060  ntotalrows = (int) mipdata->ntotalrows;
1061 
1062  /* copy value */
1063  mipdata->conshdlrusenorm = sepadata->conshdlrusenorm;
1064 
1065  /* create subscip */
1066  SCIP_CALL( SCIPcreate( &(mipdata->subscip) ) );
1067  subscip = mipdata->subscip;
1069 
1070  /* add violation constraint handler if requested */
1071  if ( sepadata->addviolconshdlr )
1072  {
1073  SCIP_CALL( SCIPincludeConshdlrViolatedCut(subscip, mipdata) );
1074  }
1075 
1076  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "sepa_cgmip separating MIP (%s)", SCIPgetProbName(scip));
1077  SCIP_CALL( SCIPcreateProb(subscip, name, NULL, NULL , NULL , NULL , NULL , NULL , NULL) );
1079 
1080  /* alloc memory for subscipdata elements */
1081  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->alpha), ncols) );
1082  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->fracalpha), ncols) );
1083  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->coltype), ncols) );
1084  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->iscomplemented), ncols) );
1085  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->isshifted), ncols) );
1086  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->ylhs), ntotalrows) );
1087  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->yrhs), ntotalrows) );
1088  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(mipdata->z), 2*ncols) );
1089 
1090  /* get temporary storage */
1091  SCIP_CALL( SCIPallocBufferArray(scip, &lhs, ntotalrows) );
1092  SCIP_CALL( SCIPallocBufferArray(scip, &rhs, ntotalrows) );
1093  SCIP_CALL( SCIPallocBufferArray(scip, &lb, ncols) );
1094  SCIP_CALL( SCIPallocBufferArray(scip, &ub, ncols) );
1095  SCIP_CALL( SCIPallocBufferArray(scip, &primsol, ncols) );
1096 
1097  /* store lhs/rhs for complementing (see below) and compute maximal nonzeros of candidate rows */
1098  for (i = 0; i < nrows; ++i)
1099  {
1100  SCIP_Real val;
1101  SCIP_ROW* row;
1102 
1103  row = rows[i];
1104  assert( row != NULL );
1105 
1106  val = SCIProwGetLhs(row) - SCIProwGetConstant(row);
1107  if ( SCIProwIsIntegral(row) )
1108  val = SCIPfeasCeil(scip, val); /* row is integral: round left hand side up */
1109  lhs[i] = val;
1110 
1111  val = SCIProwGetRhs(row) - SCIProwGetConstant(row);
1112  if ( SCIProwIsIntegral(row) )
1113  val = SCIPfeasFloor(scip, val); /* row is integral: round right hand side down */
1114  rhs[i] = val;
1115 
1116  /* skip modifiable rows and local rows, unless allowed */
1117  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
1118  continue;
1119 
1120  /* skip rows that not have been active for a longer time */
1121  if ( ! sepadata->onlyactiverows && sepadata->maxrowage > 0 && SCIProwGetAge(row) > sepadata->maxrowage )
1122  continue;
1123 
1124  /* check whether we want to skip cut produced by the CGMIP separator */
1125  if ( sepadata->onlyrankone )
1126  {
1127  if ( SCIProwGetOriginSepa(row) == sepa )
1128  continue;
1129  }
1130 
1131  /* determine maximal row size: */
1132  val = SCIPgetRowLPActivity(scip, row);
1133  if ( ! SCIPisInfinity(scip, REALABS(lhs[i])) )
1134  {
1135  if ( ! sepadata->onlyactiverows || SCIPisFeasEQ(scip, val, SCIProwGetLhs(row)) )
1136  {
1137  if ( SCIProwGetNLPNonz(row) > maxrowsize )
1138  maxrowsize = SCIProwGetNLPNonz(row);
1139  if ( SCIProwGetNLPNonz(row) < minrowsize )
1140  minrowsize = SCIProwGetNLPNonz(row);
1141  }
1142  }
1143  else
1144  {
1145  if ( ! SCIPisInfinity(scip, rhs[i]) )
1146  {
1147  if ( ! sepadata->onlyactiverows || SCIPisFeasEQ(scip, val, SCIProwGetRhs(row)) )
1148  {
1149  if ( SCIProwGetNLPNonz(row) > maxrowsize )
1150  maxrowsize = SCIProwGetNLPNonz(row);
1151  if ( SCIProwGetNLPNonz(row) < minrowsize )
1152  minrowsize = SCIProwGetNLPNonz(row);
1153  }
1154  }
1155  }
1156  }
1157  assert( maxrowsize > 0 );
1158  assert( minrowsize < INT_MAX );
1159 
1160  /* add cuts for objective function if required */
1161  if ( sepadata->useobjub )
1162  {
1163  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
1164  rhs[mipdata->nrows] = SCIPgetUpperbound(scip);
1165  assert( ! SCIPisObjIntegral(scip) || SCIPisFeasIntegral(scip, SCIPgetUpperbound(scip)) );
1166 
1167  if ( ! SCIPisInfinity(scip, SCIPgetUpperbound(scip)) && SCIPgetNObjVars(scip) > maxrowsize )
1168  maxrowsize = SCIPgetNObjVars(scip);
1169  if ( ! SCIPisInfinity(scip, SCIPgetUpperbound(scip)) && SCIPgetNObjVars(scip) < minrowsize )
1170  minrowsize = SCIPgetNObjVars(scip);
1171  }
1172  if ( sepadata->useobjlb )
1173  {
1174  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
1175 
1176  if ( SCIPisObjIntegral(scip) )
1177  lhs[mipdata->nrows] = SCIPfeasCeil(scip, SCIPgetLowerbound(scip));
1178  else
1179  lhs[mipdata->nrows] = SCIPgetLowerbound(scip);
1180 
1181  if ( ! SCIPisInfinity(scip, -SCIPgetLowerbound(scip)) && SCIPgetNObjVars(scip) > maxrowsize )
1182  maxrowsize = SCIPgetNObjVars(scip);
1183  if ( ! SCIPisInfinity(scip, -SCIPgetLowerbound(scip)) && SCIPgetNObjVars(scip) < minrowsize )
1184  minrowsize = SCIPgetNObjVars(scip);
1185  }
1186 
1187  /* store lb/ub for complementing and perform preprocessing */
1188  nshifted = 0;
1189  ncomplemented = 0;
1190  ncontconverted = 0;
1191  nintconverted = 0;
1192  nlbounds = 0;
1193  nubounds = 0;
1194  for (j = 0; j < ncols; ++j)
1195  {
1196  SCIP_COL* col;
1197  SCIP_VAR* var;
1198 
1199  col = cols[j];
1200  assert( col != NULL );
1201  var = SCIPcolGetVar(col);
1202  assert( var != NULL );
1203 
1204  primsol[j] = SCIPcolGetPrimsol(col);
1205  assert( SCIPisEQ(scip, SCIPgetVarSol(scip, var), primsol[j]) );
1206 
1207  lb[j] = SCIPvarGetLbGlobal(var);
1208  assert( SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPcolGetLb(col)) );
1209 
1210  /* if allowed, try to use stronger local bound */
1211  if ( sepadata->allowlocal && SCIPisGT(scip, SCIPvarGetLbLocal(var), lb[j]) )
1212  lb[j] = SCIPvarGetLbLocal(var);
1213 
1214  ub[j] = SCIPvarGetUbGlobal(var);
1215  assert( SCIPisEQ(scip, SCIPvarGetUbLocal(var), SCIPcolGetUb(col)) );
1216 
1217  /* if allowed, try to use stronger local bound */
1218  if ( sepadata->allowlocal && SCIPisLT(scip, SCIPvarGetUbLocal(var), ub[j]) )
1219  ub[j] = SCIPvarGetUbLocal(var);
1220 
1221  mipdata->coltype[j] = colPresent;
1222  mipdata->iscomplemented[j] = FALSE;
1223  mipdata->isshifted[j] = FALSE;
1224 
1225  /* check status of column/variable */
1226  if ( SCIPcolIsIntegral(col) )
1227  {
1228  /* integral variables taking integral values are not interesting - will be substituted out below */
1229  if ( ! SCIPisFeasIntegral(scip, primsol[j]) )
1230  {
1231  /* possibly convert fractional integral variables to take integral values */
1232  if ( sepadata->intconvert && ncols >= sepadata->intconvmin )
1233  {
1234  /* randomly convert variables */
1235  if ( ((SCIP_Real) rand())/((SCIP_Real) RAND_MAX) <= sepadata->intconvfrac )
1236  {
1237  assert( ! SCIPisInfinity(scip, ub[j]) || ! SCIPisInfinity(scip, -lb[j]) );
1238 
1239  /* if both bounds are finite, take the closer one */
1240  if ( ! SCIPisInfinity(scip, ub[j]) && ! SCIPisInfinity(scip, -lb[j]) )
1241  {
1242  assert( SCIPisFeasIntegral(scip, ub[j]) );
1243  assert( SCIPisFeasIntegral(scip, lb[j]) );
1244  assert( SCIPisFeasLT(scip, primsol[j], ub[j]) );
1245  assert( SCIPisFeasGT(scip, primsol[j], lb[j]) );
1246  if ( ub[j] - primsol[j] < primsol[j] - lb[j] )
1247  primsol[j] = ub[j];
1248  else
1249  primsol[j] = lb[j];
1250  ++nintconverted;
1251  }
1252  else
1253  {
1254  /* if only lower bound is finite */
1255  if ( ! SCIPisInfinity(scip, -lb[j]) )
1256  {
1257  assert( SCIPisFeasIntegral(scip, lb[j]) );
1258  primsol[j] = lb[j];
1259  ++nintconverted;
1260  }
1261  else
1262  {
1263  assert( ! SCIPisInfinity(scip, ub[j]) );
1264  assert( SCIPisFeasIntegral(scip, ub[j]) );
1265  primsol[j] = ub[j];
1266  ++nintconverted;
1267  }
1268  }
1269  }
1270  }
1271  }
1272 
1273  /* integral variables taking integral values are not interesting - will be substituted out below */
1274  if ( ! SCIPisFeasIntegral(scip, primsol[j]) )
1275  {
1276  /* possibly convert integral variables to be continuous */
1277  if ( sepadata->contconvert && ncols >= sepadata->contconvmin )
1278  {
1279  /* randomly convert variables */
1280  if ( ((SCIP_Real) rand())/((SCIP_Real) RAND_MAX) <= sepadata->contconvfrac )
1281  {
1282  /* preprocessing is also performed for converted columns */
1283  mipdata->coltype[j] = colConverted;
1284  ++ncontconverted;
1285  }
1286  }
1287  }
1288  }
1289  else
1290  {
1291  /* detect continuous variables, but perform preprocessing for them */
1292  mipdata->coltype[j] = colContinuous;
1293  }
1294 
1295  /* if integer variable is at its upper bound -> complementing (this also generates a 0 lower bound) */
1296  if ( mipdata->coltype[j] == colPresent && SCIPisFeasEQ(scip, primsol[j], ub[j]) )
1297  {
1298  assert( ! SCIPisInfinity(scip, ub[j]) );
1299  SCIP_CALL( transformColumn(scip, sepadata, mipdata, col, ub[j], -1.0, lhs, rhs, &(lb[j]), &(ub[j]), &(primsol[j])) );
1300  mipdata->iscomplemented[j] = TRUE;
1301  mipdata->coltype[j] = colAtUb;
1302  ++nubounds;
1303  }
1304  else
1305  {
1306  /* if a variable has a finite nonzero lower bound -> shift */
1307  if ( ! SCIPisInfinity(scip, -lb[j]) )
1308  {
1309  if ( ! SCIPisZero(scip, lb[j]) )
1310  {
1311  SCIP_CALL( transformColumn(scip, sepadata, mipdata, col, -lb[j], 1.0, lhs, rhs, &(lb[j]), &(ub[j]), &(primsol[j])) );
1312  assert( SCIPisZero(scip, lb[j]) );
1313  mipdata->isshifted[j] = TRUE;
1314  ++nshifted;
1315  }
1316 
1317  /* if integer variable is at its lower bound */
1318  if ( mipdata->coltype[j] == colPresent && SCIPisZero(scip, primsol[j]) )
1319  {
1320  mipdata->coltype[j] = colAtLb;
1321  ++nlbounds;
1322  }
1323  }
1324  else
1325  {
1326  /* lower bound is minus-infinity -> check whether upper bound is finite */
1327  if ( ! SCIPisInfinity(scip, ub[j]) )
1328  {
1329  /* complement variable */
1330  SCIP_CALL( transformColumn(scip, sepadata, mipdata, col, ub[j], -1.0, lhs, rhs, &(lb[j]), &(ub[j]), &(primsol[j])) );
1331  assert( SCIPisZero(scip, lb[j]) );
1332  mipdata->iscomplemented[j] = TRUE;
1333  ++ncomplemented;
1334 
1335  /* if integer variable is at its lower bound */
1336  if ( mipdata->coltype[j] == colPresent && SCIPisZero(scip, primsol[j]) )
1337  {
1338  mipdata->coltype[j] = colAtLb;
1339  ++nlbounds;
1340  }
1341  }
1342  }
1343  }
1344 
1345  assert( SCIPisFeasLE(scip, lb[j], primsol[j]) );
1346  assert( SCIPisFeasLE(scip, primsol[j], ub[j]) );
1347  }
1348 
1349 #ifndef NDEBUG
1350  if ( sepadata->intconvert && ncols >= sepadata->intconvmin )
1351  {
1352  SCIPdebugMessage("Converted %u fractional integral variables to have integral value.\n", nintconverted);
1353  }
1354  if ( sepadata->contconvert && ncols >= sepadata->contconvmin )
1355  {
1356  SCIPdebugMessage("Converted %u integral variables to be continuous.\n", ncontconverted);
1357  }
1358 #endif
1359  SCIPdebugMessage("original variables: %d integral, %d continuous, %u shifted, %u complemented, %u at lb, %u at ub\n",
1361  nshifted, ncomplemented, nlbounds, nubounds);
1362 
1363  /* prepare upper bound on y-variables */
1364  if ( sepadata->skipmultbounds )
1365  multvarub = SCIPinfinity(scip);
1366  else
1367  multvarub = 1.0-EPSILONVALUE;
1368 
1369  /* create artificial variables for row combinations (y-variables) */
1370  cnt = 0;
1371  for (i = 0; i < nrows; ++i)
1372  {
1373  SCIP_ROW* row;
1374 
1375  row = rows[i];
1376  assert( row != NULL );
1377 
1378  mipdata->ylhs[i] = NULL;
1379  mipdata->yrhs[i] = NULL;
1380 
1381  /* skip modifiable rows and local rows, unless allowed */
1382  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
1383  continue;
1384 
1385  /* skip rows that not have been active for a longer time */
1386  if ( ! sepadata->onlyactiverows && sepadata->maxrowage > 0 && SCIProwGetAge(row) > sepadata->maxrowage )
1387  continue;
1388 
1389  /* check whether we want to skip cut produced by the CGMIP separator */
1390  if ( sepadata->onlyrankone )
1391  {
1392  if ( SCIProwGetOriginSepa(row) == sepa )
1393  continue;
1394  }
1395 
1396  /* if we have an equation */
1397  if ( SCIPisEQ(scip, lhs[i], rhs[i]) )
1398  {
1399  SCIP_Real weight = -sepadata->objweight;
1400 
1401  assert( ! SCIPisInfinity(scip, rhs[i]) );
1402  assert( SCIPisFeasEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row)) ); /* equations should always be active */
1403  assert( SCIPisFeasEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetRhs(row)) );
1404 
1405  if ( sepadata->objweightsize )
1406  weight = - sepadata->objweight * computeObjWeightSize(SCIProwGetNLPNonz(row), minrowsize, maxrowsize);
1407 
1408  /* create two variables for each equation */
1409  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yeq1_%d", i);
1410  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->ylhs[i]), name, 0.0, multvarub,
1412  SCIP_CALL( SCIPaddVar(subscip, mipdata->ylhs[i]) );
1413  ++cnt;
1414 
1415 #ifdef SCIP_MORE_DEBUG
1416  SCIPdebugMessage("Created variable <%s> for equation <%s>.\n", name, SCIProwGetName(row));
1417 #endif
1418 
1419  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yeq2_%d", i);
1420  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->yrhs[i]), name, 0.0, multvarub,
1422  SCIP_CALL( SCIPaddVar(subscip, mipdata->yrhs[i]) );
1423  ++cnt;
1424 
1425 #ifdef SCIP_MORE_DEBUG
1426  SCIPdebugMessage("Created variable <%s> for equation <%s>.\n", name, SCIProwGetName(row));
1427 #endif
1428  }
1429  else
1430  {
1431  /* create variable for lhs of row if necessary */
1432  if ( ! SCIPisInfinity(scip, -lhs[i]) )
1433  {
1434  SCIP_Bool isactive = FALSE;
1435  SCIP_Real weight = 0.0;
1436 
1437  /* if the row is active, use objective weight equal to -sepadata->objweight */
1438  if ( SCIPisFeasEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row)) )
1439  {
1440  isactive = TRUE;
1441  if ( sepadata->objweightsize )
1442  weight = -sepadata->objweight * computeObjWeightSize(SCIProwGetNLPNonz(row), minrowsize, maxrowsize);
1443  else
1444  weight = -sepadata->objweight;
1445  }
1446 
1447  if ( ! sepadata->onlyactiverows || isactive )
1448  {
1449  /* add variable */
1450  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "ylhs_%d", i);
1451  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->ylhs[i]), name, 0.0, multvarub,
1453  SCIP_CALL( SCIPaddVar(subscip, mipdata->ylhs[i]) );
1454  ++cnt;
1455 
1456 #ifdef SCIP_MORE_DEBUG
1457  SCIPdebugMessage("Created variable <%s> for >= inequality <%s> (weight: %f).\n", name, SCIProwGetName(row), weight);
1458 #endif
1459  }
1460  }
1461 
1462  /* create variable for rhs of row if necessary */
1463  if ( ! SCIPisInfinity(scip, rhs[i]) )
1464  {
1465  SCIP_Bool isactive = FALSE;
1466  SCIP_Real weight = 0.0;
1467 
1468  /* if the row is active, use objective weight equal to -sepadata->objweight */
1469  if ( SCIPisFeasEQ(scip, SCIPgetRowLPActivity(scip, row), SCIProwGetRhs(row)) )
1470  {
1471  isactive = TRUE;
1472  if ( sepadata->objweightsize )
1473  weight = -sepadata->objweight * computeObjWeightSize(SCIProwGetNLPNonz(row), minrowsize, maxrowsize);
1474  else
1475  weight = -sepadata->objweight;
1476  }
1477 
1478  if ( ! sepadata->onlyactiverows || isactive )
1479  {
1480  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yrhs_%d", i);
1481  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->yrhs[i]), name, 0.0, multvarub,
1483  SCIP_CALL( SCIPaddVar(subscip, mipdata->yrhs[i]) );
1484  ++cnt;
1485 
1486 #ifdef SCIP_MORE_DEBUG
1487  SCIPdebugMessage("Created variable <%s> for <= inequality <%s> (weight: %f).\n", name, SCIProwGetName(row), weight);
1488 #endif
1489  }
1490  }
1491  }
1492  }
1493  assert( (int) cnt <= 2 * nrows );
1494  mipdata->n += cnt;
1495 
1496  /* create artificial variables for objective function (if required) (y-variables) */
1497  if ( sepadata->useobjub || sepadata->useobjlb )
1498  {
1499  SCIP_Real weight = 0.0;
1500 
1501  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
1502  mipdata->ylhs[mipdata->nrows] = NULL;
1503  mipdata->yrhs[mipdata->nrows] = NULL;
1504  cnt = 0;
1505 
1506  if ( sepadata->objweightsize )
1507  weight = -sepadata->objweight * computeObjWeightSize(SCIPgetNObjVars(scip), minrowsize, maxrowsize);
1508  else
1509  weight = -sepadata->objweight;
1510 
1511  /* create variable for upper objective bound if necessary */
1512  if ( sepadata->useobjub && ! SCIPisInfinity(scip, rhs[mipdata->nrows]) )
1513  {
1514  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yobjub");
1515  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->yrhs[mipdata->nrows]), name, 0.0, multvarub,
1517  SCIP_CALL( SCIPaddVar(subscip, mipdata->yrhs[mipdata->nrows]) );
1518  ++cnt;
1519 
1520 #ifdef SCIP_MORE_DEBUG
1521  SCIPdebugMessage("Created variable <%s> for upper bound on objective (weight: %f).\n", name, weight);
1522 #endif
1523  }
1524 
1525  /* create variable for lower bound objective if necessary */
1526  if ( sepadata->useobjlb && ! SCIPisInfinity(scip, -lhs[mipdata->nrows]) )
1527  {
1528  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "yobjlb");
1529  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->ylhs[mipdata->nrows]), name, 0.0, multvarub,
1531  SCIP_CALL( SCIPaddVar(subscip, mipdata->ylhs[mipdata->nrows]) );
1532  ++cnt;
1533 
1534 #ifdef SCIP_MORE_DEBUG
1535  SCIPdebugMessage("Created variable <%s> for lower bound on objective (weight: %f).\n", name, weight);
1536 #endif
1537  }
1538 
1539  assert( (int) cnt <= 2 * ntotalrows );
1540  mipdata->n += cnt;
1541  }
1542 
1543  /* create alpha, bound, and fractional variables */
1544  cnt = 0;
1545  ucnt = 0;
1546  for (j = 0; j < ncols; ++j)
1547  {
1548  mipdata->z[j] = NULL;
1549  mipdata->alpha[j] = NULL;
1550  mipdata->fracalpha[j] = NULL;
1551 
1552  if ( mipdata->coltype[j] == colPresent )
1553  {
1554  SCIP_Real obj;
1555 
1556  if ( sepadata->objlone )
1557  obj = 0.0;
1558  else
1559  obj = primsol[j];
1560 
1561  /* create alpha variables */
1562  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "alpha_%d", j);
1563  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->alpha[j]), name, -sepadata->cutcoefbnd, sepadata->cutcoefbnd, obj,
1565  SCIP_CALL( SCIPaddVar(subscip, mipdata->alpha[j]) );
1566  ++cnt;
1567 
1568  /* create fractional variables */
1569  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "f_%d", j);
1570  if ( SCIPisInfinity(scip, -lb[j]) && SCIPisInfinity(scip, ub[j]) )
1571  {
1572  /* fix fractional value to be zero for free original variables */
1573  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->fracalpha[j]), name, 0.0, 0.0, 0.0,
1575  }
1576  else
1577  {
1578  /* fractional value in [0, 1) for variables with finite bounds */
1579  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->fracalpha[j]), name, 0.0, 1.0-EPSILONVALUE, 0.0,
1581  }
1582  SCIP_CALL( SCIPaddVar(subscip, mipdata->fracalpha[j]) );
1583  ++cnt;
1584 
1585  /* create variables for upper bounds */
1586  if ( ! SCIPisInfinity(scip, ub[j]) )
1587  {
1588  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "zub_%d", j);
1589  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->z[j]), name, 0.0, multvarub,
1591  SCIP_CALL( SCIPaddVar(subscip, mipdata->z[j]) );
1592  ++ucnt;
1593  }
1594  }
1595  }
1596  assert( (int) cnt <= 2 * ncols );
1597  assert( (int) ucnt <= ncols );
1598 
1599  /* create variable for the rhs of the cut */
1600  if ( sepadata->objlone )
1601  {
1602  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->beta), "beta", -sepadata->cutcoefbnd, sepadata->cutcoefbnd, 0.0,
1604  }
1605  else
1606  {
1607  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->beta), "beta", -sepadata->cutcoefbnd, sepadata->cutcoefbnd, -1.0,
1609  }
1610  SCIP_CALL( SCIPaddVar(subscip, mipdata->beta) );
1611 
1612  /* create fractional variable for the rhs */
1613  SCIP_CALL( SCIPcreateVar(subscip, &(mipdata->fracbeta), "fracbeta", 0.0, 1.0-BETAEPSILONVALUE, 0.0,
1615  SCIP_CALL( SCIPaddVar(subscip, mipdata->fracbeta) );
1616  mipdata->n += cnt + ucnt + 2;
1617 
1618  /* get temporary storage */
1619  SCIP_CALL( SCIPallocBufferArray(scip, &consvals, (int) mipdata->n) );
1620  SCIP_CALL( SCIPallocBufferArray(scip, &consvars, (int) mipdata->n) );
1621 
1622  /* create constraints for alpha variables of CG-cut */
1623  cnt = 0;
1624  for (j = 0; j < ncols; ++j)
1625  {
1626  SCIP_ROW** colrows;
1627  SCIP_Real* colvals;
1628 
1629  /* create ordinary part for all selected variables */
1630  if ( mipdata->coltype[j] == colPresent )
1631  {
1632  SCIP_Real sigma;
1633 
1634  assert( cols[j] != NULL );
1635  colrows = SCIPcolGetRows(cols[j]);
1636  colvals = SCIPcolGetVals(cols[j]);
1637  nconsvars = 0;
1638 
1639  if ( mipdata->iscomplemented[j] )
1640  sigma = -1.0;
1641  else
1642  sigma = 1.0;
1643 
1644  /* add part for columns */
1645  for (i = 0; i < SCIPcolGetNLPNonz(cols[j]); ++i)
1646  {
1647  SCIP_ROW* row;
1648  int pos;
1649 
1650  row = colrows[i];
1651  assert( row != NULL );
1652 
1653  /* skip modifiable rows and local rows, unless allowed */
1654  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
1655  continue;
1656 
1657  pos = SCIProwGetLPPos(row);
1658  assert( 0 <= pos && pos < nrows );
1659 
1660  if ( mipdata->ylhs[pos] != NULL )
1661  {
1662  consvars[nconsvars] = mipdata->ylhs[pos];
1663  consvals[nconsvars] = -sigma * colvals[i];
1664  ++nconsvars;
1665  }
1666  if ( mipdata->yrhs[pos] != NULL )
1667  {
1668  consvars[nconsvars] = mipdata->yrhs[pos];
1669  consvals[nconsvars] = sigma * colvals[i];
1670  ++nconsvars;
1671  }
1672  assert( nconsvars <= (int) mipdata->n );
1673  }
1674  /* add part for upper bounds */
1675  if ( mipdata->z[j] != NULL )
1676  {
1677  assert( ! SCIPisInfinity(scip, ub[j]) );
1678  consvars[nconsvars] = mipdata->z[j];
1679  consvals[nconsvars] = 1.0;
1680  ++nconsvars;
1681  }
1682  assert( nconsvars <= (int) mipdata->n );
1683 
1684  /* add alpha variable */
1685  consvars[nconsvars] = mipdata->alpha[j];
1686  consvals[nconsvars] = -1.0;
1687  ++nconsvars;
1688  assert( nconsvars <= (int) mipdata->n );
1689 
1690  /* add fractional-alpha variable */
1691  consvars[nconsvars] = mipdata->fracalpha[j];
1692  consvals[nconsvars] = -1.0;
1693  ++nconsvars;
1694  assert( nconsvars <= (int) mipdata->n );
1695 
1696  /* check for lower and upper objective bounds */
1697  if ( (sepadata->useobjub || sepadata->useobjlb) && ! SCIPisZero(scip, SCIPcolGetObj(cols[j])) )
1698  {
1699  /* add lower objective bound */
1700  if ( mipdata->ylhs[mipdata->nrows] != NULL )
1701  {
1702  assert( sepadata->useobjlb );
1703  consvars[nconsvars] = mipdata->ylhs[mipdata->nrows];
1704  consvals[nconsvars] = -sigma * SCIPcolGetObj(cols[j]);
1705  ++nconsvars;
1706  }
1707 
1708  /* add upper objective bound */
1709  if ( mipdata->yrhs[mipdata->nrows] != NULL )
1710  {
1711  assert( sepadata->useobjub );
1712  consvars[nconsvars] = mipdata->yrhs[mipdata->nrows];
1713  consvals[nconsvars] = sigma * SCIPcolGetObj(cols[j]);
1714  ++nconsvars;
1715  }
1716  }
1717 
1718  /* add linear constraint */
1719  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "alpha_%d", j);
1720  SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, name, nconsvars, consvars, consvals, 0.0, 0.0,
1721  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1722  SCIP_CALL( SCIPaddCons(subscip, cons) );
1723  SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
1724  ++cnt;
1725  }
1726  /* generate part that makes sure that cut is valid for continuous variables */
1727  else if ( mipdata->coltype[j] == colContinuous || mipdata->coltype[j] == colConverted )
1728  {
1729  SCIP_Real sigma;
1730  SCIP_Real r;
1731 
1732  assert( cols[j] != NULL );
1733  colrows = SCIPcolGetRows(cols[j]);
1734  colvals = SCIPcolGetVals(cols[j]);
1735  nconsvars = 0;
1736 
1737  if ( mipdata->iscomplemented[j] )
1738  sigma = -1.0;
1739  else
1740  sigma = 1.0;
1741 
1742  /* add part for columns */
1743  for (i = 0; i < SCIPcolGetNLPNonz(cols[j]); ++i)
1744  {
1745  SCIP_ROW* row;
1746  int pos;
1747 
1748  row = colrows[i];
1749  assert( row != NULL );
1750 
1751  /* skip modifiable rows and local rows, unless allowed */
1752  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
1753  continue;
1754 
1755  pos = SCIProwGetLPPos(row);
1756  assert( 0 <= pos && pos < nrows );
1757 
1758  if ( mipdata->ylhs[pos] != NULL )
1759  {
1760  consvars[nconsvars] = mipdata->ylhs[pos];
1761  consvals[nconsvars] = -sigma * colvals[i];
1762  ++nconsvars;
1763  }
1764  if ( mipdata->yrhs[pos] != NULL )
1765  {
1766  consvars[nconsvars] = mipdata->yrhs[pos];
1767  consvals[nconsvars] = sigma * colvals[i];
1768  ++nconsvars;
1769  }
1770  assert( nconsvars <= (int) mipdata->n );
1771  }
1772 
1773  /* check for lower and upper objective bounds */
1774  if ( (sepadata->useobjub || sepadata->useobjlb) && ! SCIPisZero(scip, SCIPcolGetObj(cols[j])) )
1775  {
1776  /* add lower objective bound */
1777  if ( mipdata->ylhs[mipdata->nrows] )
1778  {
1779  assert( sepadata->useobjlb );
1780  consvars[nconsvars] = mipdata->ylhs[mipdata->nrows];
1781  consvals[nconsvars] = -sigma * SCIPcolGetObj(cols[j]);
1782  ++nconsvars;
1783  }
1784 
1785  /* add upper objective bound */
1786  if ( mipdata->yrhs[mipdata->nrows] )
1787  {
1788  assert( sepadata->useobjub );
1789  consvars[nconsvars] = mipdata->yrhs[mipdata->nrows];
1790  consvals[nconsvars] = sigma * SCIPcolGetObj(cols[j]);
1791  ++nconsvars;
1792  }
1793  }
1794 
1795  /* add linear constraint */
1796  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cont_%d", j);
1797 
1798  /* for free continuous variables require equality */
1799  r = SCIPinfinity(subscip);
1800  if ( SCIPisInfinity(scip, -lb[j]) && SCIPisInfinity(scip, ub[j]) )
1801  r = 0.0;
1802  else
1803  assert( SCIPisZero(scip, lb[j]) );
1804 
1805  SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, name, nconsvars, consvars, consvals, 0.0, r,
1806  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1807  SCIP_CALL( SCIPaddCons(subscip, cons) );
1808  SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
1809  ++cnt;
1810  }
1811  }
1812  assert( (int) cnt <= ncols );
1813  mipdata->m += cnt;
1814 
1815  /* create constraints for rhs of cut */
1816  nconsvars = 0;
1817 
1818  /* first for the rows */
1819  for (i = 0; i < nrows; ++i)
1820  {
1821  assert( rows[i] != NULL );
1822 
1823  /* skip modifiable rows and local rows, unless allowed */
1824  if ( SCIProwIsModifiable(rows[i]) || (SCIProwIsLocal(rows[i]) && !sepadata->allowlocal) )
1825  continue;
1826 
1827  /* if lhs is there */
1828  if ( mipdata->ylhs[i] != NULL && ! SCIPisZero(scip, lhs[i]) )
1829  {
1830  assert( ! SCIPisInfinity(scip, -lhs[i]) );
1831  consvars[nconsvars] = mipdata->ylhs[i];
1832  consvals[nconsvars] = -lhs[i];
1833  ++nconsvars;
1834  }
1835  /* if rhs is there */
1836  if ( mipdata->yrhs[i] != NULL && ! SCIPisZero(scip, rhs[i]) )
1837  {
1838  assert( ! SCIPisInfinity(scip, rhs[i]) );
1839  consvars[nconsvars] = mipdata->yrhs[i];
1840  consvals[nconsvars] = rhs[i];
1841  ++nconsvars;
1842  }
1843  assert( nconsvars <= (int) mipdata->n );
1844  }
1845 
1846  if ( sepadata->useobjub || sepadata->useobjlb )
1847  {
1848  /* add lower objective bound */
1849  if ( mipdata->ylhs[mipdata->nrows] != NULL && ! SCIPisZero(scip, lhs[mipdata->nrows]) )
1850  {
1851  assert( sepadata->useobjlb );
1852  assert( ! SCIPisInfinity(scip, -lhs[mipdata->nrows]) );
1853  consvars[nconsvars] = mipdata->ylhs[mipdata->nrows];
1854  consvals[nconsvars] = -lhs[mipdata->nrows];
1855  ++nconsvars;
1856  }
1857 
1858  /* add upper objective bound */
1859  if ( mipdata->yrhs[mipdata->nrows] != NULL && ! SCIPisZero(scip, rhs[mipdata->nrows]) )
1860  {
1861  assert( sepadata->useobjub );
1862  assert( ! SCIPisInfinity(scip, rhs[mipdata->nrows]) );
1863  consvars[nconsvars] = mipdata->yrhs[mipdata->nrows];
1864  consvals[nconsvars] = rhs[mipdata->nrows];
1865  ++nconsvars;
1866  }
1867  assert( nconsvars <= (int) mipdata->n );
1868  }
1869 
1870  /* next for the columns */
1871  for (j = 0; j < ncols; ++j)
1872  {
1873  /* if ub is there */
1874  if ( mipdata->z[j] != NULL && ! SCIPisZero(scip, ub[j]) )
1875  {
1876  assert( mipdata->coltype[j] == colPresent );
1877  assert( ! SCIPisInfinity(scip, ub[j]) );
1878  consvars[nconsvars] = mipdata->z[j];
1879  consvals[nconsvars] = ub[j];
1880  ++nconsvars;
1881  assert( nconsvars <= (int) mipdata->n );
1882  }
1883  }
1884  /* add beta variable */
1885  consvars[nconsvars] = mipdata->beta;
1886  consvals[nconsvars] = -1.0;
1887  ++nconsvars;
1888 
1889  /* add fractional-beta variable */
1890  consvars[nconsvars] = mipdata->fracbeta;
1891  consvals[nconsvars] = -1.0;
1892  ++nconsvars;
1893  assert( nconsvars <= (int) mipdata->n );
1894 
1895  /* add linear constraint */
1896  SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, "beta", nconsvars, consvars, consvals, 0.0, 0.0,
1897  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1898  SCIP_CALL( SCIPaddCons(subscip, cons) );
1899  SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
1900  ++mipdata->m;
1901 
1902  /* add primal separation constraint if required */
1903  if ( sepadata->primalseparation )
1904  {
1905  SCIP_SOL* bestsol;
1906  bestsol = SCIPgetBestSol(scip);
1907  if ( bestsol != NULL )
1908  {
1909  nconsvars = 0;
1910  for (j = 0; j < ncols; ++j)
1911  {
1912  if ( mipdata->alpha[j] != NULL )
1913  {
1914  SCIP_Real val;
1915  assert( mipdata->coltype[j] == colPresent );
1916 
1917  val = SCIPgetSolVal(scip, bestsol, SCIPcolGetVar(cols[j]));
1918  consvars[nconsvars] = mipdata->alpha[j];
1919  consvals[nconsvars] = val;
1920  ++nconsvars;
1921  assert( nconsvars <= (int) mipdata->n );
1922  }
1923  }
1924  consvars[nconsvars] = mipdata->beta;
1925  consvals[nconsvars] = -1.0;
1926  ++nconsvars;
1927 
1928  /* add linear constraint - allow slight deviation from equality */
1929  SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, "primalseparation", nconsvars, consvars, consvals, -EPSILONVALUE, EPSILONVALUE,
1930  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1931  SCIP_CALL( SCIPaddCons(subscip, cons) );
1932  SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
1933  ++mipdata->m;
1934  }
1935  }
1936 
1937  /* add constraint to force violated cuts if required */
1938  if ( sepadata->addviolationcons )
1939  {
1940  nconsvars = 0;
1941  for (j = 0; j < ncols; ++j)
1942  {
1943  if ( mipdata->alpha[j] != NULL )
1944  {
1945  consvars[nconsvars] = mipdata->alpha[j];
1946  consvals[nconsvars] = primsol[j];
1947  ++nconsvars;
1948  assert( nconsvars <= (int) mipdata->n );
1949  }
1950  }
1951  consvars[nconsvars] = mipdata->beta;
1952  consvals[nconsvars] = -1.0;
1953  ++nconsvars;
1954 
1955  /* add linear constraint - allow slight deviation from equality */
1956  SCIP_CALL( SCIPcreateConsLinear(subscip, &cons, "violationConstraint", nconsvars, consvars, consvals, MINEFFICACY, SCIPinfinity(subscip),
1957  TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE) );
1958  SCIP_CALL( SCIPaddCons(subscip, cons) );
1959  SCIP_CALL( SCIPreleaseCons(subscip, &cons) );
1960  ++mipdata->m;
1961  }
1962 
1963  SCIPdebugMessage("Subscip has %u vars (%d integral, %d continuous), %u conss.\n",
1964  mipdata->n, SCIPgetNIntVars(subscip), SCIPgetNContVars(subscip), mipdata->m);
1965 
1966  /* free temporary memory */
1967  SCIPfreeBufferArray(scip, &consvars);
1968  SCIPfreeBufferArray(scip, &consvals);
1969 
1970  SCIPfreeBufferArray(scip, &primsol);
1971  SCIPfreeBufferArray(scip, &lb);
1972  SCIPfreeBufferArray(scip, &ub);
1973  SCIPfreeBufferArray(scip, &rhs);
1974  SCIPfreeBufferArray(scip, &lhs);
1975 
1976  /* SCIPdebug( SCIP_CALL( SCIPprintOrigProblem(subscip, NULL, NULL, FALSE) ) ); */
1977 
1978 #ifdef SCIP_WRITEPROB
1979  {
1980  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cgsepa%s%s%s%s_%s.lp",
1981  sepadata->objlone ? "_l1" : "",
1982  sepadata->addviolationcons ? "_vc" : "",
1983  sepadata->skipmultbounds ? "_ub" : "",
1984  sepadata->primalseparation ? "_ps" : "",
1985  SCIPgetProbName(scip));
1986  SCIP_CALL( SCIPwriteOrigProblem(subscip, name, "lp", FALSE) );
1987  SCIPinfoMessage(scip, NULL, "Wrote subscip to file <%s>.\n", name);
1988  }
1989 #endif
1990 
1991  return SCIP_OKAY;
1992 }
1993 
1994 
1995 /** sets parameters for subscip */
1996 static
1998  SCIP_SEPADATA* sepadata, /**< separator data */
1999  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
2000  SCIP_Bool* success /**< if setting was successful -> stop solution otherwise */
2001  )
2002 {
2003  SCIP* subscip;
2004 
2005  assert( sepadata != NULL );
2006  assert( mipdata != NULL );
2007  assert( success != NULL );
2008 
2009  *success = TRUE;
2010 
2011  subscip = mipdata->subscip;
2012  assert( subscip != NULL );
2013 
2014  /* set objective limit, if no corresponding constraint has been added */
2015  if ( ! sepadata->addviolationcons && ! sepadata->addviolconshdlr )
2016  {
2017  SCIP_CALL( SCIPsetObjlimit(subscip, MINEFFICACY) );
2018  }
2019 
2020  /* do not abort subscip on CTRL-C */
2021  SCIP_CALL( SCIPsetBoolParam(subscip, "misc/catchctrlc", FALSE) );
2022 
2023  /* disable memory saving mode: this is likely to result in the maximal depth being reached. This is because DFS
2024  * results in a repeated branching on the alpha-variables, which often have large bounds resulting in deep levels of
2025  * the tree. */
2026  SCIP_CALL( SCIPsetRealParam(subscip, "memory/savefac", 1.0) );
2027 
2028  /* set number of solutions stored */
2029  SCIP_CALL( SCIPsetIntParam(subscip, "limits/maxsol", MAXNSOLS) );
2030 
2031  /* determine output to console */
2032 #ifdef SCIP_OUTPUT
2033  SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
2034  SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 1000) );
2035  SCIP_CALL( SCIPsetIntParam(subscip, "display/nsols/active", 2) );
2036 #else
2037  if ( sepadata->output )
2038  {
2039  SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 5) );
2040  SCIP_CALL( SCIPsetIntParam(subscip, "display/freq", 1000) );
2041  SCIP_CALL( SCIPsetIntParam(subscip, "display/nsols/active", 2) );
2042  }
2043  else
2044  {
2045  SCIP_CALL( SCIPsetIntParam(subscip, "display/verblevel", 0) );
2046  }
2047 #endif
2048 
2049  if ( sepadata->subscipfast )
2050  {
2051  /* forbid recursive call of plugins solving subMIPs (also disables CG-separation) */
2052 #ifdef SCIP_OUTPUT
2053  SCIP_CALL( SCIPsetSubscipsOff(subscip, FALSE) );
2054 #else
2055  SCIP_CALL( SCIPsetSubscipsOff(subscip, TRUE) ); /* quiet */
2056 #endif
2057  }
2058  else
2059  {
2060  /* avoid recursive call */
2061  if ( ! SCIPisParamFixed(subscip, "separating/cgmip/freq") )
2062  {
2063  SCIP_CALL( SCIPsetIntParam(subscip, "separating/cgmip/freq", -1) );
2064  }
2065  }
2066 
2067 #if 0
2069 #else
2070 
2071  /* zirounding is often successful, so allow it some more calls */
2072  if( !SCIPisParamFixed(subscip, "heuristics/zirounding/minstopncalls") )
2073  {
2074  SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/zirounding/minstopncalls", 10000) );
2075  }
2076 
2077  if ( sepadata->subscipfast )
2078  {
2079  /* set other heuristics */
2080  if( !SCIPisParamFixed(subscip, "heuristics/shifting/freq") )
2081  {
2082  SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/shifting/freq", 3) );
2083  }
2084  if( !SCIPisParamFixed(subscip, "heuristics/simplerounding/freq") )
2085  {
2086  SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/simplerounding/freq", 1) );
2087  }
2088  if( !SCIPisParamFixed(subscip, "heuristics/rounding/freq") )
2089  {
2090  SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/rounding/freq", 1) );
2091  }
2092  if( !SCIPisParamFixed(subscip, "heuristics/oneopt/freq") )
2093  {
2094  SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/oneopt/freq", 1) );
2095  }
2096 
2097  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/pscostdiving/freq", 1) ); */
2098  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/feaspump/freq", 3) ); */
2099 
2100  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/coefdiving/freq", -1) ); */
2101  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/fracdiving/freq", -1) ); */
2102  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/guideddiving/freq", -1) ); */
2103  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/linesearchdiving/freq", -1) ); */
2104  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/objpscostdiving/freq", -1) ); */
2105  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/rootsoldiving/freq", -1) ); */
2106  /* SCIP_CALL( SCIPsetIntParam(subscip, "heuristics/veclendiving/freq", -1) ); */
2107 
2108  /* use fast presolving */
2110 
2111  /* disable conflict analysis */
2112  /* SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/useprop", FALSE) ); */
2113  /* SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/useinflp", FALSE) ); */
2114  /* SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/useboundlp", FALSE) ); */
2115  /* SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/usesb", FALSE) ); */
2116  /* SCIP_CALL( SCIPsetBoolParam(subscip, "conflict/usepseudo", FALSE) ); */
2117 
2118  /* use fast separation */
2120  }
2121 #endif
2122 
2123  return SCIP_OKAY;
2124 }
2125 
2126 
2127 /** solve subscip */
2128 static
2130  SCIP* scip, /**< SCIP data structure */
2131  SCIP_SEPADATA* sepadata, /**< separator data */
2132  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
2133  SCIP_Bool* success /**< if setting was successful -> stop solution otherwise */
2134  )
2135 {
2136  SCIP* subscip;
2137  SCIP_RETCODE retcode;
2138  SCIP_STATUS status;
2139  SCIP_Real timelimit;
2140  SCIP_Real memorylimit;
2141  SCIP_Longint nodelimit;
2142 
2143  assert( scip != NULL );
2144  assert( sepadata != NULL );
2145  assert( mipdata != NULL );
2146  assert( success != NULL );
2147 
2148  *success = TRUE;
2149 
2150  subscip = mipdata->subscip;
2151 
2152  /* determine timelimit */
2153  SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
2154  if ( ! SCIPisInfinity(scip, timelimit) )
2155  timelimit -= SCIPgetSolvingTime(scip);
2156  if ( sepadata->timelimit < timelimit )
2157  timelimit = sepadata->timelimit;
2158  if ( timelimit > 0.0 )
2159  {
2160  SCIP_CALL( SCIPsetRealParam(subscip, "limits/time", timelimit) );
2161  }
2162  else
2163  {
2164  *success = FALSE;
2165  return SCIP_OKAY;
2166  }
2167 
2168  /* determine memorylimit */
2169  SCIP_CALL( SCIPgetRealParam(scip, "limits/memory", &memorylimit) );
2170  if ( sepadata->memorylimit < memorylimit )
2171  memorylimit = sepadata->memorylimit;
2172 
2173  /* substract the memory already used by the main SCIP and the estimated memory usage of external software */
2174  if ( ! SCIPisInfinity(scip, memorylimit) )
2175  {
2176  memorylimit -= SCIPgetMemUsed(scip)/1048576.0;
2177  memorylimit -= SCIPgetMemExternEstim(scip)/1048576.0;
2178  }
2179 
2180  /* set memory limit if at least twice the amount of memory is left as is currently used */
2181  if ( memorylimit > 2.0 * (SCIPgetMemUsed(scip) + SCIPgetMemExternEstim(scip)) / 1048576.0 )
2182  {
2183  SCIP_CALL( SCIPsetRealParam(subscip, "limits/memory", memorylimit) );
2184  }
2185  else
2186  {
2187  *success = FALSE;
2188  return SCIP_OKAY;
2189  }
2190 
2191  /* set nodelimit for subproblem */
2192  if ( sepadata->minnodelimit < 0 || sepadata->maxnodelimit < 0 )
2193  nodelimit = SCIP_LONGINT_MAX;
2194  else
2195  {
2196  assert( sepadata->minnodelimit >= 0 && sepadata->maxnodelimit >= 0 );
2197  nodelimit = SCIPgetNLPIterations(scip);
2198  nodelimit = MAX(sepadata->minnodelimit, nodelimit);
2199  nodelimit = MIN(sepadata->maxnodelimit, nodelimit);
2200  }
2201  assert( nodelimit >= 0 );
2202  SCIP_CALL( SCIPsetLongintParam(subscip, "limits/nodes", nodelimit) );
2203 
2204  SCIPdebugMessage("Solving sub-SCIP (time limit: %f mem limit: %f node limit: %"SCIP_LONGINT_FORMAT") ...\n", timelimit, memorylimit, nodelimit);
2205 
2206  /* check whether we want a complete solve */
2207  if ( ! sepadata->earlyterm )
2208  {
2209  retcode = SCIPsolve(subscip);
2210 
2211  /* errors in solving the subproblem should not kill the overall solving process;
2212  * hence, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. */
2213  if ( retcode != SCIP_OKAY )
2214  {
2215 #ifndef NDEBUG
2216  SCIP_CALL( retcode );
2217 #endif
2218  SCIPwarningMessage(scip, "Error while solving subproblem in CGMIP separator; sub-SCIP terminated with code <%d>\n", retcode);
2219  *success = FALSE;
2220  return SCIP_OKAY;
2221  }
2222 
2223  status = SCIPgetStatus(subscip);
2224 
2225 #ifdef SCIP_OUTPUT
2226  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2227 #else
2228  if ( sepadata->output )
2229  {
2230  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2231  }
2232 #endif
2233 
2234  /* if the solution process was terminated or the problem is infeasible (can happen because of violation constraint) */
2235  if ( status == SCIP_STATUS_TIMELIMIT || status == SCIP_STATUS_USERINTERRUPT || status == SCIP_STATUS_INFEASIBLE || status == SCIP_STATUS_INFORUNBD )
2236  {
2237  *success = FALSE;
2238  return SCIP_OKAY;
2239  }
2240 
2241  /* all other statuses except optimal are invalid */
2242  if ( status != SCIP_STATUS_OPTIMAL && status != SCIP_STATUS_NODELIMIT )
2243  {
2244  SCIPerrorMessage("Solution of subscip for CG-separation returned with invalid status %d.\n", status);
2245  return SCIP_ERROR;
2246  }
2247  }
2248  else
2249  {
2250  /* otherwise we want a heuristic solve */
2251 
2252  /* -> solve until first solution is found */
2253  SCIP_CALL( SCIPsetIntParam(subscip, "limits/bestsol", 1) );
2254  retcode = SCIPsolve(subscip);
2255  SCIP_CALL( SCIPsetIntParam(subscip, "limits/bestsol", -1) );
2256 
2257  /* errors in solving the subproblem should not kill the overall solving process;
2258  * hence, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. */
2259  if ( retcode != SCIP_OKAY )
2260  {
2261 #ifndef NDEBUG
2262  SCIP_CALL( retcode );
2263 #endif
2264  SCIPwarningMessage(scip, "Error while solving subproblem in CGMIP separator; sub-SCIP terminated with code <%d>\n", retcode);
2265  *success = FALSE;
2266  return SCIP_OKAY;
2267  }
2268 
2269  status = SCIPgetStatus(subscip);
2270 
2271  /* if the solution process was terminated or the problem is infeasible (can happen because of violation constraint) */
2272  if ( status == SCIP_STATUS_TIMELIMIT || status == SCIP_STATUS_USERINTERRUPT || status == SCIP_STATUS_NODELIMIT ||
2273  status == SCIP_STATUS_INFEASIBLE || status == SCIP_STATUS_INFORUNBD || status == SCIP_STATUS_MEMLIMIT )
2274  {
2275  /* output statistics before stopping */
2276 #ifdef SCIP_OUTPUT
2277  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2278 #else
2279  if ( sepadata->output )
2280  {
2281  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2282  }
2283 #endif
2284  *success = FALSE;
2285  return SCIP_OKAY;
2286  }
2287 
2288  /* all other statuses except optimal or bestsollimit are invalid - (problem cannot be infeasible) */
2289  if ( status != SCIP_STATUS_OPTIMAL && status != SCIP_STATUS_BESTSOLLIMIT )
2290  {
2291  SCIPerrorMessage("Solution of subscip for CG-separation returned with invalid status %d.\n", status);
2292  return SCIP_ERROR;
2293  }
2294 
2295  /* solve some more, if a feasible solution was found */
2296  if ( status == SCIP_STATUS_BESTSOLLIMIT )
2297  {
2298  SCIPdebugMessage("Continue solving separation problem ...\n");
2299 
2300  SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", STALLNODELIMIT) );
2301  retcode = SCIPsolve(subscip);
2302  SCIP_CALL( SCIPsetLongintParam(subscip, "limits/stallnodes", -1LL) );
2303 
2304  /* errors in solving the subproblem should not kill the overall solving process;
2305  * hence, the return code is caught and a warning is printed, only in debug mode, SCIP will stop. */
2306  if ( retcode != SCIP_OKAY )
2307  {
2308 #ifndef NDEBUG
2309  SCIP_CALL( retcode );
2310 #endif
2311  SCIPwarningMessage(scip, "Error while solving subproblem in CGMIP separator; sub-SCIP terminated with code <%d>\n", retcode);
2312  *success = FALSE;
2313  return SCIP_OKAY;
2314  }
2315 
2316  status = SCIPgetStatus(subscip);
2317  assert( status != SCIP_STATUS_BESTSOLLIMIT );
2318 
2319 #ifdef SCIP_OUTPUT
2320  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2321 #else
2322  if ( sepadata->output )
2323  {
2324  SCIP_CALL( SCIPprintStatistics(subscip, NULL) );
2325  }
2326 #endif
2327 
2328  /* if the solution process was terminated */
2329  if ( status == SCIP_STATUS_TIMELIMIT || status == SCIP_STATUS_USERINTERRUPT || status == SCIP_STATUS_MEMLIMIT )
2330  {
2331  *success = FALSE;
2332  return SCIP_OKAY;
2333  }
2334 
2335  /* all other statuses except optimal or bestsollimit are invalid */
2336  if ( status != SCIP_STATUS_OPTIMAL && status != SCIP_STATUS_STALLNODELIMIT && status != SCIP_STATUS_NODELIMIT )
2337  {
2338  SCIPerrorMessage("Solution of subscip for CG-separation returned with invalid status %d.\n", status);
2339  return SCIP_ERROR;
2340  }
2341  }
2342  }
2343 
2344  return SCIP_OKAY;
2345 }
2346 
2347 
2348 /** Computes cut from the given multipliers
2349  *
2350  * Note that the cut computed here in general will not be the same as the one computed with the
2351  * sub-MIP, because of numerical differences. Here, we only combine rows whose corresponding
2352  * multiplier is positive w.r.t. the feasibility tolerance. In the sub-MIP, however, the rows are
2353  * combined in any case. This makes a difference, if the coefficients in the matrix are large and
2354  * hence yield a value that is larger than the tolerance.
2355  *
2356  * Because of the transformations we have the following:
2357  *
2358  * If variable \f$x_j\f$ was complemented, we have \f$x'_j = u_j - x_j\f$. If in the transformed
2359  * system the lower bound is used, its corresponding multiplier is \f$y^T A'_j - \lfloor y^T A'_j
2360  * \rfloor\f$, which corresponds to
2361  * \f[
2362  * y^T A'_j - \lfloor y^T A'_j \rfloor = - y^T A_j - \lfloor - y^T A_j \rfloor = - y^T A_j + \lceil y^T A_j \rceil
2363  * \f]
2364  * in the original system.
2365  *
2366  * If such a variable was at its upper bound before the transformation, it is at its lower bound
2367  * afterwards. Hence, its contribution to the cut is 0.
2368  *
2369  * Note that if the original LP-solution does not satisfy some of the rows with equality the
2370  * violation of the cut might be smaller than what is computed with the reduced sub-MIP.
2371  *
2372  * Furthermore, note that if continuous variables have been shifted, the computed violated may be
2373  * different as well, because the necessary changes in the lhs/rhs are not used here anymore.
2374  *
2375  * @todo check if cut is correct if continuous variables have been shifted.
2376  */
2377 static
2379  SCIP* scip, /**< original scip */
2380  SCIP_SEPA* sepa, /**< separator */
2381  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
2382  SCIP_SEPADATA* sepadata, /**< separator data */
2383  SCIP_SOL* sol, /**< current solution for sub-MIP */
2384  SCIP_Real* cutcoefs, /**< coefficients of the cut */
2385  SCIP_Real* cutrhs, /**< rhs of the cut */
2386  SCIP_Bool* localrowsused, /**< pointer to store whether local rows were used in summation */
2387  SCIP_Bool* localboundsused, /**< pointer to store whether local bounds were used in summation */
2388  int* cutrank, /**< pointer to store the cut rank */
2389  SCIP_Bool* success /**< whether we produced a valid cut */
2390  )
2391 {
2392  SCIP* subscip;
2393  SCIP_VAR** vars;
2394  SCIP_ROW** rows;
2395  SCIP_COL** cols;
2396  SCIP_Real val;
2397  SCIP_Real maxabsweight;
2398  int nvars;
2399  int ncols;
2400  int nrows;
2401  int i;
2402  int j;
2403 
2404  assert( scip != NULL );
2405  assert( mipdata != NULL );
2406  assert( sepadata != NULL );
2407  assert( cutcoefs != NULL );
2408  assert( cutrhs != NULL );
2409  assert( localrowsused != NULL );
2410  assert( cutrank != NULL );
2411  assert( success != NULL );
2412 
2413  subscip = mipdata->subscip;
2414  assert( subscip != NULL );
2415 
2416  /* get data */
2417  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2418  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2419  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2420  assert( nrows == (int) mipdata->nrows );
2421  assert( ncols == (int) mipdata->ncols );
2422 
2423  /* initialize */
2424  *success = TRUE;
2425  *localrowsused = FALSE;
2426  *cutrank = 0;
2427  *localboundsused = FALSE;
2428  BMSclearMemoryArray(cutcoefs, nvars);
2429  *cutrhs = 0.0;
2430 
2431  /* find maximal absolute weight */
2432  maxabsweight = 0.0;
2433  for (i = 0; i < nrows; ++i)
2434  {
2435  SCIP_ROW* row;
2436  SCIP_Real absweight = 0.0;
2437 
2438  row = rows[i];
2439  assert( row != NULL );
2440 
2441  /* skip modifiable rows and local rows, unless allowed */
2442  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
2443  {
2444  assert( mipdata->ylhs[i] == NULL && mipdata->yrhs[i] == NULL );
2445  continue;
2446  }
2447 
2448  /* get weight from solution (take larger of the values of lhs/rhs) */
2449  if ( mipdata->ylhs[i] != NULL )
2450  {
2451  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[i]);
2452 
2453  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2454  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2455 
2456  if ( SCIPisFeasPositive(scip, val) )
2457  absweight = val;
2458 
2459  assert( ! sepadata->onlyrankone || SCIProwGetOriginSepa(row) != sepa );
2460  }
2461  if ( mipdata->yrhs[i] != NULL )
2462  {
2463  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[i]);
2464 
2465  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2466  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2467 
2468  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
2469  if ( SCIPisFeasGT(scip, val, absweight) )
2470  absweight = val;
2471 
2472  assert( ! sepadata->onlyrankone || SCIProwGetOriginSepa(row) != sepa );
2473  }
2474  assert( ! SCIPisNegative(scip, absweight) );
2475 
2476  if ( absweight > maxabsweight )
2477  maxabsweight = absweight;
2478  }
2479 
2480  /* get weight from objective cuts */
2481  if ( sepadata->useobjub || sepadata->useobjlb )
2482  {
2483  SCIP_Real absweight = 0.0;
2484 
2485  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
2486 
2487  if ( mipdata->ylhs[mipdata->nrows] != NULL )
2488  {
2489  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[mipdata->nrows]);
2490  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2491 
2492  if ( SCIPisFeasPositive(scip, val) )
2493  absweight = val;
2494  }
2495  if ( mipdata->yrhs[mipdata->nrows] != NULL )
2496  {
2497  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[mipdata->nrows]);
2498  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2499 
2500  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
2501  if ( SCIPisFeasGT(scip, val, absweight) )
2502  absweight = val;
2503  }
2504 
2505  if ( absweight > maxabsweight )
2506  maxabsweight = absweight;
2507  }
2508 
2509  /* calculate the row summation */
2510  for (i = 0; i < nrows; ++i)
2511  {
2512  SCIP_ROW* row;
2513  SCIP_Real weight;
2514  SCIP_Real absweight;
2515  SCIP_Bool uselhs;
2516 
2517  row = rows[i];
2518  assert( row != NULL );
2519 
2520  /* skip modifiable rows and local rows, unless allowed */
2521  if ( SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !sepadata->allowlocal) )
2522  {
2523  assert( mipdata->ylhs[i] == NULL && mipdata->yrhs[i] == NULL );
2524  continue;
2525  }
2526 
2527  /* get weight from solution */
2528  weight = 0.0;
2529  uselhs = FALSE;
2530  if ( mipdata->ylhs[i] != NULL )
2531  {
2532  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[i]);
2533  assert( ! SCIPisFeasNegative(subscip, val) );
2534 
2535  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2536  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2537 
2538  if ( SCIPisFeasPositive(scip, val) )
2539  {
2540  uselhs = TRUE;
2541  weight = -val;
2542  }
2543  }
2544  if ( mipdata->yrhs[i] != NULL )
2545  {
2546  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[i]);
2547  assert( ! SCIPisFeasNegative(subscip, val) );
2548 
2549  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2550  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2551 
2552  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
2553  if ( SCIPisFeasGT(scip, val, REALABS(weight)) )
2554  weight = val;
2555  }
2556 
2557  /* add row if weight is nonzero and lies within range */
2558  absweight = REALABS(weight);
2559  if ( ! SCIPisSumZero(scip, weight) && absweight * MAXWEIGHTRANGE >= maxabsweight )
2560  {
2561  SCIP_COL** rowcols;
2562  SCIP_Real* rowvals;
2563 
2564  rowcols = SCIProwGetCols(row);
2565  rowvals = SCIProwGetVals(row);
2566 
2567  /* add the row coefficients to the sum */
2568  for (j = 0; j < SCIProwGetNLPNonz(row); ++j)
2569  {
2570  int idx;
2571  SCIP_VAR* var;
2572 
2573  assert( rowcols[j] != NULL );
2574  var = SCIPcolGetVar(rowcols[j]);
2575 
2576  assert( var != NULL );
2577  assert( SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN );
2578  assert( SCIPvarGetCol(var) == rowcols[j] );
2579 
2580  idx = SCIPvarGetProbindex(var);
2581  assert( 0 <= idx && idx < nvars );
2582 
2583  cutcoefs[idx] += weight * rowvals[j];
2584  }
2585 
2586  /* compute rhs */
2587  if ( uselhs )
2588  {
2589  assert( ! SCIPisInfinity(scip, -SCIProwGetLhs(row)) );
2590  val = SCIProwGetLhs(row) - SCIProwGetConstant(row);
2591  if ( SCIProwIsIntegral(row) )
2592  val = SCIPfeasCeil(scip, val); /* row is integral: round left hand side up */
2593  }
2594  else
2595  {
2596  assert( ! SCIPisInfinity(scip, SCIProwGetRhs(row)) );
2597  val = SCIProwGetRhs(row) - SCIProwGetConstant(row);
2598  if ( SCIProwIsIntegral(row) )
2599  val = SCIPfeasFloor(scip, val); /* row is integral: round right hand side down */
2600  }
2601  (*cutrhs) += weight * val;
2602 
2603  *localrowsused = *localrowsused || SCIProwIsLocal(row);
2604 
2605  if ( SCIProwGetRank(row) > *cutrank )
2606  *cutrank = SCIProwGetRank(row);
2607  }
2608  }
2609  /* add 1 to cutrank */
2610  ++(*cutrank);
2611 
2612  /* get weight from objective bounds */
2613  if ( sepadata->useobjub || sepadata->useobjlb )
2614  {
2615  SCIP_Real weight = 0.0;
2616  SCIP_Bool uselhs = FALSE;
2617  SCIP_Real absweight;
2618 
2619  assert( mipdata->ntotalrows == mipdata->nrows + 1 );
2620 
2621  if ( mipdata->ylhs[mipdata->nrows] != NULL )
2622  {
2623  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[mipdata->nrows]);
2624  assert( ! SCIPisFeasNegative(subscip, val) );
2625  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2626  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2627 
2628  if ( SCIPisFeasPositive(scip, val) )
2629  {
2630  uselhs = TRUE;
2631  weight = -val;
2632  }
2633  }
2634  if ( mipdata->yrhs[mipdata->nrows] != NULL )
2635  {
2636  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[mipdata->nrows]);
2637  assert( ! SCIPisFeasNegative(subscip, val) );
2638  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2639  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2640 
2641  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
2642  if ( SCIPisFeasGT(scip, val, REALABS(weight)) )
2643  weight = val;
2644  }
2645 
2646  /* add objective row if weight is nonzero and lies within range */
2647  absweight = REALABS(weight);
2648  if ( ! SCIPisSumZero(scip, weight) && absweight * MAXWEIGHTRANGE >= maxabsweight )
2649  {
2650  SCIP_Real obj = 0.0;
2651 
2652  /* add the objective row coefficients to the sum */
2653  for (j = 0; j < ncols; ++j)
2654  {
2655  obj = SCIPcolGetObj(cols[j]);
2656  if ( ! SCIPisZero(scip, obj) )
2657  cutcoefs[j] += weight * obj;
2658  }
2659 
2660  /* compute rhs */
2661  if ( uselhs )
2662  {
2663  val = SCIPgetLowerbound(scip);
2664  assert( ! SCIPisInfinity(scip, -val) );
2665  if ( SCIPisObjIntegral(scip) )
2666  val = SCIPfeasCeil(scip, val); /* objective is integral: round left hand side up */
2667  }
2668  else
2669  {
2670  val = SCIPgetUpperbound(scip);
2671  assert( ! SCIPisInfinity(scip, val) );
2672  if ( SCIPisObjIntegral(scip) )
2673  val = SCIPfeasFloor(scip, val); /* objective is integral: round right hand side down */
2674  }
2675  (*cutrhs) += weight * val;
2676  }
2677  }
2678 
2679  /* add upper bounds */
2680  for (j = 0; j < ncols; ++j)
2681  {
2682  assert( cols[j] != NULL );
2683  if ( mipdata->z[j] != NULL )
2684  {
2685  assert( mipdata->coltype[j] == colPresent );
2686 
2687  val = SCIPgetSolVal(subscip, sol, mipdata->z[j]);
2688  assert( ! SCIPisFeasNegative(subscip, val) );
2689 
2690  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
2691  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
2692 
2693  /* if a bound has been used */
2694  if ( SCIPisSumPositive(subscip, val) )
2695  {
2696  SCIP_VAR* var;
2697  int idx;
2698 
2699  var = SCIPcolGetVar(cols[j]);
2700 
2701  assert( var != NULL );
2702  assert( SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN );
2703  assert( SCIPvarIsIntegral(var) );
2704  assert( SCIPvarGetCol(var) == cols[j] );
2705 
2706  idx = SCIPvarGetProbindex(var);
2707  assert( 0 <= idx && idx < nvars );
2708 
2709  /* check whether variable is complemented */
2710  if ( mipdata->iscomplemented[j] )
2711  {
2712  SCIP_Real lbnd;
2713  lbnd = SCIPvarGetLbGlobal(var);
2714  assert( ! SCIPisInfinity(scip, -lbnd) );
2715  assert( SCIPisIntegral(scip, lbnd) );
2716  assert( SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPcolGetLb(cols[j])) );
2717 
2718  /* variable should not be free */
2719  assert( ! SCIPisInfinity(scip, -lbnd) || ! SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) );
2720 
2721  /* if allowed, try to use stronger local bound */
2722  if ( sepadata->allowlocal && SCIPvarGetLbLocal(var) - 0.5 > lbnd )
2723  {
2724  lbnd = SCIPvarGetLbLocal(var);
2725  assert( SCIPisIntegral(scip, lbnd) );
2726  *localboundsused = TRUE;
2727  }
2728 
2729  cutcoefs[idx] -= val;
2730  *cutrhs -= lbnd * val;
2731  }
2732  else
2733  {
2734  SCIP_Real ubnd;
2735  ubnd = SCIPvarGetUbGlobal(var);
2736  assert( ! SCIPisInfinity(scip, ubnd) );
2737  assert( SCIPisIntegral(scip, ubnd) );
2738  assert( SCIPisEQ(scip, SCIPvarGetUbLocal(var), SCIPcolGetUb(cols[j])) );
2739 
2740  /* if allowed, try to use stronger local bound */
2741  if ( sepadata->allowlocal && SCIPvarGetUbLocal(var) + 0.5 < ubnd )
2742  {
2743  ubnd = SCIPvarGetUbLocal(var);
2744  assert( SCIPisIntegral(scip, ubnd) );
2745  *localboundsused = TRUE;
2746  }
2747 
2748  cutcoefs[idx] += val;
2749  *cutrhs += ubnd * val;
2750  }
2751  }
2752  }
2753  }
2754 
2755  /* check lower bounds for integral variables */
2756  for (j = 0; j < nvars; ++j)
2757  {
2758  SCIP_VAR* var;
2759  int pos;
2760 
2761  var = vars[j];
2762  assert( var != NULL );
2763  pos = SCIPcolGetLPPos(SCIPvarGetCol(var));
2764 
2765  /* a variable may have status COLUMN, but the corresponding column may not (yet) be in the LP */
2766  if ( pos >= 0 && mipdata->coltype[pos] != colContinuous && mipdata->coltype[pos] != colConverted )
2767  {
2768  assert( pos < ncols );
2769  assert( SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN );
2770  assert( SCIPvarIsIntegral(var) );
2771 
2772  /* check whether variable is complemented */
2773  if ( mipdata->iscomplemented[pos] )
2774  {
2775  assert( ! mipdata->isshifted[pos] );
2776  /* if the variable is complemented, the multiplier for the upper bound arises from the
2777  lower bound multiplier for the transformed problem - because of the minus-sign in the
2778  transformation this yields a round-up operation. */
2779  val = SCIPfeasCeil(scip, cutcoefs[j]) - cutcoefs[j];
2780  assert( ! SCIPisFeasNegative(scip, val) );
2781 
2782  /* only if variable needs to be rounded */
2783  if ( SCIPisSumPositive(scip, val) )
2784  {
2785  SCIP_Real ubnd;
2786  ubnd = SCIPvarGetUbGlobal(var);
2787  assert( ! SCIPisInfinity(scip, ubnd) );
2788  assert( SCIPisIntegral(scip, ubnd) );
2789 
2790  /* variable should not be free */
2791  assert( ! SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)) || ! SCIPisInfinity(scip, ubnd) );
2792 
2793  /* if allowed, try to use stronger local bound */
2794  if ( sepadata->allowlocal && SCIPvarGetUbLocal(var) + 0.5 < ubnd )
2795  {
2796  ubnd = SCIPvarGetUbLocal(var);
2797  assert( SCIPisIntegral(scip, ubnd) );
2798  *localboundsused = TRUE;
2799  }
2800 
2801  /* round cut coefficients, i.e., add val to cutcoefs[j] */
2802  cutcoefs[j] = SCIPfeasCeil(scip, cutcoefs[j]);
2803 
2804  /* correct rhs */
2805  if ( ! SCIPisSumZero(scip, ubnd) )
2806  *cutrhs += ubnd * val;
2807  }
2808  }
2809  else
2810  {
2811  /* compute multiplier for lower bound: */
2812  val = cutcoefs[j] - SCIPfeasFloor(scip, cutcoefs[j]);
2813  assert( ! SCIPisFeasNegative(scip, val) );
2814 
2815  /* only if variable needs to be rounded */
2816  if ( SCIPisSumPositive(scip, val) )
2817  {
2818  SCIP_Real lbnd;
2819  lbnd = SCIPvarGetLbGlobal(var);
2820  assert( ! SCIPisInfinity(scip, -lbnd) );
2821  assert( SCIPisIntegral(scip, lbnd) );
2822 
2823  /* variable should not be free */
2824  assert( ! SCIPisInfinity(scip, -lbnd) || ! SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)) );
2825 
2826  /* if allowed, try to use stronger local bound */
2827  if ( sepadata->allowlocal && SCIPvarGetLbLocal(var) - 0.5 > lbnd )
2828  {
2829  lbnd = SCIPvarGetLbLocal(var);
2830  assert( SCIPisIntegral(scip, lbnd) );
2831  *localboundsused = TRUE;
2832  }
2833 
2834  /* round cut coefficients, i.e., subtract val from cutcoefs[j] */
2835  cutcoefs[j] = SCIPfeasFloor(scip, cutcoefs[j]);
2836 
2837  /* correct rhs */
2838  if ( ! SCIPisSumZero(scip, lbnd) )
2839  *cutrhs -= lbnd * val;
2840  }
2841  }
2842  }
2843  else
2844  {
2845  /* force coefficients of all continuous variables or of variables not in the lp to zero */
2846  assert( pos == -1 || mipdata->coltype[pos] == colContinuous || mipdata->coltype[pos] == colConverted );
2847 
2848  /* check whether all coefficients for continuous or converted variables are nonnegative */
2849  if ( pos >= 0 )
2850  {
2851  if ( SCIPisNegative(scip, cutcoefs[j]) )
2852  {
2853  *success = FALSE;
2854  break;
2855  }
2856  }
2857 
2858  cutcoefs[j] = 0.0;
2859  }
2860  }
2861 
2862  /* round rhs */
2863  *cutrhs = SCIPfeasFloor(scip, *cutrhs);
2864 
2865  return SCIP_OKAY;
2866 }
2867 
2868 
2869 /** Create CG-cut directly from solution of sub-MIP */
2870 static
2872  SCIP* scip, /**< SCIP data structure */
2873  SCIP_SEPA* sepa, /**< separator */
2874  SCIP_SEPADATA* sepadata, /**< separator data */
2875  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
2876  SCIP_SOL* sol, /**< solution of sub-MIP */
2877  SCIP_Real* cutcoefs, /**< cut coefficients */
2878  SCIP_VAR** cutvars, /**< variables appearing in cut */
2879  SCIP_Real* cutvals, /**< values of variables in cut */
2880  SCIP_Real* varsolvals, /**< solution value of variables */
2881  SCIP_Real* weights, /**< weights to compute cmir cut */
2882  int* nprevrows, /**< number of previously generated rows */
2883  SCIP_ROW** prevrows, /**< previously generated rows */
2884  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
2885  unsigned int* ngen /**< number of generated cuts */
2886  )
2887 {
2888  char name[SCIP_MAXSTRLEN];
2889  SCIP_Bool cutislocal;
2890  SCIP_Bool localrowsused;
2891  SCIP_Bool localboundsused;
2892  SCIP_Real cutrhs;
2893  SCIP_Real cutact;
2894  SCIP_Bool success;
2895  SCIP_VAR** vars;
2896  int cutrank = 0;
2897  int nvars;
2898  int k;
2899 
2900  assert( scip != NULL );
2901  assert( sepadata != NULL );
2902  assert( mipdata != NULL );
2903  assert( sol != NULL );
2904  assert( cutcoefs != NULL );
2905  assert( cutvars != NULL );
2906  assert( cutvals != NULL );
2907  assert( varsolvals != NULL );
2908  assert( weights != NULL );
2909  assert( nprevrows != NULL );
2910  assert( prevrows != NULL );
2911  assert( cutoff != NULL );
2912  assert( ngen != NULL );
2913 
2914  /* get variable data */
2915  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2916 
2917  cutrhs = 0.0;
2918  localrowsused = FALSE;
2919  localboundsused = FALSE;
2920  *cutoff = FALSE;
2921  success = TRUE;
2922 
2923  /* compute coefficients */
2924  SCIP_CALL( computeCut(scip, sepa, mipdata, sepadata, sol, cutcoefs, &cutrhs, &localrowsused, &localboundsused, &cutrank, &success) );
2925  cutislocal = localrowsused || localboundsused;
2926 
2927  /* take next solution if cut was not valid */
2928  if ( ! success )
2929  {
2930  SCIPdebugMessage("cut not valid - skipping ...\n");
2931  return SCIP_OKAY;
2932  }
2933 
2934  /* compute activity */
2935  cutact = 0.0;
2936  for (k = 0; k < nvars; ++k)
2937  cutact += cutcoefs[k] * varsolvals[k];
2938 
2939  /* the following test should be treated with care because of numerical differences - see computeCut() */
2940 #if 0
2941  {
2942  /* check for correctness of computed values */
2943  SCIP* subscip;
2944  SCIP_Real obj = 0.0;
2945  SCIP_Real val;
2946  SCIP_Bool contVarShifted = FALSE;
2947  unsigned int j;
2948  SCIP_COL** cols;
2949  int ncols;
2950 
2951  subscip = mipdata->subscip;
2952  assert( subscip != NULL );
2953 
2954  SCIP_CALL( SCIPprintSol(subscip, sol, NULL, FALSE) );
2955 
2956  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2957  for (j = 0; j < mipdata->ncols; ++j)
2958  {
2959  if ( mipdata->coltype[j] == colPresent )
2960  {
2961  int idx;
2962  assert( mipdata->alpha[j] != NULL );
2963  val = SCIPgetSolVal(subscip, sol, mipdata->alpha[j]);
2964  assert( SCIPisFeasIntegral(subscip, val) );
2965  idx = SCIPvarGetProbindex(SCIPcolGetVar(cols[j]));
2966  assert( SCIPisFeasEQ(scip, val, cutcoefs[idx]) );
2967  obj += val * SCIPvarGetObj(mipdata->alpha[j]);
2968  }
2969  else
2970  {
2971  if ( (mipdata->coltype[j] == colContinuous || mipdata->coltype[j] == colConverted) && mipdata->isshifted[j] )
2972  contVarShifted = TRUE;
2973  }
2974  }
2975  assert( mipdata->beta != NULL );
2976  val = SCIPgetSolVal(subscip, sol, mipdata->beta);
2977  assert( SCIPisFeasIntegral(subscip, val) );
2978  obj += val * SCIPvarGetObj(mipdata->beta);
2979  assert( contVarShifted || SCIPisFeasEQ(scip, obj, cutact - cutrhs) );
2980  }
2981 #endif
2982 
2983  /* if successful, convert dense cut into sparse row, and add the row as a cut */
2984  if ( SCIPisFeasGT(scip, cutact, cutrhs) )
2985  {
2986  SCIP_Real cutnorm;
2987  int cutlen;
2988 
2989  /* store the cut as sparse row, calculate activity and norm of cut */
2990  SCIP_CALL( storeCutInArrays(scip, nvars, vars, cutcoefs, varsolvals, mipdata->normtype,
2991  cutvars, cutvals, &cutlen, &cutact, &cutnorm) );
2992 
2993  SCIPdebugMessage("act=%f, rhs=%f, norm=%f, eff=%f\n", cutact, cutrhs, cutnorm, (cutact - cutrhs)/cutnorm);
2994 
2995  /* if norm is 0, the cut is trivial */
2996  if ( SCIPisPositive(scip, cutnorm) )
2997  {
2998  SCIP_Bool violated = SCIPisEfficacious(scip, (cutact - cutrhs)/cutnorm);
2999 
3000  if ( violated || (sepadata->usecutpool && ! cutislocal ) )
3001  {
3002  SCIP_ROW* cut;
3003 
3004  /* create the cut */
3005  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cgcut%d_%u", SCIPgetNLPs(scip), *ngen);
3006  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, name, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, sepadata->dynamiccuts) );
3007  SCIP_CALL( SCIPaddVarsToRow(scip, cut, cutlen, cutvars, cutvals) );
3008 
3009  /* set cut rank */
3010  SCIProwChgRank(cut, cutrank);
3011 
3012  /*SCIPdebug( SCIP_CALL( SCIPprintRow(scip, cut, NULL) ) );*/
3013 
3014  /* add cut to pool */
3015  if ( ! cutislocal )
3016  {
3017  assert( violated || sepadata->usecutpool );
3018  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
3019  }
3020 
3021  /* add cut if it is violated */
3022  if ( violated )
3023  {
3024  /* check whether cut has been found before - may happend due to projection */
3025  for (k = 0; k < *nprevrows; ++k)
3026  {
3027  SCIP_Real parval;
3028 
3029  assert( prevrows[k] != NULL );
3030  parval = SCIProwGetParallelism(cut, prevrows[k], 'e');
3031  /* exit if row is parallel to existing cut and rhs is not better */
3032  if ( SCIPisEQ(scip, parval, 1.0) && SCIPisGE(scip, cutrhs, SCIProwGetRhs(prevrows[k])) )
3033  break;
3034  }
3035 
3036  /* if cut is new */
3037  if ( k >= *nprevrows )
3038  {
3039  prevrows[*nprevrows] = cut;
3040  ++(*nprevrows);
3041 
3042  SCIPdebugMessage(" -> CG-cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
3043  name, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
3044  SCIPgetCutEfficacy(scip, NULL, cut),
3045  SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
3046  SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
3047 #ifdef SCIP_DEBUG
3048  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3049 #else
3050  if ( sepadata->output )
3051  {
3052  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3053  }
3054 #endif
3055  SCIP_CALL( SCIPaddCut(scip, NULL, cut, FALSE, cutoff) );
3056  ++(*ngen);
3057  }
3058  else
3059  {
3060  SCIPdebugMessage("Cut already exists.\n");
3061  /* release the row */
3062  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3063  }
3064  }
3065  else
3066  {
3067  /* release the row */
3068  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3069  }
3070  }
3071  }
3072  }
3073 
3074  return SCIP_OKAY;
3075 }
3076 
3077 
3078 /** create CG-cut via CMIR-function */
3079 static
3081  SCIP* scip, /**< SCIP data structure */
3082  SCIP_SEPA* sepa, /**< separator */
3083  SCIP_SEPADATA* sepadata, /**< separator data */
3084  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
3085  SCIP_SOL* sol, /**< solution of sub-MIP */
3086  SCIP_Real* cutcoefs, /**< cut coefficients */
3087  SCIP_VAR** cutvars, /**< variables appearing in cut */
3088  SCIP_Real* cutvals, /**< values of variables in cut */
3089  SCIP_Real* varsolvals, /**< solution value of variables */
3090  SCIP_Real* weights, /**< weights to compute cmir cut */
3091  int* boundsfortrans, /**< bounds for cmir function of NULL */
3092  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds for cmir function or NULL */
3093  int* nprevrows, /**< number of previously generated rows */
3094  SCIP_ROW** prevrows, /**< previously generated rows */
3095  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
3096  unsigned int* ngen /**< number of generated cuts */
3097  )
3098 {
3099  char name[SCIP_MAXSTRLEN];
3100  SCIP_Longint maxdnom;
3101  SCIP_Bool cutislocal;
3102  SCIP_Real maxscale;
3103  SCIP_Real cutrhs;
3104  SCIP_Real cutact;
3105  SCIP_Bool success;
3106  SCIP_ROW** rows;
3107  SCIP_VAR** vars;
3108  SCIP* subscip;
3109  int nrows;
3110  int nvars;
3111  int k;
3112  int cutrank;
3113 
3114  assert( scip != NULL );
3115  assert( sepadata != NULL );
3116  assert( mipdata != NULL );
3117  assert( sol != NULL );
3118  assert( cutcoefs != NULL );
3119  assert( cutvars != NULL );
3120  assert( cutvals != NULL );
3121  assert( varsolvals != NULL );
3122  assert( weights != NULL );
3123  assert( nprevrows != NULL );
3124  assert( prevrows != NULL );
3125  assert( cutoff != NULL );
3126  assert( ngen != NULL );
3127 
3128  *cutoff = FALSE;
3129  subscip = mipdata->subscip;
3130  assert( subscip != NULL );
3131 
3132  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
3133  assert( nrows > 0 );
3134  assert( (int) mipdata->nrows == nrows );
3135 
3136  /* @todo more advanced settings - compare sepa_gomory.c */
3137  maxdnom = (SCIP_Longint) sepadata->cutcoefbnd+1;
3138  maxscale = 10000.0;
3139 
3140  /* get variable data */
3141  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
3142 
3143  /* generate weights */
3144  for (k = 0; k < nrows; ++k)
3145  {
3146  SCIP_Real val;
3147 
3148  weights[k] = 0;
3149  if ( mipdata->ylhs[k] != NULL )
3150  {
3151  assert( !SCIProwIsModifiable(rows[k]) && (!SCIProwIsLocal(rows[k]) || sepadata->allowlocal) );
3152 
3153  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[k]);
3154  assert( ! SCIPisFeasNegative(subscip, val) );
3155 
3156  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
3157  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
3158 
3159  if ( SCIPisFeasPositive(subscip, val) )
3160  weights[k] = -val;
3161  }
3162  if ( mipdata->yrhs[k] != NULL )
3163  {
3164  assert( !SCIProwIsModifiable(rows[k]) && (!SCIProwIsLocal(rows[k]) || sepadata->allowlocal) );
3165 
3166  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[k]);
3167  assert( ! SCIPisFeasNegative(subscip, val) );
3168 
3169  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
3170  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
3171 
3172  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
3173  if ( SCIPisFeasGT(scip, val, ABS(weights[k])) )
3174  weights[k] = val;
3175  }
3176  }
3177 
3178  /* set up data for bounds to use */
3179  if ( sepadata->cmirownbounds )
3180  {
3181  int typefortrans;
3182 
3183  assert( boundsfortrans != NULL );
3184  assert( boundtypesfortrans != NULL );
3185 
3186  if ( sepadata->allowlocal )
3187  typefortrans = -2;
3188  else
3189  typefortrans = -1;
3190 
3191  /* check all variables */
3192  for (k = 0; k < nvars; ++k)
3193  {
3194  int pos;
3195  SCIP_VAR* var;
3196 
3197  var = vars[k];
3198  assert( var != NULL );
3199  pos = SCIPcolGetLPPos(SCIPvarGetCol(var));
3200 
3201  if ( pos < 0 )
3202  continue;
3203 
3204  assert( pos < (int) mipdata->ncols );
3205  assert( SCIPvarGetStatus(var) == SCIP_VARSTATUS_COLUMN );
3206 
3207  boundsfortrans[k] = typefortrans;
3208  boundtypesfortrans[k] = SCIP_BOUNDTYPE_LOWER;
3209 
3210  if ( mipdata->coltype[pos] == colContinuous || mipdata->coltype[pos] == colConverted )
3211  {
3212  assert( SCIPvarIsIntegral(var) || mipdata->coltype[pos] != colContinuous );
3213  continue;
3214  }
3215 
3216  /* check upper bound */
3217  if ( mipdata->z[pos] != NULL && SCIPisSumPositive(subscip, SCIPgetSolVal(subscip, sol, mipdata->z[pos])) )
3218  {
3219  /* check whether variable is complemented */
3220  if ( ! mipdata->iscomplemented[pos] )
3221  boundtypesfortrans[k] = SCIP_BOUNDTYPE_UPPER;
3222  /* otherwise use lower bound */
3223  }
3224  else
3225  {
3226  /* check whether variable is complemented */
3227  if ( mipdata->iscomplemented[pos] )
3228  boundtypesfortrans[k] = SCIP_BOUNDTYPE_UPPER;
3229  /* otherwise use lower bound */
3230  }
3231  }
3232  }
3233 
3234  /* create a MIR cut using the above calculated weights */
3235  cutact = -1.0;
3236  cutrhs = -1.0;
3237  SCIP_CALL( SCIPcalcMIR(scip, NULL, BOUNDSWITCH, USEVBDS, sepadata->allowlocal, FIXINTEGRALRHS, boundsfortrans, boundtypesfortrans,
3238  (int) MAXAGGRLEN(nvars), MAXWEIGHTRANGE, MINFRAC, MAXFRAC,
3239  weights, NULL, 1.0, NULL, NULL, cutcoefs, &cutrhs, &cutact, &success, &cutislocal, &cutrank) );
3240  assert( sepadata->allowlocal || !cutislocal );
3241  SCIPdebugMessage("CMIR: success = %u, cut is%sviolated (cutact: %g, cutrhs: %g)\n", success,
3242  SCIPisFeasGT(scip, cutact, cutrhs) ? " " : " not ", cutact, cutrhs);
3243 
3244  /* if successful, convert dense cut into sparse row, and add the row as a cut */
3245  if ( success && SCIPisFeasGT(scip, cutact, cutrhs) )
3246  {
3247  SCIP_Real cutnorm;
3248  int cutlen;
3249 
3250  /* store the cut as sparse row, calculate activity and norm of cut */
3251  SCIP_CALL( storeCutInArrays(scip, nvars, vars, cutcoefs, varsolvals, mipdata->normtype,
3252  cutvars, cutvals, &cutlen, &cutact, &cutnorm) );
3253 
3254  /* only proceed if norm is positive - otherwise the cut is trivial */
3255  if ( SCIPisPositive(scip, cutnorm) )
3256  {
3257  SCIP_Bool violated;
3258 
3259  violated = SCIPisEfficacious(scip, (cutact - cutrhs)/cutnorm);
3260 
3261  /* only if the cut if violated - if it is not violated we might store non-local cuts in the pool */
3262  if ( violated || ( sepadata->usecutpool && ! cutislocal) )
3263  {
3264  SCIP_ROW* cut;
3265 
3266  /* create the cut */
3267  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cgcut%d_%u", SCIPgetNLPs(scip), *ngen);
3268  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, name, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, sepadata->dynamiccuts) );
3269  SCIP_CALL( SCIPaddVarsToRow(scip, cut, cutlen, cutvars, cutvals) );
3270  assert( success );
3271 
3272  /* set cut rank */
3273  SCIProwChgRank(cut, cutrank);
3274 
3275 #ifdef SCIP_DEBUG
3276  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3277 #else
3278  if ( sepadata->output )
3279  {
3280  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3281  }
3282 #endif
3283 
3284  /* try to scale the cut to integral values */
3285  SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
3286  maxdnom, maxscale, MAKECONTINTEGRAL, &success) );
3287 
3288  /* if the cut could be made integral */
3289  if ( success )
3290  {
3291  /* add cut to pool */
3292  if ( ! cutislocal )
3293  {
3294  assert( violated || sepadata->usecutpool );
3295  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
3296  }
3297 
3298  if ( ! SCIPisCutEfficacious(scip, NULL, cut) )
3299  {
3300  SCIPdebugMessage(" -> CG-cut <%s> no longer efficacious: act=%f, rhs=%f, norm=%f, eff=%f\n",
3301  name, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
3302  SCIPgetCutEfficacy(scip, NULL, cut));
3303 
3304  /* release the row */
3305  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3306  }
3307  else
3308  {
3309  /* check whether cut has been found before - may happend due to projection */
3310  for (k = 0; k < *nprevrows; ++k)
3311  {
3312  SCIP_Real parval;
3313 
3314  assert( prevrows[k] != NULL );
3315  parval = SCIProwGetParallelism(cut, prevrows[k], 'e');
3316  /* exit if row is parallel to existing cut and rhs is not better */
3317  if ( SCIPisEQ(scip, parval, 1.0) && SCIPisGE(scip, cutrhs, SCIProwGetRhs(prevrows[k])) )
3318  break;
3319  }
3320 
3321  /* if cut is new */
3322  if ( k >= *nprevrows )
3323  {
3324  prevrows[*nprevrows] = cut;
3325  ++(*nprevrows);
3326 
3327  SCIPdebugMessage(" -> CG-cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, rank=%d, min=%f, max=%f (range=%f)\n",
3328  name, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
3329  SCIPgetCutEfficacy(scip, NULL, cut), SCIProwGetRank(cut),
3330  SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
3331  SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
3332 #ifdef SCIP_OUTPUT
3333  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3334 #else
3335  if ( sepadata->output )
3336  {
3337  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3338  }
3339 #endif
3340  SCIP_CALL( SCIPaddCut(scip, NULL, cut, FALSE, cutoff) );
3341  ++(*ngen);
3342  }
3343  else
3344  {
3345  SCIPdebugMessage("Cut already exists.\n");
3346  /* release the row */
3347  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3348  }
3349  }
3350  }
3351  else
3352  {
3353  SCIPdebugMessage(" -> CG-cut <%s> could not be scaled to integral coefficients: act=%f, rhs=%f, norm=%f, eff=%f\n",
3354  name, cutact, cutrhs, cutnorm, SCIPgetCutEfficacy(scip, NULL, cut));
3355 
3356  /* release the row */
3357  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3358  }
3359  }
3360  }
3361  }
3362 
3363  return SCIP_OKAY;
3364 }
3365 
3366 
3367 /** create CG-cut via strong-CG-function */
3368 static
3370  SCIP* scip, /**< SCIP data structure */
3371  SCIP_SEPA* sepa, /**< separator */
3372  SCIP_SEPADATA* sepadata, /**< separator data */
3373  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
3374  SCIP_SOL* sol, /**< solution of sub-MIP */
3375  SCIP_Real* cutcoefs, /**< cut coefficients */
3376  SCIP_VAR** cutvars, /**< variables appearing in cut */
3377  SCIP_Real* cutvals, /**< values of variables in cut */
3378  SCIP_Real* varsolvals, /**< solution value of variables */
3379  SCIP_Real* weights, /**< weights to compute cmir cut */
3380  int* nprevrows, /**< number of previously generated rows */
3381  SCIP_ROW** prevrows, /**< previously generated rows */
3382  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
3383  unsigned int* ngen /**< number of generated cuts */
3384  )
3385 {
3386  char name[SCIP_MAXSTRLEN];
3387  SCIP_Longint maxdnom;
3388  SCIP_Bool cutislocal;
3389  SCIP_Real maxscale;
3390  SCIP_Real cutrhs;
3391  SCIP_Real cutact;
3392  SCIP_Bool success;
3393  SCIP_ROW** rows;
3394  SCIP_VAR** vars;
3395  SCIP* subscip;
3396  int nrows;
3397  int nvars;
3398  int k;
3399  int cutrank;
3400 
3401  assert( scip != NULL );
3402  assert( sepadata != NULL );
3403  assert( mipdata != NULL );
3404  assert( sol != NULL );
3405  assert( cutcoefs != NULL );
3406  assert( cutvars != NULL );
3407  assert( cutvals != NULL );
3408  assert( varsolvals != NULL );
3409  assert( weights != NULL );
3410  assert( nprevrows != NULL );
3411  assert( prevrows != NULL );
3412  assert( cutoff != NULL );
3413  assert( ngen != NULL );
3414 
3415  *cutoff = FALSE;
3416  subscip = mipdata->subscip;
3417  assert( subscip != NULL );
3418 
3419  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
3420  assert( nrows > 0 );
3421  assert( (int) mipdata->nrows == nrows );
3422 
3423  /* @todo more advanced settings - compare sepa_gomory.c */
3424  maxdnom = (SCIP_Longint) sepadata->cutcoefbnd + 1;
3425  maxscale = 10000.0;
3426 
3427  /* get variable data */
3428  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
3429 
3430  /* generate weights */
3431  for (k = 0; k < nrows; ++k)
3432  {
3433  SCIP_Real val;
3434 
3435  weights[k] = 0;
3436  if ( mipdata->ylhs[k] != NULL )
3437  {
3438  assert( !SCIProwIsModifiable(rows[k]) && (!SCIProwIsLocal(rows[k]) || sepadata->allowlocal) );
3439 
3440  val = SCIPgetSolVal(subscip, sol, mipdata->ylhs[k]);
3441  assert( ! SCIPisFeasNegative(subscip, val) );
3442 
3443  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
3444  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
3445 
3446  if ( SCIPisFeasPositive(subscip, val) )
3447  weights[k] = -val;
3448  }
3449  if ( mipdata->yrhs[k] != NULL )
3450  {
3451  assert( !SCIProwIsModifiable(rows[k]) && (!SCIProwIsLocal(rows[k]) || sepadata->allowlocal) );
3452 
3453  val = SCIPgetSolVal(subscip, sol, mipdata->yrhs[k]);
3454  assert( ! SCIPisFeasNegative(subscip, val) );
3455 
3456  assert( sepadata->skipmultbounds || SCIPisFeasLT(subscip, val, 1.0) );
3457  val = SCIPfrac(scip, val); /* take fractional value if variable has no upper bounds */
3458 
3459  /* in a suboptimal solution both values may be positive - take the one with larger absolute value */
3460  if ( SCIPisFeasGT(scip, val, ABS(weights[k])) )
3461  weights[k] = val;
3462  }
3463  }
3464 
3465  /* create a strong CG cut out of the weighted LP rows using the B^-1 row as weights */
3466  cutact = -1.0;
3467  cutrhs = -1.0;
3468  SCIP_CALL( SCIPcalcStrongCG(scip, BOUNDSWITCH, USEVBDS, sepadata->allowlocal, (int) MAXAGGRLEN(nvars), MAXWEIGHTRANGE, MINFRAC, MAXFRAC,
3469  weights, 1.0, cutcoefs, &cutrhs, &cutact, &success, &cutislocal, &cutrank) );
3470  assert( sepadata->allowlocal || !cutislocal );
3471  SCIPdebugMessage("Strong-CG: success = %u, cut is%sviolated (cutact: %g, cutrhs: %g)\n", success,
3472  SCIPisFeasGT(scip, cutact, cutrhs) ? " " : " not ", cutact, cutrhs);
3473 
3474  /* if successful, convert dense cut into sparse row, and add the row as a cut */
3475  if ( success && SCIPisFeasGT(scip, cutact, cutrhs) )
3476  {
3477  SCIP_Real cutnorm;
3478  int cutlen;
3479 
3480  /* store the cut as sparse row, calculate activity and norm of cut */
3481  SCIP_CALL( storeCutInArrays(scip, nvars, vars, cutcoefs, varsolvals, mipdata->normtype,
3482  cutvars, cutvals, &cutlen, &cutact, &cutnorm) );
3483 
3484  /* only proceed if norm is positive - otherwise the cut is trivial */
3485  if ( SCIPisPositive(scip, cutnorm) )
3486  {
3487  SCIP_Bool violated;
3488 
3489  violated = SCIPisEfficacious(scip, (cutact - cutrhs)/cutnorm);
3490 
3491  /* only if the cut if violated - if it is not violated we might store non-local cuts in the pool */
3492  if ( violated || ( sepadata->usecutpool && ! cutislocal) )
3493  {
3494  SCIP_ROW* cut;
3495 
3496  /* create the cut */
3497  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "cgcut%d_%u", SCIPgetNLPs(scip), *ngen);
3498  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, name, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE, sepadata->dynamiccuts) );
3499  SCIP_CALL( SCIPaddVarsToRow(scip, cut, cutlen, cutvars, cutvals) );
3500  SCIProwChgRank(cut, cutrank);
3501 
3502  assert( success );
3503 
3504 #ifdef SCIP_DEBUG
3505  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3506 #else
3507  if ( sepadata->output )
3508  {
3509  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3510  }
3511 #endif
3512 
3513  /* try to scale the cut to integral values */
3514  SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
3515  maxdnom, maxscale, MAKECONTINTEGRAL, &success) );
3516 
3517  /* if the cut could be made integral */
3518  if ( success )
3519  {
3520  /* add cut to pool */
3521  if ( ! cutislocal )
3522  {
3523  assert( violated || sepadata->usecutpool );
3524  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
3525  }
3526 
3527  if ( ! SCIPisCutEfficacious(scip, NULL, cut) )
3528  {
3529  SCIPdebugMessage(" -> CG-cut <%s> no longer efficacious: act=%f, rhs=%f, norm=%f, eff=%f rank=%d\n",
3530  name, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
3531  SCIPgetCutEfficacy(scip, NULL, cut), SCIProwGetRank(cut));
3532 
3533  /* release the row */
3534  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3535  }
3536  else
3537  {
3538  /* check whether cut has been found before - may happend due to projection */
3539  for (k = 0; k < *nprevrows; ++k)
3540  {
3541  SCIP_Real parval;
3542 
3543  assert( prevrows[k] != NULL );
3544  parval = SCIProwGetParallelism(cut, prevrows[k], 'e');
3545  /* exit if row is parallel to existing cut and rhs is not better */
3546  if ( SCIPisEQ(scip, parval, 1.0) && SCIPisGE(scip, cutrhs, SCIProwGetRhs(prevrows[k])) )
3547  break;
3548  }
3549 
3550  /* if cut is new */
3551  if ( k >= *nprevrows )
3552  {
3553  prevrows[*nprevrows] = cut;
3554  ++(*nprevrows);
3555 
3556  SCIPdebugMessage(" -> CG-cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
3557  name, SCIPgetRowLPActivity(scip, cut), SCIProwGetRhs(cut), SCIProwGetNorm(cut),
3558  SCIPgetCutEfficacy(scip, NULL, cut),
3559  SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
3560  SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
3561 
3562 #ifdef SCIP_OUTPUT
3563  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3564 #else
3565  if ( sepadata->output )
3566  {
3567  SCIP_CALL( SCIPprintRow(scip, cut, NULL) );
3568  }
3569 #endif
3570 
3571  SCIP_CALL( SCIPaddCut(scip, NULL, cut, FALSE, cutoff) );
3572  ++(*ngen);
3573  }
3574  else
3575  {
3576  SCIPdebugMessage("Cut already exists.\n");
3577  /* release the row */
3578  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3579  }
3580  }
3581  }
3582  else
3583  {
3584  SCIPdebugMessage(" -> CG-cut <%s> could not be scaled to integral coefficients: act=%f, rhs=%f, norm=%f, eff=%f\n",
3585  name, cutact, cutrhs, cutnorm, SCIPgetCutEfficacy(scip, NULL, cut));
3586 
3587  /* release the row */
3588  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
3589  }
3590  }
3591  }
3592  }
3593 
3594  return SCIP_OKAY;
3595 }
3596 
3597 
3598 /** Create CG-cuts from solutions of sub-MIP */
3599 static
3601  SCIP* scip, /**< SCIP data structure */
3602  SCIP_SEPA* sepa, /**< separator */
3603  SCIP_SEPADATA* sepadata, /**< separator data */
3604  CGMIP_MIPDATA* mipdata, /**< data for sub-MIP */
3605  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
3606  unsigned int* ngen /**< number of generated cuts */
3607  )
3608 {
3609  SCIP_BOUNDTYPE* boundtypesfortrans;
3610  SCIP_STAGE stage;
3611  SCIP_Real* varsolvals;
3612  SCIP_Real* weights;
3613  SCIP_VAR** cutvars;
3614  SCIP_Real* cutvals;
3615  SCIP_Real* cutcoefs;
3616  SCIP_ROW** prevrows;
3617  SCIP_SOL** sols;
3618  SCIP_VAR** vars;
3619  SCIP* subscip;
3620  int* boundsfortrans;
3621  int nprevrows;
3622  int ntotalrows;
3623  int nsols;
3624  int nvars;
3625  int k;
3626  int s;
3627 
3628  assert( scip != NULL );
3629  assert( sepadata != NULL );
3630  assert( mipdata != NULL );
3631  assert( cutoff != NULL );
3632  assert( ngen != NULL );
3633 
3634  subscip = mipdata->subscip;
3635  assert( subscip != NULL );
3636 
3637  *cutoff = FALSE;
3638  *ngen = 0;
3639 
3640  /* check if solving was successful and get solutions */
3641  stage = SCIPgetStage(subscip);
3642  if ( stage == SCIP_STAGE_SOLVING || stage == SCIP_STAGE_SOLVED )
3643  nsols = SCIPgetNSols(subscip);
3644  else
3645  nsols = 0;
3646 
3647  /* only if solutions have been found */
3648  if ( nsols == 0 )
3649  return SCIP_OKAY;
3650 
3651  SCIPdebugMessage("Generating CG-cuts from %d sols (cmir: %u, strong-cg: %u) ...\n", nsols, sepadata->usecmir, sepadata->usestrongcg);
3652 
3653  sols = SCIPgetSols(subscip);
3654 
3655  /* get variable data */
3656  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
3657 
3658  /* allocate temporary memory */
3659  assert(mipdata->ntotalrows <= INT_MAX);
3660  ntotalrows = (int)mipdata->ntotalrows;
3661  assert( ntotalrows >= SCIPgetNLPRows(scip) && ntotalrows <= SCIPgetNLPRows(scip) + 1 );
3662  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
3663  SCIP_CALL( SCIPallocBufferArray(scip, &varsolvals, nvars) );
3664  SCIP_CALL( SCIPallocBufferArray(scip, &cutvars, nvars) );
3665  SCIP_CALL( SCIPallocBufferArray(scip, &cutvals, nvars) );
3666  SCIP_CALL( SCIPallocBufferArray(scip, &weights, ntotalrows) );
3667  SCIP_CALL( SCIPallocBufferArray(scip, &prevrows, 2 * nsols) );
3668 
3669  /* prepare arrays for bound information, if requested */
3670  if ( sepadata->usecmir && sepadata->cmirownbounds )
3671  {
3672  SCIP_CALL( SCIPallocBufferArray(scip, &boundsfortrans, nvars) );
3673  SCIP_CALL( SCIPallocBufferArray(scip, &boundtypesfortrans, nvars) );
3674  }
3675  else
3676  {
3677  boundsfortrans = NULL;
3678  boundtypesfortrans = NULL;
3679  }
3680 
3681  /* get solution values */
3682  for (k = 0; k < nvars; ++k)
3683  {
3684  if ( SCIPvarGetStatus(vars[k]) == SCIP_VARSTATUS_COLUMN )
3685  varsolvals[k] = SCIPvarGetLPSol(vars[k]);
3686  else
3687  varsolvals[k] = 0.0;
3688  }
3689 
3690  /* loop through solutions found */
3691  nprevrows = 0;
3692  for (s = 0; s < nsols; ++s)
3693  {
3694  SCIP_SOL* sol;
3695  sol = sols[s];
3696 
3697  /* generate cuts by the C-MIR and/or Strong-CG functions */
3698  if ( sepadata->usecmir )
3699  {
3700  SCIP_CALL( createCGCutCMIR(scip, sepa, sepadata, mipdata, sol, cutcoefs, cutvars, cutvals, varsolvals, weights,
3701  boundsfortrans, boundtypesfortrans, &nprevrows, prevrows, cutoff, ngen) );
3702  }
3703 
3704  if ( sepadata->usestrongcg )
3705  {
3706  SCIP_CALL( createCGCutStrongCG(scip, sepa, sepadata, mipdata, sol, cutcoefs, cutvars, cutvals, varsolvals, weights,
3707  &nprevrows, prevrows, cutoff, ngen) );
3708  }
3709 
3710  if ( ! sepadata->usecmir && ! sepadata->usestrongcg )
3711  {
3712  SCIP_CALL( createCGCutDirect(scip, sepa, sepadata, mipdata, sol, cutcoefs, cutvars, cutvals, varsolvals, weights,
3713  &nprevrows, prevrows, cutoff, ngen) );
3714  }
3715  }
3716  assert( nprevrows <= 2 * nsols );
3717  assert( sepadata->usecmir || nprevrows <= nsols );
3718  assert( sepadata->usestrongcg || nprevrows <= nsols );
3719 
3720  /* release rows */
3721  for (k = 0; k < nprevrows; ++k)
3722  {
3723  SCIP_CALL( SCIPreleaseRow(scip, &(prevrows[k])) );
3724  }
3725 
3726  /* free temporary memory */
3727  SCIPfreeBufferArrayNull(scip, &boundsfortrans);
3728  SCIPfreeBufferArrayNull(scip, &boundtypesfortrans);
3729 
3730  SCIPfreeBufferArray(scip, &prevrows);
3731  SCIPfreeBufferArray(scip, &weights);
3732  SCIPfreeBufferArray(scip, &cutvals);
3733  SCIPfreeBufferArray(scip, &cutvars);
3734  SCIPfreeBufferArray(scip, &varsolvals);
3735  SCIPfreeBufferArray(scip, &cutcoefs);
3736 
3737  return SCIP_OKAY;
3738 }
3739 
3740 
3741 /** frees "subscip" data */
3742 static
3744  SCIP* scip, /**< SCIP data structure */
3745  SCIP_SEPA* sepa, /**< separator data */
3746  CGMIP_MIPDATA* mipdata /**< data for sub-MIP */
3747  )
3748 {
3749  SCIP_SEPADATA* sepadata;
3750  unsigned int i, j;
3751  SCIP* subscip;
3752 
3753  assert( scip != NULL );
3754  assert( sepa != NULL );
3755  assert( mipdata != NULL );
3756 
3757  /* free separator data */
3758  sepadata = SCIPsepaGetData(sepa);
3759  assert( sepadata != NULL );
3760 
3761  SCIPdebugMessage("Freeing subscip ...\n");
3762 
3763  subscip = mipdata->subscip;
3764  assert( subscip != NULL );
3765 
3766  for (j = 0; j < mipdata->ncols; ++j)
3767  {
3768  if ( mipdata->coltype[j] == colPresent )
3769  {
3770  assert( mipdata->alpha[j] != NULL );
3771  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->alpha[j])) );
3772  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->fracalpha[j])) );
3773  }
3774  }
3775  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->beta)) );
3776  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->fracbeta)) );
3777 
3778  for (i = 0; i < mipdata->nrows; ++i)
3779  {
3780  if ( mipdata->ylhs[i] != NULL )
3781  {
3782  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->ylhs[i])) );
3783  }
3784  if ( mipdata->yrhs[i] != NULL )
3785  {
3786  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->yrhs[i])) );
3787  }
3788  }
3789 
3790  if ( sepadata->useobjub || sepadata->useobjlb )
3791  {
3792  if ( mipdata->yrhs[mipdata->nrows] )
3793  {
3794  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->yrhs[mipdata->nrows])) );
3795  }
3796  if ( mipdata->ylhs[mipdata->nrows] )
3797  {
3798  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->ylhs[mipdata->nrows])) );
3799  }
3800  }
3801 
3802  for (j = 0; j < mipdata->ncols; ++j)
3803  {
3804  if ( mipdata->z[j] != NULL )
3805  {
3806  SCIP_CALL( SCIPreleaseVar(subscip, &(mipdata->z[j])) );
3807  }
3808  }
3809 
3810  if ( mipdata->subscip != NULL )
3811  {
3812  SCIP_CALL( SCIPfree(&(mipdata->subscip)) );
3813  }
3814 
3815  SCIPfreeBlockMemoryArray(scip, &(mipdata->z), 2*mipdata->ncols);
3816  SCIPfreeBlockMemoryArray(scip, &(mipdata->yrhs), mipdata->ntotalrows);
3817  SCIPfreeBlockMemoryArray(scip, &(mipdata->ylhs), mipdata->ntotalrows);
3818  SCIPfreeBlockMemoryArray(scip, &(mipdata->isshifted), mipdata->ncols);
3819  SCIPfreeBlockMemoryArray(scip, &(mipdata->iscomplemented), mipdata->ncols);
3820  SCIPfreeBlockMemoryArray(scip, &(mipdata->coltype), mipdata->ncols);
3821  SCIPfreeBlockMemoryArray(scip, &(mipdata->fracalpha), mipdata->ncols);
3822  SCIPfreeBlockMemoryArray(scip, &(mipdata->alpha), mipdata->ncols);
3823 
3824  return SCIP_OKAY;
3825 }
3826 
3827 
3828 /*
3829  * Callback methods
3830  */
3831 
3832 /** copy method for separator plugins (called when SCIP copies plugins) */
3833 static
3834 SCIP_DECL_SEPACOPY(sepaCopyCGMIP)
3835 { /*lint --e{715}*/
3836  assert( scip != NULL );
3837  assert( sepa != NULL );
3838  assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
3839 
3840  /* call inclusion method of constraint handler */
3842 
3843  return SCIP_OKAY;
3844 }
3845 
3846 
3847 /** destructor of separator to free user data (called when SCIP is exiting) */
3848 static
3849 SCIP_DECL_SEPAFREE(sepaFreeCGMIP)
3850 { /*lint --e{715}*/
3851  SCIP_SEPADATA* sepadata;
3852 
3853  assert( scip != NULL );
3854  assert( sepa != NULL );
3855  assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
3856 
3857  /* free separator data */
3858  sepadata = SCIPsepaGetData(sepa);
3859  assert( sepadata != NULL );
3860 
3861  SCIPfreeMemory(scip, &sepadata);
3862 
3863  SCIPsepaSetData(sepa, NULL);
3864 
3865  return SCIP_OKAY;
3866 }
3867 
3868 
3869 /** LP solution separation method of separator */
3870 static
3871 SCIP_DECL_SEPAEXECLP(sepaExeclpCGMIP)
3872 { /*lint --e{715}*/
3873  SCIP_SEPADATA* sepadata;
3874  CGMIP_MIPDATA* mipdata;
3875 
3876  int depth;
3877  int ncalls;
3878  int ncols;
3879  int nrows;
3880  unsigned int ngen;
3881  SCIP_Bool success;
3882  SCIP_Bool cutoff = FALSE;
3883 
3884  assert( scip != NULL );
3885  assert( sepa != NULL );
3886  assert( strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0 );
3887  assert( result != NULL );
3888 
3889  *result = SCIP_DIDNOTRUN;
3890  ngen = 0;
3891 
3892  sepadata = SCIPsepaGetData(sepa);
3893  assert(sepadata != NULL);
3894 
3895  depth = SCIPgetDepth(scip);
3896 
3897  /* only call separator, if we are not close to terminating */
3898  if ( SCIPisStopped(scip) )
3899  return SCIP_OKAY;
3900 
3901  /* only call separator up to a maximum depth */
3902  if ( sepadata->maxdepth >= 0 && depth > sepadata->maxdepth )
3903  return SCIP_OKAY;
3904 
3905  /* only call separator a given number of times at each node */
3906  ncalls = SCIPsepaGetNCallsAtNode(sepa);
3907  if ( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
3908  || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
3909  return SCIP_OKAY;
3910 
3911  /* only call separator, if an optimal LP solution is at hand */
3913  return SCIP_OKAY;
3914 
3915  /* skip separation if there are continuous variables, but only integers required */
3916  if ( SCIPgetNContVars(scip) > 0 && sepadata->onlyintvars )
3917  return SCIP_OKAY;
3918 
3919  /* only call separator, if there are fractional variables */
3920  if ( SCIPgetNLPBranchCands(scip) == 0 )
3921  return SCIP_OKAY;
3922 
3923  /* check for parameters */
3924  if ( ( sepadata->useobjub || sepadata->useobjlb ) && ( sepadata->usecmir || sepadata->usestrongcg ) )
3925  {
3927  "WARNING - sepa_cgmip: Using objective function bounds and CMIR or Strong-CG functions is useless. Turning off usage of objective function bounds.\n");
3928  SCIP_CALL( SCIPsetBoolParam(scip, "separating/cgmip/useobjub", FALSE) );
3929  SCIP_CALL( SCIPsetBoolParam(scip, "separating/cgmip/useobjlb", FALSE) );
3930  }
3931 
3932  /* get LP data */
3933  ncols = SCIPgetNLPCols(scip);
3934  nrows = SCIPgetNLPRows(scip);
3935  if ( ncols <= NCOLSTOOSMALL || nrows <= NROWSTOOSMALL )
3936  return SCIP_OKAY;
3937 
3938  /* determine whether we should run the separation based on a decision tree */
3939  if ( sepadata->decisiontree )
3940  {
3941  SCIP_Bool separate;
3942  SCIP_Real firstlptime;
3943 
3944  separate = FALSE;
3945  firstlptime = SCIPgetFirstLPTime(scip);
3946 
3947  if ( nrows <= 136 && firstlptime <= 0.05 && ncols <= 143 )
3948  separate = TRUE;
3949  else if ( nrows <= 136 && 0.05 < firstlptime && firstlptime <= 0.15 && ncols <= 143 )
3950  separate = TRUE;
3951  else if ( 136 < nrows && nrows <= 332 && ncols <= 143 )
3952  separate = TRUE;
3953  else if ( 136 < nrows && nrows <= 332 && 655 < ncols && ncols <= 1290 )
3954  separate = TRUE;
3955  else if ( 333 < nrows && nrows <= 874 && 0.15 < firstlptime && firstlptime <= 0.25 && 2614 < ncols && ncols <= 5141 )
3956  separate = TRUE;
3957  else if ( 875 < nrows && nrows <= 1676 && firstlptime <= 0.05 && 143 < ncols && ncols <= 265 )
3958  separate = TRUE;
3959  else if ( 875 < nrows && nrows <= 1676 && firstlptime <= 0.05 && 265 < ncols && ncols <= 654 )
3960  separate = TRUE;
3961  else if ( 875 < nrows && nrows <= 1676 && 0.05 < firstlptime && firstlptime <= 0.15 )
3962  separate = TRUE;
3963  else if ( 875 < nrows && nrows <= 1676 && 0.15 < firstlptime && firstlptime <= 0.25 && 1291 < ncols && ncols <= 2613 )
3964  separate = TRUE;
3965  else if ( nrows > 8146 && 0.75 < firstlptime && firstlptime <= 6.25 && 655 < ncols && ncols <= 1290 )
3966  separate = TRUE;
3967  else if ( nrows > 8146 && 0.75 < firstlptime && firstlptime <= 6.25 && 1291 < ncols && ncols <= 2613 )
3968  separate = TRUE;
3969  else if ( nrows > 8146 && firstlptime > 6.25 )
3970  separate = TRUE;
3971 
3972  if ( ! separate )
3973  {
3974  return SCIP_OKAY;
3975  }
3976  }
3977 
3978  /* preceed with separation */
3979  *result = SCIP_DIDNOTFIND;
3980 
3981  SCIPdebugMessage("separating CG-cuts via sub-MIPs: %d cols, %d rows\n", ncols, nrows);
3982 
3983  /* prepare data */
3984  SCIP_CALL( SCIPallocBlockMemory(scip, &mipdata) );
3985  mipdata->subscip = NULL;
3986  mipdata->alpha = NULL;
3987  mipdata->fracalpha = NULL;
3988  mipdata->beta = NULL;
3989  mipdata->fracbeta = NULL;
3990  mipdata->coltype = NULL;
3991  mipdata->iscomplemented = NULL;
3992  mipdata->isshifted = NULL;
3993  mipdata->ylhs = NULL;
3994  mipdata->yrhs = NULL;
3995  mipdata->z = NULL;
3996  mipdata->normtype = ' ';
3997 
3998  mipdata->conshdlrfullnorm = CONSHDLRFULLNORM;
3999  mipdata->scip = scip;
4000  mipdata->sepa = sepa;
4001  mipdata->sepadata = sepadata;
4002 
4003 
4004  /* get the type of norm to use for efficacy calculations */
4005  SCIP_CALL( SCIPgetCharParam(scip, "separating/efficacynorm", &mipdata->normtype) );
4006 
4007  /* create subscip */
4008  SCIP_CALL( createSubscip(scip, sepa, sepadata, mipdata) );
4009 
4010  /* set parameters */
4011  SCIP_CALL( subscipSetParams(sepadata, mipdata, &success) );
4012 
4013  if ( success && !SCIPisStopped(scip) )
4014  {
4015  /* solve subscip */
4016  SCIP_CALL( solveSubscip(scip, sepadata, mipdata, &success) );
4017 
4018  /* preceed if solution was successful */
4019  if ( success && ! SCIPisStopped(scip) )
4020  {
4021  SCIP_CALL( createCGCuts(scip, sepa, sepadata, mipdata, &cutoff, &ngen) );
4022  }
4023  }
4024 
4025  SCIP_CALL( freeSubscip(scip, sepa, mipdata) );
4026  SCIPfreeBlockMemory(scip, &mipdata);
4027 
4028  SCIPdebugMessage("Found %u CG-cuts.\n", ngen);
4029 
4030  if ( cutoff )
4031  *result = SCIP_CUTOFF;
4032  else if ( ngen > 0 )
4033  *result = SCIP_SEPARATED;
4034 
4035 #ifdef SCIP_OUTPUT
4036  /* SCIP_CALL( SCIPwriteLP(scip, "cuts.lp") ); */
4037  /* SCIP_CALL( SCIPwriteMIP(scip, "cuts.lp", FALSE, TRUE) ); */
4038 #endif
4039 
4040  return SCIP_OKAY;
4041 }
4042 
4043 
4044 /*
4045  * separator specific interface methods
4046  */
4047 
4048 /** creates the CGMIP MIR cut separator and includes it in SCIP */
4050  SCIP* scip /**< SCIP data structure */
4051  )
4052 {
4053  SCIP_SEPADATA* sepadata;
4054  SCIP_SEPA* sepa;
4055 
4056  /* create separator data */
4057  SCIP_CALL( SCIPallocMemory(scip, &sepadata) );
4058 
4059  sepa = NULL;
4060  /* include separator */
4062  SEPA_USESSUBSCIP, SEPA_DELAY, sepaExeclpCGMIP, NULL, sepadata) );
4063  assert(sepa != NULL);
4064 
4065  SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyCGMIP) );
4066  SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeCGMIP) );
4067 
4068  /* add separator parameters */
4069  SCIP_CALL( SCIPaddIntParam(scip,
4070  "separating/"SEPA_NAME"/maxrounds",
4071  "maximal number of cgmip separation rounds per node (-1: unlimited)",
4072  &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
4073 
4074  SCIP_CALL( SCIPaddIntParam(scip,
4075  "separating/"SEPA_NAME"/maxroundsroot",
4076  "maximal number of cgmip separation rounds in the root node (-1: unlimited)",
4077  &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
4078 
4079  SCIP_CALL( SCIPaddIntParam(scip,
4080  "separating/"SEPA_NAME"/maxdepth",
4081  "maximal depth at which the separator is applied (-1: unlimited)",
4082  &sepadata->maxdepth, FALSE, DEFAULT_MAXDEPTH, -1, INT_MAX, NULL, NULL) );
4083 
4085  "separating/"SEPA_NAME"/decisiontree",
4086  "Use decision tree to turn separation on/off?",
4087  &sepadata->decisiontree, FALSE, DEFAULT_DECISIONTREE, NULL, NULL) );
4088 
4090  "separating/"SEPA_NAME"/timelimit",
4091  "time limit for sub-MIP",
4092  &sepadata->timelimit, TRUE, DEFAULT_TIMELIMIT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
4093 
4095  "separating/"SEPA_NAME"/memorylimit",
4096  "memory limit for sub-MIP",
4097  &sepadata->memorylimit, TRUE, DEFAULT_MEMORYLIMIT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
4098 
4100  "separating/"SEPA_NAME"/minnodelimit",
4101  "minimum number of nodes considered for sub-MIP (-1: unlimited)",
4102  &sepadata->minnodelimit, FALSE, DEFAULT_MINNODELIMIT, -1LL, SCIP_LONGINT_MAX, NULL, NULL) );
4103 
4105  "separating/"SEPA_NAME"/maxnodelimit",
4106  "maximum number of nodes considered for sub-MIP (-1: unlimited)",
4107  &sepadata->maxnodelimit, FALSE, DEFAULT_MAXNODELIMIT, -1LL, SCIP_LONGINT_MAX, NULL, NULL) );
4108 
4110  "separating/"SEPA_NAME"/cutcoefbnd",
4111  "bounds on the values of the coefficients in the CG-cut",
4112  &sepadata->cutcoefbnd, TRUE, DEFAULT_CUTCOEFBND, 0.0, SCIP_REAL_MAX, NULL, NULL) );
4113 
4115  "separating/"SEPA_NAME"/onlyactiverows",
4116  "Use only active rows to generate cuts?",
4117  &sepadata->onlyactiverows, FALSE, DEFAULT_ONLYACTIVEROWS, NULL, NULL) );
4118 
4119  SCIP_CALL( SCIPaddIntParam(scip,
4120  "separating/"SEPA_NAME"/maxrowage",
4121  "maximal age of rows to consider if onlyactiverows is false",
4122  &sepadata->maxrowage, FALSE, DEFAULT_MAXROWAGE, -1, INT_MAX, NULL, NULL) );
4123 
4125  "separating/"SEPA_NAME"/onlyrankone",
4126  "Separate only rank 1 inequalities w.r.t. CG-MIP separator?",
4127  &sepadata->onlyrankone, FALSE, DEFAULT_ONLYRANKONE, NULL, NULL) );
4128 
4130  "separating/"SEPA_NAME"/onlyintvars",
4131  "Generate cuts for problems with only integer variables?",
4132  &sepadata->onlyintvars, FALSE, DEFAULT_ONLYINTVARS, NULL, NULL) );
4133 
4135  "separating/"SEPA_NAME"/allowlocal",
4136  "Allow to generate local cuts?",
4137  &sepadata->allowlocal, FALSE, DEFAULT_ALLOWLOCAL, NULL, NULL) );
4138 
4140  "separating/"SEPA_NAME"/contconvert",
4141  "Convert some integral variables to be continuous to reduce the size of the sub-MIP?",
4142  &sepadata->contconvert, FALSE, DEFAULT_CONTCONVERT, NULL, NULL) );
4143 
4145  "separating/"SEPA_NAME"/contconvfrac",
4146  "fraction of integral variables converted to be continuous (if contconvert)",
4147  &sepadata->contconvfrac, FALSE, DEFAULT_CONTCONVFRAC, 0.0, 1.0, NULL, NULL) );
4148 
4149  SCIP_CALL( SCIPaddIntParam(scip,
4150  "separating/"SEPA_NAME"/contconvmin",
4151  "minimum number of integral variables before some are converted to be continuous",
4152  &sepadata->contconvmin, FALSE, DEFAULT_CONTCONVMIN, -1, INT_MAX, NULL, NULL) );
4153 
4155  "separating/"SEPA_NAME"/intconvert",
4156  "Convert some integral variables attaining fractional values to have integral value?",
4157  &sepadata->intconvert, FALSE, DEFAULT_INTCONVERT, NULL, NULL) );
4158 
4160  "separating/"SEPA_NAME"/intconvfrac",
4161  "fraction of frac. integral variables converted to have integral value (if intconvert)",
4162  &sepadata->intconvfrac, FALSE, DEFAULT_INTCONVFRAC, 0.0, 1.0, NULL, NULL) );
4163 
4164  SCIP_CALL( SCIPaddIntParam(scip,
4165  "separating/"SEPA_NAME"/intconvmin",
4166  "minimum number of integral variables before some are converted to have integral value",
4167  &sepadata->intconvmin, FALSE, DEFAULT_INTCONVMIN, -1, INT_MAX, NULL, NULL) );
4168 
4170  "separating/"SEPA_NAME"/skipmultbounds",
4171  "Skip the upper bounds on the multipliers in the sub-MIP?",
4172  &sepadata->skipmultbounds, FALSE, DEFAULT_SKIPMULTBOUNDS, NULL, NULL) );
4173 
4175  "separating/"SEPA_NAME"/objlone",
4176  "Should the objective of the sub-MIP minimize the l1-norm of the multipliers?",
4177  &sepadata->objlone, FALSE, DEFAULT_OBJLONE, NULL, NULL) );
4178 
4180  "separating/"SEPA_NAME"/objweight",
4181  "weight used for the row combination coefficient in the sub-MIP objective",
4182  &sepadata->objweight, TRUE, DEFAULT_OBJWEIGHT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
4183 
4185  "separating/"SEPA_NAME"/objweightsize",
4186  "Weight each row by its size?",
4187  &sepadata->objweightsize, FALSE, DEFAULT_OBJWEIGHTSIZE, NULL, NULL) );
4188 
4190  "separating/"SEPA_NAME"/dynamiccuts",
4191  "should generated cuts be removed from the LP if they are no longer tight?",
4192  &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
4193 
4195  "separating/"SEPA_NAME"/usecmir",
4196  "use CMIR-generator (otherwise add cut directly)?",
4197  &sepadata->usecmir, FALSE, DEFAULT_USECMIR, NULL, NULL) );
4198 
4200  "separating/"SEPA_NAME"/usestrongcg",
4201  "use strong CG-function to strengthen cut?",
4202  &sepadata->usestrongcg, FALSE, DEFAULT_USESTRONGCG, NULL, NULL) );
4203 
4205  "separating/"SEPA_NAME"/cmirownbounds",
4206  "tell CMIR-generator which bounds to used in rounding?",
4207  &sepadata->cmirownbounds, FALSE, DEFAULT_CMIROWNBOUNDS, NULL, NULL) );
4208 
4210  "separating/"SEPA_NAME"/usecutpool",
4211  "use cutpool to store CG-cuts even if the are not efficient?",
4212  &sepadata->usecutpool, FALSE, DEFAULT_USECUTPOOL, NULL, NULL) );
4213 
4215  "separating/"SEPA_NAME"/primalseparation",
4216  "only separate cuts that are tight for the best feasible solution?",
4217  &sepadata->primalseparation, FALSE, DEFAULT_PRIMALSEPARATION, NULL, NULL) );
4218 
4220  "separating/"SEPA_NAME"/earlyterm",
4221  "terminate separation if a violated (but possibly sub-optimal) cut has been found?",
4222  &sepadata->earlyterm, FALSE, DEFAULT_EARLYTERM, NULL, NULL) );
4223 
4225  "separating/"SEPA_NAME"/addviolationcons",
4226  "add constraint to subscip that only allows violated cuts (otherwise add obj. limit)?",
4227  &sepadata->addviolationcons, FALSE, DEFAULT_ADDVIOLATIONCONS, NULL, NULL) );
4228 
4230  "separating/"SEPA_NAME"/addviolconshdlr",
4231  "add constraint handler to filter out violated cuts?",
4232  &sepadata->addviolconshdlr, FALSE, DEFAULT_ADDVIOLCONSHDLR, NULL, NULL) );
4233 
4235  "separating/"SEPA_NAME"/conshdlrusenorm",
4236  "should the violation constraint handler use the norm of a cut to check for feasibility?",
4237  &sepadata->conshdlrusenorm, FALSE, DEFAULT_CONSHDLRUSENORM, NULL, NULL) );
4238 
4240  "separating/"SEPA_NAME"/useobjub",
4241  "Use upper bound on objective function (via primal solution)?",
4242  &sepadata->useobjub, FALSE, DEFAULT_USEOBJUB, NULL, NULL) );
4243 
4245  "separating/"SEPA_NAME"/useobjlb",
4246  "Use lower bound on objective function (via primal solution)?",
4247  &sepadata->useobjlb, FALSE, DEFAULT_USEOBJLB, NULL, NULL) );
4248 
4250  "separating/"SEPA_NAME"/subscipfast",
4251  "Should the settings for the sub-MIP be optimized for speed?",
4252  &sepadata->subscipfast, FALSE, DEFAULT_SUBSCIPFAST, NULL, NULL) );
4253 
4255  "separating/"SEPA_NAME"/output",
4256  "Should information about the sub-MIP and cuts be displayed?",
4257  &sepadata->output, FALSE, DEFAULT_OUTPUT, NULL, NULL) );
4258 
4259  return SCIP_OKAY;
4260 }
4261