Scippy

SCIP

Solving Constraint Integer Programs

sepa_flowcover.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 
16 /**@file sepa_flowcover.c
17  * @brief flow cover cuts separator
18  * @author Kati Wolter
19  * @author Tobias Achterberg
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #include <assert.h>
25 #include <string.h>
26 
27 #include "scip/sepa_flowcover.h"
28 #include "scip/cons_knapsack.h"
29 #include "scip/pub_misc.h"
30 
31 
32 #define SEPA_NAME "flowcover"
33 #define SEPA_DESC "flow cover cuts separator (c-MIR approach)"
34 #define SEPA_PRIORITY -4000
35 #define SEPA_FREQ 0
36 #define SEPA_MAXBOUNDDIST 0.0
37 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
38 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
39 
40 #define DEFAULT_MAXROUNDS 5 /**< maximal number of separation rounds per node (-1: unlimited) */
41 #define DEFAULT_MAXROUNDSROOT 15 /**< maximal number of separation rounds in the root node (-1: unlimited) */
42 #define DEFAULT_MAXTRIES 100 /**< maximal number of rows to separate flow cover cuts for per separation round
43  * (-1: unlimited) */
44 #define DEFAULT_MAXTRIESROOT -1 /**< maximal number of rows to separate flow cover cuts for per separation round
45  * in the root (-1: unlimited) */
46 #define DEFAULT_MAXFAILS 50 /**< maximal number of consecutive fails to generate a cut per separation round
47  * (-1: unlimited) */
48 #define DEFAULT_MAXFAILSROOT 100 /**< maximal number of consecutive fails to generate a cut per separation round
49  * in the root (-1: unlimited) */
50 #define DEFAULT_MAXSEPACUTS 100 /**< maximal number of flow cover cuts separated per separation round */
51 #define DEFAULT_MAXSEPACUTSROOT 200 /**< maximal number of flow cover cuts separated per separation round in the root */
52 #define DEFAULT_MAXSLACK SCIP_REAL_MAX /**< maximal slack of rows to separate flow cover cuts for */
53 #define DEFAULT_MAXSLACKROOT SCIP_REAL_MAX /**< maximal slack of rows to separate flow cover cuts for in the root */
54 #define DEFAULT_SLACKSCORE 1e-03 /**< weight of slack in the scoring of the rows */
55 #define DEFAULT_MAXROWDENSITY 1.0 /**< maximal density of rows to separate flow cover cuts for */
56 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
57 #define DEFAULT_MAXTESTDELTA 10 /**< cut generation heuristic: maximal number of different deltas to try */
58 #define DEFAULT_MULTBYMINUSONE TRUE /**< should flow cover cuts be separated for 0-1 single node flow set with reversed arcs in addition? */
59 
60 #define BOUNDSWITCH 0.5
61 #define ALLOWLOCAL TRUE
62 #define DENSSCORE 1e-04
63 /*#define MAKECONTINTEGRAL FALSE*/
64 #define MINFRAC 0.01
65 #define MAXFRAC 0.95
66 #define FIXINTEGRALRHS FALSE
67 #define MAXDNOM 1000LL
68 #define MINDELTA 1e-03
69 #define MAXDELTA 1e-09
70 #define MAXSCALE 1000.0
71 #define MAXDYNPROGSPACE 1000000
72 
73 #define MAXAGGRLEN(nvars) (0.1*(nvars)+1000) /**< maximal length of base inequality */
74 #define MAXABSVBCOEF 1e+5 /**< maximal absolute coefficient in variable bounds used for snf relaxation */
75 #define MAXBOUND 1e+10 /**< maximal value of normal bounds used for snf relaxation */
76 
77 
78 /*
79  * Data structures
80  */
81 
82 /** separator data */
83 struct SCIP_SepaData
84 {
85  int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */
86  int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */
87  int maxtries; /**< maximal number of rows to separate flow cover cuts for per separation round
88  * (-1: unlimited) */
89  int maxtriesroot; /**< maximal number of rows to separate flow cover cuts for per separation round
90  * in the root (-1: unlimited) */
91  int maxfails; /**< maximal number of consecutive fails to generate a cut per separation round
92  * (-1: unlimited) */
93  int maxfailsroot; /**< maximal number of consecutive fails to generate a cut per separation round
94  * in the root (-1: unlimited) */
95  int maxsepacuts; /**< maximal number of flow cover cuts separated per separation round */
96  int maxsepacutsroot; /**< maximal number of flow cover cuts separated per separation round in the root */
97  SCIP_Real maxslack; /**< maximal slack of rows to separate flow cover cuts for */
98  SCIP_Real maxslackroot; /**< maximal slack of rows to separate flow cover cuts for in the root */
99  SCIP_Real slackscore; /**< weight of slack in the scoring of the rows */
100  SCIP_Real maxrowdensity; /**< maximal density of rows to separate flow cover cuts for */
101  SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */
102  SCIP_Bool multbyminusone; /**< should flow cover cuts be separated for 0-1 single node flow set with reversed arcs in addition? */
103  int maxtestdelta; /**< cut generation heuristic: maximal number of different deltas to try */
104 };
105 
106 
107 /*
108  * Local methods
109  */
110 
111 /** get LP solution value and index of variable lower bound (with binary variable) which is closest to the current LP
112  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
113  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
114  * given variable
115  */
116 static
117 SCIP_RETCODE getClosestVlb(
118  SCIP* scip, /**< SCIP data structure */
119  SCIP_VAR* var, /**< given active problem variable */
120  SCIP_Real bestsub, /**< closest simple upper bound of given variable */
121  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
122  SCIP_Real* rowcoefsbinary, /**< coefficient of all binary problem variables in current row */
123  SCIP_Real* varsolvals, /**< LP solution value of all active problem variables */
124  int* assoctransvars, /**< associated var in relaxed set for all vars of row; construction is not finished yet */
125  SCIP_Real* closestvlb, /**< pointer to store the LP sol value of the closest variable lower bound */
126  int* closestvlbidx /**< pointer to store the index of the closest vlb; -1 if no vlb was found */
127  )
128 {
129  int nvlbs;
130 
131  assert(scip != NULL);
132  assert(var != NULL);
133  assert(bestsub == SCIPvarGetUbGlobal(var) || bestsub == SCIPvarGetUbLocal(var)); /*lint !e777*/
134  assert(!SCIPisInfinity(scip, bestsub));
135  assert(!SCIPisZero(scip, rowcoef));
136  assert(rowcoefsbinary != NULL);
137  assert(varsolvals != NULL);
138  assert(assoctransvars != NULL);
139  assert(closestvlb != NULL);
140  assert(closestvlbidx != NULL);
141 
142  nvlbs = SCIPvarGetNVlbs(var);
143 
144  *closestvlbidx = -1;
145  *closestvlb = -SCIPinfinity(scip);
146  if( nvlbs > 0 )
147  {
148  SCIP_VAR** vlbvars;
149  SCIP_Real* vlbcoefs;
150  SCIP_Real* vlbconsts;
151  int i;
152 
153  vlbvars = SCIPvarGetVlbVars(var);
154  vlbcoefs = SCIPvarGetVlbCoefs(var);
155  vlbconsts = SCIPvarGetVlbConstants(var);
156 
157  for( i = 0; i < nvlbs; i++ )
158  {
159  SCIP_Real rowcoefbinvar;
160  SCIP_Real val1;
161  SCIP_Real val2;
162  SCIP_Real vlbsol;
163  SCIP_Bool meetscriteria;
164  int probidxbinvar;
165 
166  /* use only variable lower bounds l~_i * x_i + d_i with x_i binary which are active */
167  if( !SCIPvarIsBinary(vlbvars[i]) || !SCIPvarIsActive(vlbvars[i]) )
168  continue;
169 
170  /* check if current variable lower bound l~_i * x_i + d_i imposed on y_j meets the following criteria:
171  * (let a_j = coefficient of y_j in current row,
172  * u_j = closest simple upper bound imposed on y_j,
173  * c_i = coefficient of x_i in current row)
174  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k yet
175  * if a_j > 0:
176  * 1. u_j <= d_i
177  * 2. a_j ( u_j - d_i ) + c_i <= 0
178  * 3. a_j l~_i + c_i <= 0
179  * if a_j < 0:
180  * 1. u_j <= d_i
181  * 2. a_j ( u_j - d_i ) + c_i >= 0
182  * 3. a_j l~_i + c_i >= 0
183  */
184  probidxbinvar = SCIPvarGetProbindex(vlbvars[i]);
185  rowcoefbinvar = rowcoefsbinary[probidxbinvar];
186 
187  val1 = ( rowcoef * ( bestsub - vlbconsts[i] ) ) + rowcoefbinvar;
188  val2 = ( rowcoef * vlbcoefs[i] ) + rowcoefbinvar;
189 
190  meetscriteria = FALSE;
191  if( SCIPisPositive(scip, rowcoef) )
192  {
193  if( assoctransvars[probidxbinvar] == -1 && SCIPisFeasLE(scip, bestsub, vlbconsts[i])
194  && SCIPisFeasLE(scip, val1, 0.0) && SCIPisFeasLE(scip, val2, 0.0) )
195  meetscriteria = TRUE;
196  }
197  else
198  {
199  assert(SCIPisNegative(scip, rowcoef));
200  if( assoctransvars[probidxbinvar] == -1 && SCIPisFeasLE(scip, bestsub, vlbconsts[i])
201  && SCIPisFeasGE(scip, val1, 0.0) && SCIPisFeasGE(scip, val2, 0.0) )
202  meetscriteria = TRUE;
203  }
204 
205  /* variable lower bound does not meet criteria */
206  if( !meetscriteria )
207  continue;
208 
209  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
210  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
211  */
212  if( REALABS(vlbcoefs[i]) > MAXABSVBCOEF || SCIPisInfinity(scip, REALABS(val2)) )
213  continue;
214 
215  vlbsol = (vlbcoefs[i] * varsolvals[probidxbinvar]) + vlbconsts[i];
216  if( SCIPisGT(scip, vlbsol, *closestvlb) )
217  {
218  *closestvlb = vlbsol;
219  *closestvlbidx = i;
220  }
221  assert(*closestvlbidx >= 0);
222 
223  }
224  }
225 
226  return SCIP_OKAY;
227 }
228 
229 /** get LP solution value and index of variable upper bound (with binary variable) which is closest to the current LP
230  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
231  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
232  * given variable
233  */
234 static
235 SCIP_RETCODE getClosestVub(
236  SCIP* scip, /**< SCIP data structure */
237  SCIP_VAR* var, /**< given active problem variable */
238  SCIP_Real bestslb, /**< closest simple lower bound of given variable */
239  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
240  SCIP_Real* rowcoefsbinary, /**< coefficient of all binary problem variables in current row */
241  SCIP_Real* varsolvals, /**< LP solution value of all active problem variables */
242  int* assoctransvars, /**< associated var in relaxed set for all vars of row; construction is not finished yet */
243  SCIP_Real* closestvub, /**< pointer to store the LP sol value of the closest variable upper bound */
244  int* closestvubidx /**< pointer to store the index of the closest vub; -1 if no vub was found */
245  )
246 {
247  int nvubs;
248 
249  assert(scip != NULL);
250  assert(var != NULL);
251  assert(bestslb == SCIPvarGetLbGlobal(var) || bestslb == SCIPvarGetLbLocal(var)); /*lint !e777*/
252  assert(!SCIPisInfinity(scip, - bestslb));
253  assert(!SCIPisZero(scip, rowcoef));
254  assert(rowcoefsbinary != NULL);
255  assert(varsolvals != NULL);
256  assert(assoctransvars != NULL);
257  assert(closestvub != NULL);
258  assert(closestvubidx != NULL);
259 
260  nvubs = SCIPvarGetNVubs(var);
261 
262  *closestvubidx = -1;
263  *closestvub = SCIPinfinity(scip);
264  if( nvubs > 0 )
265  {
266  SCIP_VAR** vubvars;
267  SCIP_Real* vubcoefs;
268  SCIP_Real* vubconsts;
269  int i;
270 
271  vubvars = SCIPvarGetVubVars(var);
272  vubcoefs = SCIPvarGetVubCoefs(var);
273  vubconsts = SCIPvarGetVubConstants(var);
274 
275  for( i = 0; i < nvubs; i++ )
276  {
277  SCIP_Real rowcoefbinvar;
278  SCIP_Real val1;
279  SCIP_Real val2;
280  SCIP_Real vubsol;
281  SCIP_Bool meetscriteria;
282  int probidxbinvar;
283 
284  /* use only variable upper bound u~_i * x_i + d_i with x_i binary and which are active */
285  if( !SCIPvarIsBinary(vubvars[i]) || !SCIPvarIsActive(vubvars[i]))
286  continue;
287 
288  /* checks if current variable upper bound u~_i * x_i + d_i meets the following criteria
289  * (let a_j = coefficient of y_j in current row,
290  * l_j = closest simple lower bound imposed on y_j,
291  * c_i = coefficient of x_i in current row)
292  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k
293  * if a > 0:
294  * 1. l_j >= d_i
295  * 2. a_j ( l_i - d_i ) + c_i >= 0
296  * 3. a_j u~_i + c_i >= 0
297  * if a < 0:
298  * 1. l_j >= d_i
299  * 2. a_j ( l_j - d_i ) + c_i <= 0
300  * 3. a_j u~_i + c_i <= 0
301  */
302  probidxbinvar = SCIPvarGetProbindex(vubvars[i]);
303  rowcoefbinvar = rowcoefsbinary[probidxbinvar];
304 
305  val1 = ( rowcoef * ( bestslb - vubconsts[i] ) ) + rowcoefbinvar;
306  val2 = ( rowcoef * vubcoefs[i] ) + rowcoefbinvar;
307 
308  meetscriteria = FALSE;
309  if( SCIPisPositive(scip, rowcoef) )
310  {
311  if( assoctransvars[probidxbinvar] == -1 && SCIPisFeasGE(scip, bestslb, vubconsts[i])
312  && SCIPisFeasGE(scip, val1, 0.0) && SCIPisFeasGE(scip, val2, 0.0) )
313  meetscriteria = TRUE;
314  }
315  else
316  {
317  assert(SCIPisNegative(scip, rowcoef));
318  if( assoctransvars[probidxbinvar] == -1 && SCIPisFeasGE(scip, bestslb, vubconsts[i])
319  && SCIPisFeasLE(scip, val1, 0.0) && SCIPisFeasLE(scip, val2, 0.0) )
320  meetscriteria = TRUE;
321  }
322 
323  /* variable upper bound does not meet criteria */
324  if( !meetscriteria )
325  continue;
326 
327  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
328  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
329  */
330  if( REALABS(vubcoefs[i]) > MAXABSVBCOEF || SCIPisInfinity(scip, REALABS(val2)) )
331  continue;
332 
333  vubsol = vubcoefs[i] * varsolvals[probidxbinvar] + vubconsts[i];
334  if( SCIPisLT(scip, vubsol, *closestvub) )
335  {
336  *closestvub = vubsol;
337  *closestvubidx = i;
338  }
339  assert(*closestvubidx >= 0);
340  }
341  }
342 
343  return SCIP_OKAY;
344 }
345 
346 /** return global or local lower bound of given variable whichever is closer to the variables current LP solution value */
347 static
348 void getClosestLb(
349  SCIP* scip, /**< SCIP data structure */
350  SCIP_VAR* var, /**< active problem variable */
351  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
352  SCIP_Real* closestlb, /**< pointer to store the value of the closest variable lower bound */
353  int* closestlbtype /**< pointer to store type of closest bound; -1 if global lb, -2 otherwise */
354  )
355 {
356  assert(closestlb != NULL);
357  assert(closestlbtype != NULL);
358 
359  *closestlb = SCIPvarGetLbGlobal(var);
360  *closestlbtype = -1;
361  if( allowlocal )
362  {
363  SCIP_Real loclb;
364  loclb = SCIPvarGetLbLocal(var);
365  if( SCIPisGT(scip, loclb, *closestlb) )
366  {
367  *closestlb = loclb;
368  *closestlbtype = -2;
369  }
370  }
371 
372  /* due to numerical reasons, huge bounds are relaxed to infinite bounds; this way the bounds are not used for
373  * the construction of the 0-1 single node flow relaxation
374  */
375  if( *closestlb <= -MAXBOUND )
376  *closestlb = -SCIPinfinity(scip);
377 }
378 
379 /** return global or local upper bound of given variable whichever is closer to the variables current LP solution value */
380 static
381 void getClosestUb(
382  SCIP* scip, /**< SCIP data structure */
383  SCIP_VAR* var, /**< active problem variable */
384  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
385  SCIP_Real* closestub, /**< pointer to store the value of the closest upper bound */
386  int* closestubtype /**< pointer to store type of closest bound; -1 if global ub, -2 otherwise */
387  )
388 {
389  assert(closestub != NULL);
390  assert(closestubtype != NULL);
391 
392  *closestub = SCIPvarGetUbGlobal(var);
393  *closestubtype = -1;
394  if( allowlocal )
395  {
396  SCIP_Real locub;
397  locub = SCIPvarGetUbLocal(var);
398  if( SCIPisLT(scip, locub, *closestub) )
399  {
400  *closestub = locub;
401  *closestubtype = -2;
402  }
403  }
404 
405  /* due to numerical reasons, huge bounds are relaxed to infinite bounds; this way the bounds are not used for
406  * the construction of the 0-1 single node flow relaxation
407  */
408  if( *closestub >= MAXBOUND )
409  *closestub = SCIPinfinity(scip);
410 }
411 
412 /** construct a 0-1 single node flow relaxation (with some additional simple constraints) of a mixed integer set
413  * corresponding to the given row lhs <= a * x + const <= rhs; depending on the given values rowweight and scale
414  * the mixed integer set which should be used is defined by the mixed integer constraint
415  * a * (x,y) <= rhs - const if (rowweight = 1, scale = 1) or (rowweight = -1, scale = -1, rhs < infinity)
416  * - a * (x,y) <= - (lhs - const) if (rowweight = -1, scale = 1) or (rowweight = 1, scale = -1, lhs > -infinity)
417  * - a * (x,y) - s <= - (rhs - const) if (rowweight = 1, scale = -1, lhs = -infinity)
418  * a * (x,y) - s <= (lhs - const) if (rowweight = -1, scale = -1, rhs = infinity)
419  */
420 static
421 SCIP_RETCODE constructSNFRelaxation(
422  SCIP* scip, /**< SCIP data structure */
423  SCIP_VAR** vars, /**< active problem variables */
424  int nvars, /**< number of active problem variables */
425  SCIP_Real* varsolvals, /**< solution values of active problem variables */
426  SCIP_ROW* row, /**< given row */
427  SCIP_Real rowweight, /**< weight of given row; can be +1 or -1 */
428  SCIP_Real scale, /**< additional scaling factor for given row */
429  int* boundsfortrans, /**< pointer to store bound used for all non-binary vars of row */
430  SCIP_BOUNDTYPE* boundtypesfortrans, /**< pointer to store type of bound used for all non-binary vars of row */
431  int* assoctransvars, /**< pointer to store associated var in relaxed set for all vars of row */
432  int* transvarcoefs, /**< pointer to store coefficient of all vars in relaxed set */
433  SCIP_Real* transbinvarsolvals, /**< pointer to store sol val of bin var in vub of all vars in relaxed set */
434  SCIP_Real* transcontvarsolvals,/**< pointer to store sol val of all real vars in relaxed set */
435  SCIP_Real* transvarvubcoefs, /**< pointer to store coefficient in vub of all vars in relaxed set */
436  int* ntransvars, /**< pointer to store number of vars in relaxed set */
437  SCIP_Real* transrhs, /**< pointer to store rhs in relaxed set */
438  SCIP_Bool* success /**< pointer to store whether a relaxation was constructed */
439  )
440 {
441  SCIP_COL** nonzcols;
442  SCIP_Real* nonzcoefs;
443  SCIP_Real* rowcoefsbinary;
444  int* nonzcolsbinary;
445  int* nonzcolsnonbinary;
446  int nnonzcols;
447  int nnonzcolsbinary;
448  int nnonzcolsnonbinary;
449  int c;
450 
451  assert(scip != NULL);
452  assert(vars!= NULL);
453  assert(nvars > 0);
454  assert(varsolvals != NULL);
455  assert(row != NULL);
456  assert(rowweight == 1.0 || rowweight == -1.0);
457  assert(( rowweight == 1.0 && !SCIPisInfinity(scip, SCIProwGetRhs(row)) )
458  || ( rowweight == -1.0 && !SCIPisInfinity(scip, -SCIProwGetLhs(row)) ));
459  assert(scale == 1.0 || scale == -1.0);
460  assert(boundsfortrans != NULL);
461  assert(boundtypesfortrans != NULL);
462  assert(assoctransvars != NULL);
463  assert(transvarcoefs != NULL);
464  assert(transbinvarsolvals != NULL);
465  assert(transcontvarsolvals != NULL);
466  assert(transvarvubcoefs != NULL);
467  assert(ntransvars != NULL);
468  assert(transrhs != NULL);
469  assert(success != NULL);
470 
471  *success = FALSE;
472 
473  SCIPdebugMessage("--------------------- construction of SNF relaxation ------------------------------------\n");
474 
475  /* get nonzero columns and coefficients of row */
476  nonzcols = SCIProwGetCols(row);
477  nnonzcols = SCIProwGetNLPNonz(row);
478  nonzcoefs = SCIProwGetVals(row);
479  SCIPdebugMessage("nnonzcols = %d\n",nnonzcols);
480 
481  /* get data structures */
482  SCIP_CALL( SCIPallocBufferArray(scip, &nonzcolsbinary, nnonzcols) );
483  SCIP_CALL( SCIPallocBufferArray(scip, &nonzcolsnonbinary, nnonzcols) );
484  SCIP_CALL( SCIPallocBufferArray(scip, &rowcoefsbinary, nvars) );
485 
486  /* store nonzero columns representing binary and non-binary variables, and get active binary problem variables the
487  * coefficient in the row
488  */
489  nnonzcolsbinary = 0;
490  nnonzcolsnonbinary = 0;
491  BMSclearMemoryArray(rowcoefsbinary, nvars);
492  for( c = 0; c < nnonzcols; c++ )
493  {
494  SCIP_COL* col;
495  SCIP_VAR* var;
496 
497  col = nonzcols[c];
498  var = SCIPcolGetVar(col);
499 
500  assert(!SCIPisZero(scip, nonzcoefs[c]));
501 
502  if( SCIPvarIsBinary(var) )
503  {
504  /* saves column for binary variable */
505  nonzcolsbinary[nnonzcolsbinary] = c;
506  nnonzcolsbinary++;
507 
508  /* saves row coefficient of binary variable */
509  assert(SCIPvarGetProbindex(var) > -1 && SCIPvarGetProbindex(var) < nvars);
510  rowcoefsbinary[SCIPvarGetProbindex(var)] = rowweight * scale * nonzcoefs[c];
511  }
512  else
513  {
514  /* saves column for non-binary variable */
515  nonzcolsnonbinary[nnonzcolsnonbinary] = c;
516  nnonzcolsnonbinary++;
517  }
518  }
519  assert(nnonzcolsbinary + nnonzcolsnonbinary == nnonzcols);
520 
521  /* initialize data structures */
522  for( c = 0; c < nvars; c++ )
523  {
524  assoctransvars[c] = -1;
525  boundsfortrans[c] = -3;
526  }
527  *ntransvars = 0;
528 
529  /* initialize right hand side of constraint in 0-1 single node flow relaxation */
530  if( rowweight * scale == 1.0 && !SCIPisInfinity(scip, SCIProwGetRhs(row)) )
531  *transrhs = SCIProwGetRhs(row) - SCIProwGetConstant(row);
532  else if( rowweight * scale == 1.0 && SCIPisInfinity(scip, SCIProwGetRhs(row)) )
533  {
534  assert(rowweight == -1.0 && scale == -1.0);
535  *transrhs = SCIProwGetLhs(row) - SCIProwGetConstant(row);
536  }
537  else if( rowweight * scale == -1.0 && !SCIPisInfinity(scip, -SCIProwGetLhs(row)) )
538  *transrhs = - SCIProwGetLhs(row) + SCIProwGetConstant(row);
539  else
540  {
541  assert(rowweight == 1.0 && scale == -1.0 && SCIPisInfinity(scip, -SCIProwGetLhs(row)));
542  *transrhs = - SCIProwGetRhs(row) + SCIProwGetConstant(row);
543  }
544 
545  /* for each non-binary variable y_j in the row with nonzero row coefficient perform
546  * 1. get closest simple or variable lower bound and closest simple or variable upper bound
547  * 2. decide which bound is used to define the real variable y'_j in the 0-1 single node flow relaxation
548  * 3. construct y'_j with 0 <= y'_j <= u'_j x_j
549  * 4. store for y_j and x_j (if x_j is a binary variable in the row) that y'_j is the associated real variable
550  * in the 0-1 single node flow relaxation and for y_j the bound used to define y'_j.
551  *
552  * for each binary variable x_j in the row which has not been handled with a non-binary variable perform
553  * 1. construct y'_j with 0 <= y'_j <= u'_j x_j
554  * 2. store for x_j that y'_j is the associated real variable in the 0-1 single node flow relaxation.
555  *
556  * start with non-binary variables because a binary variable x_j which is involved in a used variable bound
557  * imposed on a non-binary variable y_j has to be handled together with the non-binary variable y_j.
558  */
559  SCIPdebugMessage("transformation for NONBINARY variables (nnonbinvars=%d):\n", nnonzcolsnonbinary);
560 
561  /* non-binary variables and binary variables contained in used variable bounds */
562  for( c = 0; c < nnonzcolsnonbinary; c++ )
563  {
564  SCIP_VAR* var;
565  SCIP_Real bestlb;
566  SCIP_Real bestub;
567  SCIP_Real bestslb;
568  SCIP_Real bestsub;
569  SCIP_Real rowcoef;
570  SCIP_Bool uselb;
571  int bestlbtype;
572  int bestubtype;
573  int bestslbtype;
574  int bestsubtype;
575  int probidx;
576 
577  bestlb = -SCIPinfinity(scip);
578  bestub = SCIPinfinity(scip);
579  bestlbtype = -3;
580  bestubtype = -3;
581 
582  var = SCIPcolGetVar(nonzcols[nonzcolsnonbinary[c]]);
583  probidx = SCIPvarGetProbindex(var);
584  rowcoef = rowweight * scale * nonzcoefs[nonzcolsnonbinary[c]];
585 
586  assert(assoctransvars[probidx] == -1);
587  assert(boundsfortrans[probidx] == -3);
588  assert(!SCIPisZero(scip, rowcoef));
589 
590  /* get closest simple lower bound and closest simple upper bound */
591  getClosestLb(scip, var, ALLOWLOCAL, &bestslb, &bestslbtype);
592  getClosestUb(scip, var, ALLOWLOCAL, &bestsub, &bestsubtype);
593 
594  SCIPdebugMessage(" %d: %g <%s, idx=%d, lp=%g, [%g(%d),%g(%d)]>:\n", c, rowcoef, SCIPvarGetName(var), probidx,
595  varsolvals[probidx], bestslb, bestslbtype, bestsub, bestsubtype);
596 
597  /* mixed integer set cannot be relaxed to 0-1 single node flow set because both simple bounds are -infinity
598  * and infinity, respectively
599  */
600  if( SCIPisInfinity(scip, -bestslb) && SCIPisInfinity(scip, bestsub) )
601  {
602  assert(!(*success));
603  goto TERMINATE;
604  }
605 
606  /* get closest lower bound that can be used to define the real variable y'_j in the 0-1 single node flow
607  * relaxation
608  */
609  if( !SCIPisInfinity(scip, bestsub) )
610  {
611  bestlb = bestslb;
612  bestlbtype = bestslbtype;
613 
614  if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
615  {
616  SCIP_Real bestvlb;
617  int bestvlbidx;
618 
619  SCIP_CALL( getClosestVlb(scip, var, bestsub, rowcoef, rowcoefsbinary, varsolvals, assoctransvars, &bestvlb, &bestvlbidx) );
620  if( SCIPisGT(scip, bestvlb, bestlb) )
621  {
622  bestlb = bestvlb;
623  bestlbtype = bestvlbidx;
624  }
625  }
626  }
627  /* get closest upper bound that can be used to define the real variable y'_j in the 0-1 single node flow
628  * relaxation
629  */
630  if( !SCIPisInfinity(scip, -bestslb) )
631  {
632  bestub = bestsub;
633  bestubtype = bestsubtype;
634 
635  if( SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
636  {
637  SCIP_Real bestvub;
638  int bestvubidx;
639 
640  SCIP_CALL( getClosestVub(scip, var, bestslb, rowcoef, rowcoefsbinary, varsolvals, assoctransvars, &bestvub, &bestvubidx) );
641  if( SCIPisLT(scip, bestvub, bestub) )
642  {
643  bestub = bestvub;
644  bestubtype = bestvubidx;
645  }
646  }
647  }
648  SCIPdebugMessage(" bestlb=%g(%d), bestub=%g(%d)\n", bestlb, bestlbtype, bestub, bestubtype);
649 
650  /* mixed integer set cannot be relaxed to 0-1 single node flow set because there are no suitable bounds
651  * to define the transformed variable y'_j
652  */
653  if( SCIPisInfinity(scip, -bestlb) && SCIPisInfinity(scip, bestub) )
654  {
655  assert(!(*success));
656  goto TERMINATE;
657  }
658 
659  /* select best upper bound if it is closer to the LP value of y_j and best lower bound otherwise and use this bound
660  * to define the real variable y'_j with 0 <= y'_j <= u'_j x_j in the 0-1 single node flow relaxation;
661  * prefer variable bounds
662  */
663  if( SCIPisEQ(scip, varsolvals[probidx], (1.0 - BOUNDSWITCH) * bestlb + BOUNDSWITCH * bestub) && bestlbtype >= 0 )
664  uselb = TRUE;
665  else if( SCIPisEQ(scip, varsolvals[probidx], (1.0 - BOUNDSWITCH) * bestlb + BOUNDSWITCH * bestub)
666  && bestubtype >= 0 )
667  uselb = FALSE;
668  else if( SCIPisLE(scip, varsolvals[probidx], (1.0 - BOUNDSWITCH) * bestlb + BOUNDSWITCH * bestub) )
669  uselb = TRUE;
670  else
671  {
672  assert(SCIPisGT(scip, varsolvals[probidx], (1.0 - BOUNDSWITCH) * bestlb + BOUNDSWITCH * bestub));
673  uselb = FALSE;
674  }
675  if( uselb )
676  {
677  /* use bestlb to define y'_j */
678 
679  assert(!SCIPisInfinity(scip, bestsub));
680  assert(!SCIPisInfinity(scip, - bestlb));
681  assert(bestsubtype == -1 || bestsubtype == -2);
682  assert(bestlbtype > -3 && bestlbtype < SCIPvarGetNVlbs(var));
683 
684  /* store for y_j that bestlb is the bound used to define y'_j and that y'_j is the associated real variable
685  * in the relaxed set
686  */
687  boundsfortrans[probidx] = bestlbtype;
688  boundtypesfortrans[probidx] = SCIP_BOUNDTYPE_LOWER;
689  assoctransvars[probidx] = *ntransvars;
690 
691  if( bestlbtype < 0 )
692  {
693  SCIP_Real val;
694  SCIP_Real contsolval;
695 
696  /* use simple lower bound in bestlb = l_j <= y_j <= u_j = bestsub to define
697  * y'_j = - a_j ( y_j - u_j ) with 0 <= y'_j <= a_j ( u_j - l_j ) x_j and x_j = 1 if a_j > 0
698  * y'_j = a_j ( y_j - u_j ) with 0 <= y'_j <= - a_j ( u_j - l_j ) x_j and x_j = 1 if a_j < 0,
699  * put j into the set
700  * N2 if a_j > 0
701  * N1 if a_j < 0
702  * and update the right hand side of the constraint in the relaxation
703  * rhs = rhs - a_j u_j
704  */
705  val = rowcoef * ( bestsub - bestlb );
706  contsolval = rowcoef * ( varsolvals[probidx] - bestsub );
707  if( SCIPisPositive(scip, rowcoef) )
708  {
709  transvarcoefs[*ntransvars] = - 1;
710  transvarvubcoefs[*ntransvars] = val;
711  transbinvarsolvals[*ntransvars] = 1.0;
712  transcontvarsolvals[*ntransvars] = - contsolval;
713  }
714  else
715  {
716  assert(SCIPisNegative(scip, rowcoef));
717  transvarcoefs[*ntransvars] = 1;
718  transvarvubcoefs[*ntransvars] = - val;
719  transbinvarsolvals[*ntransvars] = 1.0;
720  transcontvarsolvals[*ntransvars] = contsolval;
721  }
722  (*transrhs) -= (rowcoef * bestsub);
723 
724  SCIPdebugMessage(" --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
725  transvarcoefs[*ntransvars] == 1 ? "+" : "-", *ntransvars, *ntransvars, transvarvubcoefs[*ntransvars],
726  *ntransvars, *transrhs + (rowcoef * bestsub), rowcoef, bestsub, *transrhs);
727  }
728  else
729  {
730  SCIP_Real rowcoefbinary;
731  SCIP_Real varsolvalbinary;
732  SCIP_Real val;
733  SCIP_Real contsolval;
734  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(var);
735  SCIP_Real* vlbconsts = SCIPvarGetVlbConstants(var);
736  SCIP_Real* vlbcoefs = SCIPvarGetVlbCoefs(var);
737 
738  /* use variable lower bound in bestlb = l~_j x_j + d_j <= y_j <= u_j = bestsub to define
739  * y'_j = - ( a_j ( y_j - d_j ) + c_j x_j ) with 0 <= y'_j <= - ( a_j l~_j + c_j ) x_j if a_j > 0
740  * y'_j = a_j ( y_j - d_j ) + c_j x_j with 0 <= y'_j <= ( a_j l~_j + c_j ) x_j if a_j < 0,
741  * where c_j is the coefficient of x_j in the row, put j into the set
742  * N2 if a_j > 0
743  * N1 if a_j < 0
744  * and update the right hand side of the constraint in the relaxation
745  * rhs = rhs - a_j d_j
746  */
747  assert(SCIPvarIsBinary(vlbvars[bestlbtype]));
748 
749  rowcoefbinary = rowcoefsbinary[SCIPvarGetProbindex(vlbvars[bestlbtype])];
750  varsolvalbinary = varsolvals[SCIPvarGetProbindex(vlbvars[bestlbtype])];
751 
752  val = (rowcoef * vlbcoefs[bestlbtype]) + rowcoefbinary;
753  contsolval = (rowcoef * (varsolvals[probidx] - vlbconsts[bestlbtype])) + (rowcoefbinary * varsolvalbinary);
754 
755  if( SCIPisPositive(scip, rowcoef) )
756  {
757  transvarcoefs[*ntransvars] = - 1;
758  transvarvubcoefs[*ntransvars] = - val;
759  transbinvarsolvals[*ntransvars] = varsolvalbinary;
760  transcontvarsolvals[*ntransvars] = - contsolval;
761  }
762  else
763  {
764  assert(SCIPisNegative(scip, rowcoef));
765  transvarcoefs[*ntransvars] = 1;
766  transvarvubcoefs[*ntransvars] = val;
767  transbinvarsolvals[*ntransvars] = varsolvalbinary;
768  transcontvarsolvals[*ntransvars] = contsolval;
769  }
770  (*transrhs) -= (rowcoef * vlbconsts[bestlbtype]);
771 
772  /* store for x_j that y'_j is the associated real variable in the 0-1 single node flow relaxation */
773  assoctransvars[SCIPvarGetProbindex(vlbvars[bestlbtype])] = *ntransvars;
774 
775  SCIPdebugMessage(" --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=%s), rhs=%g-(%g*%g)=%g\n",
776  transvarcoefs[*ntransvars] == 1 ? "+" : "-", *ntransvars, *ntransvars, transvarvubcoefs[*ntransvars],
777  *ntransvars, SCIPvarGetName(vlbvars[bestlbtype]), *transrhs + (rowcoef * vlbconsts[bestlbtype]), rowcoef,
778  vlbconsts[bestlbtype], *transrhs );
779  }
780  }
781  else
782  {
783  /* use bestub to define y'_j */
784 
785  assert(!SCIPisInfinity(scip, bestub));
786  assert(!SCIPisInfinity(scip, - bestslb));
787  assert(bestslbtype == -1 || bestslbtype == -2);
788  assert(bestubtype > -3 && bestubtype < SCIPvarGetNVubs(var));
789 
790  /* store for y_j that bestub is the bound used to define y'_j and that y'_j is the associated real variable
791  * in the relaxed set
792  */
793  boundsfortrans[probidx] = bestubtype;
794  boundtypesfortrans[probidx] = SCIP_BOUNDTYPE_UPPER;
795  assoctransvars[probidx] = *ntransvars;
796 
797  if( bestubtype < 0 )
798  {
799  SCIP_Real val;
800  SCIP_Real contsolval;
801 
802  /* use simple upper bound in bestslb = l_j <= y_j <= u_j = bestub to define
803  * y'_j = a_j ( y_j - l_j ) with 0 <= y'_j <= a_j ( u_j - l_j ) x_j and x_j = 1 if a_j > 0
804  * y'_j = - a_j ( y_j - l_j ) with 0 <= y'_j <= - a_j ( u_j - l_j ) x_j and x_j = 1 if a_j < 0,
805  * put j into the set
806  * N1 if a_j > 0
807  * N2 if a_j < 0
808  * and update the right hand side of the constraint in the relaxation
809  * rhs = rhs - a_j l_j
810  */
811  val = rowcoef * ( bestub - bestslb );
812  contsolval = rowcoef * ( varsolvals[probidx] - bestslb );
813  if( SCIPisPositive(scip, rowcoef) )
814  {
815  transvarcoefs[*ntransvars] = 1;
816  transvarvubcoefs[*ntransvars] = val;
817  transbinvarsolvals[*ntransvars] = 1.0;
818  transcontvarsolvals[*ntransvars] = contsolval;
819  }
820  else
821  {
822  assert(SCIPisNegative(scip, rowcoef));
823  transvarcoefs[*ntransvars] = - 1;
824  transvarvubcoefs[*ntransvars] = - val;
825  transbinvarsolvals[*ntransvars] = 1.0;
826  transcontvarsolvals[*ntransvars] = - contsolval;
827  }
828  (*transrhs) -= (rowcoef * bestslb);
829 
830  SCIPdebugMessage(" --> bestub used for trans: ... %s y'_%d + ..., Y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
831  transvarcoefs[*ntransvars] == 1 ? "+" : "-", *ntransvars, *ntransvars, transvarvubcoefs[*ntransvars],
832  *ntransvars, *transrhs + (rowcoef * bestslb), rowcoef, bestslb, *transrhs);
833  }
834  else
835  {
836  SCIP_Real rowcoefbinary;
837  SCIP_Real varsolvalbinary;
838  SCIP_Real val;
839  SCIP_Real contsolval;
840 
841  SCIP_VAR** vubvars = SCIPvarGetVubVars(var);
842  SCIP_Real* vubconsts = SCIPvarGetVubConstants(var);
843  SCIP_Real* vubcoefs = SCIPvarGetVubCoefs(var);
844 
845  /* use variable upper bound in bestslb = l_j <= y_j <= u~_j x_j + d_j = bestub to define
846  * y'_j = a_j ( y_j - d_j ) + c_j x_j with 0 <= y'_j <= ( a_j u~_j + c_j ) x_j if a_j > 0
847  * y'_j = - ( a_j ( y_j - d_j ) + c_j x_j ) with 0 <= y'_j <= - ( a_j u~_j + c_j ) x_j if a_j < 0,
848  * where c_j is the coefficient of x_j in the row, put j into the set
849  * N1 if a_j > 0
850  * N2 if a_j < 0
851  * and update the right hand side of the constraint in the relaxation
852  * rhs = rhs - a_j d_j
853  */
854  assert(SCIPvarIsBinary(vubvars[bestubtype]));
855 
856  rowcoefbinary = rowcoefsbinary[SCIPvarGetProbindex(vubvars[bestubtype])];
857  varsolvalbinary = varsolvals[SCIPvarGetProbindex(vubvars[bestubtype])];
858 
859  val = ( rowcoef * vubcoefs[bestubtype] ) + rowcoefbinary;
860  contsolval = (rowcoef * (varsolvals[probidx] - vubconsts[bestubtype])) + (rowcoefbinary * varsolvalbinary);
861 
862  if( SCIPisPositive(scip, rowcoef) )
863  {
864  transvarcoefs[*ntransvars] = 1;
865  transvarvubcoefs[*ntransvars] = val;
866  transbinvarsolvals[*ntransvars] = varsolvalbinary;
867  transcontvarsolvals[*ntransvars] = contsolval;
868  }
869  else
870  {
871  assert(SCIPisNegative(scip, rowcoef));
872  transvarcoefs[*ntransvars] = - 1;
873  transvarvubcoefs[*ntransvars] = - val;
874  transbinvarsolvals[*ntransvars] = varsolvalbinary;
875  transcontvarsolvals[*ntransvars] = - contsolval;
876  }
877  (*transrhs) -= (rowcoef * vubconsts[bestubtype]);
878 
879  /* store for x_j that y'_j is the associated real variable in the 0-1 single node flow relaxation */
880  assoctransvars[SCIPvarGetProbindex(vubvars[bestubtype])] = *ntransvars;
881 
882  SCIPdebugMessage(" --> bestub used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=%s), rhs=%g-(%g*%g)=%g\n",
883  transvarcoefs[*ntransvars] == 1 ? "+" : "-", *ntransvars, *ntransvars, transvarvubcoefs[*ntransvars],
884  *ntransvars, SCIPvarGetName(vubvars[bestubtype]), *transrhs + (rowcoef * vubconsts[bestubtype]), rowcoef,
885  vubconsts[bestubtype], *transrhs);
886  }
887  }
888 
889  /* relaxing the mixed integer set to a 0-1 single node flow set was not successful because coefficient of y_j and
890  * the bounds selected for the transformation together result in an infinite variable upper bound in the 0-1 single
891  * node flow set; this can be caused by huge but finite values for the bounds or the coefficient
892  */
893  if( SCIPisInfinity(scip, transvarvubcoefs[*ntransvars]) )
894  {
895  assert(!(*success));
896  goto TERMINATE;
897  }
898 
899  assert(boundsfortrans[probidx] > -3);
900  assert(assoctransvars[probidx] >= 0 && assoctransvars[probidx] == (*ntransvars));
901  assert(transvarcoefs[*ntransvars] == 1 || transvarcoefs[*ntransvars] == - 1 );
902  assert(SCIPisFeasGE(scip, transbinvarsolvals[*ntransvars], 0.0)
903  && SCIPisFeasLE(scip, transbinvarsolvals[*ntransvars], 1.0));
904  assert(SCIPisFeasGE(scip, transvarvubcoefs[*ntransvars], 0.0)
905  && !SCIPisInfinity(scip, transvarvubcoefs[*ntransvars]));
906 
907  /* updates number of variables in transformed problem */
908  (*ntransvars)++;
909  }
910  assert(*ntransvars == nnonzcolsnonbinary);
911 
912  SCIPdebugMessage("transformation for BINARY variables (nbinvars=%d):\n", nnonzcolsbinary);
913 
914  /* binary variables not involved in used variable bounds imposed on non-binary variable */
915  for( c = 0; c < nnonzcolsbinary; c++ )
916  {
917  SCIP_VAR* var;
918  SCIP_Real rowcoef;
919  int probidx;
920  SCIP_Real val;
921  SCIP_Real contsolval;
922 
923  var = SCIPcolGetVar(nonzcols[nonzcolsbinary[c]]);
924  probidx = SCIPvarGetProbindex(var);
925  rowcoef = rowweight * scale * nonzcoefs[nonzcolsbinary[c]];
926 
927  assert(rowcoefsbinary[probidx] == rowcoef); /*lint !e777*/
928  assert(!SCIPisZero(scip, rowcoef));
929 
930  SCIPdebugMessage(" %d: %g <%s, idx=%d, lp=%g, [%g, %g]>:\n", c, rowcoef, SCIPvarGetName(var), probidx, varsolvals[probidx],
931  SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
932 
933  /* x_j has already been handled in connection with a non-binary variable */
934  if( assoctransvars[probidx] > -1 )
935  {
936  assert(assoctransvars[probidx] >= 0 && assoctransvars[probidx] <= nnonzcolsnonbinary);
937  assert(boundsfortrans[probidx] == -3);
938  SCIPdebugMessage(" --> already handled\n");
939  continue;
940  }
941 
942  assert(assoctransvars[probidx] == -1);
943  assert(boundsfortrans[probidx] == -3);
944 
945  /* store for x_j that y'_j is the associated real variable in the 0-1 single node flow relaxation */
946  assoctransvars[probidx] = *ntransvars;
947 
948  /* define
949  * y'_j = c_j x_j with 0 <= y'_j <= c_j x_j if c_j > 0
950  * y'_j = - c_j x_j with 0 <= y'_j <= - c_j x_j if c_j < 0,
951  * where c_j is the coefficient of x_j in the row and put j into the set
952  * N1 if c_j > 0
953  * N2 if c_j < 0.
954  */
955  val = rowcoef;
956  contsolval = rowcoef * varsolvals[probidx];
957  if( SCIPisPositive(scip, rowcoef) )
958  {
959  transvarcoefs[*ntransvars] = 1;
960  transvarvubcoefs[*ntransvars] = val;
961  transbinvarsolvals[*ntransvars] = varsolvals[probidx];
962  transcontvarsolvals[*ntransvars] = contsolval;
963  }
964  else
965  {
966  assert(SCIPisNegative(scip, rowcoef));
967  transvarcoefs[*ntransvars] = - 1;
968  transvarvubcoefs[*ntransvars] = - val;
969  transbinvarsolvals[*ntransvars] = varsolvals[probidx];
970  transcontvarsolvals[*ntransvars] = - contsolval;
971  }
972  assert(assoctransvars[probidx] >= 0 && assoctransvars[probidx] == (*ntransvars));
973  assert(transvarcoefs[*ntransvars] == 1 || transvarcoefs[*ntransvars] == - 1 );
974  assert(SCIPisFeasGE(scip, transbinvarsolvals[*ntransvars], 0.0)
975  && SCIPisFeasLE(scip, transbinvarsolvals[*ntransvars], 1.0));
976  assert(SCIPisFeasGE(scip, transvarvubcoefs[*ntransvars], 0.0)
977  && !SCIPisInfinity(scip, transvarvubcoefs[*ntransvars]));
978 
979  SCIPdebugMessage(" --> ... %s y'_%d + ..., y'_%d <= %g x_%d (=%s))\n",
980  transvarcoefs[*ntransvars] == 1 ? "+" : "-", *ntransvars, *ntransvars,
981  transvarvubcoefs[*ntransvars], *ntransvars, SCIPvarGetName(var) );
982 
983  /* updates number of variables in transformed problem */
984  (*ntransvars)++;
985  }
986  assert(*ntransvars >= nnonzcolsnonbinary && *ntransvars <= nnonzcols);
987 
988  /* construction was successful */
989  *success = TRUE;
990 
991 #ifdef SCIP_DEBUG
992  SCIPdebugMessage("constraint in constructed 0-1 single node flow relaxation: ");
993  for( c = 0; c < *ntransvars; c++ )
994  {
995  SCIPdebugPrintf("%s y'_%d ", transvarcoefs[c] == 1 ? "+" : "-", c);
996  }
997  SCIPdebugPrintf("<= %g\n", *transrhs);
998 #endif
999 
1000  TERMINATE:
1001  /* free data structures */
1002  SCIPfreeBufferArray(scip, &rowcoefsbinary);
1003  SCIPfreeBufferArray(scip, &nonzcolsnonbinary);
1004  SCIPfreeBufferArray(scip, &nonzcolsbinary);
1005 
1006  return SCIP_OKAY;
1007 }
1008 
1009 /** solve knapsack problem in maximization form with "<" constraint approximately by greedy; if needed, one can provide
1010  * arrays to store all selected items and all not selected items
1011  */
1012 static
1013 SCIP_RETCODE SCIPsolveKnapsackApproximatelyLT(
1014  SCIP* scip, /**< SCIP data structure */
1015  int nitems, /**< number of available items */
1016  SCIP_Real* weights, /**< item weights */
1017  SCIP_Real* profits, /**< item profits */
1018  SCIP_Real capacity, /**< capacity of knapsack */
1019  int* items, /**< item numbers */
1020  int* solitems, /**< array to store items in solution, or NULL */
1021  int* nonsolitems, /**< array to store items not in solution, or NULL */
1022  int* nsolitems, /**< pointer to store number of items in solution, or NULL */
1023  int* nnonsolitems, /**< pointer to store number of items not in solution, or NULL */
1024  SCIP_Real* solval /**< pointer to store optimal solution value, or NULL */
1025  )
1026 {
1027  SCIP_Real* tempsort;
1028  SCIP_Real solitemsweight;
1029  int j;
1030  int i;
1031 
1032  assert(weights != NULL);
1033  assert(profits != NULL);
1034  assert(SCIPisFeasGE(scip, capacity, 0.0));
1035  assert(!SCIPisInfinity(scip, capacity));
1036  assert(items != NULL);
1037  assert(nitems >= 0);
1038 
1039  if( solitems != NULL )
1040  {
1041  *nsolitems = 0;
1042  *nnonsolitems = 0;
1043  }
1044  if( solval != NULL )
1045  *solval = 0.0;
1046 
1047  /* allocate memory for temporary array used for sorting; array should contain profits divided by corresponding weights (p_1 / w_1 ... p_n / w_n )*/
1048  SCIP_CALL( SCIPallocBufferArray(scip, &tempsort, nitems) );
1049  /* initialize temporary array */
1050  for( i = nitems - 1; i >= 0; --i )
1051  {
1052  tempsort[i] = profits[i] / weights [i];
1053  }
1054 
1055  /* sort tempsort, items, weights and profits such that p_1 / w_1 >= p_2 / w_2 >= ... >= p_n / w_n */
1056  SCIPsortDownRealRealRealInt ( tempsort, weights, profits, items, nitems);
1057 
1058  /* free temporary array */
1059  SCIPfreeBufferArray(scip, &tempsort);
1060 
1061  /* select items as long as they fit into the knapsack */
1062  solitemsweight = 0.0;
1063  for( j = 0; j < nitems && SCIPisFeasLT(scip, solitemsweight + weights[j], capacity); j++ )
1064  {
1065  if( solitems != NULL )
1066  {
1067  solitems[*nsolitems] = items[j];
1068  (*nsolitems)++;
1069  }
1070  if( solval != NULL )
1071  (*solval) += profits[j];
1072  solitemsweight += weights[j];
1073  }
1074  for( ; j < nitems && solitems != NULL; j++ )
1075  {
1076  nonsolitems[*nnonsolitems] = items[j];
1077  (*nnonsolitems)++;
1078  }
1079 
1080  return SCIP_OKAY;
1081 }
1082 
1083 /** checks, whether the given scalar scales the given value to an integral number with error in the given bounds */
1084 static
1085 SCIP_Bool isIntegralScalar(
1086  SCIP_Real val, /**< value that should be scaled to an integral value */
1087  SCIP_Real scalar, /**< scalar that should be tried */
1088  SCIP_Real mindelta, /**< minimal relative allowed difference of scaled coefficient s*c and integral i */
1089  SCIP_Real maxdelta /**< maximal relative allowed difference of scaled coefficient s*c and integral i */
1090  )
1091 {
1092  SCIP_Real sval;
1093  SCIP_Real downval;
1094  SCIP_Real upval;
1095 
1096  assert(mindelta <= 0.0);
1097  assert(maxdelta >= 0.0);
1098 
1099  sval = val * scalar;
1100  downval = floor(sval);
1101  upval = ceil(sval);
1102 
1103  return (SCIPrelDiff(sval, downval) <= maxdelta || SCIPrelDiff(sval, upval) >= mindelta);
1104 }
1105 
1106 /** get integral number with error in the bounds which corresponds to given value scaled by a given scalar;
1107  * should be used in connection with isIntegralScalar()
1108  */
1109 static
1110 SCIP_Longint getIntegralVal(
1111  SCIP_Real val, /**< value that should be scaled to an integral value */
1112  SCIP_Real scalar, /**< scalar that should be tried */
1113  SCIP_Real mindelta, /**< minimal relative allowed difference of scaled coefficient s*c and integral i */
1114  SCIP_Real maxdelta /**< maximal relative allowed difference of scaled coefficient s*c and integral i */
1115  )
1116 {
1117  SCIP_Real sval;
1118  SCIP_Real upval;
1119  SCIP_Longint intval;
1120 
1121  assert(mindelta <= 0.0);
1122  assert(maxdelta >= 0.0);
1123 
1124  sval = val * scalar;
1125  upval = ceil(sval);
1126 
1127  if( SCIPrelDiff(sval, upval) >= mindelta )
1128  intval = (SCIP_Longint) upval;
1129  else
1130  intval = (SCIP_Longint) (floor(sval));
1131 
1132  return intval;
1133 }
1134 
1135 /** build the flow cover which corresponds to the given exact or approximate solution of KP^SNF; given unfinished
1136  * flow cover contains variables which have been fixed in advance
1137  */
1138 static
1139 void buildFlowCover(
1140  SCIP* scip, /**< SCIP data structure */
1141  int* coefs, /**< coefficient of all real variables in N1&N2 */
1142  SCIP_Real* vubcoefs, /**< coefficient in vub of all real variables in N1&N2 */
1143  SCIP_Real rhs, /**< right hand side of 0-1 single node flow constraint */
1144  int* solitems, /**< items in knapsack */
1145  int* nonsolitems, /**< items not in knapsack */
1146  int nsolitems, /**< number of items in knapsack */
1147  int nnonsolitems, /**< number of items not in knapsack */
1148  int* nflowcovervars, /**< pointer to store number of variables in flow cover */
1149  int* nnonflowcovervars, /**< pointer to store number of variables not in flow cover */
1150  int* flowcoverstatus, /**< pointer to store whether variable is in flow cover (+1) or not (-1) */
1151  SCIP_Real* flowcoverweight, /**< pointer to store weight of flow cover */
1152  SCIP_Real* lambda /**< pointer to store lambda */
1153 )
1154 {
1155  int j;
1156 
1157  assert(scip != NULL);
1158  assert(coefs != NULL);
1159  assert(vubcoefs != NULL);
1160  assert(solitems != NULL);
1161  assert(nonsolitems != NULL);
1162  assert(nsolitems >= 0);
1163  assert(nnonsolitems >= 0);
1164  assert(nflowcovervars != NULL && *nflowcovervars >= 0);
1165  assert(nnonflowcovervars != NULL && *nnonflowcovervars >= 0);
1166  assert(flowcoverstatus != NULL);
1167  assert(flowcoverweight != NULL);
1168  assert(lambda != NULL);
1169 
1170  /* get flowcover status for each item */
1171  for( j = 0; j < nsolitems; j++ )
1172  {
1173  /* j in N1 with z°_j = 1 => j in N1\C1 */
1174  if( coefs[solitems[j]] == 1 )
1175  {
1176  flowcoverstatus[solitems[j]] = -1;
1177  (*nnonflowcovervars)++;
1178  }
1179  /* j in N2 with z_j = 1 => j in C2 */
1180  else
1181  {
1182  assert(coefs[solitems[j]] == -1);
1183  flowcoverstatus[solitems[j]] = 1;
1184  (*nflowcovervars)++;
1185  (*flowcoverweight) -= vubcoefs[solitems[j]];
1186  }
1187  }
1188  for( j = 0; j < nnonsolitems; j++ )
1189  {
1190  /* j in N1 with z°_j = 0 => j in C1 */
1191  if( coefs[nonsolitems[j]] == 1 )
1192  {
1193  flowcoverstatus[nonsolitems[j]] = 1;
1194  (*nflowcovervars)++;
1195  (*flowcoverweight) += vubcoefs[nonsolitems[j]];
1196  }
1197  /* j in N2 with z_j = 0 => j in N2\C2 */
1198  else
1199  {
1200  assert(coefs[nonsolitems[j]] == -1);
1201  flowcoverstatus[nonsolitems[j]] = -1;
1202  (*nnonflowcovervars)++;
1203  }
1204  }
1205 
1206  /* get lambda = sum_{j in C1} u_j - sum_{j in C2} u_j - rhs */
1207  *lambda = (*flowcoverweight) - rhs;
1208 }
1209 
1210 /** get a flow cover (C1, C2) for a given 0-1 single node flow set
1211  * {(x,y) in {0,1}^n x R^n : sum_{j in N1} y_j - sum_{j in N2} y_j <= b, 0 <= y_j <= u_j x_j},
1212  * i.e., get sets C1 subset N1 and C2 subset N2 with sum_{j in C1} u_j - sum_{j in C2} u_j = b + lambda and lambda > 0
1213  */
1214 static
1215 SCIP_RETCODE getFlowCover(
1216  SCIP* scip, /**< SCIP data structure */
1217  int* coefs, /**< coefficient of all real variables in N1&N2 */
1218  SCIP_Real* solvals, /**< LP solution value of binary variable in vub of all real vars in N1&N2 */
1219  SCIP_Real* vubcoefs, /**< coefficient in vub of all real variables in N1&N2 */
1220  int nvars, /**< number of real variables in N1&N2 */
1221  SCIP_Real rhs, /**< right hand side of 0-1 single node flow constraint */
1222  int* nflowcovervars, /**< pointer to store number of variables in flow cover */
1223  int* nnonflowcovervars, /**< pointer to store number of variables not in flow cover */
1224  int* flowcoverstatus, /**< pointer to store whether variable is in flow cover (+1) or not (-1) */
1225  SCIP_Real* lambda, /**< pointer to store lambda */
1226  SCIP_Bool* found /**< pointer to store whether a cover was found */
1227  )
1228 {
1229  SCIP_Real* transprofitsint;
1230  SCIP_Real* transprofitsreal;
1231  SCIP_Real* transweightsreal;
1232  SCIP_Longint* transweightsint;
1233  int* items;
1234  int* itemsint;
1235  int* nonsolitems;
1236  int* solitems;
1237  SCIP_Real flowcoverweight;
1238  SCIP_Real flowcoverweightafterfix;
1239  SCIP_Real n1itemsweight;
1240  SCIP_Real n2itemsminweight;
1241  SCIP_Real scalar;
1242  SCIP_Real transcapacityreal;
1243 #if !defined(NDEBUG) || defined(SCIP_DEBUG)
1244  SCIP_Bool kpexact;
1245 #endif
1246  SCIP_Bool scalesuccess;
1247  SCIP_Bool transweightsrealintegral;
1248  SCIP_Longint transcapacityint;
1249  int nflowcovervarsafterfix;
1250  int nitems;
1251  int nn1items;
1252  int nnonflowcovervarsafterfix;
1253  int nnonsolitems;
1254  int nsolitems;
1255  int j;
1256 
1257  assert(scip != NULL);
1258  assert(coefs != NULL);
1259  assert(solvals != NULL);
1260  assert(vubcoefs != NULL);
1261  assert(nvars > 0);
1262  assert(nflowcovervars != NULL);
1263  assert(nnonflowcovervars != NULL);
1264  assert(flowcoverstatus != NULL);
1265  assert(lambda != NULL);
1266  assert(found != NULL);
1267 
1268  SCIPdebugMessage("--------------------- get flow cover ----------------------------------------------------\n");
1269 
1270  /* get data structures */
1271  SCIP_CALL( SCIPallocBufferArray(scip, &items, nvars) );
1272  SCIP_CALL( SCIPallocBufferArray(scip, &itemsint, nvars) );
1273  SCIP_CALL( SCIPallocBufferArray(scip, &transprofitsreal, nvars) );
1274  SCIP_CALL( SCIPallocBufferArray(scip, &transprofitsint, nvars) );
1275  SCIP_CALL( SCIPallocBufferArray(scip, &transweightsreal, nvars) );
1276  SCIP_CALL( SCIPallocBufferArray(scip, &transweightsint, nvars) );
1277  SCIP_CALL( SCIPallocBufferArray(scip, &solitems, nvars) );
1278  SCIP_CALL( SCIPallocBufferArray(scip, &nonsolitems, nvars) );
1279 
1280  BMSclearMemoryArray(flowcoverstatus, nvars);
1281  *found = FALSE;
1282  *nflowcovervars = 0;
1283  *nnonflowcovervars = 0;
1284  flowcoverweight = 0.0;
1285  nflowcovervarsafterfix = 0;
1286  nnonflowcovervarsafterfix = 0;
1287  flowcoverweightafterfix = 0.0;
1288 #if !defined(NDEBUG) || defined(SCIP_DEBUG)
1289  kpexact = FALSE;
1290 #endif
1291 
1292  /* fix some variables in advance according to the following fixing strategy
1293  * put j into N1\C1, if j in N1 and x*_j = 0,
1294  * put j into C1, if j in N1 and x*_j = 1,
1295  * put j into C2, if j in N2 and x*_j = 1,
1296  * put j into N2\C2, if j in N2 and x*_j = 0
1297  * and get the set of the remaining variables
1298  */
1299  SCIPdebugMessage("0. Fix some variables in advance:\n");
1300  nitems = 0;
1301  nn1items = 0;
1302  n1itemsweight = 0.0;
1303  n2itemsminweight = SCIP_REAL_MAX;
1304  for( j = 0; j < nvars; j++ )
1305  {
1306  assert(coefs[j] == 1 || coefs[j] == -1);
1307  assert(SCIPisFeasGE(scip, solvals[j], 0.0) && SCIPisFeasLE(scip, solvals[j], 1.0));
1308  assert(SCIPisFeasGE(scip, vubcoefs[j], 0.0));
1309 
1310  /* if u_j = 0, put j into N1\C1 and N2\C2, respectively */
1311  if( SCIPisFeasZero(scip, vubcoefs[j]) )
1312  {
1313  flowcoverstatus[j] = -1;
1314  (*nnonflowcovervars)++;
1315  continue;
1316  }
1317 
1318  /* x*_j is fractional */
1319  if( !SCIPisFeasIntegral(scip, solvals[j]) )
1320  {
1321  items[nitems] = j;
1322  nitems++;
1323  if( coefs[j] == 1 )
1324  {
1325  n1itemsweight += vubcoefs[j];
1326  nn1items++;
1327  }
1328  else
1329  n2itemsminweight = MIN(n2itemsminweight, vubcoefs[j]);
1330  }
1331  /* j is in N1 and x*_j = 0 */
1332  else if( coefs[j] == 1 && solvals[j] < 0.5 )
1333  {
1334  flowcoverstatus[j] = -1;
1335  (*nnonflowcovervars)++;
1336  SCIPdebugMessage(" <%d>: in N1-C1\n", j);
1337  }
1338  /* j is in N1 and x*_j = 1 */
1339  else if( coefs[j] == 1 && solvals[j] > 0.5 )
1340  {
1341  flowcoverstatus[j] = 1;
1342  (*nflowcovervars)++;
1343  flowcoverweight += vubcoefs[j];
1344  SCIPdebugMessage(" <%d>: in C1\n", j);
1345  }
1346  /* j is in N2 and x*_j = 1 */
1347  else if( coefs[j] == -1 && solvals[j] > 0.5 )
1348  {
1349  flowcoverstatus[j] = 1;
1350  (*nflowcovervars)++;
1351  flowcoverweight -= vubcoefs[j];
1352  SCIPdebugMessage(" <%d>: in C2\n", j);
1353  }
1354  /* j is in N2 and x*_j = 0 */
1355  else
1356  {
1357  assert(coefs[j] == -1 && solvals[j] < 0.5);
1358  flowcoverstatus[j] = -1;
1359  (*nnonflowcovervars)++;
1360  SCIPdebugMessage(" <%d>: in N2-C2\n", j);
1361  }
1362  }
1363  assert((*nflowcovervars) + (*nnonflowcovervars) + nitems == nvars);
1364  assert(nn1items >= 0);
1365 
1366  /* to find a flow cover, transform the following knapsack problem
1367  *
1368  * (KP^SNF) max sum_{j in N1} ( x*_j - 1 ) z_j + sum_{j in N2} x*_j z_j
1369  * sum_{j in N1} u_j z_j - sum_{j in N2} u_j z_j > b
1370  * z_j in {0,1} for all j in N1 & N2
1371  *
1372  * 1. to a knapsack problem in maximization form, such that all variables in the knapsack constraint have
1373  * positive weights and the constraint is a "<" constraint, by complementing all variables in N1
1374  *
1375  * (KP^SNF_rat) max sum_{j in N1} ( 1 - x*_j ) z°_j + sum_{j in N2} x*_j z_j
1376  * sum_{j in N1} u_j z°_j + sum_{j in N2} u_j z_j < - b + sum_{j in N1} u_j
1377  * z°_j in {0,1} for all j in N1
1378  * z_j in {0,1} for all j in N2,
1379  * and solve it approximately under consideration of the fixing,
1380  * or
1381  * 2. to a knapsack problem in maximization form, such that all variables in the knapsack constraint have
1382  * positive integer weights and the constraint is a "<=" constraint, by complementing all variables in N1
1383  * and multiplying the constraint by a suitable scalar C
1384  *
1385  * (KP^SNF_int) max sum_{j in N1} ( 1 - x*_j ) z°_j + sum_{j in N2} x*_j z_j
1386  * sum_{j in N1} C u_j z°_j + sum_{j in N2} C u_j z_j <= c
1387  * z°_j in {0,1} for all j in N1
1388  * z_j in {0,1} for all j in N2,
1389  * where
1390  * c = floor[ C (- b + sum_{j in N1} u_j ) ] if frac[ C (- b + sum_{j in N1} u_j ) ] > 0
1391  * c = C (- b + sum_{j in N1} u_j ) - 1 if frac[ C (- b + sum_{j in N1} u_j ) ] = 0
1392  * and solve it exactly under consideration of the fixing.
1393  */
1394  SCIPdebugMessage("1. Transform KP^SNF to KP^SNF_rat:\n");
1395  /* get weight and profit of variables in KP^SNF_rat and check, whether all weights are already integral */
1396  transweightsrealintegral = TRUE;
1397  for( j = 0; j < nitems; j++ )
1398  {
1399  transweightsreal[j] = vubcoefs[items[j]];
1400 
1401  if( !isIntegralScalar(transweightsreal[j], 1.0, -MINDELTA, MAXDELTA) )
1402  transweightsrealintegral = FALSE;
1403 
1404  if( coefs[items[j]] == 1 )
1405  {
1406  transprofitsreal[j] = 1.0 - solvals[items[j]];
1407  SCIPdebugMessage(" <%d>: j in N1: w_%d = %g, p_%d = %g %s\n", items[j], items[j], transweightsreal[j],
1408  items[j], transprofitsreal[j], SCIPisIntegral(scip, transweightsreal[j]) ? "" : " ----> NOT integral");
1409  }
1410  else
1411  {
1412  transprofitsreal[j] = solvals[items[j]];
1413  SCIPdebugMessage(" <%d>: j in N2: w_%d = %g, p_%d = %g %s\n", items[j], items[j], transweightsreal[j],
1414  items[j], transprofitsreal[j], SCIPisIntegral(scip, transweightsreal[j]) ? "" : " ----> NOT integral");
1415  }
1416  }
1417  /* get capacity of knapsack constraint in KP^SNF_rat */
1418  transcapacityreal = - rhs + flowcoverweight + n1itemsweight;
1419  SCIPdebugMessage(" transcapacity = -rhs(%g) + flowcoverweight(%g) + n1itemsweight(%g) = %g\n",
1420  rhs, flowcoverweight, n1itemsweight, transcapacityreal);
1421 
1422  /* there exists no flow cover if the capacity of knapsack constraint in KP^SNF_rat after fixing
1423  * is less than or equal to zero
1424  */
1425  if( SCIPisFeasLE(scip, transcapacityreal/10, 0.0) )
1426  {
1427  assert(!(*found));
1428  goto TERMINATE;
1429  }
1430 
1431  /* KP^SNF_rat has been solved by fixing some variables in advance */
1432  assert(nitems >= 0);
1433  if( nitems == 0)
1434  {
1435  /* get lambda = sum_{j in C1} u_j - sum_{j in C2} u_j - rhs */
1436  *lambda = flowcoverweight - rhs;
1437  *found = TRUE;
1438  goto TERMINATE;
1439  }
1440 
1441  /* Use the following strategy
1442  * solve KP^SNF_int exactly, if a suitable factor C is found and (nitems*capacity) <= MAXDYNPROGSPACE,
1443  * solve KP^SNF_rat approximately, otherwise
1444  */
1445 
1446  /* find a scaling factor C */
1447  if( transweightsrealintegral )
1448  {
1449  /* weights are already integral */
1450  scalar = 1.0;
1451  scalesuccess = TRUE;
1452  }
1453  else
1454  {
1455  scalesuccess = FALSE;
1456  SCIP_CALL( SCIPcalcIntegralScalar(transweightsreal, nitems, -MINDELTA, MAXDELTA, MAXDNOM, MAXSCALE, &scalar,
1457  &scalesuccess) );
1458  }
1459 
1460  /* initialize number of (non-)solution items, should be changed to a nonnegative number in all possible paths below */
1461  nsolitems = -1;
1462  nnonsolitems = -1;
1463 
1464  /* suitable factor C was found*/
1465  if( scalesuccess )
1466  {
1467  SCIP_Real tmp1;
1468  SCIP_Real tmp2;
1469 
1470  /* transform KP^SNF to KP^SNF_int */
1471  for( j = 0; j < nitems; ++j )
1472  {
1473  transweightsint[j] = getIntegralVal(transweightsreal[j], scalar, -MINDELTA, MAXDELTA);
1474  transprofitsint[j] = transprofitsreal[j];
1475  itemsint[j] = items[j];
1476  }
1477  if( isIntegralScalar(transcapacityreal, scalar, -MINDELTA, MAXDELTA) )
1478  {
1479  transcapacityint = getIntegralVal(transcapacityreal, scalar, -MINDELTA, MAXDELTA);
1480  transcapacityint -= 1;
1481  }
1482  else
1483  transcapacityint = (SCIP_Longint) (transcapacityreal * scalar);
1484  nflowcovervarsafterfix = *nflowcovervars;
1485  nnonflowcovervarsafterfix = *nnonflowcovervars;
1486  flowcoverweightafterfix = flowcoverweight;
1487 
1488  tmp1 = (SCIP_Real) (nitems + 1);
1489  tmp2 = (SCIP_Real) ((transcapacityint) + 1);
1490  if( transcapacityint * nitems <= MAXDYNPROGSPACE && tmp1 * tmp2 <= INT_MAX / 8.0)
1491  {
1492  SCIP_Bool success;
1493 
1494  /* solve KP^SNF_int by dynamic programming */
1495  SCIP_CALL(SCIPsolveKnapsackExactly(scip, nitems, transweightsint, transprofitsint, transcapacityint,
1496  itemsint, solitems, nonsolitems, &nsolitems, &nnonsolitems, NULL, &success));
1497 
1498  if( !success )
1499  {
1500  /* solve KP^SNF_rat approximately */
1501  SCIP_CALL(SCIPsolveKnapsackApproximatelyLT(scip, nitems, transweightsreal, transprofitsreal,
1502  transcapacityreal, items, solitems, nonsolitems, &nsolitems, &nnonsolitems, NULL));
1503  }
1504 #if !defined(NDEBUG) || defined(SCIP_DEBUG)
1505  else
1506  kpexact = TRUE;
1507 #endif
1508  }
1509  else
1510  {
1511  /* solve KP^SNF_rat approximately */
1512  SCIP_CALL(SCIPsolveKnapsackApproximatelyLT(scip, nitems, transweightsreal, transprofitsreal, transcapacityreal,
1513  items, solitems, nonsolitems, &nsolitems, &nnonsolitems, NULL));
1514  assert(!kpexact);
1515  }
1516  }
1517  else
1518  {
1519  /* solve KP^SNF_rat approximately */
1520  SCIP_CALL(SCIPsolveKnapsackApproximatelyLT(scip, nitems, transweightsreal, transprofitsreal, transcapacityreal,
1521  items, solitems, nonsolitems, &nsolitems, &nnonsolitems, NULL));
1522  assert(!kpexact);
1523  }
1524 
1525  assert(nsolitems != -1);
1526  assert(nnonsolitems != -1);
1527 
1528  /* build the flow cover from the solution of KP^SNF_rat and KP^SNF_int, respectively and the fixing */
1529  assert(*nflowcovervars + *nnonflowcovervars + nsolitems + nnonsolitems == nvars);
1530  buildFlowCover(scip, coefs, vubcoefs, rhs, solitems, nonsolitems, nsolitems, nnonsolitems, nflowcovervars,
1531  nnonflowcovervars, flowcoverstatus, &flowcoverweight, lambda);
1532  assert(*nflowcovervars + *nnonflowcovervars == nvars);
1533 
1534  /* if the found structure is not a flow cover, because of scaling, solve KP^SNF_rat approximately */
1535  if( SCIPisFeasLE(scip, *lambda, 0.0) )
1536  {
1537  assert(kpexact);
1538 
1539  /* solve KP^SNF_rat approximately */
1540  SCIP_CALL(SCIPsolveKnapsackApproximatelyLT(scip, nitems, transweightsreal, transprofitsreal, transcapacityreal,
1541  items, solitems, nonsolitems, &nsolitems, &nnonsolitems, NULL));
1542 #ifdef SCIP_DEBUG /* this time only for SCIP_DEBUG, because only then, the variable is used again */
1543  kpexact = FALSE;
1544 #endif
1545 
1546  /* build the flow cover from the solution of KP^SNF_rat and the fixing */
1547  *nflowcovervars = nflowcovervarsafterfix;
1548  *nnonflowcovervars = nnonflowcovervarsafterfix;
1549  flowcoverweight = flowcoverweightafterfix;
1550 
1551  assert(*nflowcovervars + *nnonflowcovervars + nsolitems + nnonsolitems == nvars);
1552  buildFlowCover(scip, coefs, vubcoefs, rhs, solitems, nonsolitems, nsolitems, nnonsolitems, nflowcovervars,
1553  nnonflowcovervars, flowcoverstatus, &flowcoverweight, lambda);
1554  assert(*nflowcovervars + *nnonflowcovervars == nvars);
1555  }
1556  *found = TRUE;
1557 
1558  TERMINATE:
1559  assert((!*found) || SCIPisFeasGT(scip, *lambda, 0.0));
1560 #ifdef SCIP_DEBUG
1561  if( *found )
1562  {
1563  SCIPdebugMessage("2. %s solution:\n", kpexact ? "exact" : "approximate");
1564  for( j = 0; j < nvars; j++ )
1565  {
1566  if( coefs[j] == 1 && flowcoverstatus[j] == 1 )
1567  {
1568  SCIPdebugMessage(" C1: + y_%d [u_%d = %g]\n", j, j, vubcoefs[j]);
1569  }
1570  else if( coefs[j] == -1 && flowcoverstatus[j] == 1 )
1571  {
1572  SCIPdebugMessage(" C2: - y_%d [u_%d = %g]\n", j, j, vubcoefs[j]);
1573  }
1574  }
1575  SCIPdebugMessage(" flowcoverweight(%g) = rhs(%g) + lambda(%g)\n", flowcoverweight, rhs, *lambda);
1576  }
1577 #endif
1578 
1579  /* free data structures */
1580  SCIPfreeBufferArray(scip, &nonsolitems);
1581  SCIPfreeBufferArray(scip, &solitems);
1582  SCIPfreeBufferArray(scip, &transweightsint);
1583  SCIPfreeBufferArray(scip, &transweightsreal);
1584  SCIPfreeBufferArray(scip, &transprofitsint);
1585  SCIPfreeBufferArray(scip, &transprofitsreal);
1586  SCIPfreeBufferArray(scip, &itemsint);
1587  SCIPfreeBufferArray(scip, &items);
1588 
1589  return SCIP_OKAY;
1590 }
1591 
1592 /** for a given flow cover and a given value of delta, choose L1 subset N1 \ C1 and L2 subset N2 \ C2 by comparison such that
1593  * the violation of the resulting c-MIRFCI is maximized.
1594  */
1595 static
1596 void getL1L2(
1597  SCIP* scip, /**< SCIP data structure */
1598  int ntransvars, /**< number of continuous variables in N1 & N2 */
1599  int* transvarcoefs, /**< coefficient of all continuous variables in N1 & N2 */
1600  SCIP_Real* transbinvarsolvals, /**< LP solution value of bin var in vub of all continuous vars in N1 & N2 */
1601  SCIP_Real* transcontvarsolvals,/**< LP solution value of all continuous vars in N1 & N2 */
1602  SCIP_Real* transvarvubcoefs, /**< coefficient of vub of all continuous variables in N1 & N2 */
1603  int* transvarflowcoverstatus,/**< pointer to store whether non-binary var is in L2 (2) or not (-1 or 1) */
1604  SCIP_Real delta, /**< delta */
1605  SCIP_Real lambda /**< lambda */
1606  )
1607 {
1608  SCIP_Real fbeta;
1609  SCIP_Real onedivoneminusfbeta;
1610  int j;
1611 
1612  assert(scip != NULL);
1613  assert(ntransvars >= 0);
1614  assert(transvarcoefs != NULL);
1615  assert(transbinvarsolvals != NULL);
1616  assert(transcontvarsolvals != NULL);
1617  assert(transvarvubcoefs != NULL);
1618  assert(transvarflowcoverstatus != NULL);
1619  assert(SCIPisGT(scip, delta, lambda));
1620  assert(SCIPisFeasGT(scip, lambda, 0.0));
1621 
1622  /* for beta = - lambda / delta with delta > lambda,
1623  * f_beta = beta - floor[beta] = 1 - ( lambda / delta )
1624  * 1 / ( 1 - f_beta ) = delta / lambda
1625  */
1626  fbeta = 1.0 - (lambda / delta);
1627  onedivoneminusfbeta = delta / lambda;
1628 
1629  SCIPdebugMessage(" --------------------- get L1 and L2 -----------------------------------------------------\n");
1630  SCIPdebugMessage(" L1 = { j in N1-C1 : y*_j >= ( u_j - lambda F_{f_beta}( u_j/delta) ) x*_j }\n");
1631  SCIPdebugMessage(" L2 = { j in N2-C2 : y*_j >= - lambda F_{f_beta}(- u_j/delta) x*_j }\n");
1632 
1633  /* set flowcover status of continuous variable x_j to 2, i.e., put j intp L1 and L2, respectively
1634  * if j is in N1\C1 and y*_j >= ( u_j - lambda F_{f_beta}( u_j/delta) ) x*_j
1635  * if j is in N2\C2 and y*_j >= - lambda F_{f_beta}(- u_j/delta) x*_j
1636  */
1637  for( j = 0; j < ntransvars; j++ )
1638  {
1639  SCIP_Real d;
1640  SCIP_Real downd;
1641  SCIP_Real fd;
1642  SCIP_Real mirval;
1643 
1644  assert(transvarcoefs[j] == 1 || transvarcoefs[j] == -1);
1645  assert(SCIPisGE(scip, transvarvubcoefs[j], 0.0));
1646  assert(SCIPisFeasGE(scip, transbinvarsolvals[j], 0.0) && SCIPisFeasLE(scip, transbinvarsolvals[j], 1.0));
1647 
1648  /* j in N1\C1 */
1649  if( transvarcoefs[j] == 1 && transvarflowcoverstatus[j] == -1 )
1650  {
1651  /* Let d = u_j/delta and alpha = f_beta, than the MIR function is defined as
1652  * F_{alpha}(d) = down(d) , if f_d <= alpha
1653  * F_{alpha}(d) = down(d) + (f_d - alpha)/(1 - alpha), if f_d > alpha
1654  */
1655  d = transvarvubcoefs[j] / delta;
1656  downd = SCIPfloor(scip, d);
1657  fd = d - downd;
1658  if( SCIPisSumLE(scip, fd, fbeta) )
1659  mirval = downd;
1660  else
1661  mirval = downd + ((fd - fbeta) * onedivoneminusfbeta);
1662 
1663  /* y*_j >= ( u_j - lambda F_{f_beta}(u_j/delta) ) x*_j */
1664  if( SCIPisFeasGE(scip, transcontvarsolvals[j], ( transvarvubcoefs[j] - ( lambda * mirval ) ) * transbinvarsolvals[j]) )
1665  transvarflowcoverstatus[j] = 2;
1666 
1667  SCIPdebugMessage(" <%d>: in N1-C1: %g ?>=? ( %g - %g F_{f_beta}(%g)(%g) ) %g = %g ---> fcstatus = %d\n",
1668  j, transcontvarsolvals[j], transvarvubcoefs[j], lambda, d, mirval, transbinvarsolvals[j],
1669  ( transvarvubcoefs[j] - ( lambda * mirval ) ) * transbinvarsolvals[j], transvarflowcoverstatus[j]);
1670  }
1671 
1672  /* j in N2\C2 */
1673  if( transvarcoefs[j] == -1 && transvarflowcoverstatus[j] == -1 )
1674  {
1675  /* Let d = - u_j/delta and alpha = f_beta, than the MIR function is defined as
1676  * F_{alpha}(d) = down(d) , if f_d <= alpha
1677  * F_{alpha}(d) = down(d) + (f_d - alpha)/(1 - alpha), if f_d > alpha
1678  */
1679  d = - transvarvubcoefs[j] / delta;
1680  downd = SCIPfloor(scip, d);
1681  fd = d - downd;
1682  if( SCIPisSumLE(scip, fd, fbeta) )
1683  mirval = downd;
1684  else
1685  mirval = downd + ((fd - fbeta) * onedivoneminusfbeta);
1686 
1687  /* j in N2\C2 and y*_j >= - lambda F_{f_beta}(- u_j/delta) x*_j */
1688  if( SCIPisFeasGE(scip, transcontvarsolvals[j], - ( lambda * mirval ) * transbinvarsolvals[j]) )
1689  transvarflowcoverstatus[j] = 2;
1690 
1691  SCIPdebugMessage(" <%d>: in N2-C2: %g ?>=? - %g F_{f_beta}(-%g)(%g) %g = %g ---> fcstatus = %d\n",
1692  j, transcontvarsolvals[j], lambda, d, mirval, transbinvarsolvals[j],
1693  - ( lambda * mirval ) * transbinvarsolvals[j], transvarflowcoverstatus[j]);
1694  }
1695  }
1696 }
1697 
1698 /** get for all problem variables with nonzero coefficient in current row the bound which should be used for the
1699  * substitution routine in the c-MIR routine; this bound depends on the bound used for each variable to get associated
1700  * variable in transformed problem
1701  */
1702 static
1703 SCIP_RETCODE getBoundsForSubstitution(
1704  SCIP* scip, /**< SCIP data structure */
1705  SCIP_VAR** vars, /**< active problem variables */
1706  int nvars, /**< number of active problem variables */
1707  int* boundsfortrans, /**< bound used for transformation for all non-binary vars of current row */
1708  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bound used for transform. for all non-binary vars of current row */
1709  int* assoctransvars, /**< associated var in transformed problem for all vars of current row */
1710  int* transvarcoefs, /**< coefficient of all vars in transformed problem */
1711  int* flowcoverstatus, /**< flow cover status of all non-binary vars in transformed problem;
1712  * 1 if in C1 & C2, 2 if in L2, -1 N1 \ C1 & N2 \ (C2&L2) */
1713  int ntransvars, /**< number of vars in transformed problem */
1714  int* boundsforsubst, /**< pointer to store bounds that should be used for substitution in c-mir for vars */
1715  SCIP_BOUNDTYPE* boundtypesforsubst /**< pointer to store types of bounds that should be used for substitution in c-mir */
1716  )
1717 {
1718  int j;
1719 
1720  assert(scip != NULL);
1721  assert(vars != NULL);
1722  assert(nvars >= 0);
1723  assert(boundsfortrans != NULL);
1724  assert(boundtypesfortrans != NULL);
1725  assert(assoctransvars != NULL);
1726  assert(transvarcoefs != NULL);
1727  assert(flowcoverstatus != NULL);
1728  assert(ntransvars >= 0 && ntransvars <= nvars);
1729  assert(boundsforsubst != NULL);
1730  assert(boundtypesforsubst != NULL);
1731 
1732  for( j = 0; j < nvars; j++ )
1733  {
1734  assert(SCIPvarGetProbindex(vars[j]) == j);
1735  assert(assoctransvars[j] == -1 || ( flowcoverstatus[assoctransvars[j]] == -1
1736  || flowcoverstatus[assoctransvars[j]] == 1 || flowcoverstatus[assoctransvars[j]] == 2 ));
1737 
1738  /* variable has no associated variable in transformed problem */
1739  if( assoctransvars[j] == -1 )
1740  {
1741  assert(boundsfortrans[j] == -3);
1742  boundsforsubst[j] = -3;
1743  continue;
1744  }
1745 
1746  /* binary variable */
1747  if( SCIPvarIsBinary(vars[j]) )
1748  {
1749  /* j in C1 & C2 */
1750  if( flowcoverstatus[assoctransvars[j]] == 1 )
1751  {
1752  boundsforsubst[j] = -1;
1753  boundtypesforsubst[j] = SCIP_BOUNDTYPE_UPPER;
1754  }
1755  /* j in N1\C1 & N2\C2 */
1756  else
1757  {
1758  boundsforsubst[j] = -1;
1759  boundtypesforsubst[j] = SCIP_BOUNDTYPE_LOWER;
1760  }
1761  }
1762  /* non-binary variables */
1763  else
1764  {
1765  /* j in C1 & C2 & L1 & L2 */
1766  if( flowcoverstatus[assoctransvars[j]] >= 1 )
1767  {
1768  boundsforsubst[j] = boundsfortrans[j];
1769  boundtypesforsubst[j] = boundtypesfortrans[j];
1770  }
1771  /* j in N1\C1 & N2\(C2 & L2) */
1772  else
1773  {
1774  if( boundtypesfortrans[j] == SCIP_BOUNDTYPE_UPPER )
1775  {
1776  SCIP_Real closestlb;
1777  int closestlbtype;
1778 
1779  getClosestLb(scip, vars[j], ALLOWLOCAL, &closestlb, &closestlbtype);
1780  assert(!SCIPisInfinity(scip, -closestlb));
1781  assert(closestlbtype == -1 || closestlbtype == -2);
1782 
1783  boundsforsubst[j] = closestlbtype;
1784  boundtypesforsubst[j] = SCIP_BOUNDTYPE_LOWER;
1785  }
1786  else
1787  {
1788  SCIP_Real closestub;
1789  int closestubtype;
1790 
1791  getClosestUb(scip, vars[j], ALLOWLOCAL, &closestub, &closestubtype);
1792  assert(!SCIPisInfinity(scip, closestub));
1793  assert(closestubtype == -1 || closestubtype == -2);
1794 
1795  boundsforsubst[j] = closestubtype;
1796  boundtypesforsubst[j] = SCIP_BOUNDTYPE_UPPER;
1797  }
1798  }
1799  }
1800  }
1801 
1802  return SCIP_OKAY;
1803 }
1804 
1805 /** stores nonzero elements of dense coefficient vector as sparse vector, and calculates activity and norm */
1806 static
1807 SCIP_RETCODE storeCutInArrays(
1808  SCIP* scip, /**< SCIP data structure */
1809  int nvars, /**< number of problem variables */
1810  SCIP_VAR** vars, /**< problem variables */
1811  SCIP_Real* cutcoefs, /**< dense coefficient vector */
1812  SCIP_Real* varsolvals, /**< dense variable LP solution vector */
1813  char normtype, /**< type of norm to use for efficacy norm calculation */
1814  SCIP_VAR** cutvars, /**< array to store variables of sparse cut vector */
1815  SCIP_Real* cutvals, /**< array to store coefficients of sparse cut vector */
1816  int* cutlen, /**< pointer to store number of nonzero entries in cut */
1817  SCIP_Real* cutact, /**< pointer to store activity of cut */
1818  SCIP_Real* cutnorm /**< pointer to store norm of cut vector */
1819  )
1820 {
1821  SCIP_Real val;
1822  SCIP_Real absval;
1823  SCIP_Real cutsqrnorm;
1824  SCIP_Real act;
1825  SCIP_Real norm;
1826  int len;
1827  int v;
1828 
1829  assert(nvars == 0 || cutcoefs != NULL);
1830  assert(nvars == 0 || varsolvals != NULL);
1831  assert(cutvars != NULL);
1832  assert(cutvals != NULL);
1833  assert(cutlen != NULL);
1834  assert(cutact != NULL);
1835  assert(cutnorm != NULL);
1836 
1837  len = 0;
1838  act = 0.0;
1839  norm = 0.0;
1840  switch( normtype )
1841  {
1842  case 'e':
1843  cutsqrnorm = 0.0;
1844  for( v = 0; v < nvars; ++v )
1845  {
1846  val = cutcoefs[v];
1847  if( !SCIPisZero(scip, val) )
1848  {
1849  act += val * varsolvals[v];
1850  cutsqrnorm += SQR(val);
1851  cutvars[len] = vars[v];
1852  cutvals[len] = val;
1853  len++;
1854  }
1855  }
1856  norm = SQRT(cutsqrnorm);
1857  break;
1858  case 'm':
1859  for( v = 0; v < nvars; ++v )
1860  {
1861  val = cutcoefs[v];
1862  if( !SCIPisZero(scip, val) )
1863  {
1864  act += val * varsolvals[v];
1865  absval = REALABS(val);
1866  norm = MAX(norm, absval);
1867  cutvars[len] = vars[v];
1868  cutvals[len] = val;
1869  len++;
1870  }
1871  }
1872  break;
1873  case 's':
1874  for( v = 0; v < nvars; ++v )
1875  {
1876  val = cutcoefs[v];
1877  if( !SCIPisZero(scip, val) )
1878  {
1879  act += val * varsolvals[v];
1880  norm += REALABS(val);
1881  cutvars[len] = vars[v];
1882  cutvals[len] = val;
1883  len++;
1884  }
1885  }
1886  break;
1887  case 'd':
1888  for( v = 0; v < nvars; ++v )
1889  {
1890  val = cutcoefs[v];
1891  if( !SCIPisZero(scip, val) )
1892  {
1893  act += val * varsolvals[v];
1894  norm = 1.0;
1895  cutvars[len] = vars[v];
1896  cutvals[len] = val;
1897  len++;
1898  }
1899  }
1900  break;
1901  default:
1902  SCIPerrorMessage("invalid efficacy norm parameter '%c'\n", normtype);
1903  return SCIP_INVALIDDATA;
1904  }
1905 
1906  *cutlen = len;
1907  *cutact = act;
1908  *cutnorm = norm;
1909 
1910  return SCIP_OKAY;
1911 }
1912 
1913 /** adds given cut to LP if violated */
1914 static
1915 SCIP_RETCODE addCut(
1916  SCIP* scip, /**< SCIP data structure */
1917  SCIP_SEPA* sepa, /**< separator */
1918  SCIP_SEPADATA* sepadata, /**< separator data */
1919  SCIP_VAR** vars, /**< problem variables */
1920  int nvars, /**< number of problem variables */
1921  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
1922  SCIP_Real* varsolvals, /**< solution values of active variables */
1923  SCIP_Real* cutcoefs, /**< coefficients of active variables in cut */
1924  SCIP_Real cutrhs, /**< right hand side of cut */
1925  SCIP_Bool cutislocal, /**< is the cut only locally valid? */
1926  int cutrank, /**< rank of the cut */
1927  char normtype, /**< type of norm to use for efficacy norm calculation */
1928  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
1929  int* ncuts /**< pointer to count the number of added cuts */
1930  )
1931 {
1932  SCIP_VAR** cutvars;
1933  SCIP_Real* cutvals;
1934  SCIP_Real cutact;
1935  SCIP_Real cutnorm;
1936  int cutlen;
1937  SCIP_Bool success;
1938 
1939  assert(scip != NULL);
1940  assert(varsolvals != NULL);
1941  assert(cutcoefs != NULL);
1942  assert(ncuts != NULL);
1943  assert(cutoff != NULL);
1944  assert(nvars == 0 || vars != NULL);
1945 
1946  *cutoff = FALSE;
1947 
1948  SCIPdebugMessage("--------------------- found cut ---------------------------------------------------------\n");
1949 
1950  /* gets temporary memory for storing the cut as sparse row */
1951  SCIP_CALL( SCIPallocBufferArray(scip, &cutvars, nvars) );
1952  SCIP_CALL( SCIPallocBufferArray(scip, &cutvals, nvars) );
1953 
1954  /* stores the cut as sparse row, calculates activity and norm of cut */
1955  SCIP_CALL( storeCutInArrays(scip, nvars, vars, cutcoefs, varsolvals, normtype,
1956  cutvars, cutvals, &cutlen, &cutact, &cutnorm) );
1957 
1958  if( SCIPisPositive(scip, cutnorm) && SCIPisEfficacious(scip, (cutact - cutrhs)/cutnorm) )
1959  {
1960  SCIP_ROW* cut;
1961  char cutname[SCIP_MAXSTRLEN];
1962 
1963  /* creates the cut */
1964  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "flowcover%d_%d", SCIPgetNLPs(scip), *ncuts);
1965  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs,
1966  cutislocal, FALSE, sepadata->dynamiccuts) );
1967  SCIP_CALL( SCIPaddVarsToRow(scip, cut, cutlen, cutvars, cutvals) );
1968 
1969  /* set cut rank */
1970  SCIProwChgRank(cut, cutrank);
1971 
1972  SCIPdebugMessage(" -> found potential flowcover cut <%s>: activity=%f, rhs=%f, norm=%f, eff=%f\n",
1973  cutname, cutact, cutrhs, cutnorm, SCIPgetCutEfficacy(scip, sol, cut));
1974  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, cut, NULL) ) );
1975 
1976 #if 0 /* tries to scale the cut to integral values */
1977  SCIP_CALL( SCIPmakeRowIntegral(scip, cut, -SCIPepsilon(scip), SCIPsumepsilon(scip),
1978  10, 100.0, MAKECONTINTEGRAL, &success) );
1979  if( success && !SCIPisCutEfficacious(scip, sol, cut) )
1980  {
1981  SCIPdebugMessage(" -> flowcover cut <%s> no longer efficacious: act=%f, rhs=%f, norm=%f, eff=%f\n",
1982  cutname, cutact, cutrhs, cutnorm, SCIPgetCutEfficacy(scip, sol, cut));
1983  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, cut, NULL) ) );
1984  success = FALSE;
1985  }
1986 #else
1987  success = TRUE;
1988 #endif
1989 
1990  /* if scaling was successful, adds the cut */
1991  if( success ) /*lint !e774*/ /* Boolean within 'if' always evaluates to True */
1992  {
1993  SCIPdebugMessage(" -> found flowcover cut <%s>: act=%f, rhs=%f, norm=%f, eff=%f, rank=%d, min=%f, max=%f (range=%g)\n",
1994  cutname, cutact, cutrhs, cutnorm, SCIPgetCutEfficacy(scip, sol, cut), SCIProwGetRank(cut),
1995  SCIPgetRowMinCoef(scip, cut), SCIPgetRowMaxCoef(scip, cut),
1996  SCIPgetRowMaxCoef(scip, cut)/SCIPgetRowMinCoef(scip, cut));
1997  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, cut, NULL) ) );
1998  SCIP_CALL( SCIPaddCut(scip, sol, cut, FALSE, cutoff) );
1999  if( !(*cutoff) && !cutislocal )
2000  {
2001  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
2002  }
2003  (*ncuts)++;
2004  }
2005 
2006  /* releases the row */
2007  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
2008  }
2009 
2010  /* frees temporary memory */
2011  SCIPfreeBufferArray(scip, &cutvals);
2012  SCIPfreeBufferArray(scip, &cutvars);
2013 
2014  return SCIP_OKAY;
2015 }
2016 
2017 /** calculates efficacy of the given cut */
2018 static
2019 SCIP_Real calcEfficacy(
2020  int nvars, /**< number of variables in the problem */
2021  SCIP_Real* cutcoefs, /**< dense vector of cut coefficients */
2022  SCIP_Real cutrhs, /**< right hand side of cut */
2023  SCIP_Real cutact /**< activity of cut */
2024  )
2025 {
2026  SCIP_Real sqrnorm;
2027  int i;
2028 
2029  assert(cutcoefs != NULL);
2030 
2031  sqrnorm = 0;
2032  for( i = 0; i < nvars; ++i )
2033  sqrnorm += SQR(cutcoefs[i]);
2034  sqrnorm = MAX(sqrnorm, 1e-06);
2035  return (cutact - cutrhs)/SQRT(sqrnorm);
2036 }
2037 
2038 /* generate c-MIRFCI for different sets L1 and L2 and different values of delta */
2039 static
2040 SCIP_RETCODE cutGenerationHeuristic(
2041  SCIP* scip, /**< SCIP data structure */
2042  SCIP_SEPA* sepa, /**< separator */
2043  SCIP_SEPADATA* sepadata, /**< separator data */
2044  SCIP_VAR** vars, /**< active problem variables */
2045  int nvars, /**< number of active problem variables */
2046  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2047  SCIP_Real* varsolvals, /**< solution values of active variables */
2048  SCIP_Real* rowweights, /**< weight of rows in aggregated row */
2049  SCIP_Real scalar, /**< additional scaling factor of rows in aggregation */
2050  int* boundsfortrans, /**< bound used for all non-bin vars of row */
2051  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bound used for all non-binary vars of row */
2052  int* assoctransvars, /**< associated var in relaxed set for all vars of row */
2053  int ntransvars, /**< number of real variables in N1&N2 */
2054  int* transvarcoefs, /**< coefficient of all continuous variables in N1 & N2 */
2055  SCIP_Real* transbinvarsolvals, /**< LP solution value of binary variable in vub of all real vars in N1&N2 */
2056  SCIP_Real* transcontvarsolvals,/**< LP solution value of all real vars in N1&N2 */
2057  SCIP_Real* transvarvubcoefs, /**< coefficient of vub of all continuous variables in N1 & N2 */
2058  int* transvarflowcoverstatus, /**< pointer to store whether non-binary var is in L2 (2) or not (-1 or 1) */
2059  SCIP_Real lambda, /**< lambda */
2060  char normtype, /**< type of norm to use for efficacy norm calculation */
2061  SCIP_Bool* cutoff, /**< whether a cutoff has been detected */
2062  int* ncuts /**< pointer to count the number of generated cuts */
2063  )
2064 {
2065  SCIP_BOUNDTYPE* boundtypesforsubst;
2066  SCIP_Real* candsetdelta;
2067  SCIP_Real* cutcoefs;
2068  SCIP_Real* testeddeltas;
2069  int* boundsforsubst;
2070  int* transvarflowcoverstatustmp;
2071  SCIP_Real bestdelta;
2072  SCIP_Real bestefficacy;
2073  SCIP_Real c1vubcoefsmax;
2074  SCIP_Real cutrhs;
2075  SCIP_Real cutact;
2076  SCIP_Real l2tmpvubcoefsmax;
2077  SCIP_Real nvubcoefsmax;
2078  SCIP_Bool cutislocal;
2079  SCIP_Bool success;
2080  int cutrank;
2081  int ncandsetdelta;
2082  int ntesteddeltas;
2083  int startidx;
2084  int j;
2085 
2086  assert(scip != NULL);
2087  assert(sepadata != NULL);
2088  assert(vars != NULL);
2089  assert(varsolvals != NULL);
2090  assert(rowweights != NULL);
2091  assert(boundsfortrans != NULL);
2092  assert(boundtypesfortrans != NULL);
2093  assert(assoctransvars != NULL);
2094  assert(ntransvars >=0);
2095  assert(transvarcoefs != NULL);
2096  assert(transbinvarsolvals != NULL);
2097  assert(transcontvarsolvals != NULL);
2098  assert(transvarvubcoefs != NULL);
2099  assert(transvarflowcoverstatus != NULL);
2100  assert(SCIPisFeasGT(scip, lambda, 0.0));
2101  assert(ncuts != NULL);
2102 
2103  SCIPdebugMessage("--------------------- cut generation heuristic ------------------------------------------\n");
2104 
2105  /* get data structures */
2106  SCIP_CALL( SCIPallocBufferArray(scip, &testeddeltas, ntransvars + 4) );
2107  SCIP_CALL( SCIPallocBufferArray(scip, &candsetdelta, ntransvars + 4) );
2108  SCIP_CALL( SCIPallocBufferArray(scip, &transvarflowcoverstatustmp, ntransvars) );
2109  SCIP_CALL( SCIPallocBufferArray(scip, &boundsforsubst, nvars) );
2110  SCIP_CALL( SCIPallocBufferArray(scip, &boundtypesforsubst, nvars) );
2111  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
2112 
2113  SCIPdebugMessage("1. get candidate set for the value of delta:\n");
2114 
2115  /* get candidate set for the value of delta
2116  * N* = { u_j : j in N and u_j > lambda} & { max{ u_j : j in N and u_j >= lambda } + 1, lambda + 1 };
2117  * store the following values which will be tested in any case at the beginning of the array
2118  * max{ u_j : j in N and u_j >= lambda } + 1,
2119  * lambda + 1,
2120  * max{ u_j : j in C1 and u_j > lambda },
2121  * max{ u_j : j in L~2 and u_j > lambda }
2122  */
2123  ncandsetdelta = 0;
2124  startidx = 4;
2125 
2126  /* get max{ u_j : j in C1 and u_j > lambda }, max{ u_j : j in L~2 and u_j > lambda }, and
2127  * max{ u_j : j in N and u_j >= lambda } and store { u_j : j in N and u_j > lambda }
2128  */
2129  c1vubcoefsmax = SCIP_REAL_MIN;
2130  l2tmpvubcoefsmax = SCIP_REAL_MIN;
2131  nvubcoefsmax = SCIP_REAL_MIN;
2132  for( j = 0; j < ntransvars; j++ )
2133  {
2134  SCIP_Real val;
2135 
2136  /* j is in C1 and u_j > lambda */
2137  if( transvarcoefs[j] == 1 && transvarflowcoverstatus[j] == 1
2138  && SCIPisFeasGT(scip, transvarvubcoefs[j], lambda) )
2139  c1vubcoefsmax = MAX(c1vubcoefsmax, transvarvubcoefs[j]);
2140 
2141  /* j is in L~2 and u_j > lambda */
2142  val = MIN(transvarvubcoefs[j], lambda) * transbinvarsolvals[j];
2143  if( transvarcoefs[j] == -1 && transvarflowcoverstatus[j] == -1 && SCIPisFeasLE(scip, val, transcontvarsolvals[j])
2144  && SCIPisFeasGT(scip, transvarvubcoefs[j], lambda) )
2145  l2tmpvubcoefsmax = MAX(l2tmpvubcoefsmax, transvarvubcoefs[j]);
2146 
2147  /* u_j + 1 > lambda */
2148  if( SCIPisFeasGT(scip, transvarvubcoefs[j] + 1.0, lambda) )
2149  nvubcoefsmax = MAX(nvubcoefsmax, transvarvubcoefs[j]);
2150 
2151  /* u_j > lambda */
2152  if( SCIPisFeasGT(scip, transvarvubcoefs[j], lambda) )
2153  {
2154  candsetdelta[startidx + ncandsetdelta] = transvarvubcoefs[j];
2155  ncandsetdelta++;
2156  SCIPdebugMessage(" u_j = %g\n", transvarvubcoefs[j]);
2157  }
2158  }
2159 
2160  /* store max { u_j : j in N and u_j >= lambda } + 1 */
2161  if( !SCIPisInfinity(scip, -nvubcoefsmax) )
2162  {
2163  assert(SCIPisFeasGT(scip, nvubcoefsmax + 1.0, lambda));
2164  startidx--;
2165  candsetdelta[startidx] = nvubcoefsmax + 1.0;
2166  ncandsetdelta++;
2167  SCIPdebugMessage(" max u_j+1 = %g\n", nvubcoefsmax + 1.0);
2168  }
2169 
2170  /* store lambda + 1 */
2171  startidx--;
2172  candsetdelta[startidx] = lambda + 1.0;
2173  ncandsetdelta++;
2174  SCIPdebugMessage(" lambda + 1 = %g\n", lambda + 1.0);
2175 
2176  /* store max{ u_j : j in C1 and u_j > lambda } and max{ u_j : j in C1 & L~2 and u_j > lambda } */
2177  if( !SCIPisInfinity(scip, -l2tmpvubcoefsmax) && SCIPisFeasGT(scip, l2tmpvubcoefsmax, c1vubcoefsmax) )
2178  {
2179  assert(SCIPisFeasGT(scip, l2tmpvubcoefsmax, lambda));
2180  startidx--;
2181  candsetdelta[startidx] = l2tmpvubcoefsmax;
2182  ncandsetdelta++;
2183  SCIPdebugMessage(" l2tmpvubcoefsmax = %g\n", l2tmpvubcoefsmax);
2184  }
2185  if( !SCIPisInfinity(scip, -c1vubcoefsmax) )
2186  {
2187  assert(SCIPisFeasGT(scip, c1vubcoefsmax, lambda));
2188  startidx--;
2189  candsetdelta[startidx] = c1vubcoefsmax;
2190  ncandsetdelta++;
2191  SCIPdebugMessage(" c1vubcoefsmax = %g\n", c1vubcoefsmax);
2192  }
2193 
2194  assert(startidx >= 0 && startidx <= 3);
2195  assert(startidx + ncandsetdelta >= 4);
2196  assert(ncandsetdelta >= 1 && ncandsetdelta <= ntransvars + 4);
2197 
2198  SCIPdebugMessage("2. generate c-MIRFCIs for different values of delta:\n");
2199 
2200  /* for each value of delta choose L1 subset N1\C1 and L2 subset N2\C2 by comparison, generate the
2201  * c-MIRFCI for delta, (C1, C2) and (L1, L2) and select the most efficient c-MIRFCI
2202  */
2203  ntesteddeltas = 0;
2204  bestdelta = 0.0;
2205  bestefficacy = 0.0;
2206  for( j = startidx; j < startidx + ncandsetdelta && ntesteddeltas < sepadata->maxtestdelta; j++ )
2207  {
2208  SCIP_Real delta;
2209  SCIP_Real onedivdelta;
2210  SCIP_Bool tested;
2211  int i;
2212 
2213  /* get current delta and corresponding scalar for c-MIR routine */
2214  delta = candsetdelta[j];
2215  onedivdelta = 1.0 / delta;
2216 
2217  /* do not use scaling factors (1/delta) which are too small, since this max cause numerical problems;
2218  * besides for relatively small lambda and large scaling factors (delta), we get f_beta = 1 - lambda/delta > MINFRAC
2219  */
2220  if( SCIPisZero(scip, scalar * onedivdelta) )
2221  continue;
2222 
2223  /* check, if delta was already tested */
2224  tested = FALSE;
2225  for( i = 0; i < ntesteddeltas && !tested; i++ )
2226  tested = SCIPisEQ(scip, testeddeltas[i], delta);
2227  if( tested )
2228  continue;
2229  testeddeltas[ntesteddeltas] = delta;
2230  ntesteddeltas++;
2231 
2232  SCIPdebugMessage(" delta = %g:\n", delta);
2233 
2234  /* work on copy of transvarflowcoverstatus because current choice of sets L1 and L2 will change
2235  * transvarflowcoverstatus
2236  */
2237  BMScopyMemoryArray(transvarflowcoverstatustmp, transvarflowcoverstatus, ntransvars);
2238 
2239  /* get L1 subset of N1\C1 and L2 subset of N2\C2 by comparison */
2240  getL1L2(scip, ntransvars, transvarcoefs, transbinvarsolvals, transcontvarsolvals, transvarvubcoefs,
2241  transvarflowcoverstatustmp, delta, lambda);
2242 
2243  /* get bounds for substitution in c-MIR routine for original mixed integer set;
2244  * note that the scalar has already been considered in the constructed 0-1 single node flow relaxation
2245  */
2246  SCIP_CALL( getBoundsForSubstitution(scip, vars, nvars, boundsfortrans, boundtypesfortrans, assoctransvars,
2247  transvarcoefs, transvarflowcoverstatustmp, ntransvars, boundsforsubst, boundtypesforsubst ) );
2248 
2249  /* generate c-MIRFCI for flow cover (C1,C2), L1 subset N1\C1 and L2 subset N2\C2 and delta */
2250  SCIP_CALL( SCIPcalcMIR(scip, sol, BOUNDSWITCH, TRUE, ALLOWLOCAL, FIXINTEGRALRHS, boundsforsubst, boundtypesforsubst,
2251  (int) MAXAGGRLEN(nvars), 1.0, MINFRAC, MAXFRAC, rowweights, NULL, scalar * onedivdelta, NULL, NULL, cutcoefs,
2252  &cutrhs, &cutact, &success, &cutislocal, NULL) );
2253  assert(ALLOWLOCAL || !cutislocal);
2254 
2255  /* delta leads to c-MIRFCI which is more violated */
2256  if( success )
2257  {
2258  SCIP_Real efficacy;
2259 
2260  for(i = 0; i < nvars; i++)
2261  cutcoefs[i] = lambda * cutcoefs[i];
2262  cutrhs = lambda * cutrhs;
2263  cutact = lambda * cutact;
2264 
2265  efficacy = calcEfficacy(nvars, cutcoefs, cutrhs, cutact);
2266 
2267  SCIPdebugMessage(" ---> act = %g rhs = %g eff = %g (old besteff = %g, old bestdelta=%g)\n",
2268  cutact, cutrhs, efficacy, bestefficacy, bestdelta);
2269 
2270  if( efficacy > bestefficacy )
2271  {
2272  bestdelta = delta;
2273  bestefficacy = efficacy;
2274  }
2275  }
2276  }
2277 
2278  /* delta found */
2279  if( SCIPisEfficacious(scip, bestefficacy) )
2280  {
2281  SCIP_Real onedivbestdelta;
2282  int i;
2283 
2284  assert(bestdelta != 0.0);
2285  onedivbestdelta = 1.0 / bestdelta;
2286 
2287  /* for best value of delta: get L1 subset of N1\C1 and L2 subset of N2\C2 by comparison */
2288  getL1L2(scip, ntransvars, transvarcoefs, transbinvarsolvals, transcontvarsolvals, transvarvubcoefs,
2289  transvarflowcoverstatus, bestdelta, lambda);
2290 
2291  /* for best value of delta: get bounds for substitution in c-MIR routine for original mixed integer set
2292  * note that the scalar has already been considered in the constructed 0-1 single node flow relaxation
2293  */
2294  SCIP_CALL( getBoundsForSubstitution(scip, vars, nvars, boundsfortrans, boundtypesfortrans, assoctransvars,
2295  transvarcoefs, transvarflowcoverstatus, ntransvars, boundsforsubst, boundtypesforsubst ) );
2296 
2297  /* generate c-MIRFCI for flow cover (C1,C2), L1 subset N1\C1 and L2 subset N2\C2 and bestdelta */
2298  SCIP_CALL( SCIPcalcMIR(scip, sol, BOUNDSWITCH, TRUE, ALLOWLOCAL, FIXINTEGRALRHS, boundsforsubst, boundtypesforsubst,
2299  (int) MAXAGGRLEN(nvars), 1.0, MINFRAC, MAXFRAC, rowweights, NULL, scalar * onedivbestdelta, NULL, NULL, cutcoefs,
2300  &cutrhs, &cutact, &success, &cutislocal, &cutrank) );
2301  assert(ALLOWLOCAL || !cutislocal);
2302  assert(success);
2303 
2304  for(i = 0; i < nvars; i++)
2305  cutcoefs[i] = lambda * cutcoefs[i];
2306  cutrhs = lambda * cutrhs;
2307  cutact = lambda * cutact;
2308 
2309  assert(SCIPisFeasEQ(scip, bestefficacy, calcEfficacy(nvars, cutcoefs, cutrhs, cutact)));
2310  SCIP_CALL( addCut(scip, sepa, sepadata, vars, nvars, sol, varsolvals, cutcoefs, cutrhs, cutislocal, cutrank, normtype, cutoff, ncuts) );
2311  }
2312 
2313  /* free data structures */
2314  SCIPfreeBufferArray(scip, &cutcoefs);
2315  SCIPfreeBufferArray(scip, &boundtypesforsubst);
2316  SCIPfreeBufferArray(scip, &boundsforsubst);
2317  SCIPfreeBufferArray(scip, &transvarflowcoverstatustmp);
2318  SCIPfreeBufferArray(scip, &candsetdelta);
2319  SCIPfreeBufferArray(scip, &testeddeltas);
2320 
2321  return SCIP_OKAY;
2322 }
2323 
2324 /** comparison method for scores of lp rows */
2325 static
2326 SCIP_DECL_SORTINDCOMP(rowScoreComp)
2327 { /*lint --e{715}*/
2328  SCIP_Real* rowscores = (SCIP_Real*)dataptr;
2329  SCIP_Real tmp;
2330 
2331  assert(rowscores != NULL);
2332  assert(0 <= ind1 && 0 <= ind2);
2333 
2334  tmp = rowscores[ind1] - rowscores[ind2];
2335 
2336  if( tmp < 0 )
2337  return 1;
2338  else if( tmp > 0 )
2339  return -1;
2340  else
2341  return 0;
2342 }
2343 
2344 /** search and add flowcover cuts that separate the given primal solution */
2345 static
2346 SCIP_RETCODE separateCuts(
2347  SCIP* scip, /**< SCIP data structure */
2348  SCIP_SEPA* sepa, /**< the flowcover separator */
2349  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2350  SCIP_RESULT* result /**< pointer to store the result */
2351  )
2352 {
2353  SCIP_ROW** rows;
2354  SCIP_VAR** vars;
2355  SCIP_BOUNDTYPE* boundtypesfortrans;
2356  SCIP_SEPADATA* sepadata;
2357  SCIP_Real* rowlhsscores;
2358  SCIP_Real* rowrhsscores;
2359  SCIP_Real* rowscores;
2360  SCIP_Real* rowweights;
2361  SCIP_Real* transbinvarsolvals;
2362  SCIP_Real* transcontvarsolvals;
2363  SCIP_Real* transvarvubcoefs;
2364  SCIP_Real* varsolvals;
2365  int* assoctransvars;
2366  int* boundsfortrans;
2367  int* covervars;
2368  int* noncovervars;
2369  int* roworder;
2370  int* transvarcoefs;
2371  int* transvarflowcoverstatus;
2372  SCIP_Real lambda;
2373  SCIP_Real maxslack;
2374  SCIP_Real objnorm;
2375  SCIP_Real transcapacity;
2376  SCIP_Bool transsuccess;
2377  SCIP_Bool flowcoverfound;
2378  SCIP_Bool cutoff = FALSE;
2379  char normtype;
2380  int depth;
2381  int maxfails;
2382  int maxsepacuts;
2383  int maxtries;
2384  int ncalls;
2385  int ncovervars;
2386  int nfails;
2387  int nnoncovervars;
2388  int nrows;
2389  int ntransvars;
2390  int ntries;
2391  int nvars;
2392  int ncuts;
2393  int r;
2394 
2395  assert(scip != NULL);
2396  assert(sepa != NULL);
2397  assert(result != NULL);
2398  assert(*result == SCIP_DIDNOTRUN);
2399 
2400  sepadata = SCIPsepaGetData(sepa);
2401  assert(sepadata != NULL);
2402 
2403  depth = SCIPgetDepth(scip);
2404  ncalls = SCIPsepaGetNCallsAtNode(sepa);
2405 
2406  /* only call the flow cover cuts separator a given number of times at each node */
2407  if( (depth == 0 && sepadata->maxroundsroot >= 0 && ncalls >= sepadata->maxroundsroot)
2408  || (depth > 0 && sepadata->maxrounds >= 0 && ncalls >= sepadata->maxrounds) )
2409  return SCIP_OKAY;
2410 
2411  /* get all rows and number of rows */
2412  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2413  assert(nrows == 0 || rows != NULL);
2414 
2415  /* nothing to do, if LP is empty */
2416  if( nrows == 0 )
2417  return SCIP_OKAY;
2418 
2419  /* get active problem variables */
2420  vars = SCIPgetVars(scip);
2421  nvars = SCIPgetNVars(scip);
2422  assert(nvars == 0 || vars != NULL);
2423 
2424  /* nothing to do, if problem has no variables */
2425  if( nvars == 0 )
2426  return SCIP_OKAY;
2427 
2428  /* check whether SCIP was stopped in the meantime */
2429  if( SCIPisStopped(scip) )
2430  return SCIP_OKAY;
2431 
2432  *result = SCIP_DIDNOTFIND;
2433 
2434  /* get the type of norm to use for efficacy calculations */
2435  SCIP_CALL( SCIPgetCharParam(scip, "separating/efficacynorm", &normtype) );
2436 
2437  /* get data structures */
2438  SCIP_CALL( SCIPallocBufferArray(scip, &rowlhsscores, nrows) );
2439  SCIP_CALL( SCIPallocBufferArray(scip, &rowrhsscores, nrows) );
2440  SCIP_CALL( SCIPallocBufferArray(scip, &rowscores, nrows) );
2441  SCIP_CALL( SCIPallocBufferArray(scip, &roworder, nrows) );
2442  SCIP_CALL( SCIPallocBufferArray(scip, &varsolvals, nvars) );
2443  SCIP_CALL( SCIPallocBufferArray(scip, &assoctransvars, nvars) );
2444  SCIP_CALL( SCIPallocBufferArray(scip, &transvarcoefs, nvars) );
2445  SCIP_CALL( SCIPallocBufferArray(scip, &transvarvubcoefs, nvars) );
2446  SCIP_CALL( SCIPallocBufferArray(scip, &transbinvarsolvals, nvars) );
2447  SCIP_CALL( SCIPallocBufferArray(scip, &transcontvarsolvals, nvars) );
2448  SCIP_CALL( SCIPallocBufferArray(scip, &transvarflowcoverstatus, nvars) );
2449  SCIP_CALL( SCIPallocBufferArray(scip, &boundsfortrans, nvars) );
2450  SCIP_CALL( SCIPallocBufferArray(scip, &boundtypesfortrans, nvars) );
2451  SCIP_CALL( SCIPallocBufferArray(scip, &covervars, nvars) );
2452  SCIP_CALL( SCIPallocBufferArray(scip, &noncovervars, nvars) );
2453  SCIP_CALL( SCIPallocBufferArray(scip, &rowweights, nrows) );
2454 
2455  /* get the solution values for all active variables */
2456  SCIP_CALL( SCIPgetSolVals(scip, sol, nvars, vars, varsolvals) );
2457 
2458  /* get the maximal number of cuts allowed in a separation round */
2459  if( depth == 0 )
2460  {
2461  maxtries = sepadata->maxtriesroot;
2462  maxfails = sepadata->maxfailsroot;
2463  maxsepacuts = sepadata->maxsepacutsroot;
2464  maxslack = sepadata->maxslackroot;
2465  }
2466  else
2467  {
2468  maxtries = sepadata->maxtries;
2469  maxfails = sepadata->maxfails;
2470  maxsepacuts = sepadata->maxsepacuts;
2471  maxslack = sepadata->maxslack;
2472  }
2473 
2474  /* calculate row scores for both sides of all rows, and sort rows by nonincreasing maximal score */
2475  objnorm = SCIPgetObjNorm(scip);
2476  objnorm = MAX(objnorm, 1.0);
2477  for( r = 0; r < nrows; r++ )
2478  {
2479  int nnonz;
2480 
2481  assert(SCIProwGetLPPos(rows[r]) == r);
2482 
2483  nnonz = SCIProwGetNLPNonz(rows[r]);
2484  if( nnonz == 0 )
2485  {
2486  /* ignore empty rows */
2487  rowlhsscores[r] = 0.0;
2488  rowrhsscores[r] = 0.0;
2489  }
2490  else
2491  {
2492  SCIP_Real activity;
2493  SCIP_Real lhs;
2494  SCIP_Real rhs;
2495  SCIP_Real dualsol;
2496  SCIP_Real dualscore;
2497  SCIP_Real rowdensity;
2498  SCIP_Real rownorm;
2499  SCIP_Real slack;
2500 
2501  dualsol = (sol == NULL ? SCIProwGetDualsol(rows[r]) : 1.0);
2502  activity = SCIPgetRowSolActivity(scip, rows[r], sol);
2503  lhs = SCIProwGetLhs(rows[r]);
2504  rhs = SCIProwGetRhs(rows[r]);
2505  rownorm = SCIProwGetNorm(rows[r]);
2506  rownorm = MAX(rownorm, 0.1);
2507  rowdensity = (SCIP_Real)nnonz/(SCIP_Real)nvars;
2508  assert(SCIPisPositive(scip, rownorm));
2509 
2510  slack = (activity - lhs)/rownorm;
2511  dualscore = MAX(dualsol/objnorm, 0.0001);
2512  if( !SCIPisInfinity(scip, -lhs) && SCIPisLE(scip, slack, maxslack)
2513  && (ALLOWLOCAL || !SCIProwIsLocal(rows[r])) /*lint !e506 !e774*/
2514  && rowdensity <= sepadata->maxrowdensity )
2515  {
2516  rowlhsscores[r] = dualscore + DENSSCORE * rowdensity + sepadata->slackscore * MAX(1.0 - slack, 0.0);
2517  assert(rowlhsscores[r] > 0.0);
2518  }
2519  else
2520  rowlhsscores[r] = 0.0;
2521 
2522  slack = (rhs - activity)/rownorm;
2523  dualscore = MAX(-dualsol/objnorm, 0.0001);
2524  if( !SCIPisInfinity(scip, rhs) && SCIPisLE(scip, slack, maxslack)
2525  && (ALLOWLOCAL || !SCIProwIsLocal(rows[r])) /*lint !e506 !e774*/
2526  && rowdensity <= sepadata->maxrowdensity )
2527  {
2528  rowrhsscores[r] = dualscore + DENSSCORE * rowdensity + sepadata->slackscore * MAX(1.0 - slack, 0.0);
2529  assert(rowrhsscores[r] > 0.0);
2530  }
2531  else
2532  rowrhsscores[r] = 0.0;
2533  }
2534  rowscores[r] = MAX(rowlhsscores[r], rowrhsscores[r]);
2535 
2536  roworder[r] = r;
2537  }
2538 
2539  SCIPsortInd(roworder, rowScoreComp, (void*) rowscores, nrows);
2540 #ifndef NDEBUG
2541  for( r = nrows - 1; r > 0; --r )
2542  assert(rowscores[roworder[r]] <= rowscores[roworder[r - 1]]);
2543 #endif
2544 
2545  /* initialize weights of rows for aggregation for c-mir routine */
2546  BMSclearMemoryArray(rowweights, nrows);
2547 
2548  /* try to generate a flow cover cut for each row in the LP */
2549  ncuts = 0;
2550  if( maxtries < 0 )
2551  maxtries = INT_MAX;
2552  if( maxfails < 0 )
2553  maxfails = INT_MAX;
2554  else if( depth == 0 && 2*SCIPgetNSepaRounds(scip) < maxfails )
2555  maxfails += maxfails - 2*SCIPgetNSepaRounds(scip); /* allow up to double as many fails in early separounds of root node */
2556  ntries = 0;
2557  nfails = 0;
2558  for( r = 0; r < nrows && ntries < maxtries && ncuts < maxsepacuts && rowscores[roworder[r]] > 0.0
2559  && !SCIPisStopped(scip); r++ )
2560  {
2561  SCIP_Bool wastried;
2562  int oldncuts;
2563  SCIP_Real rowact;
2564  SCIP_Real mult;
2565 
2566  oldncuts = ncuts;
2567  wastried = FALSE;
2568 
2569  /* update weight of rows for aggregation in c-MIR routine; all rows but current one have weight 0.0 */
2570  if( r > 0 )
2571  {
2572  assert(rowweights[roworder[r-1]] == 1.0 || rowweights[roworder[r-1]] == -1.0 || rowweights[roworder[r-1]] == 0.0);
2573  rowweights[roworder[r-1]] = 0.0;
2574  }
2575 
2576  /* decide which side of the row should be used */
2577  rowact = SCIPgetRowSolActivity(scip, rows[roworder[r]], sol);
2578  if( !SCIPisInfinity(scip, -SCIProwGetLhs(rows[roworder[r]]))
2579  && rowact < 0.5 * SCIProwGetLhs(rows[roworder[r]]) + 0.5 * SCIProwGetRhs(rows[roworder[r]]) )
2580  rowweights[roworder[r]] = -1.0;
2581  else if( !SCIPisInfinity(scip, SCIProwGetRhs(rows[roworder[r]])) )
2582  rowweights[roworder[r]] = 1.0;
2583  else
2584  {
2585  SCIPdebugMessage("skipping unconstrained row\n");
2586  rowweights[roworder[r]] = 0.0;
2587  continue;
2588  }
2589 
2590  SCIPdebugMessage("===================== flow cover separation for row <%s> (%d of %d) ===================== \n",
2591  SCIProwGetName(rows[roworder[r]]), r, nrows);
2592  SCIPdebug( SCIP_CALL( SCIPprintRow(scip, rows[roworder[r]], NULL) ) );
2593  SCIPdebugMessage("rowact=%g is closer to %s --> rowweight=%g\n", rowact,
2594  rowweights[roworder[r]] == 1 ? "rhs" : "lhs", rowweights[roworder[r]]);
2595 
2596  mult = 1.0;
2597  do
2598  {
2599  assert(mult == 1.0 || (mult == -1.0 && sepadata->multbyminusone));
2600 
2601  /* construct 0-1 single node flow relaxation (with some additional simple constraints) of the mixed integer set
2602  * corresponding to the current row
2603  * sum_{j in N} a_j x_j + sum_{j in M} c_j y_j + s = rhs + const or
2604  * - ( sum_{j in N} a_j x_j + sum_{j in M} c_j y_j ) + s = - ( lhs - const )
2605  * multiplied by mult in { +1, -1 }
2606  */
2607  SCIP_CALL( constructSNFRelaxation(scip, vars, nvars, varsolvals, rows[roworder[r]], rowweights[roworder[r]],
2608  mult, boundsfortrans, boundtypesfortrans, assoctransvars, transvarcoefs, transbinvarsolvals,
2609  transcontvarsolvals, transvarvubcoefs, &ntransvars, &transcapacity, &transsuccess) );
2610  if( !transsuccess )
2611  {
2612  /* transformation will fail for mult = -1, too */
2613  SCIPdebugMessage("mult=%g: no 0-1 single node flow relaxation found\n", mult);
2614  break;
2615  }
2616 
2617  flowcoverfound = FALSE;
2618 
2619  /* get a flow cover (C1, C2) for the constructed 0-1 single node flow set */
2620  SCIP_CALL( getFlowCover(scip, transvarcoefs, transbinvarsolvals, transvarvubcoefs, ntransvars, transcapacity,
2621  &ncovervars, &nnoncovervars, transvarflowcoverstatus, &lambda, &flowcoverfound) );
2622  if( !flowcoverfound )
2623  {
2624  SCIPdebugMessage("mult=%g: no flow cover found\n", mult);
2625  mult *= -1.0;
2626  continue;
2627  }
2628  assert(SCIPisFeasGT(scip, lambda, 0.0));
2629 
2630  /* generate most violated c-MIRFCI for different sets L1 and L2 and different values of delta and add it to the LP */
2631  SCIP_CALL( cutGenerationHeuristic(scip, sepa, sepadata, vars, nvars, sol, varsolvals, rowweights, mult, boundsfortrans,
2632  boundtypesfortrans, assoctransvars, ntransvars, transvarcoefs, transbinvarsolvals, transcontvarsolvals,
2633  transvarvubcoefs, transvarflowcoverstatus, lambda, normtype, &cutoff, &ncuts) );
2634 
2635  wastried = TRUE;
2636  mult *= -1.0;
2637 
2638  if ( cutoff )
2639  break;
2640  }
2641  while( sepadata->multbyminusone && mult < 0.0 );
2642 
2643  if ( cutoff )
2644  break;
2645 
2646  if( !wastried )
2647  continue;
2648 
2649  ntries++;
2650  if( ncuts == oldncuts )
2651  {
2652  nfails++;
2653  if( nfails >= maxfails )
2654  break;
2655  }
2656  else
2657  nfails = 0;
2658  }
2659 
2660  /* free data structures */
2661  SCIPfreeBufferArray(scip, &rowweights);
2662  SCIPfreeBufferArray(scip, &noncovervars);
2663  SCIPfreeBufferArray(scip, &covervars);
2664  SCIPfreeBufferArray(scip, &boundtypesfortrans);
2665  SCIPfreeBufferArray(scip, &boundsfortrans);
2666  SCIPfreeBufferArray(scip, &transvarflowcoverstatus);
2667  SCIPfreeBufferArray(scip, &transcontvarsolvals);
2668  SCIPfreeBufferArray(scip, &transbinvarsolvals);
2669  SCIPfreeBufferArray(scip, &transvarvubcoefs);
2670  SCIPfreeBufferArray(scip, &transvarcoefs);
2671  SCIPfreeBufferArray(scip, &assoctransvars);
2672  SCIPfreeBufferArray(scip, &varsolvals);
2673  SCIPfreeBufferArray(scip, &roworder);
2674  SCIPfreeBufferArray(scip, &rowscores);
2675  SCIPfreeBufferArray(scip, &rowrhsscores);
2676  SCIPfreeBufferArray(scip, &rowlhsscores);
2677 
2678  if ( cutoff )
2679  *result = SCIP_CUTOFF;
2680  else if ( ncuts > 0 )
2681  *result = SCIP_SEPARATED;
2682 
2683  return SCIP_OKAY;
2684 }
2685 
2686 
2687 /*
2688  * Callback methods of separator
2689  */
2690 
2691 /** copy method for separator plugins (called when SCIP copies plugins) */
2692 static
2693 SCIP_DECL_SEPACOPY(sepaCopyFlowcover)
2694 { /*lint --e{715}*/
2695  assert(scip != NULL);
2696  assert(sepa != NULL);
2697  assert(strcmp(SCIPsepaGetName(sepa), SEPA_NAME) == 0);
2698 
2699  /* call inclusion method of constraint handler */
2700  SCIP_CALL( SCIPincludeSepaFlowcover(scip) );
2701 
2702  return SCIP_OKAY;
2703 }
2704 
2705 /** destructor of separator to free user data (called when SCIP is exiting) */
2706 static
2707 SCIP_DECL_SEPAFREE(sepaFreeFlowcover)
2708 { /*lint --e{715}*/
2709  SCIP_SEPADATA* sepadata;
2710 
2711  /* free separator data */
2712  sepadata = SCIPsepaGetData(sepa);
2713  assert(sepadata != NULL);
2714 
2715  SCIPfreeMemory(scip, &sepadata);
2716 
2717  SCIPsepaSetData(sepa, NULL);
2718 
2719  return SCIP_OKAY;
2720 }
2721 
2722 
2723 /** LP solution separation method of separator */
2724 static
2725 SCIP_DECL_SEPAEXECLP(sepaExeclpFlowcover)
2726 { /*lint --e{715}*/
2727 
2728  *result = SCIP_DIDNOTRUN;
2729 
2730  /* only call separator, if we are not close to terminating */
2731  if( SCIPisStopped(scip) )
2732  return SCIP_OKAY;
2733 
2734  /* only call separator, if an optimal LP solution is at hand */
2735  if( SCIPgetLPSolstat(scip) != SCIP_LPSOLSTAT_OPTIMAL )
2736  return SCIP_OKAY;
2737 
2738  /* only call separator, if there are fractional variables */
2739  if( SCIPgetNLPBranchCands(scip) == 0 )
2740  return SCIP_OKAY;
2741 
2742  SCIP_CALL( separateCuts(scip, sepa, NULL, result) );
2743 
2744  return SCIP_OKAY;
2745 }
2746 
2747 
2748 /** arbitrary primal solution separation method of separator */
2749 static
2750 SCIP_DECL_SEPAEXECSOL(sepaExecsolFlowcover)
2751 { /*lint --e{715}*/
2752 
2753  *result = SCIP_DIDNOTRUN;
2754 
2755  SCIP_CALL( separateCuts(scip, sepa, sol, result) );
2756 
2757  return SCIP_OKAY;
2758 }
2759 
2760 
2761 /*
2762  * separator specific interface methods
2763  */
2764 
2765 /** creates the flowcover separator and includes it in SCIP */
2766 SCIP_RETCODE SCIPincludeSepaFlowcover(
2767  SCIP* scip /**< SCIP data structure */
2768  )
2769 {
2770  SCIP_SEPADATA* sepadata;
2771  SCIP_SEPA* sepa;
2772 
2773  /* create flowcover separator data */
2774  SCIP_CALL( SCIPallocMemory(scip, &sepadata) );
2775 
2776  /* include separator */
2777  SCIP_CALL( SCIPincludeSepaBasic(scip, &sepa, SEPA_NAME, SEPA_DESC, SEPA_PRIORITY, SEPA_FREQ, SEPA_MAXBOUNDDIST,
2778  SEPA_USESSUBSCIP, SEPA_DELAY,
2779  sepaExeclpFlowcover, sepaExecsolFlowcover,
2780  sepadata) );
2781 
2782  assert(sepa != NULL);
2783 
2784  /* set non-NULL pointers to callback methods */
2785  SCIP_CALL( SCIPsetSepaCopy(scip, sepa, sepaCopyFlowcover) );
2786  SCIP_CALL( SCIPsetSepaFree(scip, sepa, sepaFreeFlowcover) );
2787 
2788  /* add flow cover cuts separator parameters */
2789  SCIP_CALL( SCIPaddIntParam(scip,
2790  "separating/flowcover/maxrounds",
2791  "maximal number of separation rounds per node (-1: unlimited)",
2792  &sepadata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
2793  SCIP_CALL( SCIPaddIntParam(scip,
2794  "separating/flowcover/maxroundsroot",
2795  "maximal number of separation rounds in the root node (-1: unlimited)",
2796  &sepadata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
2797  SCIP_CALL( SCIPaddIntParam(scip,
2798  "separating/flowcover/maxtries",
2799  "maximal number of rows to separate flow cover cuts for per separation round (-1: unlimited)",
2800  &sepadata->maxtries, TRUE, DEFAULT_MAXTRIES, -1, INT_MAX, NULL, NULL) );
2801  SCIP_CALL( SCIPaddIntParam(scip,
2802  "separating/flowcover/maxtriesroot",
2803  "maximal number of rows to separate flow cover cuts for per separation round in the root (-1: unlimited)",
2804  &sepadata->maxtriesroot, TRUE, DEFAULT_MAXTRIESROOT, -1, INT_MAX, NULL, NULL) );
2805  SCIP_CALL( SCIPaddIntParam(scip,
2806  "separating/flowcover/maxfails",
2807  "maximal number of consecutive fails to generate a cut per separation round (-1: unlimited)",
2808  &sepadata->maxfails, TRUE, DEFAULT_MAXFAILS, -1, INT_MAX, NULL, NULL) );
2809  SCIP_CALL( SCIPaddIntParam(scip,
2810  "separating/flowcover/maxfailsroot",
2811  "maximal number of consecutive fails to generate a cut per separation round in the root (-1: unlimited)",
2812  &sepadata->maxfailsroot, TRUE, DEFAULT_MAXFAILSROOT, -1, INT_MAX, NULL, NULL) );
2813  SCIP_CALL( SCIPaddIntParam(scip,
2814  "separating/flowcover/maxsepacuts",
2815  "maximal number of flow cover cuts separated per separation round",
2816  &sepadata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
2817  SCIP_CALL( SCIPaddIntParam(scip,
2818  "separating/flowcover/maxsepacutsroot",
2819  "maximal number of flow cover cuts separated per separation round in the root",
2820  &sepadata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
2821  SCIP_CALL( SCIPaddRealParam(scip,
2822  "separating/flowcover/maxslack",
2823  "maximal slack of rows to separate flow cover cuts for",
2824  &sepadata->maxslack, TRUE, DEFAULT_MAXSLACK, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2825  SCIP_CALL( SCIPaddRealParam(scip,
2826  "separating/flowcover/maxslackroot",
2827  "maximal slack of rows to separate flow cover cuts for in the root",
2828  &sepadata->maxslackroot, TRUE, DEFAULT_MAXSLACKROOT, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2829  SCIP_CALL( SCIPaddRealParam(scip,
2830  "separating/flowcover/slackscore",
2831  "weight of slack in the scoring of the rows",
2832  &sepadata->slackscore, TRUE, DEFAULT_SLACKSCORE, 0.0, SCIP_REAL_MAX, NULL, NULL) );
2833  SCIP_CALL( SCIPaddRealParam(scip,
2834  "separating/flowcover/maxrowdensity",
2835  "maximal density of row to separate flow cover cuts for",
2836  &sepadata->maxrowdensity, TRUE, DEFAULT_MAXROWDENSITY, 0.0, 1.0, NULL, NULL) );
2837  SCIP_CALL( SCIPaddBoolParam(scip,
2838  "separating/flowcover/dynamiccuts",
2839  "should generated cuts be removed from the LP if they are no longer tight?",
2840  &sepadata->dynamiccuts, FALSE, DEFAULT_DYNAMICCUTS, NULL, NULL) );
2841  SCIP_CALL( SCIPaddBoolParam(scip,
2842  "separating/flowcover/multbyminusone",
2843  "should flow cover cuts be separated for 0-1 single node flow set with reversed arcs in addition?",
2844  &sepadata->multbyminusone, TRUE, DEFAULT_MULTBYMINUSONE, NULL, NULL) );
2845  SCIP_CALL( SCIPaddIntParam(scip,
2846  "separating/flowcover/maxtestdelta",
2847  "cut generation heuristic: maximal number of different deltas to try",
2848  &sepadata->maxtestdelta, TRUE, DEFAULT_MAXTESTDELTA, 0, INT_MAX, NULL, NULL) );
2849  return SCIP_OKAY;
2850 }
2851