Scippy

SCIP

Solving Constraint Integer Programs

cuts.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-2020 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 cuts.c
17  * @ingroup OTHER_CFILES
18  * @brief methods for aggregation of rows
19  * @author Jakob Witzig
20  * @author Leona Gottwald
21  */
22 
23 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
24 
25 #include "blockmemshell/memory.h"
26 #include "scip/cuts.h"
27 #include "scip/dbldblarith.h"
28 #include "scip/lp.h"
29 #include "scip/pub_lp.h"
30 #include "scip/pub_message.h"
31 #include "scip/pub_misc.h"
32 #include "scip/pub_misc_select.h"
33 #include "scip/pub_misc_sort.h"
34 #include "scip/pub_var.h"
35 #include "scip/scip_cut.h"
36 #include "scip/scip_lp.h"
37 #include "scip/scip_mem.h"
38 #include "scip/scip_message.h"
39 #include "scip/scip_numerics.h"
40 #include "scip/scip_prob.h"
41 #include "scip/scip_sol.h"
42 #include "scip/scip_solvingstats.h"
43 #include "scip/scip_var.h"
44 #include "scip/struct_lp.h"
45 #include "scip/struct_scip.h"
46 #include "scip/struct_set.h"
47 
48 /* =========================================== general static functions =========================================== */
49 #ifdef SCIP_DEBUG
50 static
51 void printCutQuad(
52  SCIP* scip, /**< SCIP data structure */
53  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
54  SCIP_Real* cutcoefs, /**< non-zero coefficients of cut */
55  QUAD(SCIP_Real cutrhs), /**< right hand side of the MIR row */
56  int* cutinds, /**< indices of problem variables for non-zero coefficients */
57  int cutnnz, /**< number of non-zeros in cut */
58  SCIP_Bool ignorsol,
59  SCIP_Bool islocal
60  )
61 {
62  SCIP_Real QUAD(activity);
63  SCIP_VAR** vars;
64  int i;
65 
66  assert(scip != NULL);
67  vars = SCIPgetVars(scip);
68 
69  SCIPdebugMessage("CUT:");
70  QUAD_ASSIGN(activity, 0.0);
71  for( i = 0; i < cutnnz; ++i )
72  {
73  SCIP_Real QUAD(coef);
74 
75  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[i]);
76 
77  SCIPdebugPrintf(" %+g<%s>", QUAD_TO_DBL(coef), SCIPvarGetName(vars[cutinds[i]]));
78 
79  if( !ignorsol )
80  {
81  SCIPquadprecProdQD(coef, coef, (sol == NULL ? SCIPvarGetLPSol(vars[cutinds[i]]) : SCIPgetSolVal(scip, sol, vars[cutinds[i]])));
82  }
83  else
84  {
85  if( cutcoefs[i] > 0.0 )
86  {
87  SCIPquadprecProdQD(coef, coef, (islocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]])));
88  }
89  else
90  {
91  SCIPquadprecProdQD(coef, coef, (islocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]])));
92  }
93  }
94 
95  SCIPquadprecSumQQ(activity, activity, coef);
96  }
97  SCIPdebugPrintf(" <= %.6f (activity: %g)\n", QUAD_TO_DBL(cutrhs), QUAD_TO_DBL(activity));
98 }
99 #endif
100 
101 /** macro to make sure a value is not equal to zero, i.e. NONZERO(x) != 0.0
102  * will be TRUE for every x including 0.0
103  *
104  * To avoid branches it will add 1e-100 with the same sign as x to x which will
105  * be rounded away for any sane non-zero value but will make sure the value is
106  * never exactly 0.0.
107  */
108 #define NONZERO(x) (COPYSIGN(1e-100, (x)) + (x))
109 
110 /** add a scaled row to a dense vector indexed over the problem variables and keep the
111  * index of non-zeros up-to-date
112  */
113 static
115  int*RESTRICT inds, /**< pointer to array with variable problem indices of non-zeros in variable vector */
116  SCIP_Real*RESTRICT vals, /**< array with values of variable vector */
117  int*RESTRICT nnz, /**< number of non-zeros coefficients of variable vector */
118  SCIP_ROW* row, /**< row coefficients to add to variable vector */
119  SCIP_Real scale /**< scale for adding given row to variable vector */
120  )
121 {
122  /* add up coefficients */
123  int i;
124 
125  assert(inds != NULL);
126  assert(vals != NULL);
127  assert(nnz != NULL);
128  assert(row != NULL);
129 
130  /* add the non-zeros to the aggregation row and keep non-zero index up to date */
131  for( i = 0 ; i < row->len; ++i )
132  {
133  SCIP_Real val;
134  int probindex = row->cols[i]->var_probindex;
135 
136  val = vals[probindex];
137 
138  if( val == 0.0 )
139  inds[(*nnz)++] = probindex;
140 
141  val += row->vals[i] * scale;
142 
143  /* the value must not be exactly zero due to sparsity pattern */
144  val = NONZERO(val);
145 
146  assert(val != 0.0);
147  vals[probindex] = val;
148  }
149 
150  return SCIP_OKAY;
151 }
152 
153 /** add a scaled row to a dense vector indexed over the problem variables and keep the
154  * index of non-zeros up-to-date
155  */
156 static
158  int*RESTRICT inds, /**< pointer to array with variable problem indices of non-zeros in variable vector */
159  SCIP_Real*RESTRICT vals, /**< array with values of variable vector */
160  int*RESTRICT nnz, /**< number of non-zeros coefficients of variable vector */
161  SCIP_ROW* row, /**< row coefficients to add to variable vector */
162  SCIP_Real scale /**< scale for adding given row to variable vector */
163  )
164 {
165  /* add up coefficients */
166  int i;
167 
168  assert(inds != NULL);
169  assert(vals != NULL);
170  assert(nnz != NULL);
171  assert(row != NULL);
172 
173  /* add the non-zeros to the aggregation row and keep non-zero index up to date */
174  for( i = 0 ; i < row->len; ++i )
175  {
176  SCIP_Real QUAD(val);
177  int probindex = row->cols[i]->var_probindex;
178 
179  QUAD_ARRAY_LOAD(val, vals, probindex);
180 
181  if( QUAD_HI(val) == 0.0 )
182  inds[(*nnz)++] = probindex;
183 
184  SCIPquadprecSumQD(val, val, row->vals[i] * scale);
185 
186  /* the value must not be exactly zero due to sparsity pattern */
187  QUAD_HI(val) = NONZERO(QUAD_HI(val));
188  assert(QUAD_HI(val) != 0.0);
189 
190  QUAD_ARRAY_STORE(vals, probindex, val);
191  }
192 
193  return SCIP_OKAY;
194 }
195 
196 /** calculates the cuts efficacy for the given solution */
197 static
199  SCIP* scip, /**< SCIP data structure */
200  SCIP_SOL* sol, /**< solution to calculate the efficacy for (NULL for LP solution) */
201  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
202  SCIP_Real cutrhs, /**< the right hand side of the cut */
203  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
204  int cutnnz /**< the number of non-zeros in the cut */
205  )
206 {
207  SCIP_VAR** vars;
208  SCIP_Real norm;
209  SCIP_Real activity;
210  int i;
211 
212  assert(scip != NULL);
213  assert(cutcoefs != NULL);
214  assert(cutinds != NULL);
215 
216  vars = SCIPgetVars(scip);
217 
218  activity = 0.0;
219  for( i = 0; i < cutnnz; ++i )
220  activity += cutcoefs[i] * SCIPgetSolVal(scip, sol, vars[cutinds[i]]);
221 
222  norm = SCIPgetVectorEfficacyNorm(scip, cutcoefs, cutnnz);
223  return (activity - cutrhs) / MAX(1e-6, norm);
224 }
225 
226 /** calculates the efficacy norm of the given aggregation row, which depends on the "separating/efficacynorm" parameter */
227 static
229  SCIP* scip, /**< SCIP data structure */
230  SCIP_Real* vals, /**< array of the non-zero coefficients in the vector; this is a quad precision array! */
231  int* inds, /**< array of the problem indices of variables with a non-zero coefficient in the vector */
232  int nnz /**< the number of non-zeros in the vector */
233  )
234 {
235  SCIP_Real norm = 0.0;
236  SCIP_Real QUAD(coef);
237  int i;
238 
239  assert(scip != NULL);
240  assert(scip->set != NULL);
241 
242  switch( scip->set->sepa_efficacynorm )
243  {
244  case 'e':
245  for( i = 0; i < nnz; ++i )
246  {
247  QUAD_ARRAY_LOAD(coef, vals, inds[i]);
248  norm += SQR(QUAD_TO_DBL(coef));
249  }
250  norm = SQRT(norm);
251  break;
252  case 'm':
253  for( i = 0; i < nnz; ++i )
254  {
255  SCIP_Real absval;
256  QUAD_ARRAY_LOAD(coef, vals, inds[i]);
257 
258  absval = REALABS(QUAD_TO_DBL(coef));
259  norm = MAX(norm, absval);
260  }
261  break;
262  case 's':
263  for( i = 0; i < nnz; ++i )
264  {
265  QUAD_ARRAY_LOAD(coef, vals, inds[i]);
266  norm += REALABS(QUAD_TO_DBL(coef));
267  }
268  break;
269  case 'd':
270  for( i = 0; i < nnz; ++i )
271  {
272  QUAD_ARRAY_LOAD(coef, vals, inds[i]);
273  if( !SCIPisZero(scip, QUAD_TO_DBL(coef)) )
274  {
275  norm = 1.0;
276  break;
277  }
278  }
279  break;
280  default:
281  SCIPerrorMessage("invalid efficacy norm parameter '%c'\n", scip->set->sepa_efficacynorm);
282  assert(FALSE);
283  }
284 
285  return norm;
286 }
287 
288 /** calculates the cuts efficacy for the given solution; the cut coefs are stored densely and in quad precision */
289 static
291  SCIP* scip, /**< SCIP data structure */
292  SCIP_SOL* sol, /**< solution to calculate the efficacy for (NULL for LP solution) */
293  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut; this is a quad precision array! */
294  SCIP_Real cutrhs, /**< the right hand side of the cut */
295  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
296  int cutnnz /**< the number of non-zeros in the cut */
297  )
298 {
299  SCIP_VAR** vars;
300  SCIP_Real norm;
301  SCIP_Real activity;
302  SCIP_Real QUAD(coef);
303  int i;
304 
305  assert(scip != NULL);
306  assert(cutcoefs != NULL);
307  assert(cutinds != NULL);
308  assert(scip->set != NULL);
309 
310  vars = SCIPgetVars(scip);
311 
312  activity = 0.0;
313  for( i = 0; i < cutnnz; ++i )
314  {
315  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[i]);
316  activity += QUAD_TO_DBL(coef) * SCIPgetSolVal(scip, sol, vars[cutinds[i]]);
317  }
318 
319  norm = calcEfficacyNormQuad(scip, cutcoefs, cutinds, cutnnz);
320  return (activity - cutrhs) / MAX(1e-6, norm);
321 }
322 
323 /** safely remove all items with |a_i| or |u_i - l_i)| below the given value; returns TRUE if the cut became redundant
324  * if it is a local cut, use local bounds, otherwise, using global bounds
325  * */
326 static
328  SCIP* scip, /**< SCIP data structure */
329  SCIP_Real minval, /**< minimal absolute value of coefficients that should not be removed */
330  SCIP_Bool cutislocal, /**< is the cut local? */
331  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
332  QUAD(SCIP_Real* cutrhs), /**< the right hand side of the cut */
333  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
334  int* cutnnz /**< the number of non-zeros in the cut */
335  )
336 {
337  int i;
338  SCIP_VAR** vars;
339 
340  vars = SCIPgetVars(scip);
341 
342  for( i = 0; i < *cutnnz; )
343  {
344  SCIP_Real QUAD(val);
345  SCIP_Real lb;
346  SCIP_Real ub;
347  int v;
348  SCIP_Bool isfixed;
349  v = cutinds[i];
350  QUAD_ARRAY_LOAD(val, cutcoefs, v);
351  if( cutislocal )
352  {
353  lb = SCIPvarGetLbLocal(vars[v]);
354  ub = SCIPvarGetUbLocal(vars[v]);
355  }
356  else
357  {
358  lb = SCIPvarGetLbGlobal(vars[v]);
359  ub = SCIPvarGetUbGlobal(vars[v]);
360  }
361  if( !(SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub)) && SCIPisEQ(scip, ub, lb) )
362  isfixed = TRUE;
363  else
364  isfixed = FALSE;
365 
366  if( EPSZ(QUAD_TO_DBL(val), minval) || isfixed )
367  {
368  if( REALABS(QUAD_TO_DBL(val)) > QUAD_EPSILON )
369  {
370  /* adjust left and right hand sides with max contribution */
371  if( QUAD_TO_DBL(val) < 0.0 )
372  {
373  if( SCIPisInfinity(scip, ub) )
374  return TRUE;
375  else
376  {
377  SCIPquadprecProdQD(val, val, ub);
378  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -val);
379  }
380  }
381  else
382  {
383  if( SCIPisInfinity(scip, -lb) )
384  return TRUE;
385  else
386  {
387  SCIPquadprecProdQD(val, val, lb);
388  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -val);
389  }
390  }
391  }
392 
393  QUAD_ASSIGN(val, 0.0);
394  QUAD_ARRAY_STORE(cutcoefs, v, val);
395 
396  /* remove non-zero entry */
397  --(*cutnnz);
398  cutinds[i] = cutinds[*cutnnz];
399  }
400  else
401  ++i;
402  }
403 
404  /* relax rhs to zero, if it's very close to */
405  if( QUAD_TO_DBL(*cutrhs) < 0.0 && QUAD_TO_DBL(*cutrhs) >= -SCIPepsilon(scip) )
406  QUAD_ASSIGN(*cutrhs, 0.0);
407 
408  return FALSE;
409 }
410 
411 /** safely remove all items with |a_i| or |u_i - l_i)| below the given value; returns TRUE if the cut became redundant
412  * if it is a local cut, use local bounds, otherwise, using global bounds
413  * */
414 static
416  SCIP* scip, /**< SCIP data structure */
417  SCIP_Real minval, /**< minimal absolute value of coefficients that should not be removed */
418  SCIP_Bool cutislocal, /**< is the cut local? */
419  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
420  QUAD(SCIP_Real* cutrhs), /**< the right hand side of the cut */
421  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
422  int* cutnnz /**< the number of non-zeros in the cut */
423  )
424 {
425  int i;
426  SCIP_VAR** vars;
427 
428  vars = SCIPgetVars(scip);
429 
430  /* loop over non-zeros and remove values below minval; values above QUAD_EPSILON are cancelled with their bound
431  * to avoid numerical rounding errors
432  */
433  for( i = 0; i < *cutnnz; )
434  {
435  SCIP_Real val;
436  SCIP_Real lb;
437  SCIP_Real ub;
438  int v;
439  SCIP_Bool isfixed;
440  v = cutinds[i];
441  val = cutcoefs[v];
442  if( cutislocal )
443  {
444  lb = SCIPvarGetLbLocal(vars[v]);
445  ub = SCIPvarGetUbLocal(vars[v]);
446  }
447  else
448  {
449  lb = SCIPvarGetLbGlobal(vars[v]);
450  ub = SCIPvarGetUbGlobal(vars[v]);
451  }
452  if( !(SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub)) && SCIPisEQ(scip, ub, lb) )
453  isfixed = TRUE;
454  else
455  isfixed = FALSE;
456 
457  if( EPSZ(val, minval) || isfixed )
458  {
459  if( REALABS(val) > QUAD_EPSILON )
460  {
461  /* adjust left and right hand sides with max contribution */
462  if( val < 0.0 )
463  {
464  if( SCIPisInfinity(scip, ub) )
465  return TRUE;
466  else
467  {
468  SCIPquadprecSumQD(*cutrhs, *cutrhs, -val * ub);
469  }
470  }
471  else
472  {
473  if( SCIPisInfinity(scip, -lb) )
474  return TRUE;
475  else
476  {
477  SCIPquadprecSumQD(*cutrhs, *cutrhs, -val * lb);
478  }
479  }
480  }
481 
482  cutcoefs[v] = 0.0;
483 
484  /* remove non-zero entry */
485  --(*cutnnz);
486  cutinds[i] = cutinds[*cutnnz];
487  }
488  else
489  ++i;
490  }
491 
492  /* relax rhs to zero, if it's very close to */
493  if( QUAD_TO_DBL(*cutrhs) < 0.0 && QUAD_TO_DBL(*cutrhs) >= -SCIPepsilon(scip) )
494  QUAD_ASSIGN(*cutrhs, 0.0);
495 
496  return FALSE;
497 }
498 
499 static
500 SCIP_DECL_SORTINDCOMP(compareAbsCoefsQuad)
501 {
502  SCIP_Real abscoef1;
503  SCIP_Real abscoef2;
504  SCIP_Real QUAD(coef1);
505  SCIP_Real QUAD(coef2);
506  SCIP_Real* coefs = (SCIP_Real*) dataptr;
507 
508  QUAD_ARRAY_LOAD(coef1, coefs, ind1);
509  QUAD_ARRAY_LOAD(coef2, coefs, ind2);
510 
511  abscoef1 = REALABS(QUAD_TO_DBL(coef1));
512  abscoef2 = REALABS(QUAD_TO_DBL(coef2));
513 
514  if( abscoef1 < abscoef2 )
515  return -1;
516  if( abscoef2 < abscoef1 )
517  return 1;
518 
519  return 0;
520 }
521 
522 static
523 SCIP_DECL_SORTINDCOMP(compareAbsCoefs)
524 {
525  SCIP_Real abscoef1;
526  SCIP_Real abscoef2;
527  SCIP_Real* coefs = (SCIP_Real*) dataptr;
528 
529  abscoef1 = REALABS(coefs[ind1]);
530  abscoef2 = REALABS(coefs[ind2]);
531 
532  if( abscoef1 < abscoef2 )
533  return -1;
534  if( abscoef2 < abscoef1 )
535  return 1;
536 
537  return 0;
538 }
539 
540 /** change given coefficient to new given value, adjust right hand side using the variables bound;
541  * returns TRUE if the right hand side would need to be changed to infinity and FALSE otherwise
542  */
543 static
545  SCIP* scip, /**< SCIP data structure */
546  SCIP_VAR* var, /**< variable the coefficient belongs to */
547  SCIP_Real oldcoeff, /**< old coefficient value */
548  SCIP_Real newcoeff, /**< new coefficient value */
549  SCIP_Bool cutislocal, /**< is the cut local? */
550  QUAD(SCIP_Real* cutrhs) /**< pointer to adjust right hand side of cut */
551  )
552 {
553  SCIP_Real QUAD(delta);
554  SCIPquadprecSumDD(delta, newcoeff, -oldcoeff);
555 
556  if( QUAD_TO_DBL(delta) > QUAD_EPSILON )
557  {
558  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
559  if( SCIPisInfinity(scip, ub) )
560  return TRUE;
561  else
562  {
563  SCIPquadprecProdQD(delta, delta, ub);
564  SCIPquadprecSumQQ(*cutrhs, *cutrhs, delta);
565  }
566  }
567  else if( QUAD_TO_DBL(delta) < -QUAD_EPSILON )
568  {
569  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
570  if( SCIPisInfinity(scip, -lb) )
571  return TRUE;
572  else
573  {
574  SCIPquadprecProdQD(delta, delta, lb);
575  SCIPquadprecSumQQ(*cutrhs, *cutrhs, delta);
576  }
577  }
578 
579  return FALSE;
580 }
581 
582 /** change given (quad) coefficient to new given value, adjust right hand side using the variables bound;
583  * returns TRUE if the right hand side would need to be changed to infinity and FALSE otherwise
584  */
585 static
587  SCIP* scip, /**< SCIP data structure */
588  SCIP_VAR* var, /**< variable the coefficient belongs to */
589  QUAD(SCIP_Real oldcoeff), /**< old coefficient value */
590  SCIP_Real newcoeff, /**< new coefficient value */
591  SCIP_Bool cutislocal, /**< is the cut local? */
592  QUAD(SCIP_Real* cutrhs) /**< pointer to adjust right hand side of cut */
593  )
594 {
595  SCIP_Real QUAD(delta);
596 
597  SCIPquadprecSumQD(delta, -oldcoeff, newcoeff);
598 
599  if( QUAD_TO_DBL(delta) > QUAD_EPSILON )
600  {
601  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
602  if( SCIPisInfinity(scip, ub) )
603  return TRUE;
604  else
605  {
606  SCIPquadprecProdQD(delta, delta, ub);
607  SCIPquadprecSumQQ(*cutrhs, *cutrhs, delta);
608  }
609  }
610  else if( QUAD_TO_DBL(delta) < -QUAD_EPSILON )
611  {
612  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
613  if( SCIPisInfinity(scip, -lb) )
614  return TRUE;
615  else
616  {
617  SCIPquadprecProdQD(delta, delta, lb);
618  SCIPquadprecSumQQ(*cutrhs, *cutrhs, delta);
619  }
620  }
621 
622  return FALSE;
623 }
624 
625 /** scales the cut and then tightens the coefficients of the given cut based on the maximal activity;
626  * see cons_linear.c consdataTightenCoefs() for details; the cut is given in a semi-sparse quad precision array;
627  */
628 static
630  SCIP* scip, /**< SCIP data structure */
631  SCIP_Bool cutislocal, /**< is the cut local? */
632  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
633  QUAD(SCIP_Real* cutrhs), /**< the right hand side of the cut */
634  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
635  int* cutnnz, /**< the number of non-zeros in the cut */
636  SCIP_Bool* redundant /**< whether the cut was detected to be redundant */
637  )
638 {
639  int i;
640  int nintegralvars;
641  SCIP_Bool isintegral;
642  SCIP_VAR** vars;
643  SCIP_Real QUAD(maxacttmp);
644  SCIP_Real maxact;
645  SCIP_Real maxabsintval;
646  SCIP_Real maxabscontval;
647 
648  QUAD_ASSIGN(maxacttmp, 0.0);
649 
650  vars = SCIPgetVars(scip);
651  maxabsintval = 0.0;
652  maxabscontval = 0.0;
653  nintegralvars = SCIPgetNVars(scip) - SCIPgetNContVars(scip);
654  isintegral = TRUE;
655 
656  *redundant = FALSE;
657 
658  /* compute the maximum activity and maximum absolute coefficient values for all and for integral variables in the cut */
659  for( i = 0; i < *cutnnz; ++i )
660  {
661  SCIP_Real QUAD(val);
662 
663  assert(cutinds[i] >= 0);
664  assert(vars[cutinds[i]] != NULL);
665 
666  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
667 
668  if( QUAD_TO_DBL(val) < 0.0 )
669  {
670  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
671 
672  if( SCIPisInfinity(scip, -lb) )
673  return SCIP_OKAY;
674 
675  if( cutinds[i] < nintegralvars )
676  maxabsintval = MAX(maxabsintval, -QUAD_TO_DBL(val));
677  else
678  {
679  maxabscontval = MAX(maxabscontval, -QUAD_TO_DBL(val));
680  isintegral = FALSE;
681  }
682 
683  SCIPquadprecProdQD(val, val, lb);
684 
685  SCIPquadprecSumQQ(maxacttmp, maxacttmp, val);
686  }
687  else
688  {
689  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
690 
691  if( SCIPisInfinity(scip, ub) )
692  return SCIP_OKAY;
693 
694  if( cutinds[i] < nintegralvars )
695  maxabsintval = MAX(maxabsintval, QUAD_TO_DBL(val));
696  else
697  {
698  maxabscontval = MAX(maxabscontval, QUAD_TO_DBL(val));
699  isintegral = FALSE;
700  }
701 
702  SCIPquadprecProdQD(val, val, ub);
703 
704  SCIPquadprecSumQQ(maxacttmp, maxacttmp, val);
705  }
706  }
707 
708  maxact = QUAD_TO_DBL(maxacttmp);
709 
710  /* cut is redundant in activity bounds */
711  if( SCIPisFeasLE(scip, maxact, QUAD_TO_DBL(*cutrhs)) )
712  {
713  *redundant = TRUE;
714  return SCIP_OKAY;
715  }
716 
717  /* cut is only on integral variables, try to scale to integral coefficients */
718  if( isintegral )
719  {
720  SCIP_Real equiscale;
721  SCIP_Real intscalar;
722  SCIP_Bool success;
723  SCIP_Real* intcoeffs;
724 
725  SCIP_CALL( SCIPallocBufferArray(scip, &intcoeffs, *cutnnz) );
726 
727  equiscale = 1.0 / MIN((maxact - QUAD_TO_DBL(*cutrhs)), maxabsintval);
728 
729  for( i = 0; i < *cutnnz; ++i )
730  {
731  SCIP_Real QUAD(val);
732 
733  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
734  SCIPquadprecProdQD(val, val, equiscale);
735 
736  intcoeffs[i] = QUAD_TO_DBL(val);
737  }
738 
739  SCIP_CALL( SCIPcalcIntegralScalar(intcoeffs, *cutnnz, -SCIPsumepsilon(scip), SCIPepsilon(scip),
740  (SCIP_Longint)scip->set->sepa_maxcoefratio, scip->set->sepa_maxcoefratio, &intscalar, &success) );
741 
742  SCIPfreeBufferArray(scip, &intcoeffs);
743 
744  if( success )
745  {
746  /* if successful, apply the scaling */
747  intscalar *= equiscale;
748 
749  SCIPquadprecProdQD(*cutrhs, *cutrhs, intscalar);
750 
751  for( i = 0; i < *cutnnz; )
752  {
753  SCIP_Real QUAD(val);
754  SCIP_Real intval;
755 
756  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
757  SCIPquadprecProdQD(val, val, intscalar);
758 
759  intval = SCIPround(scip, QUAD_TO_DBL(val));
760 
761  if( chgQuadCoeffWithBound(scip, vars[cutinds[i]], QUAD(val), intval, cutislocal, QUAD(cutrhs)) )
762  {
763  /* TODO maybe change the coefficient to the other value instead of discarding the cut? */
764  *redundant = TRUE;
765  return SCIP_OKAY;
766  }
767 
768  if( intval != 0.0 )
769  {
770  QUAD_ASSIGN(val, intval);
771  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], val);
772  ++i;
773  }
774  else
775  {
776  /* this must not be -0.0, otherwise the clean buffer memory is not cleared properly */
777  QUAD_ASSIGN(val, 0.0);
778  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], val);
779  --(*cutnnz);
780  cutinds[i] = cutinds[*cutnnz];
781  }
782  }
783 
784  SCIPquadprecEpsFloorQ(*cutrhs, *cutrhs, SCIPfeastol(scip)); /*lint !e666*/
785 
786  /* recompute the maximal activity after scaling to integral values */
787  QUAD_ASSIGN(maxacttmp, 0.0);
788  maxabsintval = 0.0;
789 
790  for( i = 0; i < *cutnnz; ++i )
791  {
792  SCIP_Real QUAD(val);
793 
794  assert(cutinds[i] >= 0);
795  assert(vars[cutinds[i]] != NULL);
796 
797  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
798 
799  if( QUAD_TO_DBL(val) < 0.0 )
800  {
801  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
802 
803  maxabsintval = MAX(maxabsintval, -QUAD_TO_DBL(val));
804 
805  SCIPquadprecProdQD(val, val, lb);
806 
807  SCIPquadprecSumQQ(maxacttmp, maxacttmp, val);
808  }
809  else
810  {
811  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
812 
813  maxabsintval = MAX(maxabsintval, QUAD_TO_DBL(val));
814 
815  SCIPquadprecProdQD(val, val, ub);
816 
817  SCIPquadprecSumQQ(maxacttmp, maxacttmp, val);
818  }
819  }
820 
821  maxact = QUAD_TO_DBL(maxacttmp);
822 
823  assert(EPSISINT(maxact, 1e-4));
824  maxact = SCIPround(scip, maxact);
825  QUAD_ASSIGN(maxacttmp, maxact);
826 
827  /* check again for redundancy */
828  if( SCIPisFeasLE(scip, maxact, QUAD_TO_DBL(*cutrhs)) )
829  {
830  *redundant = TRUE;
831  return SCIP_OKAY;
832  }
833  }
834  else
835  {
836  /* otherwise, apply the equilibrium scaling */
837  isintegral = FALSE;
838 
839  /* perform the scaling */
840  SCIPquadprecProdQD(maxacttmp, maxacttmp, equiscale);
841 
842  SCIPquadprecProdQD(*cutrhs, *cutrhs, equiscale);
843  maxabsintval *= equiscale;
844 
845  for( i = 0; i < *cutnnz; ++i )
846  {
847  SCIP_Real QUAD(val);
848 
849  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
850  SCIPquadprecProdQD(val, val, equiscale);
851  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], val);
852  }
853  }
854  }
855  else
856  {
857  /* cut has integer and continuous variables, so scale it to equilibrium */
858  SCIP_Real scale;
859  SCIP_Real maxabsval;
860 
861  maxabsval = maxact - QUAD_TO_DBL(*cutrhs);
862  maxabsval = MIN(maxabsval, maxabsintval);
863  maxabsval = MAX(maxabsval, maxabscontval);
864 
865  scale = 1.0 / maxabsval; /*lint !e795*/
866 
867  /* perform the scaling */
868  SCIPquadprecProdQD(maxacttmp, maxacttmp, scale);
869  maxact = QUAD_TO_DBL(maxacttmp);
870 
871  SCIPquadprecProdQD(*cutrhs, *cutrhs, scale);
872  maxabsintval *= scale;
873 
874  for( i = 0; i < *cutnnz; ++i )
875  {
876  SCIP_Real QUAD(val);
877 
878  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
879  SCIPquadprecProdQD(val, val, scale);
880  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], val);
881  }
882  }
883 
884  /* no coefficient tightening can be performed since the precondition doesn't hold for any of the variables */
885  if( SCIPisGT(scip, maxact - maxabsintval, QUAD_TO_DBL(*cutrhs)) )
886  return SCIP_OKAY;
887 
888  SCIPsortDownInd(cutinds, compareAbsCoefsQuad, (void*) cutcoefs, *cutnnz);
889 
890  /* loop over the integral variables and try to tighten the coefficients; see cons_linear for more details */
891  for( i = 0; i < *cutnnz; )
892  {
893  SCIP_Real QUAD(val);
894 
895  if( cutinds[i] >= nintegralvars )
896  {
897  ++i;
898  continue;
899  }
900 
901  QUAD_ARRAY_LOAD(val, cutcoefs, cutinds[i]);
902 
903  assert(SCIPvarIsIntegral(vars[cutinds[i]]));
904 
905  if( QUAD_TO_DBL(val) < 0.0 && SCIPisLE(scip, maxact + QUAD_TO_DBL(val), QUAD_TO_DBL(*cutrhs)) )
906  {
907  SCIP_Real QUAD(coef);
908  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
909 
910  SCIPquadprecSumQQ(coef, *cutrhs, -maxacttmp);
911 
912  if( isintegral )
913  {
914  /* if the cut is integral, the true coefficient must also be integral;
915  * thus we round it to the exact integral value
916  */
917  assert(SCIPisFeasIntegral(scip, QUAD_TO_DBL(coef)));
918  QUAD_ASSIGN(coef, SCIPround(scip, QUAD_TO_DBL(coef)));
919  }
920 
921  if( QUAD_TO_DBL(coef) > QUAD_TO_DBL(val) )
922  {
923  SCIP_Real QUAD(delta);
924  SCIP_Real QUAD(tmp);
925 
926  SCIPquadprecSumQQ(delta, -val, coef);
927  SCIPquadprecProdQD(delta, delta, lb);
928 
929  SCIPquadprecSumQQ(tmp, delta, *cutrhs);
930 
931  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
932  QUAD_TO_DBL(val), QUAD_TO_DBL(coef), QUAD_TO_DBL(*cutrhs), QUAD_TO_DBL(tmp), lb,
933  cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]));
934 
935  QUAD_ASSIGN_Q(*cutrhs, tmp);
936 
937  assert(!SCIPisPositive(scip, QUAD_TO_DBL(coef)));
938 
939  if( SCIPisNegative(scip, QUAD_TO_DBL(coef)) )
940  {
941  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
942  maxact = QUAD_TO_DBL(maxacttmp);
943  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], coef);
944  }
945  else
946  {
947  QUAD_ASSIGN(coef, 0.0);
948  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], coef);
949  --(*cutnnz);
950  cutinds[i] = cutinds[*cutnnz];
951  continue;
952  }
953  }
954  }
955  else if( QUAD_TO_DBL(val) > 0.0 && SCIPisLE(scip, maxact - QUAD_TO_DBL(val), QUAD_TO_DBL(*cutrhs)) )
956  {
957  SCIP_Real QUAD(coef);
958  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
959 
960  SCIPquadprecSumQQ(coef, maxacttmp, -*cutrhs);
961 
962  if( isintegral )
963  {
964  /* if the cut is integral, the true coefficient must also be integral;
965  * thus we round it to the exact integral value
966  */
967  assert(SCIPisFeasIntegral(scip, QUAD_TO_DBL(coef)));
968  QUAD_ASSIGN(coef, SCIPround(scip, QUAD_TO_DBL(coef)));
969  }
970 
971  if( QUAD_TO_DBL(coef) < QUAD_TO_DBL(val) )
972  {
973  SCIP_Real QUAD(delta);
974  SCIP_Real QUAD(tmp);
975 
976  SCIPquadprecSumQQ(delta, -val, coef);
977  SCIPquadprecProdQD(delta, delta, ub);
978 
979  SCIPquadprecSumQQ(tmp, delta, *cutrhs);
980 
981  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
982  QUAD_TO_DBL(val), QUAD_TO_DBL(coef), QUAD_TO_DBL(*cutrhs), QUAD_TO_DBL(tmp),
983  cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]), ub);
984 
985  QUAD_ASSIGN_Q(*cutrhs, tmp);
986 
987  assert(SCIPisGE(scip, QUAD_TO_DBL(coef), 0.0));
988  if( SCIPisPositive(scip, QUAD_TO_DBL(coef)) )
989  {
990  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
991  maxact = QUAD_TO_DBL(maxacttmp);
992  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], coef);
993  }
994  else
995  {
996  QUAD_ASSIGN(coef, 0.0);
997  QUAD_ARRAY_STORE(cutcoefs, cutinds[i], coef);
998  --(*cutnnz);
999  cutinds[i] = cutinds[*cutnnz];
1000  continue;
1001  }
1002  }
1003  }
1004  else /* due to sorting we can stop completely if the precondition was not fulfilled for this variable */
1005  break;
1006 
1007  ++i;
1008  }
1009 
1010  return SCIP_OKAY;
1011 }
1012 
1013 /** scales the cut and then tightens the coefficients of the given cut based on the maximal activity;
1014  * see cons_linear.c consdataTightenCoefs() for details; the cut is given in a semi-sparse array;
1015  */
1016 static
1018  SCIP* scip, /**< SCIP data structure */
1019  SCIP_Bool cutislocal, /**< is the cut local? */
1020  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
1021  QUAD(SCIP_Real* cutrhs), /**< the right hand side of the cut */
1022  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
1023  int* cutnnz, /**< the number of non-zeros in the cut */
1024  SCIP_Bool* redundant /**< pointer to return whtether the cut was detected to be redundant */
1025  )
1026 {
1027  int i;
1028  int nintegralvars;
1029  SCIP_Bool isintegral;
1030  SCIP_VAR** vars;
1031  SCIP_Real QUAD(maxacttmp);
1032  SCIP_Real maxact;
1033  SCIP_Real maxabsintval;
1034  SCIP_Real maxabscontval;
1035 
1036  QUAD_ASSIGN(maxacttmp, 0.0);
1037 
1038  vars = SCIPgetVars(scip);
1039  maxabsintval = 0.0;
1040  maxabscontval = 0.0;
1041  isintegral = TRUE;
1042  *redundant = FALSE;
1043  nintegralvars = SCIPgetNVars(scip) - SCIPgetNContVars(scip);
1044 
1045  for( i = 0; i < *cutnnz; ++i )
1046  {
1047  SCIP_Real val;
1048 
1049  assert(cutinds[i] >= 0);
1050  assert(vars[cutinds[i]] != NULL);
1051 
1052  val = cutcoefs[cutinds[i]];
1053 
1054  if( val < 0.0 )
1055  {
1056  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
1057 
1058  if( SCIPisInfinity(scip, -lb) )
1059  return SCIP_OKAY;
1060 
1061  if( cutinds[i] < nintegralvars )
1062  maxabsintval = MAX(maxabsintval, -val);
1063  else
1064  {
1065  maxabscontval = MAX(maxabscontval, -val);
1066  isintegral = FALSE;
1067  }
1068 
1069  SCIPquadprecSumQD(maxacttmp, maxacttmp, val * lb);
1070  }
1071  else
1072  {
1073  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
1074 
1075  if( SCIPisInfinity(scip, ub) )
1076  return SCIP_OKAY;
1077 
1078  if( cutinds[i] < nintegralvars )
1079  maxabsintval = MAX(maxabsintval, val);
1080  else
1081  {
1082  maxabscontval = MAX(maxabscontval, val);
1083  isintegral = FALSE;
1084  }
1085 
1086  SCIPquadprecSumQD(maxacttmp, maxacttmp, val * ub);
1087  }
1088  }
1089 
1090  maxact = QUAD_TO_DBL(maxacttmp);
1091 
1092  /* cut is redundant in activity bounds */
1093  if( SCIPisFeasLE(scip, maxact, QUAD_TO_DBL(*cutrhs)) )
1094  {
1095  *redundant = TRUE;
1096  return SCIP_OKAY;
1097  }
1098 
1099  /* cut is only on integral variables, try to scale to integral coefficients */
1100  if( isintegral )
1101  {
1102  SCIP_Real equiscale;
1103  SCIP_Real intscalar;
1104  SCIP_Bool success;
1105  SCIP_Real* intcoeffs;
1106 
1107  SCIP_CALL( SCIPallocBufferArray(scip, &intcoeffs, *cutnnz) );
1108 
1109  equiscale = 1.0 / MIN((maxact - QUAD_TO_DBL(*cutrhs)), maxabsintval);
1110 
1111  for( i = 0; i < *cutnnz; ++i )
1112  {
1113  SCIP_Real scaleval;
1114  SCIP_Real val;
1115 
1116  val = cutcoefs[cutinds[i]];
1117 
1118  scaleval = val * equiscale;
1119 
1120  intcoeffs[i] = scaleval;
1121  }
1122 
1123  SCIP_CALL( SCIPcalcIntegralScalar(intcoeffs, *cutnnz, -SCIPsumepsilon(scip), SCIPepsilon(scip),
1124  (SCIP_Longint)scip->set->sepa_maxcoefratio, scip->set->sepa_maxcoefratio, &intscalar, &success) );
1125 
1126  SCIPfreeBufferArray(scip, &intcoeffs);
1127 
1128  if( success )
1129  {
1130  /* if successful, apply the scaling */
1131  intscalar *= equiscale;
1132 
1133  SCIPquadprecProdQD(*cutrhs, *cutrhs, intscalar);
1134 
1135  for( i = 0; i < *cutnnz; )
1136  {
1137  SCIP_Real val;
1138  SCIP_Real intval;
1139 
1140  val = cutcoefs[cutinds[i]];
1141  val *= intscalar;
1142 
1143  intval = SCIPround(scip, val);
1144 
1145  if( chgCoeffWithBound(scip, vars[cutinds[i]], val, intval, cutislocal, QUAD(cutrhs)) )
1146  {
1147  /* TODO maybe change the coefficient to the other value instead of discarding the cut? */
1148  *redundant = TRUE;
1149  return SCIP_OKAY;
1150  }
1151 
1152  if( intval != 0.0 )
1153  {
1154  cutcoefs[cutinds[i]] = intval;
1155  ++i;
1156  }
1157  else
1158  {
1159  /* this must not be -0.0, otherwise the clean buffer memory is not cleared properly */
1160  cutcoefs[cutinds[i]] = 0.0;
1161  --(*cutnnz);
1162  cutinds[i] = cutinds[*cutnnz];
1163  }
1164  }
1165 
1166  SCIPquadprecEpsFloorQ(*cutrhs, *cutrhs, SCIPfeastol(scip)); /*lint !e666*/
1167 
1168  /* recompute the maximal activity after scaling to integral values */
1169  QUAD_ASSIGN(maxacttmp, 0.0);
1170  maxabsintval = 0.0;
1171 
1172  for( i = 0; i < *cutnnz; ++i )
1173  {
1174  SCIP_Real val;
1175 
1176  assert(cutinds[i] >= 0);
1177  assert(vars[cutinds[i]] != NULL);
1178 
1179  val = cutcoefs[cutinds[i]];
1180 
1181  if( val < 0.0 )
1182  {
1183  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
1184 
1185  maxabsintval = MAX(maxabsintval, -val);
1186 
1187  val *= lb;
1188 
1189  SCIPquadprecSumQD(maxacttmp, maxacttmp, val);
1190  }
1191  else
1192  {
1193  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
1194 
1195  maxabsintval = MAX(maxabsintval, val);
1196 
1197  val *= ub;
1198 
1199  SCIPquadprecSumQD(maxacttmp, maxacttmp, val);
1200  }
1201  }
1202 
1203  maxact = QUAD_TO_DBL(maxacttmp);
1204 
1205  assert(EPSISINT(maxact, 1e-4));
1206  maxact = SCIPround(scip, maxact);
1207  QUAD_ASSIGN(maxacttmp, maxact);
1208 
1209  /* check again for redundancy */
1210  if( SCIPisFeasLE(scip, maxact, QUAD_TO_DBL(*cutrhs)) )
1211  {
1212  *redundant = TRUE;
1213  return SCIP_OKAY;
1214  }
1215  }
1216  else
1217  {
1218  /* otherwise, apply the equilibrium scaling */
1219  isintegral = FALSE;
1220 
1221  /* perform the scaling */
1222  SCIPquadprecProdQD(maxacttmp, maxacttmp, equiscale);
1223 
1224  SCIPquadprecProdQD(*cutrhs, *cutrhs, equiscale);
1225  maxabsintval *= equiscale;
1226 
1227  for( i = 0; i < *cutnnz; ++i )
1228  cutcoefs[cutinds[i]] *= equiscale;
1229  }
1230  }
1231  else
1232  {
1233  /* cut has integer and continuous variables, so scale it to equilibrium */
1234  SCIP_Real scale;
1235  SCIP_Real maxabsval;
1236 
1237  maxabsval = maxact - QUAD_TO_DBL(*cutrhs);
1238  maxabsval = MIN(maxabsval, maxabsintval);
1239  maxabsval = MAX(maxabsval, maxabscontval);
1240 
1241  scale = 1.0 / maxabsval; /*lint !e795*/
1242 
1243  /* perform the scaling */
1244  SCIPquadprecProdQD(maxacttmp, maxacttmp, scale);
1245  maxact = QUAD_TO_DBL(maxacttmp);
1246 
1247  SCIPquadprecProdQD(*cutrhs, *cutrhs, scale);
1248  maxabsintval *= scale;
1249 
1250  for( i = 0; i < *cutnnz; ++i )
1251  cutcoefs[cutinds[i]] *= scale;
1252  }
1253 
1254  /* no coefficient tightening can be performed since the precondition doesn't hold for any of the variables */
1255  if( SCIPisGT(scip, maxact - maxabsintval, QUAD_TO_DBL(*cutrhs)) )
1256  return SCIP_OKAY;
1257 
1258  SCIPsortDownInd(cutinds, compareAbsCoefs, (void*) cutcoefs, *cutnnz);
1259 
1260  /* loop over the integral variables and try to tighten the coefficients; see cons_linear for more details */
1261  for( i = 0; i < *cutnnz; )
1262  {
1263  SCIP_Real val;
1264 
1265  if( cutinds[i] >= nintegralvars )
1266  {
1267  ++i;
1268  continue;
1269  }
1270 
1271  val = cutcoefs[cutinds[i]];
1272 
1273  assert(SCIPvarIsIntegral(vars[cutinds[i]]));
1274 
1275  if( val < 0.0 && SCIPisLE(scip, maxact + val, QUAD_TO_DBL(*cutrhs)) )
1276  {
1277  SCIP_Real QUAD(coef);
1278  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
1279 
1280  SCIPquadprecSumQQ(coef, -maxacttmp, *cutrhs);
1281 
1282  if( isintegral )
1283  {
1284  /* if the cut is integral, the true coefficient must also be integral;
1285  * thus we round it to the exact integral value
1286  */
1287  assert(SCIPisFeasIntegral(scip, QUAD_TO_DBL(coef)));
1288  QUAD_ASSIGN(coef, SCIPround(scip, QUAD_TO_DBL(coef)));
1289  }
1290 
1291  if( QUAD_TO_DBL(coef) > val )
1292  {
1293  SCIP_Real QUAD(delta);
1294  SCIP_Real QUAD(tmp);
1295 
1296  SCIPquadprecSumQD(delta, coef, -val);
1297  SCIPquadprecProdQD(delta, delta, lb);
1298 
1299  SCIPquadprecSumQQ(tmp, delta, *cutrhs);
1300  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
1301  val, QUAD_TO_DBL(coef), QUAD_TO_DBL(*cutrhs), QUAD_TO_DBL(tmp), lb,
1302  cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]));
1303 
1304  QUAD_ASSIGN_Q(*cutrhs, tmp);
1305 
1306  assert(!SCIPisPositive(scip, QUAD_TO_DBL(coef)));
1307 
1308  if( SCIPisNegative(scip, QUAD_TO_DBL(coef)) )
1309  {
1310  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
1311  maxact = QUAD_TO_DBL(maxacttmp);
1312  cutcoefs[cutinds[i]] = QUAD_TO_DBL(coef);
1313  }
1314  else
1315  {
1316  cutcoefs[cutinds[i]] = 0.0;
1317  --(*cutnnz);
1318  cutinds[i] = cutinds[*cutnnz];
1319  continue;
1320  }
1321  }
1322  }
1323  else if( val > 0.0 && SCIPisLE(scip, maxact - val, QUAD_TO_DBL(*cutrhs)) )
1324  {
1325  SCIP_Real QUAD(coef);
1326  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
1327 
1328  SCIPquadprecSumQQ(coef, maxacttmp, -*cutrhs);
1329 
1330  if( isintegral )
1331  {
1332  /* if the cut is integral, the true coefficient must also be integral;
1333  * thus we round it to the exact integral value
1334  */
1335  assert(SCIPisFeasIntegral(scip, QUAD_TO_DBL(coef)));
1336  QUAD_ASSIGN(coef, SCIPround(scip, QUAD_TO_DBL(coef)));
1337  }
1338 
1339  if( QUAD_TO_DBL(coef) < val )
1340  {
1341  SCIP_Real QUAD(delta);
1342  SCIP_Real QUAD(tmp);
1343 
1344  SCIPquadprecSumQD(delta, coef, -val);
1345  SCIPquadprecProdQD(delta, delta, ub);
1346 
1347  SCIPquadprecSumQQ(tmp, delta, *cutrhs);
1348  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
1349  val, QUAD_TO_DBL(coef), QUAD_TO_DBL(*cutrhs), QUAD_TO_DBL(tmp),
1350  cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]), ub);
1351 
1352  QUAD_ASSIGN_Q(*cutrhs, tmp);
1353 
1354  assert(! SCIPisNegative(scip, QUAD_TO_DBL(coef)));
1355 
1356  if( SCIPisPositive(scip, QUAD_TO_DBL(coef)) )
1357  {
1358  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
1359  maxact = QUAD_TO_DBL(maxacttmp);
1360  cutcoefs[cutinds[i]] = QUAD_TO_DBL(coef);
1361  }
1362  else
1363  {
1364  cutcoefs[cutinds[i]] = 0.0;
1365  --(*cutnnz);
1366  cutinds[i] = cutinds[*cutnnz];
1367  continue;
1368  }
1369  }
1370  }
1371  else /* due to sorting we can stop completely if the precondition was not fulfilled for this variable */
1372  break;
1373 
1374  ++i;
1375  }
1376 
1377  return SCIP_OKAY;
1378 }
1379 
1380 /** perform activity based coefficient tightening on the given cut; returns TRUE if the cut was detected
1381  * to be redundant due to activity bounds
1382  */
1384  SCIP* scip, /**< SCIP data structure */
1385  SCIP_Bool cutislocal, /**< is the cut local? */
1386  SCIP_Real* cutcoefs, /**< array of the non-zero coefficients in the cut */
1387  SCIP_Real* cutrhs, /**< the right hand side of the cut */
1388  int* cutinds, /**< array of the problem indices of variables with a non-zero coefficient in the cut */
1389  int* cutnnz, /**< the number of non-zeros in the cut */
1390  int* nchgcoefs /**< number of changed coefficients */
1391  )
1392 {
1393  int i;
1394  int nintegralvars;
1395  SCIP_VAR** vars;
1396  SCIP_Real* absvals;
1397  SCIP_Real QUAD(maxacttmp);
1398  SCIP_Real maxact;
1399  SCIP_Real maxabsval;
1400  SCIP_Bool redundant;
1401 
1402  assert(nchgcoefs != NULL);
1403 
1404  QUAD_ASSIGN(maxacttmp, 0.0);
1405 
1406  vars = SCIPgetVars(scip);
1407  nintegralvars = SCIPgetNVars(scip) - SCIPgetNContVars(scip);
1408  maxabsval = 0.0;
1409  SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &absvals, *cutnnz) );
1410 
1411  *nchgcoefs = 0;
1412  redundant = FALSE;
1413 
1414  for( i = 0; i < *cutnnz; ++i )
1415  {
1416  assert(cutinds[i] >= 0);
1417  assert(vars[cutinds[i]] != NULL);
1418 
1419  if( cutcoefs[i] < 0.0 )
1420  {
1421  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
1422 
1423  if( SCIPisInfinity(scip, -lb) )
1424  goto TERMINATE;
1425 
1426  if( cutinds[i] < nintegralvars )
1427  {
1428  maxabsval = MAX(maxabsval, -cutcoefs[i]);
1429  absvals[i] = -cutcoefs[i];
1430  }
1431  else
1432  {
1433  absvals[i] = 0.0;
1434  }
1435 
1436  SCIPquadprecSumQD(maxacttmp, maxacttmp, lb * cutcoefs[i]);
1437  }
1438  else
1439  {
1440  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
1441 
1442  if( SCIPisInfinity(scip, ub) )
1443  goto TERMINATE;
1444 
1445  if( cutinds[i] < nintegralvars )
1446  {
1447  maxabsval = MAX(maxabsval, cutcoefs[i]);
1448  absvals[i] = cutcoefs[i];
1449  }
1450  else
1451  {
1452  absvals[i] = 0.0;
1453  }
1454 
1455  SCIPquadprecSumQD(maxacttmp, maxacttmp, ub * cutcoefs[i]);
1456  }
1457  }
1458 
1459  maxact = QUAD_TO_DBL(maxacttmp);
1460 
1461  /* cut is redundant in activity bounds */
1462  if( SCIPisFeasLE(scip, maxact, *cutrhs) )
1463  {
1464  redundant = TRUE;
1465  goto TERMINATE;
1466  }
1467 
1468  /* no coefficient tightening can be performed since the precondition doesn't hold for any of the variables */
1469  if( SCIPisGT(scip, maxact - maxabsval, *cutrhs) )
1470  goto TERMINATE;
1471 
1472  SCIPsortDownRealRealInt(absvals, cutcoefs, cutinds, *cutnnz);
1473  SCIPfreeBufferArray(scip, &absvals);
1474 
1475  /* loop over the integral variables and try to tighten the coefficients; see cons_linear for more details */
1476  for( i = 0; i < *cutnnz;)
1477  {
1478  if( cutinds[i] >= nintegralvars )
1479  {
1480  ++i;
1481  continue;
1482  }
1483 
1484  assert(SCIPvarIsIntegral(vars[cutinds[i]]));
1485 
1486  if( cutcoefs[i] < 0.0 && SCIPisLE(scip, maxact + cutcoefs[i], *cutrhs) )
1487  {
1488  SCIP_Real coef = (*cutrhs) - maxact;
1489  SCIP_Real lb = cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]);
1490 
1491  coef = floor(coef);
1492 
1493  if( coef > cutcoefs[i] )
1494  {
1495  SCIP_Real QUAD(delta);
1496  SCIP_Real QUAD(tmp);
1497 
1498  SCIPquadprecSumDD(delta, coef, -cutcoefs[i]);
1499  SCIPquadprecProdQD(delta, delta, lb);
1500 
1501  SCIPquadprecSumQD(tmp, delta, *cutrhs);
1502  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
1503  cutcoefs[i], coef, (*cutrhs), QUAD_TO_DBL(tmp), lb,
1504  cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]));
1505 
1506  *cutrhs = QUAD_TO_DBL(tmp);
1507 
1508  assert(!SCIPisPositive(scip, coef));
1509 
1510  ++(*nchgcoefs);
1511 
1512  if( SCIPisNegative(scip, coef) )
1513  {
1514  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
1515  maxact = QUAD_TO_DBL(maxacttmp);
1516  cutcoefs[i] = coef;
1517  }
1518  else
1519  {
1520  --(*cutnnz);
1521  cutinds[i] = cutinds[*cutnnz];
1522  cutcoefs[i] = cutcoefs[*cutnnz];
1523  continue;
1524  }
1525  }
1526  }
1527  else if( cutcoefs[i] > 0.0 && SCIPisLE(scip, maxact - cutcoefs[i], *cutrhs) )
1528  {
1529  SCIP_Real coef = maxact - (*cutrhs);
1530  SCIP_Real ub = cutislocal ? SCIPvarGetUbLocal(vars[cutinds[i]]) : SCIPvarGetUbGlobal(vars[cutinds[i]]);
1531 
1532  coef = ceil(coef);
1533 
1534  if( coef < cutcoefs[i] )
1535  {
1536  SCIP_Real QUAD(delta);
1537  SCIP_Real QUAD(tmp);
1538 
1539  SCIPquadprecSumDD(delta, coef, -cutcoefs[i]);
1540  SCIPquadprecProdQD(delta, delta, ub);
1541 
1542  SCIPquadprecSumQD(tmp, delta, *cutrhs);
1543  SCIPdebugPrintf("tightened coefficient from %g to %g; rhs changed from %g to %g; the bounds are [%g,%g]\n",
1544  cutcoefs[i], coef, (*cutrhs), QUAD_TO_DBL(tmp),
1545  cutislocal ? SCIPvarGetLbLocal(vars[cutinds[i]]) : SCIPvarGetLbGlobal(vars[cutinds[i]]), ub);
1546 
1547  *cutrhs = QUAD_TO_DBL(tmp);
1548 
1549  assert(!SCIPisNegative(scip, coef));
1550 
1551  ++(*nchgcoefs);
1552 
1553  if( SCIPisPositive(scip, coef) )
1554  {
1555  SCIPquadprecSumQQ(maxacttmp, maxacttmp, delta);
1556  maxact = QUAD_TO_DBL(maxacttmp);
1557  cutcoefs[i] = coef;
1558  }
1559  else
1560  {
1561  --(*cutnnz);
1562  cutinds[i] = cutinds[*cutnnz];
1563  cutcoefs[i] = cutcoefs[*cutnnz];
1564  continue;
1565  }
1566  }
1567  }
1568  else /* due to sorting we can stop completely if the precondition was not fulfilled for this variable */
1569  break;
1570 
1571  ++i;
1572  }
1573 
1574  TERMINATE:
1575  SCIPfreeBufferArrayNull(scip, &absvals);
1576 
1577  return redundant;
1578 }
1579 
1580 /* =========================================== aggregation row =========================================== */
1581 
1582 
1583 /** create an empty aggregation row */
1585  SCIP* scip, /**< SCIP data structure */
1586  SCIP_AGGRROW** aggrrow /**< pointer to return aggregation row */
1587  )
1588 {
1589  int nvars;
1590  assert(scip != NULL);
1591  assert(aggrrow != NULL);
1592 
1593  SCIP_CALL( SCIPallocBlockMemory(scip, aggrrow) );
1594 
1595  nvars = SCIPgetNVars(scip);
1596 
1597  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*aggrrow)->vals, QUAD_ARRAY_SIZE(nvars)) );
1598  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*aggrrow)->inds, nvars) );
1599 
1600  BMSclearMemoryArray((*aggrrow)->vals, QUAD_ARRAY_SIZE(nvars));
1601 
1602  (*aggrrow)->local = FALSE;
1603  (*aggrrow)->nnz = 0;
1604  (*aggrrow)->rank = 0;
1605  QUAD_ASSIGN((*aggrrow)->rhs, 0.0);
1606  (*aggrrow)->rowsinds = NULL;
1607  (*aggrrow)->slacksign = NULL;
1608  (*aggrrow)->rowweights = NULL;
1609  (*aggrrow)->nrows = 0;
1610  (*aggrrow)->rowssize = 0;
1611 
1612  return SCIP_OKAY;
1613 }
1614 
1615 /** free a aggregation row */
1617  SCIP* scip, /**< SCIP data structure */
1618  SCIP_AGGRROW** aggrrow /**< pointer to aggregation row that should be freed */
1619  )
1620 {
1621  int nvars;
1622  assert(scip != NULL);
1623  assert(aggrrow != NULL);
1624 
1625  nvars = SCIPgetNVars(scip);
1626 
1627  SCIPfreeBlockMemoryArray(scip, &(*aggrrow)->inds, nvars);
1628  SCIPfreeBlockMemoryArray(scip, &(*aggrrow)->vals, QUAD_ARRAY_SIZE(nvars)); /*lint !e647*/
1629  SCIPfreeBlockMemoryArrayNull(scip, &(*aggrrow)->rowsinds, (*aggrrow)->rowssize);
1630  SCIPfreeBlockMemoryArrayNull(scip, &(*aggrrow)->slacksign, (*aggrrow)->rowssize);
1631  SCIPfreeBlockMemoryArrayNull(scip, &(*aggrrow)->rowweights, (*aggrrow)->rowssize);
1632  SCIPfreeBlockMemory(scip, aggrrow);
1633 }
1634 
1635 /** output aggregation row to file stream */
1637  SCIP* scip, /**< SCIP data structure */
1638  SCIP_AGGRROW* aggrrow, /**< pointer to return aggregation row */
1639  FILE* file /**< output file (or NULL for standard output) */
1640  )
1641 {
1642  SCIP_VAR** vars;
1643  SCIP_MESSAGEHDLR* messagehdlr;
1644  int i;
1645 
1646  assert(scip != NULL);
1647  assert(aggrrow != NULL);
1648 
1649  vars = SCIPgetVars(scip);
1650  assert(vars != NULL);
1651 
1652  messagehdlr = SCIPgetMessagehdlr(scip);
1653  assert(messagehdlr);
1654 
1655  /* print coefficients */
1656  if( aggrrow->nnz == 0 )
1657  SCIPmessageFPrintInfo(messagehdlr, file, "0 ");
1658 
1659  for( i = 0; i < aggrrow->nnz; ++i )
1660  {
1661  SCIP_Real QUAD(val);
1662 
1663  QUAD_ARRAY_LOAD(val, aggrrow->vals, aggrrow->inds[i]);
1664  assert(SCIPvarGetProbindex(vars[aggrrow->inds[i]]) == aggrrow->inds[i]);
1665  SCIPmessageFPrintInfo(messagehdlr, file, "%+.15g<%s> ", QUAD_TO_DBL(val), SCIPvarGetName(vars[aggrrow->inds[i]]));
1666  }
1667 
1668  /* print right hand side */
1669  SCIPmessageFPrintInfo(messagehdlr, file, "<= %.15g\n", QUAD_TO_DBL(aggrrow->rhs));
1670 }
1671 
1672 /** copy a aggregation row */
1674  SCIP* scip, /**< SCIP data structure */
1675  SCIP_AGGRROW** aggrrow, /**< pointer to return aggregation row */
1676  SCIP_AGGRROW* source /**< source aggregation row */
1677  )
1678 {
1679  int nvars;
1680 
1681  assert(scip != NULL);
1682  assert(aggrrow != NULL);
1683  assert(source != NULL);
1684 
1685  nvars = SCIPgetNVars(scip);
1686  SCIP_CALL( SCIPallocBlockMemory(scip, aggrrow) );
1687 
1688  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*aggrrow)->vals, source->vals, QUAD_ARRAY_SIZE(nvars)) );
1689  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*aggrrow)->inds, source->inds, nvars) );
1690  (*aggrrow)->nnz = source->nnz;
1691  QUAD_ASSIGN_Q((*aggrrow)->rhs, source->rhs);
1692 
1693  if( source->nrows > 0 )
1694  {
1695  assert(source->rowsinds != NULL);
1696  assert(source->slacksign != NULL);
1697  assert(source->rowweights != NULL);
1698 
1699  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*aggrrow)->rowsinds, source->rowsinds, source->nrows) );
1700  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*aggrrow)->slacksign, source->slacksign, source->nrows) );
1701  SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*aggrrow)->rowweights, source->rowweights, source->nrows) );
1702  }
1703  else
1704  {
1705  (*aggrrow)->rowsinds = NULL;
1706  (*aggrrow)->slacksign = NULL;
1707  (*aggrrow)->rowweights = NULL;
1708  }
1709 
1710  (*aggrrow)->nrows = source->nrows;
1711  (*aggrrow)->rowssize = source->nrows;
1712  (*aggrrow)->rank = source->rank;
1713  (*aggrrow)->local = source->local;
1714 
1715  return SCIP_OKAY;
1716 }
1717 
1718 /** add weighted row to aggregation row */
1720  SCIP* scip, /**< SCIP data structure */
1721  SCIP_AGGRROW* aggrrow, /**< aggregation row */
1722  SCIP_ROW* row, /**< row to add to aggregation row */
1723  SCIP_Real weight, /**< scale for adding given row to aggregation row */
1724  int sidetype /**< specify row side type (-1 = lhs, 0 = automatic, 1 = rhs) */
1725  )
1726 {
1727  int i;
1728 
1729  assert(row->lppos >= 0);
1730 
1731  /* update local flag */
1732  aggrrow->local = aggrrow->local || row->local;
1733 
1734  /* update rank */
1735  aggrrow->rank = MAX(row->rank, aggrrow->rank);
1736 
1737  {
1738  SCIP_Real sideval;
1739  SCIP_Bool uselhs;
1740 
1741  i = aggrrow->nrows++;
1742 
1743  if( aggrrow->nrows > aggrrow->rowssize )
1744  {
1745  int newsize = SCIPcalcMemGrowSize(scip, aggrrow->nrows);
1746  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->rowsinds, aggrrow->rowssize, newsize) );
1747  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->slacksign, aggrrow->rowssize, newsize) );
1748  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->rowweights, aggrrow->rowssize, newsize) );
1749  aggrrow->rowssize = newsize;
1750  }
1751  aggrrow->rowsinds[i] = SCIProwGetLPPos(row);
1752  aggrrow->rowweights[i] = weight;
1753 
1754  if ( sidetype == -1 )
1755  {
1756  assert( ! SCIPisInfinity(scip, -row->lhs) );
1757  uselhs = TRUE;
1758  }
1759  else if ( sidetype == 1 )
1760  {
1761  assert( ! SCIPisInfinity(scip, row->rhs) );
1762  uselhs = FALSE;
1763  }
1764  else
1765  {
1766  /* Automatically decide, whether we want to use the left or the right hand side of the row in the summation.
1767  * If possible, use the side that leads to a positive slack value in the summation.
1768  */
1769  if( SCIPisInfinity(scip, row->rhs) || (!SCIPisInfinity(scip, -row->lhs) && weight < 0.0) )
1770  uselhs = TRUE;
1771  else
1772  uselhs = FALSE;
1773  }
1774 
1775  if( uselhs )
1776  {
1777  aggrrow->slacksign[i] = -1;
1778  sideval = row->lhs - row->constant;
1779  if( row->integral )
1780  sideval = SCIPceil(scip, sideval); /* row is integral: round left hand side up */
1781  }
1782  else
1783  {
1784  aggrrow->slacksign[i] = +1;
1785  sideval = row->rhs - row->constant;
1786  if( row->integral )
1787  sideval = SCIPfloor(scip, sideval); /* row is integral: round right hand side up */
1788  }
1789 
1790  SCIPquadprecSumQD(aggrrow->rhs, aggrrow->rhs, weight * sideval);
1791  }
1792 
1793  /* add up coefficients */
1794  SCIP_CALL( varVecAddScaledRowCoefsQuad(aggrrow->inds, aggrrow->vals, &aggrrow->nnz, row, weight) );
1795 
1796  return SCIP_OKAY;
1797 }
1798 
1799 /** Removes a given variable @p var from position @p pos the aggregation row and updates the right-hand side according
1800  * to sign of the coefficient, i.e., rhs -= coef * bound, where bound = lb if coef >= 0 and bound = ub, otherwise.
1801  *
1802  * @note: The choice of global or local bounds depend on the validity (global or local) of the aggregation row.
1803  *
1804  * @note: The list of non-zero indices will be updated by swapping the last non-zero index to @p pos.
1805  */
1807  SCIP* scip, /**< SCIP data structure */
1808  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
1809  SCIP_VAR* var, /**< variable that should be removed */
1810  int pos, /**< position of the variable in the aggregation row */
1811  SCIP_Bool* valid /**< pointer to return whether the aggregation row is still valid */
1812  )
1813 {
1814  SCIP_Real QUAD(val);
1815  int v;
1816 
1817  assert(valid != NULL);
1818  assert(pos >= 0);
1819 
1820  v = aggrrow->inds[pos];
1821  assert(v == SCIPvarGetProbindex(var));
1822 
1823  QUAD_ARRAY_LOAD(val, aggrrow->vals, v);
1824 
1825  *valid = TRUE;
1826 
1827  /* adjust left and right hand sides with max contribution */
1828  if( QUAD_TO_DBL(val) < 0.0 )
1829  {
1830  SCIP_Real ub = aggrrow->local ? SCIPvarGetUbLocal(var) : SCIPvarGetUbGlobal(var);
1831  if( SCIPisInfinity(scip, ub) )
1832  QUAD_ASSIGN(aggrrow->rhs, SCIPinfinity(scip));
1833  else
1834  {
1835  SCIPquadprecProdQD(val, val, ub);
1836  SCIPquadprecSumQQ(aggrrow->rhs, aggrrow->rhs, -val);
1837  }
1838  }
1839  else
1840  {
1841  SCIP_Real lb = aggrrow->local ? SCIPvarGetLbLocal(var) : SCIPvarGetLbGlobal(var);
1842  if( SCIPisInfinity(scip, -lb) )
1843  QUAD_ASSIGN(aggrrow->rhs, SCIPinfinity(scip));
1844  else
1845  {
1846  SCIPquadprecProdQD(val, val, lb);
1847  SCIPquadprecSumQQ(aggrrow->rhs, aggrrow->rhs, -val);
1848  }
1849  }
1850 
1851  QUAD_ASSIGN(val, 0.0);
1852  QUAD_ARRAY_STORE(aggrrow->vals, v, val);
1853 
1854  /* remove non-zero entry */
1855  --(aggrrow->nnz);
1856  aggrrow->inds[pos] = aggrrow->inds[aggrrow->nnz];
1857 
1858  if( SCIPisInfinity(scip, QUAD_HI(aggrrow->rhs)) )
1859  *valid = FALSE;
1860 }
1861 
1862 /** add the objective function with right-hand side @p rhs and scaled by @p scale to the aggregation row */
1864  SCIP* scip, /**< SCIP data structure */
1865  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
1866  SCIP_Real rhs, /**< right-hand side of the artificial row */
1867  SCIP_Real scale /**< scalar */
1868  )
1869 {
1870  SCIP_VAR** vars;
1871  SCIP_Real QUAD(val);
1872  int nvars;
1873 
1874  assert(scip != NULL);
1875  assert(aggrrow != NULL);
1876 
1877  vars = SCIPgetVars(scip);
1878  nvars = SCIPgetNVars(scip);
1879 
1880  /* add all variables straight forward if the aggregation row is empty */
1881  if( aggrrow->nnz == 0 )
1882  {
1883  int i;
1884  for( i = 0; i < nvars; ++i )
1885  {
1886  assert(SCIPvarGetProbindex(vars[i]) == i);
1887 
1888  /* skip all variables with zero objective coefficient */
1889  if( SCIPisZero(scip, scale * SCIPvarGetObj(vars[i])) )
1890  continue;
1891 
1892  QUAD_ASSIGN(val, scale * SCIPvarGetObj(vars[i]));
1893  QUAD_ARRAY_STORE(aggrrow->vals, i, val);
1894  aggrrow->inds[aggrrow->nnz++] = i;
1895  }
1896 
1897  /* add right-hand side value */
1898  QUAD_ASSIGN(aggrrow->rhs, scale * rhs);
1899  }
1900  else
1901  {
1902  int i;
1903  /* add the non-zeros to the aggregation row and keep non-zero index up to date */
1904  for( i = 0 ; i < nvars; ++i )
1905  {
1906  assert(SCIPvarGetProbindex(vars[i]) == i);
1907 
1908  /* skip all variables with zero objective coefficient */
1909  if( SCIPisZero(scip, scale * SCIPvarGetObj(vars[i])) )
1910  continue;
1911 
1912  QUAD_ARRAY_LOAD(val, aggrrow->vals, i); /* val = aggrrow->vals[i] */
1913 
1914  if( QUAD_HI(val) == 0.0 )
1915  aggrrow->inds[aggrrow->nnz++] = i;
1916 
1917  SCIPquadprecSumQD(val, val, scale * SCIPvarGetObj(vars[i]));
1918 
1919  /* the value must not be exactly zero due to sparsity pattern */
1920  QUAD_HI(val) = NONZERO(QUAD_HI(val));
1921  assert(QUAD_HI(val) != 0.0);
1922 
1923  QUAD_ARRAY_STORE(aggrrow->vals, i, val);
1924  }
1925 
1926  /* add right-hand side value */
1927  SCIPquadprecSumQD(aggrrow->rhs, aggrrow->rhs, scale * rhs);
1928  }
1929 
1930  return SCIP_OKAY;
1931 }
1932 
1933 /** add weighted constraint to the aggregation row */
1935  SCIP* scip, /**< SCIP data structure */
1936  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
1937  int* inds, /**< variable problem indices in constraint to add to the aggregation row */
1938  SCIP_Real* vals, /**< values of constraint to add to the aggregation row */
1939  int len, /**< length of constraint to add to the aggregation row */
1940  SCIP_Real rhs, /**< right hand side of constraint to add to the aggregation row */
1941  SCIP_Real weight, /**< (positive) scale for adding given constraint to the aggregation row */
1942  int rank, /**< rank to use for given constraint */
1943  SCIP_Bool local /**< is constraint only valid locally */
1944  )
1945 {
1946  int i;
1947 
1948  assert(weight >= 0.0);
1949  assert(!SCIPisInfinity(scip, REALABS(weight * rhs)));
1950 
1951  /* update local flag */
1952  aggrrow->local = aggrrow->local || local;
1953 
1954  /* update rank */
1955  aggrrow->rank = MAX(rank, aggrrow->rank);
1956 
1957  /* add right hand side value */
1958  SCIPquadprecSumQD(aggrrow->rhs, aggrrow->rhs, weight * rhs);
1959 
1960  /* add the non-zeros to the aggregation row and keep non-zero index up to date */
1961  for( i = 0 ; i < len; ++i )
1962  {
1963  SCIP_Real QUAD(val);
1964  int probindex = inds[i];
1965 
1966  QUAD_ARRAY_LOAD(val, aggrrow->vals, probindex); /* val = aggrrow->vals[probindex] */
1967 
1968  if( QUAD_HI(val) == 0.0 )
1969  aggrrow->inds[aggrrow->nnz++] = probindex;
1970 
1971  SCIPquadprecSumQD(val, val, vals[i] * weight);
1972 
1973  /* the value must not be exactly zero due to sparsity pattern */
1974  QUAD_HI(val) = NONZERO(QUAD_HI(val));
1975  assert(QUAD_HI(val) != 0.0);
1976 
1977  QUAD_ARRAY_STORE(aggrrow->vals, probindex, val);
1978  }
1979 
1980  return SCIP_OKAY;
1981 }
1982 
1983 /** clear all entries int the aggregation row but don't free memory */
1985  SCIP_AGGRROW* aggrrow /**< the aggregation row */
1986  )
1987 {
1988  int i;
1989  SCIP_Real QUAD(tmp);
1990 
1991  QUAD_ASSIGN(tmp, 0.0);
1992 
1993  for( i = 0; i < aggrrow->nnz; ++i )
1994  {
1995  QUAD_ARRAY_STORE(aggrrow->vals, aggrrow->inds[i], tmp);
1996  }
1997 
1998  aggrrow->nnz = 0;
1999  aggrrow->nrows = 0;
2000  aggrrow->rank = 0;
2001  QUAD_ASSIGN(aggrrow->rhs, 0.0);
2002  aggrrow->local = FALSE;
2003 }
2004 
2005 /** calculates the efficacy norm of the given aggregation row, which depends on the "separating/efficacynorm" parameter
2006  *
2007  * @return the efficacy norm of the given aggregation row, which depends on the "separating/efficacynorm" parameter
2008  */
2010  SCIP* scip, /**< SCIP data structure */
2011  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2012  )
2013 {
2014  return calcEfficacyNormQuad(scip, aggrrow->vals, aggrrow->inds, aggrrow->nnz);
2015 }
2016 
2017 /** Adds one row to the aggregation row. Differs from SCIPaggrRowAddRow() by providing some additional
2018  * parameters required for SCIPaggrRowSumRows()
2019  */
2020 static
2022  SCIP* scip, /**< SCIP data structure */
2023  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2024  SCIP_ROW* row, /**< the row to add */
2025  SCIP_Real weight, /**< weight of row to add */
2026  SCIP_Bool sidetypebasis, /**< choose sidetypes of row (lhs/rhs) based on basis information? */
2027  SCIP_Bool allowlocal, /**< should local rows allowed to be used? */
2028  int negslack, /**< should negative slack variables allowed to be used? (0: no, 1: only for integral rows, 2: yes) */
2029  int maxaggrlen, /**< maximal length of aggregation row */
2030  SCIP_Bool* rowtoolong /**< is the aggregated row too long */
2031  )
2032 {
2033  SCIP_Real sideval;
2034  SCIP_Bool uselhs;
2035  int i;
2036 
2037  assert( rowtoolong != NULL );
2038  *rowtoolong = FALSE;
2039 
2040  if( SCIPisFeasZero(scip, weight) || SCIProwIsModifiable(row) || (SCIProwIsLocal(row) && !allowlocal) )
2041  {
2042  return SCIP_OKAY;
2043  }
2044 
2045  if( sidetypebasis && !SCIPisEQ(scip, SCIProwGetLhs(row), SCIProwGetRhs(row)) )
2046  {
2048 
2049  if( stat == SCIP_BASESTAT_LOWER )
2050  {
2051  assert( ! SCIPisInfinity(scip, -SCIProwGetLhs(row)) );
2052  uselhs = TRUE;
2053  }
2054  else if( stat == SCIP_BASESTAT_UPPER )
2055  {
2056  assert( ! SCIPisInfinity(scip, SCIProwGetRhs(row)) );
2057  uselhs = FALSE;
2058  }
2059  else if( SCIPisInfinity(scip, SCIProwGetRhs(row)) || (weight < 0.0 && ! SCIPisInfinity(scip, -SCIProwGetLhs(row))) )
2060  uselhs = TRUE;
2061  else
2062  uselhs = FALSE;
2063  }
2064  else if( (weight < 0.0 && !SCIPisInfinity(scip, -row->lhs)) || SCIPisInfinity(scip, row->rhs) )
2065  uselhs = TRUE;
2066  else
2067  uselhs = FALSE;
2068 
2069  if( uselhs )
2070  {
2071  assert( ! SCIPisInfinity(scip, -SCIProwGetLhs(row)) );
2072 
2073  if( weight > 0.0 && ((negslack == 0) || (negslack == 1 && !row->integral)) )
2074  return SCIP_OKAY;
2075 
2076  sideval = row->lhs - row->constant;
2077  /* row is integral? round left hand side up */
2078  if( row->integral )
2079  sideval = SCIPceil(scip, sideval);
2080  }
2081  else
2082  {
2083  assert( ! SCIPisInfinity(scip, SCIProwGetRhs(row)) );
2084 
2085  if( weight < 0.0 && ((negslack == 0) || (negslack == 1 && !row->integral)) )
2086  return SCIP_OKAY;
2087 
2088  sideval = row->rhs - row->constant;
2089  /* row is integral? round right hand side down */
2090  if( row->integral )
2091  sideval = SCIPfloor(scip, sideval);
2092  }
2093 
2094  /* add right hand side, update rank and local flag */
2095  SCIPquadprecSumQD(aggrrow->rhs, aggrrow->rhs, sideval * weight);
2096  aggrrow->rank = MAX(aggrrow->rank, row->rank);
2097  aggrrow->local = aggrrow->local || row->local;
2098 
2099  /* ensure the array for storing the row information is large enough */
2100  i = aggrrow->nrows++;
2101  if( aggrrow->nrows > aggrrow->rowssize )
2102  {
2103  int newsize = SCIPcalcMemGrowSize(scip, aggrrow->nrows);
2104  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->rowsinds, aggrrow->rowssize, newsize) );
2105  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->slacksign, aggrrow->rowssize, newsize) );
2106  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &aggrrow->rowweights, aggrrow->rowssize, newsize) );
2107  aggrrow->rowssize = newsize;
2108  }
2109 
2110  /* add information of addditional row */
2111  aggrrow->rowsinds[i] = row->lppos;
2112  aggrrow->rowweights[i] = weight;
2113  aggrrow->slacksign[i] = uselhs ? -1 : 1;
2114 
2115  /* add up coefficients */
2116  SCIP_CALL( varVecAddScaledRowCoefsQuad(aggrrow->inds, aggrrow->vals, &aggrrow->nnz, row, weight) );
2117 
2118  /* check if row is too long now */
2119  if( aggrrow->nnz > maxaggrlen )
2120  *rowtoolong = TRUE;
2121 
2122  return SCIP_OKAY;
2123 }
2124 
2125 /** aggregate rows using the given weights; the current content of the aggregation
2126  * row, \p aggrrow, gets overwritten
2127  */
2129  SCIP* scip, /**< SCIP data structure */
2130  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2131  SCIP_Real* weights, /**< row weights in row summation */
2132  int* rowinds, /**< array to store indices of non-zero entries of the weights array, or NULL */
2133  int nrowinds, /**< number of non-zero entries in weights array, -1 if rowinds is NULL */
2134  SCIP_Bool sidetypebasis, /**< choose sidetypes of row (lhs/rhs) based on basis information? */
2135  SCIP_Bool allowlocal, /**< should local rows allowed to be used? */
2136  int negslack, /**< should negative slack variables allowed to be used? (0: no, 1: only for integral rows, 2: yes) */
2137  int maxaggrlen, /**< maximal length of aggregation row */
2138  SCIP_Bool* valid /**< is the aggregation valid */
2139  )
2140 {
2141  SCIP_ROW** rows;
2142  SCIP_VAR** vars;
2143  int nrows;
2144  int nvars;
2145  int k;
2146  SCIP_Bool rowtoolong;
2147 
2148  assert( scip != NULL );
2149  assert( aggrrow != NULL );
2150  assert( valid != NULL );
2151 
2152  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
2153  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2154 
2155  SCIPaggrRowClear(aggrrow);
2156  *valid = FALSE;
2157 
2158  if( rowinds != NULL && nrowinds > -1 )
2159  {
2160  for( k = 0; k < nrowinds; ++k )
2161  {
2162  SCIP_CALL( addOneRow(scip, aggrrow, rows[rowinds[k]], weights[rowinds[k]], sidetypebasis, allowlocal, negslack, maxaggrlen, &rowtoolong) );
2163 
2164  if( rowtoolong )
2165  return SCIP_OKAY;
2166  }
2167  }
2168  else
2169  {
2170  for( k = 0; k < nrows; ++k )
2171  {
2172  if( weights[k] != 0.0 )
2173  {
2174  SCIP_CALL( addOneRow(scip, aggrrow, rows[k], weights[k], sidetypebasis, allowlocal, negslack, maxaggrlen, &rowtoolong) );
2175 
2176  if( rowtoolong )
2177  return SCIP_OKAY;
2178  }
2179  }
2180  }
2181 
2182  SCIPaggrRowRemoveZeros(scip, aggrrow, FALSE, valid);
2183 
2184  return SCIP_OKAY;
2185 }
2186 
2187 /** checks for cut redundancy and performs activity based coefficient tightening;
2188  * removes coefficients that are zero with QUAD_EPSILON tolerance and uses variable bounds
2189  * to remove small coefficients (relative to the maximum absolute coefficient)
2190  */
2191 static
2193  SCIP* scip, /**< SCIP data structure */
2194  SCIP_Bool cutislocal, /**< is the cut a local cut */
2195  int* cutinds, /**< variable problem indices of non-zeros in cut */
2196  SCIP_Real* cutcoefs, /**< non-zeros coefficients of cut */
2197  int* nnz, /**< number non-zeros coefficients of cut */
2198  SCIP_Real* cutrhs, /**< right hand side of cut */
2199  SCIP_Bool* success /**< pointer to return whether post-processing was succesful or cut is redundant */
2200  )
2201 {
2202  int i;
2203  SCIP_Bool redundant;
2204  SCIP_Real maxcoef;
2205  SCIP_Real minallowedcoef;
2206  SCIP_Real QUAD(rhs);
2207 
2208  assert(scip != NULL);
2209  assert(cutinds != NULL);
2210  assert(cutcoefs != NULL);
2211  assert(cutrhs != NULL);
2212  assert(success != NULL);
2213 
2214  *success = FALSE;
2215 
2216  QUAD_ASSIGN(rhs, *cutrhs);
2217 
2218  if( removeZeros(scip, SCIPfeastol(scip), cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz) )
2219  {
2220  /* right hand side was changed to infinity -> cut is redundant */
2221  return SCIP_OKAY;
2222  }
2223 
2224  if( *nnz == 0 )
2225  return SCIP_OKAY;
2226 
2227  SCIP_CALL( cutTightenCoefs(scip, cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz, &redundant) );
2228 
2229  if( redundant )
2230  {
2231  /* cut is redundant */
2232  return SCIP_OKAY;
2233  }
2234 
2235  maxcoef = 0.0;
2236  for( i = 0; i < *nnz; ++i )
2237  {
2238  SCIP_Real absval = REALABS(cutcoefs[cutinds[i]]);
2239  maxcoef = MAX(absval, maxcoef);
2240  }
2241 
2242  maxcoef /= scip->set->sepa_maxcoefratio;
2243  minallowedcoef = SCIPsumepsilon(scip);
2244  minallowedcoef = MAX(minallowedcoef, maxcoef);
2245 
2246  *success = ! removeZeros(scip, minallowedcoef, cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz);
2247  *cutrhs = QUAD_TO_DBL(rhs);
2248 
2249  return SCIP_OKAY;
2250 }
2251 
2252 
2253 /** checks for cut redundancy and performs activity based coefficient tightening;
2254  * removes coefficients that are zero with QUAD_EPSILON tolerance and uses variable bounds
2255  * to remove small coefficients (relative to the maximum absolute coefficient).
2256  * The cutcoefs must be a quad precision array, i.e. allocated with size
2257  * QUAD_ARRAY_SIZE(nvars) and accessed with QUAD_ARRAY_LOAD and QUAD_ARRAY_STORE
2258  * macros.
2259  */
2260 static
2262  SCIP* scip, /**< SCIP data structure */
2263  SCIP_Bool cutislocal, /**< is the cut a local cut */
2264  int* cutinds, /**< variable problem indices of non-zeros in cut */
2265  SCIP_Real* cutcoefs, /**< non-zeros coefficients of cut */
2266  int* nnz, /**< number non-zeros coefficients of cut */
2267  QUAD(SCIP_Real* cutrhs), /**< right hand side of cut */
2268  SCIP_Bool* success /**< pointer to return whether the cleanup was successful or if it is useless */
2269  )
2270 {
2271  int i;
2272  SCIP_Bool redundant;
2273  SCIP_Real maxcoef;
2274  SCIP_Real minallowedcoef;
2275 
2276  assert(scip != NULL);
2277  assert(cutinds != NULL);
2278  assert(cutcoefs != NULL);
2279  assert(QUAD_HI(cutrhs) != NULL);
2280  assert(success != NULL);
2281 
2282  *success = FALSE;
2283 
2284  if( removeZerosQuad(scip, SCIPfeastol(scip), cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz) )
2285  {
2286  /* right hand side was changed to infinity -> cut is redundant */
2287  return SCIP_OKAY;
2288  }
2289 
2290  if( *nnz == 0 )
2291  return SCIP_OKAY;
2292 
2293  SCIP_CALL( cutTightenCoefsQuad(scip, cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz, &redundant) );
2294  if( redundant )
2295  {
2296  /* cut is redundant */
2297  return SCIP_OKAY;
2298  }
2299 
2300  maxcoef = 0.0;
2301  for( i = 0; i < *nnz; ++i )
2302  {
2303  SCIP_Real abscoef;
2304  SCIP_Real QUAD(coef);
2305  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[i]); /* coef = cutcoefs[cutinds[i]] */
2306  abscoef = REALABS(QUAD_TO_DBL(coef));
2307  maxcoef = MAX(abscoef, maxcoef);
2308  }
2309 
2310  maxcoef /= scip->set->sepa_maxcoefratio;
2311  minallowedcoef = SCIPsumepsilon(scip);
2312  minallowedcoef = MAX(minallowedcoef, maxcoef);
2313 
2314  *success = ! removeZerosQuad(scip, minallowedcoef, cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz);
2315 
2316  return SCIP_OKAY;
2317 }
2318 
2319 /** removes almost zero entries from the aggregation row. */
2321  SCIP* scip, /**< SCIP datastructure */
2322  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2323  SCIP_Bool useglbbounds, /**< consider global bound although the cut is local? */
2324  SCIP_Bool* valid /**< pointer to return whether the aggregation row is still valid */
2325  )
2326 {
2327  assert(aggrrow != NULL);
2328  assert(valid != NULL);
2329 
2330  *valid = ! removeZerosQuad(scip, SCIPsumepsilon(scip), useglbbounds ? FALSE : aggrrow->local, aggrrow->vals,
2331  QUAD(&aggrrow->rhs), aggrrow->inds, &aggrrow->nnz);
2332 }
2333 
2334 /** get number of aggregated rows */
2336  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2337  )
2338 {
2339  assert(aggrrow != NULL);
2340 
2341  return aggrrow->nrows;
2342 }
2343 
2344 /** get array with lp positions of rows used in aggregation */
2346  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2347  )
2348 {
2349  assert(aggrrow != NULL);
2350  assert(aggrrow->rowsinds != NULL || aggrrow->nrows == 0);
2351 
2352  return aggrrow->rowsinds;
2353 }
2354 
2355 /** get array with weights of aggregated rows */
2357  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2358  )
2359 {
2360  assert(aggrrow != NULL);
2361  assert(aggrrow->rowweights != NULL || aggrrow->nrows == 0);
2362 
2363  return aggrrow->rowweights;
2364 }
2365 
2366 /** checks whether a given row has been added to the aggregation row */
2368  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2369  SCIP_ROW* row /**< row for which it is checked whether it has been added to the aggregation */
2370  )
2371 {
2372  int i;
2373  int rowind;
2374 
2375  assert(aggrrow != NULL);
2376  assert(row != NULL);
2377 
2378  rowind = SCIProwGetLPPos(row);
2379 
2380  for( i = 0; i < aggrrow->nrows; ++i )
2381  {
2382  if( aggrrow->rowsinds[i] == rowind )
2383  return TRUE;
2384  }
2385 
2386  return FALSE;
2387 }
2388 
2389 /** gets the array of corresponding variable problem indices for each non-zero in the aggregation row */
2391  SCIP_AGGRROW* aggrrow /**< aggregation row */
2392  )
2393 {
2394  assert(aggrrow != NULL);
2395 
2396  return aggrrow->inds;
2397 }
2398 
2399 /** gets the number of non-zeros in the aggregation row */
2401  SCIP_AGGRROW* aggrrow /**< aggregation row */
2402  )
2403 {
2404  assert(aggrrow != NULL);
2405 
2406  return aggrrow->nnz;
2407 }
2408 
2409 /** gets the rank of the aggregation row */
2411  SCIP_AGGRROW* aggrrow /**< aggregation row */
2412  )
2413 {
2414  assert(aggrrow != NULL);
2415 
2416  return aggrrow->rank;
2417 }
2418 
2419 /** checks if the aggregation row is only valid locally */
2421  SCIP_AGGRROW* aggrrow /**< aggregation row */
2422  )
2423 {
2424  assert(aggrrow != NULL);
2425 
2426  return aggrrow->local;
2427 }
2428 
2429 /** gets the right hand side of the aggregation row */
2431  SCIP_AGGRROW* aggrrow /**< aggregation row */
2432  )
2433 {
2434  assert(aggrrow != NULL);
2435 
2436  return QUAD_TO_DBL(aggrrow->rhs);
2437 }
2438 
2439 /** filters the given array of cuts to enforce a maximum parallelism constraints
2440  * for the given cut; moves filtered cuts to the end of the array and returns number of selected cuts */
2441 static
2443  SCIP_ROW* cut, /**< cut to filter orthogonality with */
2444  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2445  SCIP_Real* scores, /**< array with scores of cuts to perform selection algorithm */
2446  int ncuts, /**< number of cuts in given array */
2447  SCIP_Real goodscore, /**< threshold for the score to be considered a good cut */
2448  SCIP_Real goodmaxparall, /**< maximal parallelism for good cuts */
2449  SCIP_Real maxparall /**< maximal parallelism for all cuts that are not good */
2450  )
2451 {
2452  int i;
2453 
2454  assert( cut != NULL );
2455  assert( ncuts == 0 || cuts != NULL );
2456  assert( ncuts == 0 || scores != NULL );
2457 
2458  for( i = ncuts - 1; i >= 0; --i )
2459  {
2460  SCIP_Real thisparall;
2461  SCIP_Real thismaxparall;
2462 
2463  thisparall = SCIProwGetParallelism(cut, cuts[i], 'e');
2464  thismaxparall = scores[i] >= goodscore ? goodmaxparall : maxparall;
2465 
2466  if( thisparall > thismaxparall )
2467  {
2468  --ncuts;
2469  SCIPswapPointers((void**) &cuts[i], (void**) &cuts[ncuts]);
2470  SCIPswapReals(&scores[i], &scores[ncuts]);
2471  }
2472  }
2473 
2474  return ncuts;
2475 }
2476 
2477 /** move the cut with the highest score to the first position in the array; there must be at least one cut */
2478 static
2480  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2481  SCIP_Real* scores, /**< array with scores of cuts to perform selection algorithm */
2482  int ncuts /**< number of cuts in given array */
2483  )
2484 {
2485  int i;
2486  int bestpos;
2487  SCIP_Real bestscore;
2488 
2489  assert(ncuts > 0);
2490  assert(cuts != NULL);
2491  assert(scores != NULL);
2492 
2493  bestscore = scores[0];
2494  bestpos = 0;
2495 
2496  for( i = 1; i < ncuts; ++i )
2497  {
2498  if( scores[i] > bestscore )
2499  {
2500  bestpos = i;
2501  bestscore = scores[i];
2502  }
2503  }
2504 
2505  SCIPswapPointers((void**) &cuts[bestpos], (void**) &cuts[0]);
2506  SCIPswapReals(&scores[bestpos], &scores[0]);
2507 }
2508 
2509 /** perform a cut selection algorithm for the given array of cuts; the array is partitioned
2510  * so that the selected cuts come first and the remaining ones are at the end of the array
2511  */
2513  SCIP* scip, /**< SCIP data structure */
2514  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2515  SCIP_RANDNUMGEN* randnumgen, /**< random number generator for tie-breaking, or NULL */
2516  SCIP_Real goodscorefac, /**< factor of best score among the given cuts to consider a cut good
2517  * and filter with less strict settings of the maximum parallelism */
2518  SCIP_Real badscorefac, /**< factor of best score among the given cuts to consider a cut bad
2519  * and discard it regardless of its parallelism to other cuts */
2520  SCIP_Real goodmaxparall, /**< maximum parallelism for good cuts */
2521  SCIP_Real maxparall, /**< maximum parallelism for non-good cuts */
2522  SCIP_Real dircutoffdistweight,/**< weight of directed cutoff distance in score calculation */
2523  SCIP_Real efficacyweight, /**< weight of efficacy (shortest cutoff distance) in score calculation */
2524  SCIP_Real objparalweight, /**< weight of objective parallelism in score calculation */
2525  SCIP_Real intsupportweight, /**< weight of integral support in score calculation */
2526  int ncuts, /**< number of cuts in given array */
2527  int nforcedcuts, /**< number of forced cuts at start of given array */
2528  int maxselectedcuts, /**< maximal number of cuts to select */
2529  int* nselectedcuts /**< pointer to return number of selected cuts */
2530  )
2531 {
2532  int i;
2533  SCIP_Real* scores;
2534  SCIP_Real goodscore;
2535  SCIP_Real badscore;
2536  SCIP_Real efficacyfac;
2537  SCIP_SOL* sol;
2538 
2539  assert( nselectedcuts != NULL );
2540 
2541  /* if all cuts are forced cuts, no selection is required */
2542  if( nforcedcuts >= MIN(ncuts, maxselectedcuts) )
2543  {
2544  *nselectedcuts = nforcedcuts;
2545  return SCIP_OKAY;
2546  }
2547  *nselectedcuts = 0;
2548 
2549  SCIP_CALL( SCIPallocBufferArray(scip, &scores, ncuts) );
2550 
2551  sol = SCIPgetBestSol(scip);
2552 
2553  efficacyfac = efficacyweight;
2554 
2555  goodscore = 0.0;
2556 
2557  /* if there is an incumbent and the factor is not 0.0, compute directed cutoff distances for the incumbent */
2558  if( sol != NULL && dircutoffdistweight > 0.0 )
2559  {
2560  for( i = 0; i < ncuts; ++i )
2561  {
2562  SCIP_Real objparallelism;
2563  SCIP_Real intsupport;
2564  SCIP_Real efficacy;
2565 
2566  intsupport = intsupportweight != 0.0 ?
2567  intsupportweight * SCIProwGetNumIntCols(cuts[i], scip->set) / (SCIP_Real) SCIProwGetNNonz(cuts[i]) :
2568  0.0;
2569 
2570  objparallelism = objparalweight != 0.0 ? objparalweight * SCIProwGetObjParallelism(cuts[i], scip->set, scip->lp) : 0.0;
2571 
2572  efficacy = SCIProwGetLPEfficacy(cuts[i], scip->set, scip->stat, scip->lp);
2573 
2574  if( SCIProwIsLocal(cuts[i]) )
2575  {
2576  scores[i] = dircutoffdistweight * efficacy;
2577  }
2578  else
2579  {
2580  scores[i] = SCIProwGetLPSolCutoffDistance(cuts[i], scip->set, scip->stat, sol, scip->lp);
2581  scores[i] = dircutoffdistweight * MAX(scores[i], efficacy);
2582  }
2583 
2584  efficacy *= efficacyfac;
2585  scores[i] += objparallelism + intsupport + efficacy;
2586 
2587  /* add small term to prefer global pool cuts */
2588  scores[i] += SCIProwIsInGlobalCutpool(cuts[i]) ? 1e-4 : 0.0;
2589 
2590  if( randnumgen != NULL )
2591  {
2592  scores[i] += SCIPrandomGetReal(randnumgen, 0.0, 1e-6);
2593  }
2594 
2595  goodscore = MAX(goodscore, scores[i]);
2596  }
2597  }
2598  else
2599  {
2600  /* in case there is no solution add the directed cutoff distance weight to the efficacy weight
2601  * since the efficacy underestimates the directed cuttoff distance
2602  */
2603  efficacyfac += dircutoffdistweight;
2604  for( i = 0; i < ncuts; ++i )
2605  {
2606  SCIP_Real objparallelism;
2607  SCIP_Real intsupport;
2608  SCIP_Real efficacy;
2609 
2610  intsupport = intsupportweight > 0.0 ?
2611  intsupportweight * SCIProwGetNumIntCols(cuts[i], scip->set) / (SCIP_Real) SCIProwGetNNonz(cuts[i]) :
2612  0.0;
2613 
2614  objparallelism = objparalweight > 0.0 ? objparalweight * SCIProwGetObjParallelism(cuts[i], scip->set, scip->lp) : 0.0;
2615 
2616  efficacy = efficacyfac > 0.0 ?
2617  efficacyfac * SCIProwGetLPEfficacy(cuts[i], scip->set, scip->stat, scip->lp) :
2618  0.0;
2619 
2620  scores[i] = objparallelism + intsupport + efficacy;
2621 
2622  /* add small term to prefer global pool cuts */
2623  scores[i] += SCIProwIsInGlobalCutpool(cuts[i]) ? 1e-4 : 0.0;
2624 
2625  if( randnumgen != NULL )
2626  {
2627  scores[i] += SCIPrandomGetReal(randnumgen, 0.0, 1e-6);
2628  }
2629 
2630  goodscore = MAX(goodscore, scores[i]);
2631  }
2632  }
2633 
2634  /* compute values for filtering cuts */
2635  badscore = goodscore * badscorefac;
2636  goodscore *= goodscorefac;
2637 
2638  /* perform cut selection algorithm for the cuts */
2639  {
2640  int nnonforcedcuts;
2641  SCIP_ROW** nonforcedcuts;
2642  SCIP_Real* nonforcedscores;
2643 
2644  /* adjust pointers to the beginning of the non-forced cuts */
2645  nnonforcedcuts = ncuts - nforcedcuts;
2646  nonforcedcuts = cuts + nforcedcuts;
2647  nonforcedscores = scores + nforcedcuts;
2648 
2649  /* select the forced cuts first */
2650  *nselectedcuts = nforcedcuts;
2651  for( i = 0; i < nforcedcuts && nnonforcedcuts > 0; ++i )
2652  {
2653  nnonforcedcuts = filterWithParallelism(cuts[i], nonforcedcuts, nonforcedscores, nnonforcedcuts, goodscore, goodmaxparall, maxparall);
2654  }
2655 
2656  /* if the maximal number of cuts was exceeded after selecting the forced cuts, we can stop here */
2657  if( *nselectedcuts >= maxselectedcuts )
2658  goto TERMINATE;
2659 
2660  /* now greedily select the remaining cuts */
2661  while( nnonforcedcuts > 0 )
2662  {
2663  SCIP_ROW* selectedcut;
2664 
2665  selectBestCut(nonforcedcuts, nonforcedscores, nnonforcedcuts);
2666  selectedcut = nonforcedcuts[0];
2667 
2668  /* if the best cut of the remaining cuts is considered bad, we discard it and all remaining cuts */
2669  if( nonforcedscores[0] < badscore )
2670  goto TERMINATE;
2671 
2672  ++(*nselectedcuts);
2673 
2674  /* if the maximal number of cuts was selected, we can stop here */
2675  if( *nselectedcuts == maxselectedcuts )
2676  goto TERMINATE;
2677 
2678  /* move the pointers to the next position and filter the remaining cuts to enforce the maximum parallelism constraint */
2679  ++nonforcedcuts;
2680  ++nonforcedscores;
2681  --nnonforcedcuts;
2682 
2683  nnonforcedcuts = filterWithParallelism(selectedcut, nonforcedcuts, nonforcedscores, nnonforcedcuts, goodscore, goodmaxparall, maxparall);
2684  }
2685  }
2686 
2687  TERMINATE:
2688  SCIPfreeBufferArray(scip, &scores);
2689 
2690  return SCIP_OKAY;
2691 }
2692 
2693 /* =========================================== c-MIR =========================================== */
2694 
2695 #define MAXCMIRSCALE 1e+6 /**< maximal scaling (scale/(1-f0)) allowed in c-MIR calculations */
2696 
2697 /** finds the best lower bound of the variable to use for MIR transformation */
2698 static
2700  SCIP* scip, /**< SCIP data structure */
2701  SCIP_VAR* var, /**< problem variable */
2702  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2703  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2704  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2705  SCIP_Real* bestlb, /**< pointer to store best bound value */
2706  SCIP_Real* simplebound, /**< pointer to store simple bound value */
2707  int* bestlbtype /**< pointer to store best bound type */
2708  )
2709 {
2710  assert(bestlb != NULL);
2711  assert(bestlbtype != NULL);
2712 
2713  *bestlb = SCIPvarGetLbGlobal(var);
2714  *bestlbtype = -1;
2715 
2716  if( allowlocal )
2717  {
2718  SCIP_Real loclb;
2719 
2720  loclb = SCIPvarGetLbLocal(var);
2721  if( SCIPisGT(scip, loclb, *bestlb) )
2722  {
2723  *bestlb = loclb;
2724  *bestlbtype = -2;
2725  }
2726  }
2727 
2728  *simplebound = *bestlb;
2729 
2730  if( usevbds && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
2731  {
2732  SCIP_Real bestvlb;
2733  int bestvlbidx;
2734 
2735  SCIP_CALL( SCIPgetVarClosestVlb(scip, var, sol, &bestvlb, &bestvlbidx) );
2736  if( bestvlbidx >= 0
2737  && (bestvlb > *bestlb || (*bestlbtype < 0 && SCIPisGE(scip, bestvlb, *bestlb))) )
2738  {
2739  SCIP_VAR** vlbvars;
2740 
2741  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2742  /**@todo this check is not needed for continuous variables; but allowing all but binary variables
2743  * to be replaced by variable bounds seems to be buggy (wrong result on gesa2)
2744  */
2745  vlbvars = SCIPvarGetVlbVars(var);
2746  assert(vlbvars != NULL);
2747  if( SCIPvarGetProbindex(vlbvars[bestvlbidx]) < SCIPvarGetProbindex(var) )
2748  {
2749  *bestlb = bestvlb;
2750  *bestlbtype = bestvlbidx;
2751  }
2752  }
2753  }
2754 
2755  return SCIP_OKAY;
2756 }
2757 
2758 /** finds the best upper bound of the variable to use for MIR transformation */
2759 static
2761  SCIP* scip, /**< SCIP data structure */
2762  SCIP_VAR* var, /**< problem variable */
2763  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2764  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2765  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2766  SCIP_Real* bestub, /**< pointer to store best bound value */
2767  SCIP_Real* simplebound, /**< pointer to store simple bound */
2768  int* bestubtype /**< pointer to store best bound type */
2769  )
2770 {
2771  assert(bestub != NULL);
2772  assert(bestubtype != NULL);
2773 
2774  *bestub = SCIPvarGetUbGlobal(var);
2775  *bestubtype = -1;
2776 
2777  if( allowlocal )
2778  {
2779  SCIP_Real locub;
2780 
2781  locub = SCIPvarGetUbLocal(var);
2782  if( SCIPisLT(scip, locub, *bestub) )
2783  {
2784  *bestub = locub;
2785  *bestubtype = -2;
2786  }
2787  }
2788 
2789  *simplebound = *bestub;
2790 
2791  if( usevbds && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
2792  {
2793  SCIP_Real bestvub;
2794  int bestvubidx;
2795 
2796  SCIP_CALL( SCIPgetVarClosestVub(scip, var, sol, &bestvub, &bestvubidx) );
2797  if( bestvubidx >= 0
2798  && (bestvub < *bestub || (*bestubtype < 0 && SCIPisLE(scip, bestvub, *bestub))) )
2799  {
2800  SCIP_VAR** vubvars;
2801 
2802  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2803  /**@todo this check is not needed for continuous variables; but allowing all but binary variables
2804  * to be replaced by variable bounds seems to be buggy (wrong result on gesa2)
2805  */
2806  vubvars = SCIPvarGetVubVars(var);
2807  assert(vubvars != NULL);
2808  if( SCIPvarGetProbindex(vubvars[bestvubidx]) < SCIPvarGetProbindex(var) )
2809  {
2810  *bestub = bestvub;
2811  *bestubtype = bestvubidx;
2812  }
2813  }
2814  }
2815 
2816  return SCIP_OKAY;
2817 }
2818 
2819 /** determine the best bounds with respect to the given solution for complementing the given variable */
2820 static
2822  SCIP* scip, /**< SCIP data structure */
2823  SCIP_VAR* var, /**< variable to determine best bound for */
2824  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2825  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
2826  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2827  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2828  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
2829  SCIP_Bool ignoresol, /**< should the LP solution be ignored? (eg, apply MIR to dualray) */
2830  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
2831  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
2832  * NULL for using closest bound for all variables */
2833  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
2834  * NULL for using closest bound for all variables */
2835  SCIP_Real* bestlb, /**< pointer to store best lower bound of variable */
2836  SCIP_Real* bestub, /**< pointer to store best upper bound of variable */
2837  int* bestlbtype, /**< pointer to store type of best lower bound of variable */
2838  int* bestubtype, /**< pointer to store type of best upper bound of variable */
2839  SCIP_BOUNDTYPE* selectedbound, /**< pointer to store whether the lower bound or the upper bound should be preferred */
2840  SCIP_Bool* freevariable /**< pointer to store if this is a free variable */
2841  )
2842 {
2843  SCIP_Real simplelb;
2844  SCIP_Real simpleub;
2845  int v;
2846 
2847  v = SCIPvarGetProbindex(var);
2848 
2849  /* check if the user specified a bound to be used */
2850  if( boundsfortrans != NULL && boundsfortrans[v] > -3 )
2851  {
2852  assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || ( boundsfortrans[v] == -2 || boundsfortrans[v] == -1 ));
2853  assert(boundtypesfortrans != NULL);
2854 
2855  /* user has explicitly specified a bound to be used */
2856  if( boundtypesfortrans[v] == SCIP_BOUNDTYPE_LOWER )
2857  {
2858  /* user wants to use lower bound */
2859  *bestlbtype = boundsfortrans[v];
2860  if( *bestlbtype == -1 )
2861  *bestlb = SCIPvarGetLbGlobal(var); /* use global standard lower bound */
2862  else if( *bestlbtype == -2 )
2863  *bestlb = SCIPvarGetLbLocal(var); /* use local standard lower bound */
2864  else
2865  {
2866  SCIP_VAR** vlbvars;
2867  SCIP_Real* vlbcoefs;
2868  SCIP_Real* vlbconsts;
2869  int k;
2870 
2871  assert(!ignoresol);
2872 
2873  /* use the given variable lower bound */
2874  vlbvars = SCIPvarGetVlbVars(var);
2875  vlbcoefs = SCIPvarGetVlbCoefs(var);
2876  vlbconsts = SCIPvarGetVlbConstants(var);
2877  k = boundsfortrans[v];
2878  assert(k >= 0 && k < SCIPvarGetNVlbs(var));
2879  assert(vlbvars != NULL);
2880  assert(vlbcoefs != NULL);
2881  assert(vlbconsts != NULL);
2882 
2883  *bestlb = vlbcoefs[k] * (sol == NULL ? SCIPvarGetLPSol(vlbvars[k]) : SCIPgetSolVal(scip, sol, vlbvars[k])) + vlbconsts[k];
2884  }
2885 
2886  assert(!SCIPisInfinity(scip, - *bestlb));
2887  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2888 
2889  /* find closest upper bound in standard upper bound (and variable upper bounds for continuous variables) */
2890  SCIP_CALL( findBestUb(scip, var, sol, usevbds && fixintegralrhs, allowlocal && fixintegralrhs, bestub, &simpleub, bestubtype) );
2891  }
2892  else
2893  {
2894  assert(boundtypesfortrans[v] == SCIP_BOUNDTYPE_UPPER);
2895 
2896  /* user wants to use upper bound */
2897  *bestubtype = boundsfortrans[v];
2898  if( *bestubtype == -1 )
2899  *bestub = SCIPvarGetUbGlobal(var); /* use global standard upper bound */
2900  else if( *bestubtype == -2 )
2901  *bestub = SCIPvarGetUbLocal(var); /* use local standard upper bound */
2902  else
2903  {
2904  SCIP_VAR** vubvars;
2905  SCIP_Real* vubcoefs;
2906  SCIP_Real* vubconsts;
2907  int k;
2908 
2909  assert(!ignoresol);
2910 
2911  /* use the given variable upper bound */
2912  vubvars = SCIPvarGetVubVars(var);
2913  vubcoefs = SCIPvarGetVubCoefs(var);
2914  vubconsts = SCIPvarGetVubConstants(var);
2915  k = boundsfortrans[v];
2916  assert(k >= 0 && k < SCIPvarGetNVubs(var));
2917  assert(vubvars != NULL);
2918  assert(vubcoefs != NULL);
2919  assert(vubconsts != NULL);
2920 
2921  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2922  *bestub = vubcoefs[k] * (sol == NULL ? SCIPvarGetLPSol(vubvars[k]) : SCIPgetSolVal(scip, sol, vubvars[k])) + vubconsts[k];
2923  }
2924 
2925  assert(!SCIPisInfinity(scip, *bestub));
2926  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2927 
2928  /* find closest lower bound in standard lower bound (and variable lower bounds for continuous variables) */
2929  SCIP_CALL( findBestLb(scip, var, sol, usevbds && fixintegralrhs, allowlocal && fixintegralrhs, bestlb, &simplelb, bestlbtype) );
2930  }
2931  }
2932  else
2933  {
2934  SCIP_Real varsol;
2935 
2936  /* bound selection should be done automatically */
2937 
2938  /* find closest lower bound in standard lower bound (and variable lower bounds for continuous variables) */
2939  SCIP_CALL( findBestLb(scip, var, sol, usevbds, allowlocal, bestlb, &simplelb, bestlbtype) );
2940 
2941  /* find closest upper bound in standard upper bound (and variable upper bounds for continuous variables) */
2942  SCIP_CALL( findBestUb(scip, var, sol, usevbds, allowlocal, bestub, &simpleub, bestubtype) );
2943 
2944  /* check, if variable is free variable */
2945  if( SCIPisInfinity(scip, - *bestlb) && SCIPisInfinity(scip, *bestub) )
2946  {
2947  /* we found a free variable in the row with non-zero coefficient
2948  * -> MIR row can't be transformed in standard form
2949  */
2950  *freevariable = TRUE;
2951  return SCIP_OKAY;
2952  }
2953 
2954  if( !ignoresol )
2955  {
2956  /* select transformation bound */
2957  varsol = (sol == NULL ? SCIPvarGetLPSol(var) : SCIPgetSolVal(scip, sol, var));
2958 
2959  if( SCIPisInfinity(scip, *bestub) ) /* if there is no ub, use lb */
2960  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2961  else if( SCIPisInfinity(scip, - *bestlb) ) /* if there is no lb, use ub */
2962  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2963  else if( SCIPisLT(scip, varsol, (1.0 - boundswitch) * (*bestlb) + boundswitch * (*bestub)) )
2964  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2965  else if( SCIPisGT(scip, varsol, (1.0 - boundswitch) * (*bestlb) + boundswitch * (*bestub)) )
2966  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2967  else if( *bestlbtype == -1 ) /* prefer global standard bounds */
2968  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2969  else if( *bestubtype == -1 ) /* prefer global standard bounds */
2970  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2971  else if( ((*bestlbtype) >= 0 || (*bestubtype) >= 0) && !SCIPisEQ(scip, *bestlb - simplelb, simpleub - *bestub) )
2972  {
2973  if( *bestlb - simplelb > simpleub - *bestub )
2974  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2975  else
2976  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2977  }
2978  else if( *bestlbtype >= 0 ) /* prefer variable bounds over local bounds */
2979  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2980  else if( *bestubtype >= 0 ) /* prefer variable bounds over local bounds */
2981  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2982  else /* no decision yet? just use lower bound */
2983  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2984  }
2985  else
2986  {
2987  SCIP_Real glbub = SCIPvarGetUbGlobal(var);
2988  SCIP_Real glblb = SCIPvarGetLbGlobal(var);
2989  SCIP_Real distlb = REALABS(glblb - *bestlb);
2990  SCIP_Real distub = REALABS(glbub - *bestub);
2991 
2992  assert(!SCIPisInfinity(scip, - *bestlb) || !SCIPisInfinity(scip, *bestub));
2993 
2994  if( SCIPisInfinity(scip, - *bestlb) )
2995  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2996  else if( !SCIPisNegative(scip, *bestlb) )
2997  {
2998  if( SCIPisInfinity(scip, *bestub) )
2999  *selectedbound = SCIP_BOUNDTYPE_LOWER;
3000  else if( SCIPisZero(scip, glblb) )
3001  *selectedbound = SCIP_BOUNDTYPE_LOWER;
3002  else if( SCIPisLE(scip, distlb, distub) )
3003  *selectedbound = SCIP_BOUNDTYPE_LOWER;
3004  else
3005  *selectedbound = SCIP_BOUNDTYPE_UPPER;
3006  }
3007  else
3008  {
3009  assert(!SCIPisInfinity(scip, - *bestlb));
3010  *selectedbound = SCIP_BOUNDTYPE_LOWER;
3011  }
3012  }
3013  }
3014 
3015  return SCIP_OKAY; /*lint !e438*/
3016 }
3017 
3018 /** Transform equation \f$ a \cdot x = b; lb \leq x \leq ub \f$ into standard form
3019  * \f$ a^\prime \cdot x^\prime = b,\; 0 \leq x^\prime \leq ub' \f$.
3020  *
3021  * Transform variables (lb or ub):
3022  * \f[
3023  * \begin{array}{llll}
3024  * x^\prime_j := x_j - lb_j,& x_j = x^\prime_j + lb_j,& a^\prime_j = a_j,& \mbox{if lb is used in transformation}\\
3025  * x^\prime_j := ub_j - x_j,& x_j = ub_j - x^\prime_j,& a^\prime_j = -a_j,& \mbox{if ub is used in transformation}
3026  * \end{array}
3027  * \f]
3028  * and move the constant terms \f$ a_j\, lb_j \f$ or \f$ a_j\, ub_j \f$ to the rhs.
3029  *
3030  * Transform variables (vlb or vub):
3031  * \f[
3032  * \begin{array}{llll}
3033  * x^\prime_j := x_j - (bl_j\, zl_j + dl_j),& x_j = x^\prime_j + (bl_j\, zl_j + dl_j),& a^\prime_j = a_j,& \mbox{if vlb is used in transf.} \\
3034  * x^\prime_j := (bu_j\, zu_j + du_j) - x_j,& x_j = (bu_j\, zu_j + du_j) - x^\prime_j,& a^\prime_j = -a_j,& \mbox{if vub is used in transf.}
3035  * \end{array}
3036  * \f]
3037  * move the constant terms \f$ a_j\, dl_j \f$ or \f$ a_j\, du_j \f$ to the rhs, and update the coefficient of the VLB variable:
3038  * \f[
3039  * \begin{array}{ll}
3040  * a_{zl_j} := a_{zl_j} + a_j\, bl_j,& \mbox{or} \\
3041  * a_{zu_j} := a_{zu_j} + a_j\, bu_j &
3042  * \end{array}
3043  * \f]
3044  */
3045 static
3047  SCIP* scip, /**< SCIP data structure */
3048  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
3049  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
3050  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
3051  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
3052  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
3053  SCIP_Bool ignoresol, /**< should the LP solution be ignored? (eg, apply MIR to dualray) */
3054  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
3055  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
3056  * NULL for using closest bound for all variables */
3057  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
3058  * NULL for using closest bound for all variables */
3059  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
3060  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
3061  SCIP_Real* cutcoefs, /**< array of coefficients of cut */
3062  QUAD(SCIP_Real* cutrhs), /**< pointer to right hand side of cut */
3063  int* cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3064  int* nnz, /**< number of non-zeros in cut */
3065  int* varsign, /**< stores the sign of the transformed variable in summation */
3066  int* boundtype, /**< stores the bound used for transformed variable:
3067  * vlb/vub_idx, or -1 for global lb/ub, or -2 for local lb/ub */
3068  SCIP_Bool* freevariable, /**< stores whether a free variable was found in MIR row -> invalid summation */
3069  SCIP_Bool* localbdsused /**< pointer to store whether local bounds were used in transformation */
3070  )
3071 {
3072  SCIP_Real QUAD(tmp);
3073  SCIP_Real* bestlbs;
3074  SCIP_Real* bestubs;
3075  int* bestlbtypes;
3076  int* bestubtypes;
3077  SCIP_BOUNDTYPE* selectedbounds;
3078  int i;
3079  int aggrrowintstart;
3080  int nvars;
3081  int firstcontvar;
3082  SCIP_VAR** vars;
3083 
3084  assert(varsign != NULL);
3085  assert(boundtype != NULL);
3086  assert(freevariable != NULL);
3087  assert(localbdsused != NULL);
3088 
3089  *freevariable = FALSE;
3090  *localbdsused = FALSE;
3091 
3092  /* allocate temporary memory to store best bounds and bound types */
3093  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbs, 2*(*nnz)) );
3094  SCIP_CALL( SCIPallocBufferArray(scip, &bestubs, 2*(*nnz)) );
3095  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbtypes, 2*(*nnz)) );
3096  SCIP_CALL( SCIPallocBufferArray(scip, &bestubtypes, 2*(*nnz)) );
3097  SCIP_CALL( SCIPallocBufferArray(scip, &selectedbounds, 2*(*nnz)) );
3098 
3099  /* start with continuous variables, because using variable bounds can affect the untransformed integral
3100  * variables, and these changes have to be incorporated in the transformation of the integral variables
3101  * (continuous variables have largest problem indices!)
3102  */
3103  SCIPsortDownInt(cutinds, *nnz);
3104 
3105  vars = SCIPgetVars(scip);
3106  nvars = SCIPgetNVars(scip);
3107  firstcontvar = nvars - SCIPgetNContVars(scip);
3108 
3109  /* determine the best bounds for the continous variables */
3110  for( i = 0; i < *nnz && cutinds[i] >= firstcontvar; ++i )
3111  {
3112  SCIP_CALL( determineBestBounds(scip, vars[cutinds[i]], sol, boundswitch, usevbds, allowlocal, fixintegralrhs,
3113  ignoresol, boundsfortrans, boundtypesfortrans,
3114  bestlbs + i, bestubs + i, bestlbtypes + i, bestubtypes + i, selectedbounds + i, freevariable) );
3115 
3116  if( *freevariable )
3117  goto TERMINATE;
3118  }
3119 
3120  /* remember start of integer variables in the aggrrow */
3121  aggrrowintstart = i;
3122 
3123  /* perform bound substitution for continuous variables */
3124  for( i = 0; i < aggrrowintstart; ++i )
3125  {
3126  SCIP_Real QUAD(coef);
3127  int v = cutinds[i];
3128  SCIP_VAR* var = vars[v];
3129 
3130  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3131 
3132  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
3133  {
3134  assert(!SCIPisInfinity(scip, -bestlbs[i]));
3135 
3136  /* use lower bound as transformation bound: x'_j := x_j - lb_j */
3137  boundtype[i] = bestlbtypes[i];
3138  varsign[i] = +1;
3139 
3140  /* standard (bestlbtype < 0) or variable (bestlbtype >= 0) lower bound? */
3141  if( bestlbtypes[i] < 0 )
3142  {
3143  SCIPquadprecProdQD(tmp, coef, bestlbs[i]);
3144  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3145  *localbdsused = *localbdsused || (bestlbtypes[i] == -2);
3146  }
3147  else
3148  {
3149  SCIP_VAR** vlbvars;
3150  SCIP_Real* vlbcoefs;
3151  SCIP_Real* vlbconsts;
3152  SCIP_Real QUAD(zcoef);
3153  int zidx;
3154 
3155  vlbvars = SCIPvarGetVlbVars(var);
3156  vlbcoefs = SCIPvarGetVlbCoefs(var);
3157  vlbconsts = SCIPvarGetVlbConstants(var);
3158  assert(vlbvars != NULL);
3159  assert(vlbcoefs != NULL);
3160  assert(vlbconsts != NULL);
3161 
3162  assert(0 <= bestlbtypes[i] && bestlbtypes[i] < SCIPvarGetNVlbs(var));
3163  assert(SCIPvarIsActive(vlbvars[bestlbtypes[i]]));
3164  zidx = SCIPvarGetProbindex(vlbvars[bestlbtypes[i]]);
3165  assert(0 <= zidx && zidx < firstcontvar);
3166 
3167  SCIPquadprecProdQD(tmp, coef, vlbconsts[bestlbtypes[i]]);
3168  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3169 
3170  /* check if integral variable already exists in the row and update sparsity pattern */
3171  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3172  if( QUAD_HI(zcoef) == 0.0 )
3173  cutinds[(*nnz)++] = zidx;
3174 
3175  SCIPquadprecProdQD(coef, coef, vlbcoefs[bestlbtypes[i]]);
3176  SCIPquadprecSumQQ(zcoef, zcoef, coef);
3177  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3178  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3179  assert(QUAD_HI(zcoef) != 0.0);
3180  }
3181  }
3182  else
3183  {
3184  assert(!SCIPisInfinity(scip, bestubs[i]));
3185 
3186  /* use upper bound as transformation bound: x'_j := ub_j - x_j */
3187  boundtype[i] = bestubtypes[i];
3188  varsign[i] = -1;
3189 
3190  /* standard (bestubtype < 0) or variable (bestubtype >= 0) upper bound? */
3191  if( bestubtypes[i] < 0 )
3192  {
3193  SCIPquadprecProdQD(tmp, coef, bestubs[i]);
3194  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3195  *localbdsused = *localbdsused || (bestubtypes[i] == -2);
3196  }
3197  else
3198  {
3199  SCIP_VAR** vubvars;
3200  SCIP_Real* vubcoefs;
3201  SCIP_Real* vubconsts;
3202  SCIP_Real QUAD(zcoef);
3203  int zidx;
3204 
3205  vubvars = SCIPvarGetVubVars(var);
3206  vubcoefs = SCIPvarGetVubCoefs(var);
3207  vubconsts = SCIPvarGetVubConstants(var);
3208  assert(vubvars != NULL);
3209  assert(vubcoefs != NULL);
3210  assert(vubconsts != NULL);
3211 
3212  assert(0 <= bestubtypes[i] && bestubtypes[i] < SCIPvarGetNVubs(var));
3213  assert(SCIPvarIsActive(vubvars[bestubtypes[i]]));
3214  zidx = SCIPvarGetProbindex(vubvars[bestubtypes[i]]);
3215  assert(zidx >= 0);
3216 
3217  SCIPquadprecProdQD(tmp, coef, vubconsts[bestubtypes[i]]);
3218  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3219 
3220  /* check if integral variable already exists in the row and update sparsity pattern */
3221  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3222  if( QUAD_HI(zcoef) == 0.0 )
3223  cutinds[(*nnz)++] = zidx;
3224 
3225  SCIPquadprecProdQD(coef, coef, vubcoefs[bestubtypes[i]]);
3226  SCIPquadprecSumQQ(zcoef, zcoef, coef);
3227  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3228  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3229  assert(QUAD_HI(zcoef) != 0.0);
3230  }
3231  }
3232  }
3233 
3234  /* remove integral variables that now have a zero coefficient due to variable bound usage of continuous variables
3235  * and determine the bound to use for the integer variables that are left
3236  */
3237  while( i < *nnz )
3238  {
3239  SCIP_Real QUAD(coef);
3240  int v = cutinds[i];
3241  assert(cutinds[i] < firstcontvar);
3242 
3243  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3244 
3245  /* due to variable bound usage for the continous variables cancellation may have occurred */
3246  if( EPSZ(QUAD_TO_DBL(coef), QUAD_EPSILON) )
3247  {
3248  QUAD_ASSIGN(coef, 0.0);
3249  QUAD_ARRAY_STORE(cutcoefs, v, coef);
3250  --(*nnz);
3251  cutinds[i] = cutinds[*nnz];
3252  /* do not increase i, since last element is copied to the i-th position */
3253  continue;
3254  }
3255 
3256  /* determine the best bounds for the integral variable, usevbd can be set to FALSE here as vbds are only used for continous variables */
3257  SCIP_CALL( determineBestBounds(scip, vars[v], sol, boundswitch, FALSE, allowlocal, fixintegralrhs,
3258  ignoresol, boundsfortrans, boundtypesfortrans,
3259  bestlbs + i, bestubs + i, bestlbtypes + i, bestubtypes + i, selectedbounds + i, freevariable) );
3260 
3261  /* increase i */
3262  ++i;
3263 
3264  if( *freevariable )
3265  goto TERMINATE;
3266  }
3267 
3268  /* now perform the bound substitution on the remaining integral variables which only uses standard bounds */
3269  for( i = aggrrowintstart; i < *nnz; ++i )
3270  {
3271  SCIP_Real QUAD(coef);
3272  int v = cutinds[i];
3273 
3274  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3275 
3276  /* perform bound substitution */
3277  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
3278  {
3279  assert(!SCIPisInfinity(scip, - bestlbs[i]));
3280  assert(bestlbtypes[i] < 0);
3281 
3282  /* use lower bound as transformation bound: x'_j := x_j - lb_j */
3283  boundtype[i] = bestlbtypes[i];
3284  varsign[i] = +1;
3285 
3286  /* standard (bestlbtype < 0) or variable (bestlbtype >= 0) lower bound? */
3287  SCIPquadprecProdQD(tmp, coef, bestlbs[i]);
3288  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3289  *localbdsused = *localbdsused || (bestlbtypes[i] == -2);
3290  }
3291  else
3292  {
3293  assert(!SCIPisInfinity(scip, bestubs[i]));
3294  assert(bestubtypes[i] < 0);
3295 
3296  /* use upper bound as transformation bound: x'_j := ub_j - x_j */
3297  boundtype[i] = bestubtypes[i];
3298  varsign[i] = -1;
3299 
3300  /* standard (bestubtype < 0) or variable (bestubtype >= 0) upper bound? */
3301  SCIPquadprecProdQD(tmp, coef, bestubs[i]);
3302  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3303  *localbdsused = *localbdsused || (bestubtypes[i] == -2);
3304  }
3305  }
3306 
3307  if( fixintegralrhs )
3308  {
3309  SCIP_Real f0;
3310 
3311  /* check if rhs is fractional */
3312  f0 = EPSFRAC(QUAD_TO_DBL(*cutrhs), SCIPsumepsilon(scip));
3313  if( f0 < minfrac || f0 > maxfrac )
3314  {
3315  SCIP_Real bestviolgain;
3316  SCIP_Real bestnewf0;
3317  int besti;
3318 
3319  /* choose complementation of one variable differently such that f0 is in correct range */
3320  besti = -1;
3321  bestviolgain = -1e+100;
3322  bestnewf0 = 1.0;
3323  for( i = 0; i < *nnz; i++ )
3324  {
3325  int v;
3326  SCIP_Real QUAD(coef);
3327 
3328  v = cutinds[i];
3329  assert(0 <= v && v < nvars);
3330 
3331  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3332  assert(!EPSZ(QUAD_TO_DBL(coef), QUAD_EPSILON));
3333 
3334  if( boundtype[i] < 0
3335  && ((varsign[i] == +1 && !SCIPisInfinity(scip, bestubs[i]) && bestubtypes[i] < 0)
3336  || (varsign[i] == -1 && !SCIPisInfinity(scip, -bestlbs[i]) && bestlbtypes[i] < 0)) )
3337  {
3338  SCIP_Real fj;
3339  SCIP_Real newfj;
3340  SCIP_Real newrhs;
3341  SCIP_Real newf0;
3342  SCIP_Real solval;
3343  SCIP_Real viol;
3344  SCIP_Real newviol;
3345  SCIP_Real violgain;
3346 
3347  /* currently: a'_j = varsign * a_j -> f'_j = a'_j - floor(a'_j)
3348  * after complementation: a''_j = -varsign * a_j -> f''_j = a''_j - floor(a''_j) = 1 - f'_j
3349  * rhs'' = rhs' + varsign * a_j * (lb_j - ub_j)
3350  * cut violation from f0 and fj: f'_0 - f'_j * x'_j
3351  * after complementation: f''_0 - f''_j * x''_j
3352  *
3353  * for continuous variables, we just set f'_j = f''_j = |a'_j|
3354  */
3355  newrhs = QUAD_TO_DBL(*cutrhs) + varsign[i] * QUAD_TO_DBL(coef) * (bestlbs[i] - bestubs[i]);
3356  newf0 = EPSFRAC(newrhs, SCIPsumepsilon(scip));
3357  if( newf0 < minfrac || newf0 > maxfrac )
3358  continue;
3359  if( v >= firstcontvar )
3360  {
3361  fj = REALABS(QUAD_TO_DBL(coef));
3362  newfj = fj;
3363  }
3364  else
3365  {
3366  fj = SCIPfrac(scip, varsign[i] * QUAD_TO_DBL(coef));
3367  newfj = SCIPfrac(scip, -varsign[i] * QUAD_TO_DBL(coef));
3368  }
3369 
3370  if( !ignoresol )
3371  {
3372  solval = (sol == NULL ? SCIPvarGetLPSol(vars[v]) : SCIPgetSolVal(scip, sol, vars[v]));
3373  viol = f0 - fj * (varsign[i] == +1 ? solval - bestlbs[i] : bestubs[i] - solval);
3374  newviol = newf0 - newfj * (varsign[i] == -1 ? solval - bestlbs[i] : bestubs[i] - solval);
3375  violgain = newviol - viol;
3376  }
3377  else
3378  {
3379  /* todo: this should be done, this can improve the dualray significantly */
3380  SCIPerrorMessage("Cannot handle closest bounds with ignoring the LP solution.\n");
3381  return SCIP_INVALIDCALL;
3382  }
3383 
3384  /* prefer larger violations; for equal violations, prefer smaller f0 values since then the possibility that
3385  * we f_j > f_0 is larger and we may improve some coefficients in rounding
3386  */
3387  if( SCIPisGT(scip, violgain, bestviolgain)
3388  || (SCIPisGE(scip, violgain, bestviolgain) && newf0 < bestnewf0) )
3389  {
3390  besti = i;
3391  bestviolgain = violgain;
3392  bestnewf0 = newf0;
3393  }
3394  }
3395  }
3396 
3397  if( besti >= 0 )
3398  {
3399  SCIP_Real QUAD(coef);
3400  assert(besti < *nnz);
3401  assert(boundtype[besti] < 0);
3402  assert(!SCIPisInfinity(scip, -bestlbs[besti]));
3403  assert(!SCIPisInfinity(scip, bestubs[besti]));
3404 
3405  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[besti]);
3406  QUAD_SCALE(coef, varsign[besti]);
3407 
3408  /* switch the complementation of this variable */
3409  SCIPquadprecSumDD(tmp, bestlbs[besti], - bestubs[besti]);
3410  SCIPquadprecProdQQ(tmp, tmp, coef);
3411  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3412 
3413  if( varsign[besti] == +1 )
3414  {
3415  /* switch to upper bound */
3416  assert(bestubtypes[besti] < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
3417  boundtype[besti] = bestubtypes[besti];
3418  varsign[besti] = -1;
3419  }
3420  else
3421  {
3422  /* switch to lower bound */
3423  assert(bestlbtypes[besti] < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
3424  boundtype[besti] = bestlbtypes[besti];
3425  varsign[besti] = +1;
3426  }
3427  *localbdsused = *localbdsused || (boundtype[besti] == -2);
3428  }
3429  }
3430  }
3431 
3432  TERMINATE:
3433 
3434  /*free temporary memory */
3435  SCIPfreeBufferArray(scip, &selectedbounds);
3436  SCIPfreeBufferArray(scip, &bestubtypes);
3437  SCIPfreeBufferArray(scip, &bestlbtypes);
3438  SCIPfreeBufferArray(scip, &bestubs);
3439  SCIPfreeBufferArray(scip, &bestlbs);
3440 
3441  return SCIP_OKAY;
3442 }
3443 
3444 /** Calculate fractionalities \f$ f_0 := b - down(b), f_j := a^\prime_j - down(a^\prime_j) \f$, and derive MIR cut \f$ \tilde{a} \cdot x' \leq down(b) \f$
3445  * \f[
3446  * \begin{array}{rll}
3447  * integers :& \tilde{a}_j = down(a^\prime_j), & if \qquad f_j \leq f_0 \\
3448  * & \tilde{a}_j = down(a^\prime_j) + (f_j - f_0)/(1 - f_0),& if \qquad f_j > f_0 \\
3449  * continuous:& \tilde{a}_j = 0, & if \qquad a^\prime_j \geq 0 \\
3450  * & \tilde{a}_j = a^\prime_j/(1 - f_0), & if \qquad a^\prime_j < 0
3451  * \end{array}
3452  * \f]
3453  *
3454  * Transform inequality back to \f$ \hat{a} \cdot x \leq rhs \f$:
3455  *
3456  * (lb or ub):
3457  * \f[
3458  * \begin{array}{lllll}
3459  * x^\prime_j := x_j - lb_j,& x_j = x^\prime_j + lb_j,& a^\prime_j = a_j,& \hat{a}_j := \tilde{a}_j,& \mbox{if lb was used in transformation} \\
3460  * x^\prime_j := ub_j - x_j,& x_j = ub_j - x^\prime_j,& a^\prime_j = -a_j,& \hat{a}_j := -\tilde{a}_j,& \mbox{if ub was used in transformation}
3461  * \end{array}
3462  * \f]
3463  * and move the constant terms
3464  * \f[
3465  * \begin{array}{cl}
3466  * -\tilde{a}_j \cdot lb_j = -\hat{a}_j \cdot lb_j,& \mbox{or} \\
3467  * \tilde{a}_j \cdot ub_j = -\hat{a}_j \cdot ub_j &
3468  * \end{array}
3469  * \f]
3470  * to the rhs.
3471  *
3472  * (vlb or vub):
3473  * \f[
3474  * \begin{array}{lllll}
3475  * x^\prime_j := x_j - (bl_j \cdot zl_j + dl_j),& x_j = x^\prime_j + (bl_j\, zl_j + dl_j),& a^\prime_j = a_j,& \hat{a}_j := \tilde{a}_j,& \mbox{(vlb)} \\
3476  * x^\prime_j := (bu_j\, zu_j + du_j) - x_j,& x_j = (bu_j\, zu_j + du_j) - x^\prime_j,& a^\prime_j = -a_j,& \hat{a}_j := -\tilde{a}_j,& \mbox{(vub)}
3477  * \end{array}
3478  * \f]
3479  * move the constant terms
3480  * \f[
3481  * \begin{array}{cl}
3482  * -\tilde{a}_j\, dl_j = -\hat{a}_j\, dl_j,& \mbox{or} \\
3483  * \tilde{a}_j\, du_j = -\hat{a}_j\, du_j &
3484  * \end{array}
3485  * \f]
3486  * to the rhs, and update the VB variable coefficients:
3487  * \f[
3488  * \begin{array}{ll}
3489  * \hat{a}_{zl_j} := \hat{a}_{zl_j} - \tilde{a}_j\, bl_j = \hat{a}_{zl_j} - \hat{a}_j\, bl_j,& \mbox{or} \\
3490  * \hat{a}_{zu_j} := \hat{a}_{zu_j} + \tilde{a}_j\, bu_j = \hat{a}_{zu_j} - \hat{a}_j\, bu_j &
3491  * \end{array}
3492  * \f]
3493  */
3494 static
3496  SCIP* scip, /**< SCIP data structure */
3497  SCIP_Real*RESTRICT cutcoefs, /**< array of coefficients of cut */
3498  QUAD(SCIP_Real*RESTRICT cutrhs), /**< pointer to right hand side of cut */
3499  int*RESTRICT cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3500  int*RESTRICT nnz, /**< number of non-zeros in cut */
3501  int*RESTRICT varsign, /**< stores the sign of the transformed variable in summation */
3502  int*RESTRICT boundtype, /**< stores the bound used for transformed variable (vlb/vub_idx or -1 for lb/ub) */
3503  QUAD(SCIP_Real f0) /**< fractional value of rhs */
3504  )
3505 {
3506  SCIP_Real QUAD(tmp);
3507  SCIP_Real QUAD(onedivoneminusf0);
3508  int i;
3509  int firstcontvar;
3510  SCIP_VAR** vars;
3511  int ndelcontvars;
3512 
3513  assert(QUAD_HI(cutrhs) != NULL);
3514  assert(cutcoefs != NULL);
3515  assert(cutinds != NULL);
3516  assert(nnz != NULL);
3517  assert(boundtype != NULL);
3518  assert(varsign != NULL);
3519  assert(0.0 < QUAD_TO_DBL(f0) && QUAD_TO_DBL(f0) < 1.0);
3520 
3521  SCIPquadprecSumQD(onedivoneminusf0, -f0, 1.0);
3522  SCIPquadprecDivDQ(onedivoneminusf0, 1.0, onedivoneminusf0);
3523 
3524  /* Loop backwards to process integral variables first and be able to delete coefficients of integral variables
3525  * without destroying the ordering of the aggrrow's non-zeros.
3526  * (due to sorting in cutsTransformMIR the ordering is continuous before integral)
3527  */
3528 
3529  firstcontvar = SCIPgetNVars(scip) - SCIPgetNContVars(scip);
3530  vars = SCIPgetVars(scip);
3531 #ifndef NDEBUG
3532  /*in debug mode check that all continuous variables of the aggrrow come before the integral variables */
3533  i = 0;
3534  while( i < *nnz && cutinds[i] >= firstcontvar )
3535  ++i;
3536 
3537  while( i < *nnz )
3538  {
3539  assert(cutinds[i] < firstcontvar);
3540  ++i;
3541  }
3542 #endif
3543 
3544  for( i = *nnz - 1; i >= 0 && cutinds[i] < firstcontvar; --i )
3545  {
3546  SCIP_VAR* var;
3547  SCIP_Real QUAD(cutaj);
3548  int v;
3549 
3550  v = cutinds[i];
3551  assert(0 <= v && v < SCIPgetNVars(scip));
3552 
3553  var = vars[v];
3554  assert(var != NULL);
3555  assert(SCIPvarGetProbindex(var) == v);
3556  assert(varsign[i] == +1 || varsign[i] == -1);
3557 
3558  /* calculate the coefficient in the retransformed cut */
3559  {
3560  SCIP_Real QUAD(aj);
3561  SCIP_Real QUAD(downaj);
3562  SCIP_Real QUAD(fj);
3563 
3564  QUAD_ARRAY_LOAD(aj, cutcoefs, v);
3565  QUAD_SCALE(aj, varsign[i]);
3566 
3567  SCIPquadprecEpsFloorQ(downaj, aj, SCIPepsilon(scip)); /*lint !e666*/
3568  SCIPquadprecSumQQ(fj, aj, -downaj);
3569  assert(QUAD_TO_DBL(fj) >= -SCIPepsilon(scip) && QUAD_TO_DBL(fj) < 1.0);
3570 
3571  if( SCIPisLE(scip, QUAD_TO_DBL(fj), QUAD_TO_DBL(f0)) )
3572  {
3573  QUAD_ASSIGN_Q(cutaj, downaj);
3574  }
3575  else
3576  {
3577  SCIPquadprecSumQQ(tmp, fj, -f0);
3578  SCIPquadprecProdQQ(tmp, tmp, onedivoneminusf0);
3579  SCIPquadprecSumQQ(cutaj, tmp, downaj);
3580  }
3581 
3582  QUAD_SCALE(cutaj, varsign[i]);
3583  }
3584 
3585  /* remove zero cut coefficients from cut */
3586  if( EPSZ(QUAD_TO_DBL(cutaj), QUAD_EPSILON) )
3587  {
3588  QUAD_ASSIGN(cutaj, 0.0);
3589  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3590  --*nnz;
3591  cutinds[i] = cutinds[*nnz];
3592  continue;
3593  }
3594 
3595  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3596 
3597  /* integral var uses standard bound */
3598  assert(boundtype[i] < 0);
3599 
3600  /* move the constant term -a~_j * lb_j == -a^_j * lb_j , or a~_j * ub_j == -a^_j * ub_j to the rhs */
3601  if( varsign[i] == +1 )
3602  {
3603  /* lower bound was used */
3604  if( boundtype[i] == -1 )
3605  {
3606  assert(!SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)));
3607  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbGlobal(var));
3608  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetLbGlobal(var) */
3609  }
3610  else
3611  {
3612  assert(!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)));
3613  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbLocal(var));
3614  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetLbLocal(var) */
3615  }
3616  }
3617  else
3618  {
3619  /* upper bound was used */
3620  if( boundtype[i] == -1 )
3621  {
3622  assert(!SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)));
3623  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbGlobal(var));
3624  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetUbGlobal(var) */
3625  }
3626  else
3627  {
3628  assert(!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)));
3629  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbLocal(var));
3630  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetUbLocal(var) */
3631  }
3632  }
3633  }
3634 
3635  /* now process the continuous variables; postpone deletetion of zeros till all continuous variables have been processed */
3636  ndelcontvars = 0;
3637  while( i >= ndelcontvars )
3638  {
3639  SCIP_VAR* var;
3640  SCIP_Real QUAD(cutaj);
3641  int v;
3642 
3643  v = cutinds[i];
3644  assert(0 <= v && v < SCIPgetNVars(scip));
3645 
3646  var = vars[v];
3647  assert(var != NULL);
3648  assert(SCIPvarGetProbindex(var) == v);
3649  assert(varsign[i] == +1 || varsign[i] == -1);
3650  assert( v >= firstcontvar );
3651 
3652  /* calculate the coefficient in the retransformed cut */
3653  {
3654  SCIP_Real QUAD(aj);
3655 
3656  QUAD_ARRAY_LOAD(aj, cutcoefs, v);
3657 
3658  if( QUAD_TO_DBL(aj) * varsign[i] >= 0.0 )
3659  QUAD_ASSIGN(cutaj, 0.0);
3660  else
3661  SCIPquadprecProdQQ(cutaj, onedivoneminusf0, aj); /* cutaj = varsign[i] * aj * onedivoneminusf0; // a^_j */
3662  }
3663 
3664  /* remove zero cut coefficients from cut; move a continuous var from the beginning
3665  * to the current position, so that all integral variables stay behind the continuous
3666  * variables
3667  */
3668  if( EPSZ(QUAD_TO_DBL(cutaj), QUAD_EPSILON) )
3669  {
3670  QUAD_ASSIGN(cutaj, 0.0);
3671  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3672  cutinds[i] = cutinds[ndelcontvars];
3673  varsign[i] = varsign[ndelcontvars];
3674  boundtype[i] = boundtype[ndelcontvars];
3675  ++ndelcontvars;
3676  continue;
3677  }
3678 
3679  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3680 
3681  /* check for variable bound use */
3682  if( boundtype[i] < 0 )
3683  {
3684  /* standard bound */
3685 
3686  /* move the constant term -a~_j * lb_j == -a^_j * lb_j , or a~_j * ub_j == -a^_j * ub_j to the rhs */
3687  if( varsign[i] == +1 )
3688  {
3689  /* lower bound was used */
3690  if( boundtype[i] == -1 )
3691  {
3692  assert(!SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)));
3693  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbGlobal(var));
3694  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3695  }
3696  else
3697  {
3698  assert(!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)));
3699  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbLocal(var));
3700  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3701  }
3702  }
3703  else
3704  {
3705  /* upper bound was used */
3706  if( boundtype[i] == -1 )
3707  {
3708  assert(!SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)));
3709  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbGlobal(var));
3710  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3711  }
3712  else
3713  {
3714  assert(!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)));
3715  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbLocal(var));
3716  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3717  }
3718  }
3719  }
3720  else
3721  {
3722  SCIP_VAR** vbz;
3723  SCIP_Real* vbb;
3724  SCIP_Real* vbd;
3725  SCIP_Real QUAD(zcoef);
3726  int vbidx;
3727  int zidx;
3728 
3729  /* variable bound */
3730  vbidx = boundtype[i];
3731 
3732  /* change mirrhs and cutaj of integer variable z_j of variable bound */
3733  if( varsign[i] == +1 )
3734  {
3735  /* variable lower bound was used */
3736  assert(0 <= vbidx && vbidx < SCIPvarGetNVlbs(var));
3737  vbz = SCIPvarGetVlbVars(var);
3738  vbb = SCIPvarGetVlbCoefs(var);
3739  vbd = SCIPvarGetVlbConstants(var);
3740  }
3741  else
3742  {
3743  /* variable upper bound was used */
3744  assert(0 <= vbidx && vbidx < SCIPvarGetNVubs(var));
3745  vbz = SCIPvarGetVubVars(var);
3746  vbb = SCIPvarGetVubCoefs(var);
3747  vbd = SCIPvarGetVubConstants(var);
3748  }
3749  assert(SCIPvarIsActive(vbz[vbidx]));
3750  zidx = SCIPvarGetProbindex(vbz[vbidx]);
3751  assert(0 <= zidx && zidx < firstcontvar);
3752 
3753  SCIPquadprecProdQD(tmp, cutaj, vbd[vbidx]);
3754  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3755 
3756  SCIPquadprecProdQD(tmp, cutaj, vbb[vbidx]);
3757  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3758 
3759  /* update sparsity pattern */
3760  if( QUAD_HI(zcoef) == 0.0 )
3761  cutinds[(*nnz)++] = zidx;
3762 
3763  SCIPquadprecSumQQ(zcoef, zcoef, -tmp);
3764  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3765  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3766  assert(QUAD_HI(zcoef) != 0.0);
3767  }
3768 
3769  /* advance to next variable */
3770  --i;
3771  }
3772 
3773  /* fill the empty position due to deleted continuous variables */
3774  if( ndelcontvars > 0 )
3775  {
3776  assert(ndelcontvars <= *nnz);
3777  *nnz -= ndelcontvars;
3778  if( *nnz < ndelcontvars )
3779  {
3780  BMScopyMemoryArray(cutinds, cutinds + ndelcontvars, *nnz);
3781  }
3782  else
3783  {
3784  BMScopyMemoryArray(cutinds, cutinds + *nnz, ndelcontvars);
3785  }
3786  }
3787 
3788  return SCIP_OKAY;
3789 }
3790 
3791 /** substitute aggregated slack variables:
3792  *
3793  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
3794  * variable only appears in its own row: \f$ a^\prime_r = scale * weight[r] * slacksign[r]. \f$
3795  *
3796  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
3797  * \f[
3798  * \begin{array}{rll}
3799  * integers : & \hat{a}_r = \tilde{a}_r = down(a^\prime_r), & \mbox{if}\qquad f_r <= f0 \\
3800  * & \hat{a}_r = \tilde{a}_r = down(a^\prime_r) + (f_r - f0)/(1 - f0),& \mbox{if}\qquad f_r > f0 \\
3801  * continuous:& \hat{a}_r = \tilde{a}_r = 0, & \mbox{if}\qquad a^\prime_r >= 0 \\
3802  * & \hat{a}_r = \tilde{a}_r = a^\prime_r/(1 - f0), & \mbox{if}\qquad a^\prime_r < 0
3803  * \end{array}
3804  * \f]
3805  *
3806  * Substitute \f$ \hat{a}_r \cdot s_r \f$ by adding \f$ \hat{a}_r \f$ times the slack's definition to the cut.
3807  */
3808 static
3810  SCIP* scip, /**< SCIP data structure */
3811  SCIP_Real* weights, /**< row weights in row summation */
3812  int* slacksign, /**< stores the sign of the row's slack variable in summation */
3813  int* rowinds, /**< sparsity pattern of used rows */
3814  int nrowinds, /**< number of used rows */
3815  SCIP_Real scale, /**< additional scaling factor multiplied to all rows */
3816  SCIP_Real* cutcoefs, /**< array of coefficients of cut */
3817  QUAD(SCIP_Real* cutrhs), /**< pointer to right hand side of cut */
3818  int* cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3819  int* nnz, /**< number of non-zeros in cut */
3820  QUAD(SCIP_Real f0) /**< fractional value of rhs */
3821  )
3822 { /*lint --e{715}*/
3823  SCIP_ROW** rows;
3824  SCIP_Real QUAD(onedivoneminusf0);
3825  int i;
3826 
3827  assert(scip != NULL);
3828  assert(weights != NULL || nrowinds == 0);
3829  assert(slacksign != NULL || nrowinds == 0);
3830  assert(rowinds != NULL || nrowinds == 0);
3831  assert(scale > 0.0);
3832  assert(cutcoefs != NULL);
3833  assert(QUAD_HI(cutrhs) != NULL);
3834  assert(cutinds != NULL);
3835  assert(nnz != NULL);
3836  assert(0.0 < QUAD_TO_DBL(f0) && QUAD_TO_DBL(f0) < 1.0);
3837 
3838  SCIPquadprecSumQD(onedivoneminusf0, -f0, 1.0);
3839  SCIPquadprecDivDQ(onedivoneminusf0, 1.0, onedivoneminusf0);
3840 
3841  rows = SCIPgetLPRows(scip);
3842  for( i = 0; i < nrowinds; i++ )
3843  {
3844  SCIP_ROW* row;
3845  SCIP_Real ar;
3846  SCIP_Real downar;
3847  SCIP_Real QUAD(cutar);
3848  SCIP_Real QUAD(fr);
3849  SCIP_Real QUAD(tmp);
3850  SCIP_Real mul;
3851  int r;
3852 
3853  r = rowinds[i]; /*lint !e613*/
3854  assert(0 <= r && r < SCIPgetNLPRows(scip));
3855  assert(slacksign[i] == -1 || slacksign[i] == +1); /*lint !e613*/
3856  assert(!SCIPisZero(scip, weights[i])); /*lint !e613*/
3857 
3858  row = rows[r];
3859  assert(row != NULL);
3860  assert(row->len == 0 || row->cols != NULL);
3861  assert(row->len == 0 || row->cols_index != NULL);
3862  assert(row->len == 0 || row->vals != NULL);
3863 
3864  /* get the slack's coefficient a'_r in the aggregated row */
3865  ar = slacksign[i] * scale * weights[i]; /*lint !e613*/
3866 
3867  /* calculate slack variable's coefficient a^_r in the cut */
3868  if( row->integral
3869  && ((slacksign[i] == +1 && SCIPisFeasIntegral(scip, row->rhs - row->constant))
3870  || (slacksign[i] == -1 && SCIPisFeasIntegral(scip, row->lhs - row->constant))) ) /*lint !e613*/
3871  {
3872  /* slack variable is always integral:
3873  * a^_r = a~_r = down(a'_r) , if f_r <= f0
3874  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
3875  */
3876  downar = EPSFLOOR(ar, QUAD_EPSILON);
3877  SCIPquadprecSumDD(fr, ar, -downar);
3878  if( SCIPisLE(scip, QUAD_TO_DBL(fr), QUAD_TO_DBL(f0)) )
3879  QUAD_ASSIGN(cutar, downar);
3880  else
3881  {
3882  SCIPquadprecSumQQ(cutar, fr, -f0);
3883  SCIPquadprecProdQQ(cutar, cutar, onedivoneminusf0);
3884  SCIPquadprecSumQD(cutar, cutar, downar);
3885  }
3886  }
3887  else
3888  {
3889  /* slack variable is continuous:
3890  * a^_r = a~_r = 0 , if a'_r >= 0
3891  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
3892  */
3893  if( ar >= 0.0 )
3894  continue; /* slack can be ignored, because its coefficient is reduced to 0.0 */
3895  else
3896  SCIPquadprecProdQD(cutar, onedivoneminusf0, ar);
3897  }
3898 
3899  /* if the coefficient was reduced to zero, ignore the slack variable */
3900  if( EPSZ(QUAD_TO_DBL(cutar), QUAD_EPSILON) )
3901  continue;
3902 
3903  /* depending on the slack's sign, we have
3904  * a*x + c + s == rhs => s == - a*x - c + rhs, or a*x + c - s == lhs => s == a*x + c - lhs
3905  * substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
3906  */
3907  mul = -slacksign[i] * QUAD_TO_DBL(cutar); /*lint !e613*/
3908 
3909  /* add the slack's definition multiplied with a^_j to the cut */
3910  SCIP_CALL( varVecAddScaledRowCoefsQuad(cutinds, cutcoefs, nnz, row, mul) );
3911 
3912  /* move slack's constant to the right hand side */
3913  if( slacksign[i] == +1 ) /*lint !e613*/
3914  {
3915  SCIP_Real QUAD(rowrhs);
3916 
3917  /* a*x + c + s == rhs => s == - a*x - c + rhs: move a^_r * (rhs - c) to the right hand side */
3918  assert(!SCIPisInfinity(scip, row->rhs));
3919  SCIPquadprecSumDD(rowrhs, row->rhs, -row->constant);
3920  if( row->integral )
3921  {
3922  /* the right hand side was implicitly rounded down in row aggregation */
3923  QUAD_ASSIGN(rowrhs, SCIPfloor(scip, QUAD_TO_DBL(rowrhs)));
3924  }
3925  SCIPquadprecProdQQ(tmp, cutar, rowrhs);
3926  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3927  }
3928  else
3929  {
3930  SCIP_Real QUAD(rowlhs);
3931 
3932  /* a*x + c - s == lhs => s == a*x + c - lhs: move a^_r * (c - lhs) to the right hand side */
3933  assert(!SCIPisInfinity(scip, -row->lhs));
3934  SCIPquadprecSumDD(rowlhs, row->lhs, -row->constant);
3935  if( row->integral )
3936  {
3937  /* the left hand side was implicitly rounded up in row aggregation */
3938  QUAD_ASSIGN(rowlhs, SCIPceil(scip, QUAD_TO_DBL(rowlhs)));
3939  }
3940  SCIPquadprecProdQQ(tmp, cutar, rowlhs);
3941  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3942  }
3943  }
3944 
3945  /* relax rhs to zero, if it's very close to */
3946  if( QUAD_TO_DBL(*cutrhs) < 0.0 && QUAD_TO_DBL(*cutrhs) >= SCIPepsilon(scip) )
3947  QUAD_ASSIGN(*cutrhs, 0.0);
3948 
3949  return SCIP_OKAY;
3950 }
3951 
3952 /** calculates an MIR cut out of the weighted sum of LP rows; The weights of modifiable rows are set to 0.0, because
3953  * these rows cannot participate in an MIR cut.
3954  *
3955  * @return \ref SCIP_OKAY is returned if everything worked. Otherwise a suitable error code is passed. See \ref
3956  * SCIP_Retcode "SCIP_RETCODE" for a complete list of error codes.
3957  *
3958  * @pre This method can be called if @p scip is in one of the following stages:
3959  * - \ref SCIP_STAGE_SOLVING
3960  *
3961  * See \ref SCIP_Stage "SCIP_STAGE" for a complete list of all possible solving stages.
3962  */
3964  SCIP* scip, /**< SCIP data structure */
3965  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
3966  SCIP_Bool postprocess, /**< apply a post-processing step to the resulting cut? */
3967  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
3968  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
3969  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
3970  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
3971  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
3972  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
3973  * NULL for using closest bound for all variables */
3974  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
3975  * NULL for using closest bound for all variables */
3976  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
3977  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
3978  SCIP_Real scale, /**< additional scaling factor multiplied to the aggrrow; must be positive */
3979  SCIP_AGGRROW* aggrrow, /**< aggrrow to compute MIR cut for */
3980  SCIP_Real* cutcoefs, /**< array to store the non-zero coefficients in the cut */
3981  SCIP_Real* cutrhs, /**< pointer to store the right hand side of the cut */
3982  int* cutinds, /**< array to store the problem indices of variables with a non-zero coefficient in the cut */
3983  int* cutnnz, /**< pointer to store the number of non-zeros in the cut */
3984  SCIP_Real* cutefficacy, /**< pointer to store efficacy of cut, or NULL */
3985  int* cutrank, /**< pointer to return rank of generated cut */
3986  SCIP_Bool* cutislocal, /**< pointer to store whether the generated cut is only valid locally */
3987  SCIP_Bool* success /**< pointer to store whether the returned coefficients are a valid MIR cut */
3988  )
3989 {
3990  int i;
3991  int nvars;
3992  int* varsign;
3993  int* boundtype;
3994  SCIP_Real* tmpcoefs;
3995 
3996  SCIP_Real QUAD(rhs);
3997  SCIP_Real QUAD(downrhs);
3998  SCIP_Real QUAD(f0);
3999  SCIP_Bool freevariable;
4000  SCIP_Bool localbdsused;
4001 
4002  assert(aggrrow != NULL);
4003  assert(SCIPisPositive(scip, scale));
4004  assert(success != NULL);
4005 
4006  SCIPdebugMessage("calculating MIR cut (scale: %g)\n", scale);
4007 
4008  *success = FALSE;
4009 
4010  /* allocate temporary memory */
4011  nvars = SCIPgetNVars(scip);
4012  SCIP_CALL( SCIPallocBufferArray(scip, &varsign, nvars) );
4013  SCIP_CALL( SCIPallocBufferArray(scip, &boundtype, nvars) );
4014  SCIP_CALL( SCIPallocCleanBufferArray(scip, &tmpcoefs, QUAD_ARRAY_SIZE(nvars)) );
4015 
4016  /* initialize cut with aggregation */
4017  *cutnnz = aggrrow->nnz;
4018  *cutislocal = aggrrow->local;
4019 
4020  SCIPquadprecProdQD(rhs, aggrrow->rhs, scale);
4021 
4022  if( *cutnnz > 0 )
4023  {
4024  BMScopyMemoryArray(cutinds, aggrrow->inds, *cutnnz);
4025 
4026  for( i = 0; i < *cutnnz; ++i )
4027  {
4028  SCIP_Real QUAD(coef);
4029 
4030  int k = aggrrow->inds[i];
4031  QUAD_ARRAY_LOAD(coef, aggrrow->vals, k);
4032 
4033  SCIPquadprecProdQD(coef, coef, scale);
4034 
4035  QUAD_ARRAY_STORE(tmpcoefs, k, coef);
4036 
4037  assert(QUAD_HI(coef) != 0.0);
4038  }
4039 
4040  /* Transform equation a*x == b, lb <= x <= ub into standard form
4041  * a'*x' == b, 0 <= x' <= ub'.
4042  *
4043  * Transform variables (lb or ub):
4044  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, if lb is used in transformation
4045  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, if ub is used in transformation
4046  * and move the constant terms "a_j * lb_j" or "a_j * ub_j" to the rhs.
4047  *
4048  * Transform variables (vlb or vub):
4049  * x'_j := x_j - (bl_j * zl_j + dl_j), x_j == x'_j + (bl_j * zl_j + dl_j), a'_j == a_j, if vlb is used in transf.
4050  * x'_j := (bu_j * zu_j + du_j) - x_j, x_j == (bu_j * zu_j + du_j) - x'_j, a'_j == -a_j, if vub is used in transf.
4051  * move the constant terms "a_j * dl_j" or "a_j * du_j" to the rhs, and update the coefficient of the VLB variable:
4052  * a_{zl_j} := a_{zl_j} + a_j * bl_j, or
4053  * a_{zu_j} := a_{zu_j} + a_j * bu_j
4054  */
4055  SCIP_CALL( cutsTransformMIR(scip, sol, boundswitch, usevbds, allowlocal, fixintegralrhs, FALSE,
4056  boundsfortrans, boundtypesfortrans, minfrac, maxfrac, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, varsign, boundtype, &freevariable, &localbdsused) );
4057  assert(allowlocal || !localbdsused);
4058  *cutislocal = *cutislocal || localbdsused;
4059 
4060  if( freevariable )
4061  goto TERMINATE;
4062  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4063  }
4064 
4065  /* Calculate fractionalities f_0 := b - down(b), f_j := a'_j - down(a'_j) , and derive MIR cut
4066  * a~*x' <= down(b)
4067  * integers : a~_j = down(a'_j) , if f_j <= f_0
4068  * a~_j = down(a'_j) + (f_j - f0)/(1 - f0), if f_j > f_0
4069  * continuous: a~_j = 0 , if a'_j >= 0
4070  * a~_j = a'_j/(1 - f0) , if a'_j < 0
4071  *
4072  * Transform inequality back to a^*x <= rhs:
4073  *
4074  * (lb or ub):
4075  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, a^_j := a~_j, if lb was used in transformation
4076  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, a^_j := -a~_j, if ub was used in transformation
4077  * and move the constant terms
4078  * -a~_j * lb_j == -a^_j * lb_j, or
4079  * a~_j * ub_j == -a^_j * ub_j
4080  * to the rhs.
4081  *
4082  * (vlb or vub):
4083  * x'_j := x_j - (bl_j * zl_j + dl_j), x_j == x'_j + (bl_j * zl_j + dl_j), a'_j == a_j, a^_j := a~_j, (vlb)
4084  * x'_j := (bu_j * zu_j + du_j) - x_j, x_j == (bu_j * zu_j + du_j) - x'_j, a'_j == -a_j, a^_j := -a~_j, (vub)
4085  * move the constant terms
4086  * -a~_j * dl_j == -a^_j * dl_j, or
4087  * a~_j * du_j == -a^_j * du_j
4088  * to the rhs, and update the VB variable coefficients:
4089  * a^_{zl_j} := a^_{zl_j} - a~_j * bl_j == a^_{zl_j} - a^_j * bl_j, or
4090  * a^_{zu_j} := a^_{zu_j} + a~_j * bu_j == a^_{zu_j} - a^_j * bu_j
4091  */
4092  SCIPquadprecEpsFloorQ(downrhs, rhs, SCIPepsilon(scip)); /*lint !e666*/
4093 
4094  SCIPquadprecSumQQ(f0, rhs, -downrhs);
4095 
4096  if( QUAD_TO_DBL(f0) < minfrac || QUAD_TO_DBL(f0) > maxfrac )
4097  goto TERMINATE;
4098 
4099  /* We multiply the coefficients of the base inequality roughly by scale/(1-f0).
4100  * If this gives a scalar that is very big, we better do not generate this cut.
4101  */
4102  if( REALABS(scale)/(1.0 - QUAD_TO_DBL(f0)) > MAXCMIRSCALE )
4103  goto TERMINATE;
4104 
4105  /* renormalize f0 value */
4106  SCIPquadprecSumDD(f0, QUAD_HI(f0), QUAD_LO(f0));
4107 
4108  QUAD_ASSIGN_Q(rhs, downrhs);
4109 
4110  if( *cutnnz > 0 )
4111  {
4112  SCIP_CALL( cutsRoundMIR(scip, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, varsign, boundtype, QUAD(f0)) );
4113  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4114  }
4115 
4116  /* substitute aggregated slack variables:
4117  *
4118  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
4119  * variable only appears in its own row:
4120  * a'_r = scale * weight[r] * slacksign[r].
4121  *
4122  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
4123  * integers : a^_r = a~_r = down(a'_r) , if f_r <= f0
4124  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
4125  * continuous: a^_r = a~_r = 0 , if a'_r >= 0
4126  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
4127  *
4128  * Substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
4129  */
4130  SCIP_CALL( cutsSubstituteMIR(scip, aggrrow->rowweights, aggrrow->slacksign, aggrrow->rowsinds,
4131  aggrrow->nrows, scale, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, QUAD(f0)) );
4132  SCIPdebug( printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE) );
4133 
4134  if( postprocess )
4135  {
4136  /* remove all nearly-zero coefficients from MIR row and relax the right hand side correspondingly in order to
4137  * prevent numerical rounding errors
4138  */
4139  SCIP_CALL( postprocessCutQuad(scip, *cutislocal, cutinds, tmpcoefs, cutnnz, QUAD(&rhs), success) );
4140  }
4141  else
4142  {
4143  *success = ! removeZerosQuad(scip, SCIPsumepsilon(scip), *cutislocal, tmpcoefs, QUAD(&rhs), cutnnz, cutinds);
4144  }
4145 
4146  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4147 
4148  if( *success )
4149  {
4150  *cutrhs = QUAD_TO_DBL(rhs);
4151 
4152  /* clean tmpcoefs and go back to double precision */
4153  for( i = 0; i < *cutnnz; ++i )
4154  {
4155  SCIP_Real QUAD(coef);
4156  int j = cutinds[i];
4157 
4158  QUAD_ARRAY_LOAD(coef, tmpcoefs, j);
4159 
4160  cutcoefs[i] = QUAD_TO_DBL(coef);
4161  QUAD_ASSIGN(coef, 0.0);
4162  QUAD_ARRAY_STORE(tmpcoefs, j, coef);
4163  }
4164 
4165  if( cutefficacy != NULL )
4166  *cutefficacy = calcEfficacy(scip, sol, cutcoefs, *cutrhs, cutinds, *cutnnz);
4167 
4168  if( cutrank != NULL )
4169  *cutrank = aggrrow->rank + 1;
4170  }
4171 
4172  TERMINATE:
4173  if( !(*success) )
4174  {
4175  SCIP_Real QUAD(tmp);
4176 
4177  QUAD_ASSIGN(tmp, 0.0);
4178  for( i = 0; i < *cutnnz; ++i )
4179  {
4180  QUAD_ARRAY_STORE(tmpcoefs, cutinds[i], tmp);
4181  }
4182  }
4183 
4184  /* free temporary memory */
4185  SCIPfreeCleanBufferArray(scip, &tmpcoefs);
4186  SCIPfreeBufferArray(scip, &boundtype);
4187  SCIPfreeBufferArray(scip, &varsign);
4188 
4189  return SCIP_OKAY;
4190 }
4191 
4192 /** compute the efficacy of the MIR cut for the given values without computing the cut.
4193  * This is used for the CMIR cut generation heuristic.
4194  */
4195 static
4197  SCIP* scip, /**< SCIP datastructure */
4198  SCIP_Real*RESTRICT coefs, /**< array with coefficients in row */
4199  SCIP_Real*RESTRICT solvals, /**< solution values of variables in the row */
4200  SCIP_Real rhs, /**< right hand side of MIR cut */
4201  SCIP_Real contactivity, /**< aggregated activity of continuous variables in the row */
4202  SCIP_Real contsqrnorm, /**< squared norm of continuous variables */
4203  SCIP_Real delta, /**< delta value to compute the violation for */
4204  int nvars, /**< number of variables in the row, i.e. the size of coefs and solvals arrays */
4205  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
4206  SCIP_Real maxfrac /**< maximal fractionality of rhs to produce MIR cut for */
4207  )
4208 {
4209  int i;
4210  SCIP_Real f0pluseps;
4211  SCIP_Real f0;
4212  SCIP_Real onedivoneminusf0;
4213  SCIP_Real scale;
4214  SCIP_Real downrhs;
4215  SCIP_Real norm;
4216  SCIP_Real contscale;
4217 
4218  scale = 1.0 / delta;
4219 
4220  rhs *= scale;
4221 
4222  downrhs = SCIPfloor(scip, rhs);
4223 
4224  f0 = rhs - downrhs;
4225 
4226  if( f0 < minfrac || f0 > maxfrac )
4227  return 0.0;
4228 
4229  onedivoneminusf0 = 1.0 / (1.0 - f0);
4230 
4231  contscale = scale * onedivoneminusf0;
4232 
4233  /* We multiply the coefficients of the base inequality roughly by scale/(1-f0).
4234  * If this gives a scalar that is very big, we better do not generate this cut.
4235  */
4236  if( contscale > MAXCMIRSCALE )
4237  return 0.0;
4238 
4239  rhs = downrhs;
4240  rhs -= contscale * contactivity;
4241  norm = SQR(contscale) * contsqrnorm;
4242 
4243  assert(!SCIPisFeasZero(scip, f0));
4244  assert(!SCIPisFeasZero(scip, 1.0 - f0));
4245 
4246  f0pluseps = f0 + SCIPepsilon(scip);
4247 
4248  for( i = 0; i < nvars; ++i )
4249  {
4250  SCIP_Real floorai = floor(scale * coefs[i]);
4251  SCIP_Real fi = (scale * coefs[i]) - floorai;
4252 
4253  if( fi > f0pluseps )
4254  floorai += (fi - f0) * onedivoneminusf0;
4255 
4256  rhs -= solvals[i] * floorai;
4257  norm += SQR(floorai);
4258  }
4259 
4260  norm = SQRT(norm);
4261 
4262  return - rhs / MAX(norm, 1e-6);
4263 }
4264 
4265 /** calculates an MIR cut out of an aggregation of LP rows
4266  *
4267  * Given the aggregation, it is transformed to a mixed knapsack set via complementation (using bounds or variable bounds)
4268  * Then, different scalings of the mkset are used to generate a MIR and the best is chosen.
4269  * One of the steps of the MIR is to round the coefficients of the integer variables down,
4270  * so one would prefer to have integer coefficients for integer variables which are far away from their bounds in the
4271  * mkset.
4272  *
4273  * @return \ref SCIP_OKAY is returned if everything worked. Otherwise a suitable error code is passed. See \ref
4274  * SCIP_Retcode "SCIP_RETCODE" for a complete list of error codes.
4275  *
4276  * @pre This method can be called if @p scip is in one of the following stages:
4277  * - \ref SCIP_STAGE_SOLVING
4278  *
4279  * See \ref SCIP_Stage "SCIP_STAGE" for a complete list of all possible solving stages.
4280  */
4282  SCIP* scip, /**< SCIP data structure */
4283  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
4284  SCIP_Bool postprocess, /**< apply a post-processing step to the resulting cut? */
4285  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
4286  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
4287  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
4288  int maxtestdelta, /**< maximum number of deltas to test */
4289  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
4290  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
4291  * NULL for using closest bound for all variables */
4292  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
4293  * NULL for using closest bound for all variables */
4294  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
4295  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
4296  SCIP_AGGRROW* aggrrow, /**< aggrrow to compute MIR cut for */
4297  SCIP_Real* cutcoefs, /**< array to store the non-zero coefficients in the cut */
4298  SCIP_Real* cutrhs, /**< pointer to store the right hand side of the cut */
4299  int* cutinds, /**< array to store the problem indices of variables with a non-zero coefficient in the cut */
4300  int* cutnnz, /**< pointer to store the number of non-zeros in the cut */
4301  SCIP_Real* cutefficacy, /**< pointer to store efficacy of best cut; only cuts that are strictly better than the value of
4302  * this efficacy on input to this function are returned */
4303  int* cutrank, /**< pointer to return rank of generated cut (or NULL) */
4304  SCIP_Bool* cutislocal, /**< pointer to store whether the generated cut is only valid locally */
4305  SCIP_Bool* success /**< pointer to store whether a valid and efficacious cut was returned */
4306  )
4307 {
4308  int i;
4309  int firstcontvar;
4310  int nvars;
4311  int intstart;
4312  int ntmpcoefs;
4313  int* varsign;
4314  int* boundtype;
4315  int* mksetinds;
4316  SCIP_Real* mksetcoefs;
4317  SCIP_Real QUAD(mksetrhs);
4318  int mksetnnz;
4319  SCIP_Real* bounddist;
4320  int* bounddistpos;
4321  int nbounddist;
4322  SCIP_Real* tmpcoefs;
4323  SCIP_Real* tmpvalues;
4324  SCIP_Real* deltacands;
4325  int ndeltacands;
4326  SCIP_Real bestdelta;
4327  SCIP_Real bestefficacy;
4328  SCIP_Real maxabsmksetcoef;
4329  SCIP_VAR** vars;
4330  SCIP_Bool freevariable;
4331  SCIP_Bool localbdsused;
4332  SCIP_Real contactivity;
4333  SCIP_Real contsqrnorm;
4334 
4335  assert(aggrrow != NULL);
4336  assert(aggrrow->nrows + aggrrow->nnz >= 1);
4337  assert(success != NULL);
4338 
4339  *success = FALSE;
4340  nvars = SCIPgetNVars(scip);
4341  firstcontvar = nvars - SCIPgetNContVars(scip);
4342  vars = SCIPgetVars(scip);
4343 
4344  /* allocate temporary memory */
4345  SCIP_CALL( SCIPallocBufferArray(scip, &varsign, nvars) );
4346  SCIP_CALL( SCIPallocBufferArray(scip, &boundtype, nvars) );
4347  SCIP_CALL( SCIPallocCleanBufferArray(scip, &mksetcoefs, QUAD_ARRAY_SIZE(nvars)) );
4348  SCIP_CALL( SCIPallocBufferArray(scip, &mksetinds, nvars) );
4349  SCIP_CALL( SCIPallocBufferArray(scip, &tmpcoefs, nvars + aggrrow->nrows) );
4350  SCIP_CALL( SCIPallocBufferArray(scip, &tmpvalues, nvars + aggrrow->nrows) );
4351  SCIP_CALL( SCIPallocBufferArray(scip, &deltacands, aggrrow->nnz + 6) );
4352 
4353  /* we only compute bound distance for integer variables; we allocate an array of length aggrrow->nnz to store this, since
4354  * this is the largest number of integer variables. (in contrast to the number of total variables which can be 2 *
4355  * aggrrow->nnz variables: if all are continuous and we use variable bounds to completement, we introduce aggrrow->nnz
4356  * extra vars)
4357  */
4358  SCIP_CALL( SCIPallocBufferArray(scip, &bounddist, aggrrow->nnz) );
4359  SCIP_CALL( SCIPallocBufferArray(scip, &bounddistpos, aggrrow->nnz) );
4360 
4361  /* initialize mkset with aggregation */
4362  mksetnnz = aggrrow->nnz;
4363  QUAD_ASSIGN_Q(mksetrhs, aggrrow->rhs);
4364 
4365  BMScopyMemoryArray(mksetinds, aggrrow->inds, mksetnnz);
4366 
4367  for( i = 0; i < mksetnnz; ++i )
4368  {
4369  int j = mksetinds[i];
4370  SCIP_Real QUAD(coef);
4371  QUAD_ARRAY_LOAD(coef, aggrrow->vals, j);
4372  QUAD_ARRAY_STORE(mksetcoefs, j, coef);
4373  assert(QUAD_HI(coef) != 0.0);
4374  }
4375 
4376  *cutislocal = aggrrow->local;
4377 
4378  /* Transform equation a*x == b, lb <= x <= ub into standard form
4379  * a'*x' == b, 0 <= x' <= ub'.
4380  *
4381  * Transform variables (lb or ub):
4382  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, if lb is used in transformation
4383  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, if ub is used in transformation
4384  * and move the constant terms "a_j * lb_j" or "a_j * ub_j" to the rhs.
4385  *
4386  * Transform variables (vlb or vub):
4387  * x'_j := x_j - (bl_j * zl_j + dl_j), x_j == x'_j + (bl_j * zl_j + dl_j), a'_j == a_j, if vlb is used in transf.
4388  * x'_j := (bu_j * zu_j + du_j) - x_j, x_j == (bu_j * zu_j + du_j) - x'_j, a'_j == -a_j, if vub is used in transf.
4389  * move the constant terms "a_j * dl_j" or "a_j * du_j" to the rhs, and update the coefficient of the VLB variable:
4390  * a_{zl_j} := a_{zl_j} + a_j * bl_j, or
4391  * a_{zu_j} := a_{zu_j} + a_j * bu_j
4392  */
4393  SCIP_CALL( cutsTransformMIR(scip, sol, boundswitch, usevbds, allowlocal, FALSE, FALSE,
4394  boundsfortrans, boundtypesfortrans, minfrac, maxfrac, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, varsign, boundtype, &freevariable, &localbdsused) );
4395 
4396  assert(allowlocal || !localbdsused);
4397 
4398  if( freevariable )
4399  goto TERMINATE;
4400 
4401  SCIPdebugMessage("transformed aggrrow row:\n");
4402  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4403 
4404  /* found positions of integral variables that are strictly between their bounds */
4405  maxabsmksetcoef = -1.0;
4406  nbounddist = 0;
4407 
4408  for( i = mksetnnz - 1; i >= 0 && mksetinds[i] < firstcontvar; --i )
4409  {
4410  SCIP_VAR* var = vars[mksetinds[i]];
4411  SCIP_Real primsol = SCIPgetSolVal(scip, sol, var);
4412  SCIP_Real lb = SCIPvarGetLbLocal(var);
4413  SCIP_Real ub = SCIPvarGetUbLocal(var);
4414  SCIP_Real QUAD(coef);
4415 
4416  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4417 
4418  if( SCIPisEQ(scip, primsol, lb) || SCIPisEQ(scip, primsol, ub) )
4419  continue;
4420 
4421  bounddist[nbounddist] = MIN(ub - primsol, primsol - lb);
4422  bounddistpos[nbounddist] = i;
4423  deltacands[nbounddist] = QUAD_TO_DBL(coef);
4424  ++nbounddist;
4425  }
4426 
4427  /* no fractional variable; so abort here */
4428  if( nbounddist == 0 )
4429  goto TERMINATE;
4430 
4431  intstart = i + 1;
4432  ndeltacands = nbounddist;
4433 
4434  SCIPsortDownRealRealInt(bounddist, deltacands, bounddistpos, nbounddist);
4435 
4436  {
4437  SCIP_Real intscale;
4438  SCIP_Bool intscalesuccess;
4439 
4440  SCIP_CALL( SCIPcalcIntegralScalar(deltacands, nbounddist, -QUAD_EPSILON, SCIPsumepsilon(scip), (SCIP_Longint)10000, 10000.0, &intscale, &intscalesuccess) );
4441 
4442  if( intscalesuccess )
4443  {
4444  SCIP_Real intf0;
4445  SCIP_Real intscalerhs;
4446  SCIP_Real delta;
4447 
4448  intscalerhs = QUAD_TO_DBL(mksetrhs) * intscale;
4449  delta = 1.0 / intscale;
4450  intf0 = intscalerhs - floor(intscalerhs);
4451 
4452  if( ! SCIPisFeasIntegral(scip, intf0) )
4453  {
4454  if( intf0 < minfrac || intf0 > maxfrac )
4455  {
4456  intscale *= ceil(MAX(minfrac, (1.0 - maxfrac)) / MIN(intf0, (1.0 - intf0)));
4457  intscalerhs = QUAD_TO_DBL(mksetrhs) * intscale;
4458  delta = 1.0 / intscale;
4459  intf0 = intscalerhs - floor(intscalerhs);
4460  }
4461 
4462  if( intf0 >= minfrac && intf0 <= maxfrac )
4463  {
4464  if( ! SCIPisEQ(scip, delta, 1.0) )
4465  {
4466  deltacands[ndeltacands++] = delta;
4467  }
4468 
4469  if( intf0 < maxfrac )
4470  {
4471  SCIP_Real delta2;
4472 
4473  delta2 = 1.0 / (intscale * floor(maxfrac / intf0));
4474 
4475  if( ! SCIPisEQ(scip, delta, delta2) && ! SCIPisEQ(scip, delta2, 1.0) )
4476  {
4477  deltacands[ndeltacands++] = delta2;
4478  }
4479  }
4480  }
4481  }
4482  }
4483  }
4484 
4485  for( i = 0; i < nbounddist; ++i )
4486  {
4487  SCIP_Real absmksetcoef;
4488 
4489  absmksetcoef = REALABS(deltacands[i]);
4490  maxabsmksetcoef = MAX(absmksetcoef, maxabsmksetcoef);
4491 
4492  deltacands[i] = absmksetcoef;
4493  }
4494 
4495  /* also test 1.0 and maxabsmksetcoef + 1.0 as last delta values */
4496  if( maxabsmksetcoef != -1.0 )
4497  {
4498  deltacands[ndeltacands++] = maxabsmksetcoef + 1.0;
4499  }
4500 
4501  deltacands[ndeltacands++] = 1.0;
4502 
4503  maxtestdelta = MIN(ndeltacands, maxtestdelta);
4504 
4505  /* For each delta
4506  * Calculate fractionalities f_0 := b - down(b), f_j := a'_j - down(a'_j) , and derive MIR cut
4507  * a~*x' <= down(b)
4508  * integers : a~_j = down(a'_j) , if f_j <= f_0
4509  * a~_j = down(a'_j) + (f_j - f0)/(1 - f0), if f_j > f_0
4510  * continuous: a~_j = 0 , if a'_j >= 0
4511  * a~_j = a'_j/(1 - f0) , if a'_j < 0
4512  *
4513  * Transform inequality back to a^*x <= rhs:
4514  *
4515  * (lb or ub):
4516  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, a^_j := a~_j, if lb was used in transformation
4517  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, a^_j := -a~_j, if ub was used in transformation
4518  * and move the constant terms
4519  * -a~_j * lb_j == -a^_j * lb_j, or
4520  * a~_j * ub_j == -a^_j * ub_j
4521  * to the rhs.
4522  *
4523  * (vlb or vub):
4524  * x'_j := x_j - (bl_j * zl_j + dl_j), x_j == x'_j + (bl_j * zl_j + dl_j), a'_j == a_j, a^_j := a~_j, (vlb)
4525  * x'_j := (bu_j * zu_j + du_j) - x_j, x_j == (bu_j * zu_j + du_j) - x'_j, a'_j == -a_j, a^_j := -a~_j, (vub)
4526  * move the constant terms
4527  * -a~_j * dl_j == -a^_j * dl_j, or
4528  * a~_j * du_j == -a^_j * du_j
4529  * to the rhs, and update the VB variable coefficients:
4530  * a^_{zl_j} := a^_{zl_j} - a~_j * bl_j == a^_{zl_j} - a^_j * bl_j, or
4531  * a^_{zu_j} := a^_{zu_j} + a~_j * bu_j == a^_{zu_j} - a^_j * bu_j
4532  */
4533 
4534  ntmpcoefs = 0;
4535  for( i = intstart; i < mksetnnz; ++i )
4536  {
4537  SCIP_VAR* var;
4538  SCIP_Real solval;
4539  SCIP_Real QUAD(coef);
4540 
4541  var = vars[mksetinds[i]];
4542 
4543  /* get the soltion value of the continuous variable */
4544  solval = SCIPgetSolVal(scip, sol, var);
4545 
4546  /* now compute the solution value in the transform space considering complementation */
4547  if( boundtype[i] == -1 )
4548  {
4549  /* variable was complemented with global (simple) bound */
4550  if( varsign[i] == -1 )
4551  solval = SCIPvarGetUbGlobal(var) - solval;
4552  else
4553  solval = solval - SCIPvarGetLbGlobal(var);
4554  }
4555  else
4556  {
4557  assert(boundtype[i] == -2);
4558 
4559  /* variable was complemented with local (simple) bound */
4560  if( varsign[i] == -1 )
4561  solval = SCIPvarGetUbLocal(var) - solval;
4562  else
4563  solval = solval - SCIPvarGetLbLocal(var);
4564  }
4565 
4566  tmpvalues[ntmpcoefs] = solval;
4567  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4568  tmpcoefs[ntmpcoefs] = varsign[i] * QUAD_TO_DBL(coef);
4569  ++ntmpcoefs;
4570  }
4571 
4572  assert(ntmpcoefs == mksetnnz - intstart);
4573 
4574  contactivity = 0.0;
4575  contsqrnorm = 0.0;
4576  for( i = 0; i < intstart; ++i )
4577  {
4578  SCIP_Real solval;
4579  SCIP_Real QUAD(mksetcoef);
4580 
4581  QUAD_ARRAY_LOAD(mksetcoef, mksetcoefs, mksetinds[i]);
4582 
4583  if( varsign[i] * QUAD_TO_DBL(mksetcoef) >= 0.0 )
4584  continue;
4585 
4586  /* get the soltion value of the continuous variable */
4587  solval = SCIPgetSolVal(scip, sol, vars[mksetinds[i]]);
4588 
4589  /* now compute the solution value in the transform space considering complementation */
4590  switch( boundtype[i] )
4591  {
4592  case -1:
4593  /* variable was complemented with global (simple) bound */
4594  if( varsign[i] == -1 )
4595  solval = SCIPvarGetUbGlobal(vars[mksetinds[i]]) - solval;
4596  else
4597  solval = solval - SCIPvarGetLbGlobal(vars[mksetinds[i]]);
4598  break;
4599  case -2:
4600  /* variable was complemented with local (simple) bound */
4601  if( varsign[i] == -1 )
4602  solval = SCIPvarGetUbLocal(vars[mksetinds[i]]) - solval;
4603  else
4604  solval = solval - SCIPvarGetLbLocal(vars[mksetinds[i]]);
4605  break;
4606  default:
4607  /* variable was complemented with a variable bound */
4608  if( varsign[i] == -1 )
4609  {
4610  SCIP_Real coef;
4611  SCIP_Real constant;
4612  SCIP_Real vbdsolval;
4613 
4614  coef = SCIPvarGetVubCoefs(vars[mksetinds[i]])[boundtype[i]];
4615  constant = SCIPvarGetVubConstants(vars[mksetinds[i]])[boundtype[i]];
4616  vbdsolval = SCIPgetSolVal(scip, sol, SCIPvarGetVubVars(vars[mksetinds[i]])[boundtype[i]]);
4617 
4618  solval = (coef * vbdsolval + constant) - solval;
4619  }
4620  else
4621  {
4622  SCIP_Real coef;
4623  SCIP_Real constant;
4624  SCIP_Real vbdsolval;
4625 
4626  coef = SCIPvarGetVlbCoefs(vars[mksetinds[i]])[boundtype[i]];
4627  constant = SCIPvarGetVlbConstants(vars[mksetinds[i]])[boundtype[i]];
4628  vbdsolval = SCIPgetSolVal(scip, sol, SCIPvarGetVlbVars(vars[mksetinds[i]])[boundtype[i]]);
4629 
4630  solval = solval - (coef * vbdsolval + constant);
4631  }
4632  }
4633 
4634  contactivity += solval * (QUAD_TO_DBL(mksetcoef) * varsign[i]);
4635  contsqrnorm += QUAD_TO_DBL(mksetcoef) * QUAD_TO_DBL(mksetcoef);
4636  }
4637 
4638  {
4639  SCIP_ROW** rows;
4640 
4641  rows = SCIPgetLPRows(scip);
4642 
4643  for( i = 0; i < aggrrow->nrows; ++i )
4644  {
4645  SCIP_ROW* row;
4646  SCIP_Real slackval;
4647 
4648  row = rows[aggrrow->rowsinds[i]];
4649 
4650  if( (aggrrow->rowweights[i] * aggrrow->slacksign[i]) >= 0.0 && !row->integral )
4651  continue;
4652 
4653  /* compute solution value of slack variable */
4654  slackval = SCIPgetRowSolActivity(scip, row, sol);
4655 
4656  if( aggrrow->slacksign[i] == +1 )
4657  {
4658  /* right hand side */
4659  assert(!SCIPisInfinity(scip, row->rhs));
4660 
4661  slackval = row->rhs - slackval;
4662  }
4663  else
4664  {
4665  /* left hand side */
4666  assert(aggrrow->slacksign[i] == -1);
4667  assert(!SCIPisInfinity(scip, -row->lhs));
4668 
4669  slackval = slackval - row->lhs;
4670  }
4671 
4672  if( row->integral )
4673  {
4674  /* if row is integral add variable to tmp arrays */
4675  tmpvalues[ntmpcoefs] = slackval;
4676  tmpcoefs[ntmpcoefs] = aggrrow->rowweights[i] * aggrrow->slacksign[i];
4677  ++ntmpcoefs;
4678  }
4679  else
4680  {
4681  SCIP_Real slackcoeff = (aggrrow->rowweights[i] * aggrrow->slacksign[i]);
4682 
4683  /* otherwise add it to continuous activity */
4684  contactivity += slackval * slackcoeff;
4685  contsqrnorm += SQR(slackcoeff);
4686  }
4687  }
4688  }
4689 
4690  /* try all candidates for delta and remember best */
4691  bestdelta = SCIP_INVALID;
4692  bestefficacy = -SCIPinfinity(scip);
4693 
4694  for( i = 0; i < maxtestdelta; ++i )
4695  {
4696  int j;
4697  SCIP_Real efficacy;
4698 
4699  /* check if we have seen this value of delta before */
4700  SCIP_Bool deltaseenbefore = FALSE;
4701  for( j = 0; j < i; ++j )
4702  {
4703  if( SCIPisEQ(scip, deltacands[i], deltacands[j]) )
4704  {
4705  deltaseenbefore = TRUE;
4706  break;
4707  }
4708  }
4709 
4710  /* skip this delta value and allow one more delta value if available */
4711  if( deltaseenbefore )
4712  {
4713  maxtestdelta = MIN(maxtestdelta + 1, ndeltacands);
4714  continue;
4715  }
4716 
4717  efficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(mksetrhs), contactivity, contsqrnorm, deltacands[i], ntmpcoefs, minfrac, maxfrac);
4718 
4719  if( efficacy > bestefficacy )
4720  {
4721  bestefficacy = efficacy;
4722  bestdelta = deltacands[i];
4723  }
4724  }
4725 
4726  /* no delta was found that yielded any cut */
4727  if( bestdelta == SCIP_INVALID ) /*lint !e777*/
4728  goto TERMINATE;
4729 
4730  /* try bestdelta divided by 2, 4 and 8 */
4731  {
4732  SCIP_Real basedelta = bestdelta;
4733  for( i = 2; i <= 8 ; i *= 2 )
4734  {
4735  SCIP_Real efficacy;
4736  SCIP_Real delta;
4737 
4738  delta = basedelta / i;
4739 
4740  efficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(mksetrhs), contactivity, contsqrnorm, delta, ntmpcoefs, minfrac, maxfrac);
4741 
4742  if( efficacy > bestefficacy )
4743  {
4744  bestefficacy = efficacy;
4745  bestdelta = delta;
4746  }
4747  }
4748  }
4749 
4750  /* try to improve efficacy by switching complementation of integral variables that are not at their bounds
4751  * in order of non-increasing bound distance
4752  */
4753  for( i = 0; i < nbounddist; ++i )
4754  {
4755  int k;
4756  SCIP_Real newefficacy;
4757  SCIP_Real QUAD(newrhs);
4758  SCIP_Real bestlb;
4759  SCIP_Real bestub;
4760  SCIP_Real oldsolval;
4761  SCIP_Real simplebnd;
4762  int bestlbtype;
4763  int bestubtype;
4764 
4765  k = bounddistpos[i];
4766 
4767  SCIP_CALL( findBestLb(scip, vars[mksetinds[k]], sol, FALSE, allowlocal, &bestlb, &simplebnd, &bestlbtype) );
4768 
4769  if( SCIPisInfinity(scip, -bestlb) )
4770  continue;
4771 
4772  SCIP_CALL( findBestUb(scip, vars[mksetinds[k]], sol, FALSE, allowlocal, &bestub, &simplebnd, &bestubtype) );
4773 
4774  if( SCIPisInfinity(scip, bestub) )
4775  continue;
4776 
4777  /* switch the complementation of this variable */
4778 #ifndef NDEBUG
4779  {
4780  SCIP_Real QUAD(coef);
4781  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[k]);
4782  assert(SCIPisEQ(scip, tmpcoefs[k - intstart], varsign[k] * QUAD_TO_DBL(coef)));
4783  }
4784 #endif
4785 
4786  /* compute this: newrhs = mksetrhs + tmpcoefs[k - intstart] * (bestlb - bestub); */
4787  SCIPquadprecSumQD(newrhs, mksetrhs, tmpcoefs[k - intstart] * (bestlb - bestub));
4788  tmpcoefs[k - intstart] = -tmpcoefs[k - intstart];
4789 
4790  oldsolval = tmpvalues[k - intstart];
4791  tmpvalues[k - intstart] = varsign[k] == +1 ? bestub - SCIPgetSolVal(scip, sol, vars[mksetinds[k]]) : SCIPgetSolVal(scip, sol, vars[mksetinds[k]]) - bestlb;
4792 
4793  /* compute new violation */
4794  newefficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(newrhs), contactivity, contsqrnorm, bestdelta, ntmpcoefs, minfrac, maxfrac);
4795 
4796  /* check if violaton was increased */
4797  if( newefficacy > bestefficacy )
4798  {
4799  /* keep change of complementation */
4800  bestefficacy = newefficacy;
4801  QUAD_ASSIGN_Q(mksetrhs, newrhs);
4802 
4803  if( varsign[k] == +1 )
4804  {
4805  /* switch to upper bound */
4806  assert(bestubtype < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
4807  boundtype[k] = bestubtype;
4808  varsign[k] = -1;
4809  }
4810  else
4811  {
4812  /* switch to lower bound */
4813  assert(bestlbtype < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
4814  boundtype[k] = bestlbtype;
4815  varsign[k] = +1;
4816  }
4817 
4818  localbdsused = localbdsused || (boundtype[k] == -2);
4819  }
4820  else
4821  {
4822  /* undo the change of the complementation */
4823  tmpcoefs[k - intstart] = -tmpcoefs[k - intstart];
4824  tmpvalues[k - intstart] = oldsolval;
4825  }
4826  } /*lint !e438*/
4827 
4828  if( bestefficacy > 0.0 )
4829  {
4830  SCIP_Real mirefficacy;
4831  SCIP_Real QUAD(downrhs);
4832  SCIP_Real QUAD(f0);
4833  SCIP_Real scale;
4834 
4835  scale = 1.0 / bestdelta;
4836  SCIPquadprecProdQD(mksetrhs, mksetrhs, scale);
4837 
4838  SCIPquadprecEpsFloorQ(downrhs, mksetrhs, SCIPepsilon(scip)); /*lint !e666*/
4839  SCIPquadprecSumQQ(f0, mksetrhs, -downrhs);
4840 
4841  /* renormaliize f0 value */
4842  SCIPquadprecSumDD(f0, QUAD_HI(f0), QUAD_LO(f0));
4843 
4844  for( i = 0; i < mksetnnz; ++i )
4845  {
4846  SCIP_Real QUAD(coef);
4847 
4848  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4849  SCIPquadprecProdQD(coef, coef, scale);
4850  QUAD_ARRAY_STORE(mksetcoefs, mksetinds[i], coef);
4851  }
4852  SCIPdebugMessage("applied best scale (=%.13g):\n", scale);
4853  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4854 
4855  QUAD_ASSIGN_Q(mksetrhs, downrhs);
4856 
4857  SCIP_CALL( cutsRoundMIR(scip, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, varsign, boundtype, QUAD(f0)) );
4858 
4859  SCIPdebugMessage("rounded MIR cut:\n");
4860  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4861 
4862  /* substitute aggregated slack variables:
4863  *
4864  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
4865  * variable only appears in its own row:
4866  * a'_r = scale * weight[r] * slacksign[r].
4867  *
4868  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
4869  * integers : a^_r = a~_r = down(a'_r) , if f_r <= f0
4870  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
4871  * continuous: a^_r = a~_r = 0 , if a'_r >= 0
4872  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
4873  *
4874  * Substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
4875  */
4876  SCIP_CALL( cutsSubstituteMIR(scip, aggrrow->rowweights, aggrrow->slacksign, aggrrow->rowsinds,
4877  aggrrow->nrows, scale, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, QUAD(f0)) );
4878 
4879  SCIPdebugMessage("substituted slacks in MIR cut:\n");
4880  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4881 
4882 #ifndef NDEBUG
4883  {
4884  SCIP_Real efficacy = -QUAD_TO_DBL(mksetrhs);
4885  for( i = 0; i < mksetnnz; ++i )
4886  {
4887  SCIP_Real QUAD(coef);
4888  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4889  efficacy += QUAD_TO_DBL(coef) * SCIPgetSolVal(scip, sol, vars[mksetinds[i]]);
4890  }
4891 
4892  if( !EPSZ(SCIPrelDiff(efficacy, bestefficacy), 1e-4) )
4893  {
4894  SCIPdebugMessage("efficacy of cmir cut is different than expected efficacy: %f != %f\n", efficacy, bestefficacy);
4895  }
4896  }
4897 #endif
4898 
4899  *cutislocal = *cutislocal || localbdsused;
4900 
4901  /* remove all nearly-zero coefficients from MIR row and relax the right hand side correspondingly in order to
4902  * prevent numerical rounding errors
4903  */
4904  if( postprocess )
4905  {
4906  SCIP_CALL( postprocessCutQuad(scip, *cutislocal, mksetinds, mksetcoefs, &mksetnnz, QUAD(&mksetrhs), success) );
4907  }
4908  else
4909  {
4910  *success = ! removeZerosQuad(scip, SCIPsumepsilon(scip), *cutislocal, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz);
4911  }
4912 
4913  SCIPdebugMessage("post-processed cut (success = %s):\n", *success ? "TRUE" : "FALSE");
4914  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4915 
4916  if( *success )
4917  {
4918  mirefficacy = calcEfficacyDenseStorageQuad(scip, sol, mksetcoefs, QUAD_TO_DBL(mksetrhs), mksetinds, mksetnnz);
4919 
4920  if( SCIPisEfficacious(scip, mirefficacy) && mirefficacy > *cutefficacy )
4921  {
4922  BMScopyMemoryArray(cutinds, mksetinds, mksetnnz);
4923  for( i = 0; i < mksetnnz; ++i )
4924  {
4925  SCIP_Real QUAD(coef);
4926  int j = cutinds[i];
4927 
4928  QUAD_ARRAY_LOAD(coef, mksetcoefs, j);
4929 
4930  cutcoefs[i] = QUAD_TO_DBL(coef);
4931  QUAD_ASSIGN(coef, 0.0);
4932  QUAD_ARRAY_STORE(mksetcoefs, j, coef);
4933  }
4934  *cutnnz = mksetnnz;
4935  *cutrhs = QUAD_TO_DBL(mksetrhs);
4936  *cutefficacy = mirefficacy;
4937  if( cutrank != NULL )
4938  *cutrank = aggrrow->rank + 1;
4939  *cutislocal = *cutislocal || localbdsused;
4940  }
4941  else
4942  {
4943  *success = FALSE;
4944  }
4945  }
4946  }
4947 
4948  TERMINATE:
4949  /* if we aborted early we need to clean the mksetcoefs */
4950  if( !(*success) )
4951  {
4952  SCIP_Real QUAD(tmp);
4953  QUAD_ASSIGN(tmp, 0.0);
4954 
4955  for( i = 0; i < mksetnnz; ++i )
4956  {
4957  QUAD_ARRAY_STORE(mksetcoefs, mksetinds[i], tmp);
4958  }
4959  }
4960 
4961  /* free temporary memory */
4962  SCIPfreeBufferArray(scip, &bounddistpos);
4963  SCIPfreeBufferArray(scip, &bounddist);
4964  SCIPfreeBufferArray(scip, &deltacands);
4965  SCIPfreeBufferArray(scip, &tmpvalues);
4966  SCIPfreeBufferArray(scip, &tmpcoefs);
4967  SCIPfreeBufferArray(scip, &mksetinds);
4968  SCIPfreeCleanBufferArray(scip, &mksetcoefs);
4969  SCIPfreeBufferArray(scip, &boundtype);
4970  SCIPfreeBufferArray(scip, &varsign);
4971 
4972  return SCIP_OKAY;
4973 }
4974 
4975 /* =========================================== flow cover =========================================== */
4976 
4977 #define NO_EXACT_KNAPSACK
4978 
4979 #ifndef NO_EXACT_KNAPSACK
4980 #define MAXDNOM 1000LL
4981 #define MINDELTA 1e-03
4982 #define MAXDELTA 1e-09
4983 #define MAXSCALE 1000.0
4984 #define MAXDYNPROGSPACE 1000000
4985 #endif
4986 
4987 #define MAXABSVBCOEF 1e+5 /**< maximal absolute coefficient in variable bounds used for snf relaxation */
4988 #define MAXBOUND 1e+10 /**< maximal value of normal bounds used for snf relaxation */
4989 
4990 /** structure that contains all data required to perform the sequence independent lifting
4991  */
4992 typedef
4993 struct LiftingData
4994 {
4995  SCIP_Real* M; /**< \f$ M_0 := 0.0 \f$ and \f$ M_i := M_i-1 + m_i \f$ */
4996  SCIP_Real* m; /**< non-increasing array of variable upper bound coefficients
4997  * for all variables in \f$ C^{++} \f$ and \f$ L^- \f$,
4998  * where \f$ C = C^+ \cup C^- \f$ is the flowcover and
4999  * \f$ C^{++} := \{ j \in C^+ \mid u_j > \lambda \} \f$
5000  * \f$ L^- := \{ j \in (N^- \setminus C^-) \mid u_j > \lambda \} \f$
5001  */
5002  int r; /**< size of array m */
5003  int t; /**< index of smallest value in m that comes from a variable in \f$ C^{++} \f$ */
5004  SCIP_Real d1; /**< right hand side of single-node-flow set plus the sum of all \f$ u_j \f$ for \f$ j \in C^- \f$ */
5005  SCIP_Real d2; /**< right hand side of single-node-flow set plus the sum of all \f$ u_j \f$ for \f$ j \in N^- \f$ */
5006  SCIP_Real lambda; /**< excess of the flowcover */
5007  SCIP_Real mp; /**< smallest variable bound coefficient of variable in \f$ C^{++} (min_{j \in C++} u_j) \f$ */
5008  SCIP_Real ml; /**< \f$ ml := min(\lambda, \sum_{j \in C^+ \setminus C^{++}} u_j) \f$ */
5009 } LIFTINGDATA;
5010 
5011 /** structure that contains all the data that defines the single-node-flow relaxation of an aggregation row */
5012 typedef
5013 struct SNF_Relaxation
5014 {
5015  int* transvarcoefs; /**< coefficients of all vars in relaxed set */
5016  SCIP_Real* transbinvarsolvals; /**< sol val of bin var in vub of all vars in relaxed set */
5017  SCIP_Real* transcontvarsolvals;/**< sol val of all real vars in relaxed set */
5018  SCIP_Real* transvarvubcoefs; /**< coefficient in vub of all vars in relaxed set */
5019  int ntransvars; /**< number of vars in relaxed set */
5020  SCIP_Real transrhs; /**< rhs in relaxed set */
5021  int* origbinvars; /**< associated original binary var for all vars in relaxed set */
5022  int* origcontvars; /**< associated original continuous var for all vars in relaxed set */
5023  SCIP_Real* aggrcoefsbin; /**< aggregation coefficient of the original binary var used to define the
5024  * continuous variable in the relaxed set */
5025  SCIP_Real* aggrcoefscont; /**< aggregation coefficient of the original continous var used to define the
5026  * continuous variable in the relaxed set */
5027  SCIP_Real* aggrconstants; /**< aggregation constant used to define the continuous variable in the relaxed set */
5028 } SNF_RELAXATION;
5029 
5030 /** get solution value and index of variable lower bound (with binary variable) which is closest to the current LP
5031  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
5032  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
5033  * given variable
5034  */
5035 static
5037  SCIP* scip, /**< SCIP data structure */
5038  SCIP_VAR* var, /**< given active problem variable */
5039  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
5040  SCIP_Real* rowcoefs, /**< (dense) array of coefficients of row */
5041  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
5042  * was not used (0) or was not used but is contained in the row (-1)
5043  */
5044  SCIP_Real bestsub, /**< closest simple upper bound of given variable */
5045  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
5046  SCIP_Real* closestvlb, /**< pointer to store the LP sol value of the closest variable lower bound */
5047  int* closestvlbidx /**< pointer to store the index of the closest vlb; -1 if no vlb was found */
5048  )
5049 {
5050  int nvlbs;
5051  int nbinvars;
5052 
5053  assert(scip != NULL);
5054  assert(var != NULL);
5055  assert(bestsub == SCIPvarGetUbGlobal(var) || bestsub == SCIPvarGetUbLocal(var)); /*lint !e777*/
5056  assert(!SCIPisInfinity(scip, bestsub));
5057  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5058  assert(rowcoefs != NULL);
5059  assert(binvarused != NULL);
5060  assert(closestvlb != NULL);
5061  assert(closestvlbidx != NULL);
5062 
5063  nvlbs = SCIPvarGetNVlbs(var);
5064  nbinvars = SCIPgetNBinVars(scip);
5065 
5066  *closestvlbidx = -1;
5067  *closestvlb = -SCIPinfinity(scip);
5068  if( nvlbs > 0 )
5069  {
5070  SCIP_VAR** vlbvars;
5071  SCIP_Real* vlbcoefs;
5072  SCIP_Real* vlbconsts;
5073  int i;
5074 
5075  vlbvars = SCIPvarGetVlbVars(var);
5076  vlbcoefs = SCIPvarGetVlbCoefs(var);
5077  vlbconsts = SCIPvarGetVlbConstants(var);
5078 
5079  for( i = 0; i < nvlbs; i++ )
5080  {
5081  SCIP_Real rowcoefbinvar;
5082  SCIP_Real val1;
5083  SCIP_Real val2;
5084  SCIP_Real vlbsol;
5085  SCIP_Real rowcoefsign;
5086  int probidxbinvar;
5087 
5088  if( bestsub > vlbconsts[i] )
5089  continue;
5090 
5091  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
5092  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
5093  */
5094  if( REALABS(vlbcoefs[i]) > MAXABSVBCOEF )
5095  continue;
5096 
5097  /* use only variable lower bounds l~_i * x_i + d_i with x_i binary which are active */
5098  probidxbinvar = SCIPvarGetProbindex(vlbvars[i]);
5099 
5100  /* if the variable is not active the problem index is -1, so we cast to unsigned int before the comparison which
5101  * ensures that the problem index is between 0 and nbinvars - 1
5102  */
5103  if( (unsigned int)probidxbinvar >= (unsigned int)nbinvars )
5104  continue;
5105 
5106  assert(SCIPvarIsBinary(vlbvars[i]));
5107 
5108  /* check if current variable lower bound l~_i * x_i + d_i imposed on y_j meets the following criteria:
5109  * (let a_j = coefficient of y_j in current row,
5110  * u_j = closest simple upper bound imposed on y_j,
5111  * c_i = coefficient of x_i in current row)
5112  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k yet
5113  * if a_j > 0:
5114  * 1. u_j <= d_i
5115  * 2. a_j ( u_j - d_i ) + c_i <= 0
5116  * 3. a_j l~_i + c_i <= 0
5117  * if a_j < 0:
5118  * 1. u_j <= d_i
5119  * 2. a_j ( u_j - d_i ) + c_i >= 0
5120  * 3. a_j l~_i + c_i >= 0
5121  */
5122 
5123  /* has already been used in the SNF relaxation */
5124  if( binvarused[probidxbinvar] == 1 )
5125  continue;
5126 
5127  /* get the row coefficient */
5128  {
5129  SCIP_Real QUAD(tmp);
5130  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidxbinvar);
5131  rowcoefbinvar = QUAD_TO_DBL(tmp);
5132  }
5133  rowcoefsign = COPYSIGN(1.0, rowcoef);
5134 
5135  val2 = rowcoefsign * ((rowcoef * vlbcoefs[i]) + rowcoefbinvar);
5136 
5137  /* variable lower bound does not meet criteria */
5138  if( val2 > 0.0 || SCIPisInfinity(scip, -val2) )
5139  continue;
5140 
5141  val1 = rowcoefsign * ((rowcoef * (bestsub - vlbconsts[i])) + rowcoefbinvar);
5142 
5143  /* variable lower bound does not meet criteria */
5144  if( val1 > 0.0 )
5145  continue;
5146 
5147  vlbsol = vlbcoefs[i] * SCIPgetSolVal(scip, sol, vlbvars[i]) + vlbconsts[i];
5148  if( vlbsol > *closestvlb )
5149  {
5150  *closestvlb = vlbsol;
5151  *closestvlbidx = i;
5152  }
5153  assert(*closestvlbidx >= 0);
5154  }
5155  }
5156 
5157  return SCIP_OKAY;
5158 }
5159 
5160 /** get LP solution value and index of variable upper bound (with binary variable) which is closest to the current LP
5161  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
5162  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
5163  * given variable
5164  */
5165 static
5167  SCIP* scip, /**< SCIP data structure */
5168  SCIP_VAR* var, /**< given active problem variable */
5169  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
5170  SCIP_Real* rowcoefs, /**< (dense) array of coefficients of row */
5171  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
5172  * was not used (0) or was not used but is contained in the row (-1)
5173  */
5174  SCIP_Real bestslb, /**< closest simple lower bound of given variable */
5175  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
5176  SCIP_Real* closestvub, /**< pointer to store the LP sol value of the closest variable upper bound */
5177  int* closestvubidx /**< pointer to store the index of the closest vub; -1 if no vub was found */
5178  )
5179 {
5180  int nvubs;
5181  int nbinvars;
5182 
5183  assert(scip != NULL);
5184  assert(var != NULL);
5185  assert(bestslb == SCIPvarGetLbGlobal(var) || bestslb == SCIPvarGetLbLocal(var)); /*lint !e777*/
5186  assert(!SCIPisInfinity(scip, - bestslb));
5187  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5188  assert(rowcoefs != NULL);
5189  assert(binvarused != NULL);
5190  assert(closestvub != NULL);
5191  assert(closestvubidx != NULL);
5192 
5193  nvubs = SCIPvarGetNVubs(var);
5194  nbinvars = SCIPgetNBinVars(scip);
5195 
5196  *closestvubidx = -1;
5197  *closestvub = SCIPinfinity(scip);
5198  if( nvubs > 0 )
5199  {
5200  SCIP_VAR** vubvars;
5201  SCIP_Real* vubcoefs;
5202  SCIP_Real* vubconsts;
5203  int i;
5204 
5205  vubvars = SCIPvarGetVubVars(var);
5206  vubcoefs = SCIPvarGetVubCoefs(var);
5207  vubconsts = SCIPvarGetVubConstants(var);
5208 
5209  for( i = 0; i < nvubs; i++ )
5210  {
5211  SCIP_Real rowcoefbinvar;
5212  SCIP_Real val1;
5213  SCIP_Real val2;
5214  SCIP_Real vubsol;
5215  SCIP_Real rowcoefsign;
5216  int probidxbinvar;
5217 
5218  if( bestslb < vubconsts[i] )
5219  continue;
5220 
5221  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
5222  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
5223  */
5224  if( REALABS(vubcoefs[i]) > MAXABSVBCOEF )
5225  continue;
5226 
5227  /* use only variable upper bound u~_i * x_i + d_i with x_i binary and which are active */
5228  probidxbinvar = SCIPvarGetProbindex(vubvars[i]);
5229 
5230  /* if the variable is not active the problem index is -1, so we cast to unsigned int before the comparison which
5231  * ensures that the problem index is between 0 and nbinvars - 1
5232  */
5233  if( (unsigned int)probidxbinvar >= (unsigned int)nbinvars )
5234  continue;
5235 
5236  assert(SCIPvarIsBinary(vubvars[i]));
5237 
5238  /* checks if current variable upper bound u~_i * x_i + d_i meets the following criteria
5239  * (let a_j = coefficient of y_j in current row,
5240  * l_j = closest simple lower bound imposed on y_j,
5241  * c_i = coefficient of x_i in current row)
5242  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k
5243  * if a > 0:
5244  * 1. l_j >= d_i
5245  * 2. a_j ( l_i - d_i ) + c_i >= 0
5246  * 3. a_j u~_i + c_i >= 0
5247  * if a < 0:
5248  * 1. l_j >= d_i
5249  * 2. a_j ( l_j - d_i ) + c_i <= 0
5250  * 3. a_j u~_i + c_i <= 0
5251  */
5252 
5253  /* has already been used in the SNF relaxation */
5254  if( binvarused[probidxbinvar] == 1 )
5255  continue;
5256 
5257  /* get the row coefficient */
5258  {
5259  SCIP_Real QUAD(tmp);
5260  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidxbinvar);
5261  rowcoefbinvar = QUAD_TO_DBL(tmp);
5262  }
5263  rowcoefsign = COPYSIGN(1.0, rowcoef);
5264 
5265  val2 = rowcoefsign * ((rowcoef * vubcoefs[i]) + rowcoefbinvar);
5266 
5267  /* variable upper bound does not meet criteria */
5268  if( val2 < 0.0 || SCIPisInfinity(scip, val2) )
5269  continue;
5270 
5271  val1 = rowcoefsign * ((rowcoef * (bestslb - vubconsts[i])) + rowcoefbinvar);
5272 
5273  /* variable upper bound does not meet criteria */
5274  if( val1 < 0.0 )
5275  continue;
5276 
5277  vubsol = vubcoefs[i] * SCIPgetSolVal(scip, sol, vubvars[i]) + vubconsts[i];
5278  if( vubsol < *closestvub )
5279  {
5280  *closestvub = vubsol;
5281  *closestvubidx = i;
5282  }
5283  assert(*closestvubidx >= 0);
5284  }
5285  }
5286 
5287  return SCIP_OKAY;
5288 }
5289 
5290 /** determines the bounds to use for constructing the single-node-flow relaxation of a variable in
5291  * the given row.
5292  */
5293 static
5295  SCIP* scip, /**< SCIP data structure */
5296  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
5297  SCIP_VAR** vars, /**< array of problem variables */
5298  SCIP_Real* rowcoefs, /**< (dense) array of variable coefficients in the row */
5299  int* rowinds, /**< array with positions of non-zero values in the rowcoefs array */
5300  int varposinrow, /**< position of variable in the rowinds array for which the bounds should be determined */
5301  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
5302  * was not used (0) or was not used but is contained in the row (-1)
5303  */
5304  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
5305  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
5306  SCIP_Real* bestlb, /**< pointer to store best lower bound for transformation */
5307  SCIP_Real* bestub, /**< pointer to store best upper bound for transformation */
5308  SCIP_Real* bestslb, /**< pointer to store best simple lower bound for transformation */
5309  SCIP_Real* bestsub, /**< pointer to store best simple upper bound for transformation */
5310  int* bestlbtype, /**< pointer to store type of best lower bound */
5311  int* bestubtype, /**< pointer to store type of best upper bound */
5312  int* bestslbtype, /**< pointer to store type of best simple lower bound */
5313  int* bestsubtype, /**< pointer to store type of best simple upper bound */
5314  SCIP_BOUNDTYPE* selectedbounds, /**< pointer to store the preferred bound for the transformation */
5315  SCIP_Bool* freevariable /**< pointer to store if variable is a free variable */
5316  )
5317 {
5318  SCIP_VAR* var;
5319 
5320  SCIP_Real rowcoef;
5321  SCIP_Real solval;
5322  SCIP_Real simplebound;
5323 
5324  int probidx;
5325 
5326  bestlb[varposinrow] = -SCIPinfinity(scip);
5327  bestub[varposinrow] = SCIPinfinity(scip);
5328  bestlbtype[varposinrow] = -3;
5329  bestubtype[varposinrow] = -3;
5330 
5331  probidx = rowinds[varposinrow];
5332  var = vars[probidx];
5333  {
5334  SCIP_Real QUAD(tmp);
5335  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidx);
5336  rowcoef = QUAD_TO_DBL(tmp);
5337  }
5338 
5339  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5340 
5341  /* get closest simple lower bound and closest simple upper bound */
5342  SCIP_CALL( findBestLb(scip, var, sol, FALSE, allowlocal, &bestslb[varposinrow], &simplebound, &bestslbtype[varposinrow]) );
5343  SCIP_CALL( findBestUb(scip, var, sol, FALSE, allowlocal, &bestsub[varposinrow], &simplebound, &bestsubtype[varposinrow]) );
5344 
5345  /* do not use too large bounds */
5346  if( bestslb[varposinrow] <= -MAXBOUND )
5347  bestslb[varposinrow] = -SCIPinfinity(scip);
5348 
5349  if( bestsub[varposinrow] >= MAXBOUND )
5350  bestsub[varposinrow] = SCIPinfinity(scip);
5351 
5352  solval = SCIPgetSolVal(scip, sol, var);
5353 
5354  SCIPdebugMsg(scip, " %d: %g <%s, idx=%d, lp=%g, [%g(%d),%g(%d)]>:\n", varposinrow, rowcoef, SCIPvarGetName(var), probidx,
5355  solval, bestslb[varposinrow], bestslbtype[varposinrow], bestsub[varposinrow], bestsubtype[varposinrow]);
5356 
5357  /* mixed integer set cannot be relaxed to 0-1 single node flow set because both simple bounds are -infinity
5358  * and infinity, respectively
5359  */
5360  if( SCIPisInfinity(scip, -bestslb[varposinrow]) && SCIPisInfinity(scip, bestsub[varposinrow]) )
5361  {
5362  *freevariable = TRUE;
5363  return SCIP_OKAY;
5364  }
5365 
5366  /* get closest lower bound that can be used to define the real variable y'_j in the 0-1 single node flow
5367  * relaxation
5368  */
5369  if( !SCIPisInfinity(scip, bestsub[varposinrow]) )
5370  {
5371  bestlb[varposinrow] = bestslb[varposinrow];
5372  bestlbtype[varposinrow] = bestslbtype[varposinrow];
5373 
5375  {
5376  SCIP_Real bestvlb;
5377  int bestvlbidx;
5378 
5379  SCIP_CALL( getClosestVlb(scip, var, sol, rowcoefs, binvarused, bestsub[varposinrow], rowcoef, &bestvlb, &bestvlbidx) );
5380  if( SCIPisGT(scip, bestvlb, bestlb[varposinrow]) )
5381  {
5382  bestlb[varposinrow] = bestvlb;
5383  bestlbtype[varposinrow] = bestvlbidx;
5384  }
5385  }
5386  }
5387 
5388  /* get closest upper bound that can be used to define the real variable y'_j in the 0-1 single node flow
5389  * relaxation
5390  */
5391  if( !SCIPisInfinity(scip, -bestslb[varposinrow]) )
5392  {
5393  bestub[varposinrow] = bestsub[varposinrow];
5394  bestubtype[varposinrow] = bestsubtype[varposinrow];
5395 
5397  {
5398  SCIP_Real bestvub;
5399  int bestvubidx;
5400 
5401  SCIP_CALL( getClosestVub(scip, var, sol, rowcoefs, binvarused, bestslb[varposinrow], rowcoef, &bestvub, &bestvubidx) );
5402  if( SCIPisLT(scip, bestvub, bestub[varposinrow]) )
5403  {
5404  bestub[varposinrow] = bestvub;
5405  bestubtype[varposinrow] = bestvubidx;
5406  }
5407  }
5408  }
5409  SCIPdebugMsg(scip, " bestlb=%g(%d), bestub=%g(%d)\n", bestlb[varposinrow], bestlbtype[varposinrow], bestub[varposinrow], bestubtype[varposinrow]);
5410 
5411  /* mixed integer set cannot be relaxed to 0-1 single node flow set because there are no suitable bounds
5412  * to define the transformed variable y'_j
5413  */
5414  if( SCIPisInfinity(scip, -bestlb[varposinrow]) && SCIPisInfinity(scip, bestub[varposinrow]) )
5415  {
5416  *freevariable = TRUE;
5417  return SCIP_OKAY;
5418  }
5419 
5420  *freevariable = FALSE;
5421 
5422  /* select best upper bound if it is closer to the LP value of y_j and best lower bound otherwise and use this bound
5423  * to define the real variable y'_j with 0 <= y'_j <= u'_j x_j in the 0-1 single node flow relaxation;
5424  * prefer variable bounds
5425  */
5426  if( SCIPisEQ(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]) && bestlbtype[varposinrow] >= 0 )
5427  {
5428  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_LOWER;
5429  }
5430  else if( SCIPisEQ(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow])
5431  && bestubtype[varposinrow] >= 0 )
5432  {
5433  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_UPPER;
5434  }
5435  else if( SCIPisLE(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]) )
5436  {
5437  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_LOWER;
5438  }
5439  else
5440  {
5441  assert(SCIPisGT(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]));
5442  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_UPPER;
5443  }
5444 
5445  if( selectedbounds[varposinrow] == SCIP_BOUNDTYPE_LOWER && bestlbtype[varposinrow] >= 0 )
5446  {
5447  int vlbvarprobidx;
5448  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(var);
5449 
5450  /* mark binary variable of vlb so that it is not used for other continuous variables
5451  * by setting it's position in the aggrrow to a negative value
5452  */
5453  vlbvarprobidx = SCIPvarGetProbindex(vlbvars[bestlbtype[varposinrow]]);
5454  binvarused[vlbvarprobidx] = 1;
5455  }
5456  else if ( selectedbounds[varposinrow] == SCIP_BOUNDTYPE_UPPER && bestubtype[varposinrow] >= 0 )
5457  {
5458  int vubvarprobidx;
5459  SCIP_VAR** vubvars = SCIPvarGetVubVars(var);
5460 
5461  /* mark binary variable of vub so that it is not used for other continuous variables
5462  * by setting it's position in the aggrrow to a negative value
5463  */
5464  vubvarprobidx = SCIPvarGetProbindex(vubvars[bestubtype[varposinrow]]);
5465  binvarused[vubvarprobidx] = 1;
5466  }
5467 
5468  return SCIP_OKAY; /*lint !e438*/
5469 }
5470 
5471 /** construct a 0-1 single node flow relaxation (with some additional simple constraints) of a mixed integer set
5472  * corresponding to the given aggrrow a * x <= rhs
5473  */
5474 static
5476  SCIP* scip, /**< SCIP data structure */
5477  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5478  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
5479  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
5480  SCIP_Real* rowcoefs, /**< array of coefficients of row */
5481  QUAD(SCIP_Real rowrhs), /**< pointer to right hand side of row */
5482  int* rowinds, /**< array of variables problem indices for non-zero coefficients in row */
5483  int nnz, /**< number of non-zeros in row */
5484  SNF_RELAXATION* snf, /**< stores the sign of the transformed variable in summation */
5485  SCIP_Bool* success, /**< stores whether the transformation was valid */
5486  SCIP_Bool* localbdsused /**< pointer to store whether local bounds were used in transformation */
5487  )
5488 {
5489  SCIP_VAR** vars;
5490  int i;
5491  int nnonbinvarsrow;
5492  int8_t* binvarused;
5493  int nbinvars;
5494  SCIP_Real QUAD(transrhs);
5495 
5496  /* arrays to store the selected bound for each non-binary variable in the row */
5497  SCIP_Real* bestlb;
5498  SCIP_Real* bestub;
5499  SCIP_Real* bestslb;
5500  SCIP_Real* bestsub;
5501  int* bestlbtype;
5502  int* bestubtype;
5503  int* bestslbtype;
5504  int* bestsubtype;
5505  SCIP_BOUNDTYPE* selectedbounds;
5506 
5507  *success = FALSE;
5508 
5509  SCIPdebugMsg(scip, "--------------------- construction of SNF relaxation ------------------------------------\n");
5510 
5511  nbinvars = SCIPgetNBinVars(scip);
5512  vars = SCIPgetVars(scip);
5513 
5514  SCIP_CALL( SCIPallocBufferArray(scip, &bestlb, nnz) );
5515  SCIP_CALL( SCIPallocBufferArray(scip, &bestub, nnz) );
5516  SCIP_CALL( SCIPallocBufferArray(scip, &bestslb, nnz) );
5517  SCIP_CALL( SCIPallocBufferArray(scip, &bestsub, nnz) );
5518  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbtype, nnz) );
5519  SCIP_CALL( SCIPallocBufferArray(scip, &bestubtype, nnz) );
5520  SCIP_CALL( SCIPallocBufferArray(scip, &bestslbtype, nnz) );
5521  SCIP_CALL( SCIPallocBufferArray(scip, &bestsubtype, nnz) );
5522  SCIP_CALL( SCIPallocBufferArray(scip, &selectedbounds, nnz) );
5523 
5524  /* sort descending to have continuous variables first */
5525  SCIPsortDownInt(rowinds, nnz);
5526 
5527  /* array to store whether a binary variable is in the row (-1) or has been used (1) due to variable bound usage */
5528  SCIP_CALL( SCIPallocCleanBufferArray(scip, &binvarused, nbinvars) );
5529 
5530  for( i = nnz - 1; i >= 0 && rowinds[i] < nbinvars; --i )
5531  {
5532  int j = rowinds[i];
5533  binvarused[j] = -1;
5534  }
5535 
5536  nnonbinvarsrow = i + 1;
5537  /* determine the bounds to use for transforming the non-binary variables */
5538  for( i = 0; i < nnonbinvarsrow; ++i )
5539  {
5540  SCIP_Bool freevariable;
5541 
5542  assert(rowinds[i] >= nbinvars);
5543 
5544  SCIP_CALL( determineBoundForSNF(scip, sol, vars, rowcoefs, rowinds, i, binvarused, allowlocal, boundswitch,
5545  bestlb, bestub, bestslb, bestsub, bestlbtype, bestubtype, bestslbtype, bestsubtype, selectedbounds, &freevariable) );
5546 
5547  if( freevariable )
5548  {
5549  int j;
5550 
5551  /* clear binvarused at indices of binary variables of row */
5552  for( j = nnz - 1; j >= nnonbinvarsrow; --j )
5553  binvarused[rowinds[j]] = 0;
5554 
5555  /* clear binvarused at indices of selected variable bounds */
5556  for( j = 0; j < i; ++j )
5557  {
5558  if( selectedbounds[j] == SCIP_BOUNDTYPE_LOWER && bestlbtype[j] >= 0 )
5559  {
5560  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(vars[rowinds[j]]);
5561  binvarused[SCIPvarGetProbindex(vlbvars[bestlbtype[j]])] = 0;
5562  }
5563  else if( selectedbounds[j] == SCIP_BOUNDTYPE_UPPER && bestubtype[j] >= 0 )
5564  {
5565  SCIP_VAR** vubvars = SCIPvarGetVubVars(vars[rowinds[j]]);
5566  binvarused[SCIPvarGetProbindex(vubvars[bestubtype[j]])] = 0;
5567  }
5568  }
5569 
5570  /* terminate */
5571  goto TERMINATE;
5572  }
5573  }
5574 
5575  *localbdsused = FALSE;
5576  QUAD_ASSIGN_Q(transrhs, rowrhs);
5577  snf->ntransvars = 0;
5578 
5579  assert(snf->transvarcoefs != NULL); /* for lint */
5580  assert(snf->transvarvubcoefs != NULL);
5581  assert(snf->transbinvarsolvals != NULL);
5582  assert(snf->transcontvarsolvals != NULL);
5583  assert(snf->aggrconstants != NULL);
5584  assert(snf->aggrcoefscont != NULL);
5585  assert(snf->origcontvars != NULL);
5586  assert(snf->origbinvars != NULL);
5587  assert(snf->aggrcoefsbin != NULL);
5588 
5589  /* transform non-binary variables */
5590  for( i = 0; i < nnonbinvarsrow; ++i )
5591  {
5592  SCIP_VAR* var;
5593  SCIP_Real QUAD(rowcoef);
5594  SCIP_Real solval;
5595  int probidx;
5596 
5597  probidx = rowinds[i];
5598  var = vars[probidx];
5599  QUAD_ARRAY_LOAD(rowcoef, rowcoefs, probidx);
5600  solval = SCIPgetSolVal(scip, sol, var);
5601 
5602  assert(probidx >= nbinvars);
5603 
5604  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
5605  {
5606  /* use bestlb to define y'_j */
5607 
5608  assert(!SCIPisInfinity(scip, bestsub[i]));
5609  assert(!SCIPisInfinity(scip, - bestlb[i]));
5610  assert(bestsubtype[i] == -1 || bestsubtype[i] == -2);
5611  assert(bestlbtype[i] > -3 && bestlbtype[i] < SCIPvarGetNVlbs(var));
5612 
5613  /* store for y_j that bestlb is the bound used to define y'_j and that y'_j is the associated real variable
5614  * in the relaxed set
5615  */
5616  snf->origcontvars[snf->ntransvars] = probidx;
5617 
5618  if( bestlbtype[i] < 0 )
5619  {
5620  SCIP_Real QUAD(val);
5621  SCIP_Real QUAD(contsolval);
5622  SCIP_Real QUAD(rowcoeftimesbestsub);
5623 
5624  /* use simple lower bound in bestlb = l_j <= y_j <= u_j = bestsub to define
5625  * 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
5626  * 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,
5627  * put j into the set
5628  * N2 if a_j > 0
5629  * N1 if a_j < 0
5630  * and update the right hand side of the constraint in the relaxation
5631  * rhs = rhs - a_j u_j
5632  */
5633  SCIPquadprecSumDD(val, bestsub[i], -bestlb[i]);
5634  SCIPquadprecProdQQ(val, val, rowcoef);
5635  SCIPquadprecSumDD(contsolval, solval, -bestsub[i]);
5636  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5637 
5638  if( bestlbtype[i] == -2 || bestsubtype[i] == -2 )
5639  *localbdsused = TRUE;
5640 
5641  SCIPquadprecProdQD(rowcoeftimesbestsub, rowcoef, bestsub[i]);
5642 
5643  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5644  snf->origbinvars[snf->ntransvars] = -1;
5645  snf->aggrcoefsbin[snf->ntransvars] = 0.0;
5646 
5647  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5648  {
5649  snf->transvarcoefs[snf->ntransvars] = - 1;
5650  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5651  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5652  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5653 
5654  /* aggregation information for y'_j */
5655  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesbestsub);
5656  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5657  }
5658  else
5659  {
5660  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5661  snf->transvarcoefs[snf->ntransvars] = 1;
5662  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5663  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5664  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5665 
5666  /* aggregation information for y'_j */
5667  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesbestsub);
5668  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5669  }
5670  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesbestsub);
5671 
5672  SCIPdebugMsg(scip, " --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
5673  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5674  snf->ntransvars, QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesbestsub), QUAD_TO_DBL(rowcoef), bestsub, QUAD_TO_DBL(transrhs));
5675  }
5676  else
5677  {
5678  SCIP_Real QUAD(rowcoefbinary);
5679  SCIP_Real varsolvalbinary;
5680  SCIP_Real QUAD(val);
5681  SCIP_Real QUAD(contsolval);
5682  SCIP_Real QUAD(rowcoeftimesvlbconst);
5683  int vlbvarprobidx;
5684 
5685  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(var);
5686  SCIP_Real* vlbconsts = SCIPvarGetVlbConstants(var);
5687  SCIP_Real* vlbcoefs = SCIPvarGetVlbCoefs(var);
5688 
5689  /* use variable lower bound in bestlb = l~_j x_j + d_j <= y_j <= u_j = bestsub to define
5690  * 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
5691  * 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,
5692  * where c_j is the coefficient of x_j in the row, put j into the set
5693  * N2 if a_j > 0
5694  * N1 if a_j < 0
5695  * and update the right hand side of the constraint in the relaxation
5696  * rhs = rhs - a_j d_j
5697  */
5698 
5699  vlbvarprobidx = SCIPvarGetProbindex(vlbvars[bestlbtype[i]]);
5700  assert(binvarused[vlbvarprobidx] == 1);
5701  assert(vlbvarprobidx < nbinvars);
5702 
5703  QUAD_ARRAY_LOAD(rowcoefbinary, rowcoefs, vlbvarprobidx);
5704  varsolvalbinary = SCIPgetSolVal(scip, sol, vlbvars[bestlbtype[i]]);
5705 
5706  SCIPquadprecProdQD(val, rowcoef, vlbcoefs[bestlbtype[i]]);
5707  SCIPquadprecSumQQ(val, val, rowcoefbinary);
5708  {
5709  SCIP_Real QUAD(tmp);
5710 
5711  SCIPquadprecProdQD(tmp, rowcoefbinary, varsolvalbinary);
5712  SCIPquadprecSumDD(contsolval, solval, - vlbconsts[bestlbtype[i]]);
5713  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5714  SCIPquadprecSumQQ(contsolval, contsolval, tmp);
5715  }
5716 
5717  SCIPquadprecProdQD(rowcoeftimesvlbconst, rowcoef, vlbconsts[bestlbtype[i]]);
5718 
5719  /* clear the binvarpos array, since the variable has been processed */
5720  binvarused[vlbvarprobidx] = 0;
5721 
5722  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5723  snf->origbinvars[snf->ntransvars] = vlbvarprobidx;
5724 
5725  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5726  {
5727  snf->transvarcoefs[snf->ntransvars] = - 1;
5728  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5729  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5730  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5731 
5732  /* aggregation information for y'_j */
5733  snf->aggrcoefsbin[snf->ntransvars] = - QUAD_TO_DBL(rowcoefbinary);
5734  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5735  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesvlbconst);
5736  }
5737  else
5738  {
5739  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5740  snf->transvarcoefs[snf->ntransvars] = 1;
5741  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5742  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5743  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5744 
5745  /* aggregation information for y'_j */
5746  snf->aggrcoefsbin[snf->ntransvars] = QUAD_TO_DBL(rowcoefbinary);
5747  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5748  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesvlbconst);
5749  }
5750  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesvlbconst);
5751 
5752  SCIPdebugMsg(scip, " --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=%s), rhs=%g-(%g*%g)=%g\n",
5753  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5754  snf->ntransvars, SCIPvarGetName(vlbvars[bestlbtype[i]]), QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesvlbconst), QUAD_TO_DBL(rowcoef),
5755  vlbconsts[bestlbtype[i]], snf->transrhs );
5756  }
5757  }
5758  else
5759  {
5760  /* use bestub to define y'_j */
5761 
5762  assert(!SCIPisInfinity(scip, bestub[i]));
5763  assert(!SCIPisInfinity(scip, - bestslb[i]));
5764  assert(bestslbtype[i] == -1 || bestslbtype[i] == -2);
5765  assert(bestubtype[i] > -3 && bestubtype[i] < SCIPvarGetNVubs(var));
5766 
5767  /* store for y_j that y'_j is the associated real variable
5768  * in the relaxed set
5769  */
5770  snf->origcontvars[snf->ntransvars] = probidx;
5771 
5772  if( bestubtype[i] < 0 )
5773  {
5774  SCIP_Real QUAD(val);
5775  SCIP_Real QUAD(contsolval);
5776  SCIP_Real QUAD(rowcoeftimesbestslb);
5777 
5778  /* use simple upper bound in bestslb = l_j <= y_j <= u_j = bestub to define
5779  * 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
5780  * 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,
5781  * put j into the set
5782  * N1 if a_j > 0
5783  * N2 if a_j < 0
5784  * and update the right hand side of the constraint in the relaxation
5785  * rhs = rhs - a_j l_j
5786  */
5787  SCIPquadprecSumDD(val, bestub[i], - bestslb[i]);
5788  SCIPquadprecProdQQ(val, val, rowcoef);
5789  SCIPquadprecSumDD(contsolval, solval, - bestslb[i]);
5790  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5791 
5792  if( bestubtype[i] == -2 || bestslbtype[i] == -2 )
5793  *localbdsused = TRUE;
5794 
5795  SCIPquadprecProdQD(rowcoeftimesbestslb, rowcoef, bestslb[i]);
5796 
5797  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5798  snf->origbinvars[snf->ntransvars] = -1;
5799  snf->aggrcoefsbin[snf->ntransvars] = 0.0;
5800 
5801  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5802  {
5803  snf->transvarcoefs[snf->ntransvars] = 1;
5804  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5805  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5806  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5807 
5808  /* aggregation information for y'_j */
5809  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5810  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesbestslb);
5811  }
5812  else
5813  {
5814  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5815  snf->transvarcoefs[snf->ntransvars] = - 1;
5816  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5817  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5818  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5819 
5820  /* aggregation information for y'_j */
5821  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5822  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesbestslb);
5823  }
5824  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesbestslb);
5825 
5826  SCIPdebugMsg(scip, " --> bestub used for trans: ... %s y'_%d + ..., Y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
5827  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5828  snf->ntransvars, QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesbestslb), QUAD_TO_DBL(rowcoef), bestslb[i], QUAD_TO_DBL(transrhs));
5829  }
5830  else
5831  {
5832  SCIP_Real QUAD(rowcoefbinary);
5833  SCIP_Real varsolvalbinary;
5834  SCIP_Real QUAD(val);
5835  SCIP_Real QUAD(contsolval);
5836  SCIP_Real QUAD(rowcoeftimesvubconst);
5837  int vubvarprobidx;
5838 
5839  SCIP_VAR** vubvars = SCIPvarGetVubVars(var);
5840  SCIP_Real* vubconsts = SCIPvarGetVubConstants(var);
5841  SCIP_Real* vubcoefs = SCIPvarGetVubCoefs(var);
5842 
5843  /* use variable upper bound in bestslb = l_j <= y_j <= u~_j x_j + d_j = bestub to define
5844  * 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
5845  * 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,
5846  * where c_j is the coefficient of x_j in the row, put j into the set
5847  * N1 if a_j > 0
5848  * N2 if a_j < 0
5849  * and update the right hand side of the constraint in the relaxation
5850  * rhs = rhs - a_j d_j
5851  */
5852 
5853  vubvarprobidx = SCIPvarGetProbindex(vubvars[bestubtype[i]]);
5854  assert(binvarused[vubvarprobidx] == 1);
5855