Scippy

SCIP

Solving Constraint Integer Programs

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