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