Scippy

SCIP

Solving Constraint Integer Programs

nlhdlr_bilinear.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-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlhdlr_bilinear.c
17  * @ingroup DEFPLUGINS_NLHDLR
18  * @brief bilinear nonlinear handler
19  * @author Benjamin Mueller
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #include <string.h>
25 
26 #include "scip/nlhdlr_bilinear.h"
27 #include "scip/cons_nonlinear.h"
28 #include "scip/expr_product.h"
29 #include "scip/expr_var.h"
30 
31 /* fundamental nonlinear handler properties */
32 #define NLHDLR_NAME "bilinear"
33 #define NLHDLR_DESC "bilinear handler for expressions"
34 #define NLHDLR_DETECTPRIORITY -10 /**< it is important that the nlhdlr runs after the default nldhlr */
35 #define NLHDLR_ENFOPRIORITY -10
36 
37 #define MIN_INTERIORITY 0.01 /**< minimum interiority for a reference point for applying separation */
38 #define MIN_ABSBOUNDSIZE 0.1 /**< minimum size of variable bounds for applying separation */
39 
40 /* properties of the bilinear nlhdlr statistics table */
41 #define TABLE_NAME_BILINEAR "nlhdlr_bilinear"
42 #define TABLE_DESC_BILINEAR "bilinear nlhdlr statistics table"
43 #define TABLE_POSITION_BILINEAR 14800 /**< the position of the statistics table */
44 #define TABLE_EARLIEST_STAGE_BILINEAR SCIP_STAGE_INITSOLVE /**< output of the statistics table is only printed from this stage onwards */
45 
46 
47 /*
48  * Data structures
49  */
50 
51 /** nonlinear handler expression data */
52 struct SCIP_NlhdlrExprData
53 {
54  SCIP_Real underineqs[6]; /**< inequalities for underestimation */
55  int nunderineqs; /**< total number of inequalities for underestimation */
56  SCIP_Real overineqs[6]; /**< inequalities for overestimation */
57  int noverineqs; /**< total number of inequalities for overestimation */
58  SCIP_Longint lastnodeid; /**< id of the last node that has been used for separation */
59  int nseparoundslastnode; /**< number of separation calls of the last node */
60 };
61 
62 /** nonlinear handler data */
63 struct SCIP_NlhdlrData
64 {
65  SCIP_EXPR** exprs; /**< expressions that have been detected by the nlhdlr */
66  int nexprs; /**< total number of expression that have been detected */
67  int exprsize; /**< size of exprs array */
68  SCIP_HASHMAP* exprmap; /**< hashmap to store the position of each expression in the exprs array */
69 
70  /* parameter */
71  SCIP_Bool useinteval; /**< whether to use the interval evaluation callback of the nlhdlr */
72  SCIP_Bool usereverseprop; /**< whether to use the reverse propagation callback of the nlhdlr */
73  int maxseparoundsroot; /**< maximum number of separation rounds in the root node */
74  int maxseparounds; /**< maximum number of separation rounds in a local node */
75  int maxsepadepth; /**< maximum depth to apply separation */
76 };
77 
78 /*
79  * Local methods
80  */
81 
82 /** helper function to compute the violation of an inequality of the form xcoef * x <= ycoef * y + constant for two
83  * corner points of the domain [lbx,ubx] x [lby,uby]
84  */
85 static
87  SCIP_VAR* x, /**< first variable */
88  SCIP_VAR* y, /**< second variable */
89  SCIP_Real xcoef, /**< x-coefficient */
90  SCIP_Real ycoef, /**< y-coefficient */
91  SCIP_Real constant, /**< constant */
92  SCIP_Real* viol1, /**< buffer to store the violation of the first corner point */
93  SCIP_Real* viol2 /**< buffer to store the violation of the second corner point */
94  )
95 {
96  SCIP_Real norm;
97  assert(viol1 != NULL);
98  assert(viol2 != NULL);
99 
100  norm = SQRT(SQR(xcoef) + SQR(ycoef));
101 
102  /* inequality can be used for underestimating xy if and only if xcoef * ycoef > 0 */
103  if( xcoef * ycoef >= 0 )
104  {
105  /* violation for top-left and bottom-right corner */
106  *viol1 = MAX(0, (xcoef * SCIPvarGetLbLocal(x) - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/
107  *viol2 = MAX(0, (xcoef * SCIPvarGetUbLocal(x) - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/
108  }
109  else
110  {
111  /* violation for top-right and bottom-left corner */
112  *viol1 = MAX(0, (xcoef * SCIPvarGetUbLocal(x) - ycoef * SCIPvarGetUbLocal(y) - constant) / norm); /*lint !e666*/
113  *viol2 = MAX(0, (xcoef * SCIPvarGetLbLocal(x) - ycoef * SCIPvarGetLbLocal(y) - constant) / norm); /*lint !e666*/
114  }
115 }
116 
117 /** auxiliary function to decide whether to use inequalities for a strong relaxation of bilinear terms or not */
118 static
120  SCIP* scip, /**< SCIP data structure */
121  SCIP_VAR* x, /**< x variable */
122  SCIP_VAR* y, /**< y variable */
123  SCIP_Real refx, /**< reference point for x */
124  SCIP_Real refy /**< reference point for y */
125  )
126 {
127  SCIP_Real lbx;
128  SCIP_Real ubx;
129  SCIP_Real lby;
130  SCIP_Real uby;
131  SCIP_Real interiorityx;
132  SCIP_Real interiorityy;
133  SCIP_Real interiority;
134 
135  assert(x != NULL);
136  assert(y != NULL);
137  assert(x != y);
138 
139  /* get variable bounds */
140  lbx = SCIPvarGetLbLocal(x);
141  ubx = SCIPvarGetUbLocal(x);
142  lby = SCIPvarGetLbLocal(y);
143  uby = SCIPvarGetUbLocal(y);
144 
145  /* compute interiority */
146  interiorityx = MIN(refx-lbx, ubx-refx) / MAX(ubx-lbx, SCIPepsilon(scip)); /*lint !e666*/
147  interiorityy = MIN(refy-lby, uby-refy) / MAX(uby-lby, SCIPepsilon(scip)); /*lint !e666*/
148  interiority = 2.0*MIN(interiorityx, interiorityy);
149 
150  return ubx - lbx >= MIN_ABSBOUNDSIZE && uby - lby >= MIN_ABSBOUNDSIZE && interiority >= MIN_INTERIORITY;
151 }
152 
153 /** helper function to update the best relaxation for a bilinear term when using valid linear inequalities */
154 static
156  SCIP* scip, /**< SCIP data structure */
157  SCIP_VAR* RESTRICT x, /**< first variable */
158  SCIP_VAR* RESTRICT y, /**< second variable */
159  SCIP_Real bilincoef, /**< coefficient of the bilinear term */
160  SCIP_SIDETYPE violside, /**< side of quadratic constraint that is violated */
161  SCIP_Real refx, /**< reference point for the x variable */
162  SCIP_Real refy, /**< reference point for the y variable */
163  SCIP_Real* RESTRICT ineqs, /**< coefficients of each linear inequality; stored as triple (xcoef,ycoef,constant) */
164  int nineqs, /**< total number of inequalities */
165  SCIP_Real mccormickval, /**< value of the McCormick relaxation at the reference point */
166  SCIP_Real* RESTRICT bestcoefx, /**< pointer to update the x coefficient */
167  SCIP_Real* RESTRICT bestcoefy, /**< pointer to update the y coefficient */
168  SCIP_Real* RESTRICT bestconst, /**< pointer to update the constant */
169  SCIP_Real* RESTRICT bestval, /**< value of the best relaxation that have been found so far */
170  SCIP_Bool* success /**< buffer to store whether we found a better relaxation */
171  )
172 {
173  SCIP_Real constshift[2] = {0.0, 0.0};
174  SCIP_Real constant;
175  SCIP_Real xcoef;
176  SCIP_Real ycoef;
177  SCIP_Real lbx;
178  SCIP_Real ubx;
179  SCIP_Real lby;
180  SCIP_Real uby;
181  SCIP_Bool update;
182  SCIP_Bool overestimate;
183  int i;
184 
185  assert(x != y);
186  assert(!SCIPisZero(scip, bilincoef));
187  assert(nineqs >= 0 && nineqs <= 2);
188  assert(bestcoefx != NULL);
189  assert(bestcoefy != NULL);
190  assert(bestconst != NULL);
191  assert(bestval != NULL);
192 
193  /* no inequalities available */
194  if( nineqs == 0 )
195  return;
196  assert(ineqs != NULL);
197 
198  lbx = SCIPvarGetLbLocal(x);
199  ubx = SCIPvarGetUbLocal(x);
200  lby = SCIPvarGetLbLocal(y);
201  uby = SCIPvarGetUbLocal(y);
202  overestimate = (violside == SCIP_SIDETYPE_LEFT);
203 
204  /* check cases for which we can't compute a tighter relaxation */
205  if( SCIPisFeasLE(scip, refx, lbx) || SCIPisFeasGE(scip, refx, ubx)
206  || SCIPisFeasLE(scip, refy, lby) || SCIPisFeasGE(scip, refy, uby) )
207  return;
208 
209  /* due to the feasibility tolerances of the LP and NLP solver, it might possible that the reference point is
210  * violating the linear inequalities; to ensure that we compute a valid underestimate, we relax the linear
211  * inequality by changing its constant part
212  */
213  for( i = 0; i < nineqs; ++i )
214  {
215  constshift[i] = MAX(0.0, ineqs[3*i] * refx - ineqs[3*i+1] * refy - ineqs[3*i+2]);
216  SCIPdebugMsg(scip, "constant shift of inequality %d = %.16f\n", i, constshift[i]);
217  }
218 
219  /* try to use both inequalities */
220  if( nineqs == 2 )
221  {
222  SCIPcomputeBilinEnvelope2(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[0], ineqs[1],
223  ineqs[2] + constshift[0], ineqs[3], ineqs[4], ineqs[5] + constshift[1], &xcoef, &ycoef, &constant, &update);
224 
225  if( update )
226  {
227  SCIP_Real val = xcoef * refx + ycoef * refy + constant;
228  SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4) / (REALABS(*bestval - bilincoef * refx * refy) + 1e-4);
229  SCIP_Real absimpr = REALABS(val - (*bestval));
230 
231  /* update relaxation if possible */
232  if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
233  || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
234  {
235  *bestcoefx = xcoef;
236  *bestcoefy = ycoef;
237  *bestconst = constant;
238  *bestval = val;
239  *success = TRUE;
240  }
241  }
242  }
243 
244  /* use inequalities individually */
245  for( i = 0; i < nineqs; ++i )
246  {
247  SCIPcomputeBilinEnvelope1(scip, bilincoef, lbx, ubx, refx, lby, uby, refy, overestimate, ineqs[3*i], ineqs[3*i+1],
248  ineqs[3*i+2] + constshift[i], &xcoef, &ycoef, &constant, &update);
249 
250  if( update )
251  {
252  SCIP_Real val = xcoef * refx + ycoef * refy + constant;
253  SCIP_Real relimpr = 1.0 - (REALABS(val - bilincoef * refx * refy) + 1e-4)
254  / (REALABS(mccormickval - bilincoef * refx * refy) + 1e-4);
255  SCIP_Real absimpr = REALABS(val - (*bestval));
256 
257  /* update relaxation if possible */
258  if( relimpr > 0.05 && absimpr > 1e-3 && ((overestimate && SCIPisRelLT(scip, val, *bestval))
259  || (!overestimate && SCIPisRelGT(scip, val, *bestval))) )
260  {
261  *bestcoefx = xcoef;
262  *bestcoefy = ycoef;
263  *bestconst = constant;
264  *bestval = val;
265  *success = TRUE;
266  }
267  }
268  }
269 }
270 
271 /** helper function to determine whether a given point satisfy given inequalities */
272 static
274  SCIP* scip, /**< SCIP data structure */
275  SCIP_Real x, /**< x-coordinate */
276  SCIP_Real y, /**< y-coordinate */
277  SCIP_Real lbx, /**< lower bound of x */
278  SCIP_Real ubx, /**< upper bound of x */
279  SCIP_Real lby, /**< lower bound of y */
280  SCIP_Real uby, /**< upper bound of y */
281  SCIP_Real* ineqs, /**< inequalities of the form coefx x <= coefy y + constant */
282  int nineqs /**< total number of inequalities */
283  )
284 {
285  int i;
286 
287  assert(ineqs != NULL);
288  assert(nineqs > 0);
289 
290  /* check whether point satisfies the bounds */
291  if( SCIPisLT(scip, x, lbx) || SCIPisGT(scip, x, ubx)
292  || SCIPisLT(scip, y, lby) || SCIPisGT(scip, y, uby) )
293  return FALSE;
294 
295  /* check whether point satisfy the linear inequalities */
296  for( i = 0; i < nineqs; ++i )
297  {
298  SCIP_Real coefx = ineqs[3*i];
299  SCIP_Real coefy = ineqs[3*i+1];
300  SCIP_Real constant = ineqs[3*i+2];
301 
302  /* TODO check with an absolute comparison? */
303  if( SCIPisGT(scip, coefx*x - coefy*y - constant, 0.0) )
304  return FALSE;
305  }
306 
307  return TRUE;
308 }
309 
310 /** helper function for computing all vertices of the polytope described by the linear inequalities and the local
311  * extrema of the bilinear term along each inequality
312  *
313  * @note there are at most 22 points where the min/max can be achieved (given that there are at most 4 inequalities)
314  * - corners of [lbx,ubx]x[lby,uby] (4)
315  * - two intersection points for each inequality with the box (8)
316  * - global maximum / minimum on each inequality (4)
317  * - intersection between two inequalities (6)
318  */
319 static
321  SCIP* scip, /**< SCIP data structure */
322  SCIP_CONSHDLR* conshdlr, /**< constraint handler, if levelset == TRUE, otherwise can be NULL */
323  SCIP_EXPR* expr, /**< product expression */
324  SCIP_INTERVAL exprbounds, /**< bounds on product expression, only used if levelset == TRUE */
325  SCIP_Real* underineqs, /**< inequalities for underestimation */
326  int nunderineqs, /**< total number of inequalities for underestimation */
327  SCIP_Real* overineqs, /**< inequalities for overestimation */
328  int noverineqs, /**< total number of inequalities for overestimation */
329  SCIP_Bool levelset, /**< should the level set be considered? */
330  SCIP_Real* xs, /**< array to store x-coordinates of computed points */
331  SCIP_Real* ys, /**< array to store y-coordinates of computed points */
332  int* npoints /**< buffer to store the total number of computed points */
333  )
334 {
335  SCIP_EXPR* child1;
336  SCIP_EXPR* child2;
337  SCIP_Real ineqs[12];
338  SCIP_INTERVAL boundsx;
339  SCIP_INTERVAL boundsy;
340  SCIP_Real lbx;
341  SCIP_Real ubx;
342  SCIP_Real lby;
343  SCIP_Real uby;
344  int nineqs = 0;
345  int i;
346 
347  assert(scip != NULL);
348  assert(conshdlr != NULL || !levelset);
349  assert(expr != NULL);
350  assert(xs != NULL);
351  assert(ys != NULL);
352  assert(SCIPexprGetNChildren(expr) == 2);
353  assert(noverineqs + nunderineqs > 0);
354  assert(noverineqs + nunderineqs <= 4);
355 
356  *npoints = 0;
357 
358  /* collect inequalities */
359  for( i = 0; i < noverineqs; ++i )
360  {
361  SCIPdebugMsg(scip, "over-inequality %d: %g*x <= %g*y + %g\n", i, overineqs[3*i], overineqs[3*i+1], overineqs[3*i+2]);
362  ineqs[3*nineqs] = overineqs[3*i];
363  ineqs[3*nineqs+1] = overineqs[3*i+1];
364  ineqs[3*nineqs+2] = overineqs[3*i+2];
365  ++nineqs;
366  }
367  for( i = 0; i < nunderineqs; ++i )
368  {
369  SCIPdebugMsg(scip, "under-inequality %d: %g*x <= %g*y + %g 0\n", i, underineqs[3*i], underineqs[3*i+1], underineqs[3*i+2]);
370  ineqs[3*nineqs] = underineqs[3*i];
371  ineqs[3*nineqs+1] = underineqs[3*i+1];
372  ineqs[3*nineqs+2] = underineqs[3*i+2];
373  ++nineqs;
374  }
375  assert(nineqs == noverineqs + nunderineqs);
376 
377  /* collect children */
378  child1 = SCIPexprGetChildren(expr)[0];
379  child2 = SCIPexprGetChildren(expr)[1];
380  assert(child1 != NULL && child2 != NULL);
381  assert(child1 != child2);
382 
383  /* collect bounds of children */
384  if( !levelset )
385  {
386  /* if called from inteval, then use activity */
387  boundsx = SCIPexprGetActivity(child1);
388  boundsy = SCIPexprGetActivity(child2);
389  }
390  else
391  {
392  /* if called from reverseprop, then use bounds */
393  boundsx = SCIPgetExprBoundsNonlinear(scip, child1);
394  boundsy = SCIPgetExprBoundsNonlinear(scip, child2);
395 
396  /* if children bounds are empty, then returning with *npoints==0 is the way to go */
399  return;
400  }
401  lbx = boundsx.inf;
402  ubx = boundsx.sup;
403  lby = boundsy.inf;
404  uby = boundsy.sup;
405  SCIPdebugMsg(scip, "x = [%g,%g], y=[%g,%g]\n", lbx, ubx, lby, uby);
406 
407  /* corner points that satisfy all inequalities */
408  for( i = 0; i < 4; ++i )
409  {
410  SCIP_Real cx = i < 2 ? lbx : ubx;
411  SCIP_Real cy = (i % 2) == 0 ? lby : uby;
412 
413  SCIPdebugMsg(scip, "corner point (%g,%g) feasible? %d\n", cx, cy, isPointFeasible(scip, cx, cy, lbx, ubx, lby, uby, ineqs, nineqs));
414 
415  if( isPointFeasible(scip, cx, cy, lbx, ubx, lby, uby, ineqs, nineqs) )
416  {
417  xs[*npoints] = cx;
418  ys[*npoints] = cy;
419  ++(*npoints);
420  }
421  }
422 
423  /* intersection point of inequalities with [lbx,ubx] x [lby,uby] and extremum of xy on each inequality */
424  for( i = 0; i < nineqs; ++i )
425  {
426  SCIP_Real coefx = ineqs[3*i];
427  SCIP_Real coefy = ineqs[3*i+1];
428  SCIP_Real constant = ineqs[3*i+2];
429  SCIP_Real px[5] = {lbx, ubx, (coefy*lby + constant)/coefx, (coefy*uby + constant)/coefx, 0.0};
430  SCIP_Real py[5] = {(coefx*lbx - constant)/coefy, (coefx*ubx - constant)/coefy, lby, uby, 0.0};
431  int j;
432 
433  /* the last entry corresponds to the extremum of xy on the line */
434  py[4] = (-constant) / (2.0 * coefy);
435  px[4] = constant / (2.0 * coefx);
436 
437  for( j = 0; j < 5; ++j )
438  {
439  SCIPdebugMsg(scip, "intersection point (%g,%g) feasible? %d\n", px[j], py[j], isPointFeasible(scip, px[j], py[j], lbx, ubx, lby, uby, ineqs, nineqs));
440  if( isPointFeasible(scip, px[j], py[j], lbx, ubx, lby, uby, ineqs, nineqs) )
441  {
442  xs[*npoints] = px[j];
443  ys[*npoints] = py[j];
444  ++(*npoints);
445  }
446  }
447  }
448 
449  /* intersection point between two inequalities */
450  for( i = 0; i < nineqs - 1; ++i )
451  {
452  SCIP_Real coefx1 = ineqs[3*i];
453  SCIP_Real coefy1 = ineqs[3*i+1];
454  SCIP_Real constant1 = ineqs[3*i+2];
455  int j;
456 
457  for( j = i + 1; j < nineqs; ++j )
458  {
459  SCIP_Real coefx2 = ineqs[3*j];
460  SCIP_Real coefy2 = ineqs[3*j+1];
461  SCIP_Real constant2 = ineqs[3*j+2];
462  SCIP_Real px;
463  SCIP_Real py;
464 
465  /* no intersection point -> skip */
466  if( SCIPisZero(scip, coefx2*coefy1 - coefx1 * coefy2) )
467  continue;
468 
469  py = (constant2 * coefx1 - constant1 * coefx2)/ (coefx2 * coefy1 - coefx1 * coefy2);
470  px = (coefy1 * py + constant1) / coefx1;
471  assert(SCIPisRelEQ(scip, px, (coefy2 * py + constant2) / coefx2));
472 
473  if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
474  {
475  xs[*npoints] = px;
476  ys[*npoints] = py;
477  ++(*npoints);
478  }
479  }
480  }
481 
482  assert(*npoints <= 22);
483 
484  /* consider the intersection of the level set with
485  *
486  * 1. the boundary of the box
487  * 2. the linear inequalities
488  *
489  * this adds at most for 4 (level set curves) * 4 (inequalities) * 2 (intersection points) for all linear
490  * inequalities and 4 (level set curves) * 2 (intersection points) with the boundary of the box
491  */
492  if( !levelset )
493  return;
494 
495  /* compute intersection of level sets with the boundary */
496  for( i = 0; i < 2; ++i )
497  {
498  SCIP_Real vals[4] = {lbx, ubx, lby, uby};
499  SCIP_Real val;
500  int k;
501 
502  /* fix auxiliary variable to its lower or upper bound and consider the coefficient of the product */
503  val = (i == 0) ? exprbounds.inf : exprbounds.sup;
504  val /= SCIPgetCoefExprProduct(expr);
505 
506  for( k = 0; k < 4; ++k )
507  {
508  if( !SCIPisZero(scip, vals[k]) )
509  {
510  SCIP_Real res = val / vals[k];
511 
512  assert(SCIPisRelGE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.inf));
513  assert(SCIPisRelLE(scip, SCIPgetCoefExprProduct(expr)*res*vals[k], exprbounds.sup));
514 
515  /* fix x to lbx or ubx */
516  if( k < 2 && isPointFeasible(scip, vals[k], res, lbx, ubx, lby, uby, ineqs, nineqs) )
517  {
518  xs[*npoints] = vals[k];
519  ys[*npoints] = res;
520  ++(*npoints);
521  }
522  /* fix y to lby or uby */
523  else if( k >= 2 && isPointFeasible(scip, res, vals[k], lbx, ubx, lby, uby, ineqs, nineqs) )
524  {
525  xs[*npoints] = res;
526  ys[*npoints] = vals[k];
527  ++(*npoints);
528  }
529  }
530  }
531  }
532 
533  /* compute intersection points of level sets with the linear inequalities */
534  for( i = 0; i < nineqs; ++i )
535  {
536  SCIP_INTERVAL result;
537  SCIP_Real coefx = ineqs[3*i];
538  SCIP_Real coefy = ineqs[3*i+1];
539  SCIP_Real constant = ineqs[3*i+2];
540  SCIP_INTERVAL sqrcoef;
541  SCIP_INTERVAL lincoef;
542  SCIP_Real px;
543  SCIP_Real py;
544  int k;
545 
546  /* solve system of coefx x = coefy y + constant and X = xy which is the same as computing the solutions of
547  *
548  * (coefy / coefx) y^2 + (constant / coefx) y = inf(X) or sup(X)
549  */
550  SCIPintervalSet(&sqrcoef, coefy / coefx);
551  SCIPintervalSet(&lincoef, constant / coefx);
552 
553  for( k = 0; k < 2; ++k )
554  {
555  SCIP_INTERVAL rhs;
556  SCIP_INTERVAL ybnds;
557 
558  /* set right-hand side */
559  if( k == 0 )
560  SCIPintervalSet(&rhs, exprbounds.inf);
561  else
562  SCIPintervalSet(&rhs, exprbounds.sup);
563 
564  SCIPintervalSetBounds(&ybnds, lby, uby);
565  SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &result, sqrcoef, lincoef, rhs, ybnds);
566 
567  /* interval is empty -> no solution available */
569  continue;
570 
571  /* compute and check point */
572  py = SCIPintervalGetInf(result);
573  px = (coefy * py + constant) / coefx;
574 
575  if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
576  {
577  xs[*npoints] = px;
578  ys[*npoints] = py;
579  ++(*npoints);
580  }
581 
582  /* check for a second solution */
583  if( SCIPintervalGetInf(result) != SCIPintervalGetSup(result) ) /*lint !e777*/
584  {
585  py = SCIPintervalGetSup(result);
586  px = (coefy * py + constant) / coefx;
587 
588  if( isPointFeasible(scip, px, py, lbx, ubx, lby, uby, ineqs, nineqs) )
589  {
590  xs[*npoints] = px;
591  ys[*npoints] = py;
592  ++(*npoints);
593  }
594  }
595  }
596  }
597 
598  assert(*npoints <= 62);
599 }
600 
601 /** computes interval for a bilinear term when using at least one inequality */
602 static
604  SCIP* scip, /**< SCIP data structure */
605  SCIP_EXPR* expr, /**< product expression */
606  SCIP_Real* underineqs, /**< inequalities for underestimation */
607  int nunderineqs, /**< total number of inequalities for underestimation */
608  SCIP_Real* overineqs, /**< inequalities for overestimation */
609  int noverineqs /**< total number of inequalities for overestimation */
610  )
611 {
612  SCIP_INTERVAL interval = {0., 0.};
613  SCIP_Real xs[22];
614  SCIP_Real ys[22];
615  SCIP_Real inf;
616  SCIP_Real sup;
617  int npoints;
618  int i;
619 
620  assert(scip != NULL);
621  assert(expr != NULL);
622  assert(SCIPexprGetNChildren(expr) == 2);
623  assert(noverineqs + nunderineqs <= 4);
624 
625  /* no inequalities available -> skip computation */
626  if( noverineqs == 0 && nunderineqs == 0 )
627  {
629  return interval;
630  }
631 
632  /* x or y has empty interval -> empty */
635  {
636  SCIPintervalSetEmpty(&interval);
637  return interval;
638  }
639 
640  /* compute all feasible points (since we use levelset == FALSE, the value of interval doesn't matter) */
641  getFeasiblePointsBilinear(scip, NULL, expr, interval, underineqs, nunderineqs, overineqs,
642  noverineqs, FALSE, xs, ys, &npoints);
643 
644  /* no feasible point left -> return an empty interval */
645  if( npoints == 0 )
646  {
647  SCIPintervalSetEmpty(&interval);
648  return interval;
649  }
650 
651  /* compute the minimum and maximum over all computed points */
652  inf = xs[0] * ys[0];
653  sup = inf;
654  SCIPdebugMsg(scip, "point 0: (%g,%g) -> inf = sup = %g\n", xs[0], ys[0], inf);
655  for( i = 1; i < npoints; ++i )
656  {
657  inf = MIN(inf, xs[i] * ys[i]);
658  sup = MAX(sup, xs[i] * ys[i]);
659  SCIPdebugMsg(scip, "point %d: (%g,%g) -> inf = %g, sup = %g\n", i, xs[i], ys[i], inf, sup);
660  }
661  assert(inf <= sup);
662 
663  /* adjust infinite values */
664  inf = MAX(inf, -SCIP_INTERVAL_INFINITY);
665  sup = MIN(sup, SCIP_INTERVAL_INFINITY);
666 
667  /* multiply resulting interval with coefficient of the product expression */
668  SCIPintervalSetBounds(&interval, inf, sup);
669  if( SCIPgetCoefExprProduct(expr) != 1.0 )
671 
672  return interval;
673 }
674 
675 /** uses inequalities for bilinear terms to get stronger bounds during reverse propagation */
676 static
678  SCIP* scip, /**< SCIP data structure */
679  SCIP_CONSHDLR* conshdlr, /**< constraint handler */
680  SCIP_EXPR* expr, /**< product expression */
681  SCIP_INTERVAL exprbounds, /**< bounds on product expression */
682  SCIP_Real* underineqs, /**< inequalities for underestimation */
683  int nunderineqs, /**< total number of inequalities for underestimation */
684  SCIP_Real* overineqs, /**< inequalities for overestimation */
685  int noverineqs, /**< total number of inequalities for overestimation */
686  SCIP_INTERVAL* intervalx, /**< buffer to store the new interval for x */
687  SCIP_INTERVAL* intervaly /**< buffer to store the new interval for y */
688  )
689 {
690  SCIP_Real xs[62];
691  SCIP_Real ys[62];
692  SCIP_Real exprinf;
693  SCIP_Real exprsup;
694  SCIP_Bool first = TRUE;
695  int npoints;
696  int i;
697 
698  assert(scip != NULL);
699  assert(conshdlr != NULL);
700  assert(expr != NULL);
701  assert(intervalx != NULL);
702  assert(intervaly != NULL);
703  assert(SCIPexprGetNChildren(expr) == 2);
704 
705  assert(noverineqs + nunderineqs > 0);
706 
707  /* set intervals to be empty */
708  SCIPintervalSetEmpty(intervalx);
709  SCIPintervalSetEmpty(intervaly);
710 
711  /* compute feasible points */
712  getFeasiblePointsBilinear(scip, conshdlr, expr, exprbounds, underineqs, nunderineqs, overineqs,
713  noverineqs, TRUE, xs, ys, &npoints);
714 
715  /* no feasible points left -> problem is infeasible */
716  if( npoints == 0 )
717  return;
718 
719  /* get bounds of the product expression */
720  exprinf = exprbounds.inf;
721  exprsup = exprbounds.sup;
722 
723  /* update intervals with the computed points */
724  for( i = 0; i < npoints; ++i )
725  {
726  SCIP_Real val = SCIPgetCoefExprProduct(expr) * xs[i] * ys[i];
727 
728 #ifndef NDEBUG
729  {
734 
735  assert(nunderineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, underineqs, nunderineqs));
736  assert(noverineqs == 0 || isPointFeasible(scip, xs[i], ys[i], lbx, ubx, lby, uby, overineqs, noverineqs));
737  }
738 #endif
739 
740  /* only accept points for which the value of x*y is in the interval of the product expression
741  *
742  * NOTE: in order to consider all relevant points, we are a bit conservative here and relax the interval of
743  * the expression by SCIPfeastol()
744  */
745  if( SCIPisRelGE(scip, val, exprinf - SCIPfeastol(scip)) && SCIPisRelLE(scip, val, exprsup + SCIPfeastol(scip)) )
746  {
747  if( first )
748  {
749  SCIPintervalSet(intervalx, xs[i]);
750  SCIPintervalSet(intervaly, ys[i]);
751  first = FALSE;
752  }
753  else
754  {
755  (*intervalx).inf = MIN((*intervalx).inf, xs[i]);
756  (*intervalx).sup = MAX((*intervalx).sup, xs[i]);
757  (*intervaly).inf = MIN((*intervaly).inf, ys[i]);
758  (*intervaly).sup = MAX((*intervaly).sup, ys[i]);
759  }
760 
761  SCIPdebugMsg(scip, "consider points (%g,%g)=%g for reverse propagation\n", xs[i], ys[i], val);
762  }
763  }
764 }
765 
766 /** helper function to compute the convex envelope of a bilinear term when two linear inequalities are given; we
767  * use the same notation and formulas as in Locatelli 2016
768  */
769 static
771  SCIP* scip, /**< SCIP data structure */
772  SCIP_Real x, /**< reference point for x */
773  SCIP_Real y, /**< reference point for y */
774  SCIP_Real mi, /**< coefficient of x in the first linear inequality */
775  SCIP_Real qi, /**< constant in the first linear inequality */
776  SCIP_Real mj, /**< coefficient of x in the second linear inequality */
777  SCIP_Real qj, /**< constant in the second linear inequality */
778  SCIP_Real* RESTRICT xi, /**< buffer to store x coordinate of the first point */
779  SCIP_Real* RESTRICT yi, /**< buffer to store y coordinate of the first point */
780  SCIP_Real* RESTRICT xj, /**< buffer to store x coordinate of the second point */
781  SCIP_Real* RESTRICT yj, /**< buffer to store y coordinate of the second point */
782  SCIP_Real* RESTRICT xcoef, /**< buffer to store the x coefficient of the envelope */
783  SCIP_Real* RESTRICT ycoef, /**< buffer to store the y coefficient of the envelope */
784  SCIP_Real* RESTRICT constant /**< buffer to store the constant of the envelope */
785  )
786 {
787  SCIP_Real QUAD(xiq);
788  SCIP_Real QUAD(yiq);
789  SCIP_Real QUAD(xjq);
790  SCIP_Real QUAD(yjq);
791  SCIP_Real QUAD(xcoefq);
792  SCIP_Real QUAD(ycoefq);
793  SCIP_Real QUAD(constantq);
794  SCIP_Real QUAD(tmpq);
795 
796  assert(xi != NULL);
797  assert(yi != NULL);
798  assert(xj != NULL);
799  assert(yj != NULL);
800  assert(xcoef != NULL);
801  assert(ycoef != NULL);
802  assert(constant != NULL);
803 
804  if( SCIPisEQ(scip, mi, mj) )
805  {
806  /* xi = (x + mi * y - qi) / (2.0*mi) */
807  SCIPquadprecProdDD(xiq, mi, y);
808  SCIPquadprecSumQD(xiq, xiq, x);
809  SCIPquadprecSumQD(xiq, xiq, -qi);
810  SCIPquadprecDivQD(xiq, xiq, 2.0 * mi);
811  assert(EPSEQ((x + mi * y - qi) / (2.0*mi), QUAD_TO_DBL(xiq), 1e-3));
812 
813  /* yi = mi*(*xi) + qi */
814  SCIPquadprecProdQD(yiq, xiq, mi);
815  SCIPquadprecSumQD(yiq, yiq, qi);
816  assert(EPSEQ(mi*QUAD_TO_DBL(xiq) + qi, QUAD_TO_DBL(yiq), 1e-3));
817 
818  /* xj = (*xi) + (qi - qj)/ (2.0*mi) */
819  SCIPquadprecSumDD(xjq, qi, -qj);
820  SCIPquadprecDivQD(xjq, xjq, 2.0 * mi);
821  SCIPquadprecSumQQ(xjq, xjq, xiq);
822  assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj)/ (2.0*mi), QUAD_TO_DBL(xjq), 1e-3));
823 
824  /* yj = mj * (*xj) + qj */
825  SCIPquadprecProdQD(yjq, xjq, mj);
826  SCIPquadprecSumQD(yjq, yjq, qj);
827  assert(EPSEQ(mj * QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
828 
829  /* ycoef = (*xi) + (qi - qj) / (4.0*mi) note that this is wrong in Locatelli 2016 */
830  SCIPquadprecSumDD(ycoefq, qi, -qj);
831  SCIPquadprecDivQD(ycoefq, ycoefq, 4.0 * mi);
832  SCIPquadprecSumQQ(ycoefq, ycoefq, xiq);
833  assert(EPSEQ(QUAD_TO_DBL(xiq) + (qi - qj) / (4.0*mi), QUAD_TO_DBL(ycoefq), 1e-3));
834 
835  /* xcoef = 2.0*mi*(*xi) - mi * (*ycoef) + qi */
836  SCIPquadprecProdQD(xcoefq, xiq, 2.0 * mi);
837  SCIPquadprecProdQD(tmpq, ycoefq, -mi);
838  SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
839  SCIPquadprecSumQD(xcoefq, xcoefq, qi);
840  assert(EPSEQ(2.0*mi*QUAD_TO_DBL(xiq) - mi * QUAD_TO_DBL(ycoefq) + qi, QUAD_TO_DBL(xcoefq), 1e-3));
841 
842  /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
843  SCIPquadprecSquareQ(constantq, xjq);
844  SCIPquadprecProdQD(constantq, constantq, -mj);
845  SCIPquadprecProdQD(tmpq, ycoefq, -qj);
846  SCIPquadprecSumQQ(constantq, constantq, tmpq);
847  /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
848 
849  *xi = QUAD_TO_DBL(xiq);
850  *yi = QUAD_TO_DBL(yiq);
851  *xj = QUAD_TO_DBL(xjq);
852  *yj = QUAD_TO_DBL(yjq);
853  *ycoef = QUAD_TO_DBL(ycoefq);
854  *xcoef = QUAD_TO_DBL(xcoefq);
855  *constant = QUAD_TO_DBL(constantq);
856  }
857  else if( mi > 0.0 )
858  {
859  assert(mj > 0.0);
860 
861  /* xi = (y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)) */
862  SCIPquadprecProdDD(xiq, mi, mj);
863  SCIPquadprecSqrtQ(xiq, xiq);
864  SCIPquadprecProdQD(xiq, xiq, x);
865  SCIPquadprecSumQD(xiq, xiq, y);
866  SCIPquadprecSumQD(xiq, xiq, -qi); /* (y + SQRT(mi*mj)*x - qi) */
867  SCIPquadprecProdDD(tmpq, mi, mj);
868  SCIPquadprecSqrtQ(tmpq, tmpq);
869  SCIPquadprecSumQD(tmpq, tmpq, REALABS(mi)); /* REALABS(mi) + SQRT(mi*mj) */
870  SCIPquadprecDivQQ(xiq, xiq, tmpq);
871  assert(EPSEQ((y + SQRT(mi*mj)*x - qi) / (REALABS(mi) + SQRT(mi*mj)), QUAD_TO_DBL(xiq), 1e-3));
872 
873  /* yi = mi*(*xi) + qi */
874  SCIPquadprecProdQD(yiq, xiq, mi);
875  SCIPquadprecSumQD(yiq, yiq, qi);
876  assert(EPSEQ(mi*(QUAD_TO_DBL(xiq)) + qi, QUAD_TO_DBL(yiq), 1e-3));
877 
878  /* xj = (y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)) */
879  SCIPquadprecProdDD(xjq, mi, mj);
880  SCIPquadprecSqrtQ(xjq, xjq);
881  SCIPquadprecProdQD(xjq, xjq, x);
882  SCIPquadprecSumQD(xjq, xjq, y);
883  SCIPquadprecSumQD(xjq, xjq, -qj); /* (y + SQRT(mi*mj)*x - qj) */
884  SCIPquadprecProdDD(tmpq, mi, mj);
885  SCIPquadprecSqrtQ(tmpq, tmpq);
886  SCIPquadprecSumQD(tmpq, tmpq, REALABS(mj)); /* REALABS(mj) + SQRT(mi*mj) */
887  SCIPquadprecDivQQ(xjq, xjq, tmpq);
888  assert(EPSEQ((y + SQRT(mi*mj)*x - qj) / (REALABS(mj) + SQRT(mi*mj)), QUAD_TO_DBL(xjq), 1e-3));
889 
890  /* yj = mj*(*xj) + qj */
891  SCIPquadprecProdQD(yjq, xjq, mj);
892  SCIPquadprecSumQD(yjq, yjq, qj);
893  assert(EPSEQ(mj*QUAD_TO_DBL(xjq) + qj, QUAD_TO_DBL(yjq), 1e-3));
894 
895  /* ycoef = (2.0*mj*(*xj) + qj - 2.0*mi*(*xi) - qi) / (mj - mi) */
896  SCIPquadprecProdQD(ycoefq, xjq, 2.0 * mj);
897  SCIPquadprecSumQD(ycoefq, ycoefq, qj);
898  SCIPquadprecProdQD(tmpq, xiq, -2.0 * mi);
899  SCIPquadprecSumQQ(ycoefq, ycoefq, tmpq);
900  SCIPquadprecSumQD(ycoefq, ycoefq, -qi);
901  SCIPquadprecSumDD(tmpq, mj, -mi);
902  SCIPquadprecDivQQ(ycoefq, ycoefq, tmpq);
903  assert(EPSEQ((2.0*mj*QUAD_TO_DBL(xjq) + qj - 2.0*mi*QUAD_TO_DBL(xiq) - qi) / (mj - mi), QUAD_TO_DBL(ycoefq), 1e-3));
904 
905  /* xcoef = 2.0*mj*(*xj) + qj - mj*(*ycoef) */
906  SCIPquadprecProdQD(xcoefq, xjq, 2.0 * mj);
907  SCIPquadprecSumQD(xcoefq, xcoefq, qj);
908  SCIPquadprecProdQD(tmpq, ycoefq, -mj);
909  SCIPquadprecSumQQ(xcoefq, xcoefq, tmpq);
910  assert(EPSEQ(2.0*mj*QUAD_TO_DBL(xjq) + qj - mj*QUAD_TO_DBL(ycoefq), QUAD_TO_DBL(xcoefq), 1e-3));
911 
912  /* constant = -mj*SQR(*xj) - (*ycoef) * qj */
913  SCIPquadprecSquareQ(constantq, xjq);
914  SCIPquadprecProdQD(constantq, constantq, -mj);
915  SCIPquadprecProdQD(tmpq, ycoefq, -qj);
916  SCIPquadprecSumQQ(constantq, constantq, tmpq);
917  /* assert(EPSEQ(-mj*SQR(QUAD_TO_DBL(xjq)) - QUAD_TO_DBL(ycoefq) * qj, QUAD_TO_DBL(constantq), 1e-3)); */
918 
919  *xi = QUAD_TO_DBL(xiq);
920  *yi = QUAD_TO_DBL(yiq);
921  *xj = QUAD_TO_DBL(xjq);
922  *yj = QUAD_TO_DBL(yjq);
923  *ycoef = QUAD_TO_DBL(ycoefq);
924  *xcoef = QUAD_TO_DBL(xcoefq);
925  *constant = QUAD_TO_DBL(constantq);
926  }
927  else
928  {
929  assert(mi < 0.0 && mj < 0.0);
930 
931  /* apply variable transformation x = -x in case for overestimation */
932  computeBilinEnvelope2(scip, -x, y, -mi, qi, -mj, qj, xi, yi, xj, yj, xcoef, ycoef, constant);
933 
934  /* revert transformation; multiply cut by -1 and change -x by x */
935  *xi = -(*xi);
936  *xj = -(*xj);
937  *ycoef = -(*ycoef);
938  *constant = -(*constant);
939  }
940 }
941 
942 /** output method of statistics table to output file stream 'file' */
943 static
944 SCIP_DECL_TABLEOUTPUT(tableOutputBilinear)
945 { /*lint --e{715}*/
946  SCIP_NLHDLR* nlhdlr;
947  SCIP_NLHDLRDATA* nlhdlrdata;
948  SCIP_CONSHDLR* conshdlr;
949  SCIP_HASHMAP* hashmap;
950  SCIP_EXPRITER* it;
951  int resfound = 0;
952  int restotal = 0;
953  int c;
954 
955  conshdlr = SCIPfindConshdlr(scip, "nonlinear");
956  assert(conshdlr != NULL);
957  nlhdlr = SCIPfindNlhdlrNonlinear(conshdlr, NLHDLR_NAME);
958  assert(nlhdlr != NULL);
959  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
960  assert(nlhdlrdata != NULL);
961 
962  /* allocate memory */
963  SCIP_CALL( SCIPhashmapCreate(&hashmap, SCIPblkmem(scip), nlhdlrdata->nexprs) );
965 
966  for( c = 0; c < nlhdlrdata->nexprs; ++c )
967  {
968  assert(!SCIPhashmapExists(hashmap, nlhdlrdata->exprs[c]));
969  SCIP_CALL( SCIPhashmapInsertInt(hashmap, nlhdlrdata->exprs[c], 0) );
970  }
971 
972  /* count in how many constraints each expression is contained */
973  for( c = 0; c < SCIPconshdlrGetNConss(conshdlr); ++c )
974  {
975  SCIP_CONS* cons = SCIPconshdlrGetConss(conshdlr)[c];
976  SCIP_EXPR* expr;
977 
979 
980  for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) ) /*lint !e441*/
981  {
982  if( SCIPhashmapExists(hashmap, expr) )
983  {
984  int nuses = SCIPhashmapGetImageInt(hashmap, expr);
985  SCIP_CALL( SCIPhashmapSetImageInt(hashmap, expr, nuses + 1) );
986  }
987  }
988  }
989 
990  /* compute success ratio */
991  for( c = 0; c < nlhdlrdata->nexprs; ++c )
992  {
993  SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
994  int nuses;
995 
996  nuses = SCIPhashmapGetImageInt(hashmap, nlhdlrdata->exprs[c]);
997  assert(nuses > 0);
998 
999  nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, nlhdlrdata->exprs[c]);
1000  assert(nlhdlrexprdata != NULL);
1001 
1002  if( nlhdlrexprdata->nunderineqs > 0 || nlhdlrexprdata->noverineqs > 0 )
1003  resfound += nuses;
1004  restotal += nuses;
1005  }
1006 
1007  /* print statistics */
1008  SCIPinfoMessage(scip, file, "Bilinear Nlhdlr : %10s %10s\n", "#found", "#total");
1009  SCIPinfoMessage(scip, file, " %-17s:", "");
1010  SCIPinfoMessage(scip, file, " %10d", resfound);
1011  SCIPinfoMessage(scip, file, " %10d", restotal);
1012  SCIPinfoMessage(scip, file, "\n");
1013 
1014  /* free memory */
1015  SCIPfreeExpriter(&it);
1016  SCIPhashmapFree(&hashmap);
1017 
1018  return SCIP_OKAY;
1019 }
1020 
1021 
1022 /*
1023  * Callback methods of nonlinear handler
1024  */
1025 
1026 /** nonlinear handler copy callback */
1027 static
1028 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrBilinear)
1029 { /*lint --e{715}*/
1030  assert(targetscip != NULL);
1031  assert(sourcenlhdlr != NULL);
1032  assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1033 
1034  SCIP_CALL( SCIPincludeNlhdlrBilinear(targetscip) );
1035 
1036  return SCIP_OKAY;
1037 }
1038 
1039 /** callback to free data of handler */
1040 static
1041 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataBilinear)
1042 { /*lint --e{715}*/
1043  assert(nlhdlrdata != NULL);
1044  assert((*nlhdlrdata)->nexprs == 0);
1045 
1046  if( (*nlhdlrdata)->exprmap != NULL )
1047  {
1048  assert(SCIPhashmapGetNElements((*nlhdlrdata)->exprmap) == 0);
1049  SCIPhashmapFree(&(*nlhdlrdata)->exprmap);
1050  }
1051 
1052  SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrdata)->exprs, (*nlhdlrdata)->exprsize);
1053  SCIPfreeBlockMemory(scip, nlhdlrdata);
1054 
1055  return SCIP_OKAY;
1056 }
1057 
1058 /** callback to free expression specific data */
1059 static
1060 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataBilinear)
1061 { /*lint --e{715}*/
1062  SCIP_NLHDLRDATA* nlhdlrdata;
1063  int pos;
1064 
1065  assert(expr != NULL);
1066 
1067  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1068  assert(nlhdlrdata != NULL);
1069  assert(nlhdlrdata->nexprs > 0);
1070  assert(nlhdlrdata->exprs != NULL);
1071  assert(nlhdlrdata->exprmap != NULL);
1072  assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr));
1073 
1074  pos = SCIPhashmapGetImageInt(nlhdlrdata->exprmap, (void*)expr);
1075  assert(pos >= 0 && pos < nlhdlrdata->nexprs);
1076  assert(nlhdlrdata->exprs[pos] == expr);
1077 
1078  /* move the last expression to the free position */
1079  if( nlhdlrdata->nexprs > 0 && pos != nlhdlrdata->nexprs - 1 )
1080  {
1081  SCIP_EXPR* lastexpr = nlhdlrdata->exprs[nlhdlrdata->nexprs - 1];
1082  assert(expr != lastexpr);
1083  assert(SCIPhashmapExists(nlhdlrdata->exprmap, (void*)lastexpr));
1084 
1085  nlhdlrdata->exprs[pos] = lastexpr;
1086  nlhdlrdata->exprs[nlhdlrdata->nexprs - 1] = NULL;
1087  SCIP_CALL( SCIPhashmapSetImageInt(nlhdlrdata->exprmap, (void*)lastexpr, pos) );
1088  }
1089 
1090  /* remove expression from the nonlinear handler data */
1091  SCIP_CALL( SCIPhashmapRemove(nlhdlrdata->exprmap, (void*)expr) );
1092  SCIP_CALL( SCIPreleaseExpr(scip, &expr) );
1093  --nlhdlrdata->nexprs;
1094 
1095  /* free nonlinear handler expression data */
1096  SCIPfreeBlockMemoryNull(scip, nlhdlrexprdata);
1097 
1098  return SCIP_OKAY;
1099 }
1100 
1101 /** callback to be called in initialization */
1102 #define nlhdlrInitBilinear NULL
1103 
1104 /** callback to be called in deinitialization */
1105 static
1106 SCIP_DECL_NLHDLREXIT(nlhdlrExitBilinear)
1107 { /*lint --e{715}*/
1108  assert(SCIPnlhdlrGetData(nlhdlr) != NULL);
1109  assert(SCIPnlhdlrGetData(nlhdlr)->nexprs == 0);
1110 
1111  return SCIP_OKAY;
1112 }
1113 
1114 /** callback to detect structure in expression tree */
1115 static
1116 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectBilinear)
1117 { /*lint --e{715}*/
1118  SCIP_NLHDLRDATA* nlhdlrdata;
1119 
1120  assert(expr != NULL);
1121  assert(participating != NULL);
1122 
1123  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1124  assert(nlhdlrdata);
1125 
1126  /* only during solving will we have the extra inequalities that we rely on so much here */
1128  return SCIP_OKAY;
1129 
1130  /* check for product expressions with two children */
1131  if( SCIPisExprProduct(scip, expr) && SCIPexprGetNChildren(expr) == 2
1132  && (nlhdlrdata->exprmap == NULL || !SCIPhashmapExists(nlhdlrdata->exprmap, (void*)expr)) )
1133  {
1134  SCIP_EXPR** children;
1135  SCIP_Bool valid;
1136  int c;
1137 
1138  children = SCIPexprGetChildren(expr);
1139  assert(children != NULL);
1140 
1141  /* detection is only successful if both children will have auxiliary variable or are variables
1142  * that are not binary variables */
1143  valid = TRUE;
1144  for( c = 0; c < 2; ++c )
1145  {
1146  assert(children[c] != NULL);
1147  if( SCIPgetExprNAuxvarUsesNonlinear(children[c]) == 0 &&
1148  (!SCIPisExprVar(scip, children[c]) || SCIPvarIsBinary(SCIPgetVarExprVar(children[c]))) )
1149  {
1150  valid = FALSE;
1151  break;
1152  }
1153  }
1154 
1155  if( valid )
1156  {
1157  /* create expression data for the nonlinear handler */
1158  SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) );
1159  (*nlhdlrexprdata)->lastnodeid = -1;
1160 
1161  /* ensure that there is enough memory to store the detected expression */
1162  if( nlhdlrdata->exprsize < nlhdlrdata->nexprs + 1 )
1163  {
1164  int newsize = SCIPcalcMemGrowSize(scip, nlhdlrdata->nexprs + 1);
1165  assert(newsize > nlhdlrdata->exprsize);
1166 
1167  SCIP_CALL( SCIPreallocBlockMemoryArray(scip, &nlhdlrdata->exprs, nlhdlrdata->exprsize, newsize) );
1168  nlhdlrdata->exprsize = newsize;
1169  }
1170 
1171  /* create expression map, if not done so far */
1172  if( nlhdlrdata->exprmap == NULL )
1173  {
1174  SCIP_CALL( SCIPhashmapCreate(&nlhdlrdata->exprmap, SCIPblkmem(scip), SCIPgetNVars(scip)) );
1175  }
1176 
1177 #ifndef NDEBUG
1178  {
1179  int i;
1180 
1181  for( i = 0; i < nlhdlrdata->nexprs; ++i )
1182  assert(nlhdlrdata->exprs[i] != expr);
1183  }
1184 #endif
1185 
1186  /* add expression to nlhdlrdata and capture it */
1187  nlhdlrdata->exprs[nlhdlrdata->nexprs] = expr;
1188  SCIPcaptureExpr(expr);
1189  SCIP_CALL( SCIPhashmapInsertInt(nlhdlrdata->exprmap, (void*)expr, nlhdlrdata->nexprs) );
1190  ++nlhdlrdata->nexprs;
1191 
1192  /* tell children that we will use their auxvar and use its activity for both estimate and domain propagation */
1193  SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[0], TRUE, nlhdlrdata->useinteval
1194  || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1195  SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, children[1], TRUE, nlhdlrdata->useinteval
1196  || nlhdlrdata->usereverseprop, TRUE, TRUE) );
1197  }
1198  }
1199 
1200  if( *nlhdlrexprdata != NULL )
1201  {
1202  /* we want to join separation and domain propagation, if not disabled by parameter */
1203  *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1204  if( nlhdlrdata->useinteval || nlhdlrdata->usereverseprop )
1205  *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
1206  }
1207 
1208 #ifdef SCIP_DEBUG
1209  if( *participating )
1210  {
1211  SCIPdebugMsg(scip, "detected expr ");
1212  SCIPprintExpr(scip, expr, NULL);
1213  SCIPinfoMessage(scip, NULL, " participating: %d\n", *participating);
1214  }
1215 #endif
1216 
1217  return SCIP_OKAY;
1218 }
1219 
1220 /** auxiliary evaluation callback of nonlinear handler */
1221 static
1222 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxBilinear)
1223 { /*lint --e{715}*/
1224  SCIP_VAR* var1;
1225  SCIP_VAR* var2;
1226  SCIP_Real coef;
1227 
1228  assert(SCIPisExprProduct(scip, expr));
1229  assert(SCIPexprGetNChildren(expr) == 2);
1230 
1232  assert(var1 != NULL);
1234  assert(var2 != NULL);
1235  coef = SCIPgetCoefExprProduct(expr);
1236 
1237  *auxvalue = coef * SCIPgetSolVal(scip, sol, var1) * SCIPgetSolVal(scip, sol, var2);
1238 
1239  return SCIP_OKAY;
1240 }
1241 
1242 /** separation initialization method of a nonlinear handler (called during CONSINITLP) */
1243 #define nlhdlrInitSepaBilinear NULL
1244 
1245 /** separation deinitialization method of a nonlinear handler (called during CONSEXITSOL) */
1246 #define nlhdlrExitSepaBilinear NULL
1247 
1248 /** nonlinear handler separation callback */
1249 #define nlhdlrEnfoBilinear NULL
1250 
1251 /** nonlinear handler under/overestimation callback */
1252 static
1253 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateBilinear)
1254 { /*lint --e{715}*/
1255  SCIP_NLHDLRDATA* nlhdlrdata;
1256  SCIP_VAR* x;
1257  SCIP_VAR* y;
1258  SCIP_VAR* auxvar;
1259  SCIP_Real lincoefx = 0.0;
1260  SCIP_Real lincoefy = 0.0;
1261  SCIP_Real linconstant = 0.0;
1262  SCIP_Real refpointx;
1263  SCIP_Real refpointy;
1264  SCIP_Real violation;
1265  SCIP_Longint nodeid;
1266  SCIP_Bool mccsuccess = TRUE;
1267  SCIP_ROWPREP* rowprep;
1268 
1269  assert(rowpreps != NULL);
1270 
1271  *success = FALSE;
1272  *addedbranchscores = FALSE;
1273 
1274  /* check whether an inequality is available */
1275  if( nlhdlrexprdata->noverineqs == 0 && nlhdlrexprdata->nunderineqs == 0 )
1276  return SCIP_OKAY;
1277 
1278  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1279  assert(nlhdlrdata != NULL);
1280 
1282 
1283  /* update last node */
1284  if( nlhdlrexprdata->lastnodeid != nodeid )
1285  {
1286  nlhdlrexprdata->lastnodeid = nodeid;
1287  nlhdlrexprdata->nseparoundslastnode = 0;
1288  }
1289 
1290  /* update separation round */
1291  ++nlhdlrexprdata->nseparoundslastnode;
1292 
1293  /* check working limits */
1294  if( (SCIPgetDepth(scip) == 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparoundsroot)
1295  || (SCIPgetDepth(scip) > 0 && nlhdlrexprdata->nseparoundslastnode > nlhdlrdata->maxseparounds)
1296  || SCIPgetDepth(scip) > nlhdlrdata->maxsepadepth )
1297  return SCIP_OKAY;
1298 
1299  /* collect variables */
1301  assert(x != NULL);
1303  assert(y != NULL);
1304  auxvar = SCIPgetExprAuxVarNonlinear(expr);
1305  assert(auxvar != NULL);
1306 
1307  /* get and adjust the reference points */
1308  refpointx = MIN(MAX(SCIPgetSolVal(scip, sol, x), SCIPvarGetLbLocal(x)),SCIPvarGetUbLocal(x)); /*lint !e666*/
1309  refpointy = MIN(MAX(SCIPgetSolVal(scip, sol, y), SCIPvarGetLbLocal(y)),SCIPvarGetUbLocal(y)); /*lint !e666*/
1310  assert(SCIPisLE(scip, refpointx, SCIPvarGetUbLocal(x)) && SCIPisGE(scip, refpointx, SCIPvarGetLbLocal(x)));
1311  assert(SCIPisLE(scip, refpointy, SCIPvarGetUbLocal(y)) && SCIPisGE(scip, refpointy, SCIPvarGetLbLocal(y)));
1312 
1313  /* use McCormick inequalities to decide whether we want to separate or not */
1315  SCIPvarGetLbLocal(y), SCIPvarGetUbLocal(y), refpointy, overestimate, &lincoefx, &lincoefy, &linconstant,
1316  &mccsuccess);
1317 
1318  /* too large values in McCormick inequalities -> skip */
1319  if( !mccsuccess )
1320  return SCIP_OKAY;
1321 
1322  /* compute violation for the McCormick relaxation */
1323  violation = lincoefx * refpointx + lincoefy * refpointy + linconstant - SCIPgetSolVal(scip, sol, auxvar);
1324  if( overestimate )
1325  violation = -violation;
1326 
1327  /* only use a tighter relaxations if McCormick does not separate the reference point */
1328  if( SCIPisFeasLE(scip, violation, 0.0) && useBilinIneqs(scip, x, y, refpointx, refpointy) )
1329  {
1330  SCIP_Bool useoverestineq = SCIPgetCoefExprProduct(expr) > 0.0 ? overestimate : !overestimate;
1331  SCIP_Real mccormickval = lincoefx * refpointx + lincoefy * refpointy + linconstant;
1332  SCIP_Real* ineqs;
1333  SCIP_Real bestval;
1334  int nineqs;
1335 
1336  /* McCormick relaxation is too weak */
1337  bestval = mccormickval;
1338 
1339  /* get the inequalities that might lead to a tighter relaxation */
1340  if( useoverestineq )
1341  {
1342  ineqs = nlhdlrexprdata->overineqs;
1343  nineqs = nlhdlrexprdata->noverineqs;
1344  }
1345  else
1346  {
1347  ineqs = nlhdlrexprdata->underineqs;
1348  nineqs = nlhdlrexprdata->nunderineqs;
1349  }
1350 
1351  /* use linear inequalities to update relaxation */
1353  overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT,
1354  refpointx, refpointy, ineqs, nineqs, mccormickval,
1355  &lincoefx, &lincoefy, &linconstant, &bestval,
1356  success);
1357 
1358 #ifndef NDEBUG
1359  /* check whether cut is really valid */
1360  if( *success )
1361  {
1362  assert(!overestimate || SCIPisLE(scip, bestval, mccormickval));
1363  assert(overestimate || SCIPisGE(scip, bestval, mccormickval));
1364  }
1365 #endif
1366  }
1367 
1368  if( *success )
1369  {
1371  SCIProwprepAddConstant(rowprep, linconstant);
1372  SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, 2) );
1373  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, x, lincoefx) );
1374  SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, y, lincoefy) );
1375  SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
1376  }
1377 
1378  return SCIP_OKAY;
1379 }
1380 
1381 /** nonlinear handler interval evaluation callback */
1382 static
1383 SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalBilinear)
1384 { /*lint --e{715}*/
1385  SCIP_NLHDLRDATA* nlhdlrdata;
1386  assert(nlhdlrexprdata != NULL);
1387 
1388  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1389  assert(nlhdlrdata != NULL);
1390 
1391  if( nlhdlrdata->useinteval && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1392  {
1393  SCIP_INTERVAL tmp = intevalBilinear(scip, expr, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1394  nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs);
1395 
1396  /* intersect intervals if we have learned a tighter interval */
1397  if( SCIPisGT(scip, tmp.inf, (*interval).inf) || SCIPisLT(scip, tmp.sup, (*interval).sup) )
1398  SCIPintervalIntersect(interval, *interval, tmp);
1399  }
1400 
1401  return SCIP_OKAY;
1402 }
1403 
1404 /** nonlinear handler callback for reverse propagation */
1405 static
1406 SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropBilinear)
1407 { /*lint --e{715}*/
1408  SCIP_NLHDLRDATA* nlhdlrdata;
1409 
1410  assert(nlhdlrexprdata != NULL);
1411 
1412  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1413  assert(nlhdlrdata != NULL);
1414 
1415  if( nlhdlrdata->usereverseprop && nlhdlrexprdata->nunderineqs + nlhdlrexprdata->noverineqs > 0 )
1416  {
1417  SCIP_EXPR* childx;
1418  SCIP_EXPR* childy;
1419  SCIP_INTERVAL intervalx;
1420  SCIP_INTERVAL intervaly;
1421 
1422  assert(SCIPexprGetNChildren(expr) == 2);
1423  childx = SCIPexprGetChildren(expr)[0];
1424  childy = SCIPexprGetChildren(expr)[1];
1425  assert(childx != NULL && childy != NULL);
1426 
1429 
1430  /* compute bounds on x and y */
1431  reversePropBilinear(scip, conshdlr, expr, bounds, nlhdlrexprdata->underineqs, nlhdlrexprdata->nunderineqs,
1432  nlhdlrexprdata->overineqs, nlhdlrexprdata->noverineqs, &intervalx, &intervaly);
1433 
1434  /* tighten bounds of x */
1435  SCIPdebugMsg(scip, "try to tighten bounds of x: [%g,%g] -> [%g,%g]\n",
1437  intervalx.inf, intervalx.sup);
1438 
1439  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, SCIPexprGetChildren(expr)[0], intervalx, infeasible,
1440  nreductions) );
1441 
1442  if( !(*infeasible) )
1443  {
1444  /* tighten bounds of y */
1445  SCIPdebugMsg(scip, "try to tighten bounds of y: [%g,%g] -> [%g,%g]\n",
1447  intervaly.inf, intervaly.sup);
1449  infeasible, nreductions) );
1450  }
1451  }
1452 
1453  return SCIP_OKAY;
1454 }
1455 
1456 /*
1457  * nonlinear handler specific interface methods
1458  */
1459 
1460 /** includes bilinear nonlinear handler in nonlinear constraint handler */
1462  SCIP* scip /**< SCIP data structure */
1463  )
1464 {
1465  SCIP_NLHDLRDATA* nlhdlrdata;
1466  SCIP_NLHDLR* nlhdlr;
1467 
1468  assert(scip != NULL);
1469 
1470  /**! [SnippetIncludeNlhdlrBilinear] */
1471  /* create nonlinear handler specific data */
1472  SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) );
1473  BMSclearMemory(nlhdlrdata);
1474 
1476  NLHDLR_ENFOPRIORITY, nlhdlrDetectBilinear, nlhdlrEvalauxBilinear, nlhdlrdata) );
1477  assert(nlhdlr != NULL);
1478 
1479  SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrBilinear);
1480  SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrFreehdlrdataBilinear);
1481  SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataBilinear);
1482  SCIPnlhdlrSetInitExit(nlhdlr, nlhdlrInitBilinear, nlhdlrExitBilinear);
1484  SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalBilinear, nlhdlrReversepropBilinear);
1485 
1486  /* parameters */
1487  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/useinteval",
1488  "whether to use the interval evaluation callback of the nlhdlr",
1489  &nlhdlrdata->useinteval, FALSE, TRUE, NULL, NULL) );
1490 
1491  SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" NLHDLR_NAME "/usereverseprop",
1492  "whether to use the reverse propagation callback of the nlhdlr",
1493  &nlhdlrdata->usereverseprop, FALSE, TRUE, NULL, NULL) );
1494 
1495  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparoundsroot",
1496  "maximum number of separation rounds in the root node",
1497  &nlhdlrdata->maxseparoundsroot, FALSE, 10, 0, INT_MAX, NULL, NULL) );
1498 
1499  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxseparounds",
1500  "maximum number of separation rounds in a local node",
1501  &nlhdlrdata->maxseparounds, FALSE, 1, 0, INT_MAX, NULL, NULL) );
1502 
1503  SCIP_CALL( SCIPaddIntParam(scip, "nlhdlr/" NLHDLR_NAME "/maxsepadepth",
1504  "maximum depth to apply separation",
1505  &nlhdlrdata->maxsepadepth, FALSE, INT_MAX, 0, INT_MAX, NULL, NULL) );
1506 
1507  /* statistic table */
1508  assert(SCIPfindTable(scip, TABLE_NAME_BILINEAR) == NULL);
1510  NULL, NULL, NULL, NULL, NULL, NULL, tableOutputBilinear,
1512  /**! [SnippetIncludeNlhdlrBilinear] */
1513 
1514  return SCIP_OKAY;
1515 }
1516 
1517 /** returns an array of expressions that have been detected by the bilinear nonlinear handler */
1519  SCIP_NLHDLR* nlhdlr /**< nonlinear handler */
1520  )
1521 {
1522  SCIP_NLHDLRDATA* nlhdlrdata;
1523 
1524  assert(nlhdlr != NULL);
1525  assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1526 
1527  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1528  assert(nlhdlrdata);
1529 
1530  return nlhdlrdata->exprs;
1531 }
1532 
1533 /** returns the total number of expressions that have been detected by the bilinear nonlinear handler */
1535  SCIP_NLHDLR* nlhdlr /**< nonlinear handler */
1536  )
1537 {
1538  SCIP_NLHDLRDATA* nlhdlrdata;
1539 
1540  assert(nlhdlr != NULL);
1541  assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1542 
1543  nlhdlrdata = SCIPnlhdlrGetData(nlhdlr);
1544  assert(nlhdlrdata);
1545 
1546  return nlhdlrdata->nexprs;
1547 }
1548 
1549 /** adds a globally valid inequality of the form \f$\text{xcoef}\cdot x \leq \text{ycoef} \cdot y + \text{constant}\f$ to a product expression of the form \f$x\cdot y\f$ */
1551  SCIP* scip, /**< SCIP data structure */
1552  SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */
1553  SCIP_EXPR* expr, /**< product expression */
1554  SCIP_Real xcoef, /**< x coefficient */
1555  SCIP_Real ycoef, /**< y coefficient */
1556  SCIP_Real constant, /**< constant part */
1557  SCIP_Bool* success /**< buffer to store whether inequality has been accepted */
1558  )
1559 {
1560  SCIP_NLHDLREXPRDATA* nlhdlrexprdata;
1561  SCIP_VAR* x;
1562  SCIP_VAR* y;
1563  SCIP_Real* ineqs;
1564  SCIP_Real viol1;
1565  SCIP_Real viol2;
1566  SCIP_Bool underestimate;
1567  int nineqs;
1568  int i;
1569 
1570  assert(scip != NULL);
1571  assert(nlhdlr != NULL);
1572  assert(strcmp(SCIPnlhdlrGetName(nlhdlr), NLHDLR_NAME) == 0);
1573  assert(expr != NULL);
1574  assert(SCIPexprGetNChildren(expr) == 2);
1575  assert(xcoef != SCIP_INVALID); /*lint !e777 */
1576  assert(ycoef != SCIP_INVALID); /*lint !e777 */
1577  assert(constant != SCIP_INVALID); /*lint !e777 */
1578  assert(success != NULL);
1579 
1580  *success = FALSE;
1581 
1582  /* find nonlinear handler expression handler data */
1583  nlhdlrexprdata = SCIPgetNlhdlrExprDataNonlinear(nlhdlr, expr);
1584 
1585  if( nlhdlrexprdata == NULL )
1586  {
1587  SCIPwarningMessage(scip, "nonlinear expression data has not been found. "
1588  "Skip SCIPaddConsExprExprProductBilinearIneq()\n");
1589  return SCIP_OKAY;
1590  }
1591 
1592  /* ignore inequalities that only yield to a (possible) bound tightening */
1593  if( SCIPisFeasZero(scip, xcoef) || SCIPisFeasZero(scip, ycoef) )
1594  return SCIP_OKAY;
1595 
1596  /* collect variables */
1599  assert(x != NULL);
1600  assert(y != NULL);
1601  assert(x != y);
1602 
1603  /* normalize inequality s.t. xcoef in {-1,1} */
1604  if( !SCIPisEQ(scip, REALABS(xcoef), 1.0) )
1605  {
1606  constant /= REALABS(xcoef);
1607  ycoef /= REALABS(xcoef);
1608  xcoef /= REALABS(xcoef);
1609  }
1610 
1611  /* coefficients of the inequality determine whether the inequality can be used for under- or overestimation */
1612  underestimate = xcoef * ycoef > 0;
1613 
1614  SCIPdebugMsg(scip, "add inequality for a bilinear term: %g %s <= %g %s + %g (underestimate=%u)\n", xcoef,
1615  SCIPvarGetName(x), ycoef, SCIPvarGetName(y), constant, underestimate);
1616 
1617  /* compute violation of the inequality of the important corner points */
1618  getIneqViol(x, y, xcoef, ycoef, constant, &viol1, &viol2);
1619  SCIPdebugMsg(scip, "violations of inequality = (%g,%g)\n", viol1, viol2);
1620 
1621  /* inequality does not cutoff one of the important corner points -> skip */
1622  if( SCIPisFeasLE(scip, MAX(viol1, viol2), 0.0) )
1623  return SCIP_OKAY;
1624 
1625  if( underestimate )
1626  {
1627  ineqs = nlhdlrexprdata->underineqs;
1628  nineqs = nlhdlrexprdata->nunderineqs;
1629  }
1630  else
1631  {
1632  ineqs = nlhdlrexprdata->overineqs;
1633  nineqs = nlhdlrexprdata->noverineqs;
1634  }
1635 
1636  /* check for a duplicate */
1637  for( i = 0; i < nineqs; ++i )
1638  {
1639  if( SCIPisFeasEQ(scip, xcoef, ineqs[3*i]) && SCIPisFeasEQ(scip, ycoef, ineqs[3*i+1])
1640  && SCIPisFeasEQ(scip, constant, ineqs[3*i+2]) )
1641  {
1642  SCIPdebugMsg(scip, "inequality already found -> skip\n");
1643  return SCIP_OKAY;
1644  }
1645  }
1646 
1647  {
1648  SCIP_Real viols1[2] = {0.0, 0.0};
1649  SCIP_Real viols2[2] = {0.0, 0.0};
1650 
1651  /* compute violations of existing inequalities */
1652  for( i = 0; i < nineqs; ++i )
1653  {
1654  getIneqViol(x, y, ineqs[3*i], ineqs[3*i+1], ineqs[3*i+2], &viols1[i], &viols2[i]);
1655 
1656  /* check whether an existing inequality is dominating the candidate */
1657  if( SCIPisGE(scip, viols1[i], viol1) && SCIPisGE(scip, viols2[i], viol2) )
1658  {
1659  SCIPdebugMsg(scip, "inequality is dominated by %d -> skip\n", i);
1660  return SCIP_OKAY;
1661  }
1662 
1663  /* replace inequality if candidate is dominating it */
1664  if( SCIPisLT(scip, viols1[i], viol1) && SCIPisLT(scip, viols2[i], viol2) )
1665  {
1666  SCIPdebugMsg(scip, "inequality dominates %d -> replace\n", i);
1667  ineqs[3*i] = xcoef;
1668  ineqs[3*i+1] = ycoef;
1669  ineqs[3*i+2] = constant;
1670  *success = TRUE;
1671  }
1672  }
1673 
1674  /* inequality is not dominated by other inequalities -> add if we have less than 2 inequalities */
1675  if( nineqs < 2 )
1676  {
1677  ineqs[3*nineqs] = xcoef;
1678  ineqs[3*nineqs + 1] = ycoef;
1679  ineqs[3*nineqs + 2] = constant;
1680  *success = TRUE;
1681  SCIPdebugMsg(scip, "add inequality\n");
1682 
1683  /* increase number of inequalities */
1684  if( underestimate )
1685  ++(nlhdlrexprdata->nunderineqs);
1686  else
1687  ++(nlhdlrexprdata->noverineqs);
1688  }
1689  }
1690 
1691  if( *success )
1692  {
1693  /* With the added inequalities, we can potentially compute tighter activities for the expression,
1694  * so constraints that contain this expression should be propagated again.
1695  * We don't have a direct expression to constraint mapping, though. This call marks all expr-constraints
1696  * which include any of the variables that this expression depends on for propagation.
1697  */
1699  }
1700 
1701  return SCIP_OKAY;
1702 }
1703 
1704 /** computes coefficients of linearization of a bilinear term in a reference point */
1706  SCIP* scip, /**< SCIP data structure */
1707  SCIP_Real bilincoef, /**< coefficient of bilinear term */
1708  SCIP_Real refpointx, /**< point where to linearize first variable */
1709  SCIP_Real refpointy, /**< point where to linearize second variable */
1710  SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
1711  SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
1712  SCIP_Real* linconstant, /**< buffer to add constant of linearization */
1713  SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
1714  )
1715 {
1716  SCIP_Real constant;
1717 
1718  assert(scip != NULL);
1719  assert(lincoefx != NULL);
1720  assert(lincoefy != NULL);
1721  assert(linconstant != NULL);
1722  assert(success != NULL);
1723 
1724  if( bilincoef == 0.0 )
1725  return;
1726 
1727  if( SCIPisInfinity(scip, REALABS(refpointx)) || SCIPisInfinity(scip, REALABS(refpointy)) )
1728  {
1729  *success = FALSE;
1730  return;
1731  }
1732 
1733  /* bilincoef * x * y -> bilincoef * (refpointx * refpointy + refpointy * (x - refpointx) + refpointx * (y - refpointy))
1734  * = -bilincoef * refpointx * refpointy + bilincoef * refpointy * x + bilincoef * refpointx * y
1735  */
1736 
1737  constant = -bilincoef * refpointx * refpointy;
1738 
1739  if( SCIPisInfinity(scip, REALABS(bilincoef * refpointx)) || SCIPisInfinity(scip, REALABS(bilincoef * refpointy))
1740  || SCIPisInfinity(scip, REALABS(constant)) )
1741  {
1742  *success = FALSE;
1743  return;
1744  }
1745 
1746  *lincoefx += bilincoef * refpointy;
1747  *lincoefy += bilincoef * refpointx;
1748  *linconstant += constant;
1749 }
1750 
1751 /** computes coefficients of McCormick under- or overestimation of a bilinear term */
1753  SCIP* scip, /**< SCIP data structure */
1754  SCIP_Real bilincoef, /**< coefficient of bilinear term */
1755  SCIP_Real lbx, /**< lower bound on first variable */
1756  SCIP_Real ubx, /**< upper bound on first variable */
1757  SCIP_Real refpointx, /**< reference point for first variable */
1758  SCIP_Real lby, /**< lower bound on second variable */
1759  SCIP_Real uby, /**< upper bound on second variable */
1760  SCIP_Real refpointy, /**< reference point for second variable */
1761  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
1762  SCIP_Real* lincoefx, /**< buffer to add coefficient of first variable in linearization */
1763  SCIP_Real* lincoefy, /**< buffer to add coefficient of second variable in linearization */
1764  SCIP_Real* linconstant, /**< buffer to add constant of linearization */
1765  SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
1766  )
1767 {
1768  SCIP_Real constant;
1769  SCIP_Real coefx;
1770  SCIP_Real coefy;
1771 
1772  assert(scip != NULL);
1773  assert(!SCIPisInfinity(scip, lbx));
1774  assert(!SCIPisInfinity(scip, -ubx));
1775  assert(!SCIPisInfinity(scip, lby));
1776  assert(!SCIPisInfinity(scip, -uby));
1777  assert(SCIPisInfinity(scip, -lbx) || SCIPisLE(scip, lbx, ubx));
1778  assert(SCIPisInfinity(scip, -lby) || SCIPisLE(scip, lby, uby));
1779  assert(SCIPisInfinity(scip, -lbx) || SCIPisLE(scip, lbx, refpointx));
1780  assert(SCIPisInfinity(scip, -lby) || SCIPisLE(scip, lby, refpointy));
1781  assert(SCIPisInfinity(scip, ubx) || SCIPisGE(scip, ubx, refpointx));
1782  assert(SCIPisInfinity(scip, uby) || SCIPisGE(scip, uby, refpointy));
1783  assert(lincoefx != NULL);
1784  assert(lincoefy != NULL);
1785  assert(linconstant != NULL);
1786  assert(success != NULL);
1787 
1788  if( bilincoef == 0.0 )
1789  return;
1790 
1791  if( overestimate )
1792  bilincoef = -bilincoef;
1793 
1794  if( SCIPisRelEQ(scip, lbx, ubx) && SCIPisRelEQ(scip, lby, uby) )
1795  {
1796  /* both x and y are mostly fixed */
1797  SCIP_Real cand1;
1798  SCIP_Real cand2;
1799  SCIP_Real cand3;
1800  SCIP_Real cand4;
1801 
1802  coefx = 0.0;
1803  coefy = 0.0;
1804 
1805  /* estimate x * y by constant */
1806  cand1 = lbx * lby;
1807  cand2 = lbx * uby;
1808  cand3 = ubx * lby;
1809  cand4 = ubx * uby;
1810 
1811  /* take most conservative value for underestimator */
1812  if( bilincoef < 0.0 )
1813  constant = bilincoef * MAX( MAX(cand1, cand2), MAX(cand3, cand4) );
1814  else
1815  constant = bilincoef * MIN( MIN(cand1, cand2), MIN(cand3, cand4) );
1816  }
1817  else if( bilincoef > 0.0 )
1818  {
1819  /* either x or y is not fixed and coef > 0.0 */
1820  if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, -lby) &&
1821  (SCIPisInfinity(scip, ubx) || SCIPisInfinity(scip, uby)
1822  || (uby - refpointy) * (ubx - refpointx) >= (refpointy - lby) * (refpointx - lbx)) )
1823  {
1824  if( SCIPisRelEQ(scip, lbx, ubx) )
1825  {
1826  /* x*y = lbx * y + (x-lbx) * y >= lbx * y + (x-lbx) * lby >= lbx * y + min{(ubx-lbx) * lby, 0 * lby} */
1827  coefx = 0.0;
1828  coefy = bilincoef * lbx;
1829  constant = bilincoef * (lby < 0.0 ? (ubx-lbx) * lby : 0.0);
1830  }
1831  else if( SCIPisRelEQ(scip, lby, uby) )
1832  {
1833  /* x*y = lby * x + (y-lby) * x >= lby * x + (y-lby) * lbx >= lby * x + min{(uby-lby) * lbx, 0 * lbx} */
1834  coefx = bilincoef * lby;
1835  coefy = 0.0;
1836  constant = bilincoef * (lbx < 0.0 ? (uby-lby) * lbx : 0.0);
1837  }
1838  else
1839  {
1840  coefx = bilincoef * lby;
1841  coefy = bilincoef * lbx;
1842  constant = -bilincoef * lbx * lby;
1843  }
1844  }
1845  else if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, uby) )
1846  {
1847  if( SCIPisRelEQ(scip, lbx, ubx) )
1848  {
1849  /* x*y = ubx * y + (x-ubx) * y >= ubx * y + (x-ubx) * uby >= ubx * y + min{(lbx-ubx) * uby, 0 * uby} */
1850  coefx = 0.0;
1851  coefy = bilincoef * ubx;
1852  constant = bilincoef * (uby > 0.0 ? (lbx - ubx) * uby : 0.0);
1853  }
1854  else if( SCIPisRelEQ(scip, lby, uby) )
1855  {
1856  /* x*y = uby * x + (y-uby) * x >= uby * x + (y-uby) * ubx >= uby * x + min{(lby-uby) * ubx, 0 * ubx} */
1857  coefx = bilincoef * uby;
1858  coefy = 0.0;
1859  constant = bilincoef * (ubx > 0.0 ? (lby - uby) * ubx : 0.0);
1860  }
1861  else
1862  {
1863  coefx = bilincoef * uby;
1864  coefy = bilincoef * ubx;
1865  constant = -bilincoef * ubx * uby;
1866  }
1867  }
1868  else
1869  {
1870  *success = FALSE;
1871  return;
1872  }
1873  }
1874  else
1875  {
1876  /* either x or y is not fixed and coef < 0.0 */
1877  if( !SCIPisInfinity(scip, ubx) && !SCIPisInfinity(scip, -lby) &&
1878  (SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, uby)
1879  || (ubx - lbx) * (refpointy - lby) <= (uby - lby) * (refpointx - lbx)) )
1880  {
1881  if( SCIPisRelEQ(scip, lbx, ubx) )
1882  {
1883  /* x*y = ubx * y + (x-ubx) * y <= ubx * y + (x-ubx) * lby <= ubx * y + max{(lbx-ubx) * lby, 0 * lby} */
1884  coefx = 0.0;
1885  coefy = bilincoef * ubx;
1886  constant = bilincoef * (lby < 0.0 ? (lbx - ubx) * lby : 0.0);
1887  }
1888  else if( SCIPisRelEQ(scip, lby, uby) )
1889  {
1890  /* x*y = lby * x + (y-lby) * x <= lby * x + (y-lby) * ubx <= lby * x + max{(uby-lby) * ubx, 0 * ubx} */
1891  coefx = bilincoef * lby;
1892  coefy = 0.0;
1893  constant = bilincoef * (ubx > 0.0 ? (uby - lby) * ubx : 0.0);
1894  }
1895  else
1896  {
1897  coefx = bilincoef * lby;
1898  coefy = bilincoef * ubx;
1899  constant = -bilincoef * ubx * lby;
1900  }
1901  }
1902  else if( !SCIPisInfinity(scip, -lbx) && !SCIPisInfinity(scip, uby) )
1903  {
1904  if( SCIPisRelEQ(scip, lbx, ubx) )
1905  {
1906  /* x*y = lbx * y + (x-lbx) * y <= lbx * y + (x-lbx) * uby <= lbx * y + max{(ubx-lbx) * uby, 0 * uby} */
1907  coefx = 0.0;
1908  coefy = bilincoef * lbx;
1909  constant = bilincoef * (uby > 0.0 ? (ubx - lbx) * uby : 0.0);
1910  }
1911  else if( SCIPisRelEQ(scip, lby, uby) )
1912  {
1913  /* x*y = uby * x + (y-uby) * x <= uby * x + (y-uby) * lbx <= uby * x + max{(lby-uby) * lbx, 0 * lbx} */
1914  coefx = bilincoef * uby;
1915  coefy = 0.0;
1916  constant = bilincoef * (lbx < 0.0 ? (lby - uby) * lbx : 0.0);
1917  }
1918  else
1919  {
1920  coefx = bilincoef * uby;
1921  coefy = bilincoef * lbx;
1922  constant = -bilincoef * lbx * uby;
1923  }
1924  }
1925  else
1926  {
1927  *success = FALSE;
1928  return;
1929  }
1930  }
1931 
1932  if( SCIPisInfinity(scip, REALABS(coefx)) || SCIPisInfinity(scip, REALABS(coefy))
1933  || SCIPisInfinity(scip, REALABS(constant)) )
1934  {
1935  *success = FALSE;
1936  return;
1937  }
1938 
1939  if( overestimate )
1940  {
1941  coefx = -coefx;
1942  coefy = -coefy;
1943  constant = -constant;
1944  }
1945 
1946  SCIPdebugMsg(scip, "%.15g * x[%.15g,%.15g] * y[%.15g,%.15g] %c= %.15g * x %+.15g * y %+.15g\n", bilincoef, lbx, ubx,
1947  lby, uby, overestimate ? '<' : '>', coefx, coefy, constant);
1948 
1949  *lincoefx += coefx;
1950  *lincoefy += coefy;
1951  *linconstant += constant;
1952 }
1953 
1954 /** computes coefficients of linearization of a bilinear term in a reference point when given a linear inequality
1955  * involving only the variables of the bilinear term
1956  *
1957  * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
1958  * by Marco Locatelli
1959  */
1961  SCIP* scip, /**< SCIP data structure */
1962  SCIP_Real bilincoef, /**< coefficient of bilinear term */
1963  SCIP_Real lbx, /**< lower bound on first variable */
1964  SCIP_Real ubx, /**< upper bound on first variable */
1965  SCIP_Real refpointx, /**< reference point for first variable */
1966  SCIP_Real lby, /**< lower bound on second variable */
1967  SCIP_Real uby, /**< upper bound on second variable */
1968  SCIP_Real refpointy, /**< reference point for second variable */
1969  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
1970  SCIP_Real xcoef, /**< x coefficient of linear inequality; must be in {-1,0,1} */
1971  SCIP_Real ycoef, /**< y coefficient of linear inequality */
1972  SCIP_Real constant, /**< constant of linear inequality */
1973  SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
1974  SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
1975  SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
1976  SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
1977  )
1978 {
1979  SCIP_Real xs[2] = {lbx, ubx};
1980  SCIP_Real ys[2] = {lby, uby};
1981  SCIP_Real minx;
1982  SCIP_Real maxx;
1983  SCIP_Real miny;
1984  SCIP_Real maxy;
1985  SCIP_Real QUAD(lincoefyq);
1986  SCIP_Real QUAD(lincoefxq);
1987  SCIP_Real QUAD(linconstantq);
1988  SCIP_Real QUAD(denomq);
1989  SCIP_Real QUAD(mjq);
1990  SCIP_Real QUAD(qjq);
1991  SCIP_Real QUAD(xjq);
1992  SCIP_Real QUAD(yjq);
1993  SCIP_Real QUAD(tmpq);
1994  SCIP_Real vx;
1995  SCIP_Real vy;
1996  int n;
1997  int i;
1998 
1999  assert(scip != NULL);
2000  assert(!SCIPisInfinity(scip, lbx));
2001  assert(!SCIPisInfinity(scip, -ubx));
2002  assert(!SCIPisInfinity(scip, lby));
2003  assert(!SCIPisInfinity(scip, -uby));
2004  assert(SCIPisLE(scip, lbx, ubx));
2005  assert(SCIPisLE(scip, lby, uby));
2006  assert(SCIPisLE(scip, lbx, refpointx));
2007  assert(SCIPisGE(scip, ubx, refpointx));
2008  assert(SCIPisLE(scip, lby, refpointy));
2009  assert(SCIPisGE(scip, uby, refpointy));
2010  assert(lincoefx != NULL);
2011  assert(lincoefy != NULL);
2012  assert(linconstant != NULL);
2013  assert(success != NULL);
2014  assert(xcoef == 0.0 || xcoef == -1.0 || xcoef == 1.0); /*lint !e777*/
2015  assert(ycoef != SCIP_INVALID && ycoef != 0.0); /*lint !e777*/
2016  assert(constant != SCIP_INVALID); /*lint !e777*/
2017 
2018  *success = FALSE;
2019  *lincoefx = SCIP_INVALID;
2020  *lincoefy = SCIP_INVALID;
2021  *linconstant = SCIP_INVALID;
2022 
2023  /* reference point does not satisfy linear inequality */
2024  if( SCIPisFeasGT(scip, xcoef * refpointx - ycoef * refpointy - constant, 0.0) )
2025  return;
2026 
2027  /* compute minimal and maximal bounds on x and y for accepting the reference point */
2028  minx = lbx + 0.01 * (ubx-lbx);
2029  maxx = ubx - 0.01 * (ubx-lbx);
2030  miny = lby + 0.01 * (uby-lby);
2031  maxy = uby - 0.01 * (uby-lby);
2032 
2033  /* check whether the reference point is in [minx,maxx]x[miny,maxy] */
2034  if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2035  || SCIPisLE(scip, refpointy, miny) || SCIPisGE(scip, refpointy, maxy) )
2036  return;
2037 
2038  /* always consider xy without the bilinear coefficient */
2039  if( bilincoef < 0.0 )
2040  overestimate = !overestimate;
2041 
2042  /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2043  /* mj = xcoef / ycoef */
2044  SCIPquadprecDivDD(mjq, xcoef, ycoef);
2045 
2046  /* qj = -constant / ycoef */
2047  SCIPquadprecDivDD(qjq, -constant, ycoef);
2048 
2049  /* mj > 0 => underestimate; mj < 0 => overestimate */
2050  if( SCIPisNegative(scip, QUAD_TO_DBL(mjq)) != overestimate )
2051  return;
2052 
2053  /* get the corner point that satisfies the linear inequality xcoef*x <= ycoef*y + constant */
2054  if( !overestimate )
2055  {
2056  ys[0] = uby;
2057  ys[1] = lby;
2058  }
2059 
2060  vx = SCIP_INVALID;
2061  vy = SCIP_INVALID;
2062  n = 0;
2063  for( i = 0; i < 2; ++i )
2064  {
2065  SCIP_Real activity = xcoef * xs[i] - ycoef * ys[i] - constant;
2066  if( SCIPisLE(scip, activity, 0.0) )
2067  {
2068  /* corner point is satisfies inequality */
2069  vx = xs[i];
2070  vy = ys[i];
2071  }
2072  else if( SCIPisFeasGT(scip, activity, 0.0) )
2073  /* corner point is clearly cut off */
2074  ++n;
2075  }
2076 
2077  /* skip if no corner point satisfies the inequality or if no corner point is cut off
2078  * (that is, all corner points satisfy the inequality almost [1e-9..1e-6]) */
2079  if( n != 1 || vx == SCIP_INVALID || vy == SCIP_INVALID ) /*lint !e777*/
2080  return;
2081 
2082  /* denom = mj*(refpointx - vx) + vy - refpointy */
2083  SCIPquadprecSumDD(denomq, refpointx, -vx); /* refpoint - vx */
2084  SCIPquadprecProdQQ(denomq, denomq, mjq); /* mj * (refpoint - vx) */
2085  SCIPquadprecSumQD(denomq, denomq, vy); /* mj * (refpoint - vx) + vy */
2086  SCIPquadprecSumQD(denomq, denomq, -refpointy); /* mj * (refpoint - vx) + vy - refpointy */
2087 
2088  if( SCIPisZero(scip, QUAD_TO_DBL(denomq)) )
2089  return;
2090 
2091  /* (xj,yj) is the projection onto the line xcoef*x = ycoef*y + constant */
2092  /* xj = (refpointx*(vy - qj) - vx*(refpointy - qj)) / denom */
2093  SCIPquadprecProdQD(xjq, qjq, -1.0); /* - qj */
2094  SCIPquadprecSumQD(xjq, xjq, vy); /* vy - qj */
2095  SCIPquadprecProdQD(xjq, xjq, refpointx); /* refpointx * (vy - qj) */
2096  SCIPquadprecProdQD(tmpq, qjq, -1.0); /* - qj */
2097  SCIPquadprecSumQD(tmpq, tmpq, refpointy); /* refpointy - qj */
2098  SCIPquadprecProdQD(tmpq, tmpq, -vx); /* - vx * (refpointy - qj) */
2099  SCIPquadprecSumQQ(xjq, xjq, tmpq); /* refpointx * (vy - qj) - vx * (refpointy - qj) */
2100  SCIPquadprecDivQQ(xjq, xjq, denomq); /* (refpointx * (vy - qj) - vx * (refpointy - qj)) / denom */
2101 
2102  /* yj = mj * xj + qj */
2103  SCIPquadprecProdQQ(yjq, mjq, xjq);
2104  SCIPquadprecSumQQ(yjq, yjq, qjq);
2105 
2106  assert(SCIPisFeasEQ(scip, xcoef*QUAD_TO_DBL(xjq) - ycoef*QUAD_TO_DBL(yjq) - constant, 0.0));
2107 
2108  /* check whether the projection is in [minx,maxx] x [miny,maxy]; this avoids numerical difficulties when the
2109  * projection is close to the variable bounds
2110  */
2111  if( SCIPisLE(scip, QUAD_TO_DBL(xjq), minx) || SCIPisGE(scip, QUAD_TO_DBL(xjq), maxx)
2112  || SCIPisLE(scip, QUAD_TO_DBL(yjq), miny) || SCIPisGE(scip, QUAD_TO_DBL(yjq), maxy) )
2113  return;
2114 
2115  assert(vy - QUAD_TO_DBL(mjq)*vx - QUAD_TO_DBL(qjq) != 0.0);
2116 
2117  /* lincoefy = (mj*SQR(xj) - 2.0*mj*vx*xj - qj*vx + vx*vy) / (vy - mj*vx - qj) */
2118  SCIPquadprecSquareQ(lincoefyq, xjq); /* xj^2 */
2119  SCIPquadprecProdQQ(lincoefyq, lincoefyq, mjq); /* mj * xj^2 */
2120  SCIPquadprecProdQQ(tmpq, mjq, xjq); /* mj * xj */
2121  SCIPquadprecProdQD(tmpq, tmpq, -2.0 * vx); /* -2 * vx * mj * xj */
2122  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj */
2123  SCIPquadprecProdQD(tmpq, qjq, -vx); /* -qj * vx */
2124  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx */
2125  SCIPquadprecProdDD(tmpq, vx, vy); /* vx * vy */
2126  SCIPquadprecSumQQ(lincoefyq, lincoefyq, tmpq); /* mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy */
2127  SCIPquadprecProdQD(tmpq, mjq, vx); /* mj * vx */
2128  SCIPquadprecSumQD(tmpq, tmpq, -vy); /* -vy + mj * vx */
2129  SCIPquadprecSumQQ(tmpq, tmpq, qjq); /* -vy + mj * vx + qj */
2130  QUAD_SCALE(tmpq, -1.0); /* vy - mj * vx - qj */
2131  SCIPquadprecDivQQ(lincoefyq, lincoefyq, tmpq); /* (mj * xj^2 -2 * vx * mj * xj -qj * vx + vx * vy) / (vy - mj * vx - qj) */
2132 
2133  /* lincoefx = 2.0*mj*xj + qj - mj*(*lincoefy) */
2134  SCIPquadprecProdQQ(lincoefxq, mjq, xjq); /* mj * xj */
2135  QUAD_SCALE(lincoefxq, 2.0); /* 2 * mj * xj */
2136  SCIPquadprecSumQQ(lincoefxq, lincoefxq, qjq); /* 2 * mj * xj + qj */
2137  SCIPquadprecProdQQ(tmpq, mjq, lincoefyq); /* mj * lincoefy */
2138  QUAD_SCALE(tmpq, -1.0); /* - mj * lincoefy */
2139  SCIPquadprecSumQQ(lincoefxq, lincoefxq, tmpq); /* 2 * mj * xj + qj - mj * lincoefy */
2140 
2141  /* linconstant = -mj*SQR(xj) - (*lincoefy)*qj */
2142  SCIPquadprecSquareQ(linconstantq, xjq); /* xj^2 */
2143  SCIPquadprecProdQQ(linconstantq, linconstantq, mjq); /* mj * xj^2 */
2144  QUAD_SCALE(linconstantq, -1.0); /* - mj * xj^2 */
2145  SCIPquadprecProdQQ(tmpq, lincoefyq, qjq); /* lincoefy * qj */
2146  QUAD_SCALE(tmpq, -1.0); /* - lincoefy * qj */
2147  SCIPquadprecSumQQ(linconstantq, linconstantq, tmpq); /* - mj * xj^2 - lincoefy * qj */
2148 
2149  /* consider the bilinear coefficient */
2150  SCIPquadprecProdQD(lincoefxq, lincoefxq, bilincoef);
2151  SCIPquadprecProdQD(lincoefyq, lincoefyq, bilincoef);
2152  SCIPquadprecProdQD(linconstantq, linconstantq, bilincoef);
2153  *lincoefx = QUAD_TO_DBL(lincoefxq);
2154  *lincoefy = QUAD_TO_DBL(lincoefyq);
2155  *linconstant = QUAD_TO_DBL(linconstantq);
2156 
2157  /* cut needs to be tight at (vx,vy) and (xj,yj); otherwise we consider the cut to be numerically bad */
2158  *success = SCIPisFeasEQ(scip, (*lincoefx)*vx + (*lincoefy)*vy + (*linconstant), bilincoef*vx*vy)
2159  && SCIPisFeasEQ(scip, (*lincoefx)*QUAD_TO_DBL(xjq) + (*lincoefy)*QUAD_TO_DBL(yjq) + (*linconstant),
2160  bilincoef*QUAD_TO_DBL(xjq)*QUAD_TO_DBL(yjq));
2161 
2162 #ifndef NDEBUG
2163  {
2164  SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2165 
2166  /* cut needs to under- or overestimate the bilinear term at the reference point */
2167  if( bilincoef < 0.0 )
2168  overestimate = !overestimate;
2169 
2170  if( overestimate )
2171  assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2172  else
2173  assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2174  }
2175 #endif
2176 }
2177 
2178 /** computes coefficients of linearization of a bilinear term in a reference point when given two linear inequalities
2179  * involving only the variables of the bilinear term
2180  *
2181  * @note the formulas are extracted from "Convex envelopes of bivariate functions through the solution of KKT systems"
2182  * by Marco Locatelli
2183  */
2185  SCIP* scip, /**< SCIP data structure */
2186  SCIP_Real bilincoef, /**< coefficient of bilinear term */
2187  SCIP_Real lbx, /**< lower bound on first variable */
2188  SCIP_Real ubx, /**< upper bound on first variable */
2189  SCIP_Real refpointx, /**< reference point for first variable */
2190  SCIP_Real lby, /**< lower bound on second variable */
2191  SCIP_Real uby, /**< upper bound on second variable */
2192  SCIP_Real refpointy, /**< reference point for second variable */
2193  SCIP_Bool overestimate, /**< whether to compute an overestimator instead of an underestimator */
2194  SCIP_Real xcoef1, /**< x coefficient of linear inequality; must be in {-1,0,1} */
2195  SCIP_Real ycoef1, /**< y coefficient of linear inequality */
2196  SCIP_Real constant1, /**< constant of linear inequality */
2197  SCIP_Real xcoef2, /**< x coefficient of linear inequality; must be in {-1,0,1} */
2198  SCIP_Real ycoef2, /**< y coefficient of linear inequality */
2199  SCIP_Real constant2, /**< constant of linear inequality */
2200  SCIP_Real* RESTRICT lincoefx, /**< buffer to store coefficient of first variable in linearization */
2201  SCIP_Real* RESTRICT lincoefy, /**< buffer to store coefficient of second variable in linearization */
2202  SCIP_Real* RESTRICT linconstant, /**< buffer to store constant of linearization */
2203  SCIP_Bool* RESTRICT success /**< buffer to store whether linearization was successful */
2204  )
2205 {
2206  SCIP_Real mi, mj, qi, qj, xi, xj, yi, yj;
2207  SCIP_Real xcoef, ycoef, constant;
2208  SCIP_Real minx, maxx, miny, maxy;
2209 
2210  assert(scip != NULL);
2211  assert(!SCIPisInfinity(scip, lbx));
2212  assert(!SCIPisInfinity(scip, -ubx));
2213  assert(!SCIPisInfinity(scip, lby));
2214  assert(!SCIPisInfinity(scip, -uby));
2215  assert(SCIPisLE(scip, lbx, ubx));
2216  assert(SCIPisLE(scip, lby, uby));
2217  assert(SCIPisLE(scip, lbx, refpointx));
2218  assert(SCIPisGE(scip, ubx, refpointx));
2219  assert(SCIPisLE(scip, lby, refpointy));
2220  assert(SCIPisGE(scip, uby, refpointy));
2221  assert(lincoefx != NULL);
2222  assert(lincoefy != NULL);
2223  assert(linconstant != NULL);
2224  assert(success != NULL);
2225  assert(xcoef1 != 0.0 && xcoef1 != SCIP_INVALID); /*lint !e777*/
2226  assert(ycoef1 != SCIP_INVALID && ycoef1 != 0.0); /*lint !e777*/
2227  assert(constant1 != SCIP_INVALID); /*lint !e777*/
2228  assert(xcoef2 != 0.0 && xcoef2 != SCIP_INVALID); /*lint !e777*/
2229  assert(ycoef2 != SCIP_INVALID && ycoef2 != 0.0); /*lint !e777*/
2230  assert(constant2 != SCIP_INVALID); /*lint !e777*/
2231 
2232  *success = FALSE;
2233  *lincoefx = SCIP_INVALID;
2234  *lincoefy = SCIP_INVALID;
2235  *linconstant = SCIP_INVALID;
2236 
2237  /* reference point does not satisfy linear inequalities */
2238  if( SCIPisFeasGT(scip, xcoef1 * refpointx - ycoef1 * refpointy - constant1, 0.0)
2239  || SCIPisFeasGT(scip, xcoef2 * refpointx - ycoef2 * refpointy - constant2, 0.0) )
2240  return;
2241 
2242  /* compute minimal and maximal bounds on x and y for accepting the reference point */
2243  minx = lbx + 0.01 * (ubx-lbx);
2244  maxx = ubx - 0.01 * (ubx-lbx);
2245  miny = lby + 0.01 * (uby-lby);
2246  maxy = uby - 0.01 * (uby-lby);
2247 
2248  /* check the reference point is in the interior of the domain */
2249  if( SCIPisLE(scip, refpointx, minx) || SCIPisGE(scip, refpointx, maxx)
2250  || SCIPisLE(scip, refpointy, miny) || SCIPisFeasGE(scip, refpointy, maxy) )
2251  return;
2252 
2253  /* the sign of the x-coefficients of the two inequalities must be different; otherwise the convex or concave
2254  * envelope can be computed via SCIPcomputeBilinEnvelope1 for each inequality separately
2255  */
2256  if( (xcoef1 > 0) == (xcoef2 > 0) )
2257  return;
2258 
2259  /* always consider xy without the bilinear coefficient */
2260  if( bilincoef < 0.0 )
2261  overestimate = !overestimate;
2262 
2263  /* we use same notation as in "Convex envelopes of bivariate functions through the solution of KKT systems", 2016 */
2264  mi = xcoef1 / ycoef1;
2265  qi = -constant1 / ycoef1;
2266  mj = xcoef2 / ycoef2;
2267  qj = -constant2 / ycoef2;
2268 
2269  /* mi, mj > 0 => underestimate; mi, mj < 0 => overestimate */
2270  if( SCIPisNegative(scip, mi) != overestimate || SCIPisNegative(scip, mj) != overestimate )
2271  return;
2272 
2273  /* compute cut according to Locatelli 2016 */
2274  computeBilinEnvelope2(scip, refpointx, refpointy, mi, qi, mj, qj, &xi, &yi, &xj, &yj, &xcoef, &ycoef, &constant);
2275  assert(SCIPisRelEQ(scip, mi*xi + qi, yi));
2276  assert(SCIPisRelEQ(scip, mj*xj + qj, yj));
2277 
2278  /* it might happen that (xi,yi) = (xj,yj) if the two lines intersect */
2279  if( SCIPisEQ(scip, xi, xj) && SCIPisEQ(scip, yi, yj) )
2280  return;
2281 
2282  /* check whether projected points are in the interior */
2283  if( SCIPisLE(scip, xi, minx) || SCIPisGE(scip, xi, maxx) || SCIPisLE(scip, yi, miny) || SCIPisGE(scip, yi, maxy) )
2284  return;
2285  if( SCIPisLE(scip, xj, minx) || SCIPisGE(scip, xj, maxx) || SCIPisLE(scip, yj, miny) || SCIPisGE(scip, yj, maxy) )
2286  return;
2287 
2288  *lincoefx = bilincoef * xcoef;
2289  *lincoefy = bilincoef * ycoef;
2290  *linconstant = bilincoef * constant;
2291 
2292  /* cut needs to be tight at (vx,vy) and (xj,yj) */
2293  *success = SCIPisFeasEQ(scip, (*lincoefx)*xi + (*lincoefy)*yi + (*linconstant), bilincoef*xi*yi)
2294  && SCIPisFeasEQ(scip, (*lincoefx)*xj + (*lincoefy)*yj + (*linconstant), bilincoef*xj*yj);
2295 
2296 #ifndef NDEBUG
2297  {
2298  SCIP_Real activity = (*lincoefx)*refpointx + (*lincoefy)*refpointy + (*linconstant);
2299 
2300  /* cut needs to under- or overestimate the bilinear term at the reference point */
2301  if( bilincoef < 0.0 )
2302  overestimate = !overestimate;
2303 
2304  if( overestimate )
2305  assert(SCIPisFeasGE(scip, activity, bilincoef*refpointx*refpointy));
2306  else
2307  assert(SCIPisFeasLE(scip, activity, bilincoef*refpointx*refpointy));
2308  }
2309 #endif
2310 }
void SCIPintervalSetEntire(SCIP_Real infinity, SCIP_INTERVAL *resultant)
SCIP_Bool SCIPisFeasZero(SCIP *scip, SCIP_Real val)
#define SCIPreallocBlockMemoryArray(scip, ptr, oldnum, newnum)
Definition: scip_mem.h:90
static void computeBilinEnvelope2(SCIP *scip, SCIP_Real x, SCIP_Real y, SCIP_Real mi, SCIP_Real qi, SCIP_Real mj, SCIP_Real qj, SCIP_Real *RESTRICT xi, SCIP_Real *RESTRICT yi, SCIP_Real *RESTRICT xj, SCIP_Real *RESTRICT yj, SCIP_Real *RESTRICT xcoef, SCIP_Real *RESTRICT ycoef, SCIP_Real *RESTRICT constant)
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:491
SCIP_Real SCIPfeastol(SCIP *scip)
SCIP_RETCODE SCIPincludeTable(SCIP *scip, const char *name, const char *desc, SCIP_Bool active, SCIP_DECL_TABLECOPY((*tablecopy)), SCIP_DECL_TABLEFREE((*tablefree)), SCIP_DECL_TABLEINIT((*tableinit)), SCIP_DECL_TABLEEXIT((*tableexit)), SCIP_DECL_TABLEINITSOL((*tableinitsol)), SCIP_DECL_TABLEEXITSOL((*tableexitsol)), SCIP_DECL_TABLEOUTPUT((*tableoutput)), SCIP_TABLEDATA *tabledata, int position, SCIP_STAGE earlieststage)
Definition: scip_table.c:47
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_NODE * SCIPgetCurrentNode(SCIP *scip)
Definition: scip_tree.c:82
SCIP_STAGE SCIPgetStage(SCIP *scip)
Definition: scip_general.c:356
#define nlhdlrExitSepaBilinear
static SCIP_DECL_NLHDLRDETECT(nlhdlrDetectBilinear)
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition: scip_expr.c:1476
SCIP_RETCODE SCIPhashmapSetImageInt(SCIP_HASHMAP *hashmap, void *origin, int image)
Definition: misc.c:3297
SCIP_Bool SCIPisRelEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPquadprecSumQD(r, a, b)
Definition: dbldblarith.h:53
unsigned int SCIPgetExprNAuxvarUsesNonlinear(SCIP_EXPR *expr)
static void getIneqViol(SCIP_VAR *x, SCIP_VAR *y, SCIP_Real xcoef, SCIP_Real ycoef, SCIP_Real constant, SCIP_Real *viol1, SCIP_Real *viol2)
SCIP_CONSHDLR * SCIPfindConshdlr(SCIP *scip, const char *name)
Definition: scip_cons.c:877
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3798
void SCIPcomputeBilinEnvelope2(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real xcoef1, SCIP_Real ycoef1, SCIP_Real constant1, SCIP_Real xcoef2, SCIP_Real ycoef2, SCIP_Real constant2, SCIP_Real *RESTRICT lincoefx, SCIP_Real *RESTRICT lincoefy, SCIP_Real *RESTRICT linconstant, SCIP_Bool *RESTRICT success)
void SCIPnlhdlrSetFreeExprData(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRFREEEXPRDATA((*freeexprdata)))
Definition: nlhdlr.c:85
#define MIN_INTERIORITY
int SCIPcalcMemGrowSize(SCIP *scip, int num)
Definition: scip_mem.c:130
SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:17966
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPvarIsBinary(SCIP_VAR *var)
Definition: var.c:17431
SCIP_TABLE * SCIPfindTable(SCIP *scip, const char *name)
Definition: scip_table.c:85
SCIP_RETCODE SCIPregisterExprUsageNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Bool useauxvar, SCIP_Bool useactivityforprop, SCIP_Bool useactivityforsepabelow, SCIP_Bool useactivityforsepaabove)
#define SCIPquadprecSquareQ(r, a)
Definition: dbldblarith.h:59
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPquadprecProdQQ(r, a, b)
Definition: dbldblarith.h:57
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4547
#define FALSE
Definition: def.h:87
SCIP_RETCODE SCIPhashmapCreate(SCIP_HASHMAP **hashmap, BMS_BLKMEM *blkmem, int mapsize)
Definition: misc.c:3014
SCIP_RETCODE SCIPmarkExprPropagateNonlinear(SCIP *scip, SCIP_EXPR *expr)
#define EPSEQ(x, y, eps)
Definition: def.h:202
void SCIPnlhdlrSetFreeHdlrData(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRFREEHDLRDATA((*freehdlrdata)))
Definition: nlhdlr.c:74
SCIP_NLHDLREXPRDATA * SCIPgetNlhdlrExprDataNonlinear(SCIP_NLHDLR *nlhdlr, SCIP_EXPR *expr)
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
void SCIPintervalSolveUnivariateQuadExpression(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL sqrcoeff, SCIP_INTERVAL lincoeff, SCIP_INTERVAL rhs, SCIP_INTERVAL xbnds)
#define TRUE
Definition: def.h:86
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_NLHDLR * SCIPfindNlhdlrNonlinear(SCIP_CONSHDLR *conshdlr, const char *name)
SCIP_RETCODE SCIPhashmapInsertInt(SCIP_HASHMAP *hashmap, void *origin, int image)
Definition: misc.c:3132
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition: expr.c:3954
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
#define NLHDLR_NAME
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:99
static SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateBilinear)
#define QUAD_SCALE(x, a)
Definition: dbldblarith.h:41
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void SCIPcaptureExpr(SCIP_EXPR *expr)
Definition: scip_expr.c:1399
static SCIP_Bool useBilinIneqs(SCIP *scip, SCIP_VAR *x, SCIP_VAR *y, SCIP_Real refx, SCIP_Real refy)
variable expression handler
SCIP_RETCODE SCIPensureRowprepSize(SCIP *scip, SCIP_ROWPREP *rowprep, int size)
Definition: misc_rowprep.c:855
#define nlhdlrInitSepaBilinear
void SCIPwarningMessage(SCIP *scip, const char *formatstr,...)
Definition: scip_message.c:111
#define SCIPdebugMsg
Definition: scip_message.h:69
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:74
#define SCIPquadprecDivQD(r, a, b)
Definition: dbldblarith.h:56
SCIP_EXPR * SCIPexpriterGetCurrent(SCIP_EXPRITER *iterator)
Definition: expriter.c:673
SCIP_VAR ** x
Definition: circlepacking.c:54
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_Real SCIPepsilon(SCIP *scip)
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1454
#define QUAD_TO_DBL(x)
Definition: dbldblarith.h:40
bilinear nonlinear handler
SCIP_Bool SCIPisRelGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisRelLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPhashmapExists(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3363
#define NLHDLR_DETECTPRIORITY
static SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrFreehdlrdataBilinear)
SCIP_Longint SCIPnodeGetNumber(SCIP_NODE *node)
Definition: tree.c:7432
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3808
SCIP_Real inf
Definition: intervalarith.h:46
SCIP_Bool SCIPisRelGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
const char * SCIPnlhdlrGetName(SCIP_NLHDLR *nlhdlr)
Definition: nlhdlr.c:141
SCIP_VAR * SCIPgetVarExprVar(SCIP_EXPR *expr)
Definition: expr_var.c:403
static SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrBilinear)
static SCIP_DECL_NLHDLREXIT(nlhdlrExitBilinear)
SCIP_RETCODE SCIPincludeNlhdlrNonlinear(SCIP *scip, SCIP_NLHDLR **nlhdlr, const char *name, const char *desc, int detectpriority, int enfopriority, SCIP_DECL_NLHDLRDETECT((*detect)), SCIP_DECL_NLHDLREVALAUX((*evalaux)), SCIP_NLHDLRDATA *nlhdlrdata)
BMS_BLKMEM * SCIPblkmem(SCIP *scip)
Definition: scip_mem.c:48
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
#define QUAD(x)
Definition: dbldblarith.h:38
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17251
void SCIPhashmapFree(SCIP_HASHMAP **hashmap)
Definition: misc.c:3048
SCIP_EXPR * SCIPgetExprNonlinear(SCIP_CONS *cons)
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalMulScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define REALABS(x)
Definition: def.h:201
void SCIPintervalIntersect(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
#define TABLE_DESC_BILINEAR
#define NLHDLR_DESC
#define SCIP_CALL(x)
Definition: def.h:384
SCIP_Bool SCIPisFeasGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxBilinear)
SCIP_Real sup
Definition: intervalarith.h:47
#define SCIPfreeBlockMemoryNull(scip, ptr)
Definition: scip_mem.h:100
SCIP_Bool SCIPisFeasLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define SCIPquadprecProdDD(r, a, b)
Definition: dbldblarith.h:49
#define SCIP_INTERVAL_INFINITY
Definition: def.h:199
static void reversePropBilinear(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_EXPR *expr, SCIP_INTERVAL exprbounds, SCIP_Real *underineqs, int nunderineqs, SCIP_Real *overineqs, int noverineqs, SCIP_INTERVAL *intervalx, SCIP_INTERVAL *intervaly)
int SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4590
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2300
SCIP_Real SCIPgetCoefExprProduct(SCIP_EXPR *expr)
static void updateBilinearRelaxation(SCIP *scip, SCIP_VAR *RESTRICT x, SCIP_VAR *RESTRICT y, SCIP_Real bilincoef, SCIP_SIDETYPE violside, SCIP_Real refx, SCIP_Real refy, SCIP_Real *RESTRICT ineqs, int nineqs, SCIP_Real mccormickval, SCIP_Real *RESTRICT bestcoefx, SCIP_Real *RESTRICT bestcoefy, SCIP_Real *RESTRICT bestconst, SCIP_Real *RESTRICT bestval, SCIP_Bool *success)
#define SCIPquadprecSqrtQ(r, a)
Definition: dbldblarith.h:62
void SCIPnlhdlrSetSepa(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINITSEPA((*initsepa)), SCIP_DECL_NLHDLRENFO((*enfo)), SCIP_DECL_NLHDLRESTIMATE((*estimate)), SCIP_DECL_NLHDLREXITSEPA((*exitsepa)))
Definition: nlhdlr.c:123
#define SCIP_Bool
Definition: def.h:84
#define TABLE_EARLIEST_STAGE_BILINEAR
static SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropBilinear)
void SCIProwprepAddConstant(SCIP_ROWPREP *rowprep, SCIP_Real constant)
Definition: misc_rowprep.c:728
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:661
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1407
constraint handler for nonlinear constraints specified by algebraic expressions
#define MAX(x, y)
Definition: tclique_def.h:83
#define TABLE_NAME_BILINEAR
SCIP_EXPR * SCIPexpriterGetNext(SCIP_EXPRITER *iterator)
Definition: expriter.c:848
void SCIPnlhdlrSetProp(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINTEVAL((*inteval)), SCIP_DECL_NLHDLRREVERSEPROP((*reverseprop)))
Definition: nlhdlr.c:110
#define SCIPquadprecProdQD(r, a, b)
Definition: dbldblarith.h:54
void SCIPaddBilinMcCormick(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real *lincoefx, SCIP_Real *lincoefy, SCIP_Real *linconstant, SCIP_Bool *success)
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
#define BMSclearMemory(ptr)
Definition: memory.h:122
#define SCIPquadprecSumQQ(r, a, b)
Definition: dbldblarith.h:58
SCIP_RETCODE SCIPtightenExprIntervalNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_INTERVAL newbounds, SCIP_Bool *cutoff, int *ntightenings)
SCIP_EXPR ** SCIPgetExprsBilinear(SCIP_NLHDLR *nlhdlr)
int SCIPgetNVars(SCIP *scip)
Definition: scip_prob.c:1990
void SCIPfreeExpriter(SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2314
product expression handler
SCIP_RETCODE SCIPaddIneqBilinear(SCIP *scip, SCIP_NLHDLR *nlhdlr, SCIP_EXPR *expr, SCIP_Real xcoef, SCIP_Real ycoef, SCIP_Real constant, SCIP_Bool *success)
static SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalBilinear)
#define SCIPquadprecDivDD(r, a, b)
Definition: dbldblarith.h:52
SCIP_RETCODE SCIPhashmapRemove(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3379
int SCIPhashmapGetNElements(SCIP_HASHMAP *hashmap)
Definition: misc.c:3473
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_Bool isPointFeasible(SCIP *scip, SCIP_Real x, SCIP_Real y, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real lby, SCIP_Real uby, SCIP_Real *ineqs, int nineqs)
static void getFeasiblePointsBilinear(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_EXPR *expr, SCIP_INTERVAL exprbounds, SCIP_Real *underineqs, int nunderineqs, SCIP_Real *overineqs, int noverineqs, SCIP_Bool levelset, SCIP_Real *xs, SCIP_Real *ys, int *npoints)
#define TABLE_POSITION_BILINEAR
#define SCIP_NLHDLR_METHOD_ACTIVITY
Definition: type_nlhdlr.h:45
#define nlhdlrEnfoBilinear
static SCIP_INTERVAL intevalBilinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Real *underineqs, int nunderineqs, SCIP_Real *overineqs, int noverineqs)
void SCIPaddBilinLinearization(SCIP *scip, SCIP_Real bilincoef, SCIP_Real refpointx, SCIP_Real refpointy, SCIP_Real *lincoefx, SCIP_Real *lincoefy, SCIP_Real *linconstant, SCIP_Bool *success)
SCIP_Bool SCIPisExprVar(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1421
SCIP_RETCODE SCIPincludeNlhdlrBilinear(SCIP *scip)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPnlhdlrSetInitExit(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINIT((*init)), SCIP_DECL_NLHDLREXIT((*exit_)))
Definition: nlhdlr.c:97
int SCIPgetNExprsBilinear(SCIP_NLHDLR *nlhdlr)
#define SCIP_Real
Definition: def.h:177
SCIP_VAR ** y
Definition: circlepacking.c:55
#define SCIP_INVALID
Definition: def.h:197
#define SCIPquadprecSumDD(r, a, b)
Definition: dbldblarith.h:51
SCIP_NLHDLRDATA * SCIPnlhdlrGetData(SCIP_NLHDLR *nlhdlr)
Definition: nlhdlr.c:191
#define SCIP_Longint
Definition: def.h:162
struct SCIP_NlhdlrExprData SCIP_NLHDLREXPRDATA
Definition: type_nlhdlr.h:404
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
Definition: misc_rowprep.c:538
struct SCIP_NlhdlrData SCIP_NLHDLRDATA
Definition: type_nlhdlr.h:403
#define MIN_ABSBOUNDSIZE
void SCIPcomputeBilinEnvelope1(SCIP *scip, SCIP_Real bilincoef, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real refpointx, SCIP_Real lby, SCIP_Real uby, SCIP_Real refpointy, SCIP_Bool overestimate, SCIP_Real xcoef, SCIP_Real ycoef, SCIP_Real constant, SCIP_Real *RESTRICT lincoefx, SCIP_Real *RESTRICT lincoefy, SCIP_Real *RESTRICT linconstant, SCIP_Bool *RESTRICT success)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17976
SCIP_RETCODE SCIPaddRowprepTerm(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR *var, SCIP_Real coef)
Definition: misc_rowprep.c:881
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:102
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
#define SCIPallocClearBlockMemory(scip, ptr)
Definition: scip_mem.h:82
SCIP_Bool SCIPexpriterIsEnd(SCIP_EXPRITER *iterator)
Definition: expriter.c:959
SCIPallocBlockMemory(scip, subsol))
#define NLHDLR_ENFOPRIORITY
int SCIPhashmapGetImageInt(SCIP_HASHMAP *hashmap, void *origin)
Definition: misc.c:3221
static SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataBilinear)
SCIP_INTERVAL SCIPgetExprBoundsNonlinear(SCIP *scip, SCIP_EXPR *expr)
#define SCIP_NLHDLR_METHOD_SEPABOTH
Definition: type_nlhdlr.h:44
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1352
#define SCIPquadprecDivQQ(r, a, b)
Definition: dbldblarith.h:60
static SCIP_DECL_TABLEOUTPUT(tableOutputBilinear)
SCIP_Bool SCIPisRelLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPsetPtrarrayVal(SCIP *scip, SCIP_PTRARRAY *ptrarray, int idx, void *val)
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:48
#define nlhdlrInitBilinear
void SCIPnlhdlrSetCopyHdlr(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRCOPYHDLR((*copy)))
Definition: nlhdlr.c:63
enum SCIP_SideType SCIP_SIDETYPE
Definition: type_lp.h:58