Scippy

SCIP

Solving Constraint Integer Programs

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