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
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
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
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
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
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
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
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
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
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
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
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  */
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 */
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 */
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 */
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 */
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 */
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  */
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 */
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 */
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 */
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  */
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
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  {
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  */
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  if( weights[k] != 0.0 )
2138  {
2139  SCIP_CALL( addOneRow(scip, aggrrow, rows[k], weights[k], sidetypebasis, allowlocal, negslack, maxaggrlen, &rowtoolong) );
2140 
2141  if( rowtoolong )
2142  return SCIP_OKAY;
2143  }
2144  }
2145  }
2146 
2147  SCIPaggrRowRemoveZeros(scip, aggrrow, valid);
2148 
2149  return SCIP_OKAY;
2150 }
2151 
2152 /** checks for cut redundancy and performs activity based coefficient tightening;
2153  * removes coefficients that are zero with QUAD_EPSILON tolerance and uses variable bounds
2154  * to remove small coefficients (relative to the maximum absolute coefficient)
2155  */
2156 static
2158  SCIP* scip, /**< SCIP data structure */
2159  SCIP_Bool cutislocal, /**< is the cut a local cut */
2160  int* cutinds, /**< variable problem indices of non-zeros in cut */
2161  SCIP_Real* cutcoefs, /**< non-zeros coefficients of cut */
2162  int* nnz, /**< number non-zeros coefficients of cut */
2163  SCIP_Real* cutrhs, /**< right hand side of cut */
2164  SCIP_Bool* success /**< pointer to return whether post-processing was succesful or cut is redundant */
2165  )
2166 {
2167  int i;
2168  SCIP_Bool redundant;
2169  SCIP_Real maxcoef;
2170  SCIP_Real minallowedcoef;
2171  SCIP_Real QUAD(rhs);
2172 
2173  assert(scip != NULL);
2174  assert(cutinds != NULL);
2175  assert(cutcoefs != NULL);
2176  assert(cutrhs != NULL);
2177  assert(success != NULL);
2178 
2179  *success = FALSE;
2180 
2181  QUAD_ASSIGN(rhs, *cutrhs);
2182 
2183  if( removeZeros(scip, SCIPfeastol(scip), cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz) )
2184  {
2185  /* right hand side was changed to infinity -> cut is redundant */
2186  return SCIP_OKAY;
2187  }
2188 
2189  if( *nnz == 0 )
2190  return SCIP_OKAY;
2191 
2192  SCIP_CALL( cutTightenCoefs(scip, cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz, &redundant) );
2193 
2194  if( redundant )
2195  {
2196  /* cut is redundant */
2197  return SCIP_OKAY;
2198  }
2199 
2200  maxcoef = 0.0;
2201  for( i = 0; i < *nnz; ++i )
2202  {
2203  SCIP_Real absval = REALABS(cutcoefs[cutinds[i]]);
2204  maxcoef = MAX(absval, maxcoef);
2205  }
2206 
2207  maxcoef /= scip->set->sepa_maxcoefratio;
2208  minallowedcoef = SCIPsumepsilon(scip);
2209  minallowedcoef = MAX(minallowedcoef, maxcoef);
2210 
2211  *success = ! removeZeros(scip, minallowedcoef, cutislocal, cutcoefs, QUAD(&rhs), cutinds, nnz);
2212  *cutrhs = QUAD_TO_DBL(rhs);
2213 
2214  return SCIP_OKAY;
2215 }
2216 
2217 
2218 /** checks for cut redundancy and performs activity based coefficient tightening;
2219  * removes coefficients that are zero with QUAD_EPSILON tolerance and uses variable bounds
2220  * to remove small coefficients (relative to the maximum absolute coefficient).
2221  * The cutcoefs must be a quad precision array, i.e. allocated with size
2222  * QUAD_ARRAY_SIZE(nvars) and accessed with QUAD_ARRAY_LOAD and QUAD_ARRAY_STORE
2223  * macros.
2224  */
2225 static
2227  SCIP* scip, /**< SCIP data structure */
2228  SCIP_Bool cutislocal, /**< is the cut a local cut */
2229  int* cutinds, /**< variable problem indices of non-zeros in cut */
2230  SCIP_Real* cutcoefs, /**< non-zeros coefficients of cut */
2231  int* nnz, /**< number non-zeros coefficients of cut */
2232  QUAD(SCIP_Real* cutrhs), /**< right hand side of cut */
2233  SCIP_Bool* success /**< pointer to return whether the cleanup was successful or if it is useless */
2234  )
2235 {
2236  int i;
2237  SCIP_Bool redundant;
2238  SCIP_Real maxcoef;
2239  SCIP_Real minallowedcoef;
2240 
2241  assert(scip != NULL);
2242  assert(cutinds != NULL);
2243  assert(cutcoefs != NULL);
2244  assert(QUAD_HI(cutrhs) != NULL);
2245  assert(success != NULL);
2246 
2247  *success = FALSE;
2248 
2249  if( removeZerosQuad(scip, SCIPfeastol(scip), cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz) )
2250  {
2251  /* right hand side was changed to infinity -> cut is redundant */
2252  return SCIP_OKAY;
2253  }
2254 
2255  if( *nnz == 0 )
2256  return SCIP_OKAY;
2257 
2258  SCIP_CALL( cutTightenCoefsQuad(scip, cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz, &redundant) );
2259  if( redundant )
2260  {
2261  /* cut is redundant */
2262  return SCIP_OKAY;
2263  }
2264 
2265  maxcoef = 0.0;
2266  for( i = 0; i < *nnz; ++i )
2267  {
2268  SCIP_Real abscoef;
2269  SCIP_Real QUAD(coef);
2270  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[i]); /* coef = cutcoefs[cutinds[i]] */
2271  abscoef = REALABS(QUAD_TO_DBL(coef));
2272  maxcoef = MAX(abscoef, maxcoef);
2273  }
2274 
2275  maxcoef /= scip->set->sepa_maxcoefratio;
2276  minallowedcoef = SCIPsumepsilon(scip);
2277  minallowedcoef = MAX(minallowedcoef, maxcoef);
2278 
2279  *success = ! removeZerosQuad(scip, minallowedcoef, cutislocal, cutcoefs, QUAD(cutrhs), cutinds, nnz);
2280 
2281  return SCIP_OKAY;
2282 }
2283 
2284 /** removes almost zero entries from the aggregation row. */
2286  SCIP* scip, /**< SCIP datastructure */
2287  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2288  SCIP_Bool* valid /**< pointer to return whether the aggregation row is still valid */
2289  )
2290 {
2291  assert(aggrrow != NULL);
2292  assert(valid != NULL);
2293 
2294  *valid = ! removeZerosQuad(scip, SCIPsumepsilon(scip), aggrrow->local, aggrrow->vals, QUAD(&aggrrow->rhs), aggrrow->inds, &aggrrow->nnz);
2295 }
2296 
2297 /** get number of aggregated rows */
2299  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2300  )
2301 {
2302  assert(aggrrow != NULL);
2303 
2304  return aggrrow->nrows;
2305 }
2306 
2307 /** get array with lp positions of rows used in aggregation */
2309  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2310  )
2311 {
2312  assert(aggrrow != NULL);
2313  assert(aggrrow->rowsinds != NULL || aggrrow->nrows == 0);
2314 
2315  return aggrrow->rowsinds;
2316 }
2317 
2318 /** get array with weights of aggregated rows */
2320  SCIP_AGGRROW* aggrrow /**< the aggregation row */
2321  )
2322 {
2323  assert(aggrrow != NULL);
2324  assert(aggrrow->rowweights != NULL || aggrrow->nrows == 0);
2325 
2326  return aggrrow->rowweights;
2327 }
2328 
2329 /** checks whether a given row has been added to the aggregation row */
2331  SCIP_AGGRROW* aggrrow, /**< the aggregation row */
2332  SCIP_ROW* row /**< row for which it is checked whether it has been added to the aggregation */
2333  )
2334 {
2335  int i;
2336  int rowind;
2337 
2338  assert(aggrrow != NULL);
2339  assert(row != NULL);
2340 
2341  rowind = SCIProwGetLPPos(row);
2342 
2343  for( i = 0; i < aggrrow->nrows; ++i )
2344  {
2345  if( aggrrow->rowsinds[i] == rowind )
2346  return TRUE;
2347  }
2348 
2349  return FALSE;
2350 }
2351 
2352 /** gets the array of corresponding variable problem indices for each non-zero in the aggregation row */
2354  SCIP_AGGRROW* aggrrow /**< aggregation row */
2355  )
2356 {
2357  assert(aggrrow != NULL);
2358 
2359  return aggrrow->inds;
2360 }
2361 
2362 /** gets the number of non-zeros in the aggregation row */
2364  SCIP_AGGRROW* aggrrow /**< aggregation row */
2365  )
2366 {
2367  assert(aggrrow != NULL);
2368 
2369  return aggrrow->nnz;
2370 }
2371 
2372 /** gets the rank of the aggregation row */
2374  SCIP_AGGRROW* aggrrow /**< aggregation row */
2375  )
2376 {
2377  assert(aggrrow != NULL);
2378 
2379  return aggrrow->rank;
2380 }
2381 
2382 /** checks if the aggregation row is only valid locally */
2384  SCIP_AGGRROW* aggrrow /**< aggregation row */
2385  )
2386 {
2387  assert(aggrrow != NULL);
2388 
2389  return aggrrow->local;
2390 }
2391 
2392 /** gets the right hand side of the aggregation row */
2394  SCIP_AGGRROW* aggrrow /**< aggregation row */
2395  )
2396 {
2397  assert(aggrrow != NULL);
2398 
2399  return QUAD_TO_DBL(aggrrow->rhs);
2400 }
2401 
2402 /** filters the given array of cuts to enforce a maximum parallelism constraints
2403  * for the given cut; moves filtered cuts to the end of the array and returns number of selected cuts */
2404 static
2406  SCIP_ROW* cut, /**< cut to filter orthogonality with */
2407  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2408  SCIP_Real* scores, /**< array with scores of cuts to perform selection algorithm */
2409  int ncuts, /**< number of cuts in given array */
2410  SCIP_Real goodscore, /**< threshold for the score to be considered a good cut */
2411  SCIP_Real goodmaxparall, /**< maximal parallelism for good cuts */
2412  SCIP_Real maxparall /**< maximal parallelism for all cuts that are not good */
2413  )
2414 {
2415  int i;
2416 
2417  assert( cut != NULL );
2418  assert( ncuts == 0 || cuts != NULL );
2419  assert( ncuts == 0 || scores != NULL );
2420 
2421  for( i = ncuts - 1; i >= 0; --i )
2422  {
2423  SCIP_Real thisparall;
2424  SCIP_Real thismaxparall;
2425 
2426  thisparall = SCIProwGetParallelism(cut, cuts[i], 'e');
2427  thismaxparall = scores[i] >= goodscore ? goodmaxparall : maxparall;
2428 
2429  if( thisparall > thismaxparall )
2430  {
2431  --ncuts;
2432  SCIPswapPointers((void**) &cuts[i], (void**) &cuts[ncuts]);
2433  SCIPswapReals(&scores[i], &scores[ncuts]);
2434  }
2435  }
2436 
2437  return ncuts;
2438 }
2439 
2440 /** move the cut with the highest score to the first position in the array; there must be at least one cut */
2441 static
2443  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2444  SCIP_Real* scores, /**< array with scores of cuts to perform selection algorithm */
2445  int ncuts /**< number of cuts in given array */
2446  )
2447 {
2448  int i;
2449  int bestpos;
2450  SCIP_Real bestscore;
2451 
2452  assert(ncuts > 0);
2453  assert(cuts != NULL);
2454  assert(scores != NULL);
2455 
2456  bestscore = scores[0];
2457  bestpos = 0;
2458 
2459  for( i = 1; i < ncuts; ++i )
2460  {
2461  if( scores[i] > bestscore )
2462  {
2463  bestpos = i;
2464  bestscore = scores[i];
2465  }
2466  }
2467 
2468  SCIPswapPointers((void**) &cuts[bestpos], (void**) &cuts[0]);
2469  SCIPswapReals(&scores[bestpos], &scores[0]);
2470 }
2471 
2472 /** perform a cut selection algorithm for the given array of cuts; the array is partitioned
2473  * so that the selected cuts come first and the remaining ones are at the end of the array
2474  */
2476  SCIP* scip, /**< SCIP data structure */
2477  SCIP_ROW** cuts, /**< array with cuts to perform selection algorithm */
2478  SCIP_RANDNUMGEN* randnumgen, /**< random number generator for tie-breaking, or NULL */
2479  SCIP_Real goodscorefac, /**< factor of best score among the given cuts to consider a cut good
2480  * and filter with less strict settings of the maximum parallelism */
2481  SCIP_Real badscorefac, /**< factor of best score among the given cuts to consider a cut bad
2482  * and discard it regardless of its parallelism to other cuts */
2483  SCIP_Real goodmaxparall, /**< maximum parallelism for good cuts */
2484  SCIP_Real maxparall, /**< maximum parallelism for non-good cuts */
2485  SCIP_Real dircutoffdistweight,/**< weight of directed cutoff distance in score calculation */
2486  SCIP_Real efficacyweight, /**< weight of efficacy (shortest cutoff distance) in score calculation */
2487  SCIP_Real objparalweight, /**< weight of objective parallelism in score calculation */
2488  SCIP_Real intsupportweight, /**< weight of integral support in score calculation */
2489  int ncuts, /**< number of cuts in given array */
2490  int nforcedcuts, /**< number of forced cuts at start of given array */
2491  int maxselectedcuts, /**< maximal number of cuts to select */
2492  int* nselectedcuts /**< pointer to return number of selected cuts */
2493  )
2494 {
2495  int i;
2496  SCIP_Real* scores;
2497  SCIP_Real goodscore;
2498  SCIP_Real badscore;
2499  SCIP_Real efficacyfac;
2500  SCIP_SOL* sol;
2501 
2502  assert( nselectedcuts != NULL );
2503 
2504  /* if all cuts are forced cuts, no selection is required */
2505  if( nforcedcuts >= MIN(ncuts, maxselectedcuts) )
2506  {
2507  *nselectedcuts = nforcedcuts;
2508  return SCIP_OKAY;
2509  }
2510  *nselectedcuts = 0;
2511 
2512  SCIP_CALL( SCIPallocBufferArray(scip, &scores, ncuts) );
2513 
2514  sol = SCIPgetBestSol(scip);
2515 
2516  efficacyfac = efficacyweight;
2517 
2518  goodscore = 0.0;
2519 
2520  /* if there is an incumbent and the factor is not 0.0, compute directed cutoff distances for the incumbent */
2521  if( sol != NULL && dircutoffdistweight > 0.0 )
2522  {
2523  for( i = 0; i < ncuts; ++i )
2524  {
2525  SCIP_Real objparallelism;
2526  SCIP_Real intsupport;
2527  SCIP_Real efficacy;
2528 
2529  intsupport = intsupportweight != 0.0 ?
2530  intsupportweight * SCIProwGetNumIntCols(cuts[i], scip->set) / (SCIP_Real) SCIProwGetNNonz(cuts[i]) :
2531  0.0;
2532 
2533  objparallelism = objparalweight != 0.0 ? objparalweight * SCIProwGetObjParallelism(cuts[i], scip->set, scip->lp) : 0.0;
2534 
2535  efficacy = SCIProwGetLPEfficacy(cuts[i], scip->set, scip->stat, scip->lp);
2536 
2537  if( SCIProwIsLocal(cuts[i]) )
2538  {
2539  scores[i] = dircutoffdistweight * efficacy;
2540  }
2541  else
2542  {
2543  scores[i] = SCIProwGetLPSolCutoffDistance(cuts[i], scip->set, scip->stat, sol, scip->lp);
2544  scores[i] = dircutoffdistweight * MAX(scores[i], efficacy);
2545  }
2546 
2547  efficacy *= efficacyfac;
2548  scores[i] += objparallelism + intsupport + efficacy;
2549 
2550  /* add small term to prefer global pool cuts */
2551  scores[i] += SCIProwIsInGlobalCutpool(cuts[i]) ? 1e-4 : 0.0;
2552 
2553  if( randnumgen != NULL )
2554  {
2555  scores[i] += SCIPrandomGetReal(randnumgen, 0.0, 1e-6);
2556  }
2557 
2558  goodscore = MAX(goodscore, scores[i]);
2559  }
2560  }
2561  else
2562  {
2563  efficacyfac += objparalweight;
2564  for( i = 0; i < ncuts; ++i )
2565  {
2566  SCIP_Real objparallelism;
2567  SCIP_Real intsupport;
2568  SCIP_Real efficacy;
2569 
2570  intsupport = intsupportweight > 0.0 ?
2571  intsupportweight * SCIProwGetNumIntCols(cuts[i], scip->set) / (SCIP_Real) SCIProwGetNNonz(cuts[i]) :
2572  0.0;
2573 
2574  objparallelism = objparalweight > 0.0 ? objparalweight * SCIProwGetObjParallelism(cuts[i], scip->set, scip->lp) : 0.0;
2575 
2576  efficacy = efficacyfac > 0.0 ?
2577  efficacyfac * SCIProwGetLPEfficacy(cuts[i], scip->set, scip->stat, scip->lp) :
2578  0.0;
2579 
2580  scores[i] = objparallelism + intsupport + efficacy;
2581 
2582  /* add small term to prefer global pool cuts */
2583  scores[i] += SCIProwIsInGlobalCutpool(cuts[i]) ? 1e-4 : 0.0;
2584 
2585  if( randnumgen != NULL )
2586  {
2587  scores[i] += SCIPrandomGetReal(randnumgen, 0.0, 1e-6);
2588  }
2589 
2590  goodscore = MAX(goodscore, scores[i]);
2591  }
2592  }
2593 
2594  /* compute values for filtering cuts */
2595  badscore = goodscore * badscorefac;
2596  goodscore *= goodscorefac;
2597 
2598  /* perform cut selection algorithm for the cuts */
2599  {
2600  int nnonforcedcuts;
2601  SCIP_ROW** nonforcedcuts;
2602  SCIP_Real* nonforcedscores;
2603 
2604  /* adjust pointers to the beginning of the non-forced cuts */
2605  nnonforcedcuts = ncuts - nforcedcuts;
2606  nonforcedcuts = cuts + nforcedcuts;
2607  nonforcedscores = scores + nforcedcuts;
2608 
2609  /* select the forced cuts first */
2610  *nselectedcuts = nforcedcuts;
2611  for( i = 0; i < nforcedcuts && nnonforcedcuts > 0; ++i )
2612  {
2613  nnonforcedcuts = filterWithParallelism(cuts[i], nonforcedcuts, nonforcedscores, nnonforcedcuts, goodscore, goodmaxparall, maxparall);
2614  }
2615 
2616  /* if the maximal number of cuts was exceeded after selecting the forced cuts, we can stop here */
2617  if( *nselectedcuts >= maxselectedcuts )
2618  goto TERMINATE;
2619 
2620  /* now greedily select the remaining cuts */
2621  while( nnonforcedcuts > 0 )
2622  {
2623  SCIP_ROW* selectedcut;
2624 
2625  selectBestCut(nonforcedcuts, nonforcedscores, nnonforcedcuts);
2626  selectedcut = nonforcedcuts[0];
2627 
2628  /* if the best cut of the remaining cuts is considered bad, we discard it and all remaining cuts */
2629  if( nonforcedscores[0] < badscore )
2630  goto TERMINATE;
2631 
2632  ++(*nselectedcuts);
2633 
2634  /* if the maximal number of cuts was selected, we can stop here */
2635  if( *nselectedcuts == maxselectedcuts )
2636  goto TERMINATE;
2637 
2638  /* move the pointers to the next position and filter the remaining cuts to enforce the maximum parallelism constraint */
2639  ++nonforcedcuts;
2640  ++nonforcedscores;
2641  --nnonforcedcuts;
2642 
2643  nnonforcedcuts = filterWithParallelism(selectedcut, nonforcedcuts, nonforcedscores, nnonforcedcuts, goodscore, goodmaxparall, maxparall);
2644  }
2645  }
2646 
2647  TERMINATE:
2648  SCIPfreeBufferArray(scip, &scores);
2649 
2650  return SCIP_OKAY;
2651 }
2652 
2653 /* =========================================== c-MIR =========================================== */
2654 
2655 #define MAXCMIRSCALE 1e+6 /**< maximal scaling (scale/(1-f0)) allowed in c-MIR calculations */
2656 
2657 /** finds the best lower bound of the variable to use for MIR transformation */
2658 static
2660  SCIP* scip, /**< SCIP data structure */
2661  SCIP_VAR* var, /**< problem variable */
2662  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2663  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2664  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2665  SCIP_Real* bestlb, /**< pointer to store best bound value */
2666  int* bestlbtype /**< pointer to store best bound type */
2667  )
2668 {
2669  assert(bestlb != NULL);
2670  assert(bestlbtype != NULL);
2671 
2672  *bestlb = SCIPvarGetLbGlobal(var);
2673  *bestlbtype = -1;
2674 
2675  if( allowlocal )
2676  {
2677  SCIP_Real loclb;
2678 
2679  loclb = SCIPvarGetLbLocal(var);
2680  if( SCIPisGT(scip, loclb, *bestlb) )
2681  {
2682  *bestlb = loclb;
2683  *bestlbtype = -2;
2684  }
2685  }
2686 
2687  if( usevbds && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
2688  {
2689  SCIP_Real bestvlb;
2690  int bestvlbidx;
2691 
2692  SCIP_CALL( SCIPgetVarClosestVlb(scip, var, sol, &bestvlb, &bestvlbidx) );
2693  if( bestvlbidx >= 0
2694  && (bestvlb > *bestlb || (*bestlbtype < 0 && SCIPisGE(scip, bestvlb, *bestlb))) )
2695  {
2696  SCIP_VAR** vlbvars;
2697 
2698  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2699  /**@todo this check is not needed for continuous variables; but allowing all but binary variables
2700  * to be replaced by variable bounds seems to be buggy (wrong result on gesa2)
2701  */
2702  vlbvars = SCIPvarGetVlbVars(var);
2703  assert(vlbvars != NULL);
2704  if( SCIPvarGetProbindex(vlbvars[bestvlbidx]) < SCIPvarGetProbindex(var) )
2705  {
2706  *bestlb = bestvlb;
2707  *bestlbtype = bestvlbidx;
2708  }
2709  }
2710  }
2711 
2712  return SCIP_OKAY;
2713 }
2714 
2715 /** finds the best upper bound of the variable to use for MIR transformation */
2716 static
2718  SCIP* scip, /**< SCIP data structure */
2719  SCIP_VAR* var, /**< problem variable */
2720  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2721  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2722  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2723  SCIP_Real* bestub, /**< pointer to store best bound value */
2724  int* bestubtype /**< pointer to store best bound type */
2725  )
2726 {
2727  assert(bestub != NULL);
2728  assert(bestubtype != NULL);
2729 
2730  *bestub = SCIPvarGetUbGlobal(var);
2731  *bestubtype = -1;
2732 
2733  if( allowlocal )
2734  {
2735  SCIP_Real locub;
2736 
2737  locub = SCIPvarGetUbLocal(var);
2738  if( SCIPisLT(scip, locub, *bestub) )
2739  {
2740  *bestub = locub;
2741  *bestubtype = -2;
2742  }
2743  }
2744 
2745  if( usevbds && SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS )
2746  {
2747  SCIP_Real bestvub;
2748  int bestvubidx;
2749 
2750  SCIP_CALL( SCIPgetVarClosestVub(scip, var, sol, &bestvub, &bestvubidx) );
2751  if( bestvubidx >= 0
2752  && (bestvub < *bestub || (*bestubtype < 0 && SCIPisLE(scip, bestvub, *bestub))) )
2753  {
2754  SCIP_VAR** vubvars;
2755 
2756  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2757  /**@todo this check is not needed for continuous variables; but allowing all but binary variables
2758  * to be replaced by variable bounds seems to be buggy (wrong result on gesa2)
2759  */
2760  vubvars = SCIPvarGetVubVars(var);
2761  assert(vubvars != NULL);
2762  if( SCIPvarGetProbindex(vubvars[bestvubidx]) < SCIPvarGetProbindex(var) )
2763  {
2764  *bestub = bestvub;
2765  *bestubtype = bestvubidx;
2766  }
2767  }
2768  }
2769 
2770  return SCIP_OKAY;
2771 }
2772 
2773 /** determine the best bounds with respect to the given solution for complementing the given variable */
2774 static
2776  SCIP* scip, /**< SCIP data structure */
2777  SCIP_VAR* var, /**< variable to determine best bound for */
2778  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2779  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
2780  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2781  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2782  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
2783  SCIP_Bool ignoresol, /**< should the LP solution be ignored? (eg, apply MIR to dualray) */
2784  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
2785  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
2786  * NULL for using closest bound for all variables */
2787  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
2788  * NULL for using closest bound for all variables */
2789  SCIP_Real* bestlb, /**< pointer to store best lower bound of variable */
2790  SCIP_Real* bestub, /**< pointer to store best upper bound of variable */
2791  int* bestlbtype, /**< pointer to store type of best lower bound of variable */
2792  int* bestubtype, /**< pointer to store type of best upper bound of variable */
2793  SCIP_BOUNDTYPE* selectedbound, /**< pointer to store whether the lower bound or the upper bound should be preferred */
2794  SCIP_Bool* freevariable /**< pointer to store if this is a free variable */
2795  )
2796 {
2797  int v;
2798 
2799  v = SCIPvarGetProbindex(var);
2800 
2801  /* check if the user specified a bound to be used */
2802  if( boundsfortrans != NULL && boundsfortrans[v] > -3 )
2803  {
2804  assert(SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS || ( boundsfortrans[v] == -2 || boundsfortrans[v] == -1 ));
2805 
2806  /* user has explicitly specified a bound to be used */
2807  if( boundtypesfortrans[v] == SCIP_BOUNDTYPE_LOWER )
2808  {
2809  /* user wants to use lower bound */
2810  *bestlbtype = boundsfortrans[v];
2811  if( *bestlbtype == -1 )
2812  *bestlb = SCIPvarGetLbGlobal(var); /* use global standard lower bound */
2813  else if( *bestlbtype == -2 )
2814  *bestlb = SCIPvarGetLbLocal(var); /* use local standard lower bound */
2815  else
2816  {
2817  SCIP_VAR** vlbvars;
2818  SCIP_Real* vlbcoefs;
2819  SCIP_Real* vlbconsts;
2820  int k;
2821 
2822  assert(!ignoresol);
2823 
2824  /* use the given variable lower bound */
2825  vlbvars = SCIPvarGetVlbVars(var);
2826  vlbcoefs = SCIPvarGetVlbCoefs(var);
2827  vlbconsts = SCIPvarGetVlbConstants(var);
2828  k = boundsfortrans[v];
2829  assert(k >= 0 && k < SCIPvarGetNVlbs(var));
2830  assert(vlbvars != NULL);
2831  assert(vlbcoefs != NULL);
2832  assert(vlbconsts != NULL);
2833 
2834  *bestlb = vlbcoefs[k] * (sol == NULL ? SCIPvarGetLPSol(vlbvars[k]) : SCIPgetSolVal(scip, sol, vlbvars[k])) + vlbconsts[k];
2835  }
2836 
2837  assert(!SCIPisInfinity(scip, - *bestlb));
2838  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2839 
2840  /* find closest upper bound in standard upper bound (and variable upper bounds for continuous variables) */
2841  SCIP_CALL( findBestUb(scip, var, sol, usevbds && fixintegralrhs, allowlocal && fixintegralrhs, bestub, bestubtype) );
2842  }
2843  else
2844  {
2845  assert(boundtypesfortrans[v] == SCIP_BOUNDTYPE_UPPER);
2846 
2847  /* user wants to use upper bound */
2848  *bestubtype = boundsfortrans[v];
2849  if( *bestubtype == -1 )
2850  *bestub = SCIPvarGetUbGlobal(var); /* use global standard upper bound */
2851  else if( *bestubtype == -2 )
2852  *bestub = SCIPvarGetUbLocal(var); /* use local standard upper bound */
2853  else
2854  {
2855  SCIP_VAR** vubvars;
2856  SCIP_Real* vubcoefs;
2857  SCIP_Real* vubconsts;
2858  int k;
2859 
2860  assert(!ignoresol);
2861 
2862  /* use the given variable upper bound */
2863  vubvars = SCIPvarGetVubVars(var);
2864  vubcoefs = SCIPvarGetVubCoefs(var);
2865  vubconsts = SCIPvarGetVubConstants(var);
2866  k = boundsfortrans[v];
2867  assert(k >= 0 && k < SCIPvarGetNVubs(var));
2868  assert(vubvars != NULL);
2869  assert(vubcoefs != NULL);
2870  assert(vubconsts != NULL);
2871 
2872  /* we have to avoid cyclic variable bound usage, so we enforce to use only variable bounds variables of smaller index */
2873  *bestub = vubcoefs[k] * (sol == NULL ? SCIPvarGetLPSol(vubvars[k]) : SCIPgetSolVal(scip, sol, vubvars[k])) + vubconsts[k];
2874  }
2875 
2876  assert(!SCIPisInfinity(scip, *bestub));
2877  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2878 
2879  /* find closest lower bound in standard lower bound (and variable lower bounds for continuous variables) */
2880  SCIP_CALL( findBestLb(scip, var, sol, usevbds && fixintegralrhs, allowlocal && fixintegralrhs, bestlb, bestlbtype) );
2881  }
2882  }
2883  else
2884  {
2885  SCIP_Real varsol;
2886 
2887  /* bound selection should be done automatically */
2888 
2889  /* find closest lower bound in standard lower bound (and variable lower bounds for continuous variables) */
2890  SCIP_CALL( findBestLb(scip, var, sol, usevbds, allowlocal, bestlb, bestlbtype) );
2891 
2892  /* find closest upper bound in standard upper bound (and variable upper bounds for continuous variables) */
2893  SCIP_CALL( findBestUb(scip, var, sol, usevbds, allowlocal, bestub, bestubtype) );
2894 
2895  /* check, if variable is free variable */
2896  if( SCIPisInfinity(scip, - *bestlb) && SCIPisInfinity(scip, *bestub) )
2897  {
2898  /* we found a free variable in the row with non-zero coefficient
2899  * -> MIR row can't be transformed in standard form
2900  */
2901  *freevariable = TRUE;
2902  return SCIP_OKAY;
2903  }
2904 
2905  if( !ignoresol )
2906  {
2907  /* select transformation bound */
2908  varsol = (sol == NULL ? SCIPvarGetLPSol(var) : SCIPgetSolVal(scip, sol, var));
2909 
2910  if( SCIPisInfinity(scip, *bestub) ) /* if there is no ub, use lb */
2911  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2912  else if( SCIPisInfinity(scip, - *bestlb) ) /* if there is no lb, use ub */
2913  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2914  else if( SCIPisLT(scip, varsol, (1.0 - boundswitch) * (*bestlb) + boundswitch * (*bestub)) )
2915  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2916  else if( SCIPisGT(scip, varsol, (1.0 - boundswitch) * (*bestlb) + boundswitch * (*bestub)) )
2917  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2918  else if( *bestlbtype == -1 ) /* prefer global standard bounds */
2919  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2920  else if( *bestubtype == -1 ) /* prefer global standard bounds */
2921  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2922  else if( *bestlbtype >= 0 ) /* prefer variable bounds over local bounds */
2923  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2924  else if( *bestubtype >= 0 ) /* prefer variable bounds over local bounds */
2925  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2926  else /* no decision yet? just use lower bound */
2927  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2928  }
2929  else
2930  {
2931  SCIP_Real glbub = SCIPvarGetUbGlobal(var);
2932  SCIP_Real glblb = SCIPvarGetLbGlobal(var);
2933  SCIP_Real distlb = REALABS(glblb - *bestlb);
2934  SCIP_Real distub = REALABS(glbub - *bestub);
2935 
2936  assert(!SCIPisInfinity(scip, - *bestlb) || !SCIPisInfinity(scip, *bestub));
2937 
2938  if( SCIPisInfinity(scip, - *bestlb) )
2939  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2940  else if( !SCIPisNegative(scip, *bestlb) )
2941  {
2942  if( SCIPisInfinity(scip, *bestub) )
2943  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2944  else if( SCIPisZero(scip, glblb) )
2945  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2946  else if( SCIPisLE(scip, distlb, distub) )
2947  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2948  else
2949  *selectedbound = SCIP_BOUNDTYPE_UPPER;
2950  }
2951  else
2952  {
2953  assert(!SCIPisInfinity(scip, - *bestlb));
2954  *selectedbound = SCIP_BOUNDTYPE_LOWER;
2955  }
2956  }
2957  }
2958 
2959  return SCIP_OKAY;
2960 }
2961 
2962 /** Transform equation \f$ a \cdot x = b; lb \leq x \leq ub \f$ into standard form
2963  * \f$ a^\prime \cdot x^\prime = b,\; 0 \leq x^\prime \leq ub' \f$.
2964  *
2965  * Transform variables (lb or ub):
2966  * \f[
2967  * \begin{array}{llll}
2968  * 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}\\
2969  * 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}
2970  * \end{array}
2971  * \f]
2972  * and move the constant terms \f$ a_j\, lb_j \f$ or \f$ a_j\, ub_j \f$ to the rhs.
2973  *
2974  * Transform variables (vlb or vub):
2975  * \f[
2976  * \begin{array}{llll}
2977  * 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.} \\
2978  * 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.}
2979  * \end{array}
2980  * \f]
2981  * 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:
2982  * \f[
2983  * \begin{array}{ll}
2984  * a_{zl_j} := a_{zl_j} + a_j\, bl_j,& \mbox{or} \\
2985  * a_{zu_j} := a_{zu_j} + a_j\, bu_j &
2986  * \end{array}
2987  * \f]
2988  */
2989 static
2991  SCIP* scip, /**< SCIP data structure */
2992  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
2993  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
2994  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
2995  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
2996  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
2997  SCIP_Bool ignoresol, /**< should the LP solution be ignored? (eg, apply MIR to dualray) */
2998  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
2999  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
3000  * NULL for using closest bound for all variables */
3001  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
3002  * NULL for using closest bound for all variables */
3003  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
3004  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
3005  SCIP_Real* cutcoefs, /**< array of coefficients of cut */
3006  QUAD(SCIP_Real* cutrhs), /**< pointer to right hand side of cut */
3007  int* cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3008  int* nnz, /**< number of non-zeros in cut */
3009  int* varsign, /**< stores the sign of the transformed variable in summation */
3010  int* boundtype, /**< stores the bound used for transformed variable:
3011  * vlb/vub_idx, or -1 for global lb/ub, or -2 for local lb/ub */
3012  SCIP_Bool* freevariable, /**< stores whether a free variable was found in MIR row -> invalid summation */
3013  SCIP_Bool* localbdsused /**< pointer to store whether local bounds were used in transformation */
3014  )
3015 {
3016  SCIP_Real QUAD(tmp);
3017  SCIP_Real* bestlbs;
3018  SCIP_Real* bestubs;
3019  int* bestlbtypes;
3020  int* bestubtypes;
3021  SCIP_BOUNDTYPE* selectedbounds;
3022  int i;
3023  int aggrrowintstart;
3024  int nvars;
3025  int firstcontvar;
3026  SCIP_VAR** vars;
3027 
3028  assert(varsign != NULL);
3029  assert(boundtype != NULL);
3030  assert(freevariable != NULL);
3031  assert(localbdsused != NULL);
3032 
3033  *freevariable = FALSE;
3034  *localbdsused = FALSE;
3035 
3036  /* allocate temporary memory to store best bounds and bound types */
3037  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbs, 2*(*nnz)) );
3038  SCIP_CALL( SCIPallocBufferArray(scip, &bestubs, 2*(*nnz)) );
3039  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbtypes, 2*(*nnz)) );
3040  SCIP_CALL( SCIPallocBufferArray(scip, &bestubtypes, 2*(*nnz)) );
3041  SCIP_CALL( SCIPallocBufferArray(scip, &selectedbounds, 2*(*nnz)) );
3042 
3043  /* start with continuous variables, because using variable bounds can affect the untransformed integral
3044  * variables, and these changes have to be incorporated in the transformation of the integral variables
3045  * (continuous variables have largest problem indices!)
3046  */
3047  SCIPsortDownInt(cutinds, *nnz);
3048 
3049  vars = SCIPgetVars(scip);
3050  nvars = SCIPgetNVars(scip);
3051  firstcontvar = nvars - SCIPgetNContVars(scip);
3052 
3053  /* determine the best bounds for the continous variables */
3054  for( i = 0; i < *nnz && cutinds[i] >= firstcontvar; ++i )
3055  {
3056  SCIP_CALL( determineBestBounds(scip, vars[cutinds[i]], sol, boundswitch, usevbds, allowlocal, fixintegralrhs,
3057  ignoresol, boundsfortrans, boundtypesfortrans,
3058  bestlbs + i, bestubs + i, bestlbtypes + i, bestubtypes + i, selectedbounds + i, freevariable) );
3059 
3060  if( *freevariable )
3061  goto TERMINATE;
3062  }
3063 
3064  /* remember start of integer variables in the aggrrow */
3065  aggrrowintstart = i;
3066 
3067  /* perform bound substitution for continuous variables */
3068  for( i = 0; i < aggrrowintstart; ++i )
3069  {
3070  SCIP_Real QUAD(coef);
3071  int v = cutinds[i];
3072  SCIP_VAR* var = vars[v];
3073 
3074  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3075 
3076  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
3077  {
3078  assert(!SCIPisInfinity(scip, -bestlbs[i]));
3079 
3080  /* use lower bound as transformation bound: x'_j := x_j - lb_j */
3081  boundtype[i] = bestlbtypes[i];
3082  varsign[i] = +1;
3083 
3084  /* standard (bestlbtype < 0) or variable (bestlbtype >= 0) lower bound? */
3085  if( bestlbtypes[i] < 0 )
3086  {
3087  SCIPquadprecProdQD(tmp, coef, bestlbs[i]);
3088  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3089  *localbdsused = *localbdsused || (bestlbtypes[i] == -2);
3090  }
3091  else
3092  {
3093  SCIP_VAR** vlbvars;
3094  SCIP_Real* vlbcoefs;
3095  SCIP_Real* vlbconsts;
3096  SCIP_Real QUAD(zcoef);
3097  int zidx;
3098 
3099  vlbvars = SCIPvarGetVlbVars(var);
3100  vlbcoefs = SCIPvarGetVlbCoefs(var);
3101  vlbconsts = SCIPvarGetVlbConstants(var);
3102  assert(vlbvars != NULL);
3103  assert(vlbcoefs != NULL);
3104  assert(vlbconsts != NULL);
3105 
3106  assert(0 <= bestlbtypes[i] && bestlbtypes[i] < SCIPvarGetNVlbs(var));
3107  assert(SCIPvarIsActive(vlbvars[bestlbtypes[i]]));
3108  zidx = SCIPvarGetProbindex(vlbvars[bestlbtypes[i]]);
3109  assert(0 <= zidx && zidx < firstcontvar);
3110 
3111  SCIPquadprecProdQD(tmp, coef, vlbconsts[bestlbtypes[i]]);
3112  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3113 
3114  /* check if integral variable already exists in the row and update sparsity pattern */
3115  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3116  if( QUAD_HI(zcoef) == 0.0 )
3117  cutinds[(*nnz)++] = zidx;
3118 
3119  SCIPquadprecProdQD(coef, coef, vlbcoefs[bestlbtypes[i]]);
3120  SCIPquadprecSumQQ(zcoef, zcoef, coef);
3121  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3122  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3123  assert(QUAD_HI(zcoef) != 0.0);
3124  }
3125  }
3126  else
3127  {
3128  assert(!SCIPisInfinity(scip, bestubs[i]));
3129 
3130  /* use upper bound as transformation bound: x'_j := ub_j - x_j */
3131  boundtype[i] = bestubtypes[i];
3132  varsign[i] = -1;
3133 
3134  /* standard (bestubtype < 0) or variable (bestubtype >= 0) upper bound? */
3135  if( bestubtypes[i] < 0 )
3136  {
3137  SCIPquadprecProdQD(tmp, coef, bestubs[i]);
3138  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3139  *localbdsused = *localbdsused || (bestubtypes[i] == -2);
3140  }
3141  else
3142  {
3143  SCIP_VAR** vubvars;
3144  SCIP_Real* vubcoefs;
3145  SCIP_Real* vubconsts;
3146  SCIP_Real QUAD(zcoef);
3147  int zidx;
3148 
3149  vubvars = SCIPvarGetVubVars(var);
3150  vubcoefs = SCIPvarGetVubCoefs(var);
3151  vubconsts = SCIPvarGetVubConstants(var);
3152  assert(vubvars != NULL);
3153  assert(vubcoefs != NULL);
3154  assert(vubconsts != NULL);
3155 
3156  assert(0 <= bestubtypes[i] && bestubtypes[i] < SCIPvarGetNVubs(var));
3157  assert(SCIPvarIsActive(vubvars[bestubtypes[i]]));
3158  zidx = SCIPvarGetProbindex(vubvars[bestubtypes[i]]);
3159  assert(zidx >= 0);
3160 
3161  SCIPquadprecProdQD(tmp, coef, vubconsts[bestubtypes[i]]);
3162  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3163 
3164  /* check if integral variable already exists in the row and update sparsity pattern */
3165  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3166  if( QUAD_HI(zcoef) == 0.0 )
3167  cutinds[(*nnz)++] = zidx;
3168 
3169  SCIPquadprecProdQD(coef, coef, vubcoefs[bestubtypes[i]]);
3170  SCIPquadprecSumQQ(zcoef, zcoef, coef);
3171  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3172  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3173  assert(QUAD_HI(zcoef) != 0.0);
3174  }
3175  }
3176  }
3177 
3178  /* remove integral variables that now have a zero coefficient due to variable bound usage of continuous variables
3179  * and determine the bound to use for the integer variables that are left
3180  */
3181  while( i < *nnz )
3182  {
3183  SCIP_Real QUAD(coef);
3184  int v = cutinds[i];
3185  assert(cutinds[i] < firstcontvar);
3186 
3187  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3188 
3189  /* due to variable bound usage for the continous variables cancellation may have occurred */
3190  if( EPSZ(QUAD_TO_DBL(coef), QUAD_EPSILON) )
3191  {
3192  QUAD_ASSIGN(coef, 0.0);
3193  QUAD_ARRAY_STORE(cutcoefs, v, coef);
3194  --(*nnz);
3195  cutinds[i] = cutinds[*nnz];
3196  /* do not increase i, since last element is copied to the i-th position */
3197  continue;
3198  }
3199 
3200  /* determine the best bounds for the integral variable, usevbd can be set to FALSE here as vbds are only used for continous variables */
3201  SCIP_CALL( determineBestBounds(scip, vars[v], sol, boundswitch, FALSE, allowlocal, fixintegralrhs,
3202  ignoresol, boundsfortrans, boundtypesfortrans,
3203  bestlbs + i, bestubs + i, bestlbtypes + i, bestubtypes + i, selectedbounds + i, freevariable) );
3204 
3205  /* increase i */
3206  ++i;
3207 
3208  if( *freevariable )
3209  goto TERMINATE;
3210  }
3211 
3212  /* now perform the bound substitution on the remaining integral variables which only uses standard bounds */
3213  for( i = aggrrowintstart; i < *nnz; ++i )
3214  {
3215  SCIP_Real QUAD(coef);
3216  int v = cutinds[i];
3217 
3218  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3219 
3220  /* perform bound substitution */
3221  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
3222  {
3223  assert(!SCIPisInfinity(scip, - bestlbs[i]));
3224  assert(bestlbtypes[i] < 0);
3225 
3226  /* use lower bound as transformation bound: x'_j := x_j - lb_j */
3227  boundtype[i] = bestlbtypes[i];
3228  varsign[i] = +1;
3229 
3230  /* standard (bestlbtype < 0) or variable (bestlbtype >= 0) lower bound? */
3231  SCIPquadprecProdQD(tmp, coef, bestlbs[i]);
3232  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3233  *localbdsused = *localbdsused || (bestlbtypes[i] == -2);
3234  }
3235  else
3236  {
3237  assert(!SCIPisInfinity(scip, bestubs[i]));
3238  assert(bestubtypes[i] < 0);
3239 
3240  /* use upper bound as transformation bound: x'_j := ub_j - x_j */
3241  boundtype[i] = bestubtypes[i];
3242  varsign[i] = -1;
3243 
3244  /* standard (bestubtype < 0) or variable (bestubtype >= 0) upper bound? */
3245  SCIPquadprecProdQD(tmp, coef, bestubs[i]);
3246  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3247  *localbdsused = *localbdsused || (bestubtypes[i] == -2);
3248  }
3249  }
3250 
3251  if( fixintegralrhs )
3252  {
3253  SCIP_Real f0;
3254 
3255  /* check if rhs is fractional */
3256  f0 = EPSFRAC(QUAD_TO_DBL(*cutrhs), SCIPsumepsilon(scip));
3257  if( f0 < minfrac || f0 > maxfrac )
3258  {
3259  SCIP_Real bestviolgain;
3260  SCIP_Real bestnewf0;
3261  int besti;
3262 
3263  /* choose complementation of one variable differently such that f0 is in correct range */
3264  besti = -1;
3265  bestviolgain = -1e+100;
3266  bestnewf0 = 1.0;
3267  for( i = 0; i < *nnz; i++ )
3268  {
3269  int v;
3270  SCIP_Real QUAD(coef);
3271 
3272  v = cutinds[i];
3273  assert(0 <= v && v < nvars);
3274 
3275  QUAD_ARRAY_LOAD(coef, cutcoefs, v);
3276  assert(!EPSZ(QUAD_TO_DBL(coef), QUAD_EPSILON));
3277 
3278  if( boundtype[i] < 0
3279  && ((varsign[i] == +1 && !SCIPisInfinity(scip, bestubs[i]) && bestubtypes[i] < 0)
3280  || (varsign[i] == -1 && !SCIPisInfinity(scip, -bestlbs[i]) && bestlbtypes[i] < 0)) )
3281  {
3282  SCIP_Real fj;
3283  SCIP_Real newfj;
3284  SCIP_Real newrhs;
3285  SCIP_Real newf0;
3286  SCIP_Real solval;
3287  SCIP_Real viol;
3288  SCIP_Real newviol;
3289  SCIP_Real violgain;
3290 
3291  /* currently: a'_j = varsign * a_j -> f'_j = a'_j - floor(a'_j)
3292  * after complementation: a''_j = -varsign * a_j -> f''_j = a''_j - floor(a''_j) = 1 - f'_j
3293  * rhs'' = rhs' + varsign * a_j * (lb_j - ub_j)
3294  * cut violation from f0 and fj: f'_0 - f'_j * x'_j
3295  * after complementation: f''_0 - f''_j * x''_j
3296  *
3297  * for continuous variables, we just set f'_j = f''_j = |a'_j|
3298  */
3299  newrhs = QUAD_TO_DBL(*cutrhs) + varsign[i] * QUAD_TO_DBL(coef) * (bestlbs[i] - bestubs[i]);
3300  newf0 = EPSFRAC(newrhs, SCIPsumepsilon(scip));
3301  if( newf0 < minfrac || newf0 > maxfrac )
3302  continue;
3303  if( v >= firstcontvar )
3304  {
3305  fj = REALABS(QUAD_TO_DBL(coef));
3306  newfj = fj;
3307  }
3308  else
3309  {
3310  fj = SCIPfrac(scip, varsign[i] * QUAD_TO_DBL(coef));
3311  newfj = SCIPfrac(scip, -varsign[i] * QUAD_TO_DBL(coef));
3312  }
3313 
3314  if( !ignoresol )
3315  {
3316  solval = (sol == NULL ? SCIPvarGetLPSol(vars[v]) : SCIPgetSolVal(scip, sol, vars[v]));
3317  viol = f0 - fj * (varsign[i] == +1 ? solval - bestlbs[i] : bestubs[i] - solval);
3318  newviol = newf0 - newfj * (varsign[i] == -1 ? solval - bestlbs[i] : bestubs[i] - solval);
3319  violgain = newviol - viol;
3320  }
3321  else
3322  {
3323  /* todo: this should be done, this can improve the dualray significantly */
3324  SCIPerrorMessage("Cannot handle closest bounds with ignoring the LP solution.\n");
3325  return SCIP_INVALIDCALL;
3326  }
3327 
3328  /* prefer larger violations; for equal violations, prefer smaller f0 values since then the possibility that
3329  * we f_j > f_0 is larger and we may improve some coefficients in rounding
3330  */
3331  if( SCIPisGT(scip, violgain, bestviolgain)
3332  || (SCIPisGE(scip, violgain, bestviolgain) && newf0 < bestnewf0) )
3333  {
3334  besti = i;
3335  bestviolgain = violgain;
3336  bestnewf0 = newf0;
3337  }
3338  }
3339  }
3340 
3341  if( besti >= 0 )
3342  {
3343  SCIP_Real QUAD(coef);
3344  assert(besti < *nnz);
3345  assert(boundtype[besti] < 0);
3346  assert(!SCIPisInfinity(scip, -bestlbs[besti]));
3347  assert(!SCIPisInfinity(scip, bestubs[besti]));
3348 
3349  QUAD_ARRAY_LOAD(coef, cutcoefs, cutinds[besti]);
3350  QUAD_SCALE(coef, varsign[besti]);
3351 
3352  /* switch the complementation of this variable */
3353  SCIPquadprecSumDD(tmp, bestlbs[besti], - bestubs[besti]);
3354  SCIPquadprecProdQQ(tmp, tmp, coef);
3355  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3356 
3357  if( varsign[besti] == +1 )
3358  {
3359  /* switch to upper bound */
3360  assert(bestubtypes[besti] < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
3361  boundtype[besti] = bestubtypes[besti];
3362  varsign[besti] = -1;
3363  }
3364  else
3365  {
3366  /* switch to lower bound */
3367  assert(bestlbtypes[besti] < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
3368  boundtype[besti] = bestlbtypes[besti];
3369  varsign[besti] = +1;
3370  }
3371  *localbdsused = *localbdsused || (boundtype[besti] == -2);
3372  }
3373  }
3374  }
3375 
3376  TERMINATE:
3377 
3378  /*free temporary memory */
3379  SCIPfreeBufferArray(scip, &selectedbounds);
3380  SCIPfreeBufferArray(scip, &bestubtypes);
3381  SCIPfreeBufferArray(scip, &bestlbtypes);
3382  SCIPfreeBufferArray(scip, &bestubs);
3383  SCIPfreeBufferArray(scip, &bestlbs);
3384 
3385  return SCIP_OKAY;
3386 }
3387 
3388 /** 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$
3389  * \f[
3390  * \begin{array}{rll}
3391  * integers :& \tilde{a}_j = down(a^\prime_j), & if \qquad f_j \leq f_0 \\
3392  * & \tilde{a}_j = down(a^\prime_j) + (f_j - f_0)/(1 - f_0),& if \qquad f_j > f_0 \\
3393  * continuous:& \tilde{a}_j = 0, & if \qquad a^\prime_j \geq 0 \\
3394  * & \tilde{a}_j = a^\prime_j/(1 - f_0), & if \qquad a^\prime_j < 0
3395  * \end{array}
3396  * \f]
3397  *
3398  * Transform inequality back to \f$ \hat{a} \cdot x \leq rhs \f$:
3399  *
3400  * (lb or ub):
3401  * \f[
3402  * \begin{array}{lllll}
3403  * 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} \\
3404  * 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}
3405  * \end{array}
3406  * \f]
3407  * and move the constant terms
3408  * \f[
3409  * \begin{array}{cl}
3410  * -\tilde{a}_j \cdot lb_j = -\hat{a}_j \cdot lb_j,& \mbox{or} \\
3411  * \tilde{a}_j \cdot ub_j = -\hat{a}_j \cdot ub_j &
3412  * \end{array}
3413  * \f]
3414  * to the rhs.
3415  *
3416  * (vlb or vub):
3417  * \f[
3418  * \begin{array}{lllll}
3419  * 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)} \\
3420  * 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)}
3421  * \end{array}
3422  * \f]
3423  * move the constant terms
3424  * \f[
3425  * \begin{array}{cl}
3426  * -\tilde{a}_j\, dl_j = -\hat{a}_j\, dl_j,& \mbox{or} \\
3427  * \tilde{a}_j\, du_j = -\hat{a}_j\, du_j &
3428  * \end{array}
3429  * \f]
3430  * to the rhs, and update the VB variable coefficients:
3431  * \f[
3432  * \begin{array}{ll}
3433  * \hat{a}_{zl_j} := \hat{a}_{zl_j} - \tilde{a}_j\, bl_j = \hat{a}_{zl_j} - \hat{a}_j\, bl_j,& \mbox{or} \\
3434  * \hat{a}_{zu_j} := \hat{a}_{zu_j} + \tilde{a}_j\, bu_j = \hat{a}_{zu_j} - \hat{a}_j\, bu_j &
3435  * \end{array}
3436  * \f]
3437  */
3438 static
3440  SCIP* scip, /**< SCIP data structure */
3441  SCIP_Real*RESTRICT cutcoefs, /**< array of coefficients of cut */
3442  QUAD(SCIP_Real*RESTRICT cutrhs), /**< pointer to right hand side of cut */
3443  int*RESTRICT cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3444  int*RESTRICT nnz, /**< number of non-zeros in cut */
3445  int*RESTRICT varsign, /**< stores the sign of the transformed variable in summation */
3446  int*RESTRICT boundtype, /**< stores the bound used for transformed variable (vlb/vub_idx or -1 for lb/ub) */
3447  QUAD(SCIP_Real f0) /**< fractional value of rhs */
3448  )
3449 {
3450  SCIP_Real QUAD(tmp);
3451  SCIP_Real QUAD(onedivoneminusf0);
3452  int i;
3453  int firstcontvar;
3454  SCIP_VAR** vars;
3455  int ndelcontvars;
3456 
3457  assert(QUAD_HI(cutrhs) != NULL);
3458  assert(cutcoefs != NULL);
3459  assert(cutinds != NULL);
3460  assert(nnz != NULL);
3461  assert(boundtype != NULL);
3462  assert(varsign != NULL);
3463  assert(0.0 < QUAD_TO_DBL(f0) && QUAD_TO_DBL(f0) < 1.0);
3464 
3465  SCIPquadprecSumQD(onedivoneminusf0, -f0, 1.0);
3466  SCIPquadprecDivDQ(onedivoneminusf0, 1.0, onedivoneminusf0);
3467 
3468  /* Loop backwards to process integral variables first and be able to delete coefficients of integral variables
3469  * without destroying the ordering of the aggrrow's non-zeros.
3470  * (due to sorting in cutsTransformMIR the ordering is continuous before integral)
3471  */
3472 
3473  firstcontvar = SCIPgetNVars(scip) - SCIPgetNContVars(scip);
3474  vars = SCIPgetVars(scip);
3475 #ifndef NDEBUG
3476  /*in debug mode check that all continuous variables of the aggrrow come before the integral variables */
3477  i = 0;
3478  while( i < *nnz && cutinds[i] >= firstcontvar )
3479  ++i;
3480 
3481  while( i < *nnz )
3482  {
3483  assert(cutinds[i] < firstcontvar);
3484  ++i;
3485  }
3486 #endif
3487 
3488  for( i = *nnz - 1; i >= 0 && cutinds[i] < firstcontvar; --i )
3489  {
3490  SCIP_VAR* var;
3491  SCIP_Real QUAD(cutaj);
3492  int v;
3493 
3494  v = cutinds[i];
3495  assert(0 <= v && v < SCIPgetNVars(scip));
3496 
3497  var = vars[v];
3498  assert(var != NULL);
3499  assert(SCIPvarGetProbindex(var) == v);
3500  assert(varsign[i] == +1 || varsign[i] == -1);
3501 
3502  /* calculate the coefficient in the retransformed cut */
3503  {
3504  SCIP_Real QUAD(aj);
3505  SCIP_Real QUAD(downaj);
3506  SCIP_Real QUAD(fj);
3507 
3508  QUAD_ARRAY_LOAD(aj, cutcoefs, v);
3509  QUAD_SCALE(aj, varsign[i]);
3510 
3511  SCIPquadprecEpsFloorQ(downaj, aj, SCIPepsilon(scip)); /*lint !e666*/
3512  SCIPquadprecSumQQ(fj, aj, -downaj);
3513  assert(QUAD_TO_DBL(fj) >= -SCIPepsilon(scip) && QUAD_TO_DBL(fj) < 1.0);
3514 
3515  if( SCIPisLE(scip, QUAD_TO_DBL(fj), QUAD_TO_DBL(f0)) )
3516  {
3517  QUAD_ASSIGN_Q(cutaj, downaj);
3518  }
3519  else
3520  {
3521  SCIPquadprecSumQQ(tmp, fj, -f0);
3522  SCIPquadprecProdQQ(tmp, tmp, onedivoneminusf0);
3523  SCIPquadprecSumQQ(cutaj, tmp, downaj);
3524  }
3525 
3526  QUAD_SCALE(cutaj, varsign[i]);
3527  }
3528 
3529  /* remove zero cut coefficients from cut */
3530  if( EPSZ(QUAD_TO_DBL(cutaj), QUAD_EPSILON) )
3531  {
3532  QUAD_ASSIGN(cutaj, 0.0);
3533  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3534  --*nnz;
3535  cutinds[i] = cutinds[*nnz];
3536  continue;
3537  }
3538 
3539  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3540 
3541  /* integral var uses standard bound */
3542  assert(boundtype[i] < 0);
3543 
3544  /* move the constant term -a~_j * lb_j == -a^_j * lb_j , or a~_j * ub_j == -a^_j * ub_j to the rhs */
3545  if( varsign[i] == +1 )
3546  {
3547  /* lower bound was used */
3548  if( boundtype[i] == -1 )
3549  {
3550  assert(!SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)));
3551  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbGlobal(var));
3552  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetLbGlobal(var) */
3553  }
3554  else
3555  {
3556  assert(!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)));
3557  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbLocal(var));
3558  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetLbLocal(var) */
3559  }
3560  }
3561  else
3562  {
3563  /* upper bound was used */
3564  if( boundtype[i] == -1 )
3565  {
3566  assert(!SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)));
3567  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbGlobal(var));
3568  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetUbGlobal(var) */
3569  }
3570  else
3571  {
3572  assert(!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)));
3573  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbLocal(var));
3574  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp); /* rhs += cutaj * SCIPvarGetUbLocal(var) */
3575  }
3576  }
3577  }
3578 
3579  /* now process the continuous variables; postpone deletetion of zeros till all continuous variables have been processed */
3580  ndelcontvars = 0;
3581  while( i >= ndelcontvars )
3582  {
3583  SCIP_VAR* var;
3584  SCIP_Real QUAD(cutaj);
3585  int v;
3586 
3587  v = cutinds[i];
3588  assert(0 <= v && v < SCIPgetNVars(scip));
3589 
3590  var = vars[v];
3591  assert(var != NULL);
3592  assert(SCIPvarGetProbindex(var) == v);
3593  assert(varsign[i] == +1 || varsign[i] == -1);
3594  assert( v >= firstcontvar );
3595 
3596  /* calculate the coefficient in the retransformed cut */
3597  {
3598  SCIP_Real QUAD(aj);
3599 
3600  QUAD_ARRAY_LOAD(aj, cutcoefs, v);
3601 
3602  if( QUAD_TO_DBL(aj) * varsign[i] >= 0.0 )
3603  QUAD_ASSIGN(cutaj, 0.0);
3604  else
3605  SCIPquadprecProdQQ(cutaj, onedivoneminusf0, aj); /* cutaj = varsign[i] * aj * onedivoneminusf0; // a^_j */
3606  }
3607 
3608  /* remove zero cut coefficients from cut; move a continuous var from the beginning
3609  * to the current position, so that all integral variables stay behind the continuous
3610  * variables
3611  */
3612  if( EPSZ(QUAD_TO_DBL(cutaj), QUAD_EPSILON) )
3613  {
3614  QUAD_ASSIGN(cutaj, 0.0);
3615  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3616  cutinds[i] = cutinds[ndelcontvars];
3617  varsign[i] = varsign[ndelcontvars];
3618  boundtype[i] = boundtype[ndelcontvars];
3619  ++ndelcontvars;
3620  continue;
3621  }
3622 
3623  QUAD_ARRAY_STORE(cutcoefs, v, cutaj);
3624 
3625  /* check for variable bound use */
3626  if( boundtype[i] < 0 )
3627  {
3628  /* standard bound */
3629 
3630  /* move the constant term -a~_j * lb_j == -a^_j * lb_j , or a~_j * ub_j == -a^_j * ub_j to the rhs */
3631  if( varsign[i] == +1 )
3632  {
3633  /* lower bound was used */
3634  if( boundtype[i] == -1 )
3635  {
3636  assert(!SCIPisInfinity(scip, -SCIPvarGetLbGlobal(var)));
3637  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbGlobal(var));
3638  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3639  }
3640  else
3641  {
3642  assert(!SCIPisInfinity(scip, -SCIPvarGetLbLocal(var)));
3643  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetLbLocal(var));
3644  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3645  }
3646  }
3647  else
3648  {
3649  /* upper bound was used */
3650  if( boundtype[i] == -1 )
3651  {
3652  assert(!SCIPisInfinity(scip, SCIPvarGetUbGlobal(var)));
3653  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbGlobal(var));
3654  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3655  }
3656  else
3657  {
3658  assert(!SCIPisInfinity(scip, SCIPvarGetUbLocal(var)));
3659  SCIPquadprecProdQD(tmp, cutaj, SCIPvarGetUbLocal(var));
3660  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3661  }
3662  }
3663  }
3664  else
3665  {
3666  SCIP_VAR** vbz;
3667  SCIP_Real* vbb;
3668  SCIP_Real* vbd;
3669  SCIP_Real QUAD(zcoef);
3670  int vbidx;
3671  int zidx;
3672 
3673  /* variable bound */
3674  vbidx = boundtype[i];
3675 
3676  /* change mirrhs and cutaj of integer variable z_j of variable bound */
3677  if( varsign[i] == +1 )
3678  {
3679  /* variable lower bound was used */
3680  assert(0 <= vbidx && vbidx < SCIPvarGetNVlbs(var));
3681  vbz = SCIPvarGetVlbVars(var);
3682  vbb = SCIPvarGetVlbCoefs(var);
3683  vbd = SCIPvarGetVlbConstants(var);
3684  }
3685  else
3686  {
3687  /* variable upper bound was used */
3688  assert(0 <= vbidx && vbidx < SCIPvarGetNVubs(var));
3689  vbz = SCIPvarGetVubVars(var);
3690  vbb = SCIPvarGetVubCoefs(var);
3691  vbd = SCIPvarGetVubConstants(var);
3692  }
3693  assert(SCIPvarIsActive(vbz[vbidx]));
3694  zidx = SCIPvarGetProbindex(vbz[vbidx]);
3695  assert(0 <= zidx && zidx < firstcontvar);
3696 
3697  SCIPquadprecProdQD(tmp, cutaj, vbd[vbidx]);
3698  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3699 
3700  SCIPquadprecProdQD(tmp, cutaj, vbb[vbidx]);
3701  QUAD_ARRAY_LOAD(zcoef, cutcoefs, zidx);
3702 
3703  /* update sparsity pattern */
3704  if( QUAD_HI(zcoef) == 0.0 )
3705  cutinds[(*nnz)++] = zidx;
3706 
3707  SCIPquadprecSumQQ(zcoef, zcoef, -tmp);
3708  QUAD_HI(zcoef) = NONZERO(QUAD_HI(zcoef));
3709  QUAD_ARRAY_STORE(cutcoefs, zidx, zcoef);
3710  assert(QUAD_HI(zcoef) != 0.0);
3711  }
3712 
3713  /* advance to next variable */
3714  --i;
3715  }
3716 
3717  /* fill the empty position due to deleted continuous variables */
3718  if( ndelcontvars > 0 )
3719  {
3720  assert(ndelcontvars <= *nnz);
3721  *nnz -= ndelcontvars;
3722  if( *nnz < ndelcontvars )
3723  {
3724  BMScopyMemoryArray(cutinds, cutinds + ndelcontvars, *nnz);
3725  }
3726  else
3727  {
3728  BMScopyMemoryArray(cutinds, cutinds + *nnz, ndelcontvars);
3729  }
3730  }
3731 
3732  return SCIP_OKAY;
3733 }
3734 
3735 /** substitute aggregated slack variables:
3736  *
3737  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
3738  * variable only appears in its own row: \f$ a^\prime_r = scale * weight[r] * slacksign[r]. \f$
3739  *
3740  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
3741  * \f[
3742  * \begin{array}{rll}
3743  * integers : & \hat{a}_r = \tilde{a}_r = down(a^\prime_r), & \mbox{if}\qquad f_r <= f0 \\
3744  * & \hat{a}_r = \tilde{a}_r = down(a^\prime_r) + (f_r - f0)/(1 - f0),& \mbox{if}\qquad f_r > f0 \\
3745  * continuous:& \hat{a}_r = \tilde{a}_r = 0, & \mbox{if}\qquad a^\prime_r >= 0 \\
3746  * & \hat{a}_r = \tilde{a}_r = a^\prime_r/(1 - f0), & \mbox{if}\qquad a^\prime_r < 0
3747  * \end{array}
3748  * \f]
3749  *
3750  * Substitute \f$ \hat{a}_r \cdot s_r \f$ by adding \f$ \hat{a}_r \f$ times the slack's definition to the cut.
3751  */
3752 static
3754  SCIP* scip, /**< SCIP data structure */
3755  SCIP_Real* weights, /**< row weights in row summation */
3756  int* slacksign, /**< stores the sign of the row's slack variable in summation */
3757  int* rowinds, /**< sparsity pattern of used rows */
3758  int nrowinds, /**< number of used rows */
3759  SCIP_Real scale, /**< additional scaling factor multiplied to all rows */
3760  SCIP_Real* cutcoefs, /**< array of coefficients of cut */
3761  QUAD(SCIP_Real* cutrhs), /**< pointer to right hand side of cut */
3762  int* cutinds, /**< array of variables problem indices for non-zero coefficients in cut */
3763  int* nnz, /**< number of non-zeros in cut */
3764  QUAD(SCIP_Real f0) /**< fractional value of rhs */
3765  )
3766 { /*lint --e{715}*/
3767  SCIP_ROW** rows;
3768  SCIP_Real QUAD(onedivoneminusf0);
3769  int i;
3770 
3771  assert(scip != NULL);
3772  assert(weights != NULL || nrowinds == 0);
3773  assert(slacksign != NULL || nrowinds == 0);
3774  assert(rowinds != NULL || nrowinds == 0);
3775  assert(scale > 0.0);
3776  assert(cutcoefs != NULL);
3777  assert(QUAD_HI(cutrhs) != NULL);
3778  assert(cutinds != NULL);
3779  assert(nnz != NULL);
3780  assert(0.0 < QUAD_TO_DBL(f0) && QUAD_TO_DBL(f0) < 1.0);
3781 
3782  SCIPquadprecSumQD(onedivoneminusf0, -f0, 1.0);
3783  SCIPquadprecDivDQ(onedivoneminusf0, 1.0, onedivoneminusf0);
3784 
3785  rows = SCIPgetLPRows(scip);
3786  for( i = 0; i < nrowinds; i++ )
3787  {
3788  SCIP_ROW* row;
3789  SCIP_Real ar;
3790  SCIP_Real downar;
3791  SCIP_Real QUAD(cutar);
3792  SCIP_Real QUAD(fr);
3793  SCIP_Real QUAD(tmp);
3794  SCIP_Real mul;
3795  int r;
3796 
3797  r = rowinds[i]; /*lint !e613*/
3798  assert(0 <= r && r < SCIPgetNLPRows(scip));
3799  assert(slacksign[i] == -1 || slacksign[i] == +1); /*lint !e613*/
3800  assert(!SCIPisZero(scip, weights[i])); /*lint !e613*/
3801 
3802  row = rows[r];
3803  assert(row != NULL);
3804  assert(row->len == 0 || row->cols != NULL);
3805  assert(row->len == 0 || row->cols_index != NULL);
3806  assert(row->len == 0 || row->vals != NULL);
3807 
3808  /* get the slack's coefficient a'_r in the aggregated row */
3809  ar = slacksign[i] * scale * weights[i]; /*lint !e613*/
3810 
3811  /* calculate slack variable's coefficient a^_r in the cut */
3812  if( row->integral
3813  && ((slacksign[i] == +1 && SCIPisFeasIntegral(scip, row->rhs - row->constant))
3814  || (slacksign[i] == -1 && SCIPisFeasIntegral(scip, row->lhs - row->constant))) ) /*lint !e613*/
3815  {
3816  /* slack variable is always integral:
3817  * a^_r = a~_r = down(a'_r) , if f_r <= f0
3818  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
3819  */
3820  downar = EPSFLOOR(ar, QUAD_EPSILON);
3821  SCIPquadprecSumDD(fr, ar, -downar);
3822  if( SCIPisLE(scip, QUAD_TO_DBL(fr), QUAD_TO_DBL(f0)) )
3823  QUAD_ASSIGN(cutar, downar);
3824  else
3825  {
3826  SCIPquadprecSumQQ(cutar, fr, -f0);
3827  SCIPquadprecProdQQ(cutar, cutar, onedivoneminusf0);
3828  SCIPquadprecSumQD(cutar, cutar, downar);
3829  }
3830  }
3831  else
3832  {
3833  /* slack variable is continuous:
3834  * a^_r = a~_r = 0 , if a'_r >= 0
3835  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
3836  */
3837  if( ar >= 0.0 )
3838  continue; /* slack can be ignored, because its coefficient is reduced to 0.0 */
3839  else
3840  SCIPquadprecProdQD(cutar, onedivoneminusf0, ar);
3841  }
3842 
3843  /* if the coefficient was reduced to zero, ignore the slack variable */
3844  if( EPSZ(QUAD_TO_DBL(cutar), QUAD_EPSILON) )
3845  continue;
3846 
3847  /* depending on the slack's sign, we have
3848  * a*x + c + s == rhs => s == - a*x - c + rhs, or a*x + c - s == lhs => s == a*x + c - lhs
3849  * substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
3850  */
3851  mul = -slacksign[i] * QUAD_TO_DBL(cutar); /*lint !e613*/
3852 
3853  /* add the slack's definition multiplied with a^_j to the cut */
3854  SCIP_CALL( varVecAddScaledRowCoefsQuad(cutinds, cutcoefs, nnz, row, mul) );
3855 
3856  /* move slack's constant to the right hand side */
3857  if( slacksign[i] == +1 ) /*lint !e613*/
3858  {
3859  SCIP_Real QUAD(rowrhs);
3860 
3861  /* a*x + c + s == rhs => s == - a*x - c + rhs: move a^_r * (rhs - c) to the right hand side */
3862  assert(!SCIPisInfinity(scip, row->rhs));
3863  SCIPquadprecSumDD(rowrhs, row->rhs, -row->constant);
3864  if( row->integral )
3865  {
3866  /* the right hand side was implicitly rounded down in row aggregation */
3867  QUAD_ASSIGN(rowrhs, SCIPfloor(scip, QUAD_TO_DBL(rowrhs)));
3868  }
3869  SCIPquadprecProdQQ(tmp, cutar, rowrhs);
3870  SCIPquadprecSumQQ(*cutrhs, *cutrhs, -tmp);
3871  }
3872  else
3873  {
3874  SCIP_Real QUAD(rowlhs);
3875 
3876  /* a*x + c - s == lhs => s == a*x + c - lhs: move a^_r * (c - lhs) to the right hand side */
3877  assert(!SCIPisInfinity(scip, -row->lhs));
3878  SCIPquadprecSumDD(rowlhs, row->lhs, -row->constant);
3879  if( row->integral )
3880  {
3881  /* the left hand side was implicitly rounded up in row aggregation */
3882  QUAD_ASSIGN(rowlhs, SCIPceil(scip, QUAD_TO_DBL(rowlhs)));
3883  }
3884  SCIPquadprecProdQQ(tmp, cutar, rowlhs);
3885  SCIPquadprecSumQQ(*cutrhs, *cutrhs, tmp);
3886  }
3887  }
3888 
3889  /* relax rhs to zero, if it's very close to */
3890  if( QUAD_TO_DBL(*cutrhs) < 0.0 && QUAD_TO_DBL(*cutrhs) >= SCIPepsilon(scip) )
3891  QUAD_ASSIGN(*cutrhs, 0.0);
3892 
3893  return SCIP_OKAY;
3894 }
3895 
3896 /** calculates an MIR cut out of the weighted sum of LP rows; The weights of modifiable rows are set to 0.0, because
3897  * these rows cannot participate in an MIR cut.
3898  *
3899  * @return \ref SCIP_OKAY is returned if everything worked. Otherwise a suitable error code is passed. See \ref
3900  * SCIP_Retcode "SCIP_RETCODE" for a complete list of error codes.
3901  *
3902  * @pre This method can be called if @p scip is in one of the following stages:
3903  * - \ref SCIP_STAGE_SOLVING
3904  *
3905  * See \ref SCIP_Stage "SCIP_STAGE" for a complete list of all possible solving stages.
3906  */
3908  SCIP* scip, /**< SCIP data structure */
3909  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
3910  SCIP_Bool postprocess, /**< apply a post-processing step to the resulting cut? */
3911  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
3912  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
3913  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
3914  SCIP_Bool fixintegralrhs, /**< should complementation tried to be adjusted such that rhs gets fractional? */
3915  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
3916  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
3917  * NULL for using closest bound for all variables */
3918  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
3919  * NULL for using closest bound for all variables */
3920  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
3921  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
3922  SCIP_Real scale, /**< additional scaling factor multiplied to the aggrrow; must be positive */
3923  SCIP_AGGRROW* aggrrow, /**< aggrrow to compute MIR cut for */
3924  SCIP_Real* cutcoefs, /**< array to store the non-zero coefficients in the cut */
3925  SCIP_Real* cutrhs, /**< pointer to store the right hand side of the cut */
3926  int* cutinds, /**< array to store the problem indices of variables with a non-zero coefficient in the cut */
3927  int* cutnnz, /**< pointer to store the number of non-zeros in the cut */
3928  SCIP_Real* cutefficacy, /**< pointer to store efficacy of cut, or NULL */
3929  int* cutrank, /**< pointer to return rank of generated cut */
3930  SCIP_Bool* cutislocal, /**< pointer to store whether the generated cut is only valid locally */
3931  SCIP_Bool* success /**< pointer to store whether the returned coefficients are a valid MIR cut */
3932  )
3933 {
3934  int i;
3935  int nvars;
3936  int* varsign;
3937  int* boundtype;
3938  SCIP_Real* tmpcoefs;
3939 
3940  SCIP_Real QUAD(rhs);
3941  SCIP_Real QUAD(downrhs);
3942  SCIP_Real QUAD(f0);
3943  SCIP_Bool freevariable;
3944  SCIP_Bool localbdsused;
3945 
3946  assert(aggrrow != NULL);
3947  assert(SCIPisPositive(scip, scale));
3948  assert(success != NULL);
3949 
3950  SCIPdebugMessage("calculating MIR cut (scale: %g)\n", scale);
3951 
3952  *success = FALSE;
3953 
3954  /* allocate temporary memory */
3955  nvars = SCIPgetNVars(scip);
3956  SCIP_CALL( SCIPallocBufferArray(scip, &varsign, nvars) );
3957  SCIP_CALL( SCIPallocBufferArray(scip, &boundtype, nvars) );
3958  SCIP_CALL( SCIPallocCleanBufferArray(scip, &tmpcoefs, QUAD_ARRAY_SIZE(nvars)) );
3959 
3960  /* initialize cut with aggregation */
3961  *cutnnz = aggrrow->nnz;
3962  *cutislocal = aggrrow->local;
3963 
3964  SCIPquadprecProdQD(rhs, aggrrow->rhs, scale);
3965 
3966  if( *cutnnz > 0 )
3967  {
3968  BMScopyMemoryArray(cutinds, aggrrow->inds, *cutnnz);
3969 
3970  for( i = 0; i < *cutnnz; ++i )
3971  {
3972  SCIP_Real QUAD(coef);
3973 
3974  int k = aggrrow->inds[i];
3975  QUAD_ARRAY_LOAD(coef, aggrrow->vals, k);
3976 
3977  SCIPquadprecProdQD(coef, coef, scale);
3978 
3979  QUAD_ARRAY_STORE(tmpcoefs, k, coef);
3980 
3981  assert(QUAD_HI(coef) != 0.0);
3982  }
3983 
3984  /* Transform equation a*x == b, lb <= x <= ub into standard form
3985  * a'*x' == b, 0 <= x' <= ub'.
3986  *
3987  * Transform variables (lb or ub):
3988  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, if lb is used in transformation
3989  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, if ub is used in transformation
3990  * and move the constant terms "a_j * lb_j" or "a_j * ub_j" to the rhs.
3991  *
3992  * Transform variables (vlb or vub):
3993  * 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.
3994  * 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.
3995  * move the constant terms "a_j * dl_j" or "a_j * du_j" to the rhs, and update the coefficient of the VLB variable:
3996  * a_{zl_j} := a_{zl_j} + a_j * bl_j, or
3997  * a_{zu_j} := a_{zu_j} + a_j * bu_j
3998  */
3999  SCIP_CALL( cutsTransformMIR(scip, sol, boundswitch, usevbds, allowlocal, fixintegralrhs, FALSE,
4000  boundsfortrans, boundtypesfortrans, minfrac, maxfrac, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, varsign, boundtype, &freevariable, &localbdsused) );
4001  assert(allowlocal || !localbdsused);
4002  *cutislocal = *cutislocal || localbdsused;
4003 
4004  if( freevariable )
4005  goto TERMINATE;
4006  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4007  }
4008 
4009  /* Calculate fractionalities f_0 := b - down(b), f_j := a'_j - down(a'_j) , and derive MIR cut
4010  * a~*x' <= down(b)
4011  * integers : a~_j = down(a'_j) , if f_j <= f_0
4012  * a~_j = down(a'_j) + (f_j - f0)/(1 - f0), if f_j > f_0
4013  * continuous: a~_j = 0 , if a'_j >= 0
4014  * a~_j = a'_j/(1 - f0) , if a'_j < 0
4015  *
4016  * Transform inequality back to a^*x <= rhs:
4017  *
4018  * (lb or ub):
4019  * 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
4020  * 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
4021  * and move the constant terms
4022  * -a~_j * lb_j == -a^_j * lb_j, or
4023  * a~_j * ub_j == -a^_j * ub_j
4024  * to the rhs.
4025  *
4026  * (vlb or vub):
4027  * 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)
4028  * 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)
4029  * move the constant terms
4030  * -a~_j * dl_j == -a^_j * dl_j, or
4031  * a~_j * du_j == -a^_j * du_j
4032  * to the rhs, and update the VB variable coefficients:
4033  * a^_{zl_j} := a^_{zl_j} - a~_j * bl_j == a^_{zl_j} - a^_j * bl_j, or
4034  * a^_{zu_j} := a^_{zu_j} + a~_j * bu_j == a^_{zu_j} - a^_j * bu_j
4035  */
4036  SCIPquadprecEpsFloorQ(downrhs, rhs, SCIPepsilon(scip)); /*lint !e666*/
4037 
4038  SCIPquadprecSumQQ(f0, rhs, -downrhs);
4039 
4040  if( QUAD_TO_DBL(f0) < minfrac || QUAD_TO_DBL(f0) > maxfrac )
4041  goto TERMINATE;
4042 
4043  /* We multiply the coefficients of the base inequality roughly by scale/(1-f0).
4044  * If this gives a scalar that is very big, we better do not generate this cut.
4045  */
4046  if( REALABS(scale)/(1.0 - QUAD_TO_DBL(f0)) > MAXCMIRSCALE )
4047  goto TERMINATE;
4048 
4049  /* renormalize f0 value */
4050  SCIPquadprecSumDD(f0, QUAD_HI(f0), QUAD_LO(f0));
4051 
4052  QUAD_ASSIGN_Q(rhs, downrhs);
4053 
4054  if( *cutnnz > 0 )
4055  {
4056  SCIP_CALL( cutsRoundMIR(scip, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, varsign, boundtype, QUAD(f0)) );
4057  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4058  }
4059 
4060  /* substitute aggregated slack variables:
4061  *
4062  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
4063  * variable only appears in its own row:
4064  * a'_r = scale * weight[r] * slacksign[r].
4065  *
4066  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
4067  * integers : a^_r = a~_r = down(a'_r) , if f_r <= f0
4068  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
4069  * continuous: a^_r = a~_r = 0 , if a'_r >= 0
4070  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
4071  *
4072  * Substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
4073  */
4074  SCIP_CALL( cutsSubstituteMIR(scip, aggrrow->rowweights, aggrrow->slacksign, aggrrow->rowsinds,
4075  aggrrow->nrows, scale, tmpcoefs, QUAD(&rhs), cutinds, cutnnz, QUAD(f0)) );
4076  SCIPdebug( printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE) );
4077 
4078  if( postprocess )
4079  {
4080  /* remove all nearly-zero coefficients from MIR row and relax the right hand side correspondingly in order to
4081  * prevent numerical rounding errors
4082  */
4083  SCIP_CALL( postprocessCutQuad(scip, *cutislocal, cutinds, tmpcoefs, cutnnz, QUAD(&rhs), success) );
4084  }
4085  else
4086  {
4087  *success = ! removeZerosQuad(scip, SCIPsumepsilon(scip), *cutislocal, tmpcoefs, QUAD(&rhs), cutnnz, cutinds);
4088  }
4089 
4090  SCIPdebug(printCutQuad(scip, sol, tmpcoefs, QUAD(rhs), cutinds, *cutnnz, FALSE, FALSE));
4091 
4092  if( *success )
4093  {
4094  *cutrhs = QUAD_TO_DBL(rhs);
4095 
4096  /* clean tmpcoefs and go back to double precision */
4097  for( i = 0; i < *cutnnz; ++i )
4098  {
4099  SCIP_Real QUAD(coef);
4100  int j = cutinds[i];
4101 
4102  QUAD_ARRAY_LOAD(coef, tmpcoefs, j);
4103 
4104  cutcoefs[i] = QUAD_TO_DBL(coef);
4105  QUAD_ASSIGN(coef, 0.0);
4106  QUAD_ARRAY_STORE(tmpcoefs, j, coef);
4107  }
4108 
4109  if( cutefficacy != NULL )
4110  *cutefficacy = calcEfficacy(scip, sol, cutcoefs, *cutrhs, cutinds, *cutnnz);
4111 
4112  if( cutrank != NULL )
4113  *cutrank = aggrrow->rank + 1;
4114  }
4115 
4116  TERMINATE:
4117  if( !(*success) )
4118  {
4119  SCIP_Real QUAD(tmp);
4120 
4121  QUAD_ASSIGN(tmp, 0.0);
4122  for( i = 0; i < *cutnnz; ++i )
4123  {
4124  QUAD_ARRAY_STORE(tmpcoefs, cutinds[i], tmp);
4125  }
4126  }
4127 
4128  /* free temporary memory */
4129  SCIPfreeCleanBufferArray(scip, &tmpcoefs);
4130  SCIPfreeBufferArray(scip, &boundtype);
4131  SCIPfreeBufferArray(scip, &varsign);
4132 
4133  return SCIP_OKAY;
4134 }
4135 
4136 /** compute the efficacy of the MIR cut for the given values without computing the cut.
4137  * This is used for the CMIR cut generation heuristic.
4138  */
4139 static
4141  SCIP* scip, /**< SCIP datastructure */
4142  SCIP_Real*RESTRICT coefs, /**< array with coefficients in row */
4143  SCIP_Real*RESTRICT solvals, /**< solution values of variables in the row */
4144  SCIP_Real rhs, /**< right hand side of MIR cut */
4145  SCIP_Real contactivity, /**< aggregated activity of continuous variables in the row */
4146  SCIP_Real contsqrnorm, /**< squared norm of continuous variables */
4147  SCIP_Real delta, /**< delta value to compute the violation for */
4148  int nvars, /**< number of variables in the row, i.e. the size of coefs and solvals arrays */
4149  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
4150  SCIP_Real maxfrac /**< maximal fractionality of rhs to produce MIR cut for */
4151  )
4152 {
4153  int i;
4154  SCIP_Real f0pluseps;
4155  SCIP_Real f0;
4156  SCIP_Real onedivoneminusf0;
4157  SCIP_Real scale;
4158  SCIP_Real downrhs;
4159  SCIP_Real norm;
4160  SCIP_Real contscale;
4161 
4162  scale = 1.0 / delta;
4163 
4164  rhs *= scale;
4165 
4166  downrhs = SCIPfloor(scip, rhs);
4167 
4168  f0 = rhs - downrhs;
4169 
4170  if( f0 < minfrac || f0 > maxfrac )
4171  return 0.0;
4172 
4173  onedivoneminusf0 = 1.0 / (1.0 - f0);
4174 
4175  contscale = scale * onedivoneminusf0;
4176 
4177  /* We multiply the coefficients of the base inequality roughly by scale/(1-f0).
4178  * If this gives a scalar that is very big, we better do not generate this cut.
4179  */
4180  if( contscale > MAXCMIRSCALE )
4181  return 0.0;
4182 
4183  rhs = downrhs;
4184  rhs -= contscale * contactivity;
4185  norm = SQR(contscale) * contsqrnorm;
4186 
4187  assert(!SCIPisFeasZero(scip, f0));
4188  assert(!SCIPisFeasZero(scip, 1.0 - f0));
4189 
4190  f0pluseps = f0 + SCIPepsilon(scip);
4191 
4192  for( i = 0; i < nvars; ++i )
4193  {
4194  SCIP_Real floorai = floor(scale * coefs[i]);
4195  SCIP_Real fi = (scale * coefs[i]) - floorai;
4196 
4197  if( fi > f0pluseps )
4198  floorai += (fi - f0) * onedivoneminusf0;
4199 
4200  rhs -= solvals[i] * floorai;
4201  norm += SQR(floorai);
4202  }
4203 
4204  norm = SQRT(norm);
4205 
4206  return - rhs / MAX(norm, 1e-6);
4207 }
4208 
4209 /** calculates an MIR cut out of an aggregation of LP rows
4210  *
4211  * Given the aggregation, it is transformed to a mixed knapsack set via complementation (using bounds or variable bounds)
4212  * Then, different scalings of the mkset are used to generate a MIR and the best is chosen.
4213  * One of the steps of the MIR is to round the coefficients of the integer variables down,
4214  * so one would prefer to have integer coefficients for integer variables which are far away from their bounds in the
4215  * mkset.
4216  *
4217  * @return \ref SCIP_OKAY is returned if everything worked. Otherwise a suitable error code is passed. See \ref
4218  * SCIP_Retcode "SCIP_RETCODE" for a complete list of error codes.
4219  *
4220  * @pre This method can be called if @p scip is in one of the following stages:
4221  * - \ref SCIP_STAGE_SOLVING
4222  *
4223  * See \ref SCIP_Stage "SCIP_STAGE" for a complete list of all possible solving stages.
4224  */
4226  SCIP* scip, /**< SCIP data structure */
4227  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
4228  SCIP_Bool postprocess, /**< apply a post-processing step to the resulting cut? */
4229  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
4230  SCIP_Bool usevbds, /**< should variable bounds be used in bound transformation? */
4231  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
4232  int maxtestdelta, /**< maximum number of deltas to test */
4233  int* boundsfortrans, /**< bounds that should be used for transformed variables: vlb_idx/vub_idx,
4234  * -1 for global lb/ub, -2 for local lb/ub, or -3 for using closest bound;
4235  * NULL for using closest bound for all variables */
4236  SCIP_BOUNDTYPE* boundtypesfortrans, /**< type of bounds that should be used for transformed variables;
4237  * NULL for using closest bound for all variables */
4238  SCIP_Real minfrac, /**< minimal fractionality of rhs to produce MIR cut for */
4239  SCIP_Real maxfrac, /**< maximal fractionality of rhs to produce MIR cut for */
4240  SCIP_AGGRROW* aggrrow, /**< aggrrow to compute MIR cut for */
4241  SCIP_Real* cutcoefs, /**< array to store the non-zero coefficients in the cut */
4242  SCIP_Real* cutrhs, /**< pointer to store the right hand side of the cut */
4243  int* cutinds, /**< array to store the problem indices of variables with a non-zero coefficient in the cut */
4244  int* cutnnz, /**< pointer to store the number of non-zeros in the cut */
4245  SCIP_Real* cutefficacy, /**< pointer to store efficacy of best cut; only cuts that are strictly better than the value of
4246  * this efficacy on input to this function are returned */
4247  int* cutrank, /**< pointer to return rank of generated cut (or NULL) */
4248  SCIP_Bool* cutislocal, /**< pointer to store whether the generated cut is only valid locally */
4249  SCIP_Bool* success /**< pointer to store whether a valid and efficacious cut was returned */
4250  )
4251 {
4252  int i;
4253  int firstcontvar;
4254  int nvars;
4255  int intstart;
4256  int ntmpcoefs;
4257  int* varsign;
4258  int* boundtype;
4259  int* mksetinds;
4260  SCIP_Real* mksetcoefs;
4261  SCIP_Real QUAD(mksetrhs);
4262  int mksetnnz;
4263  SCIP_Real* bounddist;
4264  int* bounddistpos;
4265  int nbounddist;
4266  SCIP_Real* tmpcoefs;
4267  SCIP_Real* tmpvalues;
4268  SCIP_Real* deltacands;
4269  int ndeltacands;
4270  SCIP_Real bestdelta;
4271  SCIP_Real bestefficacy;
4272  SCIP_Real maxabsmksetcoef;
4273  SCIP_VAR** vars;
4274  SCIP_Bool freevariable;
4275  SCIP_Bool localbdsused;
4276  SCIP_Real contactivity;
4277  SCIP_Real contsqrnorm;
4278 
4279  assert(aggrrow != NULL);
4280  assert(aggrrow->nrows + aggrrow->nnz >= 1);
4281  assert(success != NULL);
4282 
4283  *success = FALSE;
4284  nvars = SCIPgetNVars(scip);
4285  firstcontvar = nvars - SCIPgetNContVars(scip);
4286  vars = SCIPgetVars(scip);
4287 
4288  /* allocate temporary memory */
4289  SCIP_CALL( SCIPallocBufferArray(scip, &varsign, nvars) );
4290  SCIP_CALL( SCIPallocBufferArray(scip, &boundtype, nvars) );
4291  SCIP_CALL( SCIPallocCleanBufferArray(scip, &mksetcoefs, QUAD_ARRAY_SIZE(nvars)) );
4292  SCIP_CALL( SCIPallocBufferArray(scip, &mksetinds, nvars) );
4293  SCIP_CALL( SCIPallocBufferArray(scip, &tmpcoefs, nvars + aggrrow->nrows) );
4294  SCIP_CALL( SCIPallocBufferArray(scip, &tmpvalues, nvars + aggrrow->nrows) );
4295  SCIP_CALL( SCIPallocBufferArray(scip, &deltacands, aggrrow->nnz + 6) );
4296 
4297  /* we only compute bound distance for integer variables; we allocate an array of length aggrrow->nnz to store this, since
4298  * this is the largest number of integer variables. (in contrast to the number of total variables which can be 2 *
4299  * aggrrow->nnz variables: if all are continuous and we use variable bounds to completement, we introduce aggrrow->nnz
4300  * extra vars)
4301  */
4302  SCIP_CALL( SCIPallocBufferArray(scip, &bounddist, aggrrow->nnz) );
4303  SCIP_CALL( SCIPallocBufferArray(scip, &bounddistpos, aggrrow->nnz) );
4304 
4305  /* initialize mkset with aggregation */
4306  mksetnnz = aggrrow->nnz;
4307  QUAD_ASSIGN_Q(mksetrhs, aggrrow->rhs);
4308 
4309  BMScopyMemoryArray(mksetinds, aggrrow->inds, mksetnnz);
4310 
4311  for( i = 0; i < mksetnnz; ++i )
4312  {
4313  int j = mksetinds[i];
4314  SCIP_Real QUAD(coef);
4315  QUAD_ARRAY_LOAD(coef, aggrrow->vals, j);
4316  QUAD_ARRAY_STORE(mksetcoefs, j, coef);
4317  assert(QUAD_HI(coef) != 0.0);
4318  }
4319 
4320  *cutislocal = aggrrow->local;
4321 
4322  /* Transform equation a*x == b, lb <= x <= ub into standard form
4323  * a'*x' == b, 0 <= x' <= ub'.
4324  *
4325  * Transform variables (lb or ub):
4326  * x'_j := x_j - lb_j, x_j == x'_j + lb_j, a'_j == a_j, if lb is used in transformation
4327  * x'_j := ub_j - x_j, x_j == ub_j - x'_j, a'_j == -a_j, if ub is used in transformation
4328  * and move the constant terms "a_j * lb_j" or "a_j * ub_j" to the rhs.
4329  *
4330  * Transform variables (vlb or vub):
4331  * 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.
4332  * 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.
4333  * move the constant terms "a_j * dl_j" or "a_j * du_j" to the rhs, and update the coefficient of the VLB variable:
4334  * a_{zl_j} := a_{zl_j} + a_j * bl_j, or
4335  * a_{zu_j} := a_{zu_j} + a_j * bu_j
4336  */
4337  SCIP_CALL( cutsTransformMIR(scip, sol, boundswitch, usevbds, allowlocal, FALSE, FALSE,
4338  boundsfortrans, boundtypesfortrans, minfrac, maxfrac, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, varsign, boundtype, &freevariable, &localbdsused) );
4339 
4340  assert(allowlocal || !localbdsused);
4341 
4342  if( freevariable )
4343  goto TERMINATE;
4344 
4345  SCIPdebugMessage("transformed aggrrow row:\n");
4346  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4347 
4348  /* found positions of integral variables that are strictly between their bounds */
4349  maxabsmksetcoef = -1.0;
4350  nbounddist = 0;
4351 
4352  for( i = mksetnnz - 1; i >= 0 && mksetinds[i] < firstcontvar; --i )
4353  {
4354  SCIP_VAR* var = vars[mksetinds[i]];
4355  SCIP_Real primsol = SCIPgetSolVal(scip, sol, var);
4356  SCIP_Real lb = SCIPvarGetLbLocal(var);
4357  SCIP_Real ub = SCIPvarGetUbLocal(var);
4358  SCIP_Real QUAD(coef);
4359 
4360  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4361 
4362  if( SCIPisEQ(scip, primsol, lb) || SCIPisEQ(scip, primsol, ub) )
4363  continue;
4364 
4365  bounddist[nbounddist] = MIN(ub - primsol, primsol - lb);
4366  bounddistpos[nbounddist] = i;
4367  deltacands[nbounddist] = QUAD_TO_DBL(coef);
4368  ++nbounddist;
4369  }
4370 
4371  /* no fractional variable; so abort here */
4372  if( nbounddist == 0 )
4373  goto TERMINATE;
4374 
4375  intstart = i + 1;
4376  ndeltacands = nbounddist;
4377 
4378  SCIPsortDownRealRealInt(bounddist, deltacands, bounddistpos, nbounddist);
4379 
4380  {
4381  SCIP_Real intscale;
4382  SCIP_Bool intscalesuccess;
4383 
4384  SCIP_CALL( SCIPcalcIntegralScalar(deltacands, nbounddist, -QUAD_EPSILON, SCIPsumepsilon(scip), (SCIP_Longint)10000, 10000.0, &intscale, &intscalesuccess) );
4385 
4386  if( intscalesuccess )
4387  {
4388  SCIP_Real intf0;
4389  SCIP_Real intscalerhs;
4390  SCIP_Real delta;
4391 
4392  intscalerhs = QUAD_TO_DBL(mksetrhs) * intscale;
4393  delta = 1.0 / intscale;
4394  intf0 = intscalerhs - floor(intscalerhs);
4395 
4396  if( ! SCIPisFeasIntegral(scip, intf0) )
4397  {
4398  if( intf0 < minfrac || intf0 > maxfrac )
4399  {
4400  intscale *= ceil(MAX(minfrac, (1.0 - maxfrac)) / MIN(intf0, (1.0 - intf0)));
4401  intscalerhs = QUAD_TO_DBL(mksetrhs) * intscale;
4402  delta = 1.0 / intscale;
4403  intf0 = intscalerhs - floor(intscalerhs);
4404  }
4405 
4406  if( intf0 >= minfrac && intf0 <= maxfrac )
4407  {
4408  if( ! SCIPisEQ(scip, delta, 1.0) )
4409  {
4410  deltacands[ndeltacands++] = delta;
4411  }
4412 
4413  if( intf0 < maxfrac )
4414  {
4415  SCIP_Real delta2;
4416 
4417  delta2 = 1.0 / (intscale * floor(maxfrac / intf0));
4418 
4419  if( ! SCIPisEQ(scip, delta, delta2) && ! SCIPisEQ(scip, delta2, 1.0) )
4420  {
4421  deltacands[ndeltacands++] = delta2;
4422  }
4423  }
4424  }
4425  }
4426  }
4427  }
4428 
4429  for( i = 0; i < nbounddist; ++i )
4430  {
4431  SCIP_Real absmksetcoef;
4432 
4433  absmksetcoef = REALABS(deltacands[i]);
4434  maxabsmksetcoef = MAX(absmksetcoef, maxabsmksetcoef);
4435 
4436  deltacands[i] = absmksetcoef;
4437  }
4438 
4439  /* also test 1.0 and maxabsmksetcoef + 1.0 as last delta values */
4440  if( maxabsmksetcoef != -1.0 )
4441  {
4442  deltacands[ndeltacands++] = maxabsmksetcoef + 1.0;
4443  }
4444 
4445  deltacands[ndeltacands++] = 1.0;
4446 
4447  maxtestdelta = MIN(ndeltacands, maxtestdelta);
4448 
4449  /* For each delta
4450  * Calculate fractionalities f_0 := b - down(b), f_j := a'_j - down(a'_j) , and derive MIR cut
4451  * a~*x' <= down(b)
4452  * integers : a~_j = down(a'_j) , if f_j <= f_0
4453  * a~_j = down(a'_j) + (f_j - f0)/(1 - f0), if f_j > f_0
4454  * continuous: a~_j = 0 , if a'_j >= 0
4455  * a~_j = a'_j/(1 - f0) , if a'_j < 0
4456  *
4457  * Transform inequality back to a^*x <= rhs:
4458  *
4459  * (lb or ub):
4460  * 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
4461  * 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
4462  * and move the constant terms
4463  * -a~_j * lb_j == -a^_j * lb_j, or
4464  * a~_j * ub_j == -a^_j * ub_j
4465  * to the rhs.
4466  *
4467  * (vlb or vub):
4468  * 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)
4469  * 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)
4470  * move the constant terms
4471  * -a~_j * dl_j == -a^_j * dl_j, or
4472  * a~_j * du_j == -a^_j * du_j
4473  * to the rhs, and update the VB variable coefficients:
4474  * a^_{zl_j} := a^_{zl_j} - a~_j * bl_j == a^_{zl_j} - a^_j * bl_j, or
4475  * a^_{zu_j} := a^_{zu_j} + a~_j * bu_j == a^_{zu_j} - a^_j * bu_j
4476  */
4477 
4478  ntmpcoefs = 0;
4479  for( i = intstart; i < mksetnnz; ++i )
4480  {
4481  SCIP_VAR* var;
4482  SCIP_Real solval;
4483  SCIP_Real QUAD(coef);
4484 
4485  var = vars[mksetinds[i]];
4486 
4487  /* get the soltion value of the continuous variable */
4488  solval = SCIPgetSolVal(scip, sol, var);
4489 
4490  /* now compute the solution value in the transform space considering complementation */
4491  if( boundtype[i] == -1 )
4492  {
4493  /* variable was complemented with global (simple) bound */
4494  if( varsign[i] == -1 )
4495  solval = SCIPvarGetUbGlobal(var) - solval;
4496  else
4497  solval = solval - SCIPvarGetLbGlobal(var);
4498  }
4499  else
4500  {
4501  assert(boundtype[i] == -2);
4502 
4503  /* variable was complemented with local (simple) bound */
4504  if( varsign[i] == -1 )
4505  solval = SCIPvarGetUbLocal(var) - solval;
4506  else
4507  solval = solval - SCIPvarGetLbLocal(var);
4508  }
4509 
4510  tmpvalues[ntmpcoefs] = solval;
4511  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4512  tmpcoefs[ntmpcoefs] = varsign[i] * QUAD_TO_DBL(coef);
4513  ++ntmpcoefs;
4514  }
4515 
4516  assert(ntmpcoefs == mksetnnz - intstart);
4517 
4518  contactivity = 0.0;
4519  contsqrnorm = 0.0;
4520  for( i = 0; i < intstart; ++i )
4521  {
4522  SCIP_Real solval;
4523  SCIP_Real QUAD(mksetcoef);
4524 
4525  QUAD_ARRAY_LOAD(mksetcoef, mksetcoefs, mksetinds[i]);
4526 
4527  if( varsign[i] * QUAD_TO_DBL(mksetcoef) >= 0.0 )
4528  continue;
4529 
4530  /* get the soltion value of the continuous variable */
4531  solval = SCIPgetSolVal(scip, sol, vars[mksetinds[i]]);
4532 
4533  /* now compute the solution value in the transform space considering complementation */
4534  switch( boundtype[i] )
4535  {
4536  case -1:
4537  /* variable was complemented with global (simple) bound */
4538  if( varsign[i] == -1 )
4539  solval = SCIPvarGetUbGlobal(vars[mksetinds[i]]) - solval;
4540  else
4541  solval = solval - SCIPvarGetLbGlobal(vars[mksetinds[i]]);
4542  break;
4543  case -2:
4544  /* variable was complemented with local (simple) bound */
4545  if( varsign[i] == -1 )
4546  solval = SCIPvarGetUbLocal(vars[mksetinds[i]]) - solval;
4547  else
4548  solval = solval - SCIPvarGetLbLocal(vars[mksetinds[i]]);
4549  break;
4550  default:
4551  /* variable was complemented with a variable bound */
4552  if( varsign[i] == -1 )
4553  {
4554  SCIP_Real coef;
4555  SCIP_Real constant;
4556  SCIP_Real vbdsolval;
4557 
4558  coef = SCIPvarGetVubCoefs(vars[mksetinds[i]])[boundtype[i]];
4559  constant = SCIPvarGetVubConstants(vars[mksetinds[i]])[boundtype[i]];
4560  vbdsolval = SCIPgetSolVal(scip, sol, SCIPvarGetVubVars(vars[mksetinds[i]])[boundtype[i]]);
4561 
4562  solval = (coef * vbdsolval + constant) - solval;
4563  }
4564  else
4565  {
4566  SCIP_Real coef;
4567  SCIP_Real constant;
4568  SCIP_Real vbdsolval;
4569 
4570  coef = SCIPvarGetVlbCoefs(vars[mksetinds[i]])[boundtype[i]];
4571  constant = SCIPvarGetVlbConstants(vars[mksetinds[i]])[boundtype[i]];
4572  vbdsolval = SCIPgetSolVal(scip, sol, SCIPvarGetVlbVars(vars[mksetinds[i]])[boundtype[i]]);
4573 
4574  solval = solval - (coef * vbdsolval + constant);
4575  }
4576  }
4577 
4578  contactivity += solval * (QUAD_TO_DBL(mksetcoef) * varsign[i]);
4579  contsqrnorm += QUAD_TO_DBL(mksetcoef) * QUAD_TO_DBL(mksetcoef);
4580  }
4581 
4582  {
4583  SCIP_ROW** rows;
4584 
4585  rows = SCIPgetLPRows(scip);
4586 
4587  for( i = 0; i < aggrrow->nrows; ++i )
4588  {
4589  SCIP_ROW* row;
4590  SCIP_Real slackval;
4591 
4592  row = rows[aggrrow->rowsinds[i]];
4593 
4594  if( (aggrrow->rowweights[i] * aggrrow->slacksign[i]) >= 0.0 && !row->integral )
4595  continue;
4596 
4597  /* compute solution value of slack variable */
4598  slackval = SCIPgetRowSolActivity(scip, row, sol);
4599 
4600  if( aggrrow->slacksign[i] == +1 )
4601  {
4602  /* right hand side */
4603  assert(!SCIPisInfinity(scip, row->rhs));
4604 
4605  slackval = row->rhs - slackval;
4606  }
4607  else
4608  {
4609  /* left hand side */
4610  assert(aggrrow->slacksign[i] == -1);
4611  assert(!SCIPisInfinity(scip, -row->lhs));
4612 
4613  slackval = slackval - row->lhs;
4614  }
4615 
4616  if( row->integral )
4617  {
4618  /* if row is integral add variable to tmp arrays */
4619  tmpvalues[ntmpcoefs] = slackval;
4620  tmpcoefs[ntmpcoefs] = aggrrow->rowweights[i] * aggrrow->slacksign[i];
4621  ++ntmpcoefs;
4622  }
4623  else
4624  {
4625  SCIP_Real slackcoeff = (aggrrow->rowweights[i] * aggrrow->slacksign[i]);
4626 
4627  /* otherwise add it to continuous activity */
4628  contactivity += slackval * slackcoeff;
4629  contsqrnorm += SQR(slackcoeff);
4630  }
4631  }
4632  }
4633 
4634  /* try all candidates for delta and remember best */
4635  bestdelta = SCIP_INVALID;
4636  bestefficacy = -SCIPinfinity(scip);
4637 
4638  for( i = 0; i < maxtestdelta; ++i )
4639  {
4640  int j;
4641  SCIP_Real efficacy;
4642 
4643  /* check if we have seen this value of delta before */
4644  SCIP_Bool deltaseenbefore = FALSE;
4645  for( j = 0; j < i; ++j )
4646  {
4647  if( SCIPisEQ(scip, deltacands[i], deltacands[j]) )
4648  {
4649  deltaseenbefore = TRUE;
4650  break;
4651  }
4652  }
4653 
4654  /* skip this delta value and allow one more delta value if available */
4655  if( deltaseenbefore )
4656  {
4657  maxtestdelta = MIN(maxtestdelta + 1, ndeltacands);
4658  continue;
4659  }
4660 
4661  efficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(mksetrhs), contactivity, contsqrnorm, deltacands[i], ntmpcoefs, minfrac, maxfrac);
4662 
4663  if( efficacy > bestefficacy )
4664  {
4665  bestefficacy = efficacy;
4666  bestdelta = deltacands[i];
4667  }
4668  }
4669 
4670  /* no delta was found that yielded any cut */
4671  if( bestdelta == SCIP_INVALID ) /*lint !e777*/
4672  goto TERMINATE;
4673 
4674  /* try bestdelta divided by 2, 4 and 8 */
4675  {
4676  SCIP_Real basedelta = bestdelta;
4677  for( i = 2; i <= 8 ; i *= 2 )
4678  {
4679  SCIP_Real efficacy;
4680  SCIP_Real delta;
4681 
4682  delta = basedelta / i;
4683 
4684  efficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(mksetrhs), contactivity, contsqrnorm, delta, ntmpcoefs, minfrac, maxfrac);
4685 
4686  if( efficacy > bestefficacy )
4687  {
4688  bestefficacy = efficacy;
4689  bestdelta = delta;
4690  }
4691  }
4692  }
4693 
4694  /* try to improve efficacy by switching complementation of integral variables that are not at their bounds
4695  * in order of non-increasing bound distance
4696  */
4697  for( i = 0; i < nbounddist; ++i )
4698  {
4699  int k;
4700  SCIP_Real newefficacy;
4701  SCIP_Real QUAD(newrhs);
4702  SCIP_Real bestlb;
4703  SCIP_Real bestub;
4704  SCIP_Real oldsolval;
4705  int bestlbtype;
4706  int bestubtype;
4707 
4708  k = bounddistpos[i];
4709 
4710  SCIP_CALL( findBestLb(scip, vars[mksetinds[k]], sol, FALSE, allowlocal, &bestlb, &bestlbtype) );
4711 
4712  if( SCIPisInfinity(scip, -bestlb) )
4713  continue;
4714 
4715  SCIP_CALL( findBestUb(scip, vars[mksetinds[k]], sol, FALSE, allowlocal, &bestub, &bestubtype) );
4716 
4717  if( SCIPisInfinity(scip, bestub) )
4718  continue;
4719 
4720  /* switch the complementation of this variable */
4721 #ifndef NDEBUG
4722  {
4723  SCIP_Real QUAD(coef);
4724  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[k]);
4725  assert(SCIPisEQ(scip, tmpcoefs[k - intstart], varsign[k] * QUAD_TO_DBL(coef)));
4726  }
4727 #endif
4728 
4729  /* compute this: newrhs = mksetrhs + tmpcoefs[k - intstart] * (bestlb - bestub); */
4730  SCIPquadprecSumQD(newrhs, mksetrhs, tmpcoefs[k - intstart] * (bestlb - bestub));
4731  tmpcoefs[k - intstart] = -tmpcoefs[k - intstart];
4732 
4733  oldsolval = tmpvalues[k - intstart];
4734  tmpvalues[k - intstart] = varsign[k] == +1 ? bestub - SCIPgetSolVal(scip, sol, vars[mksetinds[k]]) : SCIPgetSolVal(scip, sol, vars[mksetinds[k]]) - bestlb;
4735 
4736  /* compute new violation */
4737  newefficacy = computeMIREfficacy(scip, tmpcoefs, tmpvalues, QUAD_TO_DBL(newrhs), contactivity, contsqrnorm, bestdelta, ntmpcoefs, minfrac, maxfrac);
4738 
4739  /* check if violaton was increased */
4740  if( newefficacy > bestefficacy )
4741  {
4742  /* keep change of complementation */
4743  bestefficacy = newefficacy;
4744  QUAD_ASSIGN_Q(mksetrhs, newrhs);
4745 
4746  if( varsign[k] == +1 )
4747  {
4748  /* switch to upper bound */
4749  assert(bestubtype < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
4750  boundtype[k] = bestubtype;
4751  varsign[k] = -1;
4752  }
4753  else
4754  {
4755  /* switch to lower bound */
4756  assert(bestlbtype < 0); /* cannot switch to a variable bound (would lead to further coef updates) */
4757  boundtype[k] = bestlbtype;
4758  varsign[k] = +1;
4759  }
4760 
4761  localbdsused = localbdsused || (boundtype[k] == -2);
4762  }
4763  else
4764  {
4765  /* undo the change of the complementation */
4766  tmpcoefs[k - intstart] = -tmpcoefs[k - intstart];
4767  tmpvalues[k - intstart] = oldsolval;
4768  }
4769  }
4770 
4771  if( bestefficacy > 0.0 )
4772  {
4773  SCIP_Real mirefficacy;
4774  SCIP_Real QUAD(downrhs);
4775  SCIP_Real QUAD(f0);
4776  SCIP_Real scale;
4777 
4778  scale = 1.0 / bestdelta;
4779  SCIPquadprecProdQD(mksetrhs, mksetrhs, scale);
4780 
4781  SCIPquadprecEpsFloorQ(downrhs, mksetrhs, SCIPepsilon(scip)); /*lint !e666*/
4782  SCIPquadprecSumQQ(f0, mksetrhs, -downrhs);
4783 
4784  /* renormaliize f0 value */
4785  SCIPquadprecSumDD(f0, QUAD_HI(f0), QUAD_LO(f0));
4786 
4787  for( i = 0; i < mksetnnz; ++i )
4788  {
4789  SCIP_Real QUAD(coef);
4790 
4791  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4792  SCIPquadprecProdQD(coef, coef, scale);
4793  QUAD_ARRAY_STORE(mksetcoefs, mksetinds[i], coef);
4794  }
4795  SCIPdebugMessage("applied best scale (=%.13g):\n", scale);
4796  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4797 
4798  QUAD_ASSIGN_Q(mksetrhs, downrhs);
4799 
4800  SCIP_CALL( cutsRoundMIR(scip, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, varsign, boundtype, QUAD(f0)) );
4801 
4802  SCIPdebugMessage("rounded MIR cut:\n");
4803  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4804 
4805  /* substitute aggregated slack variables:
4806  *
4807  * The coefficient of the slack variable s_r is equal to the row's weight times the slack's sign, because the slack
4808  * variable only appears in its own row:
4809  * a'_r = scale * weight[r] * slacksign[r].
4810  *
4811  * Depending on the slacks type (integral or continuous), its coefficient in the cut calculates as follows:
4812  * integers : a^_r = a~_r = down(a'_r) , if f_r <= f0
4813  * a^_r = a~_r = down(a'_r) + (f_r - f0)/(1 - f0), if f_r > f0
4814  * continuous: a^_r = a~_r = 0 , if a'_r >= 0
4815  * a^_r = a~_r = a'_r/(1 - f0) , if a'_r < 0
4816  *
4817  * Substitute a^_r * s_r by adding a^_r times the slack's definition to the cut.
4818  */
4819  SCIP_CALL( cutsSubstituteMIR(scip, aggrrow->rowweights, aggrrow->slacksign, aggrrow->rowsinds,
4820  aggrrow->nrows, scale, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz, QUAD(f0)) );
4821 
4822  SCIPdebugMessage("substituted slacks in MIR cut:\n");
4823  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4824 
4825 #ifndef NDEBUG
4826  {
4827  SCIP_Real efficacy = -QUAD_TO_DBL(mksetrhs);
4828  for( i = 0; i < mksetnnz; ++i )
4829  {
4830  SCIP_Real QUAD(coef);
4831  QUAD_ARRAY_LOAD(coef, mksetcoefs, mksetinds[i]);
4832  efficacy += QUAD_TO_DBL(coef) * SCIPgetSolVal(scip, sol, vars[mksetinds[i]]);
4833  }
4834 
4835  if( !EPSZ(SCIPrelDiff(efficacy, bestefficacy), 1e-4) )
4836  {
4837  SCIPdebugMessage("efficacy of cmir cut is different than expected efficacy: %f != %f\n", efficacy, bestefficacy);
4838  }
4839  }
4840 #endif
4841 
4842  *cutislocal = *cutislocal || localbdsused;
4843 
4844  /* remove all nearly-zero coefficients from MIR row and relax the right hand side correspondingly in order to
4845  * prevent numerical rounding errors
4846  */
4847  if( postprocess )
4848  {
4849  SCIP_CALL( postprocessCutQuad(scip, *cutislocal, mksetinds, mksetcoefs, &mksetnnz, QUAD(&mksetrhs), success) );
4850  }
4851  else
4852  {
4853  *success = ! removeZerosQuad(scip, SCIPsumepsilon(scip), *cutislocal, mksetcoefs, QUAD(&mksetrhs), mksetinds, &mksetnnz);
4854  }
4855 
4856  SCIPdebugMessage("post-processed cut (success = %s):\n", *success ? "TRUE" : "FALSE");
4857  SCIPdebug(printCutQuad(scip, sol, mksetcoefs, QUAD(mksetrhs), mksetinds, mksetnnz, FALSE, FALSE));
4858 
4859  if( *success )
4860  {
4861  mirefficacy = calcEfficacyDenseStorageQuad(scip, sol, mksetcoefs, QUAD_TO_DBL(mksetrhs), mksetinds, mksetnnz);
4862 
4863  if( SCIPisEfficacious(scip, mirefficacy) && mirefficacy > *cutefficacy )
4864  {
4865  BMScopyMemoryArray(cutinds, mksetinds, mksetnnz);
4866  for( i = 0; i < mksetnnz; ++i )
4867  {
4868  SCIP_Real QUAD(coef);
4869  int j = cutinds[i];
4870 
4871  QUAD_ARRAY_LOAD(coef, mksetcoefs, j);
4872 
4873  cutcoefs[i] = QUAD_TO_DBL(coef);
4874  QUAD_ASSIGN(coef, 0.0);
4875  QUAD_ARRAY_STORE(mksetcoefs, j, coef);
4876  }
4877  *cutnnz = mksetnnz;
4878  *cutrhs = QUAD_TO_DBL(mksetrhs);
4879  *cutefficacy = mirefficacy;
4880  if( cutrank != NULL )
4881  *cutrank = aggrrow->rank + 1;
4882  *cutislocal = *cutislocal || localbdsused;
4883  }
4884  else
4885  {
4886  *success = FALSE;
4887  }
4888  }
4889  }
4890 
4891  TERMINATE:
4892  /* if we aborted early we need to clean the mksetcoefs */
4893  if( !(*success) )
4894  {
4895  SCIP_Real QUAD(tmp);
4896  QUAD_ASSIGN(tmp, 0.0);
4897 
4898  for( i = 0; i < mksetnnz; ++i )
4899  {
4900  QUAD_ARRAY_STORE(mksetcoefs, mksetinds[i], tmp);
4901  }
4902  }
4903 
4904  /* free temporary memory */
4905  SCIPfreeBufferArray(scip, &bounddistpos);
4906  SCIPfreeBufferArray(scip, &bounddist);
4907  SCIPfreeBufferArray(scip, &deltacands);
4908  SCIPfreeBufferArray(scip, &tmpvalues);
4909  SCIPfreeBufferArray(scip, &tmpcoefs);
4910  SCIPfreeBufferArray(scip, &mksetinds);
4911  SCIPfreeCleanBufferArray(scip, &mksetcoefs);
4912  SCIPfreeBufferArray(scip, &boundtype);
4913  SCIPfreeBufferArray(scip, &varsign);
4914 
4915  return SCIP_OKAY;
4916 }
4917 
4918 /* =========================================== flow cover =========================================== */
4919 
4920 #define NO_EXACT_KNAPSACK
4921 
4922 #ifndef NO_EXACT_KNAPSACK
4923 #define MAXDNOM 1000LL
4924 #define MINDELTA 1e-03
4925 #define MAXDELTA 1e-09
4926 #define MAXSCALE 1000.0
4927 #define MAXDYNPROGSPACE 1000000
4928 #endif
4929 
4930 #define MAXABSVBCOEF 1e+5 /**< maximal absolute coefficient in variable bounds used for snf relaxation */
4931 #define MAXBOUND 1e+10 /**< maximal value of normal bounds used for snf relaxation */
4932 
4933 /** structure that contains all data required to perform the sequence independent lifting
4934  */
4935 typedef
4936 struct LiftingData
4937 {
4938  SCIP_Real* M; /**< \f$ M_0 := 0.0 \f$ and \f$ M_i := M_i-1 + m_i \f$ */
4939  SCIP_Real* m; /**< non-increasing array of variable upper bound coefficients
4940  * for all variables in \f$ C^{++} \f$ and \f$ L^- \f$,
4941  * where \f$ C = C^+ \cup C^- \f$ is the flowcover and
4942  * \f$ C^{++} := \{ j \in C^+ \mid u_j > \lambda \} \f$
4943  * \f$ L^- := \{ j \in (N^- \setminus C^-) \mid u_j > \lambda \} \f$
4944  */
4945  int r; /**< size of array m */
4946  int t; /**< index of smallest value in m that comes from a variable in \f$ C^{++} \f$ */
4947  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$ */
4948  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$ */
4949  SCIP_Real lambda; /**< excess of the flowcover */
4950  SCIP_Real mp; /**< smallest variable bound coefficient of variable in \f$ C^{++} (min_{j \in C++} u_j) \f$ */
4951  SCIP_Real ml; /**< \f$ ml := min(\lambda, \sum_{j \in C^+ \setminus C^{++}} u_j) \f$ */
4952 } LIFTINGDATA;
4953 
4954 /** structure that contains all the data that defines the single-node-flow relaxation of an aggregation row */
4955 typedef
4956 struct SNF_Relaxation
4957 {
4958  int* transvarcoefs; /**< coefficients of all vars in relaxed set */
4959  SCIP_Real* transbinvarsolvals; /**< sol val of bin var in vub of all vars in relaxed set */
4960  SCIP_Real* transcontvarsolvals;/**< sol val of all real vars in relaxed set */
4961  SCIP_Real* transvarvubcoefs; /**< coefficient in vub of all vars in relaxed set */
4962  int ntransvars; /**< number of vars in relaxed set */
4963  SCIP_Real transrhs; /**< rhs in relaxed set */
4964  int* origbinvars; /**< associated original binary var for all vars in relaxed set */
4965  int* origcontvars; /**< associated original continuous var for all vars in relaxed set */
4966  SCIP_Real* aggrcoefsbin; /**< aggregation coefficient of the original binary var used to define the
4967  * continuous variable in the relaxed set */
4968  SCIP_Real* aggrcoefscont; /**< aggregation coefficient of the original continous var used to define the
4969  * continuous variable in the relaxed set */
4970  SCIP_Real* aggrconstants; /**< aggregation constant used to define the continuous variable in the relaxed set */
4971 } SNF_RELAXATION;
4972 
4973 /** get solution value and index of variable lower bound (with binary variable) which is closest to the current LP
4974  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
4975  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
4976  * given variable
4977  */
4978 static
4980  SCIP* scip, /**< SCIP data structure */
4981  SCIP_VAR* var, /**< given active problem variable */
4982  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
4983  SCIP_Real* rowcoefs, /**< (dense) array of coefficients of row */
4984  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
4985  * was not used (0) or was not used but is contained in the row (-1)
4986  */
4987  SCIP_Real bestsub, /**< closest simple upper bound of given variable */
4988  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
4989  SCIP_Real* closestvlb, /**< pointer to store the LP sol value of the closest variable lower bound */
4990  int* closestvlbidx /**< pointer to store the index of the closest vlb; -1 if no vlb was found */
4991  )
4992 {
4993  int nvlbs;
4994  int nbinvars;
4995 
4996  assert(scip != NULL);
4997  assert(var != NULL);
4998  assert(bestsub == SCIPvarGetUbGlobal(var) || bestsub == SCIPvarGetUbLocal(var)); /*lint !e777*/
4999  assert(!SCIPisInfinity(scip, bestsub));
5000  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5001  assert(rowcoefs != NULL);
5002  assert(binvarused != NULL);
5003  assert(closestvlb != NULL);
5004  assert(closestvlbidx != NULL);
5005 
5006  nvlbs = SCIPvarGetNVlbs(var);
5007  nbinvars = SCIPgetNBinVars(scip);
5008 
5009  *closestvlbidx = -1;
5010  *closestvlb = -SCIPinfinity(scip);
5011  if( nvlbs > 0 )
5012  {
5013  SCIP_VAR** vlbvars;
5014  SCIP_Real* vlbcoefs;
5015  SCIP_Real* vlbconsts;
5016  int i;
5017 
5018  vlbvars = SCIPvarGetVlbVars(var);
5019  vlbcoefs = SCIPvarGetVlbCoefs(var);
5020  vlbconsts = SCIPvarGetVlbConstants(var);
5021 
5022  for( i = 0; i < nvlbs; i++ )
5023  {
5024  SCIP_Real rowcoefbinvar;
5025  SCIP_Real val1;
5026  SCIP_Real val2;
5027  SCIP_Real vlbsol;
5028  SCIP_Real rowcoefsign;
5029  int probidxbinvar;
5030 
5031  if( bestsub > vlbconsts[i] )
5032  continue;
5033 
5034  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
5035  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
5036  */
5037  if( REALABS(vlbcoefs[i]) > MAXABSVBCOEF )
5038  continue;
5039 
5040  /* use only variable lower bounds l~_i * x_i + d_i with x_i binary which are active */
5041  probidxbinvar = SCIPvarGetProbindex(vlbvars[i]);
5042 
5043  /* if the variable is not active the problem index is -1, so we cast to unsigned int before the comparison which
5044  * ensures that the problem index is between 0 and nbinvars - 1
5045  */
5046  if( (unsigned int)probidxbinvar >= (unsigned int)nbinvars )
5047  continue;
5048 
5049  assert(SCIPvarIsBinary(vlbvars[i]));
5050 
5051  /* check if current variable lower bound l~_i * x_i + d_i imposed on y_j meets the following criteria:
5052  * (let a_j = coefficient of y_j in current row,
5053  * u_j = closest simple upper bound imposed on y_j,
5054  * c_i = coefficient of x_i in current row)
5055  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k yet
5056  * if a_j > 0:
5057  * 1. u_j <= d_i
5058  * 2. a_j ( u_j - d_i ) + c_i <= 0
5059  * 3. a_j l~_i + c_i <= 0
5060  * if a_j < 0:
5061  * 1. u_j <= d_i
5062  * 2. a_j ( u_j - d_i ) + c_i >= 0
5063  * 3. a_j l~_i + c_i >= 0
5064  */
5065 
5066  /* has already been used in the SNF relaxation */
5067  if( binvarused[probidxbinvar] == 1 )
5068  continue;
5069 
5070  /* get the row coefficient */
5071  {
5072  SCIP_Real QUAD(tmp);
5073  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidxbinvar);
5074  rowcoefbinvar = QUAD_TO_DBL(tmp);
5075  }
5076  rowcoefsign = COPYSIGN(1.0, rowcoef);
5077 
5078  val2 = rowcoefsign * ((rowcoef * vlbcoefs[i]) + rowcoefbinvar);
5079 
5080  /* variable lower bound does not meet criteria */
5081  if( val2 > 0.0 || SCIPisInfinity(scip, -val2) )
5082  continue;
5083 
5084  val1 = rowcoefsign * ((rowcoef * (bestsub - vlbconsts[i])) + rowcoefbinvar);
5085 
5086  /* variable lower bound does not meet criteria */
5087  if( val1 > 0.0 )
5088  continue;
5089 
5090  vlbsol = vlbcoefs[i] * SCIPgetSolVal(scip, sol, vlbvars[i]) + vlbconsts[i];
5091  if( vlbsol > *closestvlb )
5092  {
5093  *closestvlb = vlbsol;
5094  *closestvlbidx = i;
5095  }
5096  assert(*closestvlbidx >= 0);
5097  }
5098  }
5099 
5100  return SCIP_OKAY;
5101 }
5102 
5103 /** get LP solution value and index of variable upper bound (with binary variable) which is closest to the current LP
5104  * solution value of a given variable; candidates have to meet certain criteria in order to ensure the nonnegativity
5105  * of the variable upper bound imposed on the real variable in the 0-1 single node flow relaxation associated with the
5106  * given variable
5107  */
5108 static
5110  SCIP* scip, /**< SCIP data structure */
5111  SCIP_VAR* var, /**< given active problem variable */
5112  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
5113  SCIP_Real* rowcoefs, /**< (dense) array of coefficients of row */
5114  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
5115  * was not used (0) or was not used but is contained in the row (-1)
5116  */
5117  SCIP_Real bestslb, /**< closest simple lower bound of given variable */
5118  SCIP_Real rowcoef, /**< coefficient of given variable in current row */
5119  SCIP_Real* closestvub, /**< pointer to store the LP sol value of the closest variable upper bound */
5120  int* closestvubidx /**< pointer to store the index of the closest vub; -1 if no vub was found */
5121  )
5122 {
5123  int nvubs;
5124  int nbinvars;
5125 
5126  assert(scip != NULL);
5127  assert(var != NULL);
5128  assert(bestslb == SCIPvarGetLbGlobal(var) || bestslb == SCIPvarGetLbLocal(var)); /*lint !e777*/
5129  assert(!SCIPisInfinity(scip, - bestslb));
5130  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5131  assert(rowcoefs != NULL);
5132  assert(binvarused != NULL);
5133  assert(closestvub != NULL);
5134  assert(closestvubidx != NULL);
5135 
5136  nvubs = SCIPvarGetNVubs(var);
5137  nbinvars = SCIPgetNBinVars(scip);
5138 
5139  *closestvubidx = -1;
5140  *closestvub = SCIPinfinity(scip);
5141  if( nvubs > 0 )
5142  {
5143  SCIP_VAR** vubvars;
5144  SCIP_Real* vubcoefs;
5145  SCIP_Real* vubconsts;
5146  int i;
5147 
5148  vubvars = SCIPvarGetVubVars(var);
5149  vubcoefs = SCIPvarGetVubCoefs(var);
5150  vubconsts = SCIPvarGetVubConstants(var);
5151 
5152  for( i = 0; i < nvubs; i++ )
5153  {
5154  SCIP_Real rowcoefbinvar;
5155  SCIP_Real val1;
5156  SCIP_Real val2;
5157  SCIP_Real vubsol;
5158  SCIP_Real rowcoefsign;
5159  int probidxbinvar;
5160 
5161  if( bestslb < vubconsts[i] )
5162  continue;
5163 
5164  /* for numerical reasons, ignore variable bounds with large absolute coefficient and
5165  * those which lead to an infinite variable bound coefficient (val2) in snf relaxation
5166  */
5167  if( REALABS(vubcoefs[i]) > MAXABSVBCOEF )
5168  continue;
5169 
5170  /* use only variable upper bound u~_i * x_i + d_i with x_i binary and which are active */
5171  probidxbinvar = SCIPvarGetProbindex(vubvars[i]);
5172 
5173  /* if the variable is not active the problem index is -1, so we cast to unsigned int before the comparison which
5174  * ensures that the problem index is between 0 and nbinvars - 1
5175  */
5176  if( (unsigned int)probidxbinvar >= (unsigned int)nbinvars )
5177  continue;
5178 
5179  assert(SCIPvarIsBinary(vubvars[i]));
5180 
5181  /* checks if current variable upper bound u~_i * x_i + d_i meets the following criteria
5182  * (let a_j = coefficient of y_j in current row,
5183  * l_j = closest simple lower bound imposed on y_j,
5184  * c_i = coefficient of x_i in current row)
5185  * 0. no other non-binary variable y_k has used a variable bound with x_i to get transformed variable y'_k
5186  * if a > 0:
5187  * 1. l_j >= d_i
5188  * 2. a_j ( l_i - d_i ) + c_i >= 0
5189  * 3. a_j u~_i + c_i >= 0
5190  * if a < 0:
5191  * 1. l_j >= d_i
5192  * 2. a_j ( l_j - d_i ) + c_i <= 0
5193  * 3. a_j u~_i + c_i <= 0
5194  */
5195 
5196  /* has already been used in the SNF relaxation */
5197  if( binvarused[probidxbinvar] == 1 )
5198  continue;
5199 
5200  /* get the row coefficient */
5201  {
5202  SCIP_Real QUAD(tmp);
5203  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidxbinvar);
5204  rowcoefbinvar = QUAD_TO_DBL(tmp);
5205  }
5206  rowcoefsign = COPYSIGN(1.0, rowcoef);
5207 
5208  val2 = rowcoefsign * ((rowcoef * vubcoefs[i]) + rowcoefbinvar);
5209 
5210  /* variable upper bound does not meet criteria */
5211  if( val2 < 0.0 || SCIPisInfinity(scip, val2) )
5212  continue;
5213 
5214  val1 = rowcoefsign * ((rowcoef * (bestslb - vubconsts[i])) + rowcoefbinvar);
5215 
5216  /* variable upper bound does not meet criteria */
5217  if( val1 < 0.0 )
5218  continue;
5219 
5220  vubsol = vubcoefs[i] * SCIPgetSolVal(scip, sol, vubvars[i]) + vubconsts[i];
5221  if( vubsol < *closestvub )
5222  {
5223  *closestvub = vubsol;
5224  *closestvubidx = i;
5225  }
5226  assert(*closestvubidx >= 0);
5227  }
5228  }
5229 
5230  return SCIP_OKAY;
5231 }
5232 
5233 /** determines the bounds to use for constructing the single-node-flow relaxation of a variable in
5234  * the given row.
5235  */
5236 static
5238  SCIP* scip, /**< SCIP data structure */
5239  SCIP_SOL* sol, /**< solution to use for variable bound; NULL for LP solution */
5240  SCIP_VAR** vars, /**< array of problem variables */
5241  SCIP_Real* rowcoefs, /**< (dense) array of variable coefficients in the row */
5242  int* rowinds, /**< array with positions of non-zero values in the rowcoefs array */
5243  int varposinrow, /**< position of variable in the rowinds array for which the bounds should be determined */
5244  int8_t* binvarused, /**< array that stores if a binary variable was already used (+1)
5245  * was not used (0) or was not used but is contained in the row (-1)
5246  */
5247  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
5248  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
5249  SCIP_Real* bestlb, /**< pointer to store best lower bound for transformation */
5250  SCIP_Real* bestub, /**< pointer to store best upper bound for transformation */
5251  SCIP_Real* bestslb, /**< pointer to store best simple lower bound for transformation */
5252  SCIP_Real* bestsub, /**< pointer to store best simple upper bound for transformation */
5253  int* bestlbtype, /**< pointer to store type of best lower bound */
5254  int* bestubtype, /**< pointer to store type of best upper bound */
5255  int* bestslbtype, /**< pointer to store type of best simple lower bound */
5256  int* bestsubtype, /**< pointer to store type of best simple upper bound */
5257  SCIP_BOUNDTYPE* selectedbounds, /**< pointer to store the preferred bound for the transformation */
5258  SCIP_Bool* freevariable /**< pointer to store if variable is a free variable */
5259  )
5260 {
5261  SCIP_VAR* var;
5262 
5263  SCIP_Real rowcoef;
5264  SCIP_Real solval;
5265 
5266  int probidx;
5267 
5268  bestlb[varposinrow] = -SCIPinfinity(scip);
5269  bestub[varposinrow] = SCIPinfinity(scip);
5270  bestlbtype[varposinrow] = -3;
5271  bestubtype[varposinrow] = -3;
5272 
5273  probidx = rowinds[varposinrow];
5274  var = vars[probidx];
5275  {
5276  SCIP_Real QUAD(tmp);
5277  QUAD_ARRAY_LOAD(tmp, rowcoefs, probidx);
5278  rowcoef = QUAD_TO_DBL(tmp);
5279  }
5280 
5281  assert(!EPSZ(rowcoef, QUAD_EPSILON));
5282 
5283  /* get closest simple lower bound and closest simple upper bound */
5284  SCIP_CALL( findBestLb(scip, var, sol, FALSE, allowlocal, &bestslb[varposinrow], &bestslbtype[varposinrow]) );
5285  SCIP_CALL( findBestUb(scip, var, sol, FALSE, allowlocal, &bestsub[varposinrow], &bestsubtype[varposinrow]) );
5286 
5287  /* do not use too large bounds */
5288  if( bestslb[varposinrow] <= -MAXBOUND )
5289  bestslb[varposinrow] = -SCIPinfinity(scip);
5290 
5291  if( bestsub[varposinrow] >= MAXBOUND )
5292  bestsub[varposinrow] = SCIPinfinity(scip);
5293 
5294  solval = SCIPgetSolVal(scip, sol, var);
5295 
5296  SCIPdebugMsg(scip, " %d: %g <%s, idx=%d, lp=%g, [%g(%d),%g(%d)]>:\n", varposinrow, rowcoef, SCIPvarGetName(var), probidx,
5297  solval, bestslb[varposinrow], bestslbtype[varposinrow], bestsub[varposinrow], bestsubtype[varposinrow]);
5298 
5299  /* mixed integer set cannot be relaxed to 0-1 single node flow set because both simple bounds are -infinity
5300  * and infinity, respectively
5301  */
5302  if( SCIPisInfinity(scip, -bestslb[varposinrow]) && SCIPisInfinity(scip, bestsub[varposinrow]) )
5303  {
5304  *freevariable = TRUE;
5305  return SCIP_OKAY;
5306  }
5307 
5308  /* get closest lower bound that can be used to define the real variable y'_j in the 0-1 single node flow
5309  * relaxation
5310  */
5311  if( !SCIPisInfinity(scip, bestsub[varposinrow]) )
5312  {
5313  bestlb[varposinrow] = bestslb[varposinrow];
5314  bestlbtype[varposinrow] = bestslbtype[varposinrow];
5315 
5317  {
5318  SCIP_Real bestvlb;
5319  int bestvlbidx;
5320 
5321  SCIP_CALL( getClosestVlb(scip, var, sol, rowcoefs, binvarused, bestsub[varposinrow], rowcoef, &bestvlb, &bestvlbidx) );
5322  if( SCIPisGT(scip, bestvlb, bestlb[varposinrow]) )
5323  {
5324  bestlb[varposinrow] = bestvlb;
5325  bestlbtype[varposinrow] = bestvlbidx;
5326  }
5327  }
5328  }
5329 
5330  /* get closest upper bound that can be used to define the real variable y'_j in the 0-1 single node flow
5331  * relaxation
5332  */
5333  if( !SCIPisInfinity(scip, -bestslb[varposinrow]) )
5334  {
5335  bestub[varposinrow] = bestsub[varposinrow];
5336  bestubtype[varposinrow] = bestsubtype[varposinrow];
5337 
5339  {
5340  SCIP_Real bestvub;
5341  int bestvubidx;
5342 
5343  SCIP_CALL( getClosestVub(scip, var, sol, rowcoefs, binvarused, bestslb[varposinrow], rowcoef, &bestvub, &bestvubidx) );
5344  if( SCIPisLT(scip, bestvub, bestub[varposinrow]) )
5345  {
5346  bestub[varposinrow] = bestvub;
5347  bestubtype[varposinrow] = bestvubidx;
5348  }
5349  }
5350  }
5351  SCIPdebugMsg(scip, " bestlb=%g(%d), bestub=%g(%d)\n", bestlb[varposinrow], bestlbtype[varposinrow], bestub[varposinrow], bestubtype[varposinrow]);
5352 
5353  /* mixed integer set cannot be relaxed to 0-1 single node flow set because there are no suitable bounds
5354  * to define the transformed variable y'_j
5355  */
5356  if( SCIPisInfinity(scip, -bestlb[varposinrow]) && SCIPisInfinity(scip, bestub[varposinrow]) )
5357  {
5358  *freevariable = TRUE;
5359  return SCIP_OKAY;
5360  }
5361 
5362  *freevariable = FALSE;
5363 
5364  /* select best upper bound if it is closer to the LP value of y_j and best lower bound otherwise and use this bound
5365  * to define the real variable y'_j with 0 <= y'_j <= u'_j x_j in the 0-1 single node flow relaxation;
5366  * prefer variable bounds
5367  */
5368  if( SCIPisEQ(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]) && bestlbtype[varposinrow] >= 0 )
5369  {
5370  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_LOWER;
5371  }
5372  else if( SCIPisEQ(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow])
5373  && bestubtype[varposinrow] >= 0 )
5374  {
5375  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_UPPER;
5376  }
5377  else if( SCIPisLE(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]) )
5378  {
5379  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_LOWER;
5380  }
5381  else
5382  {
5383  assert(SCIPisGT(scip, solval, (1.0 - boundswitch) * bestlb[varposinrow] + boundswitch * bestub[varposinrow]));
5384  selectedbounds[varposinrow] = SCIP_BOUNDTYPE_UPPER;
5385  }
5386 
5387  if( selectedbounds[varposinrow] == SCIP_BOUNDTYPE_LOWER && bestlbtype[varposinrow] >= 0 )
5388  {
5389  int vlbvarprobidx;
5390  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(var);
5391 
5392  /* mark binary variable of vlb so that it is not used for other continuous variables
5393  * by setting it's position in the aggrrow to a negative value
5394  */
5395  vlbvarprobidx = SCIPvarGetProbindex(vlbvars[bestlbtype[varposinrow]]);
5396  binvarused[vlbvarprobidx] = 1;
5397  }
5398  else if ( selectedbounds[varposinrow] == SCIP_BOUNDTYPE_UPPER && bestubtype[varposinrow] >= 0 )
5399  {
5400  int vubvarprobidx;
5401  SCIP_VAR** vubvars = SCIPvarGetVubVars(var);
5402 
5403  /* mark binary variable of vub so that it is not used for other continuous variables
5404  * by setting it's position in the aggrrow to a negative value
5405  */
5406  vubvarprobidx = SCIPvarGetProbindex(vubvars[bestubtype[varposinrow]]);
5407  binvarused[vubvarprobidx] = 1;
5408  }
5409 
5410  return SCIP_OKAY;
5411 }
5412 
5413 /** construct a 0-1 single node flow relaxation (with some additional simple constraints) of a mixed integer set
5414  * corresponding to the given aggrrow a * x <= rhs
5415  */
5416 static
5418  SCIP* scip, /**< SCIP data structure */
5419  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5420  SCIP_Real boundswitch, /**< fraction of domain up to which lower bound is used in transformation */
5421  SCIP_Bool allowlocal, /**< should local information allowed to be used, resulting in a local cut? */
5422  SCIP_Real* rowcoefs, /**< array of coefficients of row */
5423  QUAD(SCIP_Real rowrhs), /**< pointer to right hand side of row */
5424  int* rowinds, /**< array of variables problem indices for non-zero coefficients in row */
5425  int nnz, /**< number of non-zeros in row */
5426  SNF_RELAXATION* snf, /**< stores the sign of the transformed variable in summation */
5427  SCIP_Bool* success, /**< stores whether the transformation was valid */
5428  SCIP_Bool* localbdsused /**< pointer to store whether local bounds were used in transformation */
5429  )
5430 {
5431  SCIP_VAR** vars;
5432  int i;
5433  int nnonbinvarsrow;
5434  int8_t* binvarused;
5435  int nbinvars;
5436  SCIP_Real QUAD(transrhs);
5437 
5438  /* arrays to store the selected bound for each non-binary variable in the row */
5439  SCIP_Real* bestlb;
5440  SCIP_Real* bestub;
5441  SCIP_Real* bestslb;
5442  SCIP_Real* bestsub;
5443  int* bestlbtype;
5444  int* bestubtype;
5445  int* bestslbtype;
5446  int* bestsubtype;
5447  SCIP_BOUNDTYPE* selectedbounds;
5448 
5449  *success = FALSE;
5450 
5451  SCIPdebugMsg(scip, "--------------------- construction of SNF relaxation ------------------------------------\n");
5452 
5453  nbinvars = SCIPgetNBinVars(scip);
5454  vars = SCIPgetVars(scip);
5455 
5456  SCIP_CALL( SCIPallocBufferArray(scip, &bestlb, nnz) );
5457  SCIP_CALL( SCIPallocBufferArray(scip, &bestub, nnz) );
5458  SCIP_CALL( SCIPallocBufferArray(scip, &bestslb, nnz) );
5459  SCIP_CALL( SCIPallocBufferArray(scip, &bestsub, nnz) );
5460  SCIP_CALL( SCIPallocBufferArray(scip, &bestlbtype, nnz) );
5461  SCIP_CALL( SCIPallocBufferArray(scip, &bestubtype, nnz) );
5462  SCIP_CALL( SCIPallocBufferArray(scip, &bestslbtype, nnz) );
5463  SCIP_CALL( SCIPallocBufferArray(scip, &bestsubtype, nnz) );
5464  SCIP_CALL( SCIPallocBufferArray(scip, &selectedbounds, nnz) );
5465 
5466  /* sort descending to have continuous variables first */
5467  SCIPsortDownInt(rowinds, nnz);
5468 
5469  /* array to store whether a binary variable is in the row (-1) or has been used (1) due to variable bound usage */
5470  SCIP_CALL( SCIPallocCleanBufferArray(scip, &binvarused, nbinvars) );
5471 
5472  for( i = nnz - 1; i >= 0 && rowinds[i] < nbinvars; --i )
5473  {
5474  int j = rowinds[i];
5475  binvarused[j] = -1;
5476  }
5477 
5478  nnonbinvarsrow = i + 1;
5479  /* determine the bounds to use for transforming the non-binary variables */
5480  for( i = 0; i < nnonbinvarsrow; ++i )
5481  {
5482  SCIP_Bool freevariable;
5483 
5484  assert(rowinds[i] >= nbinvars);
5485 
5486  SCIP_CALL( determineBoundForSNF(scip, sol, vars, rowcoefs, rowinds, i, binvarused, allowlocal, boundswitch,
5487  bestlb, bestub, bestslb, bestsub, bestlbtype, bestubtype, bestslbtype, bestsubtype, selectedbounds, &freevariable) );
5488 
5489  if( freevariable )
5490  {
5491  int j;
5492 
5493  /* clear binvarused at indices of binary variables of row */
5494  for( j = nnz - 1; j >= nnonbinvarsrow; --j )
5495  binvarused[rowinds[j]] = 0;
5496 
5497  /* clear binvarused at indices of selected variable bounds */
5498  for( j = 0; j < i; ++j )
5499  {
5500  if( selectedbounds[j] == SCIP_BOUNDTYPE_LOWER && bestlbtype[j] >= 0 )
5501  {
5502  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(vars[rowinds[j]]);
5503  binvarused[SCIPvarGetProbindex(vlbvars[bestlbtype[j]])] = 0;
5504  }
5505  else if( selectedbounds[j] == SCIP_BOUNDTYPE_UPPER && bestubtype[j] >= 0 )
5506  {
5507  SCIP_VAR** vubvars = SCIPvarGetVubVars(vars[rowinds[j]]);
5508  binvarused[SCIPvarGetProbindex(vubvars[bestubtype[j]])] = 0;
5509  }
5510  }
5511 
5512  /* terminate */
5513  goto TERMINATE;
5514  }
5515  }
5516 
5517  *localbdsused = FALSE;
5518  QUAD_ASSIGN_Q(transrhs, rowrhs);
5519  snf->ntransvars = 0;
5520 
5521  /* transform non-binary variables */
5522  for( i = 0; i < nnonbinvarsrow; ++i )
5523  {
5524  SCIP_VAR* var;
5525  SCIP_Real QUAD(rowcoef);
5526  SCIP_Real solval;
5527  int probidx;
5528 
5529  probidx = rowinds[i];
5530  var = vars[probidx];
5531  QUAD_ARRAY_LOAD(rowcoef, rowcoefs, probidx);
5532  solval = SCIPgetSolVal(scip, sol, var);
5533 
5534  assert(probidx >= nbinvars);
5535 
5536  if( selectedbounds[i] == SCIP_BOUNDTYPE_LOWER )
5537  {
5538  /* use bestlb to define y'_j */
5539 
5540  assert(!SCIPisInfinity(scip, bestsub[i]));
5541  assert(!SCIPisInfinity(scip, - bestlb[i]));
5542  assert(bestsubtype[i] == -1 || bestsubtype[i] == -2);
5543  assert(bestlbtype[i] > -3 && bestlbtype[i] < SCIPvarGetNVlbs(var));
5544 
5545  /* store for y_j that bestlb is the bound used to define y'_j and that y'_j is the associated real variable
5546  * in the relaxed set
5547  */
5548  snf->origcontvars[snf->ntransvars] = probidx;
5549 
5550  if( bestlbtype[i] < 0 )
5551  {
5552  SCIP_Real QUAD(val);
5553  SCIP_Real QUAD(contsolval);
5554  SCIP_Real QUAD(rowcoeftimesbestsub);
5555 
5556  /* use simple lower bound in bestlb = l_j <= y_j <= u_j = bestsub to define
5557  * 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
5558  * 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,
5559  * put j into the set
5560  * N2 if a_j > 0
5561  * N1 if a_j < 0
5562  * and update the right hand side of the constraint in the relaxation
5563  * rhs = rhs - a_j u_j
5564  */
5565  SCIPquadprecSumDD(val, bestsub[i], -bestlb[i]);
5566  SCIPquadprecProdQQ(val, val, rowcoef);
5567  SCIPquadprecSumDD(contsolval, solval, -bestsub[i]);
5568  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5569 
5570  if( bestlbtype[i] == -2 || bestsubtype[i] == -2 )
5571  *localbdsused = TRUE;
5572 
5573  SCIPquadprecProdQD(rowcoeftimesbestsub, rowcoef, bestsub[i]);
5574 
5575  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5576  snf->origbinvars[snf->ntransvars] = -1;
5577  snf->aggrcoefsbin[snf->ntransvars] = 0.0;
5578 
5579  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5580  {
5581  snf->transvarcoefs[snf->ntransvars] = - 1;
5582  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5583  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5584  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5585 
5586  /* aggregation information for y'_j */
5587  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesbestsub);
5588  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5589  }
5590  else
5591  {
5592  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5593  snf->transvarcoefs[snf->ntransvars] = 1;
5594  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5595  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5596  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5597 
5598  /* aggregation information for y'_j */
5599  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesbestsub);
5600  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5601  }
5602  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesbestsub);
5603 
5604  SCIPdebugMsg(scip, " --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
5605  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5606  snf->ntransvars, QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesbestsub), QUAD_TO_DBL(rowcoef), bestsub, QUAD_TO_DBL(transrhs));
5607  }
5608  else
5609  {
5610  SCIP_Real QUAD(rowcoefbinary);
5611  SCIP_Real varsolvalbinary;
5612  SCIP_Real QUAD(val);
5613  SCIP_Real QUAD(contsolval);
5614  SCIP_Real QUAD(rowcoeftimesvlbconst);
5615  int vlbvarprobidx;
5616 
5617  SCIP_VAR** vlbvars = SCIPvarGetVlbVars(var);
5618  SCIP_Real* vlbconsts = SCIPvarGetVlbConstants(var);
5619  SCIP_Real* vlbcoefs = SCIPvarGetVlbCoefs(var);
5620 
5621  /* use variable lower bound in bestlb = l~_j x_j + d_j <= y_j <= u_j = bestsub to define
5622  * 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
5623  * 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,
5624  * where c_j is the coefficient of x_j in the row, put j into the set
5625  * N2 if a_j > 0
5626  * N1 if a_j < 0
5627  * and update the right hand side of the constraint in the relaxation
5628  * rhs = rhs - a_j d_j
5629  */
5630 
5631  vlbvarprobidx = SCIPvarGetProbindex(vlbvars[bestlbtype[i]]);
5632  assert(binvarused[vlbvarprobidx] == 1);
5633  assert(vlbvarprobidx < nbinvars);
5634 
5635  QUAD_ARRAY_LOAD(rowcoefbinary, rowcoefs, vlbvarprobidx);
5636  varsolvalbinary = SCIPgetSolVal(scip, sol, vlbvars[bestlbtype[i]]);
5637 
5638  SCIPquadprecProdQD(val, rowcoef, vlbcoefs[bestlbtype[i]]);
5639  SCIPquadprecSumQQ(val, val, rowcoefbinary);
5640  {
5641  SCIP_Real QUAD(tmp);
5642 
5643  SCIPquadprecProdQD(tmp, rowcoefbinary, varsolvalbinary);
5644  SCIPquadprecSumDD(contsolval, solval, - vlbconsts[bestlbtype[i]]);
5645  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5646  SCIPquadprecSumQQ(contsolval, contsolval, tmp);
5647  }
5648 
5649  SCIPquadprecProdQD(rowcoeftimesvlbconst, rowcoef, vlbconsts[bestlbtype[i]]);
5650 
5651  /* clear the binvarpos array, since the variable has been processed */
5652  binvarused[vlbvarprobidx] = 0;
5653 
5654  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5655  snf->origbinvars[snf->ntransvars] = vlbvarprobidx;
5656 
5657  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5658  {
5659  snf->transvarcoefs[snf->ntransvars] = - 1;
5660  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5661  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5662  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5663 
5664  /* aggregation information for y'_j */
5665  snf->aggrcoefsbin[snf->ntransvars] = - QUAD_TO_DBL(rowcoefbinary);
5666  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5667  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesvlbconst);
5668  }
5669  else
5670  {
5671  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5672  snf->transvarcoefs[snf->ntransvars] = 1;
5673  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5674  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5675  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5676 
5677  /* aggregation information for y'_j */
5678  snf->aggrcoefsbin[snf->ntransvars] = QUAD_TO_DBL(rowcoefbinary);
5679  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5680  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesvlbconst);
5681  }
5682  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesvlbconst);
5683 
5684  SCIPdebugMsg(scip, " --> bestlb used for trans: ... %s y'_%d + ..., y'_%d <= %g x_%d (=%s), rhs=%g-(%g*%g)=%g\n",
5685  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5686  snf->ntransvars, SCIPvarGetName(vlbvars[bestlbtype[i]]), QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesvlbconst), QUAD_TO_DBL(rowcoef),
5687  vlbconsts[bestlbtype[i]], snf->transrhs );
5688  }
5689  }
5690  else
5691  {
5692  /* use bestub to define y'_j */
5693 
5694  assert(!SCIPisInfinity(scip, bestub[i]));
5695  assert(!SCIPisInfinity(scip, - bestslb[i]));
5696  assert(bestslbtype[i] == -1 || bestslbtype[i] == -2);
5697  assert(bestubtype[i] > -3 && bestubtype[i] < SCIPvarGetNVubs(var));
5698 
5699  /* store for y_j that y'_j is the associated real variable
5700  * in the relaxed set
5701  */
5702  snf->origcontvars[snf->ntransvars] = probidx;
5703 
5704  if( bestubtype[i] < 0 )
5705  {
5706  SCIP_Real QUAD(val);
5707  SCIP_Real QUAD(contsolval);
5708  SCIP_Real QUAD(rowcoeftimesbestslb);
5709 
5710  /* use simple upper bound in bestslb = l_j <= y_j <= u_j = bestub to define
5711  * 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
5712  * 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,
5713  * put j into the set
5714  * N1 if a_j > 0
5715  * N2 if a_j < 0
5716  * and update the right hand side of the constraint in the relaxation
5717  * rhs = rhs - a_j l_j
5718  */
5719  SCIPquadprecSumDD(val, bestub[i], - bestslb[i]);
5720  SCIPquadprecProdQQ(val, val, rowcoef);
5721  SCIPquadprecSumDD(contsolval, solval, - bestslb[i]);
5722  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5723 
5724  if( bestubtype[i] == -2 || bestslbtype[i] == -2 )
5725  *localbdsused = TRUE;
5726 
5727  SCIPquadprecProdQD(rowcoeftimesbestslb, rowcoef, bestslb[i]);
5728 
5729  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5730  snf->origbinvars[snf->ntransvars] = -1;
5731  snf->aggrcoefsbin[snf->ntransvars] = 0.0;
5732 
5733  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5734  {
5735  snf->transvarcoefs[snf->ntransvars] = 1;
5736  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5737  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5738  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5739 
5740  /* aggregation information for y'_j */
5741  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5742  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesbestslb);
5743  }
5744  else
5745  {
5746  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5747  snf->transvarcoefs[snf->ntransvars] = - 1;
5748  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5749  snf->transbinvarsolvals[snf->ntransvars] = 1.0;
5750  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5751 
5752  /* aggregation information for y'_j */
5753  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5754  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesbestslb);
5755  }
5756  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesbestslb);
5757 
5758  SCIPdebugMsg(scip, " --> bestub used for trans: ... %s y'_%d + ..., Y'_%d <= %g x_%d (=1), rhs=%g-(%g*%g)=%g\n",
5759  snf->transvarcoefs[snf->ntransvars] == 1 ? "+" : "-", snf->ntransvars, snf->ntransvars, snf->transvarvubcoefs[snf->ntransvars],
5760  snf->ntransvars, QUAD_TO_DBL(transrhs) + QUAD_TO_DBL(rowcoeftimesbestslb), QUAD_TO_DBL(rowcoef), bestslb[i], QUAD_TO_DBL(transrhs));
5761  }
5762  else
5763  {
5764  SCIP_Real QUAD(rowcoefbinary);
5765  SCIP_Real varsolvalbinary;
5766  SCIP_Real QUAD(val);
5767  SCIP_Real QUAD(contsolval);
5768  SCIP_Real QUAD(rowcoeftimesvubconst);
5769  int vubvarprobidx;
5770 
5771  SCIP_VAR** vubvars = SCIPvarGetVubVars(var);
5772  SCIP_Real* vubconsts = SCIPvarGetVubConstants(var);
5773  SCIP_Real* vubcoefs = SCIPvarGetVubCoefs(var);
5774 
5775  /* use variable upper bound in bestslb = l_j <= y_j <= u~_j x_j + d_j = bestub to define
5776  * 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
5777  * 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,
5778  * where c_j is the coefficient of x_j in the row, put j into the set
5779  * N1 if a_j > 0
5780  * N2 if a_j < 0
5781  * and update the right hand side of the constraint in the relaxation
5782  * rhs = rhs - a_j d_j
5783  */
5784 
5785  vubvarprobidx = SCIPvarGetProbindex(vubvars[bestubtype[i]]);
5786  assert(binvarused[vubvarprobidx] == 1);
5787  assert(vubvarprobidx < nbinvars);
5788 
5789  QUAD_ARRAY_LOAD(rowcoefbinary, rowcoefs, vubvarprobidx);
5790  varsolvalbinary = SCIPgetSolVal(scip, sol, vubvars[bestubtype[i]]);
5791 
5792  /* clear the binvarpos array, since the variable has been processed */
5793  binvarused[vubvarprobidx] = 0;
5794 
5795  SCIPquadprecProdQD(val, rowcoef, vubcoefs[bestubtype[i]]);
5796  SCIPquadprecSumQQ(val, val, rowcoefbinary);
5797  {
5798  SCIP_Real QUAD(tmp);
5799  SCIPquadprecProdQD(tmp, rowcoefbinary, varsolvalbinary);
5800  SCIPquadprecSumDD(contsolval, solval, - vubconsts[bestubtype[i]]);
5801  SCIPquadprecProdQQ(contsolval, contsolval, rowcoef);
5802  SCIPquadprecSumQQ(contsolval, contsolval, tmp);
5803  }
5804 
5805  SCIPquadprecProdQD(rowcoeftimesvubconst, rowcoef, vubconsts[bestubtype[i]]);
5806  /* store aggregation information for y'_j for transforming cuts for the SNF relaxation back to the problem variables later */
5807  snf->origbinvars[snf->ntransvars] = vubvarprobidx;
5808 
5809  if( QUAD_TO_DBL(rowcoef) > QUAD_EPSILON )
5810  {
5811  snf->transvarcoefs[snf->ntransvars] = 1;
5812  snf->transvarvubcoefs[snf->ntransvars] = QUAD_TO_DBL(val);
5813  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5814  snf->transcontvarsolvals[snf->ntransvars] = QUAD_TO_DBL(contsolval);
5815 
5816  /* aggregation information for y'_j */
5817  snf->aggrcoefsbin[snf->ntransvars] = QUAD_TO_DBL(rowcoefbinary);
5818  snf->aggrcoefscont[snf->ntransvars] = QUAD_TO_DBL(rowcoef);
5819  snf->aggrconstants[snf->ntransvars] = - QUAD_TO_DBL(rowcoeftimesvubconst);
5820  }
5821  else
5822  {
5823  assert(QUAD_TO_DBL(rowcoef) < QUAD_EPSILON);
5824  snf->transvarcoefs[snf->ntransvars] = - 1;
5825  snf->transvarvubcoefs[snf->ntransvars] = - QUAD_TO_DBL(val);
5826  snf->transbinvarsolvals[snf->ntransvars] = varsolvalbinary;
5827  snf->transcontvarsolvals[snf->ntransvars] = - QUAD_TO_DBL(contsolval);
5828 
5829  /* aggregation information for y'_j */
5830  snf->aggrcoefsbin[snf->ntransvars] = - QUAD_TO_DBL(rowcoefbinary);
5831  snf->aggrcoefscont[snf->ntransvars] = - QUAD_TO_DBL(rowcoef);
5832  snf->aggrconstants[snf->ntransvars] = QUAD_TO_DBL(rowcoeftimesvubconst);
5833  }
5834  SCIPquadprecSumQQ(transrhs, transrhs, -rowcoeftimesvubconst);
5835 
5836  /* store for x_j that y'_j is the associated real variable in the 0-1 single node flow relaxation