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 ≥ 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
70 SCIP_RETCODE exprdataCreate(
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
109 SCIP_RETCODE exprdataFree(
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
131 void transformExpr(
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
170 SCIP_RETCODE detectExpr(
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 */
319 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, xexpr,
320 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0,
321 xexpr == yexpr,
322 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0,
323 SCIPgetExprNAuxvarUsesNonlinear(expr) > 0) );
324
325 if( xexpr != yexpr && SCIPgetExprNAuxvarUsesNonlinear(expr) > 0 )
326 {
327 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, yexpr, TRUE, FALSE, TRUE, TRUE) );
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
345 SCIP_INTERVAL intEvalQuotient(
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 */
363 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) )
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 {
376 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result);
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("interval, a);
401 SCIPintervalDivScalar(SCIP_INTERVAL_INFINITY, "interval, 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, "interval, numinterval, denominterval);
417 SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, "interval, 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
429 SCIP_INTERVAL reversepropQuotient(
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 */
444 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bnds) )
445 return result;
446
447 /* substract constant from bounds of the expression */
448 SCIPintervalSubScalar(SCIP_INTERVAL_INFINITY, &bnds, bnds, e);
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 {
453 SCIPintervalSetEntire(SCIP_INTERVAL_INFINITY, &result);
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, "ient, 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
491 SCIP_RETCODE createRowprep(
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
532 SCIP_RETCODE estimateUnivariate(
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
624 SCIP_RETCODE estimateUnivariateQuotient(
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 */
657 SCIPintervalSetBounds(&bnd,
658 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(x)),
659 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(x)));
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
703 void hcGradCut(
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 ≤/≥ 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 ≤ 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 ≤ z, which results in a valid underestimator for x / y over the set
740 * {(x,y) | lbz ≤ x / y ≤ ubz}. Note that overestimating/underestimating the bilinear term with McCormick
741 * results in an underestimator/overestimator for x / y.
742 *
743 * 2. lbx ≥ 0 or ubx ≤ 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
752 SCIP_RETCODE estimateBivariate(
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
911 SCIP_RETCODE estimateBivariateQuotient(
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 */
956 SCIPintervalSetBounds(&bnd,
957 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[0])),
958 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(vars[0])));
959 SCIP_CALL( SCIPevalExprActivity(scip, xexpr) );
960 SCIPintervalIntersectEps(&bnd, SCIPepsilon(scip), bnd, SCIPexprGetActivity(xexpr));
961 lbx = bnd.inf;
962 ubx = bnd.sup;
963
964 SCIPintervalSetBounds(&bnd,
965 -infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, -SCIPvarGetLbGlobal(vars[1])),
966 infty2infty(SCIPinfinity(scip), SCIP_INTERVAL_INFINITY, SCIPvarGetUbGlobal(vars[1])));
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 */
991 SCIP_CALL( estimateBivariate(scip,
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
1135 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) );
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 */
1260 SCIP_RETCODE SCIPincludeNlhdlrQuotient(
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
|
(1) Event alloc_arg: |
"SCIPincludeNlhdlrNonlinear" allocates memory that is stored into "nlhdlr". [details] |
|
(2) Event cond_true: |
Condition "(_restat_ = SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, "quotient", "nonlinear handler for quotient expressions", 20, 20, nlhdlrDetectQuotient, nlhdlrEvalauxQuotient, nlhdlrdata)) != SCIP_OKAY", taking true branch. |
|
(3) Event leaked_storage: |
Variable "nlhdlr" going out of scope leaks the storage it points to. |
1272 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, NLHDLR_NAME, NLHDLR_DESC, NLHDLR_DETECTPRIORITY,
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 }
1283