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_exp.c 26 * @ingroup DEFPLUGINS_EXPR 27 * @brief exponential expression handler 28 * @author Stefan Vigerske 29 * @author Benjamin Mueller 30 * @author Ksenia Bestuzheva 31 */ 32 33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 34 35 #define _USE_MATH_DEFINES /* to get M_E on Windows */ /*lint !750 */ 36 37 #include <string.h> 38 #include <math.h> 39 40 #include "scip/expr_exp.h" 41 #include "scip/expr_value.h" 42 43 #define EXPRHDLR_NAME "exp" 44 #define EXPRHDLR_DESC "exponential expression" 45 #define EXPRHDLR_PRECEDENCE 85000 46 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(10181.0) 47 48 /* 49 * Data structures 50 */ 51 52 /* 53 * Local methods 54 */ 55 56 /** computes coefficients of secant of an exponential term */ 57 static 58 void addExpSecant( 59 SCIP* scip, /**< SCIP data structure */ 60 SCIP_Real lb, /**< lower bound on variable */ 61 SCIP_Real ub, /**< upper bound on variable */ 62 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */ 63 SCIP_Real* linconstant, /**< buffer to add constant of secant */ 64 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */ 65 ) 66 { 67 SCIP_Real coef; 68 SCIP_Real constant; 69 70 assert(scip != NULL); 71 assert(!SCIPisInfinity(scip, lb)); 72 assert(!SCIPisInfinity(scip, -ub)); 73 assert(SCIPisLE(scip, lb, ub)); 74 assert(lincoef != NULL); 75 assert(linconstant != NULL); 76 assert(success != NULL); 77 78 if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) ) 79 { 80 /* unboundedness */ 81 *success = FALSE; 82 return; 83 } 84 85 /* if lb and ub are too close use a safe secant */ 86 if( SCIPisEQ(scip, lb, ub) ) 87 { 88 coef = 0.0; 89 constant = exp(ub); 90 } 91 else 92 { 93 coef = (exp(ub) - exp(lb)) / (ub - lb); 94 constant = exp(ub) - coef * ub; 95 } 96 97 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) ) 98 { 99 *success = FALSE; 100 return; 101 } 102 103 *lincoef += coef; 104 *linconstant += constant; 105 } 106 107 /** computes coefficients of linearization of an exponential term in a reference point */ 108 static 109 void addExpLinearization( 110 SCIP* scip, /**< SCIP data structure */ 111 SCIP_Real refpoint, /**< point for which to compute value of linearization */ 112 SCIP_Bool isint, /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */ 113 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */ 114 SCIP_Real* linconstant, /**< buffer to add constant of secant */ 115 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */ 116 ) 117 { 118 SCIP_Real constant; 119 SCIP_Real coef; 120 121 assert(scip != NULL); 122 assert(lincoef != NULL); 123 assert(linconstant != NULL); 124 assert(success != NULL); 125 126 if( SCIPisInfinity(scip, REALABS(refpoint)) ) 127 { 128 *success = FALSE; 129 return; 130 } 131 132 if( !isint || SCIPisIntegral(scip, refpoint) ) 133 { 134 coef = exp(refpoint); 135 constant = exp(refpoint) * (1.0 - refpoint); 136 } 137 else 138 { 139 /* exp(x) -> secant between f=floor(refpoint) and f+1 = ((e-1)*e^f) * x + e^f - f * ((e-1)*e^f) */ 140 SCIP_Real f; 141 142 f = SCIPfloor(scip, refpoint); 143 144 coef = (M_E - 1.0) * exp(f); 145 constant = exp(f) - f * coef; 146 } 147 148 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) ) 149 { 150 *success = FALSE; 151 return; 152 } 153 154 *lincoef += coef; 155 *linconstant += constant; 156 } 157 158 /* 159 * Callback methods of expression handler 160 */ 161 162 /** simplifies an exp expression 163 * 164 * Evaluates the exponential function when its child is a value expression. 165 * 166 * TODO: exp(log(*)) = * 167 */ 168 static 169 SCIP_DECL_EXPRSIMPLIFY(simplifyExp) 170 { 171 SCIP_EXPR* child; 172 173 assert(scip != NULL); 174 assert(expr != NULL); 175 assert(simplifiedexpr != NULL); 176 assert(SCIPexprGetNChildren(expr) == 1); 177 178 child = SCIPexprGetChildren(expr)[0]; 179 assert(child != NULL); 180 181 /**! [SnippetExprSimplifyExp] */ 182 /* check for value expression */ 183 if( SCIPisExprValue(scip, child) ) 184 { 185 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, exp(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) ); 186 } 187 else 188 { 189 *simplifiedexpr = expr; 190 191 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */ 192 SCIPcaptureExpr(*simplifiedexpr); 193 } 194 /**! [SnippetExprSimplifyExp] */ 195 196 return SCIP_OKAY; 197 } 198 199 /** expression handler copy callback */ 200 static 201 SCIP_DECL_EXPRCOPYHDLR(copyhdlrExp) 202 { /*lint --e{715}*/ 203 SCIP_CALL( SCIPincludeExprhdlrExp(scip) ); 204 205 return SCIP_OKAY; 206 } 207 208 /** expression data copy callback */ 209 static 210 SCIP_DECL_EXPRCOPYDATA(copydataExp) 211 { /*lint --e{715}*/ 212 assert(targetexprdata != NULL); 213 assert(sourceexpr != NULL); 214 assert(SCIPexprGetData(sourceexpr) == NULL); 215 216 *targetexprdata = NULL; 217 218 return SCIP_OKAY; 219 } 220 221 /** expression data free callback */ 222 static 223 SCIP_DECL_EXPRFREEDATA(freedataExp) 224 { /*lint --e{715}*/ 225 assert(expr != NULL); 226 227 SCIPexprSetData(expr, NULL); 228 229 return SCIP_OKAY; 230 } 231 232 /** expression parse callback */ 233 static 234 SCIP_DECL_EXPRPARSE(parseExp) 235 { /*lint --e{715}*/ 236 SCIP_EXPR* childexpr; 237 238 assert(expr != NULL); 239 240 /**! [SnippetExprParseExp] */ 241 /* parse child expression from remaining string */ 242 SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) ); 243 assert(childexpr != NULL); 244 245 /* create exponential expression */ 246 SCIP_CALL( SCIPcreateExprExp(scip, expr, childexpr, ownercreate, ownercreatedata) ); 247 assert(*expr != NULL); 248 249 /* release child expression since it has been captured by the exponential expression */ 250 SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) ); 251 252 *success = TRUE; 253 /**! [SnippetExprParseExp] */ 254 255 return SCIP_OKAY; 256 } 257 258 /** expression point evaluation callback */ 259 static 260 SCIP_DECL_EXPREVAL(evalExp) 261 { /*lint --e{715}*/ 262 assert(expr != NULL); 263 assert(SCIPexprGetData(expr) == NULL); 264 assert(SCIPexprGetNChildren(expr) == 1); 265 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/ 266 267 *val = exp(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0])); 268 269 return SCIP_OKAY; 270 } 271 272 /** expression derivative evaluation callback */ 273 static 274 SCIP_DECL_EXPRBWDIFF(bwdiffExp) 275 { /*lint --e{715}*/ 276 assert(expr != NULL); 277 assert(childidx == 0); 278 assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[0])); 279 assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/ 280 281 *val = SCIPexprGetEvalValue(expr); 282 283 return SCIP_OKAY; 284 } 285 286 /** expression interval evaluation callback */ 287 static 288 SCIP_DECL_EXPRINTEVAL(intevalExp) 289 { /*lint --e{715}*/ 290 SCIP_INTERVAL childinterval; 291 292 assert(expr != NULL); 293 assert(SCIPexprGetData(expr) == NULL); 294 assert(SCIPexprGetNChildren(expr) == 1); 295 296 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]); 297 298 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) ) 299 SCIPintervalSetEmpty(interval); 300 else 301 SCIPintervalExp(SCIP_INTERVAL_INFINITY, interval, childinterval); 302 303 return SCIP_OKAY; 304 } 305 306 /** expression estimator callback */ 307 static 308 SCIP_DECL_EXPRESTIMATE(estimateExp) 309 { /*lint --e{715}*/ 310 assert(scip != NULL); 311 assert(expr != NULL); 312 assert(SCIPexprGetNChildren(expr) == 1); 313 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 314 assert(coefs != NULL); 315 assert(constant != NULL); 316 assert(islocal != NULL); 317 assert(branchcand != NULL); 318 assert(*branchcand == TRUE); 319 assert(success != NULL); 320 321 *success = TRUE; 322 *coefs = 0.0; 323 *constant = 0.0; 324 325 if( overestimate ) 326 { 327 addExpSecant(scip, localbounds[0].inf, localbounds[0].sup, coefs, constant, success); 328 *islocal = TRUE; /* secants are only valid locally */ 329 } 330 else 331 { 332 addExpLinearization(scip, refpoint[0], SCIPexprIsIntegral(SCIPexprGetChildren(expr)[0]), coefs, constant, success); 333 *islocal = FALSE; /* linearization are globally valid */ 334 *branchcand = FALSE; 335 } 336 337 return SCIP_OKAY; 338 } 339 340 /** initital estimates callback for an exponential expression */ 341 static 342 SCIP_DECL_EXPRINITESTIMATES(initestimatesExp) 343 { 344 SCIP_Real refpointsunder[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID}; 345 SCIP_Bool overest[4] = {FALSE, FALSE, FALSE, TRUE}; 346 SCIP_EXPR* child; 347 SCIP_Real lb; 348 SCIP_Real ub; 349 SCIP_Bool success; 350 int i; 351 352 assert(scip != NULL); 353 assert(expr != NULL); 354 assert(SCIPexprGetNChildren(expr) == 1); 355 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 356 357 /* get expression data */ 358 child = SCIPexprGetChildren(expr)[0]; 359 assert(child != NULL); 360 361 lb = bounds[0].inf; 362 ub = bounds[0].sup; 363 364 if( !overestimate ) 365 { 366 SCIP_Real lbfinite; 367 SCIP_Real ubfinite; 368 369 /* make bounds finite */ 370 lbfinite = SCIPisInfinity(scip, -lb) ? MIN(-5.0, ub - 0.1 * REALABS(ub)) : lb; /*lint !e666*/ 371 ubfinite = SCIPisInfinity(scip, ub) ? MAX( 3.0, lb + 0.1 * REALABS(lb)) : ub; /*lint !e666*/ 372 373 refpointsunder[0] = (7.0 * lbfinite + ubfinite) / 8.0; 374 refpointsunder[1] = (lbfinite + ubfinite) / 2.0; 375 refpointsunder[2] = (lbfinite + 7.0 * ubfinite) / 8.0; 376 } 377 378 *nreturned = 0; 379 for( i = 0; i < 4; ++i ) 380 { 381 if( !overest[i] && overestimate ) 382 continue; 383 384 if( overest[i] && (!overestimate || SCIPisInfinity(scip, ub) || SCIPisInfinity(scip, -lb)) ) 385 continue; 386 387 assert(overest[i] || (SCIPisLE(scip, refpointsunder[i], ub) && SCIPisGE(scip, refpointsunder[i], lb))); /*lint !e661*/ 388 389 coefs[*nreturned][0] = 0.0; 390 constant[*nreturned] = 0.0; 391 392 success = TRUE; 393 394 if( !overest[i] ) 395 { 396 assert(i < 3); 397 /* coverity[overrun] */ 398 addExpLinearization(scip, refpointsunder[i], SCIPexprIsIntegral(child), coefs[*nreturned], &constant[*nreturned], &success); /*lint !e661*/ 399 } 400 else 401 addExpSecant(scip, lb, ub, coefs[*nreturned], &constant[*nreturned], &success); 402 403 if( success ) 404 ++*nreturned; 405 } 406 407 return SCIP_OKAY; 408 } 409 410 /** expression reverse propagation callback */ 411 static 412 SCIP_DECL_EXPRREVERSEPROP(reversepropExp) 413 { /*lint --e{715}*/ 414 assert(scip != NULL); 415 assert(expr != NULL); 416 assert(SCIPexprGetNChildren(expr) == 1); 417 assert(SCIPintervalGetInf(bounds) >= 0.0); 418 419 if( SCIPintervalGetSup(bounds) <= 0.0 ) 420 { 421 *infeasible = TRUE; 422 return SCIP_OKAY; 423 } 424 425 /* f = exp(c0) -> c0 = log(f) */ 426 SCIPintervalLog(SCIP_INTERVAL_INFINITY, childrenbounds, bounds); 427 428 return SCIP_OKAY; 429 } 430 431 /** expression hash callback */ 432 static 433 SCIP_DECL_EXPRHASH(hashExp) 434 { /*lint --e{715}*/ 435 assert(scip != NULL); 436 assert(expr != NULL); 437 assert(SCIPexprGetNChildren(expr) == 1); 438 assert(hashkey != NULL); 439 assert(childrenhashes != NULL); 440 441 *hashkey = EXPRHDLR_HASHKEY; 442 *hashkey ^= childrenhashes[0]; 443 444 return SCIP_OKAY; 445 } 446 447 /** expression curvature detection callback */ 448 static 449 SCIP_DECL_EXPRCURVATURE(curvatureExp) 450 { /*lint --e{715}*/ 451 assert(scip != NULL); 452 assert(expr != NULL); 453 assert(childcurv != NULL); 454 assert(success != NULL); 455 assert(SCIPexprGetNChildren(expr) == 1); 456 457 /* expression is convex if child is convex; expression cannot be concave or linear */ 458 if( exprcurvature == SCIP_EXPRCURV_CONVEX ) 459 { 460 *success = TRUE; 461 *childcurv = SCIP_EXPRCURV_CONVEX; 462 } 463 else 464 *success = FALSE; 465 466 return SCIP_OKAY; 467 } 468 469 /** expression monotonicity detection callback */ 470 static 471 SCIP_DECL_EXPRMONOTONICITY(monotonicityExp) 472 { /*lint --e{715}*/ 473 assert(scip != NULL); 474 assert(expr != NULL); 475 assert(result != NULL); 476 assert(childidx == 0); 477 478 *result = SCIP_MONOTONE_INC; 479 480 return SCIP_OKAY; 481 } 482 483 /** creates the handler for exponential expressions and includes it into SCIP */ 484 SCIP_RETCODE SCIPincludeExprhdlrExp( 485 SCIP* scip /**< SCIP data structure */ 486 ) 487 { 488 SCIP_EXPRHDLR* exprhdlr; 489 490 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, 491 EXPRHDLR_PRECEDENCE, evalExp, NULL) ); 492 assert(exprhdlr != NULL); 493 494 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrExp, NULL); 495 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataExp, freedataExp); 496 SCIPexprhdlrSetSimplify(exprhdlr, simplifyExp); 497 SCIPexprhdlrSetParse(exprhdlr, parseExp); 498 SCIPexprhdlrSetIntEval(exprhdlr, intevalExp); 499 SCIPexprhdlrSetEstimate(exprhdlr, initestimatesExp, estimateExp); 500 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropExp); 501 SCIPexprhdlrSetHash(exprhdlr, hashExp); 502 SCIPexprhdlrSetDiff(exprhdlr, bwdiffExp, NULL, NULL); 503 SCIPexprhdlrSetCurvature(exprhdlr, curvatureExp); 504 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityExp); 505 506 return SCIP_OKAY; 507 } 508 509 /** creates an exponential expression */ 510 SCIP_RETCODE SCIPcreateExprExp( 511 SCIP* scip, /**< SCIP data structure */ 512 SCIP_EXPR** expr, /**< pointer where to store expression */ 513 SCIP_EXPR* child, /**< single child */ 514 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 515 void* ownercreatedata /**< data to pass to ownercreate */ 516 ) 517 { 518 assert(expr != NULL); 519 assert(child != NULL); 520 assert(SCIPfindExprhdlr(scip, EXPRHDLR_NAME) != NULL); 521 522 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, EXPRHDLR_NAME), NULL, 1, &child, ownercreate, ownercreatedata) ); 523 524 return SCIP_OKAY; 525 } 526 527 /** indicates whether expression is of exp-type */ /*lint -e{715}*/ 528 SCIP_Bool SCIPisExprExp( 529 SCIP* scip, /**< SCIP data structure */ 530 SCIP_EXPR* expr /**< expression */ 531 ) 532 { /*lint --e{715}*/ 533 assert(expr != NULL); 534 535 return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0; 536 } 537