Scippy

SCIP

Solving Constraint Integer Programs

nlhdlr_quotient.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file nlhdlr_quotient.c
17  * @ingroup DEFPLUGINS_NLHDLR
18  * @brief quotient nonlinear handler
19  * @author Benjamin Mueller
20  * @author Fabian Wegscheider
21  *
22  * @todo implement INITSEPA
23  * @todo use the convex envelope for x/y described in Tawarmalani and Sahinidis (2002) if y has a finite upper bound
24  */
25 
26 #include <string.h>
27 
28 #include "scip/nlhdlr_quotient.h"
29 #include "scip/cons_nonlinear.h"
30 #include "scip/pub_misc_rowprep.h"
31 #include "scip/nlhdlr.h"
32 #include "scip/nlhdlr_bilinear.h"
33 
34 /* fundamental nonlinear handler properties */
35 #define NLHDLR_NAME "quotient"
36 #define NLHDLR_DESC "nonlinear handler for quotient expressions"
37 #define NLHDLR_DETECTPRIORITY 20
38 #define NLHDLR_ENFOPRIORITY 20
39 
40 /** translate from one value of infinity to another
41  *
42  * if val is &ge; infty1, then give infty2, else give val
43  */
44 #define infty2infty(infty1, infty2, val) ((val) >= (infty1) ? (infty2) : (val))
45 
46 /*lint -e666*/
47 
48 /*
49  * Data structures
50  */
51 
52 /** nonlinear handler expression data */
53 struct SCIP_NlhdlrExprData
54 {
55  SCIP_EXPR* numexpr; /**< expression of the numerator */
56  SCIP_Real numcoef; /**< coefficient of the numerator */
57  SCIP_Real numconst; /**< constant of the numerator */
58  SCIP_EXPR* denomexpr; /**< expression of the denominator */
59  SCIP_Real denomcoef; /**< coefficient of the denominator */
60  SCIP_Real denomconst; /**< constant of the denominator */
61  SCIP_Real constant; /**< constant */
62 };
63 
64 /*
65  * Local methods
66  */
67 
68 /** helper method to create nonlinear handler expression data */
69 static
71  SCIP* scip, /**< SCIP data structure */
72  SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< nonlinear handler expression data */
73  SCIP_EXPR* numexpr, /**< expression of the numerator */
74  SCIP_Real numcoef, /**< coefficient of the numerator */
75  SCIP_Real numconst, /**< constant of the numerator */
76  SCIP_EXPR* denomexpr, /**< expression of the denominator */
77  SCIP_Real denomcoef, /**< coefficient of the denominator */
78  SCIP_Real denomconst, /**< constant of the denominator */
79  SCIP_Real constant /**< constant */
80  )
81 {
82  assert(nlhdlrexprdata != NULL);
83  assert(numexpr != NULL);
84  assert(denomexpr != NULL);
85  assert(!SCIPisZero(scip, numcoef));
86  assert(!SCIPisZero(scip, denomcoef));
87 
88  /* allocate memory */
89  SCIP_CALL( SCIPallocBlockMemory(scip, nlhdlrexprdata) );
90 
91  /* store values */
92  (*nlhdlrexprdata)->numexpr = numexpr;
93  (*nlhdlrexprdata)->numcoef = numcoef;
94  (*nlhdlrexprdata)->numconst = numconst;
95  (*nlhdlrexprdata)->denomexpr = denomexpr;
96  (*nlhdlrexprdata)->denomcoef = denomcoef;
97  (*nlhdlrexprdata)->denomconst = denomconst;
98  (*nlhdlrexprdata)->constant = constant;
99 
100  /* capture expressions */
101  SCIPcaptureExpr(numexpr);
102  SCIPcaptureExpr(denomexpr);
103 
104  return SCIP_OKAY;
105 }
106 
107 /** helper method to free nonlinear handler expression data */
108 static
110  SCIP* scip, /**< SCIP data structure */
111  SCIP_NLHDLREXPRDATA** nlhdlrexprdata /**< nonlinear handler expression data */
112  )
113 {
114  assert(nlhdlrexprdata != NULL);
115  assert(*nlhdlrexprdata != NULL);
116  assert((*nlhdlrexprdata)->numexpr != NULL);
117  assert((*nlhdlrexprdata)->denomexpr != NULL);
118 
119  /* release expressions */
120  SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->denomexpr) );
121  SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->numexpr) );
122 
123  /* free expression data of nonlinear handler */
124  SCIPfreeBlockMemory(scip, nlhdlrexprdata);
125 
126  return SCIP_OKAY;
127 }
128 
129 /** helper method to transform an expression g(x) into a*f(x) + b */
130 static
132  SCIP* scip, /**< SCIP data structure */
133  SCIP_EXPR* expr, /**< expression */
134  SCIP_EXPR** target, /**< pointer to store the expression f(x) */
135  SCIP_Real* coef, /**< pointer to store the coefficient */
136  SCIP_Real* constant /**< pointer to store the constant */
137  )
138 {
139  assert(expr != NULL);
140  assert(target != NULL);
141  assert(coef != NULL);
142  assert(constant != NULL);
143 
144  /* expression is a sum with one child */
145  if( SCIPisExprSum(scip, expr) && SCIPexprGetNChildren(expr) == 1 )
146  {
147  *target = SCIPexprGetChildren(expr)[0];
148  *coef = SCIPgetCoefsExprSum(expr)[0];
149  *constant = SCIPgetConstantExprSum(expr);
150  }
151  else /* otherwise return 1 * f(x) + 0 */
152  {
153  *target = expr;
154  *coef = 1.0;
155  *constant = 0.0;
156  }
157 }
158 
159 /** helper method to detect an expression of the form (a*x + b) / (c*y + d) + e
160  *
161  * Due to the expansion of products, there are two types of expressions that can be detected:
162  *
163  * 1. prod(f(x), pow(g(y),-1))
164  * 2. sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1))
165  *
166  * @todo At the moment quotients like xy / z are not detected, because they are turned into a product expression
167  * with three children, i.e., x * y * (1 / z).
168  */
169 static
171  SCIP* scip, /**< SCIP data structure */
172  SCIP_EXPR* expr, /**< expression */
173  SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nonlinear handler expression data */
174  SCIP_Bool* success /**< pointer to store whether nonlinear handler should be called for this expression */
175  )
176 {
177  SCIP_EXPR** children;
178  SCIP_EXPR* denomexpr = NULL;
179  SCIP_EXPR* numexpr = NULL;
180  SCIP_EXPR* xexpr = NULL;
181  SCIP_EXPR* yexpr = NULL;
182  SCIP_Real a, b, c, d, e;
183  SCIP_Real nomfac = 1.0;
184  SCIP_Real numconst = 0.0;
185 
186  assert(scip != NULL);
187  assert(expr != NULL);
188 
189  *success = FALSE;
190  a = 0.0;
191  b = 0.0;
192  c = 0.0;
193  d = 0.0;
194  e = 0.0;
195 
196  /* possible structures only have two children */
197  if( SCIPexprGetNChildren(expr) != 2 )
198  return SCIP_OKAY;
199 
200  /* expression must be either a product or a sum */
201  if( !SCIPisExprProduct(scip, expr) && !SCIPisExprSum(scip, expr) )
202  return SCIP_OKAY;
203 
204  children = SCIPexprGetChildren(expr);
205  assert(children != NULL);
206 
207  /* case: prod(f(x), pow(g(y),-1)) */
208  if( SCIPisExprProduct(scip, expr) )
209  {
210  if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0 ) /*lint !e777*/
211  {
212  denomexpr = SCIPexprGetChildren(children[0])[0];
213  numexpr = children[1];
214  }
215  else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0 ) /*lint !e777*/
216  {
217  denomexpr = SCIPexprGetChildren(children[1])[0];
218  numexpr = children[0];
219  }
220 
221  /* remember to scale the numerator by the coefficient stored in the product expression */
222  nomfac = SCIPgetCoefExprProduct(expr);
223  }
224  /* case: sum(prod(f(x),pow(g(y),-1)), pow(g(y),-1)) */
225  else
226  {
227  SCIP_Real* sumcoefs;
228 
229  assert(SCIPisExprSum(scip, expr));
230  sumcoefs = SCIPgetCoefsExprSum(expr);
231 
232  /* children[0] is 1/g(y) and children[1] is a product of f(x) and 1/g(y) */
233  if( SCIPisExprPower(scip, children[0]) && SCIPgetExponentExprPow(children[0]) == -1.0
234  && SCIPisExprProduct(scip, children[1]) && SCIPexprGetNChildren(children[1]) == 2 ) /* lint !e777 */
235  {
236  SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[1]);
237 
238  if( children[0] == SCIPexprGetChildren(children[1])[0] )
239  {
240  denomexpr = SCIPexprGetChildren(children[0])[0];
241  numexpr = SCIPexprGetChildren(children[1])[1];
242  }
243  else if( children[0] == SCIPexprGetChildren(children[1])[1] )
244  {
245  denomexpr = SCIPexprGetChildren(children[0])[0];
246  numexpr = SCIPexprGetChildren(children[1])[0];
247  }
248 
249  /* remember scalar and constant for numerator */
250  nomfac = sumcoefs[1] * prodcoef;
251  numconst = sumcoefs[0];
252  }
253  /* children[1] is 1/g(y) and children[0] is a product of f(x) and 1/g(y) */
254  else if( SCIPisExprPower(scip, children[1]) && SCIPgetExponentExprPow(children[1]) == -1.0
255  && SCIPisExprProduct(scip, children[0]) && SCIPexprGetNChildren(children[0]) == 2 ) /* lint !e777 */
256  {
257  SCIP_Real prodcoef = SCIPgetCoefExprProduct(children[0]);
258 
259  if( children[1] == SCIPexprGetChildren(children[0])[0] )
260  {
261  denomexpr = SCIPexprGetChildren(children[1])[0];
262  numexpr = SCIPexprGetChildren(children[0])[1];
263  }
264  else if( children[1] == SCIPexprGetChildren(children[0])[1] )
265  {
266  denomexpr = SCIPexprGetChildren(children[1])[0];
267  numexpr = SCIPexprGetChildren(children[0])[0];
268  }
269 
270  /* remember scalar and constant for numerator */
271  nomfac = sumcoefs[0] * prodcoef;
272  numconst = sumcoefs[1];
273  }
274 
275  /* remember the constant of the sum expression */
276  e = SCIPgetConstantExprSum(expr);
277  }
278 
279  if( denomexpr != NULL && numexpr != NULL )
280  {
281  /* transform numerator and denominator to detect structures like (a * f(x) + b) / (c * f(x) + d) */
282  transformExpr(scip, numexpr, &xexpr, &a, &b);
283  transformExpr(scip, denomexpr, &yexpr, &c, &d);
284 
285  SCIPdebugMsg(scip, "detected numerator (%g * %p + %g) and denominator (%g * %p + %g)\n", a, (void*)xexpr, b,
286  c, (void*)yexpr, d);
287 
288  /* detection is only be successful if the expression of the numerator an denominator are the same
289  * (so boundtightening can be stronger than default) or we are going to provide estimators (there will be an auxvar)
290  */
291  *success = (xexpr == yexpr) || (SCIPgetExprNAuxvarUsesNonlinear(expr) > 0);
292 
293 #ifdef SCIP_DEBUG
294  SCIPinfoMessage(scip, NULL, "Expression for numerator: ");
295  SCIP_CALL( SCIPprintExpr(scip, xexpr, NULL) );
296  SCIPinfoMessage(scip, NULL, "\nExpression for denominator: ");
297  SCIP_CALL( SCIPprintExpr(scip, yexpr, NULL) );
298  SCIPinfoMessage(scip, NULL, "\n");
299 #endif
300  }
301 
302  /* register usage of xexpr and yexpr
303  * create nonlinear handler expression data
304  */
305  if( *success )
306  {
307  assert(xexpr != NULL);
308  assert(xexpr != NULL);
309  assert(a != 0.0);
310  assert(c != 0.0);
311 
312  assert(SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 || xexpr == yexpr);
313 
314  /* request auxiliary variables for xexpr and yexpr if we will estimate
315  * mark that the bounds of the expression are important to construct the estimators
316  * (TODO check the curvature of the univariate quotient, as bounds may actually not be used)
317  * if univariate, then we also do inteval and reverseprop, so mark that the activities will be used for inteval
318  */
321  xexpr == yexpr,
323  SCIPgetExprNAuxvarUsesNonlinear(expr) > 0) );
324 
325  if( xexpr != yexpr && SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 )
326  {
328  }
329 
330  a = nomfac * a;
331  b = nomfac * b + numconst;
332 
333  SCIPdebug( SCIP_CALL( SCIPprintExpr(scip, expr, NULL) ); )
334  SCIPdebug( SCIPinfoMessage(scip, NULL, "\n") );
335  SCIPdebugMsg(scip, "detected quotient expression (%g * %p + %g) / (%g * %p + %g) + %g\n", a, (void*)xexpr,
336  b, c, (void*)yexpr, d, e);
337  SCIP_CALL( exprdataCreate(scip, nlhdlrexprdata, xexpr, a, b, yexpr, c, d, e) );
338  }
339 
340  return SCIP_OKAY;
341 }
342 
343 /** helper method to compute interval for (a x + b) / (c x + d) + e */
344 static
346  SCIP* scip, /**< SCIP data structure */
347  SCIP_INTERVAL bnds, /**< bounds on x */
348  SCIP_Real a, /**< coefficient in numerator */
349  SCIP_Real b, /**< constant in numerator */
350  SCIP_Real c, /**< coefficient in denominator */
351  SCIP_Real d, /**< constant in denominator */
352  SCIP_Real e /**< constant */
353  )
354 {
355  SCIP_INTERVAL result;
356  SCIP_INTERVAL denominterval;
357  SCIP_INTERVAL numinterval;
358  int i;
359 
360  assert(scip != NULL);
361 
362  /* return empty interval if the domain of x is empty */
364  {
365  SCIPintervalSetEmpty(&result);
366  return result;
367  }
368 
369  /* compute bounds for denominator */
370  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, bnds, c);
371  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d);
372 
373  /* there is no useful interval if 0 is in the interior of the interval of the denominator */
374  if( SCIPintervalGetInf(denominterval) < 0.0 && SCIPintervalGetSup(denominterval) > 0.0 )
375  {
377  return result;
378  }
379 
380  /* a d = b c implies that f(x) = b / d + e, i.e., f is constant */
381  if( a*d - b*c == 0.0 )
382  {
383  SCIPintervalSet(&result, b / d + e);
384  return result;
385  }
386 
387  /*
388  * evaluate for [x.inf,x.inf] and [x.sup,x.sup] independently
389  */
390  SCIPintervalSetEmpty(&result);
391 
392  for( i = 0; i < 2; ++i )
393  {
394  SCIP_INTERVAL quotinterval;
395  SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup;
396 
397  /* set the resulting interval to a / c if the bounds is infinite */
398  if( SCIPisInfinity(scip, REALABS(val)) )
399  {
400  SCIPintervalSet(&quotinterval, a);
401  SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, &quotinterval, quotinterval, c);
402  }
403  else
404  {
405  /* a x' + b */
406  SCIPintervalSet(&numinterval, val);
407  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, a);
408  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numinterval, numinterval, b);
409 
410  /* c x' + d */
411  SCIPintervalSet(&denominterval, val);
412  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, c);
413  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominterval, denominterval, d);
414 
415  /* (a x' + b) / (c x' + d) + e */
416  SCIPintervalDiv(SCIP_INTERVAL_INFINITY, &quotinterval, numinterval, denominterval);
417  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &quotinterval, quotinterval, e);
418  }
419 
420  /* unify with the resulting interval */
421  SCIPintervalUnify(&result, result, quotinterval);
422  }
423 
424  return result;
425 }
426 
427 /** helper method to compute reverse propagation for (a x + b) / (c x + d) + e */
428 static
430  SCIP_INTERVAL bnds, /**< bounds on (a x + b) / (c x + d) + e */
431  SCIP_Real a, /**< coefficient in numerator */
432  SCIP_Real b, /**< constant in numerator */
433  SCIP_Real c, /**< coefficient in denominator */
434  SCIP_Real d, /**< constant in denominator */
435  SCIP_Real e /**< constant */
436  )
437 {
438  SCIP_INTERVAL result;
439  int i;
440 
441  SCIPintervalSetEmpty(&result);
442 
443  /* return empty interval if the domain of the expression is empty */
445  return result;
446 
447  /* substract constant from bounds of the expression */
449 
450  /* if the expression is constant or the limit lies inside the domain, nothing can be propagated */
451  if( a*d - b*c == 0.0 || (bnds.inf < a / c && bnds.sup > a / c) )
452  {
454  return result;
455  }
456 
457  /* compute bounds for [x.inf,x.inf] and [x.sup,x.sup] independently */
458  for( i = 0; i < 2; ++i )
459  {
460  SCIP_INTERVAL denominator;
461  SCIP_INTERVAL numerator;
462  SCIP_INTERVAL quotient;
463  SCIP_Real val = (i == 0) ? bnds.inf : bnds.sup;
464 
465  /* (d * x' - b) */
466  SCIPintervalSet(&numerator, d);
467  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, val);
468  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &numerator, numerator, -b);
469 
470  /* (a - c * x') */
471  SCIPintervalSet(&denominator, -c);
472  SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, val);
473  SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &denominator, denominator, a);
474 
475  /* (d * x' - b) / (a - c * x') */
476  SCIPintervalDiv(SCIP_INTERVAL_INFINITY, &quotient, numerator, denominator);
477 
478  /* unify with the resulting interval */
479  SCIPintervalUnify(&result, result, quotient);
480  }
481 
482  return result;
483 }
484 
485 /** adds data to given rowprep; the generated estimator is always locally valid
486  *
487  * @note the constant is moved to the left- or right-hand side
488  * @note other than the name of this function may indicate, it does not create a rowprep
489  */
490 static
492  SCIP* scip, /**< SCIP data structure */
493  SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */
494  SCIP_VAR** vars, /**< variables */
495  SCIP_Real* coefs, /**< coefficients */
496  SCIP_Real constant, /**< constant */
497  int nlinvars /**< total number of variables */
498  )
499 {
500  assert(scip != NULL);
501  assert(rowprep != NULL);
502  assert(coefs != NULL);
503  assert(vars != NULL);
504 
505  /* create rowprep */
506  SCIProwprepAddSide(rowprep, -constant);
507  SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlinvars + 1) );
508 
509  /* add coefficients */
510  SCIP_CALL( SCIPaddRowprepTerms(scip, rowprep, nlinvars, vars, coefs) );
511 
512  return SCIP_OKAY;
513 }
514 
515 /** computes an estimator at a given point for the univariate case (ax + b) / (cx + d) + e
516  *
517  * Depending on the reference point, the estimator is a tangent or a secant on the graph.
518  * It depends on whether we are under- or overestimating, whether we are on the left or
519  * on the right side of the singularity at -d/c, and whether it is the monotone increasing
520  * (ad - bc > 0) or decreasing part (ad - bc < 0). Together, there are 8 cases:
521  *
522  * - mon. incr. + overestimate + left hand side --> secant
523  * - mon. incr. + overestimate + right hand side --> tangent
524  * - mon. incr. + understimate + left hand side --> tangent
525  * - mon. incr. + understimate + right hand side --> secant
526  * - mon. decr. + overestimate + left hand side --> tangent
527  * - mon. decr. + overestimate + right hand side --> secant
528  * - mon. decr. + understimate + left hand side --> secant
529  * - mon. decr. + understimate + right hand side --> tangent
530  */
531 static
533  SCIP* scip, /**< SCIP data structure */
534  SCIP_Real lbx, /**< local lower bound of x */
535  SCIP_Real ubx, /**< local upper bound of x */
536  SCIP_Real gllbx, /**< global lower bound of x */
537  SCIP_Real glubx, /**< global upper bound of x */
538  SCIP_Real solx, /**< solution value of x */
539  SCIP_Real a, /**< coefficient in numerator */
540  SCIP_Real b, /**< constant in numerator */
541  SCIP_Real c, /**< coefficient in denominator */
542  SCIP_Real d, /**< constant in denominator */
543  SCIP_Real e, /**< constant */
544  SCIP_Real* coef, /**< pointer to store the coefficient */
545  SCIP_Real* constant, /**< pointer to store the constant */
546  SCIP_Bool overestimate, /**< whether the expression should be overestimated */
547  SCIP_Bool* local, /**< pointer to store whether the estimate is locally valid */
548  SCIP_Bool* branchinguseful, /**< pointer to store whether branching on the expression would improve the estimator */
549  SCIP_Bool* success /**< buffer to store whether separation was successful */
550  )
551 {
552  SCIP_Real singularity;
553  SCIP_Bool isinleftpart;
554  SCIP_Bool monincreasing;
555 
556  assert(lbx <= solx && solx <= ubx);
557  assert(coef != NULL);
558  assert(constant != NULL);
559  assert(local != NULL);
560  assert(branchinguseful != NULL);
561  assert(success != NULL);
562 
563  *branchinguseful = TRUE;
564  *success = FALSE;
565  *coef = 0.0;
566  *constant = 0.0;
567  singularity = -d / c;
568 
569  /* estimate is globally valid if local and global bounds are equal */
570  *local = gllbx != lbx || glubx != ubx; /*lint !e777*/
571 
572  /* if 0 is in the denom interval, estimation is not possible */
573  if( SCIPisLE(scip, lbx, singularity) && SCIPisGE(scip, ubx, singularity) )
574  return SCIP_OKAY;
575 
576  isinleftpart = (ubx < singularity);
577  monincreasing = (a * d - b * c > 0.0);
578 
579  /* this encodes the 8 cases explained above */
580  if( monincreasing == (overestimate == isinleftpart) )
581  {
582  SCIP_Real lbeval;
583  SCIP_Real ubeval;
584 
585  /* if one of the bounds is infinite, secant cannot be computed */
586  if( SCIPisInfinity(scip, -lbx) || SCIPisInfinity(scip, ubx) )
587  return SCIP_OKAY;
588 
589  lbeval = (a * lbx + b) / (c * lbx + d) + e;
590  ubeval = (a * ubx + b) / (c * ubx + d) + e;
591 
592  /* compute coefficient and constant of linear estimator */
593  *coef = (ubeval - lbeval) / (ubx - lbx);
594  *constant = ubeval - (*coef) * ubx;
595  }
596  else
597  {
598  SCIP_Real soleval;
599 
600  soleval = (a * solx + b) / (c * solx + d) + e;
601 
602  /* compute coefficient and constant of linear estimator */
603  *coef = (a * d - b * c) / SQR(d + c * solx);
604  *constant = soleval - (*coef) * solx;
605 
606  /* gradient cuts are globally valid if the singularity is not in [gllbx,glubx] */
607  *local = SCIPisLE(scip, gllbx, singularity) && SCIPisGE(scip, glubx, singularity);
608 
609  /* branching will not improve the convexification via tangent cuts */
610  *branchinguseful = FALSE;
611  }
612 
613  /* avoid huge values in the cut */
614  if( SCIPisHugeValue(scip, REALABS(*coef)) || SCIPisHugeValue(scip, REALABS(*constant)) )
615  return SCIP_OKAY;
616 
617  *success = TRUE;
618 
619  return SCIP_OKAY;
620 }
621 
622 /** helper method to compute estimator for the univariate case; the estimator is stored in a given rowprep */
623 static
625  SCIP* scip, /**< SCIP data structure */
626  SCIP_SOL* sol, /**< solution point (or NULL for the LP solution) */
627  SCIP_EXPR* xexpr, /**< argument expression */
628  SCIP_Real a, /**< coefficient in numerator */
629  SCIP_Real b, /**< constant in numerator */
630  SCIP_Real c, /**< coefficient in denominator */
631  SCIP_Real d, /**< constant in denominator */
632  SCIP_Real e, /**< constant */
633  SCIP_Bool overestimate, /**< whether the expression should be overestimated */
634  SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */
635  SCIP_Bool* branchinguseful, /**< pointer to store whether branching on the expression would improve the estimator */
636  SCIP_Bool* success /**< buffer to store whether separation was successful */
637  )
638 {
639  SCIP_VAR* x;
640  SCIP_Real constant;
641  SCIP_Real coef;
642  SCIP_Real gllbx;
643  SCIP_Real glubx;
644  SCIP_Real lbx;
645  SCIP_Real ubx;
646  SCIP_Real solx;
647  SCIP_Bool local;
648  SCIP_INTERVAL bnd;
649 
650  assert(rowprep != NULL);
651  assert(branchinguseful != NULL);
652  assert(success != NULL);
653 
654  x = SCIPgetExprAuxVarNonlinear(xexpr);
655 
656  /* get local bounds on xexpr */
660  SCIP_CALL( SCIPevalExprActivity(scip, xexpr) );
661  SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr));
662  lbx = bnd.inf;
663  ubx = bnd.sup;
664 
665  /* check whether variable has been fixed or has empty interval */
666  if( SCIPisEQ(scip, lbx, ubx) || ubx < lbx )
667  {
668  *success = FALSE;
669  return SCIP_OKAY;
670  }
671 
672  /* get global variable bounds */
673  gllbx = SCIPvarGetLbGlobal(x);
674  glubx = SCIPvarGetUbGlobal(x);
675 
676  /* get and adjust solution value */
677  solx = SCIPgetSolVal(scip, sol, x);
678  solx = MIN(MAX(solx, lbx), ubx);
679 
680  /* compute an estimator */
681  SCIP_CALL( estimateUnivariate(scip, lbx, ubx, gllbx, glubx, solx, a, b, c, d, e, &coef, &constant, overestimate, &local, branchinguseful, success) );
682 
683  /* add estimator to rowprep, if successful */
684  if( *success )
685  {
686  (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%lld", SCIPvarGetName(x), SCIPgetNLPs(scip));
687  SCIP_CALL( createRowprep(scip, rowprep, &x, &coef, constant, 1) );
688  SCIProwprepSetLocal(rowprep, local);
689  }
690 
691  return SCIP_OKAY;
692 }
693 
694 /** helper method to compute a gradient cut for
695  * \f[
696  * h^c(x,y) := \frac{1}{y} \left(\frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx}} + \sqrt{\text{ubx}}}\right)^2
697  * \f]
698  * at a given reference point
699  *
700  * See Zamora and Grossmann (1988) for more details.
701  */
702 static
704  SCIP_Real lbx, /**< lower bound of x */
705  SCIP_Real ubx, /**< upper bound of x */
706  SCIP_Real solx, /**< solution value of x */
707  SCIP_Real soly, /**< solution value of y */
708  SCIP_Real* coefx, /**< pointer to store the coefficient of x */
709  SCIP_Real* coefy, /**< pointer to store the coefficient of y */
710  SCIP_Real* constant /**< pointer to store the constant */
711  )
712 {
713  SCIP_Real tmp1;
714  SCIP_Real tmp2;
715 
716  assert(lbx >= 0.0);
717  assert(lbx <= ubx);
718  assert(soly > 0.0);
719  assert(coefx != NULL);
720  assert(coefy != NULL);
721  assert(constant != NULL);
722 
723  tmp1 = SQRT(lbx * ubx) + solx;
724  tmp2 = SQR(SQRT(lbx) + SQRT(ubx)) * soly; /*lint !e666*/
725  assert(tmp2 > 0.0);
726 
727  *coefx = 2.0 * tmp1 / tmp2;
728  *coefy = -SQR(tmp1) / (tmp2 * soly);
729  *constant = 2.0 * SQRT(lbx * ubx) * tmp1 / tmp2;
730 }
731 
732 /** computes an over- or underestimator at a given point for the bivariate case x/y &le;/&ge; z
733  *
734  * There are the following cases for y > 0:
735  *
736  * 1. lbx < 0 < ubx:
737  * Rewrite x / y = z as x = y * z and use McCormick to compute a valid inequality of the form
738  * x = y * z &le; a * y + b * z + c. Note that b > 0 because of y > 0. The inequality is then transformed
739  * to x / b - a/b * y - c/b &le; z, which results in a valid underestimator for x / y over the set
740  * {(x,y) | lbz &le; x / y &le; ubz}. Note that overestimating/underestimating the bilinear term with McCormick
741  * results in an underestimator/overestimator for x / y.
742  *
743  * 2. lbx &ge; 0 or ubx &le; 0:
744  * - overestimation: use \f$z \leq \frac{1}{\text{lby}\cdot\text{uby}} \min(\text{uby}\cdot x - \text{lbx}\cdot y + \text{lbx}\cdot\text{lby}, \text{lby}\cdot x - \text{ubx}\cdot y + \text{ubx}\cdot\text{uby})\f$
745  * - underestimation: use \f$z \geq x/y \geq \frac{1}{y} \frac{x + \sqrt{\text{lbx}\cdot\text{ubx}}}{\sqrt{\text{lbx} + \sqrt{\text{ubx}}}}\f$ and build gradient cut
746  *
747  * If y < 0, swap and negate its bounds and compute the respective opposite estimator (and negate it).
748  *
749  * If 0 is in the interval of y, nothing is possible.
750  */
751 static
753  SCIP* scip, /**< SCIP data structure */
754  SCIP_Real lbx, /**< lower bound of x */
755  SCIP_Real ubx, /**< upper bound of x */
756  SCIP_Real lby, /**< lower bound of y */
757  SCIP_Real uby, /**< upper bound of y */
758  SCIP_Real lbz, /**< lower bound of z */
759  SCIP_Real ubz, /**< lower bound of z */
760  SCIP_Real solx, /**< reference point for x */
761  SCIP_Real soly, /**< reference point for y */
762  SCIP_Real solz, /**< reference point for z */
763  SCIP_Bool overestimate, /**< whether the expression should be overestimated */
764  SCIP_Real* coefx, /**< pointer to store the x coefficient */
765  SCIP_Real* coefy, /**< pointer to store the y coefficient */
766  SCIP_Real* constant, /**< pointer to store the constant */
767  SCIP_Bool* branchingusefulx, /**< pointer to store whether branching on x would improve the estimator */
768  SCIP_Bool* branchingusefuly, /**< pointer to store whether branching on y would improve the estimator */
769  SCIP_Bool* success /**< buffer to store whether computing the estimator was successful */
770  )
771 {
772  SCIP_Bool negatedx = FALSE;
773  SCIP_Bool negatedy = FALSE;
774 
775  assert(lbx <= solx && solx <= ubx);
776  assert(lby <= soly && soly <= uby);
777  assert(lbz <= solz && solz <= ubz);
778  assert(coefx != NULL);
779  assert(coefy != NULL);
780  assert(constant != NULL);
781  assert(branchingusefulx != NULL);
782  assert(branchingusefuly != NULL);
783  assert(success != NULL);
784 
785  *branchingusefulx = TRUE;
786  *branchingusefuly = TRUE;
787  *success = TRUE;
788  *coefx = 0.0;
789  *coefy = 0.0;
790  *constant = 0.0;
791 
792  /* if 0 is in [lby,uby], then it is not possible to compute an estimator */
793  if( SCIPisLE(scip, lby, 0.0) && SCIPisGE(scip, uby, 0.0) )
794  {
795  *success = FALSE;
796  return SCIP_OKAY;
797  }
798 
799  /* negate bounds of y if it is not positive */
800  if( uby < 0.0 )
801  {
802  SCIP_Real tmp = uby;
803 
804  uby = -lby;
805  lby = -tmp;
806  soly = -soly;
807  negatedy = TRUE;
808  overestimate = !overestimate;
809  }
810 
811  /* case 1: 0 is in the interior of [lbx,ubx] */
812  if( lbx < 0.0 && 0.0 < ubx )
813  {
814  SCIP_Real mccoefy = 0.0;
815  SCIP_Real mccoefaux = 0.0;
816  SCIP_Real mcconst = 0.0;
817 
818  /* as explained in the description of this method, overestimating/underestimating the bilinear term results in an
819  * underestimator/overestimator for x / y
820  */
821  SCIPaddBilinMcCormick(scip, 1.0, lbz, ubz, solz, lby, uby, soly, !overestimate, &mccoefaux, &mccoefy, &mcconst,
822  success);
823  assert(mccoefaux >= 0.0);
824 
825  if( !(*success) )
826  return SCIP_OKAY;
827 
828  /* resulting estimator is x/b - a/b * y - c/b, where a*y + b*z + c is the estimator for y*z */
829  *coefx = 1.0 / mccoefaux;
830  *coefy = -mccoefy / mccoefaux;
831  *constant = -mcconst / mccoefaux;
832  }
833  /* case 2: 0 is not in the interior of [lbx,ubx] */
834  else
835  {
836  /* negate bounds of x if it is negative */
837  if( ubx <= 0.0 )
838  {
839  SCIP_Real tmp = ubx;
840 
841  ubx = -lbx;
842  lbx = -tmp;
843  solx = -solx;
844  negatedx = TRUE;
845  overestimate = !overestimate;
846  }
847 
848  /* case 2a */
849  if( overestimate )
850  {
851  /* check where the minimum is attained */
852  if( uby * solx - lbx * soly + lbx * lby <= lby * solx - ubx * soly + ubx * uby )
853  {
854  *coefx = 1.0 / lby;
855  *coefy = -lbx / (lby * uby);
856  *constant = lbx / uby;
857  }
858  else
859  {
860  *coefx = 1.0 / uby;
861  *coefy = -ubx / (lby * uby);
862  *constant = ubx / lby;
863  }
864  }
865  /* case 2b */
866  else
867  {
868  /* compute gradient cut for h^c(x,y) at (solx,soly) */
869  hcGradCut(lbx, ubx, solx, soly, coefx, coefy, constant);
870 
871  /* estimator is independent of the bounds of y */
872  *branchingusefuly = FALSE;
873  }
874  }
875 
876  /* reverse negations of x and y in the resulting estimator */
877  if( negatedx )
878  *coefx = -(*coefx);
879  if( negatedy )
880  *coefy = -(*coefy);
881 
882  /* if exactly one variable has been negated, then we have computed an underestimate/overestimate for the negated
883  * expression, which results in an overestimate/underestimate for the original expression
884  */
885  if( negatedx != negatedy )
886  {
887  *coefx = -(*coefx);
888  *coefy = -(*coefy);
889  *constant = -(*constant);
890  }
891 
892  /* avoid huge values in the estimator */
893  if( SCIPisHugeValue(scip, REALABS(*coefx)) || SCIPisHugeValue(scip, REALABS(*coefy))
894  || SCIPisHugeValue(scip, REALABS(*constant)) )
895  {
896  *success = FALSE;
897  return SCIP_OKAY;
898  }
899 
900  return SCIP_OKAY;
901 }
902 
903 /** construct an estimator for a quotient expression of the form (ax + b) / (cy + d) + e
904  *
905  * The resulting estimator is stored in a rowprep.
906  *
907  * The method first computes an estimator for x' / y' with x := ax + b and y := cy + d
908  * and then transforms this estimator to one for the quotient (ax + b) / (cy + d) + e.
909  */
910 static
912  SCIP* scip, /**< SCIP data structure */
913  SCIP_EXPR* xexpr, /**< numerator expression */
914  SCIP_EXPR* yexpr, /**< denominator expression */
915  SCIP_VAR* auxvar, /**< auxiliary variable */
916  SCIP_SOL* sol, /**< solution point (or NULL for the LP solution) */
917  SCIP_Real a, /**< coefficient of numerator */
918  SCIP_Real b, /**< constant of numerator */
919  SCIP_Real c, /**< coefficient of denominator */
920  SCIP_Real d, /**< constant of denominator */
921  SCIP_Real e, /**< constant term */
922  SCIP_Bool overestimate, /**< whether the expression should be overestimated */
923  SCIP_ROWPREP* rowprep, /**< a rowprep where to store the estimator */
924  SCIP_Bool* branchingusefulx, /**< pointer to store whether branching on x would improve the estimator */
925  SCIP_Bool* branchingusefuly, /**< pointer to store whether branching on y would improve the estimator */
926  SCIP_Bool* success /**< buffer to store whether separation was successful */
927  )
928 {
929  SCIP_VAR* vars[2];
930  SCIP_Real coefs[2] = {0.0, 0.0};
931  SCIP_Real constant = 0.0;
932  SCIP_Real solx;
933  SCIP_Real soly;
934  SCIP_Real solz;
935  SCIP_Real lbx;
936  SCIP_Real ubx;
937  SCIP_Real lby;
938  SCIP_Real uby;
939  SCIP_Real lbz;
940  SCIP_Real ubz;
941  SCIP_INTERVAL bnd;
942 
943  assert(xexpr != NULL);
944  assert(yexpr != NULL);
945  assert(xexpr != yexpr);
946  assert(auxvar != NULL);
947  assert(rowprep != NULL);
948  assert(branchingusefulx != NULL);
949  assert(branchingusefuly != NULL);
950  assert(success != NULL);
951 
952  vars[0] = SCIPgetExprAuxVarNonlinear(xexpr);
953  vars[1] = SCIPgetExprAuxVarNonlinear(yexpr);
954 
955  /* get bounds for x, y, and z */
959  SCIP_CALL( SCIPevalExprActivity(scip, xexpr) );
960  SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr));
961  lbx = bnd.inf;
962  ubx = bnd.sup;
963 
967  SCIP_CALL( SCIPevalExprActivity(scip, yexpr) );
968  SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(yexpr));
969  lby = bnd.inf;
970  uby = bnd.sup;
971 
972  lbz = SCIPvarGetLbLocal(auxvar);
973  ubz = SCIPvarGetUbLocal(auxvar);
974 
975  /* check whether one of the variables has been fixed or has empty domain */
976  if( SCIPisEQ(scip, lbx, ubx) || SCIPisEQ(scip, lby, uby) || ubx < lbx || uby < lby )
977  {
978  *success = FALSE;
979  return SCIP_OKAY;
980  }
981 
982  /* get and adjust solution values */
983  solx = SCIPgetSolVal(scip, sol, vars[0]);
984  soly = SCIPgetSolVal(scip, sol, vars[1]);
985  solz = SCIPgetSolVal(scip, sol, auxvar);
986  solx = MIN(MAX(solx, lbx), ubx);
987  soly = MIN(MAX(soly, lby), uby);
988  solz = MIN(MAX(solz, lbz), ubz);
989 
990  /* compute an estimator */
992  MIN(a * lbx, a * ubx) + b, MAX(a * lbx, a * ubx) + b, /* bounds of x' */
993  MIN(c * lby, c * uby) + d, MAX(c * lby, c * uby) + d, /* bounds of y' */
994  lbz, ubz, a * solx + b, c * soly + d, solz, overestimate, &coefs[0], &coefs[1], &constant,
995  branchingusefulx, branchingusefuly, success) );
996 
997  /* add estimator to rowprep, if successful */
998  if( *success )
999  {
1000  /* transform estimator Ax' + By'+ C = A(ax + b) + B (cy + d) + C = (Aa) x + (Bc) y + (C + Ab + Bd);
1001  * add the constant e separately
1002  */
1003  constant += coefs[0] * b + coefs[1] * d + e;
1004  coefs[0] *= a;
1005  coefs[1] *= c;
1006 
1007  /* prepare rowprep */
1008  (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "quot_%s_%s_%lld", SCIPvarGetName(vars[0]), SCIPvarGetName(vars[1]),
1009  SCIPgetNLPs(scip));
1010  SCIP_CALL( createRowprep(scip, rowprep, vars, coefs, constant, 2) );
1011  }
1012 
1013  return SCIP_OKAY;
1014 }
1015 
1016 /*
1017  * Callback methods of nonlinear handler
1018  */
1019 
1020 /** nonlinear handler copy callback */
1021 static
1022 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuotient)
1023 { /*lint --e{715}*/
1024  assert(targetscip != NULL);
1025  assert(sourcenlhdlr != NULL);
1026  assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), NLHDLR_NAME) == 0);
1027 
1028  SCIP_CALL( SCIPincludeNlhdlrQuotient(targetscip) );
1029 
1030  return SCIP_OKAY;
1031 }
1032 
1033 
1034 /** callback to free expression specific data */
1035 static
1036 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataQuotient)
1037 { /*lint --e{715}*/
1038  assert(nlhdlrexprdata != NULL);
1039  assert(*nlhdlrexprdata != NULL);
1040 
1041  /* free expression data of nonlinear handler */
1042  SCIP_CALL( exprdataFree(scip, nlhdlrexprdata) );
1043 
1044  return SCIP_OKAY;
1045 }
1046 
1047 
1048 /** callback to detect structure in expression tree */
1049 static
1050 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuotient)
1051 { /*lint --e{715}*/
1052  SCIP_Bool success;
1053 
1054  assert(nlhdlrexprdata != NULL);
1055 
1056  /* call detection routine */
1057  SCIP_CALL( detectExpr(scip, expr, nlhdlrexprdata, &success) );
1058 
1059  if( success )
1060  {
1061  if( SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 )
1062  *participating = SCIP_NLHDLR_METHOD_SEPABOTH;
1063 
1064  if( (*nlhdlrexprdata)->numexpr == (*nlhdlrexprdata)->denomexpr )
1065  {
1066  /* if univariate, then we also do inteval and reverseprop */
1067  *participating |= SCIP_NLHDLR_METHOD_ACTIVITY;
1068 
1069  /* if univariate, then all our methods are enforcing */
1070  *enforcing |= *participating;
1071  }
1072  }
1073 
1074  return SCIP_OKAY;
1075 }
1076 
1077 
1078 /** auxiliary evaluation callback of nonlinear handler */
1079 static
1080 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuotient)
1081 { /*lint --e{715}*/
1082  SCIP_VAR* auxvarx;
1083  SCIP_VAR* auxvary;
1084  SCIP_Real solvalx;
1085  SCIP_Real solvaly;
1086  SCIP_Real nomval;
1087  SCIP_Real denomval;
1088 
1089  assert(expr != NULL);
1090  assert(auxvalue != NULL);
1091 
1092  /**! [SnippetNlhdlrEvalauxQuotient] */
1093  /* get auxiliary variables */
1094  auxvarx = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->numexpr);
1095  auxvary = SCIPgetExprAuxVarNonlinear(nlhdlrexprdata->denomexpr);
1096  assert(auxvarx != NULL);
1097  assert(auxvary != NULL);
1098 
1099  /* get solution values of the auxiliary variables */
1100  solvalx = SCIPgetSolVal(scip, sol, auxvarx);
1101  solvaly = SCIPgetSolVal(scip, sol, auxvary);
1102 
1103  /* evaluate expression w.r.t. the values of the auxiliary variables */
1104  nomval = nlhdlrexprdata->numcoef * solvalx + nlhdlrexprdata->numconst;
1105  denomval = nlhdlrexprdata->denomcoef * solvaly + nlhdlrexprdata->denomconst;
1106 
1107  /* return SCIP_INVALID if the denominator evaluates to zero */
1108  *auxvalue = (denomval != 0.0) ? nlhdlrexprdata->constant + nomval / denomval : SCIP_INVALID;
1109  /**! [SnippetNlhdlrEvalauxQuotient] */
1110 
1111  return SCIP_OKAY;
1112 }
1113 
1114 
1115 /** nonlinear handler under/overestimation callback
1116  *
1117  * @todo which of the paramters did I not use, but have to be taken into consideration?
1118 */
1119 static
1120 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateQuotient)
1121 { /*lint --e{715}*/
1122  SCIP_Bool branchingusefulx = FALSE;
1123  SCIP_Bool branchingusefuly = FALSE;
1124  SCIP_ROWPREP* rowprep;
1125 
1126  assert(nlhdlr != NULL);
1127  assert(expr != NULL);
1128  assert(nlhdlrexprdata != NULL);
1129  assert(rowpreps != NULL);
1130 
1131  /** ![SnippetNlhdlrEstimateQuotient] */
1132  *addedbranchscores = FALSE;
1133  *success = FALSE;
1134 
1136 
1137  if( nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr )
1138  {
1139  /* univariate case */
1140  SCIP_CALL( estimateUnivariateQuotient(scip, sol, nlhdlrexprdata->numexpr, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1141  nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant, overestimate, rowprep,
1142  &branchingusefulx, success) );
1143  }
1144  else
1145  {
1146  /* bivariate case */
1147  SCIP_CALL( estimateBivariateQuotient(scip, nlhdlrexprdata->numexpr, nlhdlrexprdata->denomexpr, SCIPgetExprAuxVarNonlinear(expr), sol,
1148  nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst, nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst,
1149  nlhdlrexprdata->constant, overestimate, rowprep,
1150  &branchingusefulx, &branchingusefuly, success) );
1151  }
1152 
1153  if( *success )
1154  {
1155  SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) );
1156  }
1157  else
1158  {
1159  SCIPfreeRowprep(scip, &rowprep);
1160  }
1161 
1162  /* add branching scores if requested */
1163  if( addbranchscores )
1164  {
1165  SCIP_EXPR* exprs[2];
1166  SCIP_Real violation;
1167  int nexprs = 0;
1168 
1169  if( branchingusefulx )
1170  exprs[nexprs++] = nlhdlrexprdata->numexpr;
1171  if( branchingusefuly )
1172  exprs[nexprs++] = nlhdlrexprdata->denomexpr;
1173 
1174  /* compute violation w.r.t. the auxiliary variable(s) */
1175 #ifndef BRSCORE_ABSVIOL
1176  SCIP_CALL( SCIPgetExprRelAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) );
1177 #else
1178  SCIP_CALL( SCIPgetExprAbsAuxViolationNonlinear(scip, expr, auxvalue, sol, &violation, NULL, NULL) );
1179 #endif
1180  assert(violation > 0.0); /* there should be a violation if we were called to enforce */
1181 
1182  SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, exprs, nexprs, violation, sol, addedbranchscores) );
1183  }
1184  /** ![SnippetNlhdlrEstimateQuotient] */
1185 
1186  return SCIP_OKAY;
1187 }
1188 
1189 
1190 /** nonlinear handler interval evaluation callback */
1191 static
1192 SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuotient)
1193 { /*lint --e{715}*/
1194  SCIP_INTERVAL bnds;
1195 
1196  assert(nlhdlrexprdata != NULL);
1197  assert(nlhdlrexprdata->numexpr != NULL);
1198  assert(nlhdlrexprdata->denomexpr != NULL);
1199 
1200  /* it is not possible to compute tighter intervals if both expressions are different
1201  * we should not be called in this case, as we haven't said we would participate in this activity in detect
1202  */
1203  assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr);
1204 
1205  /**! [SnippetNlhdlrIntevalQuotient] */
1206  /* get activity of the numerator (= denominator) expression */
1207  bnds = SCIPexprGetActivity(nlhdlrexprdata->numexpr);
1208 
1209  /* call interval evaluation for the univariate quotient expression */
1210  *interval = intEvalQuotient(scip, bnds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1211  nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant);
1212  /**! [SnippetNlhdlrIntevalQuotient] */
1213 
1214  return SCIP_OKAY;
1215 }
1216 
1217 
1218 /** nonlinear handler callback for reverse propagation */
1219 static
1220 SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuotient)
1221 { /*lint --e{715}*/
1222  SCIP_INTERVAL result;
1223 
1224  assert(nlhdlrexprdata != NULL);
1225  assert(nlhdlrexprdata->numexpr != NULL);
1226  assert(nlhdlrexprdata->denomexpr != NULL);
1227 
1228  /* it is not possible to compute tighter intervals if both expressions are different
1229  * we should not be called in this case, as we haven't said we would participate in this activity in detect
1230  */
1231  assert(nlhdlrexprdata->numexpr == nlhdlrexprdata->denomexpr);
1232 
1233  SCIPdebugMsg(scip, "call reverse propagation for expression (%g %p + %g) / (%g %p + %g) + %g bounds [%g,%g]\n",
1234  nlhdlrexprdata->numcoef, (void*)nlhdlrexprdata->numexpr, nlhdlrexprdata->numconst,
1235  nlhdlrexprdata->denomcoef, (void*)nlhdlrexprdata->denomexpr, nlhdlrexprdata->denomconst,
1236  nlhdlrexprdata->constant, bounds.inf, bounds.sup);
1237 
1238  /* call reverse propagation */
1239  /**! [SnippetNlhdlrReversepropQuotient] */
1240  result = reversepropQuotient(bounds, nlhdlrexprdata->numcoef, nlhdlrexprdata->numconst,
1241  nlhdlrexprdata->denomcoef, nlhdlrexprdata->denomconst, nlhdlrexprdata->constant);
1242 
1243  SCIPdebugMsg(scip, "try to tighten bounds of %p: [%g,%g] -> [%g,%g]\n",
1244  (void*)nlhdlrexprdata->numexpr, SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).inf,
1245  SCIPgetExprBoundsNonlinear(scip, nlhdlrexprdata->numexpr).sup, result.inf, result.sup);
1246 
1247  /* tighten bounds of the expression */
1248  SCIP_CALL( SCIPtightenExprIntervalNonlinear(scip, nlhdlrexprdata->numexpr, result, infeasible, nreductions) );
1249  /**! [SnippetNlhdlrReversepropQuotient] */
1250 
1251  return SCIP_OKAY;
1252 }
1253 
1254 
1255 /*
1256  * nonlinear handler specific interface methods
1257  */
1258 
1259 /** includes quotient nonlinear handler in nonlinear constraint handler */
1261  SCIP* scip /**< SCIP data structure */
1262  )
1263 {
1264  SCIP_NLHDLRDATA* nlhdlrdata;
1265  SCIP_NLHDLR* nlhdlr;
1266 
1267  assert(scip != NULL);
1268 
1269  /* create nonlinear handler data */
1270  nlhdlrdata = NULL;
1271 
1273  NLHDLR_ENFOPRIORITY, nlhdlrDetectQuotient, nlhdlrEvalauxQuotient, nlhdlrdata) );
1274  assert(nlhdlr != NULL);
1275 
1276  SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrQuotient);
1277  SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrFreeExprDataQuotient);
1278  SCIPnlhdlrSetSepa(nlhdlr, NULL, NULL, nlhdlrEstimateQuotient, NULL);
1279  SCIPnlhdlrSetProp(nlhdlr, nlhdlrIntevalQuotient, nlhdlrReversepropQuotient);
1280 
1281  return SCIP_OKAY;
1282 }
void SCIPintervalSetEntire(SCIP_Real infinity, SCIP_INTERVAL *resultant)
void SCIPintervalUnify(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition: scip_expr.c:1476
unsigned int SCIPgetExprNAuxvarUsesNonlinear(SCIP_EXPR *expr)
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1706
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3798
static SCIP_DECL_NLHDLRDETECT(nlhdlrDetectQuotient)
SCIP_Real SCIPvarGetLbGlobal(SCIP_VAR *var)
Definition: var.c:17910
#define SCIP_MAXSTRLEN
Definition: def.h:293
void SCIPnlhdlrSetFreeExprData(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRFREEEXPRDATA((*freeexprdata)))
Definition: nlhdlr.c:85
SCIP_Real SCIPvarGetLbLocal(SCIP_VAR *var)
Definition: var.c:17966
SCIP_Real SCIPgetConstantExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1172
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPregisterExprUsageNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Bool useauxvar, SCIP_Bool useactivityforprop, SCIP_Bool useactivityforsepabelow, SCIP_Bool useactivityforsepaabove)
static SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrQuotient)
#define FALSE
Definition: def.h:87
char * SCIProwprepGetName(SCIP_ROWPREP *rowprep)
Definition: misc_rowprep.c:664
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10755
void SCIPintervalIntersectEps(SCIP_INTERVAL *resultant, SCIP_Real eps, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3343
#define TRUE
Definition: def.h:86
#define SCIPdebug(x)
Definition: pub_message.h:84
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
static void transformExpr(SCIP *scip, SCIP_EXPR *expr, SCIP_EXPR **target, SCIP_Real *coef, SCIP_Real *constant)
static SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateQuotient)
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition: expr.c:3954
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:99
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void SCIPcaptureExpr(SCIP_EXPR *expr)
Definition: scip_expr.c:1399
SCIP_RETCODE SCIPensureRowprepSize(SCIP *scip, SCIP_ROWPREP *rowprep, int size)
Definition: misc_rowprep.c:855
#define SCIPdebugMsg
Definition: scip_message.h:69
SCIP_VAR ** x
Definition: circlepacking.c:54
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
static SCIP_INTERVAL intEvalQuotient(SCIP *scip, SCIP_INTERVAL bnds, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e)
static SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrFreeExprDataQuotient)
SCIP_Real SCIPepsilon(SCIP *scip)
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1454
bilinear nonlinear handler
#define NLHDLR_DETECTPRIORITY
static void hcGradCut(SCIP_Real lbx, SCIP_Real ubx, SCIP_Real solx, SCIP_Real soly, SCIP_Real *coefx, SCIP_Real *coefy, SCIP_Real *constant)
void SCIPintervalDiv(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
static SCIP_RETCODE estimateUnivariate(SCIP *scip, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real gllbx, SCIP_Real glubx, SCIP_Real solx, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e, SCIP_Real *coef, SCIP_Real *constant, SCIP_Bool overestimate, SCIP_Bool *local, SCIP_Bool *branchinguseful, SCIP_Bool *success)
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3808
#define NLHDLR_ENFOPRIORITY
SCIP_Real SCIPvarGetUbGlobal(SCIP_VAR *var)
Definition: var.c:17920
SCIP_Real inf
Definition: intervalarith.h:46
SCIP_Real * SCIPgetCoefsExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1157
static SCIP_RETCODE exprdataFree(SCIP *scip, SCIP_NLHDLREXPRDATA **nlhdlrexprdata)
static SCIP_RETCODE detectExpr(SCIP *scip, SCIP_EXPR *expr, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_Bool *success)
SCIP_RETCODE SCIPincludeNlhdlrQuotient(SCIP *scip)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
const char * SCIPnlhdlrGetName(SCIP_NLHDLR *nlhdlr)
Definition: nlhdlr.c:141
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)
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
const char * SCIPvarGetName(SCIP_VAR *var)
Definition: var.c:17251
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalMulScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define REALABS(x)
Definition: def.h:201
SCIP_Bool SCIPisExprSum(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1443
void SCIPfreeRowprep(SCIP *scip, SCIP_ROWPREP **rowprep)
Definition: misc_rowprep.c:558
#define SCIP_CALL(x)
Definition: def.h:384
void SCIProwprepAddSide(SCIP_ROWPREP *rowprep, SCIP_Real side)
Definition: misc_rowprep.c:714
static SCIP_RETCODE estimateBivariate(SCIP *scip, SCIP_Real lbx, SCIP_Real ubx, SCIP_Real lby, SCIP_Real uby, SCIP_Real lbz, SCIP_Real ubz, SCIP_Real solx, SCIP_Real soly, SCIP_Real solz, SCIP_Bool overestimate, SCIP_Real *coefx, SCIP_Real *coefy, SCIP_Real *constant, SCIP_Bool *branchingusefulx, SCIP_Bool *branchingusefuly, SCIP_Bool *success)
SCIP_Real sup
Definition: intervalarith.h:47
#define infty2infty(infty1, infty2, val)
quotient nonlinear handler
#define SCIP_INTERVAL_INFINITY
Definition: def.h:199
SCIP_Real SCIPgetCoefExprProduct(SCIP_EXPR *expr)
SCIP_Bool SCIPisHugeValue(SCIP *scip, SCIP_Real val)
void SCIPnlhdlrSetSepa(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINITSEPA((*initsepa)), SCIP_DECL_NLHDLRENFO((*enfo)), SCIP_DECL_NLHDLRESTIMATE((*estimate)), SCIP_DECL_NLHDLREXITSEPA((*exitsepa)))
Definition: nlhdlr.c:123
#define SCIP_Bool
Definition: def.h:84
static SCIP_RETCODE estimateUnivariateQuotient(SCIP *scip, SCIP_SOL *sol, SCIP_EXPR *xexpr, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e, SCIP_Bool overestimate, SCIP_ROWPREP *rowprep, SCIP_Bool *branchinguseful, SCIP_Bool *success)
static SCIP_DECL_NLHDLRINTEVAL(nlhdlrIntevalQuotient)
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1407
constraint handler for nonlinear constraints specified by algebraic expressions
#define MAX(x, y)
Definition: tclique_def.h:83
void SCIPnlhdlrSetProp(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRINTEVAL((*inteval)), SCIP_DECL_NLHDLRREVERSEPROP((*reverseprop)))
Definition: nlhdlr.c:110
static SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalauxQuotient)
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)
void SCIPintervalSubScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1465
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
static SCIP_RETCODE createRowprep(SCIP *scip, SCIP_ROWPREP *rowprep, SCIP_VAR **vars, SCIP_Real *coefs, SCIP_Real constant, int nlinvars)
static SCIP_RETCODE exprdataCreate(SCIP *scip, SCIP_NLHDLREXPRDATA **nlhdlrexprdata, SCIP_EXPR *numexpr, SCIP_Real numcoef, SCIP_Real numconst, SCIP_EXPR *denomexpr, SCIP_Real denomcoef, SCIP_Real denomconst, SCIP_Real constant)
SCIP_RETCODE SCIPtightenExprIntervalNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_INTERVAL newbounds, SCIP_Bool *cutoff, int *ntightenings)
static SCIP_INTERVAL reversepropQuotient(SCIP_INTERVAL bnds, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e)
void SCIPintervalDivScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_VAR ** b
Definition: circlepacking.c:56
#define NLHDLR_DESC
#define SCIP_NLHDLR_METHOD_ACTIVITY
Definition: type_nlhdlr.h:45
SCIP_VAR * a
Definition: circlepacking.c:57
SCIP_RETCODE SCIPgetExprAbsAuxViolationNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Real auxvalue, SCIP_SOL *sol, SCIP_Real *viol, SCIP_Bool *violunder, SCIP_Bool *violover)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
SCIP_RETCODE SCIPaddExprsViolScoreNonlinear(SCIP *scip, SCIP_EXPR **exprs, int nexprs, SCIP_Real violscore, SCIP_SOL *sol, SCIP_Bool *success)
#define SCIP_Real
Definition: def.h:177
SCIP_RETCODE SCIPgetExprRelAuxViolationNonlinear(SCIP *scip, SCIP_EXPR *expr, SCIP_Real auxvalue, SCIP_SOL *sol, SCIP_Real *viol, SCIP_Bool *violunder, SCIP_Bool *violover)
#define SCIP_INVALID
Definition: def.h:197
static SCIP_RETCODE estimateBivariateQuotient(SCIP *scip, SCIP_EXPR *xexpr, SCIP_EXPR *yexpr, SCIP_VAR *auxvar, SCIP_SOL *sol, SCIP_Real a, SCIP_Real b, SCIP_Real c, SCIP_Real d, SCIP_Real e, SCIP_Bool overestimate, SCIP_ROWPREP *rowprep, SCIP_Bool *branchingusefulx, SCIP_Bool *branchingusefuly, SCIP_Bool *success)
struct SCIP_NlhdlrExprData SCIP_NLHDLREXPRDATA
Definition: type_nlhdlr.h:404
SCIP_RETCODE SCIPcreateRowprep(SCIP *scip, SCIP_ROWPREP **rowprep, SCIP_SIDETYPE sidetype, SCIP_Bool local)
Definition: misc_rowprep.c:538
struct SCIP_NlhdlrData SCIP_NLHDLRDATA
Definition: type_nlhdlr.h:403
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17976
SCIP_VAR * SCIPgetExprAuxVarNonlinear(SCIP_EXPR *expr)
private functions of nonlinear handlers of nonlinear constraints
SCIPallocBlockMemory(scip, subsol))
void SCIProwprepSetLocal(SCIP_ROWPREP *rowprep, SCIP_Bool islocal)
Definition: misc_rowprep.c:748
SCIP_RETCODE SCIPaddRowprepTerms(SCIP *scip, SCIP_ROWPREP *rowprep, int nvars, SCIP_VAR **vars, SCIP_Real *coefs)
Definition: misc_rowprep.c:906
SCIP_Longint SCIPgetNLPs(SCIP *scip)
SCIP_INTERVAL SCIPgetExprBoundsNonlinear(SCIP *scip, SCIP_EXPR *expr)
#define SCIP_NLHDLR_METHOD_SEPABOTH
Definition: type_nlhdlr.h:44
SCIP_Real SCIPgetSolVal(SCIP *scip, SCIP_SOL *sol, SCIP_VAR *var)
Definition: scip_sol.c:1352
#define NLHDLR_NAME
preparation of a linear inequality to become a SCIP_ROW
SCIP_RETCODE SCIPsetPtrarrayVal(SCIP *scip, SCIP_PTRARRAY *ptrarray, int idx, void *val)
static SCIP_DECL_NLHDLRREVERSEPROP(nlhdlrReversepropQuotient)
void SCIPnlhdlrSetCopyHdlr(SCIP_NLHDLR *nlhdlr, SCIP_DECL_NLHDLRCOPYHDLR((*copy)))
Definition: nlhdlr.c:63