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