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