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