SCIP Doxygen Documentation
 
Loading...
Searching...
No Matches
expr_pow.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 expr_pow.c
26 * @ingroup DEFPLUGINS_EXPR
27 * @brief power expression handler
28 * @author Benjamin Mueller
29 * @author Ksenia Bestuzheva
30 *
31 * @todo signpower for exponent < 1 ?
32 */
33
34/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
35
36/*lint --e{835}*/
37/*lint -e777*/
38
39#include <string.h>
40
41#include "scip/expr_pow.h"
42#include "scip/pub_expr.h"
43#include "scip/expr_value.h"
44#include "scip/expr_product.h"
45#include "scip/expr_sum.h"
46#include "scip/expr_exp.h"
47#include "scip/expr_abs.h"
48
49#define POWEXPRHDLR_NAME "pow"
50#define POWEXPRHDLR_DESC "power expression"
51#define POWEXPRHDLR_PRECEDENCE 55000
52#define POWEXPRHDLR_HASHKEY SCIPcalcFibHash(21163.0)
53
54#define SIGNPOWEXPRHDLR_NAME "signpower"
55#define SIGNPOWEXPRHDLR_DESC "signed power expression"
56#define SIGNPOWEXPRHDLR_PRECEDENCE 56000
57#define SIGNPOWEXPRHDLR_HASHKEY SCIPcalcFibHash(21163.1)
58
59#define INITLPMAXPOWVAL 1e+06 /**< maximal allowed absolute value of power expression at bound,
60 * used for adjusting bounds in the convex case in initestimates */
61
62/*
63 * Data structures
64 */
65
66/** sign of a value (-1 or +1)
67 *
68 * 0.0 has sign +1 here (shouldn't matter, though)
69 */
70#define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
71
72#define SIGNPOW_ROOTS_KNOWN 10 /**< up to which (integer) exponents precomputed roots have been stored */
73
74/** The positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 is needed in separation.
75 * Here we store these roots for small integer values of n.
76 */
77static
79 -1.0, /* no root for n=0 */
80 -1.0, /* no root for n=1 */
81 0.41421356237309504880, /* root for n=2 (-1+sqrt(2)) */
82 0.5, /* root for n=3 */
83 0.56042566045031785945, /* root for n=4 */
84 0.60582958618826802099, /* root for n=5 */
85 0.64146546982884663257, /* root for n=6 */
86 0.67033204760309682774, /* root for n=7 */
87 0.69428385661425826738, /* root for n=8 */
88 0.71453772716733489700, /* root for n=9 */
89 0.73192937842370733350 /* root for n=10 */
90};
91
92/** expression handler data */
93struct SCIP_ExprhdlrData
94{
95 SCIP_Real minzerodistance; /**< minimal distance from zero to enforce for child in bound tightening */
96 SCIP_Bool warnedonpole; /**< whether we warned on enforcing a minimal distance from zero for child */
97};
98
99/** expression data */
100struct SCIP_ExprData
101{
102 SCIP_Real exponent; /**< exponent */
103 SCIP_Real root; /**< positive root of (n-1) y^n + n y^(n-1) - 1, or
104 SCIP_INVALID if not computed yet */
105};
106
107/*
108 * Local methods
109 */
110
111/** computes positive root of the polynomial (n-1) y^n + n y^(n-1) - 1 for n > 1 */
112static
114 SCIP* scip, /**< SCIP data structure */
115 SCIP_Real* root, /**< buffer where to store computed root */
116 SCIP_Real exponent /**< exponent n */
117 )
118{
119 SCIP_Real polyval;
120 SCIP_Real gradval;
121 int iter;
122
123 assert(scip != NULL);
124 assert(exponent > 1.0);
125 assert(root != NULL);
126
127 /* lookup for popular integer exponent */
128 if( SCIPisIntegral(scip, exponent) && exponent-0.5 < SIGNPOW_ROOTS_KNOWN )
129 {
130 *root = signpow_roots[(int)SCIPfloor(scip, exponent+0.5)];
131 return SCIP_OKAY;
132 }
133
134 /* lookup for weymouth exponent */
135 if( SCIPisEQ(scip, exponent, 1.852) )
136 {
137 *root = 0.39821689389382575186;
138 return SCIP_OKAY;
139 }
140
141 /* search for a positive root of (n-1) y^n + n y^(n-1) - 1
142 * use the closest precomputed root as starting value
143 */
144 if( exponent >= SIGNPOW_ROOTS_KNOWN )
146 else if( exponent <= 2.0 )
147 *root = signpow_roots[2];
148 else
149 *root = signpow_roots[(int)SCIPfloor(scip, exponent)];
150
151 for( iter = 0; iter < 1000; ++iter )
152 {
153 polyval = (exponent - 1.0) * pow(*root, exponent) + exponent * pow(*root, exponent - 1.0) - 1.0;
154 if( fabs(polyval) < 1e-12 && SCIPisZero(scip, polyval) )
155 break;
156
157 /* gradient of (n-1) y^n + n y^(n-1) - 1 is n(n-1)y^(n-1) + n(n-1)y^(n-2) */
158 gradval = (exponent - 1.0) * exponent * (pow(*root, exponent - 1.0) + pow(*root, exponent - 2.0));
159 if( SCIPisZero(scip, gradval) )
160 break;
161
162 /* update root by adding -polyval/gradval (Newton's method) */
163 *root -= polyval / gradval;
164 if( *root < 0.0 )
165 *root = 0.0;
166 }
167
168 if( !SCIPisZero(scip, polyval) )
169 {
170 SCIPerrorMessage("failed to compute root for exponent %g\n", exponent);
171 return SCIP_ERROR;
172 }
173 SCIPdebugMsg(scip, "root for %g is %.20g, certainty = %g\n", exponent, *root, polyval);
174 /* @todo cache root value for other expressions (an exponent seldom comes alone)?? (they are actually really fast to compute...) */
175
176 return SCIP_OKAY;
177}
178
179/** computes negative root of the polynomial (n-1) y^n - n y^(n-1) + 1 for n < -1 */
180static
182 SCIP* scip, /**< SCIP data structure */
183 SCIP_Real* root, /**< buffer where to store computed root */
184 SCIP_Real exponent /**< exponent n */
185 )
186{
187 SCIP_Real polyval;
188 SCIP_Real gradval;
189 int iter;
190
191 assert(scip != NULL);
192 assert(exponent < -1.0);
193 assert(root != NULL);
194
195 *root = -2.0; /* that's the solution for n=-2 */
196
197 for( iter = 0; iter < 1000; ++iter )
198 {
199 polyval = (exponent - 1.0) * pow(*root, exponent) - exponent * pow(*root, exponent - 1.0) + 1.0;
200 if( fabs(polyval) < 1e-12 && SCIPisZero(scip, polyval) )
201 break;
202
203 /* gradient of (n-1) y^n - n y^(n-1) + 1 is n(n-1)y^(n-1) - n(n-1)y^(n-2) */
204 gradval = (exponent - 1.0) * exponent * (pow(*root, exponent - 1.0) - pow(*root, exponent - 2.0));
205 if( SCIPisZero(scip, gradval) )
206 break;
207
208 /* update root by adding -polyval/gradval (Newton's method) */
209 *root -= polyval / gradval;
210 if( *root >= 0.0 )
211 *root = -1;
212 }
213
214 if( !SCIPisZero(scip, polyval) )
215 {
216 SCIPerrorMessage("failed to compute root for exponent %g\n", exponent);
217 return SCIP_ERROR;
218 }
219 SCIPdebugMsg(scip, "root for %g is %.20g, certainty = %g\n", exponent, *root, polyval);
220 /* @todo cache root value for other expressions (an exponent seldom comes alone)?? (they are actually really fast to compute...) */
221
222 return SCIP_OKAY;
223}
224
225/** creates expression data */
226static
228 SCIP* scip, /**< SCIP data structure */
229 SCIP_EXPRDATA** exprdata, /**< pointer where to store expression data */
230 SCIP_Real exponent /**< exponent of the power expression */
231 )
232{
233 assert(exprdata != NULL);
234
235 SCIP_CALL( SCIPallocBlockMemory(scip, exprdata) );
236
237 (*exprdata)->exponent = exponent;
238 (*exprdata)->root = SCIP_INVALID;
239
240 return SCIP_OKAY;
241}
242
243/** computes a tangent at a reference point by linearization
244 *
245 * for a normal power, linearization in xref is xref^exponent + exponent * xref^(exponent-1) (x - xref)
246 * = (1-exponent) * xref^exponent + exponent * xref^(exponent-1) * x
247 *
248 * for a signpower, linearization is the same if xref is positive
249 * for xref negative it is -(-xref)^exponent + exponent * (-xref)^(exponent-1) (x-xref)
250 * = (1-exponent) * (-xref)^(exponent-1) * xref + exponent * (-xref)^(exponent-1) * x
251 */
252static
254 SCIP* scip, /**< SCIP data structure */
255 SCIP_Bool signpower, /**< are we signpower or normal power */
256 SCIP_Real exponent, /**< exponent */
257 SCIP_Real xref, /**< reference point where to linearize */
258 SCIP_Real* constant, /**< buffer to store constant term of secant */
259 SCIP_Real* slope, /**< buffer to store slope of secant */
260 SCIP_Bool* success /**< buffer to store whether secant could be computed */
261 )
262{
263 SCIP_Real xrefpow;
264
265 assert(scip != NULL);
266 assert(constant != NULL);
267 assert(slope != NULL);
268 assert(success != NULL);
269 assert(xref != 0.0 || exponent > 0.0);
270 /* non-integral exponent -> reference point must be >= 0 or we do signpower */
271 assert(EPSISINT(exponent, 0.0) || signpower || !SCIPisNegative(scip, xref));
272
273 /* TODO power is not differentiable at 0.0 for exponent < 0
274 * should we forbid here that xref > 0, do something smart here, or just return success=FALSE?
275 */
276 /* assert(exponent >= 1.0 || xref > 0.0); */
277
278 if( !EPSISINT(exponent, 0.0) && !signpower && xref < 0.0 )
279 xref = 0.0;
280
281 xrefpow = pow(signpower ? REALABS(xref) : xref, exponent - 1.0);
282
283 /* if huge xref and/or exponent too large, then pow may overflow */
284 if( !SCIPisFinite(xrefpow) )
285 {
286 *success = FALSE;
287 return;
288 }
289
290 *constant = (1.0 - exponent) * xrefpow * xref;
291 *slope = exponent * xrefpow;
292 *success = TRUE;
293}
294
295/** computes a secant between lower and upper bound
296 *
297 * secant is xlb^exponent + (xub^exponent - xlb^exponent) / (xub - xlb) * (x - xlb)
298 * = xlb^exponent - slope * xlb + slope * x with slope = (xub^exponent - xlb^exponent) / (xub - xlb)
299 * same if signpower
300 */
301static
303 SCIP* scip, /**< SCIP data structure */
304 SCIP_Bool signpower, /**< are we signpower or normal power */
305 SCIP_Real exponent, /**< exponent */
306 SCIP_Real xlb, /**< lower bound on x */
307 SCIP_Real xub, /**< upper bound on x */
308 SCIP_Real* constant, /**< buffer to store constant term of secant */
309 SCIP_Real* slope, /**< buffer to store slope of secant */
310 SCIP_Bool* success /**< buffer to store whether secant could be computed */
311 )
312{
313 assert(scip != NULL);
314 assert(constant != NULL);
315 assert(slope != NULL);
316 assert(success != NULL);
317 assert(xlb >= 0.0 || EPSISINT(exponent, 0.0) || signpower);
318 assert(xub >= 0.0 || EPSISINT(exponent, 0.0) || signpower);
319 assert(exponent != 1.0);
320
321 *success = FALSE;
322
323 /* infinite bounds will not work */
325 return;
326
327 /* first handle some special cases */
328 if( xlb == xub )
329 {
330 /* usually taken care of in separatePointPow already, but we might be called with different bounds here,
331 * e.g., when handling odd or signed power
332 */
333 *slope = 0.0;
334 *constant = pow(xlb, exponent);
335 }
336 else if( EPSISINT(exponent / 2.0, 0.0) && !signpower && xub > 0.1 && SCIPisFeasEQ(scip, xlb, -xub) )
337 {
338 /* for normal power with even exponents with xlb ~ -xub the slope would be very close to 0
339 * since xub^n - xlb^n is prone to cancellation here, we omit computing this secant (it's probably useless)
340 * unless the bounds are close to 0 as well (xub <= 0.1 in the "if" above)
341 * or we have exactly xlb=-xub, where we can return a clean 0.0 (though it's probably useless)
342 */
343 if( xlb == -xub )
344 {
345 *slope = 0.0;
346 *constant = pow(xlb, exponent);
347 }
348 else
349 {
350 return;
351 }
352 }
353 else if( xlb == 0.0 && exponent > 0.0 )
354 {
355 assert(xub >= 0.0);
356 *slope = pow(xub, exponent-1.0);
357 *constant = 0.0;
358 }
359 else if( xub == 0.0 && exponent > 0.0 )
360 {
361 /* normal pow: slope = - xlb^exponent / (-xlb) = xlb^(exponent-1)
362 * signpower: slope = (-xlb)^exponent / (-xlb) = (-xlb)^(exponent-1)
363 */
364 assert(xlb <= 0.0); /* so signpower or exponent is integral */
365 if( signpower )
366 *slope = pow(-xlb, exponent-1.0);
367 else
368 *slope = pow(xlb, exponent-1.0);
369 *constant = 0.0;
370 }
371 else if( SCIPisEQ(scip, xlb, xub) && (!signpower || xlb >= 0.0 || xub <= 0.0) )
372 {
373 /* Computing the slope as (xub^n - xlb^n)/(xub-xlb) can lead to cancellation.
374 * To avoid this, we replace xub^n by a Taylor expansion of pow at xlb:
375 * xub^n = xlb^n + n xlb^(n-1) (xub-xlb) + 0.5 n*(n-1) xlb^(n-2) (xub-xlb)^2 + 1/6 n*(n-1)*(n-2) xi^(n-3) (xub-xlb)^3 for some xlb < xi < xub
376 * Dropping the last term, the slope is (with an error of O((xub-xlb)^2) = 1e-18)
377 * n*xlb^(n-1) + 0.5 n*(n-1) xlb^(n-2)*(xub-xlb)
378 * = n*xlb^(n-1) (1 - 0.5*(n-1)) + 0.5 n*(n-1) xlb^(n-2)*xub
379 * = 0.5*n*((3-n)*xlb^(n-1) + (n-1) xlb^(n-2)*xub)
380 *
381 * test n=2: 0.5*2*((3-2)*xlb + (2-1) 1*xub) = xlb + xub ok
382 * n=3: 0.5*3*((3-3)*xlb + (3-1) xlb*xub) = 3*xlb*xub ~ xlb^2 + xlb*xub + xub^2 ok
383 *
384 * The constant is
385 * xlb^n - 0.5*n*((3-n) xlb^(n-1) + (n-1) xlb^(n-2)*xub) * xlb
386 * = xlb^n - 0.5*n*(3-n) xlb^n - 0.5*n*(n-1) xlb^(n-1)*xub
387 * = (1-0.5*n*(3-n)) xlb^n - 0.5 n*(n-1) xlb^(n-1) xub
388 *
389 * test n=2: (1-0.5*2*(3-2)) xlb^2 - 0.5 2*(2-1) xlb xub = -xlb*xub
390 * old formula: xlb^2 - (xlb+xub) * xlb = -xlb*xub ok
391 *
392 * For signpower with xub <= 0, we can negate xlb and xub:
393 * slope: (sign(xub)|xub|^n - sign(xlb)*|xlb|^n) / (xub-xlb) = -((-xub)^n - (-xlb)^n) / (xub - xlb) = ((-xub)^n - (-xlb)^n) / (-xub - (-xlb))
394 * constant: sign(xlb)|xlb|^n + slope * (xub - xlb) = -((-xlb)^n - slope * (xub - xlb)) = -((-xlb)^n + slope * ((-xub) - (-xlb)))
395 */
396 SCIP_Real xlb_n; /* xlb^n */
397 SCIP_Real xlb_n1; /* xlb^(n-1) */
398 SCIP_Real xlb_n2; /* xlb^(n-2) */
399
400 if( signpower && xub <= 0.0 )
401 {
402 xlb *= -1.0;
403 xub *= -1.0;
404 }
405
406 xlb_n = pow(xlb, exponent);
407 xlb_n1 = pow(xlb, exponent - 1.0);
408 xlb_n2 = pow(xlb, exponent - 2.0);
409
410 *slope = 0.5*exponent * ((3.0-exponent) * xlb_n1 + (exponent-1.0) * xlb_n2 * xub);
411 *constant = (1.0 - 0.5*exponent*(3.0-exponent)) * xlb_n - 0.5*exponent*(exponent-1.0) * xlb_n1 * xub;
412
413 if( signpower && xub <= 0.0 )
414 *constant *= -1.0;
415 }
416 else
417 {
418 SCIP_Real lbval;
419 SCIP_Real ubval;
420
421 if( signpower )
422 lbval = SIGN(xlb) * pow(REALABS(xlb), exponent);
423 else
424 lbval = pow(xlb, exponent);
425 if( !SCIPisFinite(lbval) )
426 return;
427
428 if( signpower )
429 ubval = SIGN(xub) * pow(REALABS(xub), exponent);
430 else
431 ubval = pow(xub, exponent);
432 if( !SCIPisFinite(ubval) )
433 return;
434
435 /* we still can have bad numerics when xlb^exponent and xub^exponent are very close, but xlb and xub are not
436 * for now, only check that things did not cancel out completely
437 */
438 if( lbval == ubval )
439 return;
440
441 *slope = (ubval - lbval) / (xub - xlb);
442 *constant = lbval - *slope * xlb;
443 }
444
445 /* check whether we had overflows */
446 if( !SCIPisFinite(*slope) || !SCIPisFinite(*constant) )
447 return;
448
449 *success = TRUE;
450}
451
452/** Separation for parabola
453 *
454 * - even positive powers: x^2, x^4, x^6 with x arbitrary, or
455 * - positive powers > 1: x^1.5, x^2.5 with x >= 0
456 <pre>
457 100 +--------------------------------------------------------------------+
458 |* + + + *|
459 90 |** x**2 ********|
460 | * * |
461 80 |-+* *+-|
462 | ** ** |
463 70 |-+ * * +-|
464 | ** ** |
465 60 |-+ * * +-|
466 | ** ** |
467 50 |-+ * * +-|
468 | ** ** |
469 40 |-+ * * +-|
470 | ** ** |
471 30 |-+ ** ** +-|
472 | ** ** |
473 20 |-+ ** ** +-|
474 | *** *** |
475 10 |-+ *** *** +-|
476 | + ***** + ***** + |
477 0 +--------------------------------------------------------------------+
478 -10 -5 0 5 10
479 </pre>
480 */
481static
483 SCIP* scip, /**< SCIP data structure */
484 SCIP_Real exponent, /**< exponent */
485 SCIP_Bool overestimate, /**< should the power be overestimated? */
486 SCIP_Real xlb, /**< lower bound on x */
487 SCIP_Real xub, /**< upper bound on x */
488 SCIP_Real xref, /**< reference point (where to linearize) */
489 SCIP_Real* constant, /**< buffer to store constant term of estimator */
490 SCIP_Real* slope, /**< buffer to store slope of estimator */
491 SCIP_Bool* islocal, /**< buffer to store whether estimator only locally valid, that is,
492 * it depends on given bounds */
493 SCIP_Bool* success /**< buffer to store whether estimator could be computed */
494 )
495{
496 assert(scip != NULL);
497 assert(constant != NULL);
498 assert(slope != NULL);
499 assert(islocal != NULL);
500 assert(success != NULL);
501 assert((exponent >= 0.0 && EPSISINT(exponent/2.0, 0.0)) || (exponent > 1.0 && xlb >= 0.0));
502
503 if( !overestimate )
504 {
505 computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
506 *islocal = FALSE;
507 }
508 else
509 {
510 /* overestimation -> secant */
511 computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
512 *islocal = TRUE;
513 }
514}
515
516
517/** Separation for signpower
518 *
519 * - odd positive powers, x^3, x^5, x^7
520 * - sign(x)|x|^n for n > 1
521 * - lower bound on x is negative (otherwise one should use separation for parabola)
522 <pre>
523 100 +--------------------------------------------------------------------+
524 | + + + **|
525 | x*abs(x) ******* |
526 | ** |
527 | ** |
528 50 |-+ *** +-|
529 | *** |
530 | *** |
531 | ***** |
532 | ***** |
533 0 |-+ **************** +-|
534 | ***** |
535 | ***** |
536 | *** |
537 | *** |
538 -50 |-+ *** +-|
539 | ** |
540 | ** |
541 | ** |
542 |** + + + |
543 -100 +--------------------------------------------------------------------+
544 -10 -5 0 5 10
545 </pre>
546 */
547static
549 SCIP* scip, /**< SCIP data structure */
550 SCIP_Real exponent, /**< exponent */
551 SCIP_Real root, /**< positive root of the polynomial (n-1) y^n + n y^(n-1) - 1,
552 * if xubglobal > 0 */
553 SCIP_Bool overestimate, /**< should the power be overestimated? */
554 SCIP_Real xlb, /**< lower bound on x, assumed to be non-positive */
555 SCIP_Real xub, /**< upper bound on x */
556 SCIP_Real xref, /**< reference point (where to linearize) */
557 SCIP_Real xlbglobal, /**< global lower bound on x */
558 SCIP_Real xubglobal, /**< global upper bound on x */
559 SCIP_Real* constant, /**< buffer to store constant term of estimator */
560 SCIP_Real* slope, /**< buffer to store slope of estimator */
561 SCIP_Bool* islocal, /**< buffer to store whether estimator only locally valid, that is,
562 * it depends on given bounds */
563 SCIP_Bool* branchcand, /**< buffer to indicate whether estimator would improve by branching
564 * on it */
565 SCIP_Bool* success /**< buffer to store whether estimator could be computed */
566 )
567{
568 assert(scip != NULL);
569 assert(constant != NULL);
570 assert(slope != NULL);
571 assert(islocal != NULL);
572 assert(branchcand != NULL);
573 assert(*branchcand == TRUE); /* the default */
574 assert(success != NULL);
575 assert(exponent >= 1.0);
576 assert(xlb < 0.0); /* otherwise estimateParabola should have been called */
577 assert(xubglobal <= 0.0 || (root > 0.0 && root < 1.0));
578
579 *success = FALSE;
580
581 if( !SCIPisPositive(scip, xub) )
582 {
583 /* easy case */
584 if( !overestimate )
585 {
586 /* underestimator is secant */
587 computeSecant(scip, TRUE, exponent, xlb, xub, constant, slope, success);
588 *islocal = TRUE;
589 }
590 else
591 {
592 /* overestimator is tangent */
593
594 /* we must linearize left of 0 */
595 if( xref > 0.0 )
596 xref = 0.0;
597
598 computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
599
600 /* if global upper bound is > 0, then the tangent is only valid locally if the reference point is right of
601 * -root*xubglobal
602 */
604
605 /* tangent doesn't move after branching */
606 *branchcand = FALSE;
607 }
608 }
609 else
610 {
611 SCIP_Real c;
612
613 if( !overestimate )
614 {
615 /* compute the special point which decides between secant and tangent */
616 c = -xlb * root;
617
618 if( xref < c )
619 {
620 /* underestimator is secant between xlb and c */
621 computeSecant(scip, TRUE, exponent, xlb, c, constant, slope, success);
622 *islocal = TRUE;
623 }
624 else
625 {
626 /* underestimator is tangent */
627 computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
628
629 /* if reference point is left of -root*xlbglobal (c w.r.t. global bounds),
630 * then tangent is not valid w.r.t. global bounds
631 */
632 *islocal = xref < -root * xlbglobal;
633
634 /* tangent doesn't move after branching */
635 *branchcand = FALSE;
636 }
637 }
638 else
639 {
640 /* compute the special point which decides between secant and tangent */
641 c = -xub * root;
642
643 if( xref <= c )
644 {
645 /* overestimator is tangent */
646 computeTangent(scip, TRUE, exponent, xref, constant, slope, success);
647
648 /* if reference point is right of -root*xubglobal (c w.r.t. global bounds),
649 * then tangent is not valid w.r.t. global bounds
650 */
651 *islocal = xref > -root * xubglobal;
652
653 /* tangent doesn't move after branching */
654 *branchcand = FALSE;
655 }
656 else
657 {
658 /* overestimator is secant */
659 computeSecant(scip, TRUE, exponent, c, xub, constant, slope, success);
660 *islocal = TRUE;
661 }
662 }
663 }
664}
665
666/** Separation for positive hyperbola
667 *
668 * - x^-2, x^-4 with x arbitrary
669 * - x^-0.5, x^-1, x^-1.5, x^-3, x^-5 with x >= 0
670 <pre>
671 5 +----------------------------------------------------------------------+
672 | + * +* + |
673 | * * x**(-2) ******* |
674 4 |-+ * * +-|
675 | * * |
676 | * * |
677 | * * |
678 3 |-+ * * +-|
679 | * * |
680 | * * |
681 2 |-+ * * +-|
682 | * * |
683 | * * |
684 1 |-+ * * +-|
685 | * * |
686 | ** ** |
687 | ********** ********** |
688 0 |******************* *******************|
689 | |
690 | + + + |
691 -1 +----------------------------------------------------------------------+
692 -10 -5 0 5 10
693 </pre>
694 */
695static
697 SCIP* scip, /**< SCIP data structure */
698 SCIP_Real exponent, /**< exponent */
699 SCIP_Real root, /**< negative root of the polynomial (n-1) y^n - n y^(n-1) + 1,
700 * if x has mixed sign (w.r.t. global bounds?) and underestimating */
701 SCIP_Bool overestimate, /**< should the power be overestimated? */
702 SCIP_Real xlb, /**< lower bound on x */
703 SCIP_Real xub, /**< upper bound on x */
704 SCIP_Real xref, /**< reference point (where to linearize) */
705 SCIP_Real xlbglobal, /**< global lower bound on x */
706 SCIP_Real xubglobal, /**< global upper bound on x */
707 SCIP_Real* constant, /**< buffer to store constant term of estimator */
708 SCIP_Real* slope, /**< buffer to store slope of estimator */
709 SCIP_Bool* islocal, /**< buffer to store whether estimator only locally valid, that is,
710 * it depends on given bounds */
711 SCIP_Bool* branchcand, /**< buffer to indicate whether estimator would improve by branching
712 * on it */
713 SCIP_Bool* success /**< buffer to store whether estimator could be computed */
714 )
715{
716 assert(scip != NULL);
717 assert(constant != NULL);
718 assert(slope != NULL);
719 assert(islocal != NULL);
720 assert(branchcand != NULL);
721 assert(*branchcand == TRUE); /* the default */
722 assert(success != NULL);
723 assert(exponent < 0.0);
724 assert(EPSISINT(exponent/2.0, 0.0) || xlb >= 0.0);
725
726 *success = FALSE;
727
728 if( !overestimate )
729 {
730 if( xlb >= 0.0 || xub <= 0.0 )
731 {
732 /* underestimate and fixed sign -> tangent */
733
734 /* make sure xref has the same sign as xlb,xub */
735 if( xref < 0.0 && xlb >= 0.0 )
736 xref = xlb;
737 else if( xref > 0.0 && xub <= 0.0 )
738 xref = xub;
739
740 if( SCIPisZero(scip, xref) )
741 {
742 /* estimator would need to have an (essentially) infinite scope
743 * first try to make up a better refpoint
744 */
745 if( xub > 0.0 )
746 {
747 /* thus xlb >= 0.0; stay close to xlb (probably = 0) */
748 if( !SCIPisInfinity(scip, xub) )
749 xref = 0.9 * xlb + 0.1 * xub;
750 else
751 xref = 0.1;
752 }
753 else
754 {
755 /* xub <= 0.0; stay close to xub (probably = 0) */
756 if( !SCIPisInfinity(scip, -xlb) )
757 xref = 0.1 * xlb + 0.9 * xub;
758 else
759 xref = 0.1;
760 }
761
762 /* if still close to 0, then also bounds are close to 0, then just give up */
763 if( SCIPisZero(scip, xref) )
764 return;
765 }
766
767 computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
768
769 /* tangent will not change if branching on x (even if only locally valid, see checks below) */
770 *branchcand = FALSE;
771
772 if( EPSISINT(exponent/2.0, 0.0) )
773 {
774 /* for even exponents (as in the picture):
775 * if x has fixed sign globally, then our tangent is also globally valid
776 * however, if x has mixed sign, then it depends on the constellation between reference point and global
777 * bounds, whether the tangent is globally valid (see also the longer discussion for the mixed-sign
778 * underestimator below )
779 */
780 if( xref > 0.0 && xlbglobal < 0.0 )
781 {
782 assert(xubglobal > 0.0); /* since xref > 0.0 */
783 assert(root < 0.0); /* root needs to be given */
784 /* if on right side, then tangent is only locally valid if xref is too much to the left */
785 *islocal = xref < xlbglobal * root;
786 }
787 else if( xref < 0.0 && xubglobal > 0.0 )
788 {
789 assert(xlbglobal < 0.0); /* since xref < 0.0 */
790 assert(root < 0.0); /* root needs to be given */
791 /* if on left side, then tangent is only locally valid if xref is too much to the right */
792 *islocal = xref > xubglobal * root;
793 }
794 else
795 *islocal = FALSE;
796 }
797 else
798 {
799 /* for odd exponents, the tangent is only locally valid if the sign of x is not fixed globally */
800 *islocal = xlbglobal * xubglobal < 0.0;
801 }
802 }
803 else
804 {
805 /* underestimate but mixed sign */
806 if( SCIPisInfinity(scip, -xlb) )
807 {
808 if( SCIPisInfinity(scip, xub) )
809 {
810 /* underestimator is constant 0, but that is globally valid */
811 *constant = 0.0;
812 *slope = 0.0;
813 *islocal = FALSE;
814 *success = TRUE;
815 return;
816 }
817
818 /* switch sign of x (mirror on ordinate) to make left bound finite and use its estimator */
819 estimateHyperbolaPositive(scip, exponent, root, overestimate, -xub, -xlb, -xref, -xubglobal, -xlbglobal,
820 constant, slope, islocal, branchcand, success);
821 if( *success )
822 *slope = -*slope;
823 }
824 else
825 {
826 /* The convex envelope of x^exponent for x in [xlb, infinity] is a line (secant) between xlb and some positive
827 * coordinate xhat, and x^exponent for x > xhat.
828 * Further, on [xlb,xub] with xub < xhat, the convex envelope is the secant between xlb and xub.
829 *
830 * To find xhat, consider the affine-linear function l(x) = xlb^n + c * (x - xlb) where n = exponent
831 * we look for a value of x such that f(x) and l(x) coincide and such that l(x) will be tangent to f(x) on that
832 * point, that is
833 * xhat > 0 such that f(xhat) = l(xhat) and f'(xhat) = l'(xhat)
834 * => xhat^n = xlb^n + c * (xhat - xlb) and n * xhat^(n-1) = c
835 * => xhat^n = xlb^n + n * xhat^n - n * xhat^(n-1) * xlb
836 * => 0 = xlb^n + (n-1) * xhat^n - n * xhat^(n-1) * xlb
837 *
838 * Divide by xlb^n, one gets a polynomial that looks very much like the one for signpower, but a sign is
839 * different (since this is *not signed* power):
840 * 0 = 1 + (n-1) * y^n - n * y^(n-1) where y = xhat/xlb
841 *
842 * The solution y < 0 (because xlb < 0 and we want xhat > 0) is what we expect to be given as "root".
843 */
844 assert(root < 0.0); /* root needs to be given */
845 if( xref <= xlb * root )
846 {
847 /* If the reference point is left of xhat (=xlb*root), then we can take the
848 * secant between xlb and root*xlb (= tangent at root*xlb).
849 * However, if xub < root*xlb, then we can tilt the estimator to be the secant between xlb and xub.
850 */
851 computeSecant(scip, FALSE, exponent, xlb, MIN(xlb * root, xub), constant, slope, success);
852 *islocal = TRUE;
853 }
854 else
855 {
856 /* If reference point is right of xhat, then take the tangent at xref.
857 * This will still be underestimating for x in [xlb,0], too.
858 * The tangent is globally valid, if we had also generated w.r.t. global bounds.
859 */
860 computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
861 *islocal = xref < xlbglobal * root;
862 *branchcand = FALSE;
863 }
864 }
865 }
866 }
867 else
868 {
869 /* overestimate and mixed sign -> pole is within domain -> cannot overestimate */
870 if( xlb < 0.0 && xub > 0.0 )
871 return;
872
873 /* overestimate and fixed sign -> secant */
874 computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
875 *islocal = TRUE;
876 }
877
878}
879
880/** Separation for mixed-sign hyperbola
881 *
882 * - x^-1, x^-3, x^-5 without x >= 0 (either x arbitrary or x negative)
883 <pre>
884 +----------------------------------------------------------------------+
885 | + * + |
886 4 |-+ * x**(-1) *******-|
887 | * |
888 | * |
889 | * |
890 2 |-+ * +-|
891 | * |
892 | ** |
893 | ********* |
894 0 |********************* *********************|
895 | ********* |
896 | ** |
897 | * |
898 -2 |-+ * +-|
899 | * |
900 | * |
901 | * |
902 -4 |-+ * +-|
903 | + *+ + |
904 +----------------------------------------------------------------------+
905 -10 -5 0 5 10
906 </pre>
907 */
908static
910 SCIP* scip, /**< SCIP data structure */
911 SCIP_Real exponent, /**< exponent */
912 SCIP_Bool overestimate, /**< should the power be overestimated? */
913 SCIP_Real xlb, /**< lower bound on x */
914 SCIP_Real xub, /**< upper bound on x */
915 SCIP_Real xref, /**< reference point (where to linearize) */
916 SCIP_Real xlbglobal, /**< global lower bound on x */
917 SCIP_Real xubglobal, /**< global upper bound on x */
918 SCIP_Real* constant, /**< buffer to store constant term of estimator */
919 SCIP_Real* slope, /**< buffer to store slope of estimator */
920 SCIP_Bool* islocal, /**< buffer to store whether estimator only locally valid, that is,
921 it depends on given bounds */
922 SCIP_Bool* branchcand, /**< buffer to indicate whether estimator would improve by branching
923 on it */
924 SCIP_Bool* success /**< buffer to store whether estimator could be computed */
925 )
926{
927 assert(scip != NULL);
928 assert(constant != NULL);
929 assert(slope != NULL);
930 assert(islocal != NULL);
931 assert(branchcand != NULL);
932 assert(*branchcand == TRUE); /* the default */
933 assert(success != NULL);
934 assert(exponent < 0.0);
935 assert(EPSISINT((exponent-1.0)/2.0, 0.0));
936 assert(xlb < 0.0);
937
938 *success = FALSE;
939
940 if( xub <= 0.0 )
941 {
942 /* x is negative */
943 if( !overestimate )
944 {
945 /* underestimation -> secant */
946 computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
947 *islocal = TRUE;
948 }
949 else if( !SCIPisZero(scip, xlb/10.0) )
950 {
951 /* overestimation -> tangent */
952
953 /* need to linearize left of 0 */
954 if( xref > 0.0 )
955 xref = 0.0;
956
957 if( SCIPisZero(scip, xref) )
958 {
959 /* if xref is very close to 0.0, then slope would be infinite
960 * try to move closer to lower bound (if xlb < -10*eps)
961 */
962 if( !SCIPisInfinity(scip, -xlb) )
963 xref = 0.1*xlb + 0.9*xub;
964 else
965 xref = 0.1;
966 }
967
968 computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
969
970 /* if x does not have a fixed sign globally, then our tangent is not globally valid
971 * (power is not convex on global domain)
972 */
973 *islocal = xlbglobal * xubglobal < 0.0;
974
975 /* tangent doesn't move by branching */
976 *branchcand = FALSE;
977 }
978 /* else: xlb is very close to zero, xub is <= 0, so slope would be infinite
979 * (for any reference point in [xlb, xub]) -> do not estimate
980 */
981 }
982 /* else: x has mixed sign -> pole is within domain -> cannot estimate */
983}
984
985/** Separation for roots with exponent in [0,1]
986 *
987 * - x^0.5 with x >= 0
988 <pre>
989 8 +----------------------------------------------------------------------+
990 | + + + + |
991 7 |-+ x**0.5 ********|
992 | *********|
993 | ******** |
994 6 |-+ ******** +-|
995 | ****** |
996 5 |-+ ****** +-|
997 | ****** |
998 | ***** |
999 4 |-+ **** +-|
1000 | ***** |
1001 3 |-+ **** +-|
1002 | *** |
1003 | *** |
1004 2 |-+ ** +-|
1005 | ** |
1006 1 |** +-|
1007 |* |
1008 |* + + + + |
1009 0 +----------------------------------------------------------------------+
1010 0 10 20 30 40 50
1011 </pre>
1012 */
1013static
1015 SCIP* scip, /**< SCIP data structure */
1016 SCIP_Real exponent, /**< exponent */
1017 SCIP_Bool overestimate, /**< should the power be overestimated? */
1018 SCIP_Real xlb, /**< lower bound on x */
1019 SCIP_Real xub, /**< upper bound on x */
1020 SCIP_Real xref, /**< reference point (where to linearize) */
1021 SCIP_Real* constant, /**< buffer to store constant term of estimator */
1022 SCIP_Real* slope, /**< buffer to store slope of estimator */
1023 SCIP_Bool* islocal, /**< buffer to store whether estimator only locally valid, that is,
1024 * it depends on given bounds */
1025 SCIP_Bool* success /**< buffer to store whether estimator could be computed */
1026 )
1027{
1028 assert(scip != NULL);
1029 assert(constant != NULL);
1030 assert(slope != NULL);
1031 assert(islocal != NULL);
1032 assert(success != NULL);
1033 assert(exponent > 0.0);
1034 assert(exponent < 1.0);
1035 assert(xlb >= 0.0);
1036
1037 if( !overestimate )
1038 {
1039 /* underestimate -> secant */
1040 computeSecant(scip, FALSE, exponent, xlb, xub, constant, slope, success);
1041 *islocal = TRUE;
1042 }
1043 else
1044 {
1045 /* overestimate -> tangent */
1046
1047 /* need to linearize right of 0 */
1048 if( xref < 0.0 )
1049 xref = 0.0;
1050
1051 if( SCIPisZero(scip, xref) )
1052 {
1053 if( SCIPisZero(scip, xub) )
1054 {
1055 *success = FALSE;
1056 *islocal = FALSE;
1057 return;
1058 }
1059
1060 /* if xref is 0 (then xlb=0 probably), then slope is infinite, then try to move away from 0 */
1061 if( xub < 0.2 )
1062 xref = 0.5 * xlb + 0.5 * xub;
1063 else
1064 xref = 0.1;
1065 }
1066
1067 computeTangent(scip, FALSE, exponent, xref, constant, slope, success);
1068 *islocal = FALSE;
1069 }
1070}
1071
1072/** builds an estimator for a power function */
1073static
1075 SCIP* scip, /**< SCIP data structure */
1076 SCIP_EXPRDATA* exprdata, /**< expression data */
1077 SCIP_Bool overestimate, /**< is this an overestimator? */
1078 SCIP_Real childlb, /**< local lower bound on the child */
1079 SCIP_Real childub, /**< local upper bound on the child */
1080 SCIP_Real childglb, /**< global lower bound on the child */
1081 SCIP_Real childgub, /**< global upper bound on the child */
1082 SCIP_Bool childintegral, /**< whether child is integral */
1083 SCIP_Real refpoint, /**< reference point */
1084 SCIP_Real exponent, /**< esponent */
1085 SCIP_Real* coef, /**< pointer to store the coefficient of the estimator */
1086 SCIP_Real* constant, /**< pointer to store the constant of the estimator */
1087 SCIP_Bool* success, /**< pointer to store whether the estimator was built successfully */
1088 SCIP_Bool* islocal, /**< pointer to store whether the estimator is valid w.r.t. local bounds
1089 * only */
1090 SCIP_Bool* branchcand /**< pointer to indicate whether to consider child for branching
1091 * (initialized to TRUE) */
1092 )
1093{
1094 SCIP_Bool isinteger;
1095 SCIP_Bool iseven;
1096
1097 assert(scip != NULL);
1098 assert(exprdata != NULL);
1099 assert(coef != NULL);
1100 assert(constant != NULL);
1101 assert(success != NULL);
1102 assert(islocal != NULL);
1103 assert(branchcand != NULL);
1104
1105 isinteger = EPSISINT(exponent, 0.0);
1106 iseven = isinteger && EPSISINT(exponent / 2.0, 0.0);
1107
1108 if( exponent == 2.0 )
1109 {
1110 /* important special case: quadratic case */
1111 /* initialize, because SCIPaddSquareXyz only adds to existing values */
1112 *success = TRUE;
1113 *coef = 0.0;
1114 *constant = 0.0;
1115
1116 if( overestimate )
1117 {
1118 SCIPaddSquareSecant(scip, 1.0, childlb, childub, coef, constant, success);
1119 *islocal = TRUE; /* secants are only valid locally */
1120 }
1121 else
1122 {
1124 *islocal = FALSE; /* linearizations are globally valid */
1125 *branchcand = FALSE; /* there is no improvement due to branching */
1126 }
1127 }
1128 else if( exponent > 0.0 && iseven )
1129 {
1130 estimateParabola(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1131 /* if estimate is locally valid, then we computed a secant and so branching can improve it */
1132 *branchcand = *islocal;
1133 }
1134 else if( exponent > 1.0 && childlb >= 0.0 )
1135 {
1136 /* make sure we linearize in convex region (if we will linearize) */
1137 if( refpoint < 0.0 )
1138 refpoint = 0.0;
1139
1140 estimateParabola(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1141
1142 /* if estimate is locally valid, then we computed a secant and so branching can improve it */
1143 *branchcand = *islocal;
1144
1145 /* if odd power, then check whether tangent on parabola is also globally valid, that is reference point is
1146 * right of -root*global-lower-bound
1147 */
1148 if( !*islocal && !iseven && childglb < 0.0 )
1149 {
1151 *islocal = TRUE;
1152 else
1153 {
1154 if( exprdata->root == SCIP_INVALID )
1155 {
1156 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1157 }
1158 *islocal = refpoint < exprdata->root * (-childglb);
1159 }
1160 }
1161 }
1162 else if( exponent > 1.0 ) /* and !iseven && childlb < 0.0 due to previous if */
1163 {
1164 /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
1165 if( exprdata->root == SCIP_INVALID && childgub > 0.0 )
1166 {
1167 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1168 }
1169 estimateSignedpower(scip, exponent, exprdata->root, overestimate, childlb, childub, refpoint,
1170 -childglb, childgub, constant, coef, islocal, branchcand, success);
1171 }
1172 else if( exponent < 0.0 && (iseven || childlb >= 0.0) )
1173 {
1174 /* compute root if not known yet; only needed if mixed sign (globally) and iseven */
1175 if( exprdata->root == SCIP_INVALID && iseven )
1176 {
1177 SCIP_CALL( computeHyperbolaRoot(scip, &exprdata->root, exponent) );
1178 }
1179 estimateHyperbolaPositive(scip, exponent, exprdata->root, overestimate, childlb, childub, refpoint,
1180 childglb, childgub, constant, coef, islocal, branchcand, success);
1181 }
1182 else if( exponent < 0.0 )
1183 {
1184 assert(!iseven); /* should hold due to previous if */
1185 assert(childlb < 0.0); /* should hold due to previous if */
1186 assert(isinteger); /* should hold because childlb < 0.0 (same as assert above) */
1187
1189 constant, coef, islocal, branchcand, success);
1190 }
1191 else
1192 {
1193 assert(exponent < 1.0); /* the only case that should be left */
1194 assert(exponent > 0.0); /* should hold due to previous if */
1195
1196 estimateRoot(scip, exponent, overestimate, childlb, childub, refpoint, constant, coef, islocal, success);
1197
1198 /* if estimate is locally valid, then we computed a secant, and so branching can improve it */
1199 *branchcand = *islocal;
1200 }
1201
1202 return SCIP_OKAY;
1203}
1204
1205/** fills an array of reference points for estimating on the convex side */
1206static
1208 SCIP* scip, /**< SCIP data structure */
1209 SCIP_Real exponent, /**< exponent of the power expression */
1210 SCIP_Real lb, /**< lower bound on the child variable */
1211 SCIP_Real ub, /**< upper bound on the child variable */
1212 SCIP_Real* refpoints /**< array to store the reference points */
1213 )
1214{
1215 SCIP_Real maxabsbnd;
1216
1217 assert(refpoints != NULL);
1218
1219 maxabsbnd = pow(INITLPMAXPOWVAL, 1 / exponent);
1220
1221 /* make sure the absolute values of bounds are not too large */
1222 if( ub > -maxabsbnd )
1223 lb = MAX(lb, -maxabsbnd);
1224 if( lb < maxabsbnd )
1225 ub = MIN(ub, maxabsbnd);
1226
1227 /* in the case when ub < -maxabsbnd or lb > maxabsbnd, we still want to at least make bounds finite */
1228 if( SCIPisInfinity(scip, -lb) )
1229 lb = MIN(-10.0, ub - 0.1*REALABS(ub)); /*lint !e666 */
1230 if( SCIPisInfinity(scip, ub) )
1231 ub = MAX( 10.0, lb + 0.1*REALABS(lb)); /*lint !e666 */
1232
1233 refpoints[0] = (7.0 * lb + ub) / 8.0;
1234 refpoints[1] = (lb + ub) / 2.0;
1235 refpoints[2] = (lb + 7.0 * ub) / 8.0;
1236}
1237
1238/** fills an array of reference points for sign(x)*abs(x)^n or x^n (n odd), where x has mixed signs
1239 *
1240 * The reference points are: the lower and upper bounds (one for secant and one for tangent);
1241 * and for the second tangent, the point on the convex part of the function between the point
1242 * deciding between tangent and secant, and the corresponding bound
1243 */
1244static
1246 SCIP* scip, /**< SCIP data structure */
1247 SCIP_EXPRDATA* exprdata, /**< expression data */
1248 SCIP_Real lb, /**< lower bound on the child variable */
1249 SCIP_Real ub, /**< upper bound on the child variable */
1250 SCIP_Real exponent, /**< exponent */
1251 SCIP_Bool underestimate, /**< are the refpoints for an underestimator */
1252 SCIP_Real* refpoints /**< array to store the reference points */
1253 )
1254{
1255 assert(refpoints != NULL);
1256
1257 if( (underestimate && SCIPisInfinity(scip, -lb)) || (!underestimate && SCIPisInfinity(scip, ub)) )
1258 return SCIP_OKAY;
1259
1260 if( exprdata->root == SCIP_INVALID )
1261 {
1262 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
1263 }
1264
1265 /* make bounds finite (due to a previous if, only one can be infinite here) */
1266 if( SCIPisInfinity(scip, -lb) )
1267 lb = -ub * exprdata->root - 1.0;
1268 else if( SCIPisInfinity(scip, ub) )
1269 ub = -lb * exprdata->root + 1.0;
1270
1271 if( underestimate )
1272 {
1273 /* secant point */
1274 refpoints[0] = lb;
1275
1276 /* tangent points, depending on the special point */
1277 if( -lb * exprdata->root < ub - 2.0 )
1278 refpoints[2] = ub;
1279 if( -lb * exprdata->root < ub - 4.0 )
1280 refpoints[1] = (-lb * exprdata->root + ub) / 2.0;
1281 }
1282
1283 if( !underestimate )
1284 {
1285 /* secant point */
1286 refpoints[2] = ub;
1287
1288 /* tangent points, depending on the special point */
1289 if( -ub * exprdata->root > lb + 2.0 )
1290 refpoints[0] = lb;
1291 if( -ub * exprdata->root > lb + 4.0 )
1292 refpoints[1] = (lb - ub * exprdata->root) / 2.0;
1293 }
1294
1295 return SCIP_OKAY;
1296}
1297
1298/** choose reference points for adding initestimates cuts for a power expression */
1299static
1301 SCIP* scip, /**< SCIP data structure */
1302 SCIP_EXPRDATA* exprdata, /**< expression data */
1303 SCIP_Real lb, /**< lower bound on the child variable */
1304 SCIP_Real ub, /**< upper bound on the child variable */
1305 SCIP_Real* refpointsunder, /**< array to store reference points for underestimators */
1306 SCIP_Real* refpointsover, /**< array to store reference points for overestimators */
1307 SCIP_Bool underestimate, /**< whether refpoints for underestimation are needed */
1308 SCIP_Bool overestimate /**< whether refpoints for overestimation are needed */
1309 )
1310{
1311 SCIP_Bool convex;
1312 SCIP_Bool concave;
1313 SCIP_Bool mixedsign;
1314 SCIP_Bool even;
1315 SCIP_Real exponent;
1316
1317 assert(scip != NULL);
1318 assert(exprdata != NULL);
1320
1321 exponent = exprdata->exponent;
1322 even = EPSISINT(exponent, 0.0) && EPSISINT(exponent / 2.0, 0.0);
1323
1324 convex = FALSE;
1325 concave = FALSE;
1326 mixedsign = lb < 0.0 && ub > 0.0;
1327
1328 /* convex case:
1329 * - parabola with an even degree or positive domain
1330 * - hyperbola with a positive domain
1331 * - even hyperbola with a negative domain
1332 */
1333 if( (exponent > 1.0 && (lb >= 0 || even)) || (exponent < 0.0 && lb >= 0) || (exponent < 0.0 && even && ub <= 0.0) )
1334 convex = TRUE;
1335 /* concave case:
1336 * - parabola or hyperbola with a negative domain and (due to previous if) an uneven degree
1337 * - root
1338 */
1339 else if( ub <= 0 || (exponent > 0.0 && exponent < 1.0) )
1340 concave = TRUE;
1341
1342 if( underestimate )
1343 {
1344 if( convex )
1345 addTangentRefpoints(scip, exponent, lb, ub, refpointsunder);
1346 else if( (concave && !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub)) ||
1347 (exponent < 0.0 && even && mixedsign) ) /* concave with finite bounds or mixed even hyperbola */
1348 {
1349 /* for secant, refpoint doesn't matter, but we add it to signal that the corresponding cut should be created */
1350 refpointsunder[0] = (lb + ub) / 2.0;
1351 }
1352 else if( exponent > 1.0 && !even && mixedsign ) /* mixed signpower */
1353 {
1354 SCIP_CALL( addSignpowerRefpoints(scip, exprdata, lb, ub, exponent, TRUE, refpointsunder) );
1355 }
1356 else /* mixed odd hyperbola or an infinite bound */
1357 assert((exponent < 0.0 && !even && mixedsign) || SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub));
1358 }
1359
1360 if( overestimate )
1361 {
1362 if( convex && !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) )
1363 refpointsover[0] = (lb + ub) / 2.0;
1364 else if( concave )
1365 addTangentRefpoints(scip, exponent, lb, ub, refpointsover);
1366 else if( exponent > 1.0 && !even && mixedsign ) /* mixed signpower */
1367 {
1368 SCIP_CALL( addSignpowerRefpoints(scip, exprdata, lb, ub, exponent, FALSE, refpointsover) );
1369 }
1370 else /* mixed hyperbola or an infinite bound */
1371 assert((exponent < 0.0 && mixedsign) || SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub));
1372 }
1373
1374 return SCIP_OKAY;
1375}
1376
1377/*
1378 * Callback methods of expression handler
1379 */
1380
1381/** compares two power expressions
1382 *
1383 * the order of two power (normal or signed) is base_1^expo_1 < base_2^expo_2 if and only if
1384 * base_1 < base2 or, base_1 = base_2 and expo_1 < expo_2
1385 */
1386static
1388{ /*lint --e{715}*/
1389 SCIP_Real expo1;
1390 SCIP_Real expo2;
1391 int compareresult;
1392
1393 /**! [SnippetExprComparePow] */
1395 if( compareresult != 0 )
1396 return compareresult;
1397
1400
1401 return expo1 == expo2 ? 0 : expo1 < expo2 ? -1 : 1;
1402 /**! [SnippetExprComparePow] */
1403}
1404
1405/** simplifies a pow expression
1406 *
1407 * Evaluates the power function when its child is a value expression
1408 */
1409static
1411{ /*lint --e{715}*/
1412 SCIP_EXPR* base;
1413 SCIP_Real exponent;
1414
1415 assert(scip != NULL);
1416 assert(expr != NULL);
1418 assert(SCIPexprGetNChildren(expr) == 1);
1419
1420 base = SCIPexprGetChildren(expr)[0];
1421 assert(base != NULL);
1422
1423 exponent = SCIPgetExponentExprPow(expr);
1424
1425 SCIPdebugPrintf("[simplifyPow] simplifying power with expo %g\n", exponent);
1426
1427 /* enforces POW1 */
1428 if( exponent == 0.0 )
1429 {
1430 SCIPdebugPrintf("[simplifyPow] POW1\n");
1431 /* TODO: more checks? */
1434 return SCIP_OKAY;
1435 }
1436
1437 /* enforces POW2 */
1438 if( exponent == 1.0 )
1439 {
1440 SCIPdebugPrintf("[simplifyPow] POW2\n");
1443 return SCIP_OKAY;
1444 }
1445
1446 /* enforces POW3 */
1447 if( SCIPisExprValue(scip, base) )
1448 {
1449 SCIP_Real baseval;
1450
1451 SCIPdebugPrintf("[simplifyPow] POW3\n");
1453
1454 /* the assert below was failing on st_e35 for baseval=-1e-15 and fractional exponent
1455 * in the subNLP heuristic; I assume that this was because baseval was evaluated after
1456 * variable fixings and that there were just floating-point inaccuracies and 0 was meant,
1457 * so I treat -1e-15 as 0 here
1458 */
1459 if( baseval < 0.0 && fmod(exponent, 1.0) != 0.0 && baseval > -SCIPepsilon(scip) )
1460 baseval = 0.0;
1461
1462 /* TODO check if those are all important asserts */
1463 assert(baseval >= 0.0 || fmod(exponent, 1.0) == 0.0);
1464 assert(baseval != 0.0 || exponent != 0.0);
1465
1466 if( baseval != 0.0 || exponent > 0.0 )
1467 {
1469 return SCIP_OKAY;
1470 }
1471 }
1472
1473 /* enforces POW11 (exp(x)^n = exp(n*x)) */
1474 if( SCIPisExprExp(scip, base) )
1475 {
1476 SCIP_EXPR* child;
1477 SCIP_EXPR* prod;
1480
1481 SCIPdebugPrintf("[simplifyPow] POW11\n");
1482 child = SCIPexprGetChildren(base)[0];
1483
1484 /* multiply child of exponential with exponent */
1486
1487 /* simplify product */
1490
1491 /* create exponential with new child */
1494
1495 /* the final simplified expression is the simplification of the just created exponential */
1498
1499 return SCIP_OKAY;
1500 }
1501
1502 /* enforces POW10 */
1503 if( SCIPisExprVar(scip, base) )
1504 {
1506
1507 SCIPdebugPrintf("[simplifyPow] POW10\n");
1509
1510 assert(basevar != NULL);
1511
1512 /* TODO: if exponent is negative, we could fix the binary variable to 1. However, this is a bit tricky because
1513 * variables can not be tighten in EXITPRE, where the simplify is also called
1514 */
1515 if( SCIPvarIsBinary(basevar) && exponent > 0.0 )
1516 {
1519 return SCIP_OKAY;
1520 }
1521 }
1522
1523 if( EPSISINT(exponent, 0.0) )
1524 {
1525 SCIP_EXPR* aux;
1527
1528 /* enforces POW5
1529 * given (pow n (prod 1.0 expr_1 ... expr_k)) we distribute the exponent:
1530 * -> (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1531 * notes: - since base is simplified and its coefficient is 1.0 (SP8)
1532 * - n is an integer (excluding 1 and 0; see POW1-2 above)
1533 */
1535 {
1537 int i;
1538
1539 /* create empty product */
1541
1542 for( i = 0; i < SCIPexprGetNChildren(base); ++i )
1543 {
1544 /* create (pow n expr_i) and simplify */
1547 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1548
1549 /* append (pow n expr_i) to product */
1552 }
1553
1554 /* simplify (prod 1.0 (pow n expr_1) ... (pow n expr_k))
1555 * this calls simplifyProduct directly, since we know its children are simplified */
1558 return SCIP_OKAY;
1559 }
1560
1561 /* enforces POW6
1562 * given (pow n (sum 0.0 coef expr)) we can move `pow` inside `sum`:
1563 * (pow n (sum 0.0 coef expr) ) -> (sum 0.0 coef^n (pow n expr))
1564 * notes: - since base is simplified and its constant is 0, then coef != 1.0 (SS7)
1565 * - n is an integer (excluding 1 and 0; see POW1-2 above)
1566 */
1568 {
1569 SCIP_Real newcoef;
1570
1571 SCIPdebugPrintf("[simplifyPow] seeing a sum with one term, exponent %g\n", exponent);
1572
1573 /* assert SS7 holds */
1574 assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
1575
1576 /* create (pow n expr) and simplify it
1577 * note: we call simplifyPow directly, since we know that `expr` is simplified */
1578 newcoef = pow(SCIPgetCoefsExprSum(base)[0], exponent);
1581 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1582
1583 /* create (sum (pow n expr)) and simplify it
1584 * this calls simplifySum directly, since we know its children are simplified */
1587 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1589 return SCIP_OKAY;
1590 }
1591
1592 /* enforces POW7
1593 * (const + sum alpha_i expr_i)^2 = sum alpha_i^2 expr_i^2
1594 * + sum_{j < i} 2*alpha_i alpha_j expr_i expr_j
1595 * + sum const alpha_i expr_i
1596 * TODO: put some limits on the number of children of the sum being expanded
1597 */
1598 if( SCIPisExprSum(scip, base) && exponent == 2.0 )
1599 {
1600 int i;
1601 int nchildren;
1605 SCIP_Real* coefs;
1606 SCIP_Real constant;
1607
1608 SCIPdebugPrintf("[simplifyPow] expanding sum^%g\n", exponent);
1609
1610 nchildren = SCIPexprGetNChildren(base);
1611 nexpandedchildren = nchildren * (nchildren + 1) / 2 + nchildren;
1614
1615 for( i = 0; i < nchildren; ++i )
1616 {
1617 int j;
1621
1622 /* create and simplify expr_i * expr_j */
1623 for( j = 0; j < i; ++j )
1624 {
1626 coefs[i*(i+1)/2 + j] = 2 * SCIPgetCoefsExprSum(base)[i] * SCIPgetCoefsExprSum(base)[j];
1627
1629 ownercreatedata) );
1631 ownercreate, ownercreatedata) ); /* this calls simplifyProduct */
1633 }
1634 /* create and simplify expr_i * expr_i */
1636 coefs[i*(i+1)/2 + i] = SCIPgetCoefsExprSum(base)[i] * SCIPgetCoefsExprSum(base)[i];
1637
1639 ownercreatedata) );
1641 ownercreatedata) ); /* this calls simplifyProduct */
1643 }
1644 /* create const * alpha_i expr_i */
1645 for( i = 0; i < nchildren; ++i )
1646 {
1649 }
1650
1651 constant = SCIPgetConstantExprSum(base);
1652 constant *= constant;
1653 /* create sum of all the above and simplify it with simplifySum since all of its children are simplified! */
1657 ownercreatedata) ); /* this calls simplifySum */
1658
1659 /* release everything */
1661 /* release the *created* expanded children */
1662 for( i = 0; i < nexpandedchildren - nchildren; ++i )
1663 {
1665 }
1667 SCIPfreeBufferArray(scip, &coefs);
1668
1669 return SCIP_OKAY;
1670 }
1671 }
1672 else
1673 {
1674 /* enforces POW9
1675 *
1676 * FIXME code of POW6 is very similar
1677 */
1678 if( SCIPexprGetNChildren(base) == 1
1680 && SCIPgetConstantExprSum(base) == 0.0
1681 && SCIPgetCoefsExprSum(base)[0] >= 0.0 )
1682 {
1684 SCIP_EXPR* aux;
1685 SCIP_Real newcoef;
1686
1687 SCIPdebugPrintf("[simplifyPow] seeing a sum with one term, exponent %g\n", exponent);
1688 /* assert SS7 holds */
1689 assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
1690
1691 /* create (pow n expr) and simplify it
1692 * note: we call simplifyPow directly, since we know that `expr` is simplified */
1694 ownercreatedata) );
1696 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1697
1698 /* create (sum (pow n expr)) and simplify it
1699 * this calls simplifySum directly, since we know its child is simplified! */
1700 newcoef = pow(SCIPgetCoefsExprSum(base)[0], exponent);
1702 ownercreatedata) );
1704 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1706
1707 return SCIP_OKAY;
1708 }
1709 }
1710
1711 /* enforces POW8
1712 * given (pow n (pow expo expr)) we distribute the exponent:
1713 * -> (pow n*expo expr)
1714 * notes: n is not 1 or 0, see POW1-2 above
1715 */
1716 if( SCIPisExprPower(scip, base) )
1717 {
1718 SCIP_Real newexponent;
1719 SCIP_Real baseexponent;
1720
1722 newexponent = baseexponent * exponent;
1723
1724 /* some checks (see POW8 definition in scip_expr.h) to make sure we don't loose an
1725 * implicit SCIPexprGetChildren(base)[0] >= 0 constraint
1726 *
1727 * if newexponent is fractional, then we will still need expr >= 0
1728 * if both exponents were integer, then we never required and will not require expr >= 0
1729 * if base exponent was an even integer, then we did not require expr >= 0
1730 * (but may need to use |expr|^newexponent)
1731 */
1732 if( !EPSISINT(newexponent, 0.0) ||
1733 (EPSISINT(baseexponent, 0.0) && EPSISINT(exponent, 0.0)) ||
1734 (EPSISINT(baseexponent, 0.0) && ((int)baseexponent) % 2 == 0) )
1735 {
1736 SCIP_EXPR* aux;
1737
1738 if( EPSISINT(baseexponent, 0.0) && ((int)baseexponent) % 2 == 0 &&
1739 (!EPSISINT(newexponent, 0.0) || ((int)newexponent) % 2 == 1) )
1740 {
1741 /* If base exponent was even integer and new exponent will be fractional,
1742 * then simplify to |expr|^newexponent to allow eval for expr < 0.
1743 * If base exponent was even integer and new exponent will be odd integer,
1744 * then simplify to |expr|^newexponent to preserve value for expr < 0.
1745 */
1747
1750 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1753 }
1754 else
1755 {
1757 ownercreatedata) );
1758 }
1759
1761 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
1762
1763 return SCIP_OKAY;
1764 }
1765 }
1766
1767 SCIPdebugPrintf("[simplifyPow] power is simplified\n");
1768 *simplifiedexpr = expr;
1769
1770 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
1772
1773 return SCIP_OKAY;
1774}
1775
1776/** expression handler copy callback */
1777static
1779{ /*lint --e{715}*/
1781
1782 return SCIP_OKAY;
1783}
1784
1785/** expression handler free callback */
1786static
1788{ /*lint --e{715}*/
1791
1793
1794 return SCIP_OKAY;
1795}
1796
1797/** expression data copy callback */
1798static
1815
1816/** expression data free callback */
1817static
1819{ /*lint --e{715}*/
1820 SCIP_EXPRDATA* exprdata;
1821
1822 assert(expr != NULL);
1823
1824 exprdata = SCIPexprGetData(expr);
1825 assert(exprdata != NULL);
1826
1827 SCIPfreeBlockMemory(scip, &exprdata);
1828 SCIPexprSetData(expr, NULL);
1829
1830 return SCIP_OKAY;
1831}
1832
1833/** expression print callback */
1834/** @todo: use precedence for better printing */
1835static
1837{ /*lint --e{715}*/
1838 assert(expr != NULL);
1839
1840 /**! [SnippetExprPrintPow] */
1841 switch( stage )
1842 {
1844 {
1845 /* print function with opening parenthesis */
1846 SCIPinfoMessage(scip, file, "(");
1847 break;
1848 }
1849
1851 {
1852 assert(currentchild == 0);
1853 break;
1854 }
1855
1857 {
1858 SCIP_Real exponent = SCIPgetExponentExprPow(expr);
1859
1860 /* print closing parenthesis */
1861 if( exponent >= 0.0 )
1862 SCIPinfoMessage(scip, file, ")^%g", exponent);
1863 else
1864 SCIPinfoMessage(scip, file, ")^(%g)", exponent);
1865
1866 break;
1867 }
1868
1870 default:
1871 break;
1872 }
1873 /**! [SnippetExprPrintPow] */
1874
1875 return SCIP_OKAY;
1876}
1877
1878/** expression point evaluation callback */
1879static
1881{ /*lint --e{715}*/
1882 SCIP_Real exponent;
1883 SCIP_Real base;
1884
1885 assert(expr != NULL);
1886 assert(SCIPexprGetNChildren(expr) == 1);
1888
1889 exponent = SCIPgetExponentExprPow(expr);
1891
1892 *val = pow(base, exponent);
1893
1894 /* if there is a domain, pole, or range error, pow() should return some kind of NaN, infinity, or HUGE_VAL
1895 * we could also work with floating point exceptions or errno, but I am not sure this would be thread-safe
1896 */
1897 if( !SCIPisFinite(*val) || *val == HUGE_VAL || *val == -HUGE_VAL )
1898 *val = SCIP_INVALID;
1899
1900 return SCIP_OKAY;
1901}
1902
1903/** derivative evaluation callback
1904 *
1905 * computes <gradient, children.dot>
1906 * if expr is child^p, then computes
1907 * p child^(p-1) dot(child)
1908 */
1909static
1911{ /*lint --e{715}*/
1912 SCIP_EXPR* child;
1913 SCIP_Real exponent;
1914
1915 assert(expr != NULL);
1916 assert(SCIPexprGetData(expr) != NULL);
1918
1919 child = SCIPexprGetChildren(expr)[0];
1920 assert(child != NULL);
1921 assert(!SCIPisExprValue(scip, child));
1923
1924 exponent = SCIPgetExponentExprPow(expr);
1925 assert(exponent != 0.0);
1926
1927 /* x^exponent is not differentiable for x = 0 and exponent in ]0,1[ */
1928 if( exponent > 0.0 && exponent < 1.0 && SCIPexprGetEvalValue(child) == 0.0 )
1929 *dot = SCIP_INVALID;
1930 else
1931 *dot = exponent * pow(SCIPexprGetEvalValue(child), exponent - 1.0) * SCIPexprGetDot(child);
1932
1933 return SCIP_OKAY;
1934}
1935
1936/** expression backward forward derivative evaluation callback
1937 *
1938 * computes partial/partial child ( <gradient, children.dot> )
1939 * if expr is child^n, then computes
1940 * n * (n - 1) child^(n-2) dot(child)
1941 */
1942static
1944{ /*lint --e{715}*/
1945 SCIP_EXPR* child;
1946 SCIP_Real exponent;
1947
1948 assert(expr != NULL);
1949 assert(SCIPexprGetData(expr) != NULL);
1951 assert(childidx == 0);
1952
1953 child = SCIPexprGetChildren(expr)[0];
1954 assert(child != NULL);
1955 assert(!SCIPisExprValue(scip, child));
1957
1958 exponent = SCIPgetExponentExprPow(expr);
1959 assert(exponent != 0.0);
1960
1961 /* x^exponent is not twice differentiable for x = 0 and exponent in ]0,1[ u ]1,2[ */
1962 if( exponent > 0.0 && exponent < 2.0 && SCIPexprGetEvalValue(child) == 0.0 && exponent != 1.0 )
1963 *bardot = SCIP_INVALID;
1964 else
1965 *bardot = exponent * (exponent - 1.0) * pow(SCIPexprGetEvalValue(child), exponent - 2.0) * SCIPexprGetDot(child);
1966
1967 return SCIP_OKAY;
1968}
1969
1970/** expression derivative evaluation callback */
1971static
1973{ /*lint --e{715}*/
1974 SCIP_EXPR* child;
1975 SCIP_Real childval;
1976 SCIP_Real exponent;
1977
1978 assert(expr != NULL);
1979 assert(SCIPexprGetData(expr) != NULL);
1980 assert(childidx == 0);
1981
1982 child = SCIPexprGetChildren(expr)[0];
1983 assert(child != NULL);
1984 assert(!SCIPisExprValue(scip, child));
1985
1988
1989 exponent = SCIPgetExponentExprPow(expr);
1990 assert(exponent != 0.0);
1991
1992 /* x^exponent is not differentiable for x = 0 and exponent in ]0,1[ */
1993 if( exponent > 0.0 && exponent < 1.0 && childval == 0.0 )
1994 *val = SCIP_INVALID;
1995 else
1996 *val = exponent * pow(childval, exponent - 1.0);
1997
1998 return SCIP_OKAY;
1999}
2000
2001/** expression interval evaluation callback */
2002static
2004{ /*lint --e{715}*/
2006 SCIP_Real exponent;
2007
2008 assert(expr != NULL);
2009 assert(SCIPexprGetNChildren(expr) == 1);
2010
2012
2013 exponent = SCIPgetExponentExprPow(expr);
2014
2015 if( exponent < 0.0 )
2016 {
2020
2021 if( exprhdlrdata->minzerodistance > 0.0 )
2022 {
2023 /* avoid small interval around 0 if possible, see also reversepropPow */
2024 if( childinterval.inf > -exprhdlrdata->minzerodistance && childinterval.inf < exprhdlrdata->minzerodistance )
2025 {
2026 if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2027 {
2028 SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2029 "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2030 exponent, childinterval.inf, exprhdlrdata->minzerodistance);
2031 SCIPinfoMessage(scip, NULL, "Expression: ");
2032 SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2033 SCIPinfoMessage(scip, NULL, "\n");
2034 exprhdlrdata->warnedonpole = TRUE;
2035 }
2036 childinterval.inf = exprhdlrdata->minzerodistance;
2037 }
2038 else if( childinterval.sup < exprhdlrdata->minzerodistance
2039 && childinterval.sup > -exprhdlrdata->minzerodistance )
2040 {
2041 if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2042 {
2043 SCIPinfoMessage(scip, NULL, "Changing upper bound for child of pow(.,%g) from %g to %g.\n"
2044 "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2045 exponent, childinterval.sup, -exprhdlrdata->minzerodistance);
2046 SCIPinfoMessage(scip, NULL, "Expression: ");
2047 SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2048 SCIPinfoMessage(scip, NULL, "\n");
2049 exprhdlrdata->warnedonpole = TRUE;
2050 }
2051 childinterval.sup = -exprhdlrdata->minzerodistance;
2052 }
2053 }
2054 }
2055
2057 {
2059 return SCIP_OKAY;
2060 }
2061
2063
2064 /* make sure 0^negative is an empty interval, as some other codes do not handle intervals like [inf,inf] well
2065 * TODO maybe change SCIPintervalPowerScalar?
2066 */
2067 if( exponent < 0.0 && childinterval.inf == 0.0 && childinterval.sup == 0.0 )
2069
2070 return SCIP_OKAY;
2071}
2072
2073/** expression estimator callback */
2074static
2076{ /*lint --e{715}*/
2077 SCIP_EXPRDATA* exprdata;
2078 SCIP_EXPR* child;
2079 SCIP_Real childlb;
2080 SCIP_Real childub;
2081 SCIP_Real exponent;
2082 SCIP_Bool isinteger;
2083
2084 assert(scip != NULL);
2085 assert(expr != NULL);
2086 assert(SCIPexprGetNChildren(expr) == 1);
2088 assert(refpoint != NULL);
2089 assert(coefs != NULL);
2090 assert(constant != NULL);
2091 assert(islocal != NULL);
2092 assert(branchcand != NULL);
2093 assert(*branchcand == TRUE); /* the default */
2094 assert(success != NULL);
2095
2096 *success = FALSE;
2097
2098 /* get aux variables: we over- or underestimate childvar^exponent */
2099 child = SCIPexprGetChildren(expr)[0];
2100 assert(child != NULL);
2101
2102 SCIPdebugMsg(scip, "%sestimation of x^%g at x=%.15g\n",
2103 overestimate ? "over" : "under", SCIPgetExponentExprPow(expr), *refpoint);
2104
2105 /* we can not generate a cut at +/- infinity */
2107 return SCIP_OKAY;
2108
2109 childlb = localbounds[0].inf;
2110 childub = localbounds[0].sup;
2111
2112 exprdata = SCIPexprGetData(expr);
2113 exponent = exprdata->exponent;
2114 assert(exponent != 1.0 && exponent != 0.0); /* this should have been simplified */
2115
2116 /* if child is constant, then return a constant estimator
2117 * this can help with small infeasibilities if boundtightening is relaxing bounds too much
2118 */
2119 if( childlb == childub )
2120 {
2121 *coefs = 0.0;
2122 *constant = pow(childlb, exponent);
2123 *success = TRUE;
2124 *islocal = globalbounds[0].inf != globalbounds[0].sup;
2125 *branchcand = FALSE;
2126 return SCIP_OKAY;
2127 }
2128
2129 isinteger = EPSISINT(exponent, 0.0);
2130
2131 /* if exponent is not integral, then child must be non-negative */
2132 if( !isinteger && childlb < 0.0 )
2133 {
2134 /* somewhere we should have tightened the bound on x, but small tightening are not always applied by SCIP
2135 * it is ok to do this tightening here, but let's assert that we were close to 0.0 already
2136 */
2138 childlb = 0.0;
2139 }
2140 assert(isinteger || childlb >= 0.0);
2141
2142 SCIP_CALL( buildPowEstimator(scip, exprdata, overestimate, childlb, childub, globalbounds[0].inf,
2143 globalbounds[0].sup, SCIPexprIsIntegral(child), MAX(childlb, *refpoint), exponent, coefs,
2144 constant, success, islocal, branchcand) );
2145
2146 return SCIP_OKAY;
2147}
2148
2149/** expression reverse propagaton callback */
2150static
2152{ /*lint --e{715}*/
2153 SCIP_INTERVAL child;
2155 SCIP_Real exponent;
2156
2157 assert(scip != NULL);
2158 assert(expr != NULL);
2159 assert(SCIPexprGetNChildren(expr) == 1);
2160
2161 exponent = SCIPgetExponentExprPow(expr);
2162 child = childrenbounds[0];
2163
2164 SCIPdebugMsg(scip, "reverseprop x^%g in [%.15g,%.15g], x = [%.15g,%.15g]", exponent, bounds.inf, bounds.sup,
2165 child.inf, child.sup);
2166
2168 {
2169 *infeasible = TRUE;
2170 return SCIP_OKAY;
2171 }
2172
2174 {
2175 /* if exponent is not integral, then make sure that child is non-negative */
2176 if( !EPSISINT(exponent, 0.0) && child.inf < 0.0 )
2177 {
2178 SCIPintervalSetBounds(&interval, 0.0, child.sup);
2179 }
2180 else
2181 {
2182 SCIPdebugMsgPrint(scip, "-> no improvement\n");
2183 return SCIP_OKAY;
2184 }
2185 }
2186 else
2187 {
2188 /* f = pow(c0, alpha) -> c0 = pow(f, 1/alpha) */
2190 }
2191
2192 if( exponent < 0.0 )
2193 {
2195
2198
2199 if( exprhdlrdata->minzerodistance > 0.0 )
2200 {
2201 /* push lower bound from >= -epsilon to >= epsilon to avoid pole at 0 (domain error)
2202 * push upper bound from <= epsilon to <= -epsilon to avoid pole at 0 (domain error)
2203 * this can lead to a cutoff if domain would otherwise be very close around 0
2204 */
2205 if( interval.inf > -exprhdlrdata->minzerodistance && interval.inf < exprhdlrdata->minzerodistance )
2206 {
2207 if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2208 {
2209 SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2210 "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2211 exponent, interval.inf, exprhdlrdata->minzerodistance);
2212 SCIPinfoMessage(scip, NULL, "Expression: ");
2213 SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2214 SCIPinfoMessage(scip, NULL, "\n");
2215 exprhdlrdata->warnedonpole = TRUE;
2216 }
2217 interval.inf = exprhdlrdata->minzerodistance;
2218 }
2219 else if( interval.sup < exprhdlrdata->minzerodistance && interval.sup > -exprhdlrdata->minzerodistance )
2220 {
2221 if( !exprhdlrdata->warnedonpole && SCIPgetVerbLevel(scip) > SCIP_VERBLEVEL_NONE )
2222 {
2223 SCIPinfoMessage(scip, NULL, "Changing lower bound for child of pow(.,%g) from %g to %g.\n"
2224 "Check your model formulation or use option expr/" POWEXPRHDLR_NAME "/minzerodistance to avoid this warning.\n",
2225 exponent, interval.sup, -exprhdlrdata->minzerodistance);
2226 SCIPinfoMessage(scip, NULL, "Expression: ");
2227 SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2228 SCIPinfoMessage(scip, NULL, "\n");
2229 exprhdlrdata->warnedonpole = TRUE;
2230 }
2231 interval.sup = -exprhdlrdata->minzerodistance;
2232 }
2233 }
2234 }
2235
2236 SCIPdebugMsgPrint(scip, " -> [%.15g,%.15g]\n", interval.inf, interval.sup);
2237
2239
2240 return SCIP_OKAY;
2241}
2242
2243/** initial estimates callback for a power expression */
2244static
2246{
2247 SCIP_EXPRDATA* exprdata;
2248 SCIP_EXPR* child;
2249 SCIP_Real childlb;
2250 SCIP_Real childub;
2251 SCIP_Real exponent;
2252 SCIP_Bool isinteger;
2253 SCIP_Bool branchcand;
2254 SCIP_Bool success;
2255 SCIP_Bool islocal;
2258 SCIP_Bool overest[6] = {FALSE, FALSE, FALSE, TRUE, TRUE, TRUE};
2259 int i;
2260
2261 assert(scip != NULL);
2262 assert(expr != NULL);
2263
2264 child = SCIPexprGetChildren(expr)[0];
2265 assert(child != NULL);
2266
2267 childlb = bounds[0].inf;
2268 childub = bounds[0].sup;
2269
2270 /* if child is essentially constant, then there should be no point in separation */
2271 if( SCIPisEQ(scip, childlb, childub) )
2272 {
2273 SCIPdebugMsg(scip, "skip initestimates as child seems essentially fixed [%.15g,%.15g]\n", childlb, childub);
2274 return SCIP_OKAY;
2275 }
2276
2277 exprdata = SCIPexprGetData(expr);
2278 exponent = exprdata->exponent;
2279 assert(exponent != 1.0 && exponent != 0.0); /* this should have been simplified */
2280
2281 isinteger = EPSISINT(exponent, 0.0);
2282
2283 /* if exponent is not integral, then child must be non-negative */
2284 if( !isinteger && childlb < 0.0 )
2285 {
2286 /* somewhere we should have tightened the bound on x, but small tightening are not always applied by SCIP
2287 * it is ok to do this tightening here, but let's assert that we were close to 0.0 already
2288 */
2290 childlb = 0.0;
2291 }
2292 assert(isinteger || childlb >= 0.0);
2293
2294 /* TODO simplify to get 3 refpoints for either under- or overestimate */
2296 overestimate) );
2297
2298 for( i = 0; i < 6 && *nreturned < SCIP_EXPR_MAXINITESTIMATES; ++i )
2299 {
2300 SCIP_Real refpoint;
2301
2302 if( (overest[i] && !overestimate) || (!overest[i] && overestimate) )
2303 continue;
2304
2305 assert(overest[i] || i < 3); /* make sure that no out-of-bounds array access will be attempted */
2306 refpoint = overest[i] ? refpointsover[i % 3] : refpointsunder[i % 3];
2307
2308 if( refpoint == SCIP_INVALID )
2309 continue;
2310
2312
2313 branchcand = TRUE;
2315 SCIPexprIsIntegral(child), refpoint, exponent, coefs[*nreturned], &constant[*nreturned],
2316 &success, &islocal, &branchcand) );
2317
2318 if( success )
2319 {
2320 SCIPdebugMsg(scip, "initestimate x^%g for base in [%g,%g] at ref=%g, over:%u -> %g*x+%g\n", exponent,
2321 childlb, childub, refpoint, overest[i], coefs[*nreturned][0], constant[*nreturned]);
2322 ++*nreturned;
2323 }
2324 }
2325
2326 return SCIP_OKAY;
2327}
2328
2329/** expression hash callback */
2330static
2332{ /*lint --e{715}*/
2333 assert(scip != NULL);
2334 assert(expr != NULL);
2335 assert(SCIPexprGetNChildren(expr) == 1);
2336 assert(hashkey != NULL);
2338
2339 /* TODO include exponent into hashkey */
2341 *hashkey ^= childrenhashes[0];
2342
2343 return SCIP_OKAY;
2344}
2345
2346/** expression curvature detection callback */
2347static
2349{ /*lint --e{715}*/
2350 SCIP_EXPR* child;
2352 SCIP_Real exponent;
2353
2354 assert(scip != NULL);
2355 assert(expr != NULL);
2357 assert(childcurv != NULL);
2358 assert(success != NULL);
2359 assert(SCIPexprGetNChildren(expr) == 1);
2360
2361 exponent = SCIPgetExponentExprPow(expr);
2362 child = SCIPexprGetChildren(expr)[0];
2363 assert(child != NULL);
2364
2367
2369 /* SCIPexprcurvPowerInv return unknown actually means that curv cannot be obtained */
2371
2372 return SCIP_OKAY;
2373}
2374
2375/** expression monotonicity detection callback */
2376static
2378{ /*lint --e{715}*/
2380 SCIP_Real exponent;
2381 SCIP_Real inf;
2382 SCIP_Real sup;
2383 SCIP_Bool expisint;
2384
2385 assert(scip != NULL);
2386 assert(expr != NULL);
2387 assert(result != NULL);
2388 assert(SCIPexprGetNChildren(expr) == 1);
2389 assert(childidx == 0);
2390
2391 assert(SCIPexprGetChildren(expr)[0] != NULL);
2394
2398 exponent = SCIPgetExponentExprPow(expr);
2399 expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
2400
2401 if( expisint )
2402 {
2403 SCIP_Bool expisodd = ceil(exponent/2) != exponent/2;
2404
2405 if( expisodd )
2406 {
2407 /* x^1, x^3, ... */
2408 if( exponent >= 0.0 )
2410
2411 /* ..., x^-3, x^-1 are decreasing if 0 is not in ]inf,sup[ */
2412 else if( inf >= 0.0 || sup <= 0.0 )
2414 }
2415 /* ..., x^-4, x^-2, x^2, x^4, ... */
2416 else
2417 {
2418 /* function is not monotone if 0 is in ]inf,sup[ */
2419 if( inf >= 0.0 )
2420 *result = exponent >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
2421 else if( sup <= 0.0 )
2422 *result = exponent >= 0.0 ? SCIP_MONOTONE_DEC : SCIP_MONOTONE_INC;
2423 }
2424 }
2425 else
2426 {
2427 /* note that the expression is not defined for negative input values
2428 *
2429 * - increasing iff exponent >= 0
2430 * - decreasing iff exponent <= 0
2431 */
2432 *result = exponent >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
2433 }
2434
2435 return SCIP_OKAY;
2436}
2437
2438/** expression integrality detection callback */
2439static
2441{ /*lint --e{715}*/
2442 SCIP_EXPR* child;
2443 SCIP_Real exponent;
2444 SCIP_Bool expisint;
2445
2446 assert(scip != NULL);
2447 assert(expr != NULL);
2448 assert(isintegral != NULL);
2449 assert(SCIPexprGetNChildren(expr) == 1);
2450
2451 *isintegral = FALSE;
2452
2453 child = SCIPexprGetChildren(expr)[0];
2454 assert(child != NULL);
2455
2456 /* expression can not be integral if child is not */
2457 if( !SCIPexprIsIntegral(child) )
2458 return SCIP_OKAY;
2459
2460 exponent = SCIPgetExponentExprPow(expr);
2461 assert(exponent != 0.0);
2462 expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
2463
2464 /* expression is integral if and only if exponent non-negative and integral */
2465 *isintegral = expisint && exponent >= 0.0;
2466
2467 return SCIP_OKAY;
2468}
2469
2470/** simplifies a signpower expression
2471 */
2472static
2474{ /*lint --e{715}*/
2475 SCIP_EXPR* base;
2476 SCIP_Real exponent;
2477
2478 assert(scip != NULL);
2479 assert(expr != NULL);
2481 assert(SCIPexprGetNChildren(expr) == 1);
2482
2483 base = SCIPexprGetChildren(expr)[0];
2484 assert(base != NULL);
2485
2486 exponent = SCIPgetExponentExprPow(expr);
2487 SCIPdebugPrintf("[simplifySignpower] simplifying power with expo %g\n", exponent);
2488 assert(exponent >= 1.0);
2489
2490 /* enforces SPOW2 */
2491 if( exponent == 1.0 )
2492 {
2493 SCIPdebugPrintf("[simplifySignpower] POW2\n");
2496 return SCIP_OKAY;
2497 }
2498
2499 /* enforces SPOW3 */
2500 if( SCIPisExprValue(scip, base) )
2501 {
2502 SCIP_Real baseval;
2503
2504 SCIPdebugPrintf("[simplifySignpower] POW3\n");
2506
2509
2510 return SCIP_OKAY;
2511 }
2512
2513 /* enforces SPOW11 (exp(x)^n = exp(n*x))
2514 * since exp() is always nonnegative, we can treat signpower as normal power here
2515 */
2516 if( SCIPisExprExp(scip, base) )
2517 {
2518 SCIP_EXPR* child;
2519 SCIP_EXPR* prod;
2522
2523 SCIPdebugPrintf("[simplifySignpower] POW11\n");
2524 child = SCIPexprGetChildren(base)[0];
2525
2526 /* multiply child of exponential with exponent */
2528
2529 /* simplify product */
2532
2533 /* create exponential with new child */
2536
2537 /* the final simplified expression is the simplification of the just created exponential */
2540
2541 return SCIP_OKAY;
2542 }
2543
2544 /* enforces SPOW6 */
2545 if( EPSISINT(exponent, 0.0) && ((int)exponent) % 2 == 1 )
2546 {
2547 SCIP_EXPR* aux;
2548
2549 /* we do not just change the expression data of expression to say it is a normal power, since, at the moment,
2550 * simplify identifies that expressions changed by checking that the pointer of the input expression is
2551 * different from the returned (simplified) expression
2552 */
2554
2556 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2557
2558 return SCIP_OKAY;
2559 }
2560
2561 /* enforces SPOW10 */
2562 if( SCIPisExprVar(scip, base) )
2563 {
2565
2566 SCIPdebugPrintf("[simplifySignpower] POW10\n");
2568
2569 assert(basevar != NULL);
2570
2572 {
2575 return SCIP_OKAY;
2576 }
2577 }
2578
2579 /* TODO if( SCIPisExprSignpower(scip, base) ... */
2580
2581 /* enforces SPOW8
2582 * given (signpow n (pow expo expr)) we distribute the exponent:
2583 * -> (signpow n*expo expr) for even n (i.e., sign(x^n) * |x|^n = 1 * x^n)
2584 * notes: n is an even integer (see SPOW6 above)
2585 * FIXME: doesn't this extend to any exponent?
2586 * If (pow expo expr) can be negative, it should mean that (-1)^expo = -1
2587 * then (signpow n (pow expo expr)) = sign(expr^expo) * |expr^expo|^n
2588 * then sign(expr^expo) = sign(expr) and |expr^expo| = |expr|^expo and so
2589 * (signpow n (pow expo expr)) = sign(expr^expo) * |expr^expo|^n = sign(expr) * |expr|^(expo*n) = signpow n*expo expr
2590 */
2591 if( EPSISINT(exponent, 0.0) && SCIPisExprPower(scip, base) )
2592 {
2593 SCIP_EXPR* aux;
2594 SCIP_Real newexponent;
2595
2596 assert(((int)exponent) % 2 == 0 ); /* odd case should have been handled by SPOW6 */
2597
2602
2603 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2604
2605 return SCIP_OKAY;
2606 }
2607
2608 /* enforces SPOW9 */
2609 if( SCIPisExprSum(scip, base)
2610 && SCIPexprGetNChildren(base) == 1
2611 && SCIPgetConstantExprSum(base) == 0.0 )
2612 {
2614 SCIP_EXPR* aux;
2615 SCIP_Real newcoef;
2616
2617 SCIPdebugPrintf("[simplifySignpower] seeing a sum with one term, exponent %g\n", exponent);
2618 /* assert SS7 holds */
2619 assert(SCIPgetCoefsExprSum(base)[0] != 1.0);
2620
2621 /* create (signpow n expr) and simplify it
2622 * note: we call simplifySignpower directly, since we know that `expr` is simplified */
2627 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2628
2629 /* create (sum (signpow n expr)) and simplify it
2630 * this calls simplifySum directly, since we know its child is simplified */
2633 SCIP_CALL( SCIPreleaseExpr(scip, &aux) );
2635 return SCIP_OKAY;
2636 }
2637
2638 SCIPdebugPrintf("[simplifySignpower] signpower is simplified\n");
2639 *simplifiedexpr = expr;
2640
2641 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
2643
2644 return SCIP_OKAY;
2645}
2646
2647/** expression handler copy callback */
2648static
2650{ /*lint --e{715}*/
2652
2653 return SCIP_OKAY;
2654}
2655
2656/** expression print callback */
2657static
2659{ /*lint --e{715}*/
2660 assert(expr != NULL);
2661
2662 switch( stage )
2663 {
2665 {
2666 SCIPinfoMessage(scip, file, "signpower(");
2667 break;
2668 }
2669
2671 {
2672 assert(currentchild == 0);
2673 break;
2674 }
2675
2677 {
2678 SCIPinfoMessage(scip, file, ",%g)", SCIPgetExponentExprPow(expr));
2679 break;
2680 }
2681
2683 default:
2684 break;
2685 }
2686
2687 return SCIP_OKAY;
2688}
2689
2690/** expression parse callback */
2691static
2693{ /*lint --e{715}*/
2695 SCIP_Real exponent;
2696
2697 assert(expr != NULL);
2698
2699 /**! [SnippetExprParseSignpower] */
2700 /* parse child expression string */
2702 assert(childexpr != NULL);
2703
2704 string = *endstring;
2705 while( *string == ' ' )
2706 ++string;
2707
2708 if( *string != ',' )
2709 {
2710 SCIPerrorMessage("Expected comma after first argument of signpower().\n");
2711 return SCIP_READERROR;
2712 }
2713 ++string;
2714
2715 if( !SCIPparseReal(scip, string, &exponent, (char**)endstring) )
2716 {
2717 SCIPerrorMessage("Expected numeric exponent for second argument of signpower().\n");
2718 return SCIP_READERROR;
2719 }
2720
2721 if( exponent <= 1.0 || SCIPisInfinity(scip, exponent) )
2722 {
2723 SCIPerrorMessage("Expected finite exponent >= 1.0 for signpower().\n");
2724 return SCIP_READERROR;
2725 }
2726
2727 /* create signpower expression */
2729 assert(*expr != NULL);
2730
2731 /* release child expression since it has been captured by the signpower expression */
2733
2734 *success = TRUE;
2735 /**! [SnippetExprParseSignpower] */
2736
2737 return SCIP_OKAY;
2738}
2739
2740/** expression point evaluation callback */
2741static
2743{ /*lint --e{715}*/
2744 SCIP_Real exponent;
2745 SCIP_Real base;
2746
2747 assert(expr != NULL);
2748 assert(SCIPexprGetNChildren(expr) == 1);
2750
2751 exponent = SCIPgetExponentExprPow(expr);
2753
2754 *val = SIGN(base) * pow(REALABS(base), exponent);
2755
2756 /* if there is a range error, pow() should return some kind of infinity, or HUGE_VAL
2757 * we could also work with floating point exceptions or errno, but I am not sure this would be thread-safe
2758 */
2759 if( !SCIPisFinite(*val) || *val == HUGE_VAL || *val == -HUGE_VAL )
2760 *val = SCIP_INVALID;
2761
2762 return SCIP_OKAY;
2763}
2764
2765/** expression derivative evaluation callback */
2766static
2768{ /*lint --e{715}*/
2769 SCIP_EXPR* child;
2770 SCIP_Real childval;
2771 SCIP_Real exponent;
2772
2773 assert(expr != NULL);
2774 assert(SCIPexprGetData(expr) != NULL);
2775 assert(childidx == 0);
2776
2777 child = SCIPexprGetChildren(expr)[0];
2778 assert(child != NULL);
2779 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
2780
2783
2784 exponent = SCIPgetExponentExprPow(expr);
2785 assert(exponent >= 1.0);
2786
2787 *val = exponent * pow(REALABS(childval), exponent - 1.0);
2788
2789 return SCIP_OKAY;
2790}
2791
2792/** expression interval evaluation callback */
2793static
2812
2813/** expression estimator callback */
2814static
2816{ /*lint --e{715}*/
2817 SCIP_EXPRDATA* exprdata;
2818 SCIP_Real childlb;
2819 SCIP_Real childub;
2820 SCIP_Real childglb;
2821 SCIP_Real childgub;
2822 SCIP_Real exponent;
2823
2824 assert(scip != NULL);
2825 assert(expr != NULL);
2826 assert(SCIPexprGetNChildren(expr) == 1);
2828 assert(refpoint != NULL);
2829 assert(coefs != NULL);
2830 assert(constant != NULL);
2831 assert(islocal != NULL);
2832 assert(branchcand != NULL);
2833 assert(*branchcand == TRUE); /* the default */
2834 assert(success != NULL);
2835
2836 *success = FALSE;
2837
2838 SCIPdebugMsg(scip, "%sestimation of signed x^%g at x=%g\n", overestimate ? "over" : "under",
2840
2841 /* we can not generate a cut at +/- infinity */
2843 return SCIP_OKAY;
2844
2845 childlb = localbounds[0].inf;
2846 childub = localbounds[0].sup;
2847
2848 childglb = globalbounds[0].inf;
2849 childgub = globalbounds[0].sup;
2850
2851 exprdata = SCIPexprGetData(expr);
2852 exponent = exprdata->exponent;
2853 assert(exponent > 1.0); /* exponent == 1 should have been simplified */
2854
2855 /* if child is constant, then return a constant estimator
2856 * this can help with small infeasibilities if boundtightening is relaxing bounds too much
2857 */
2858 if( childlb == childub )
2859 {
2860 *coefs = 0.0;
2861 *constant = SIGN(childlb)*pow(REALABS(childlb), exponent);
2862 *success = TRUE;
2864 *branchcand = FALSE;
2865 return SCIP_OKAY;
2866 }
2867
2868 if( childlb >= 0.0 )
2869 {
2870 estimateParabola(scip, exponent, overestimate, childlb, childub, MAX(0.0, *refpoint), constant, coefs,
2871 islocal, success);
2872
2873 *branchcand = *islocal;
2874
2875 /* if odd or signed power, then check whether tangent on parabola is also globally valid, that is
2876 * reference point is right of -root*global-lower-bound
2877 */
2878 if( !*islocal && childglb < 0.0 )
2879 {
2881 *islocal = TRUE;
2882 else
2883 {
2884 if( exprdata->root == SCIP_INVALID )
2885 {
2886 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2887 }
2888 *islocal = *refpoint < exprdata->root * (-childglb);
2889 }
2890 }
2891 }
2892 else /* and childlb < 0.0 due to previous if */
2893 {
2894 /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
2895 if( exprdata->root == SCIP_INVALID && childgub > 0.0 )
2896 {
2897 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2898 }
2899 estimateSignedpower(scip, exponent, exprdata->root, overestimate, childlb, childub, *refpoint,
2900 childglb, childgub, constant, coefs, islocal, branchcand, success);
2901 }
2902
2903 return SCIP_OKAY;
2904}
2905
2906/** initial estimates callback for a signpower expression */
2907static
2909{
2910 SCIP_EXPRDATA* exprdata;
2911 SCIP_Real childlb;
2912 SCIP_Real childub;
2913 SCIP_Real exponent;
2914 SCIP_Bool branchcand;
2915 SCIP_Bool success;
2916 SCIP_Bool islocal;
2919 SCIP_Bool overest[6] = {FALSE, FALSE, FALSE, TRUE, TRUE, TRUE};
2920 SCIP_Real refpoint;
2921 int i;
2922
2923 assert(scip != NULL);
2924 assert(expr != NULL);
2925 assert(SCIPexprGetNChildren(expr) == 1);
2927
2928 childlb = bounds[0].inf;
2929 childub = bounds[0].sup;
2930
2931 /* if child is essentially constant, then there should be no point in separation */
2932 if( SCIPisEQ(scip, childlb, childub) )
2933 {
2934 SCIPdebugMsg(scip, "skip initestimates as child seems essentially fixed [%.15g,%.15g]\n", childlb, childub);
2935 return SCIP_OKAY;
2936 }
2937
2938 exprdata = SCIPexprGetData(expr);
2939 exponent = exprdata->exponent;
2940 assert(exponent > 1.0); /* this should have been simplified */
2941
2942 if( childlb >= 0.0 )
2943 {
2944 if( !overestimate )
2946 if( overestimate && !SCIPisInfinity(scip, childub) )
2947 refpointsover[0] = (childlb + childub) / 2.0;
2948 }
2949 else if( childub <= 0.0 )
2950 {
2951 if( !overestimate && !SCIPisInfinity(scip, -childlb) )
2952 refpointsunder[0] = (childlb + childub) / 2.0;
2953 if( overestimate )
2955 }
2956 else
2957 {
2958 SCIP_CALL( addSignpowerRefpoints(scip, exprdata, childlb, childub, exponent, !overestimate, refpointsunder) );
2959 }
2960
2961 /* add cuts for all refpoints */
2962 for( i = 0; i < 6 && *nreturned < SCIP_EXPR_MAXINITESTIMATES; ++i )
2963 {
2964 if( (overest[i] && !overestimate) || (!overest[i] && overestimate) )
2965 continue;
2966
2967 assert(overest[i] || i < 3); /* make sure that no out-of-bounds array access will be attempted */
2968 refpoint = overest[i] ? refpointsover[i % 3] : refpointsunder[i % 3];
2969 if( refpoint == SCIP_INVALID )
2970 continue;
2972
2973 if( childlb >= 0 )
2974 {
2975 estimateParabola(scip, exponent, overest[i], childlb, childub, refpoint, &constant[*nreturned], coefs[*nreturned],
2976 &islocal, &success);
2977 }
2978 else
2979 {
2980 /* compute root if not known yet; only needed if mixed sign (global child ub > 0) */
2981 if( exprdata->root == SCIP_INVALID && childub > 0.0 )
2982 {
2983 SCIP_CALL( computeSignpowerRoot(scip, &exprdata->root, exponent) );
2984 }
2985 branchcand = TRUE;
2986 estimateSignedpower(scip, exponent, exprdata->root, overest[i], childlb, childub, refpoint,
2987 childlb, childub, &constant[*nreturned], coefs[*nreturned], &islocal,
2988 &branchcand, &success);
2989 }
2990
2991 if( success )
2992 ++*nreturned;
2993 }
2994
2995 return SCIP_OKAY;
2996}
2997
2998/** expression reverse propagaton callback */
2999static
3001{ /*lint --e{715}*/
3004 SCIP_Real exponent;
3005
3006 assert(scip != NULL);
3007 assert(expr != NULL);
3008 assert(SCIPexprGetNChildren(expr) == 1);
3009
3010 exponent = SCIPgetExponentExprPow(expr);
3011
3012 SCIPdebugMsg(scip, "reverseprop signpow(x,%g) in [%.15g,%.15g]", exponent, bounds.inf, bounds.sup);
3013
3015 {
3016 SCIPdebugMsgPrint(scip, "-> no improvement\n");
3017 return SCIP_OKAY;
3018 }
3019
3020 /* f = pow(c0, alpha) -> c0 = pow(f, 1/alpha) */
3021 SCIPintervalSet(&exprecip, exponent);
3023 if( exprecip.inf == exprecip.sup )
3024 {
3026 }
3027 else
3028 {
3033 }
3034
3035 SCIPdebugMsgPrint(scip, " -> [%.15g,%.15g]\n", interval.inf, interval.sup);
3036
3038
3039 return SCIP_OKAY;
3040}
3041
3042/** expression hash callback */
3043static
3045{ /*lint --e{715}*/
3046 assert(scip != NULL);
3047 assert(expr != NULL);
3048 assert(SCIPexprGetNChildren(expr) == 1);
3049 assert(hashkey != NULL);
3051
3052 /* TODO include exponent into hashkey */
3054 *hashkey ^= childrenhashes[0];
3055
3056 return SCIP_OKAY;
3057}
3058
3059/** expression curvature detection callback */
3060static
3062{ /*lint --e{715}*/
3063 SCIP_EXPR* child;
3065
3066 assert(scip != NULL);
3067 assert(expr != NULL);
3069 assert(childcurv != NULL);
3070 assert(success != NULL);
3071 assert(SCIPexprGetNChildren(expr) == 1);
3072
3073 child = SCIPexprGetChildren(expr)[0];
3074 assert(child != NULL);
3075
3078
3080 {
3081 /* signpower is only convex if argument is convex and non-negative */
3083 *success = childinterval.inf >= 0.0;
3084 }
3086 {
3087 /* signpower is only concave if argument is concave and non-positive */
3089 *success = childinterval.sup <= 0.0;
3090 }
3091 else
3092 *success = FALSE;
3093
3094 return SCIP_OKAY;
3095}
3096
3097/** expression monotonicity detection callback */
3098static
3100{ /*lint --e{715}*/
3101 assert(scip != NULL);
3102 assert(expr != NULL);
3103 assert(result != NULL);
3104
3106 return SCIP_OKAY;
3107}
3108
3109/** creates the handler for power expression and includes it into SCIP */
3111 SCIP* scip /**< SCIP data structure */
3112 )
3113{
3114 SCIP_EXPRHDLR* exprhdlr;
3116
3118
3121 assert(exprhdlr != NULL);
3122
3126 SCIPexprhdlrSetPrint(exprhdlr, printPow);
3130 SCIPexprhdlrSetHash(exprhdlr, hashPow);
3136
3137 SCIP_CALL( SCIPaddRealParam(scip, "expr/" POWEXPRHDLR_NAME "/minzerodistance",
3138 "minimal distance from zero to enforce for child in bound tightening",
3139 &exprhdlrdata->minzerodistance, FALSE, SCIPepsilon(scip), 0.0, 1.0, NULL, NULL) );
3140
3141 return SCIP_OKAY;
3142}
3143
3144/** creates the handler for signed power expression and includes it into SCIP */
3172
3173/** creates a power expression */
3175 SCIP* scip, /**< SCIP data structure */
3176 SCIP_EXPR** expr, /**< pointer where to store expression */
3177 SCIP_EXPR* child, /**< single child */
3178 SCIP_Real exponent, /**< exponent of the power expression */
3179 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
3180 void* ownercreatedata /**< data to pass to ownercreate */
3181 )
3182{
3183 SCIP_EXPRDATA* exprdata;
3184
3185 assert(expr != NULL);
3186 assert(child != NULL);
3187
3188 SCIP_CALL( createData(scip, &exprdata, exponent) );
3189 assert(exprdata != NULL);
3190
3191 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrPower(scip), exprdata, 1, &child, ownercreate,
3192 ownercreatedata) );
3193
3194 return SCIP_OKAY;
3195}
3196
3197/** creates a signpower expression */
3199 SCIP* scip, /**< SCIP data structure */
3200 SCIP_EXPR** expr, /**< pointer where to store expression */
3201 SCIP_EXPR* child, /**< single child */
3202 SCIP_Real exponent, /**< exponent of the power expression */
3203 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
3204 void* ownercreatedata /**< data to pass to ownercreate */
3205 )
3206{
3207 SCIP_EXPRDATA* exprdata;
3208
3209 assert(expr != NULL);
3210 assert(child != NULL);
3212
3213 SCIP_CALL( createData(scip, &exprdata, exponent) );
3214 assert(exprdata != NULL);
3215
3218
3219 return SCIP_OKAY;
3220}
3221
3222/** indicates whether expression is of signpower-type */ /*lint -e{715}*/
3224 SCIP* scip, /**< SCIP data structure */
3225 SCIP_EXPR* expr /**< expression */
3226 )
3227{ /*lint --e{715}*/
3228 assert(expr != NULL);
3229
3231}
3232
3233/** computes coefficients of linearization of a square term in a reference point */
3235 SCIP* scip, /**< SCIP data structure */
3236 SCIP_Real sqrcoef, /**< coefficient of square term */
3237 SCIP_Real refpoint, /**< point where to linearize */
3238 SCIP_Bool isint, /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */
3239 SCIP_Real* lincoef, /**< buffer to add coefficient of linearization */
3240 SCIP_Real* linconstant, /**< buffer to add constant of linearization */
3241 SCIP_Bool* success /**< buffer to set to FALSE if linearization has failed due to large numbers */
3242 )
3243{
3244 assert(scip != NULL);
3245 assert(lincoef != NULL);
3247 assert(success != NULL);
3248
3249 if( sqrcoef == 0.0 )
3250 return;
3251
3253 {
3254 *success = FALSE;
3255 return;
3256 }
3257
3258 if( !isint || SCIPisIntegral(scip, refpoint) )
3259 {
3260 SCIP_Real tmp;
3261
3262 /* sqrcoef * x^2 -> tangent in refpoint = sqrcoef * 2 * refpoint * (x - refpoint) */
3263
3264 tmp = sqrcoef * refpoint;
3265
3266 if( SCIPisInfinity(scip, 2.0 * REALABS(tmp)) )
3267 {
3268 *success = FALSE;
3269 return;
3270 }
3271
3272 *lincoef += 2.0 * tmp;
3273 tmp *= refpoint;
3274 *linconstant -= tmp;
3275 }
3276 else
3277 {
3278 /* sqrcoef * x^2 -> secant between f=floor(refpoint) and f+1 = sqrcoef * (f^2 + ((f+1)^2 - f^2) * (x-f))
3279 * = sqrcoef * (-f*(f+1) + (2*f+1)*x)
3280 */
3281 SCIP_Real f;
3282 SCIP_Real coef;
3283 SCIP_Real constant;
3284
3285 f = SCIPfloor(scip, refpoint);
3286
3287 coef = sqrcoef * (2.0 * f + 1.0);
3288 constant = -sqrcoef * f * (f + 1.0);
3289
3290 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
3291 {
3292 *success = FALSE;
3293 return;
3294 }
3295
3296 *lincoef += coef;
3297 *linconstant += constant;
3298 }
3299}
3300
3301/** computes coefficients of secant of a square term */
3303 SCIP* scip, /**< SCIP data structure */
3304 SCIP_Real sqrcoef, /**< coefficient of square term */
3305 SCIP_Real lb, /**< lower bound on variable */
3306 SCIP_Real ub, /**< upper bound on variable */
3307 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */
3308 SCIP_Real* linconstant, /**< buffer to add constant of secant */
3309 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
3310 )
3311{
3312 SCIP_Real coef;
3313 SCIP_Real constant;
3314
3315 assert(scip != NULL);
3316 assert(!SCIPisInfinity(scip, lb));
3317 assert(!SCIPisInfinity(scip, -ub));
3318 assert(SCIPisLE(scip, lb, ub));
3319 assert(lincoef != NULL);
3321 assert(success != NULL);
3322
3323 if( sqrcoef == 0.0 )
3324 return;
3325
3326 if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
3327 {
3328 /* unboundedness */
3329 *success = FALSE;
3330 return;
3331 }
3332
3333 /* sqrcoef * x^2 -> sqrcoef * (lb * lb + (ub*ub - lb*lb)/(ub-lb) * (x-lb)) = sqrcoef * (lb*lb + (ub+lb)*(x-lb))
3334 * = sqrcoef * ((lb+ub)*x - lb*ub)
3335 */
3336 coef = sqrcoef * (lb + ub);
3337 constant = -sqrcoef * lb * ub;
3338 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
3339 {
3340 *success = FALSE;
3341 return;
3342 }
3343
3344 *lincoef += coef;
3345 *linconstant += constant;
3346}
3347
3348/* from pub_expr.h */
3349
3350/** gets the exponent of a power or signed power expression */ /*lint -e{715}*/
3352 SCIP_EXPR* expr /**< expression */
3353 )
3354{
3355 SCIP_EXPRDATA* exprdata;
3356
3357 assert(expr != NULL);
3358
3359 exprdata = SCIPexprGetData(expr);
3360 assert(exprdata != NULL);
3361
3362 return exprdata->exponent;
3363}
#define EPSISINT(x, eps)
Definition def.h:223
#define SCIP_INVALID
Definition def.h:206
#define SCIP_INTERVAL_INFINITY
Definition def.h:208
#define TRUE
Definition def.h:95
#define FALSE
Definition def.h:96
#define REALABS(x)
Definition def.h:210
#define SCIP_CALL(x)
Definition def.h:388
absolute expression handler
exponential expression handler
#define SIGNPOWEXPRHDLR_PRECEDENCE
Definition expr_pow.c:56
#define POWEXPRHDLR_PRECEDENCE
Definition expr_pow.c:51
void SCIPaddSquareLinearization(SCIP *scip, SCIP_Real sqrcoef, SCIP_Real refpoint, SCIP_Bool isint, SCIP_Real *lincoef, SCIP_Real *linconstant, SCIP_Bool *success)
Definition expr_pow.c:3234
static SCIP_RETCODE chooseRefpointsPow(SCIP *scip, SCIP_EXPRDATA *exprdata, SCIP_Real lb, SCIP_Real ub, SCIP_Real *refpointsunder, SCIP_Real *refpointsover, SCIP_Bool underestimate, SCIP_Bool overestimate)
Definition expr_pow.c:1300
static void estimateSignedpower(SCIP *scip, SCIP_Real exponent, SCIP_Real root, SCIP_Bool overestimate, SCIP_Real xlb, SCIP_Real xub, SCIP_Real xref, SCIP_Real xlbglobal, SCIP_Real xubglobal, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *islocal, SCIP_Bool *branchcand, SCIP_Bool *success)
Definition expr_pow.c:548
static void computeTangent(SCIP *scip, SCIP_Bool signpower, SCIP_Real exponent, SCIP_Real xref, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *success)
Definition expr_pow.c:253
#define SIGNPOWEXPRHDLR_NAME
Definition expr_pow.c:54
static void estimateRoot(SCIP *scip, SCIP_Real exponent, SCIP_Bool overestimate, SCIP_Real xlb, SCIP_Real xub, SCIP_Real xref, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *islocal, SCIP_Bool *success)
Definition expr_pow.c:1014
#define SIGNPOWEXPRHDLR_HASHKEY
Definition expr_pow.c:57
#define POWEXPRHDLR_NAME
Definition expr_pow.c:49
static SCIP_RETCODE buildPowEstimator(SCIP *scip, SCIP_EXPRDATA *exprdata, SCIP_Bool overestimate, SCIP_Real childlb, SCIP_Real childub, SCIP_Real childglb, SCIP_Real childgub, SCIP_Bool childintegral, SCIP_Real refpoint, SCIP_Real exponent, SCIP_Real *coef, SCIP_Real *constant, SCIP_Bool *success, SCIP_Bool *islocal, SCIP_Bool *branchcand)
Definition expr_pow.c:1074
static void estimateParabola(SCIP *scip, SCIP_Real exponent, SCIP_Bool overestimate, SCIP_Real xlb, SCIP_Real xub, SCIP_Real xref, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *islocal, SCIP_Bool *success)
Definition expr_pow.c:482
#define SIGN(x)
Definition expr_pow.c:70
static SCIP_RETCODE createData(SCIP *scip, SCIP_EXPRDATA **exprdata, SCIP_Real exponent)
Definition expr_pow.c:227
void SCIPaddSquareSecant(SCIP *scip, SCIP_Real sqrcoef, SCIP_Real lb, SCIP_Real ub, SCIP_Real *lincoef, SCIP_Real *linconstant, SCIP_Bool *success)
Definition expr_pow.c:3302
static SCIP_RETCODE addSignpowerRefpoints(SCIP *scip, SCIP_EXPRDATA *exprdata, SCIP_Real lb, SCIP_Real ub, SCIP_Real exponent, SCIP_Bool underestimate, SCIP_Real *refpoints)
Definition expr_pow.c:1245
#define SIGNPOW_ROOTS_KNOWN
Definition expr_pow.c:72
static void estimateHyperbolaPositive(SCIP *scip, SCIP_Real exponent, SCIP_Real root, SCIP_Bool overestimate, SCIP_Real xlb, SCIP_Real xub, SCIP_Real xref, SCIP_Real xlbglobal, SCIP_Real xubglobal, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *islocal, SCIP_Bool *branchcand, SCIP_Bool *success)
Definition expr_pow.c:696
#define POWEXPRHDLR_HASHKEY
Definition expr_pow.c:52
static SCIP_RETCODE computeSignpowerRoot(SCIP *scip, SCIP_Real *root, SCIP_Real exponent)
Definition expr_pow.c:113
static void estimateHyperbolaMixed(SCIP *scip, SCIP_Real exponent, SCIP_Bool overestimate, SCIP_Real xlb, SCIP_Real xub, SCIP_Real xref, SCIP_Real xlbglobal, SCIP_Real xubglobal, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *islocal, SCIP_Bool *branchcand, SCIP_Bool *success)
Definition expr_pow.c:909
#define SIGNPOWEXPRHDLR_DESC
Definition expr_pow.c:55
#define INITLPMAXPOWVAL
Definition expr_pow.c:59
#define POWEXPRHDLR_DESC
Definition expr_pow.c:50
static void addTangentRefpoints(SCIP *scip, SCIP_Real exponent, SCIP_Real lb, SCIP_Real ub, SCIP_Real *refpoints)
Definition expr_pow.c:1207
static SCIP_RETCODE computeHyperbolaRoot(SCIP *scip, SCIP_Real *root, SCIP_Real exponent)
Definition expr_pow.c:181
static SCIP_Real signpow_roots[SIGNPOW_ROOTS_KNOWN+1]
Definition expr_pow.c:78
static void computeSecant(SCIP *scip, SCIP_Bool signpower, SCIP_Real exponent, SCIP_Real xlb, SCIP_Real xub, SCIP_Real *constant, SCIP_Real *slope, SCIP_Bool *success)
Definition expr_pow.c:302
power and signed power expression handlers
product expression handler
sum expression handler
constant value expression handler
SCIP_RETCODE SCIPcreateExprProduct(SCIP *scip, SCIP_EXPR **expr, int nchildren, SCIP_EXPR **children, SCIP_Real coefficient, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
SCIP_RETCODE SCIPcreateExprAbs(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_abs.c:528
SCIP_RETCODE SCIPcreateExprSignpower(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_Real exponent, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_pow.c:3198
SCIP_Bool SCIPisExprExp(SCIP *scip, SCIP_EXPR *expr)
Definition expr_exp.c:528
SCIP_RETCODE SCIPcreateExprExp(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_exp.c:510
SCIP_Bool SCIPisExprSignpower(SCIP *scip, SCIP_EXPR *expr)
Definition expr_pow.c:3223
SCIP_RETCODE SCIPcreateExprSum(SCIP *scip, SCIP_EXPR **expr, int nchildren, SCIP_EXPR **children, SCIP_Real *coefficients, SCIP_Real constant, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_sum.c:1079
SCIP_RETCODE SCIPcreateExprValue(SCIP *scip, SCIP_EXPR **expr, SCIP_Real value, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_value.c:270
SCIP_RETCODE SCIPcreateExprPow(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_Real exponent, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition expr_pow.c:3174
SCIP_RETCODE SCIPincludeExprhdlrSignpower(SCIP *scip)
Definition expr_pow.c:3145
SCIP_RETCODE SCIPincludeExprhdlrPow(SCIP *scip)
Definition expr_pow.c:3110
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
SCIP_VERBLEVEL SCIPgetVerbLevel(SCIP *scip)
#define SCIPdebugMsgPrint
#define SCIPdebugMsg
#define SCIPisFinite(x)
Definition pub_misc.h:1901
SCIP_RETCODE SCIPaddRealParam(SCIP *scip, const char *name, const char *desc, SCIP_Real *valueptr, SCIP_Bool isadvanced, SCIP_Real defaultvalue, SCIP_Real minvalue, SCIP_Real maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition scip_param.c:139
const char * SCIPexprhdlrGetName(SCIP_EXPRHDLR *exprhdlr)
Definition expr.c:532
void SCIPexprhdlrSetCompare(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:460
void SCIPexprhdlrSetIntegrality(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:438
void SCIPexprhdlrSetCurvature(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:416
void SCIPexprhdlrSetParse(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:405
SCIP_EXPRHDLRDATA * SCIPexprhdlrGetData(SCIP_EXPRHDLR *exprhdlr)
Definition expr.c:562
void SCIPexprhdlrSetIntEval(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:486
void SCIPexprhdlrSetMonotonicity(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:427
void SCIPexprhdlrSetReverseProp(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:508
void SCIPexprhdlrSetHash(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:449
SCIP_RETCODE SCIPincludeExprhdlr(SCIP *scip, SCIP_EXPRHDLR **exprhdlr, const char *name, const char *desc, unsigned int precedence, SCIP_DECL_EXPREVAL((*eval)), SCIP_EXPRHDLRDATA *data)
Definition scip_expr.c:814
void SCIPexprhdlrSetSimplify(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:497
void SCIPexprhdlrSetDiff(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRBWDIFF((*bwdiff)), SCIP_DECL_EXPRFWDIFF((*fwdiff)),)
Definition expr.c:471
void SCIPexprhdlrSetCopyFreeHdlr(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCOPYHDLR((*copyhdlr)),)
Definition expr.c:368
SCIP_EXPRHDLR * SCIPgetExprhdlrPower(SCIP *scip)
Definition scip_expr.c:915
void SCIPexprhdlrSetPrint(SCIP_EXPRHDLR *exprhdlr,)
Definition expr.c:394
SCIP_EXPRHDLR * SCIPfindExprhdlr(SCIP *scip, const char *name)
Definition scip_expr.c:859
void SCIPexprhdlrSetCopyFreeData(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCOPYDATA((*copydata)),)
Definition expr.c:381
void SCIPexprhdlrSetEstimate(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRINITESTIMATES((*initestimates)),)
Definition expr.c:519
SCIP_RETCODE SCIPcreateExpr(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPRHDLR *exprhdlr, SCIP_EXPRDATA *exprdata, int nchildren, SCIP_EXPR **children, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition scip_expr.c:964
SCIP_RETCODE SCIPappendExprChild(SCIP *scip, SCIP_EXPR *expr, SCIP_EXPR *child)
Definition scip_expr.c:1220
void SCIPexprSetData(SCIP_EXPR *expr, SCIP_EXPRDATA *exprdata)
Definition expr.c:3849
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition expr.c:3801
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition expr_pow.c:3351
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1454
SCIP_EXPRCURV SCIPexprcurvPowerInv(SCIP_INTERVAL basebounds, SCIP_Real exponent, SCIP_EXPRCURV powercurv)
Definition exprcurv.c:208
SCIP_Bool SCIPisExprSum(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1443
SCIP_Bool SCIPexprIsIntegral(SCIP_EXPR *expr)
Definition expr.c:4020
SCIP_Real * SCIPgetCoefsExprSum(SCIP_EXPR *expr)
Definition expr_sum.c:1166
SCIP_Bool SCIPisExprValue(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1432
int SCIPcompareExpr(SCIP *scip, SCIP_EXPR *expr1, SCIP_EXPR *expr2)
Definition scip_expr.c:1724
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition scip_expr.c:1407
SCIP_Real SCIPexprGetDot(SCIP_EXPR *expr)
Definition expr.c:3915
SCIP_Bool SCIPisExprVar(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1421
SCIP_EXPRDATA * SCIPexprGetData(SCIP_EXPR *expr)
Definition expr.c:3834
SCIP_RETCODE SCIPparseExpr(SCIP *scip, SCIP_EXPR **expr, const char *exprstr, const char **finalpos, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition scip_expr.c:1370
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition scip_expr.c:1476
SCIP_Real SCIPgetValueExprValue(SCIP_EXPR *expr)
Definition expr_value.c:294
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1465
SCIP_Real SCIPexprGetEvalValue(SCIP_EXPR *expr)
Definition expr.c:3875
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition expr.c:3811
SCIP_Real SCIPgetConstantExprSum(SCIP_EXPR *expr)
Definition expr_sum.c:1181
SCIP_VAR * SCIPgetVarExprVar(SCIP_EXPR *expr)
Definition expr_var.c:416
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition expr.c:3957
void SCIPcaptureExpr(SCIP_EXPR *expr)
Definition scip_expr.c:1399
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition scip_expr.c:1707
SCIP_EXPRHDLR * SCIPexprGetHdlr(SCIP_EXPR *expr)
Definition expr.c:3824
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
SCIP_Bool SCIPintervalIsEntire(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSignPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
void SCIPintervalUnify(SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_INTERVAL operand2)
void SCIPintervalSet(SCIP_INTERVAL *resultant, SCIP_Real value)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalPowerScalarInverse(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL basedomain, SCIP_Real exponent, SCIP_INTERVAL image)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPintervalReciprocal(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalPowerScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
#define SCIPallocClearBlockMemory(scip, ptr)
Definition scip_mem.h:91
#define SCIPallocBufferArray(scip, ptr, num)
Definition scip_mem.h:124
#define SCIPfreeBufferArray(scip, ptr)
Definition scip_mem.h:136
#define SCIPfreeBlockMemory(scip, ptr)
Definition scip_mem.h:108
#define SCIPallocBlockMemory(scip, ptr)
Definition scip_mem.h:89
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisIntegral(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisPositive(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisFeasZero(SCIP *scip, SCIP_Real val)
SCIP_Real SCIPfloor(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisNegative(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Real SCIPepsilon(SCIP *scip)
SCIP_Bool SCIPparseReal(SCIP *scip, const char *str, SCIP_Real *value, char **endptr)
SCIP_Bool SCIPvarIsBinary(SCIP_VAR *var)
Definition var.c:17421
return SCIP_OKAY
int c
assert(minobj< SCIPgetCutoffbound(scip))
#define NULL
Definition lpi_spx1.cpp:161
public functions to work with algebraic expressions
#define SCIPerrorMessage
Definition pub_message.h:64
#define SCIPdebugPrintf
Definition pub_message.h:99
SCIP_Real sup
SCIP_Real inf
#define MAX(x, y)
Definition tclique_def.h:92
#define SCIP_DECL_EXPR_OWNERCREATE(x)
Definition type_expr.h:140
#define SCIP_DECL_EXPRREVERSEPROP(x)
Definition type_expr.h:654
#define SCIP_DECL_EXPRINITESTIMATES(x)
Definition type_expr.h:605
#define SCIP_DECL_EXPRBWFWDIFF(x)
Definition type_expr.h:517
#define SCIP_DECL_EXPRCURVATURE(x)
Definition type_expr.h:337
struct SCIP_ExprhdlrData SCIP_EXPRHDLRDATA
Definition type_expr.h:192
struct SCIP_ExprData SCIP_EXPRDATA
Definition type_expr.h:53
#define SCIP_DECL_EXPRFREEDATA(x)
Definition type_expr.h:265
@ SCIP_EXPRCURV_CONVEX
Definition type_expr.h:60
@ SCIP_EXPRCURV_UNKNOWN
Definition type_expr.h:59
@ SCIP_EXPRCURV_CONCAVE
Definition type_expr.h:61
#define SCIP_EXPR_MAXINITESTIMATES
Definition type_expr.h:195
#define SCIP_DECL_EXPRPARSE(x)
Definition type_expr.h:309
#define SCIP_DECL_EXPRBWDIFF(x)
Definition type_expr.h:446
#define SCIP_DECL_EXPRINTEVAL(x)
Definition type_expr.h:536
#define SCIP_DECL_EXPRMONOTONICITY(x)
Definition type_expr.h:355
#define SCIP_EXPRITER_VISITINGCHILD
Definition type_expr.h:677
@ SCIP_MONOTONE_UNKNOWN
Definition type_expr.h:68
@ SCIP_MONOTONE_INC
Definition type_expr.h:69
@ SCIP_MONOTONE_DEC
Definition type_expr.h:70
#define SCIP_DECL_EXPRCOMPARE(x)
Definition type_expr.h:407
#define SCIP_DECL_EXPRSIMPLIFY(x)
Definition type_expr.h:629
#define SCIP_DECL_EXPREVAL(x)
Definition type_expr.h:423
#define SCIP_DECL_EXPRFWDIFF(x)
Definition type_expr.h:477
#define SCIP_DECL_EXPRHASH(x)
Definition type_expr.h:388
#define SCIP_DECL_EXPRCOPYHDLR(x)
Definition type_expr.h:207
#define SCIP_DECL_EXPRPRINT(x)
Definition type_expr.h:286
#define SCIP_DECL_EXPRFREEHDLR(x)
Definition type_expr.h:221
#define SCIP_DECL_EXPRINTEGRALITY(x)
Definition type_expr.h:372
#define SCIP_EXPRITER_VISITEDCHILD
Definition type_expr.h:678
#define SCIP_DECL_EXPRCOPYDATA(x)
Definition type_expr.h:246
#define SCIP_EXPRITER_LEAVEEXPR
Definition type_expr.h:679
#define SCIP_DECL_EXPRESTIMATE(x)
Definition type_expr.h:572
#define SCIP_EXPRITER_ENTEREXPR
Definition type_expr.h:676
@ SCIP_VERBLEVEL_NONE
@ SCIP_READERROR
@ SCIP_ERROR
enum SCIP_Retcode SCIP_RETCODE