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