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_abs.c 26 * @ingroup DEFPLUGINS_EXPR 27 * @brief absolute expression handler 28 * @author Stefan Vigerske 29 * @author Benjamin Mueller 30 */ 31 32 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 33 34 #include <string.h> 35 36 #include "scip/expr_value.h" 37 #include "scip/expr_abs.h" 38 #include "scip/expr.h" 39 40 #define EXPRHDLR_NAME "abs" 41 #define EXPRHDLR_DESC "absolute value expression" 42 #define EXPRHDLR_PRECEDENCE 70000 43 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(7187.0) 44 45 /* 46 * Data structures 47 */ 48 49 /* 50 * Local methods 51 */ 52 53 /** computes both tangent underestimates and secant */ 54 static 55 SCIP_RETCODE computeCutsAbs( 56 SCIP* scip, /**< SCIP data structure */ 57 SCIP_INTERVAL bounds, /**< bounds of child */ 58 SCIP_Bool overestimate, /**< whether the expression shall be overestimated or underestimated */ 59 SCIP_Real** coefs, /**< buffer to store coefficients of computed estimators */ 60 SCIP_Real* constant, /**< buffer to store constant of computed estimators */ 61 int* nreturned /**< buffer to store number of estimators that have been computed */ 62 ) 63 { 64 assert(scip != NULL); 65 66 *nreturned = 0; 67 68 /**! [SnippetExprInitestimatesAbs] */ 69 if( !overestimate ) 70 { 71 /* compute left tangent -x <= z */ 72 coefs[*nreturned][0] = -1.0; 73 constant[*nreturned] = 0.0; 74 (*nreturned)++; 75 76 /* compute right tangent x <= z */ 77 coefs[*nreturned][0] = 1.0; 78 constant[*nreturned] = 0.0; 79 (*nreturned)++; 80 } 81 82 /* compute secant */ 83 if( overestimate ) 84 { 85 SCIP_Real lb; 86 SCIP_Real ub; 87 88 lb = bounds.inf; 89 ub = bounds.sup; 90 91 /* it does not make sense to add a cut if child variable is unbounded or fixed */ 92 if( !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) && !SCIPisEQ(scip, lb, ub) ) 93 { 94 if( !SCIPisPositive(scip, ub) ) 95 { 96 /* z = -x, so add z <= -x here (-x <= z is the underestimator that is added above) */ 97 coefs[*nreturned][0] = -1.0; 98 constant[*nreturned] = 0.0; 99 (*nreturned)++; 100 } 101 else if( !SCIPisNegative(scip, lb) ) 102 { 103 /* z = x, so add z <= x here (x <= z is the underestimator that is added above) */ 104 coefs[*nreturned][0] = 1.0; 105 constant[*nreturned] = 0.0; 106 (*nreturned)++; 107 } 108 else 109 { 110 /* z = abs(x), x still has mixed sign */ 111 SCIP_Real alpha; 112 113 /* let alpha = (|ub|-|lb|) / (ub-lb) then the resulting secant looks like 114 * 115 * z - |ub| <= alpha * (x - ub) <=> z <= alpha * x + |ub| - alpha * ub 116 */ 117 alpha = (REALABS(ub) - REALABS(lb)) / (ub - lb); 118 119 coefs[*nreturned][0] = alpha; 120 constant[*nreturned] = REALABS(ub) - alpha * ub; 121 (*nreturned)++; 122 } 123 } 124 } 125 /**! [SnippetExprInitestimatesAbs] */ 126 127 return SCIP_OKAY; 128 } 129 130 131 /* 132 * Callback methods of expression handler 133 */ 134 135 /** simplifies an abs expression 136 * 137 * Evaluates the absolute value function when its child is a value expression. 138 * 139 * TODO: abs(*) = * if * >= 0 or - * if * < 0 140 */ 141 static 142 SCIP_DECL_EXPRSIMPLIFY(simplifyAbs) 143 { /*lint --e{715}*/ 144 SCIP_EXPR* child; 145 146 assert(scip != NULL); 147 assert(expr != NULL); 148 assert(simplifiedexpr != NULL); 149 assert(SCIPexprGetNChildren(expr) == 1); 150 151 child = SCIPexprGetChildren(expr)[0]; 152 assert(child != NULL); 153 154 /* check for value expression */ 155 if( SCIPisExprValue(scip, child) ) 156 { 157 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, REALABS(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) ); 158 } 159 else 160 { 161 *simplifiedexpr = expr; 162 163 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */ 164 SCIPcaptureExpr(*simplifiedexpr); 165 } 166 167 return SCIP_OKAY; 168 } 169 170 static 171 SCIP_DECL_EXPRCOPYHDLR(copyhdlrAbs) 172 { /*lint --e{715}*/ 173 SCIP_CALL( SCIPincludeExprhdlrAbs(scip) ); 174 175 return SCIP_OKAY; 176 } 177 178 static 179 SCIP_DECL_EXPRPARSE(parseAbs) 180 { /*lint --e{715}*/ 181 SCIP_EXPR* childexpr; 182 183 assert(expr != NULL); 184 185 /* parse child expression from remaining string */ 186 SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) ); 187 assert(childexpr != NULL); 188 189 /* create absolute expression */ 190 SCIP_CALL( SCIPcreateExprAbs(scip, expr, childexpr, ownercreate, ownercreatedata) ); 191 assert(*expr != NULL); 192 193 /* release child expression since it has been captured by the absolute expression */ 194 SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) ); 195 196 *success = TRUE; 197 198 return SCIP_OKAY; 199 } 200 201 /** expression point evaluation callback */ 202 static 203 SCIP_DECL_EXPREVAL(evalAbs) 204 { /*lint --e{715}*/ 205 assert(expr != NULL); 206 assert(SCIPexprGetNChildren(expr) == 1); 207 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/ 208 209 *val = REALABS(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0])); 210 211 return SCIP_OKAY; 212 } 213 214 215 /** expression derivative evaluation callback */ 216 static 217 SCIP_DECL_EXPRBWDIFF(bwdiffAbs) 218 { /*lint --e{715}*/ 219 SCIP_EXPR* child; 220 221 assert(expr != NULL); 222 assert(childidx == 0); 223 assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/ 224 225 child = SCIPexprGetChildren(expr)[0]; 226 assert(child != NULL); 227 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0); 228 229 *val = (SCIPexprGetEvalValue(child) >= 0.0) ? 1.0 : -1.0; 230 231 return SCIP_OKAY; 232 } 233 234 /** expression interval evaluation callback */ 235 static 236 SCIP_DECL_EXPRINTEVAL(intevalAbs) 237 { /*lint --e{715}*/ 238 SCIP_INTERVAL childinterval; 239 240 assert(expr != NULL); 241 assert(SCIPexprGetNChildren(expr) == 1); 242 243 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]); 244 245 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) ) 246 SCIPintervalSetEmpty(interval); 247 else 248 SCIPintervalAbs(SCIP_INTERVAL_INFINITY, interval, childinterval); 249 250 return SCIP_OKAY; 251 } 252 253 /** expression estimator callback */ 254 static 255 SCIP_DECL_EXPRESTIMATE(estimateAbs) 256 { /*lint --e{715}*/ 257 assert(scip != NULL); 258 assert(expr != NULL); 259 assert(SCIPexprGetNChildren(expr) == 1); 260 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 261 assert(coefs != NULL); 262 assert(constant != NULL); 263 assert(islocal != NULL); 264 assert(branchcand != NULL); 265 assert(*branchcand == TRUE); 266 assert(success != NULL); 267 268 SCIPdebugMsg(scip, "%sestimate |child| over locdom=[%g,%g] glbdom=[%g,%g]\n", overestimate ? "over" : "under", 269 localbounds[0].inf, localbounds[0].sup, globalbounds[0].inf, globalbounds[0].sup); 270 271 /**! [SnippetExprEstimateAbs] */ 272 if( !overestimate ) 273 { 274 *constant = 0.0; 275 276 if( refpoint[0] <= 0.0 ) 277 *coefs = -1.0; 278 else 279 *coefs = 1.0; 280 281 *islocal = FALSE; 282 *branchcand = FALSE; 283 } 284 else 285 { 286 /* overestimator */ 287 SCIP_Real lb; 288 SCIP_Real ub; 289 290 lb = localbounds[0].inf; 291 ub = localbounds[0].sup; 292 293 if( !SCIPisPositive(scip, ub) ) 294 { 295 /* |x| = -x */ 296 *coefs = -1.0; 297 *constant = 0.0; 298 *islocal = SCIPisPositive(scip, globalbounds[0].sup); 299 *branchcand = FALSE; 300 } 301 else if( !SCIPisNegative(scip, lb) ) 302 { 303 /* |x| = x */ 304 *coefs = 1.0; 305 *constant = 0.0; 306 *islocal = SCIPisNegative(scip, globalbounds[0].inf); 307 *branchcand = FALSE; 308 } 309 else if( !SCIPisRelEQ(scip, lb, -ub) ) 310 { 311 /* |x| with x having mixed sign and ub+lb does not cancel out -> secant */ 312 SCIP_Real alpha; 313 314 assert(lb < 0.0); 315 assert(ub > 0.0); 316 317 /* let alpha = (|ub|-|lb|) / (ub-lb) = (ub+lb)/(ub-lb) 318 * then the resulting secant is -lb + alpha * (x - lb) = -lb - alpha*lb + alpha*x 319 */ 320 alpha = (ub + lb) / (ub - lb); 321 322 *coefs = alpha; 323 *constant = -lb - alpha * lb; 324 *islocal = TRUE; 325 } 326 else if( lb == -ub ) /*lint !e777*/ 327 { 328 /* alpha = 0 */ 329 *coefs = 0.0; 330 *constant = -lb; 331 *islocal = TRUE; 332 } 333 else 334 { 335 *success = FALSE; 336 return SCIP_OKAY; 337 } 338 } 339 /**! [SnippetExprEstimateAbs] */ 340 341 SCIPdebugMsg(scip, "-> %g * <child> %+g, local=%u branchcand=%u\n", *coefs, *constant, *islocal, *branchcand); 342 343 *success = TRUE; 344 345 return SCIP_OKAY; 346 } 347 348 /** expression estimate initialization callback */ 349 static 350 SCIP_DECL_EXPRINITESTIMATES(initEstimatesAbs) 351 { /*lint --e{715}*/ 352 assert(expr != NULL); 353 assert(SCIPexprGetNChildren(expr) == 1); 354 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 355 356 /* compute initial cuts */ 357 SCIP_CALL( computeCutsAbs(scip, bounds[0], overestimate, coefs, constant, nreturned) ); 358 359 return SCIP_OKAY; 360 } 361 362 /** expression reverse propagation callback */ 363 static 364 SCIP_DECL_EXPRREVERSEPROP(reversepropAbs) 365 { /*lint --e{715}*/ 366 SCIP_INTERVAL childbounds; 367 SCIP_INTERVAL left; 368 SCIP_INTERVAL right; 369 370 assert(scip != NULL); 371 assert(expr != NULL); 372 assert(SCIPexprGetNChildren(expr) == 1); 373 assert(bounds.inf >= 0.0); /* bounds should have been intersected with activity, which is >= 0 */ 374 375 /**! [SnippetExprReversepropAbs] */ 376 /* abs(x) in I -> x \in (-I \cup I) \cap bounds(x) */ 377 right = bounds; /* I */ 378 SCIPintervalSetBounds(&left, -right.sup, -right.inf); /* -I */ 379 380 childbounds = childrenbounds[0]; 381 /* childbounds can be empty here already, but that should work fine here */ 382 383 SCIPintervalIntersect(&left, left, childbounds); /* -I \cap bounds(x), could become empty */ 384 SCIPintervalIntersect(&right, right, childbounds); /* I \cap bounds(x), could become empty */ 385 386 /* compute smallest interval containing (-I \cap bounds(x)) \cup (I \cap bounds(x)) = (-I \cup I) \cap bounds(x) 387 * this works also if left or right is empty 388 */ 389 SCIPintervalUnify(&childrenbounds[0], left, right); 390 /**! [SnippetExprReversepropAbs] */ 391 392 return SCIP_OKAY; 393 } 394 395 /** expression hash callback */ 396 static 397 SCIP_DECL_EXPRHASH(hashAbs) 398 { /*lint --e{715}*/ 399 assert(scip != NULL); 400 assert(expr != NULL); 401 assert(SCIPexprGetNChildren(expr) == 1); 402 assert(hashkey != NULL); 403 assert(childrenhashes != NULL); 404 405 *hashkey = EXPRHDLR_HASHKEY; 406 *hashkey ^= childrenhashes[0]; 407 408 return SCIP_OKAY; 409 } 410 411 /** expression curvature detection callback */ 412 static 413 SCIP_DECL_EXPRCURVATURE(curvatureAbs) 414 { /*lint --e{715}*/ 415 SCIP_EXPR* child; 416 SCIP_INTERVAL childbounds; 417 SCIP_Real childinf; 418 SCIP_Real childsup; 419 420 assert(scip != NULL); 421 assert(expr != NULL); 422 assert(exprcurvature != SCIP_EXPRCURV_UNKNOWN); 423 assert(success != NULL); 424 assert(childcurv != NULL); 425 assert(SCIPexprGetNChildren(expr) == 1); 426 427 child = SCIPexprGetChildren(expr)[0]; 428 assert(child != NULL); 429 430 /**! [SnippetExprCurvatureAbs] */ 431 /* expression is |child|, get domain of child */ 432 SCIP_CALL( SCIPevalExprActivity(scip, child) ); 433 childbounds = SCIPexprGetActivity(child); 434 childinf = SCIPintervalGetInf(childbounds); 435 childsup = SCIPintervalGetSup(childbounds); 436 437 *success = TRUE; 438 if( childinf >= 0.0 ) /* |f(x)| = f(x) */ 439 childcurv[0] = exprcurvature; 440 else if( childsup <= 0.0 ) /* |f(x)| = -f(x) */ 441 childcurv[0] = SCIPexprcurvNegate(exprcurvature); 442 else if( exprcurvature == SCIP_EXPRCURV_CONVEX ) /* |f(x)|, f mixed sign, is convex if f is linear */ 443 childcurv[0] = SCIP_EXPRCURV_LINEAR; 444 else /* |f(x)|, f mixed sign, is never concave nor linear */ 445 *success = FALSE; 446 /**! [SnippetExprCurvatureAbs] */ 447 448 return SCIP_OKAY; 449 } 450 451 /** expression monotonicity detection callback */ 452 static 453 SCIP_DECL_EXPRMONOTONICITY(monotonicityAbs) 454 { /*lint --e{715}*/ 455 SCIP_EXPR* child; 456 SCIP_INTERVAL childbounds; 457 458 assert(scip != NULL); 459 assert(expr != NULL); 460 assert(result != NULL); 461 assert(childidx == 0); 462 463 child = SCIPexprGetChildren(expr)[0]; 464 assert(child != NULL); 465 466 /**! [SnippetExprMonotonicityAbs] */ 467 SCIP_CALL( SCIPevalExprActivity(scip, child) ); 468 childbounds = SCIPexprGetActivity(child); 469 470 if( childbounds.sup <= 0.0 ) 471 *result = SCIP_MONOTONE_DEC; 472 else if( childbounds.inf >= 0.0 ) 473 *result = SCIP_MONOTONE_INC; 474 else 475 *result = SCIP_MONOTONE_UNKNOWN; 476 /**! [SnippetExprMonotonicityAbs] */ 477 478 return SCIP_OKAY; 479 } 480 481 /** expression integrality detection callback */ 482 static 483 SCIP_DECL_EXPRINTEGRALITY(integralityAbs) 484 { /*lint --e{715}*/ 485 SCIP_EXPR* child; 486 487 assert(scip != NULL); 488 assert(expr != NULL); 489 assert(isintegral != NULL); 490 assert(SCIPexprGetNChildren(expr) == 1); 491 492 child = SCIPexprGetChildren(expr)[0]; 493 assert(child != NULL); 494 495 *isintegral = SCIPexprIsIntegral(child); 496 497 return SCIP_OKAY; 498 } 499 500 501 /** creates the handler for absolute expression and includes it into SCIP */ 502 SCIP_RETCODE SCIPincludeExprhdlrAbs( 503 SCIP* scip /**< SCIP data structure */ 504 ) 505 { 506 SCIP_EXPRHDLR* exprhdlr; 507 508 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, 509 EXPRHDLR_PRECEDENCE, evalAbs, NULL) ); 510 assert(exprhdlr != NULL); 511 512 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrAbs, NULL); 513 SCIPexprhdlrSetSimplify(exprhdlr, simplifyAbs); 514 SCIPexprhdlrSetParse(exprhdlr, parseAbs); 515 SCIPexprhdlrSetIntEval(exprhdlr, intevalAbs); 516 SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesAbs, estimateAbs); 517 SCIPexprhdlrSetHash(exprhdlr, hashAbs); 518 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropAbs); 519 SCIPexprhdlrSetDiff(exprhdlr, bwdiffAbs, NULL, NULL); 520 SCIPexprhdlrSetCurvature(exprhdlr, curvatureAbs); 521 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityAbs); 522 SCIPexprhdlrSetIntegrality(exprhdlr, integralityAbs); 523 524 return SCIP_OKAY; 525 } 526 527 /** creates an absolute value expression */ 528 SCIP_RETCODE SCIPcreateExprAbs( 529 SCIP* scip, /**< SCIP data structure */ 530 SCIP_EXPR** expr, /**< pointer where to store expression */ 531 SCIP_EXPR* child, /**< single child */ 532 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 533 void* ownercreatedata /**< data to pass to ownercreate */ 534 ) 535 { 536 assert(expr != NULL); 537 assert(child != NULL); 538 assert(SCIPfindExprhdlr(scip, EXPRHDLR_NAME) != NULL); 539 540 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, EXPRHDLR_NAME), NULL, 1, &child, ownercreate, ownercreatedata) ); 541 542 return SCIP_OKAY; 543 } 544 545 /** indicates whether expression is of abs-type */ /*lint -e{715}*/ 546 SCIP_Bool SCIPisExprAbs( 547 SCIP* scip, /**< SCIP data structure */ 548 SCIP_EXPR* expr /**< expression */ 549 ) 550 { /*lint --e{715}*/ 551 assert(expr != NULL); 552 553 return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0; 554 } 555