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