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 nlhdlr_convex.c 26 * @ingroup DEFPLUGINS_NLHDLR 27 * @brief nonlinear handlers for convex and concave expressions 28 * @author Benjamin Mueller 29 * @author Stefan Vigerske 30 */ 31 32 #include <string.h> 33 34 #include "scip/nlhdlr_convex.h" 35 #include "scip/pub_nlhdlr.h" 36 #include "scip/scip_expr.h" 37 #include "scip/cons_nonlinear.h" 38 #include "scip/expr_var.h" 39 #include "scip/expr_abs.h" 40 #include "scip/pub_misc_rowprep.h" 41 #include "scip/dbldblarith.h" 42 43 /* fundamental nonlinear handler properties */ 44 #define CONVEX_NLHDLR_NAME "convex" 45 #define CONVEX_NLHDLR_DESC "handler that identifies and estimates convex expressions" 46 #define CONVEX_NLHDLR_DETECTPRIORITY 50 47 #define CONVEX_NLHDLR_ENFOPRIORITY 50 48 49 #define CONCAVE_NLHDLR_NAME "concave" 50 #define CONCAVE_NLHDLR_DESC "handler that identifies and estimates concave expressions" 51 #define CONCAVE_NLHDLR_DETECTPRIORITY 40 52 #define CONCAVE_NLHDLR_ENFOPRIORITY 40 53 54 #define DEFAULT_DETECTSUM FALSE 55 #define DEFAULT_EXTENDEDFORM TRUE 56 #define DEFAULT_CVXQUADRATIC_CONVEX TRUE 57 #define DEFAULT_CVXQUADRATIC_CONCAVE FALSE 58 #define DEFAULT_CVXSIGNOMIAL TRUE 59 #define DEFAULT_CVXPRODCOMP TRUE 60 #define DEFAULT_HANDLETRIVIAL FALSE 61 #define DEFAULT_MAXPERTURB 0.01 62 63 #define INITLPMAXVARVAL 1000.0 /**< maximal absolute value of variable for still generating a linearization cut at that point in initlp */ 64 #define RANDNUMINITSEED 220802 /**< initial seed for random number generator for point perturbation */ 65 66 /*lint -e440*/ 67 /*lint -e441*/ 68 /*lint -e666*/ 69 /*lint -e777*/ 70 71 /* 72 * Data structures 73 */ 74 75 /** nonlinear handler expression data */ 76 struct SCIP_NlhdlrExprData 77 { 78 SCIP_EXPR* nlexpr; /**< expression (copy) for which this nlhdlr estimates */ 79 SCIP_HASHMAP* nlexpr2origexpr; /**< mapping of our copied expression to original expression */ 80 81 int nleafs; /**< number of distinct leafs of nlexpr, i.e., number of distinct (auxiliary) variables handled */ 82 SCIP_EXPR** leafexprs; /**< distinct leaf expressions (excluding value-expressions), thus variables */ 83 }; 84 85 /** nonlinear handler data */ 86 struct SCIP_NlhdlrData 87 { 88 SCIP_Bool isnlhdlrconvex; /**< whether this data is used for the convex nlhdlr (TRUE) or the concave one (FALSE) */ 89 SCIP_SOL* evalsol; /**< solution used for evaluating expression in a different point, 90 e.g., for facet computation of vertex-polyhedral function */ 91 SCIP_RANDNUMGEN* randnumgen; /**< random number generator used to perturb reference point in estimateGradient() */ 92 93 /* parameters */ 94 SCIP_Bool detectsum; /**< whether to run detection when the root of an expression is a non-quadratic sum */ 95 SCIP_Bool extendedform; /**< whether to create extended formulations instead of looking for maximal possible subexpression */ 96 SCIP_Real maxperturb; /**< maximal perturbation of non-differentiable reference point */ 97 98 /* advanced parameters (maybe remove some day) */ 99 SCIP_Bool cvxquadratic; /**< whether to use convexity check on quadratics */ 100 SCIP_Bool cvxsignomial; /**< whether to use convexity check on signomials */ 101 SCIP_Bool cvxprodcomp; /**< whether to use convexity check on product composition f(h)*h */ 102 SCIP_Bool handletrivial; /**< whether to handle trivial expressions, i.e., those where all children are variables */ 103 }; 104 105 /** data struct to be be passed on to vertexpoly-evalfunction (see SCIPcomputeFacetVertexPolyhedralNonlinear) */ 106 typedef struct 107 { 108 SCIP_NLHDLREXPRDATA* nlhdlrexprdata; 109 SCIP_SOL* evalsol; 110 SCIP* scip; 111 } VERTEXPOLYFUN_EVALDATA; 112 113 /** stack used in constructExpr to store expressions that need to be investigated ("to do list") */ 114 typedef struct 115 { 116 SCIP_EXPR** stack; /**< stack elements */ 117 int stacksize; /**< allocated space (in number of pointers) */ 118 int stackpos; /**< position of top element of stack */ 119 } EXPRSTACK; 120 121 #define DECL_CURVCHECK(x) SCIP_RETCODE x( \ 122 SCIP* scip, /**< SCIP data structure */ \ 123 SCIP_EXPR* nlexpr, /**< nlhdlr-expr to check */ \ 124 SCIP_Bool isrootexpr, /**< whether nlexpr is the root from where detection has been started */ \ 125 EXPRSTACK* stack, /**< stack where to add generated leafs */ \ 126 SCIP_HASHMAP* nlexpr2origexpr, /**< mapping from our expression copy to original expression */ \ 127 SCIP_NLHDLRDATA* nlhdlrdata, /**< data of nlhdlr */ \ 128 SCIP_HASHMAP* assumevarfixed, /**< hashmap containing variables that should be assumed to be fixed, or NULL */ \ 129 SCIP_Bool* success /**< whether we found something */ \ 130 ) 131 132 /* 133 * static methods 134 */ 135 136 /** create nlhdlr-expression 137 * 138 * does not create children, i.e., assumes that this will be a leaf 139 */ 140 static 141 SCIP_RETCODE nlhdlrExprCreate( 142 SCIP* scip, /**< SCIP data structure */ 143 SCIP_HASHMAP* nlexpr2origexpr, /**< mapping from copied to original expression */ 144 SCIP_EXPR** nlhdlrexpr, /**< buffer to store created expr */ 145 SCIP_EXPR* origexpr, /**< original expression to be copied */ 146 SCIP_EXPRCURV curv /**< curvature to achieve */ 147 ) 148 { 149 assert(scip != NULL); 150 assert(nlexpr2origexpr != NULL); 151 assert(nlhdlrexpr != NULL); 152 assert(origexpr != NULL); 153 154 if( SCIPexprGetNChildren(origexpr) == 0 ) 155 { 156 /* for leaves, do not copy */ 157 *nlhdlrexpr = origexpr; 158 SCIPcaptureExpr(*nlhdlrexpr); 159 if( !SCIPhashmapExists(nlexpr2origexpr, (void*)*nlhdlrexpr) ) 160 { 161 SCIP_CALL( SCIPhashmapInsert(nlexpr2origexpr, (void*)*nlhdlrexpr, (void*)origexpr) ); 162 } 163 return SCIP_OKAY; 164 } 165 166 /* create copy of expression, but without children */ 167 SCIP_CALL( SCIPduplicateExprShallow(scip, origexpr, nlhdlrexpr, NULL, NULL) ); 168 assert(*nlhdlrexpr != NULL); /* copies within the same SCIP must always work */ 169 170 /* store the curvature we want to get in the curvature flag of the copied expression 171 * it's a bit of a misuse, but once we are done with everything, this is actually correct 172 */ 173 SCIPexprSetCurvature(*nlhdlrexpr, curv); 174 175 /* remember which the original expression was */ 176 SCIP_CALL( SCIPhashmapInsert(nlexpr2origexpr, (void*)*nlhdlrexpr, (void*)origexpr) ); 177 178 return SCIP_OKAY; 179 } 180 181 /** expand nlhdlr-expression by adding children according to original expression */ 182 static 183 SCIP_RETCODE nlhdlrExprGrowChildren( 184 SCIP* scip, /**< SCIP data structure */ 185 SCIP_HASHMAP* nlexpr2origexpr, /**< mapping from copied to original expression */ 186 SCIP_EXPR* nlhdlrexpr, /**< expression for which to create children */ 187 SCIP_EXPRCURV* childrencurv /**< curvature required for children, or NULL if to set to UNKNOWN */ 188 ) 189 { 190 SCIP_EXPR* origexpr; 191 SCIP_EXPR* child; 192 int nchildren; 193 int i; 194 195 assert(scip != NULL); 196 assert(nlhdlrexpr != NULL); 197 assert(SCIPexprGetNChildren(nlhdlrexpr) == 0); 198 199 origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlhdlrexpr); 200 201 nchildren = SCIPexprGetNChildren(origexpr); 202 if( nchildren == 0 ) 203 return SCIP_OKAY; 204 205 for( i = 0; i < nchildren; ++i ) 206 { 207 SCIP_CALL( nlhdlrExprCreate(scip, nlexpr2origexpr, &child, SCIPexprGetChildren(origexpr)[i], 208 childrencurv != NULL ? childrencurv[i] : SCIP_EXPRCURV_UNKNOWN) ); 209 SCIP_CALL( SCIPappendExprChild(scip, nlhdlrexpr, child) ); 210 /* append captures child, so we can release the capture from nlhdlrExprCreate */ 211 SCIP_CALL( SCIPreleaseExpr(scip, &child) ); 212 } 213 214 assert(SCIPexprGetNChildren(nlhdlrexpr) == SCIPexprGetNChildren(origexpr)); 215 216 return SCIP_OKAY; 217 } 218 219 /** evaluate expression at solution w.r.t. auxiliary variables */ 220 static 221 SCIP_DECL_VERTEXPOLYFUN(nlhdlrExprEvalConcave) 222 { 223 VERTEXPOLYFUN_EVALDATA* evaldata = (VERTEXPOLYFUN_EVALDATA*)funcdata; 224 int i; 225 226 assert(args != NULL); 227 assert(nargs == evaldata->nlhdlrexprdata->nleafs); 228 assert(evaldata != NULL); 229 230 #ifdef SCIP_MORE_DEBUG 231 SCIPdebugMsg(evaldata->scip, "eval vertexpolyfun at\n"); 232 #endif 233 for( i = 0; i < nargs; ++i ) 234 { 235 #ifdef SCIP_MORE_DEBUG 236 SCIPdebugMsg(evaldata->scip, " <%s> = %g\n", 237 SCIPvarGetName(SCIPgetVarExprVar(evaldata->nlhdlrexprdata->leafexprs[i])), args[i]); 238 #endif 239 SCIP_CALL_ABORT( SCIPsetSolVal(evaldata->scip, evaldata->evalsol, 240 SCIPgetVarExprVar(evaldata->nlhdlrexprdata->leafexprs[i]), args[i]) ); 241 } 242 243 SCIP_CALL_ABORT( SCIPevalExpr(evaldata->scip, evaldata->nlhdlrexprdata->nlexpr, evaldata->evalsol, 0L) ); 244 245 return SCIPexprGetEvalValue(evaldata->nlhdlrexprdata->nlexpr); 246 } 247 248 /** initialize expression stack */ 249 static 250 SCIP_RETCODE exprstackInit( 251 SCIP* scip, /**< SCIP data structure */ 252 EXPRSTACK* exprstack, /**< stack to initialize */ 253 int initsize /**< initial size */ 254 ) 255 { 256 assert(scip != NULL); 257 assert(exprstack != NULL); 258 assert(initsize > 0); 259 260 SCIP_CALL( SCIPallocBufferArray(scip, &exprstack->stack, initsize) ); 261 exprstack->stacksize = initsize; 262 exprstack->stackpos = -1; 263 264 return SCIP_OKAY; 265 } 266 267 /** free expression stack */ 268 static 269 void exprstackFree( 270 SCIP* scip, /**< SCIP data structure */ 271 EXPRSTACK* exprstack /**< free expression stack */ 272 ) 273 { 274 assert(scip != NULL); 275 assert(exprstack != NULL); 276 277 SCIPfreeBufferArray(scip, &exprstack->stack); 278 } 279 280 /** add expressions to expression stack */ 281 static 282 SCIP_RETCODE exprstackPush( 283 SCIP* scip, /**< SCIP data structure */ 284 EXPRSTACK* exprstack, /**< expression stack */ 285 int nexprs, /**< number of expressions to push */ 286 SCIP_EXPR** exprs /**< expressions to push */ 287 ) 288 { 289 assert(scip != NULL); 290 assert(exprstack != NULL); 291 292 if( nexprs == 0 ) 293 return SCIP_OKAY; 294 295 assert(exprs != NULL); 296 297 if( exprstack->stackpos+1 + nexprs > exprstack->stacksize ) /*lint !e644*/ 298 { 299 exprstack->stacksize = SCIPcalcMemGrowSize(scip, exprstack->stackpos+1 + nexprs); /*lint !e644*/ 300 SCIP_CALL( SCIPreallocBufferArray(scip, &exprstack->stack, exprstack->stacksize) ); 301 } 302 303 memcpy(exprstack->stack + (exprstack->stackpos+1), exprs, nexprs * sizeof(SCIP_EXPR*)); /*lint !e679*/ /*lint !e737*/ 304 exprstack->stackpos += nexprs; 305 306 return SCIP_OKAY; 307 } 308 309 /** gives expression from top of expression stack and removes it from stack */ 310 static 311 SCIP_EXPR* exprstackPop( 312 EXPRSTACK* exprstack /**< expression stack */ 313 ) 314 { 315 assert(exprstack != NULL); 316 assert(exprstack->stackpos >= 0); 317 318 return exprstack->stack[exprstack->stackpos--]; 319 } 320 321 /** indicate whether expression stack is empty */ 322 static 323 SCIP_Bool exprstackIsEmpty( 324 EXPRSTACK* exprstack /**< expression stack */ 325 ) 326 { 327 assert(exprstack != NULL); 328 329 return exprstack->stackpos < 0; 330 } 331 332 /** looks whether given expression is (proper) quadratic and has a given curvature 333 * 334 * If having a given curvature, currently require all arguments of quadratic to be linear. 335 * Hence, not using this for a simple square term, as curvCheckExprhdlr may provide a better condition on argument curvature then. 336 * Also we wouldn't do anything useful for a single bilinear term. 337 * Thus, run on sum's only. 338 */ 339 static 340 DECL_CURVCHECK(curvCheckQuadratic) 341 { /*lint --e{715}*/ 342 SCIP_EXPR* expr; 343 SCIP_EXPRCURV presentcurv; 344 SCIP_EXPRCURV wantedcurv; 345 SCIP_HASHSET* lonelysquares = NULL; 346 SCIP_Bool isquadratic; 347 int nbilinexprs; 348 int nquadexprs; 349 int i; 350 351 assert(nlexpr != NULL); 352 assert(stack != NULL); 353 assert(nlexpr2origexpr != NULL); 354 assert(success != NULL); 355 356 *success = FALSE; 357 358 if( !nlhdlrdata->cvxquadratic ) 359 return SCIP_OKAY; 360 361 if( !SCIPisExprSum(scip, nlexpr) ) 362 return SCIP_OKAY; 363 364 wantedcurv = SCIPexprGetCurvature(nlexpr); 365 if( wantedcurv == SCIP_EXPRCURV_LINEAR ) 366 return SCIP_OKAY; 367 assert(wantedcurv == SCIP_EXPRCURV_CONVEX || wantedcurv == SCIP_EXPRCURV_CONCAVE); 368 369 expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr); 370 assert(expr != NULL); 371 372 /* check whether quadratic */ 373 SCIP_CALL( SCIPcheckExprQuadratic(scip, expr, &isquadratic) ); 374 375 /* if not quadratic, then give up here */ 376 if( !isquadratic ) 377 return SCIP_OKAY; 378 379 SCIPexprGetQuadraticData(expr, NULL, NULL, NULL, NULL, &nquadexprs, &nbilinexprs, NULL, NULL); 380 381 /* if only single square term (+linear), then give up here (let curvCheckExprhdlr handle this) */ 382 if( nquadexprs <= 1 ) 383 return SCIP_OKAY; 384 385 /* if root expression is only sum of squares (+linear) and detectsum is disabled, then give up here, too */ 386 if( isrootexpr && !nlhdlrdata->detectsum && nbilinexprs == 0 ) 387 return SCIP_OKAY; 388 389 /* get curvature of quadratic 390 * TODO as we know what curvature we want, we could first do some simple checks like computing xQx for a random x 391 */ 392 SCIP_CALL( SCIPcomputeExprQuadraticCurvature(scip, expr, &presentcurv, assumevarfixed, FALSE) ); 393 394 /* if not having desired curvature, return */ 395 if( presentcurv != wantedcurv ) 396 return SCIP_OKAY; 397 398 *success = TRUE; 399 400 if( !nlhdlrdata->detectsum ) 401 { 402 /* first step towards block-decomposition of quadratic term: 403 * collect all square-expressions (in original expr) which have no adjacent bilinear term 404 * we will treat these x^2 as linear, i.e., add an auxvar for them, so x^2 maybe linearized 405 * more efficiently (in particular if x is discrete) 406 */ 407 SCIP_CALL( SCIPhashsetCreate(&lonelysquares, SCIPblkmem(scip), nquadexprs) ); 408 for( i = 0; i < nquadexprs; ++i ) 409 { 410 int nadjbilin; 411 SCIP_EXPR* sqrexpr; 412 413 SCIPexprGetQuadraticQuadTerm(expr, i, NULL, NULL, NULL, &nadjbilin, NULL, &sqrexpr); 414 if( nadjbilin == 0 ) 415 { 416 assert(sqrexpr != NULL); 417 SCIP_CALL( SCIPhashsetInsert(lonelysquares, SCIPblkmem(scip), (void*)sqrexpr) ); 418 } 419 } 420 } 421 422 /* add immediate children to nlexpr */ 423 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) ); 424 assert(SCIPexprGetNChildren(nlexpr) == SCIPexprGetNChildren(expr)); 425 426 /* put children that are not square or product on stack 427 * grow child for children that are square or product and put this child on stack 428 * require all children to be linear 429 */ 430 for( i = 0; i < SCIPexprGetNChildren(nlexpr); ++i ) 431 { 432 SCIP_EXPR* child; 433 SCIP_EXPRCURV curvlinear[2] = { SCIP_EXPRCURV_LINEAR, SCIP_EXPRCURV_LINEAR }; 434 435 child = SCIPexprGetChildren(nlexpr)[i]; 436 assert(child != NULL); 437 438 assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)child) == SCIPexprGetChildren(expr)[i]); 439 440 if( SCIPisExprPower(scip, child) && SCIPgetExponentExprPow(child) == 2.0 && 441 (lonelysquares == NULL || !SCIPhashsetExists(lonelysquares, SCIPexprGetChildren(expr)[i])) ) 442 { 443 /* square term that isn't lonely, i.e., orig-version of child is a square-expr and nadjbilin>0 */ 444 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, curvlinear) ); 445 assert(SCIPexprGetNChildren(child) == 1); 446 SCIP_CALL( exprstackPush(scip, stack, 1, SCIPexprGetChildren(child)) ); 447 } 448 else if( SCIPisExprProduct(scip, child) && SCIPexprGetNChildren(SCIPexprGetChildren(expr)[i]) == 2 ) 449 /* using original version of child here as NChildren(child)==0 atm */ 450 { 451 /* bilinear term */ 452 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, curvlinear) ); 453 assert(SCIPexprGetNChildren(child) == 2); 454 SCIP_CALL( exprstackPush(scip, stack, 2, SCIPexprGetChildren(child)) ); 455 } 456 else 457 { 458 /* linear term (or term to be considered as linear) or lonely square term 459 * if we want extended formulations, then require linearity, so an auxvar will be introduced if it is nonlinear 460 * if we do not want extended formulations, then the term needs to have curvature "wantedcurv" 461 * thus, if the coef is negative, then the child needs to have the curvature opposite to "wantedcurv" 462 */ 463 if( nlhdlrdata->extendedform ) 464 SCIPexprSetCurvature(child, SCIP_EXPRCURV_LINEAR); 465 else 466 SCIPexprSetCurvature(child, SCIPexprcurvMultiply(SCIPgetCoefsExprSum(nlexpr)[i], wantedcurv)); 467 SCIP_CALL( exprstackPush(scip, stack, 1, &child) ); 468 } 469 } 470 471 if( lonelysquares != NULL ) 472 SCIPhashsetFree(&lonelysquares, SCIPblkmem(scip)); 473 474 return SCIP_OKAY; 475 } 476 477 /** looks whether top of given expression looks like a signomial that can have a given curvature 478 * 479 * e.g., sqrt(x)*sqrt(y) is convex if x,y >= 0 and x and y are convex 480 * 481 * unfortunately, doesn't work for tls, because i) it's originally sqrt(x*y), and ii) it is expanded into some sqrt(z*y+y); 482 * but works for cvxnonsep_nsig 483 */ 484 static 485 DECL_CURVCHECK(curvCheckSignomial) 486 { /*lint --e{715}*/ 487 SCIP_EXPR* expr; 488 SCIP_EXPR* child; 489 SCIP_Real* exponents; 490 SCIP_INTERVAL* bounds; 491 SCIP_EXPRCURV* curv; 492 int nfactors; 493 int i; 494 495 assert(nlexpr != NULL); 496 assert(stack != NULL); 497 assert(nlexpr2origexpr != NULL); 498 assert(success != NULL); 499 500 *success = FALSE; 501 502 if( !nlhdlrdata->cvxsignomial ) 503 return SCIP_OKAY; 504 505 if( !SCIPisExprProduct(scip, nlexpr) ) 506 return SCIP_OKAY; 507 508 expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr); 509 assert(expr != NULL); 510 511 nfactors = SCIPexprGetNChildren(expr); 512 if( nfactors <= 1 ) /* boooring */ 513 return SCIP_OKAY; 514 515 SCIP_CALL( SCIPallocBufferArray(scip, &exponents, nfactors) ); 516 SCIP_CALL( SCIPallocBufferArray(scip, &bounds, nfactors) ); 517 SCIP_CALL( SCIPallocBufferArray(scip, &curv, nfactors) ); 518 519 for( i = 0; i < nfactors; ++i ) 520 { 521 child = SCIPexprGetChildren(expr)[i]; 522 assert(child != NULL); 523 524 if( !SCIPisExprPower(scip, child) ) 525 { 526 exponents[i] = 1.0; 527 SCIP_CALL( SCIPevalExprActivity(scip, child) ); 528 bounds[i] = SCIPexprGetActivity(child); 529 } 530 else 531 { 532 exponents[i] = SCIPgetExponentExprPow(child); 533 SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(child)[0]) ); 534 bounds[i] = SCIPexprGetActivity(SCIPexprGetChildren(child)[0]); 535 } 536 } 537 538 if( !SCIPexprcurvMonomialInv(SCIPexprcurvMultiply(SCIPgetCoefExprProduct(expr), SCIPexprGetCurvature(nlexpr)), 539 nfactors, exponents, bounds, curv) ) 540 goto TERMINATE; 541 542 /* add immediate children to nlexpr 543 * some entries in curv actually apply to arguments of pow's, will correct this next 544 */ 545 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, curv) ); 546 assert(SCIPexprGetNChildren(nlexpr) == nfactors); 547 548 /* put children that are not power on stack 549 * grow child for children that are power and put this child on stack 550 * if extendedform, then require children to be linear 551 * unless they are linear, an auxvar will be introduced for them and thus they will be handled as var here 552 */ 553 for( i = 0; i < nfactors; ++i ) 554 { 555 child = SCIPexprGetChildren(nlexpr)[i]; 556 assert(child != NULL); 557 558 if( SCIPisExprPower(scip, child) ) 559 { 560 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, child, &curv[i]) ); 561 assert(SCIPexprGetNChildren(child) == 1); 562 child = SCIPexprGetChildren(child)[0]; 563 } 564 assert(SCIPexprGetNChildren(child) == 0); 565 566 if( nlhdlrdata->extendedform ) 567 { 568 SCIPexprSetCurvature(child, SCIP_EXPRCURV_LINEAR); 569 #ifdef SCIP_DEBUG 570 SCIPinfoMessage(scip, NULL, "Extendedform: Require linearity for "); 571 SCIPprintExpr(scip, child, NULL); 572 SCIPinfoMessage(scip, NULL, "\n"); 573 #endif 574 } 575 576 SCIP_CALL( exprstackPush(scip, stack, 1, &child) ); 577 } 578 579 *success = TRUE; 580 581 TERMINATE: 582 SCIPfreeBufferArray(scip, &curv); 583 SCIPfreeBufferArray(scip, &bounds); 584 SCIPfreeBufferArray(scip, &exponents); 585 586 return SCIP_OKAY; 587 } 588 589 /** looks for \f$f(c h(x)+d) h(x) \cdot \text{constant}\f$ and tries to conclude conditions on curvature 590 * 591 * Assume \f$h\f$ is univariate: 592 * - First derivative is \f$f'(c h + d) c h' h + f(c h + d) h'\f$. 593 * - Second derivative is \f{align}{&f''(c h + d) c h' c h' h + f'(c h + d) (c h'' h + c h' h') + f'(c h + d) c h' h' + f(c h + d) h'' \\ 594 * =& f''(c h + d) c^2 h'^2 h + f'(c h + d) c h'' h + 2 f'(c h + d) c h'^2 + f(c h + d) h''.\f} 595 * Remove always positive factors leaves \f[f''(c h + d) h,\quad f'(c h + d) c h'' h,\quad f'(c h + d) c,\quad f(c h + d) h''.\f] 596 * For convexity we want all these terms to be nonnegative. For concavity we want all of them to be nonpositive. 597 * Note, that in each term either both \f$f'(c h + d)\f$ and \f$c\f$ occur, or none of them. 598 * - Thus, \f$f(c h(x) + d)h(x)\f$ is convex if \f$cf\f$ is monotonically increasing \f$(c f' \geq 0)\f$ and either 599 * - \f$f\f$ is convex \f$(f'' \geq 0)\f$ and \f$h\f$ is nonnegative \f$(h \geq 0)\f$ and \f$h\f$ is convex \f$(h'' \geq 0)\f$ and [\f$f\f$ is nonnegative \f$(f \geq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$], or 600 * - \f$f\f$ is concave \f$(f'' \leq 0)\f$ and \f$h\f$ is nonpositive \f$(h \leq 0)\f$ and \f$h\f$ is concave \f$(h'' \leq 0)\f$ and [\f$f\f$ is nonpositive \f$(f \leq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$]. 601 * - Further, \f$f(c h(x) + d)h(x)\f$ is concave if \f$cf\f$ is monotonically decreasing \f$(c f' \leq 0)\f$ and either 602 * - f is convex \f$(f'' \geq 0)\f$ and \f$h\f$ is nonpositive \f$(h \leq 0)\f$ and \f$h\f$ is concave \f$(h'' \leq 0)\f$ and [\f$f\f$ is nonnegative \f$(f \geq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$], or 603 * - f is concave \f$(f'' \leq 0)\f$ and \f$h\f$ is nonnegative \f$(h >= 0)\f$ and \f$h\f$ is convex \f$(h'' \geq 0)\f$ and [\f$f\f$ is nonpositive \f$(f \leq 0)\f$ or \f$h\f$ is linear \f$(h''=0)\f$]. 604 * 605 * This should hold also for multivariate and linear \f$h\f$, as things are invariant under linear transformations. 606 * Similar to signomial, I'll assume that this will also hold for other multivariate \f$h\f$ (someone has a formal proof?). 607 */ 608 static 609 DECL_CURVCHECK(curvCheckProductComposite) 610 { /*lint --e{715}*/ 611 SCIP_EXPR* expr; 612 SCIP_EXPR* f; 613 SCIP_EXPR* h = NULL; 614 SCIP_Real c = 0.0; 615 SCIP_EXPR* ch = NULL; /* c * h */ 616 SCIP_Real d; 617 SCIP_INTERVAL fbounds; 618 SCIP_INTERVAL hbounds; 619 SCIP_MONOTONE fmonotonicity; 620 SCIP_EXPRCURV desiredcurv; 621 SCIP_EXPRCURV hcurv; 622 SCIP_EXPRCURV dummy; 623 int fidx; 624 625 assert(nlexpr != NULL); 626 assert(stack != NULL); 627 assert(nlexpr2origexpr != NULL); 628 assert(success != NULL); 629 630 *success = FALSE; 631 632 if( !nlhdlrdata->cvxprodcomp ) 633 return SCIP_OKAY; 634 635 if( !SCIPisExprProduct(scip, nlexpr) ) 636 return SCIP_OKAY; 637 638 expr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr); 639 assert(expr != NULL); 640 641 if( SCIPexprGetNChildren(expr) != 2 ) 642 return SCIP_OKAY; 643 644 /* check whether we have f(c * h(x)) * h(x) or h(x) * f(c * h(x)) */ 645 for( fidx = 0; fidx <= 1; ++fidx ) 646 { 647 f = SCIPexprGetChildren(expr)[fidx]; 648 649 if( SCIPexprGetNChildren(f) != 1 ) 650 continue; 651 652 ch = SCIPexprGetChildren(f)[0]; 653 c = 1.0; 654 h = ch; 655 656 /* check whether ch is of the form c*h(x), then switch h to child ch */ 657 if( SCIPisExprSum(scip, ch) && SCIPexprGetNChildren(ch) == 1 ) 658 { 659 c = SCIPgetCoefsExprSum(ch)[0]; 660 h = SCIPexprGetChildren(ch)[0]; 661 assert(c != 1.0 || SCIPgetConstantExprSum(ch) != 0.0); /* we could handle this, but it should have been simplified away */ 662 } 663 664 #ifndef NLHDLR_CONVEX_UNITTEST 665 /* can assume that duplicate subexpressions have been identified and comparing pointer is sufficient */ 666 if( SCIPexprGetChildren(expr)[1-fidx] == h ) 667 #else 668 /* called from unittest -> duplicate subexpressions were not identified -> compare more expensively */ 669 if( SCIPcompareExpr(scip, SCIPexprGetChildren(expr)[1-fidx], h) == 0 ) 670 #endif 671 break; 672 } 673 if( fidx == 2 ) 674 return SCIP_OKAY; 675 676 /* constant of c*h(x)+d */ 677 d = h != ch ? SCIPgetConstantExprSum(ch) : 0.0; 678 679 #ifdef SCIP_MORE_DEBUG 680 SCIPinfoMessage(scip, NULL, "f(c*h+d)*h with f = %s, c = %g, d = %g, h = ", SCIPexprhdlrGetName(SCIPexprGetHdlr(f)), c, d); 681 SCIPprintExpr(scip, h, NULL); 682 SCIPinfoMessage(scip, NULL, "\n"); 683 #endif 684 685 assert(c != 0.0); 686 687 SCIP_CALL( SCIPevalExprActivity(scip, f) ); 688 SCIP_CALL( SCIPevalExprActivity(scip, h) ); 689 fbounds = SCIPexprGetActivity(f); 690 hbounds = SCIPexprGetActivity(h); 691 692 /* if h has mixed sign, then cannot conclude anything */ 693 if( hbounds.inf < 0.0 && hbounds.sup > 0.0 ) 694 return SCIP_OKAY; 695 696 /* If we have some convex or concave x*abs(c*x+d), then gradients at x=-d/c may be very wrong due to 697 * rounding errors and non-differentiability of abs() at zero (#3411). Therefore, we skip handling 698 * such expression in this nonlinear handler when one of the bounds of c*x+d is very close to zero. 699 * (If zero is in between the bounds of c*x+d, then the composition wouldn't be regarded as convex/concave anyway.) 700 */ 701 if( SCIPisExprAbs(scip, f) && (SCIPisZero(scip, c*hbounds.inf+d) || SCIPisZero(scip, c*hbounds.sup+d)) ) 702 return SCIP_OKAY; 703 704 SCIP_CALL( SCIPcallExprMonotonicity(scip, f, 0, &fmonotonicity) ); 705 706 /* if f is not monotone, then cannot conclude anything */ 707 if( fmonotonicity == SCIP_MONOTONE_UNKNOWN ) 708 return SCIP_OKAY; 709 710 /* curvature we want to achieve (negate if product has negative coef) */ 711 desiredcurv = SCIPexprcurvMultiply(SCIPgetCoefExprProduct(nlexpr), SCIPexprGetCurvature(nlexpr)); 712 713 /* now check the conditions as stated above */ 714 if( desiredcurv == SCIP_EXPRCURV_CONVEX ) 715 { 716 /* f(c h(x)+d)h(x) is convex if c*f is monotonically increasing (c f' >= 0) and either 717 * - f is convex (f'' >= 0) and h is nonnegative (h >= 0) and h is convex (h'' >= 0) and [f is nonnegative (f >= 0) or h is linear (h''=0)], or 718 * - f is concave (f'' <= 0) and h is nonpositive (h <= 0) and h is concave (h'' <= 0) and [f is nonpositive (f <= 0) or h is linear (h''=0)] 719 * as the curvature requirements on f are on f only and not the composition f(h), we can ignore the requirements returned by SCIPcallExprCurvature (last arg) 720 */ 721 if( (c > 0.0 && fmonotonicity != SCIP_MONOTONE_INC) || (c < 0.0 && fmonotonicity != SCIP_MONOTONE_DEC) ) 722 return SCIP_OKAY; 723 724 /* check whether f can be convex (h>=0) or concave (h<=0), resp., and derive requirements for h */ 725 if( hbounds.inf >= 0 ) 726 { 727 SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONVEX, success, &dummy) ); 728 729 /* now h also needs to be convex; and if f < 0, then h actually needs to be linear */ 730 if( fbounds.inf < 0.0 ) 731 hcurv = SCIP_EXPRCURV_LINEAR; 732 else 733 hcurv = SCIP_EXPRCURV_CONVEX; 734 } 735 else 736 { 737 SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONCAVE, success, &dummy) ); 738 739 /* now h also needs to be concave; and if f > 0, then h actually needs to be linear */ 740 if( fbounds.sup > 0.0 ) 741 hcurv = SCIP_EXPRCURV_LINEAR; 742 else 743 hcurv = SCIP_EXPRCURV_CONCAVE; 744 } 745 } 746 else 747 { 748 /* f(c h(x)+d)*h(x) is concave if c*f is monotonically decreasing (c f' <= 0) and either 749 * - f is convex (f'' >= 0) and h is nonpositive (h <= 0) and h is concave (h'' <= 0) and [f is nonnegative (f >= 0) or h is linear (h''=0)], or 750 * - f is concave (f'' <= 0) and h is nonnegative (h >= 0) and h is convex (h'' >= 0) and [f is nonpositive (f <= 0) or h is linear (h''=0)] 751 * as the curvature requirements on f are on f only and not the composition f(h), we can ignore the requirements returned by SCIPcallExprCurvature (last arg) 752 */ 753 if( (c > 0.0 && fmonotonicity != SCIP_MONOTONE_DEC) || (c < 0.0 && fmonotonicity != SCIP_MONOTONE_INC) ) 754 return SCIP_OKAY; 755 756 /* check whether f can be convex (h<=0) or concave (h>=0), resp., and derive requirements for h */ 757 if( hbounds.sup <= 0 ) 758 { 759 SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONVEX, success, &dummy) ); 760 761 /* now h also needs to be concave; and if f < 0, then h actually needs to be linear */ 762 if( fbounds.inf < 0.0 ) 763 hcurv = SCIP_EXPRCURV_LINEAR; 764 else 765 hcurv = SCIP_EXPRCURV_CONCAVE; 766 } 767 else 768 { 769 SCIP_CALL( SCIPcallExprCurvature(scip, f, SCIP_EXPRCURV_CONCAVE, success, &dummy) ); 770 771 /* now h also needs to be convex; and if f > 0, then h actually needs to be linear */ 772 if( fbounds.sup > 0.0 ) 773 hcurv = SCIP_EXPRCURV_LINEAR; 774 else 775 hcurv = SCIP_EXPRCURV_CONVEX; 776 } 777 } 778 779 if( !*success ) 780 return SCIP_OKAY; 781 782 /* add immediate children (f and ch) to nlexpr; we set required curvature for h further below */ 783 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) ); 784 assert(SCIPexprGetNChildren(nlexpr) == 2); 785 786 /* copy of f (and h) should have same child position in nlexpr as f (and h) has on expr (resp) */ 787 assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(nlexpr)[fidx]) == (void*)f); 788 #ifndef NLHDLR_CONVEX_UNITTEST 789 assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(nlexpr)[1-fidx]) == (void*)h); 790 #endif 791 /* push this h onto stack for further checking */ 792 SCIP_CALL( exprstackPush(scip, stack, 1, &(SCIPexprGetChildren(nlexpr)[1-fidx])) ); 793 794 /* if we prefer extended formulations, then we always want h() to be linear */ 795 if( nlhdlrdata->extendedform ) 796 hcurv = SCIP_EXPRCURV_LINEAR; 797 798 /* h-child of product should have curvature hcurv */ 799 SCIPexprSetCurvature(SCIPexprGetChildren(nlexpr)[1-fidx], hcurv); 800 801 if( h != ch ) 802 { 803 /* add copy of ch as child to copy of f */ 804 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, SCIPexprGetChildren(nlexpr)[fidx], NULL) ); 805 assert(SCIPexprGetNChildren(SCIPexprGetChildren(nlexpr)[fidx]) == 1); 806 assert(SCIPhashmapGetImage(nlexpr2origexpr, (void*)SCIPexprGetChildren(SCIPexprGetChildren(nlexpr)[fidx])[0]) == (void*)ch); 807 808 /* add copy of h (created above as child of product) as child in copy of ch */ 809 SCIP_CALL( SCIPappendExprChild(scip, 810 SCIPexprGetChildren(SCIPexprGetChildren(nlexpr)[fidx])[0] /* copy of ch */, 811 SCIPexprGetChildren(nlexpr)[1-fidx] /* copy of h */) ); 812 } 813 else 814 { 815 /* add copy of h (created above as child of product) as child in copy of f */ 816 SCIP_CALL( SCIPappendExprChild(scip, 817 SCIPexprGetChildren(nlexpr)[fidx] /* copy of f */, 818 SCIPexprGetChildren(nlexpr)[1-fidx] /* copy of h */) ); 819 } 820 821 return SCIP_OKAY; 822 } 823 824 /** use expression handlers curvature callback to check whether given curvature can be achieved */ 825 static 826 DECL_CURVCHECK(curvCheckExprhdlr) 827 { /*lint --e{715}*/ 828 SCIP_EXPR* origexpr; 829 int nchildren; 830 SCIP_EXPRCURV* childcurv; 831 832 assert(nlexpr != NULL); 833 assert(stack != NULL); 834 assert(nlexpr2origexpr != NULL); 835 assert(success != NULL); 836 837 origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, nlexpr); 838 assert(origexpr != NULL); 839 nchildren = SCIPexprGetNChildren(origexpr); 840 841 if( nchildren == 0 ) 842 { 843 /* if originally no children, then should be var or value, which should have every curvature, 844 * so should always be success 845 */ 846 SCIP_CALL( SCIPcallExprCurvature(scip, origexpr, SCIPexprGetCurvature(nlexpr), success, NULL) ); 847 assert(*success); 848 849 return SCIP_OKAY; 850 } 851 852 /* ignore sums if > 1 children 853 * NOTE: this means that for something like 1+f(x), even if f is a trivial convex expression, we would handle 1+f(x) 854 * with this nlhdlr, instead of formulating this as 1+z and handling z=f(x) with the default nlhdlr, i.e., the exprhdlr 855 * today, I prefer handling this here, as it avoids introducing an extra auxiliary variable 856 */ 857 if( isrootexpr && !nlhdlrdata->detectsum && SCIPisExprSum(scip, nlexpr) && nchildren > 1 ) 858 return SCIP_OKAY; 859 860 SCIP_CALL( SCIPallocBufferArray(scip, &childcurv, nchildren) ); 861 862 /* check whether and under which conditions origexpr can have desired curvature */ 863 SCIP_CALL( SCIPcallExprCurvature(scip, origexpr, SCIPexprGetCurvature(nlexpr), success, childcurv) ); 864 #ifdef SCIP_MORE_DEBUG 865 SCIPprintExpr(scip, origexpr, NULL); 866 SCIPinfoMessage(scip, NULL, " is %s? %d\n", SCIPexprcurvGetName(SCIPexprGetCurvature(nlexpr)), *success); 867 #endif 868 if( !*success ) 869 goto TERMINATE; 870 871 /* if origexpr can have curvature curv, then don't treat it as leaf, but include its children */ 872 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, childcurv) ); 873 assert(SCIPexprGetChildren(nlexpr) != NULL); 874 assert(SCIPexprGetNChildren(nlexpr) == nchildren); 875 876 /* If we prefer extended formulations, then require all children to be linear. 877 * Unless they are, auxvars will be introduced and they will be handles as variables, which can be an 878 * advantage in the context of extended formulations. 879 */ 880 if( nlhdlrdata->extendedform ) 881 { 882 int i; 883 for( i = 0; i < nchildren; ++i ) 884 SCIPexprSetCurvature(SCIPexprGetChildren(nlexpr)[i], SCIP_EXPRCURV_LINEAR); 885 #ifdef SCIP_DEBUG 886 SCIPinfoMessage(scip, NULL, "require linearity for children of "); 887 SCIPprintExpr(scip, origexpr, NULL); 888 SCIPinfoMessage(scip, NULL, "\n"); 889 #endif 890 } 891 892 /* add children expressions to to-do list (stack) */ 893 SCIP_CALL( exprstackPush(scip, stack, nchildren, SCIPexprGetChildren(nlexpr)) ); 894 895 TERMINATE: 896 SCIPfreeBufferArray(scip, &childcurv); 897 898 return SCIP_OKAY; 899 } 900 901 /** curvature check and expression-growing methods 902 * 903 * some day this could be plugins added by users at runtime, but for now we have a fixed list here 904 * @note curvCheckExprhdlr should be last 905 */ 906 static DECL_CURVCHECK((*CURVCHECKS[])) = { curvCheckProductComposite, curvCheckSignomial, curvCheckQuadratic, curvCheckExprhdlr }; 907 /** number of curvcheck methods */ 908 static const int NCURVCHECKS = sizeof(CURVCHECKS) / sizeof(void*); 909 910 /** checks whether expression is a sum with more than one child and each child being a variable or going to be a variable if `expr` is a nlhdlr-specific copy 911 * 912 * Within constructExpr(), we can have an expression of any type which is a copy of an original expression, 913 * but without children. At the end of constructExpr() (after the loop with the stack), these expressions 914 * will remain as leafs and will eventually be turned into variables in collectLeafs(). Thus, we treat 915 * every child that has no children as if it were a variable. Theoretically, there is still the possibility 916 * that it could be a constant (value-expression), but simplify should have removed these. 917 */ 918 static 919 SCIP_Bool exprIsMultivarLinear( 920 SCIP* scip, /**< SCIP data structure */ 921 SCIP_EXPR* expr /**< expression to check */ 922 ) 923 { 924 int nchildren; 925 int c; 926 927 assert(expr != NULL); 928 929 if( !SCIPisExprSum(scip, expr) ) 930 return FALSE; 931 932 nchildren = SCIPexprGetNChildren(expr); 933 if( nchildren <= 1 ) 934 return FALSE; 935 936 for( c = 0; c < nchildren; ++c ) 937 /*if( !SCIPisExprVar(scip, SCIPexprGetChildren(expr)[c]) ) */ 938 if( SCIPexprGetNChildren(SCIPexprGetChildren(expr)[c]) > 0 ) 939 return FALSE; 940 941 return TRUE; 942 } 943 944 /** constructs a subexpression (as nlhdlr-expression) of maximal size that has a given curvature 945 * 946 * If the curvature cannot be achieved for an expression in the original expression graph, 947 * then this expression becomes a leaf in the nlhdlr-expression. 948 * 949 * Sets `*rootnlexpr` to NULL if failed. 950 */ 951 static 952 SCIP_RETCODE constructExpr( 953 SCIP* scip, /**< SCIP data structure */ 954 SCIP_NLHDLRDATA* nlhdlrdata, /**< nonlinear handler data */ 955 SCIP_EXPR** rootnlexpr, /**< buffer to store created expression */ 956 SCIP_HASHMAP* nlexpr2origexpr, /**< mapping from our expression copy to original expression */ 957 int* nleafs, /**< number of leafs in constructed expression */ 958 SCIP_EXPR* rootexpr, /**< expression */ 959 SCIP_EXPRCURV curv, /**< curvature to achieve */ 960 SCIP_HASHMAP* assumevarfixed, /**< hashmap containing variables that should be assumed to be fixed, or NULL */ 961 SCIP_Bool assumecurvature, /**< whether to assume that desired curvature is given (skips curvature checks) */ 962 SCIP_Bool* curvsuccess /**< pointer to store whether the curvature could be achieved 963 * w.r.t. the original variables (might be NULL) */ 964 ) 965 { 966 SCIP_EXPR* nlexpr; 967 EXPRSTACK stack; /* to do list: expressions where to check whether they can have the desired curvature when taking their children into account */ 968 int oldstackpos; 969 SCIP_Bool isrootexpr = TRUE; 970 971 assert(scip != NULL); 972 assert(nlhdlrdata != NULL); 973 assert(rootnlexpr != NULL); 974 assert(nlexpr2origexpr != NULL); 975 assert(nleafs != NULL); 976 assert(rootexpr != NULL); 977 assert(curv == SCIP_EXPRCURV_CONVEX || curv == SCIP_EXPRCURV_CONCAVE); 978 979 /* create root expression */ 980 SCIP_CALL( nlhdlrExprCreate(scip, nlexpr2origexpr, rootnlexpr, rootexpr, curv) ); 981 982 *nleafs = 0; 983 if( curvsuccess != NULL ) 984 *curvsuccess = TRUE; 985 986 SCIP_CALL( exprstackInit(scip, &stack, 20) ); 987 SCIP_CALL( exprstackPush(scip, &stack, 1, rootnlexpr) ); 988 while( !exprstackIsEmpty(&stack) ) 989 { 990 /* take expression from stack */ 991 nlexpr = exprstackPop(&stack); 992 assert(nlexpr != NULL); 993 assert(SCIPexprGetNChildren(nlexpr) == 0); 994 995 oldstackpos = stack.stackpos; 996 if( nlhdlrdata->isnlhdlrconvex && !SCIPexprhdlrHasBwdiff(SCIPexprGetHdlr(nlexpr)) ) 997 { 998 /* if bwdiff is not implemented, then we could not generate cuts in the convex nlhdlr, so "stop" (treat nlexpr as variable) */ 999 } 1000 else if( !nlhdlrdata->isnlhdlrconvex && exprIsMultivarLinear(scip, (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr)) ) 1001 { 1002 /* if we are in the concave handler, we would like to treat linear multivariate subexpressions by a new auxvar always, 1003 * e.g., handle log(x+y) as log(z), z=x+y, because the estimation problem will be smaller then without making the estimator worse 1004 * (cons_nonlinear does this, too) 1005 * this check takes care of this when x and y are original variables 1006 * however, it isn't unlikely that we will have sums that become linear after we add auxvars for some children 1007 * this will be handled in a postprocessing below 1008 * for now, the check is performed on the original expression since there is not enough information in nlexpr yet 1009 */ 1010 #ifdef SCIP_MORE_DEBUG 1011 SCIPprintExpr(scip, SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr), NULL); 1012 SCIPinfoMessage(scip, NULL, "... is a multivariate linear sum that we'll treat as auxvar\n"); 1013 #endif 1014 } 1015 else if( SCIPexprGetCurvature(nlexpr) != SCIP_EXPRCURV_UNKNOWN && !assumecurvature ) 1016 { 1017 /* if we are here, either convexity or concavity is required; try to check for this curvature */ 1018 SCIP_Bool success; 1019 int method; 1020 1021 /* try through curvature check methods until one succeeds */ 1022 for( method = 0; method < NCURVCHECKS; ++method ) 1023 { 1024 SCIP_CALL( CURVCHECKS[method](scip, nlexpr, isrootexpr, &stack, nlexpr2origexpr, nlhdlrdata, assumevarfixed, &success) ); 1025 if( success ) 1026 break; 1027 } 1028 } 1029 else 1030 { 1031 /* if we don't care about curvature in this subtree anymore (very unlikely), 1032 * or we are told to assume that the desired curvature is present (assumecurvature==TRUE), 1033 * then only continue iterating this subtree to assemble leaf expressions 1034 */ 1035 SCIP_CALL( nlhdlrExprGrowChildren(scip, nlexpr2origexpr, nlexpr, NULL) ); 1036 1037 /* add children expressions, if any, to to-do list (stack) */ 1038 SCIP_CALL( exprstackPush(scip, &stack, SCIPexprGetNChildren(nlexpr), SCIPexprGetChildren(nlexpr)) ); 1039 } 1040 assert(stack.stackpos >= oldstackpos); /* none of the methods above should have removed something from the stack */ 1041 1042 isrootexpr = FALSE; 1043 1044 /* if nothing was added, then none of the successors of nlexpr were added to the stack 1045 * this is either because nlexpr was already a variable or value expressions, thus a leaf, 1046 * or because the desired curvature could not be achieved, so it will be handled as variables, thus a leaf 1047 */ 1048 if( stack.stackpos == oldstackpos ) 1049 { 1050 ++*nleafs; 1051 1052 /* check whether the new leaf is not an original variable (or constant) */ 1053 if( curvsuccess != NULL && !SCIPisExprVar(scip, nlexpr) && !SCIPisExprValue(scip, nlexpr) ) 1054 *curvsuccess = FALSE; 1055 } 1056 } 1057 1058 exprstackFree(scip, &stack); 1059 1060 if( !nlhdlrdata->isnlhdlrconvex && *rootnlexpr != NULL ) 1061 { 1062 /* remove multivariate linear subexpressions, that is, change some f(z1+z2) into f(z3) (z3=z1+z2 will be done by nlhdlr_default) 1063 * this handles the case that was not covered by the above check, which could recognize f(x+y) for x, y original variables 1064 */ 1065 SCIP_EXPRITER* it; 1066 1067 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 1068 SCIP_CALL( SCIPexpriterInit(it, *rootnlexpr, SCIP_EXPRITER_DFS, FALSE) ); 1069 SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD); 1070 1071 while( !SCIPexpriterIsEnd(it) ) 1072 { 1073 SCIP_EXPR* child; 1074 1075 child = SCIPexpriterGetChildExprDFS(it); 1076 assert(child != NULL); 1077 1078 /* We want to change some f(x+y+z) into just f(), where f is the expression the iterator points to 1079 * and x+y+z is child. A child of a child, e.g., z, may not be a variable yet (these are added in collectLeafs later), 1080 * but an expression of some nonlinear type without children. 1081 */ 1082 if( exprIsMultivarLinear(scip, child) ) 1083 { 1084 /* turn child (x+y+z) into a sum without children 1085 * collectLeafs() should then replace this by an auxvar 1086 */ 1087 #ifdef SCIP_MORE_DEBUG 1088 SCIPprintExpr(scip, child, NULL); 1089 SCIPinfoMessage(scip, NULL, "... is a multivariate linear sum that we'll treat as auxvar instead (postprocess)\n"); 1090 #endif 1091 1092 /* TODO remove children from nlexpr2origexpr ? 1093 * should also do this if they are not used somewhere else; we could check nuses for this 1094 * however, it shouldn't matter to have some stray entries in the hashmap either 1095 */ 1096 SCIP_CALL( SCIPremoveExprChildren(scip, child) ); 1097 assert(SCIPexprGetNChildren(child) == 0); 1098 1099 (void) SCIPexpriterSkipDFS(it); 1100 } 1101 else 1102 { 1103 (void) SCIPexpriterGetNext(it); 1104 } 1105 } 1106 1107 SCIPfreeExpriter(&it); 1108 } 1109 1110 if( *rootnlexpr != NULL ) 1111 { 1112 SCIP_Bool istrivial = TRUE; 1113 1114 /* if handletrivial is enabled, then only require that rootnlexpr itself has required curvature (so has children; see below) and 1115 * that we are not a trivial sum (because the previous implementation of this nlhdlr didn't allow this, either) 1116 */ 1117 if( !nlhdlrdata->handletrivial || SCIPisExprSum(scip, *rootnlexpr) ) 1118 { 1119 /* if all children do not have children, i.e., are variables, or will be replaced by auxvars, then free 1120 * also if rootnlexpr has no children, then free 1121 */ 1122 int i; 1123 for( i = 0; i < SCIPexprGetNChildren(*rootnlexpr); ++i ) 1124 { 1125 if( SCIPexprGetNChildren(SCIPexprGetChildren(*rootnlexpr)[i]) > 0 ) 1126 { 1127 istrivial = FALSE; 1128 break; 1129 } 1130 } 1131 } 1132 else if( SCIPexprGetNChildren(*rootnlexpr) > 0 ) /* if handletrivial, then just require children */ 1133 istrivial = FALSE; 1134 1135 if( istrivial ) 1136 { 1137 SCIP_CALL( SCIPreleaseExpr(scip, rootnlexpr) ); 1138 } 1139 } 1140 1141 return SCIP_OKAY; 1142 } 1143 1144 /** collects (non-value) leaf expressions and ensure that they correspond to a variable (original or auxiliary) 1145 * 1146 * For children where we could not achieve the desired curvature, get the auxvar and replace the child by a 1147 * var-expression that points to this auxvar. 1148 * Collect all leaf expressions (if not a value-expression) and index them. 1149 */ 1150 static 1151 SCIP_RETCODE collectLeafs( 1152 SCIP* scip, /**< SCIP data structure */ 1153 SCIP_NLHDLREXPRDATA* nlhdlrexprdata /**< nlhdlr expression data */ 1154 ) 1155 { 1156 SCIP_EXPRITER* it; 1157 SCIP_EXPR* nlexpr; 1158 SCIP_HASHMAP* leaf2index; 1159 int i; 1160 1161 assert(nlhdlrexprdata != NULL); 1162 assert(nlhdlrexprdata->nlexpr != NULL); 1163 assert(nlhdlrexprdata->nlexpr2origexpr != NULL); 1164 /* nleafs should be the upper bound on the number of variables given by constructExpr 1165 * leafexprs should be NULL, as this is what we want to setup here 1166 */ 1167 assert(nlhdlrexprdata->nleafs > 0); 1168 assert(nlhdlrexprdata->leafexprs == NULL); 1169 1170 /* collect all auxvars and collect all variables */ 1171 SCIP_CALL( SCIPhashmapCreate(&leaf2index, SCIPblkmem(scip), nlhdlrexprdata->nleafs) ); 1172 nlhdlrexprdata->nleafs = 0; /* we start a new count, this time skipping value-expressions */ 1173 1174 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 1175 SCIP_CALL( SCIPexpriterInit(it, nlhdlrexprdata->nlexpr, SCIP_EXPRITER_DFS, FALSE) ); 1176 SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD); 1177 1178 for( nlexpr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); nlexpr = SCIPexpriterGetNext(it) ) 1179 { 1180 SCIP_EXPR* child; 1181 SCIP_EXPR* origexpr; 1182 1183 assert(nlexpr != NULL); 1184 1185 child = SCIPexpriterGetChildExprDFS(it); 1186 1187 /* if the to-be-visited child has children, then it doesn't need to be replaced by a new expression (representing the auxvar) */ 1188 if( SCIPexprGetNChildren(child) > 0 ) 1189 continue; 1190 1191 origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)child); 1192 assert(origexpr != NULL); 1193 1194 if( SCIPexprGetNChildren(origexpr) > 0 ) 1195 { 1196 SCIP_EXPR* newchild; 1197 int childidx; 1198 SCIP_VAR* var; 1199 1200 /* having a child that had children in original but not in copy means that we could not achieve the desired curvature 1201 * thus, replace by a new child that points to the auxvar of the original expression 1202 * we registered in createNlhdlrExprData that we need an auxvar, so it should exist now 1203 */ 1204 var = SCIPgetExprAuxVarNonlinear(origexpr); 1205 assert(var != NULL); 1206 1207 SCIP_CALL( SCIPcreateExprVar(scip, &newchild, var, NULL, NULL) ); /* this captures newchild once */ 1208 1209 childidx = SCIPexpriterGetChildIdxDFS(it); 1210 SCIP_CALL( SCIPreplaceExprChild(scip, nlexpr, childidx, newchild) ); /* this captures newchild again */ 1211 1212 /* do not remove child->origexpr from hashmap, as child may appear again due to common subexprs 1213 * (created by curvCheckProductComposite, for example) 1214 * if it doesn't reappear, though, but the memory address is reused, we need to make sure it 1215 * points to the right origexpr 1216 */ 1217 /* SCIP_CALL( SCIPhashmapRemove(nlexpr2origexpr, (void*)child) ); */ 1218 SCIP_CALL( SCIPhashmapSetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)newchild, (void*)origexpr) ); 1219 1220 if( !SCIPhashmapExists(leaf2index, (void*)newchild) ) 1221 { 1222 /* new leaf -> new index and remember in hashmap */ 1223 SCIP_CALL( SCIPhashmapInsertInt(leaf2index, (void*)newchild, nlhdlrexprdata->nleafs++) ); 1224 } 1225 1226 child = newchild; 1227 SCIP_CALL( SCIPreleaseExpr(scip, &newchild) ); /* because it was captured by both create and replace */ 1228 } 1229 else if( SCIPisExprVar(scip, child) ) 1230 { 1231 /* if variable, then add to hashmap, if not already there */ 1232 if( !SCIPhashmapExists(leaf2index, (void*)child) ) 1233 { 1234 SCIP_CALL( SCIPhashmapInsertInt(leaf2index, (void*)child, nlhdlrexprdata->nleafs++) ); 1235 } 1236 } 1237 /* else: it's probably a value-expression, nothing to do */ 1238 1239 /* update integrality flag for future leaf expressions: convex nlhdlr may use this information */ 1240 SCIP_CALL( SCIPcomputeExprIntegrality(scip, child) ); 1241 } 1242 assert(nlhdlrexprdata->nleafs > 0); 1243 1244 SCIPfreeExpriter(&it); 1245 1246 /* assemble auxvars array */ 1247 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(nlhdlrexprdata->leafexprs), nlhdlrexprdata->nleafs) ); 1248 for( i = 0; i < SCIPhashmapGetNEntries(leaf2index); ++i ) 1249 { 1250 SCIP_HASHMAPENTRY* entry; 1251 SCIP_EXPR* leaf; 1252 int idx; 1253 1254 entry = SCIPhashmapGetEntry(leaf2index, i); 1255 if( entry == NULL ) 1256 continue; 1257 1258 leaf = (SCIP_EXPR*) SCIPhashmapEntryGetOrigin(entry); 1259 assert(leaf != NULL); 1260 assert(SCIPisExprVar(scip, leaf)); 1261 1262 idx = SCIPhashmapEntryGetImageInt(entry); 1263 assert(idx >= 0); 1264 assert(idx < nlhdlrexprdata->nleafs); 1265 1266 nlhdlrexprdata->leafexprs[idx] = leaf; 1267 1268 SCIPdebugMsg(scip, "leaf %d: <%s>\n", idx, SCIPvarGetName(SCIPgetVarExprVar(leaf))); 1269 } 1270 1271 SCIPhashmapFree(&leaf2index); 1272 1273 return SCIP_OKAY; 1274 } 1275 1276 /** creates nonlinear handler expression data structure and registers expr usage */ 1277 static 1278 SCIP_RETCODE createNlhdlrExprData( 1279 SCIP* scip, /**< SCIP data structure */ 1280 SCIP_NLHDLRDATA* nlhdlrdata, /**< nlhdlr data */ 1281 SCIP_NLHDLREXPRDATA** nlhdlrexprdata, /**< pointer to store nlhdlr expression data */ 1282 SCIP_EXPR* expr, /**< original expression */ 1283 SCIP_EXPR* nlexpr, /**< our copy of expression */ 1284 SCIP_HASHMAP* nlexpr2origexpr, /**< mapping of expression copy to original */ 1285 int nleafs, /**< number of leafs as counted by constructExpr */ 1286 SCIP_NLHDLR_METHOD participating /**< the enfo methods in which we plan to participate */ 1287 ) 1288 { 1289 SCIP_EXPRITER* it; 1290 SCIP_Bool usingaux; 1291 1292 assert(scip != NULL); 1293 assert(expr != NULL); 1294 assert(nlhdlrexprdata != NULL); 1295 assert(*nlhdlrexprdata == NULL); 1296 assert(nlexpr != NULL); 1297 assert(nlexpr2origexpr != NULL); 1298 1299 assert(SCIPexprGetNChildren(nlexpr) > 0); 1300 assert(SCIPexprGetChildren(nlexpr) != NULL); 1301 1302 SCIP_CALL( SCIPallocClearBlockMemory(scip, nlhdlrexprdata) ); 1303 (*nlhdlrexprdata)->nlexpr = nlexpr; 1304 (*nlhdlrexprdata)->nlexpr2origexpr = nlexpr2origexpr; 1305 (*nlhdlrexprdata)->nleafs = nleafs; 1306 1307 usingaux = FALSE; 1308 1309 SCIP_CALL( SCIPcreateExpriter(scip, &it) ); 1310 SCIP_CALL( SCIPexpriterInit(it, nlexpr, SCIP_EXPRITER_DFS, FALSE) ); 1311 SCIPexpriterSetStagesDFS(it, SCIP_EXPRITER_VISITINGCHILD); 1312 1313 for( ; !SCIPexpriterIsEnd(it); (void) SCIPexpriterGetNext(it) ) 1314 { 1315 SCIP_EXPR* child; 1316 SCIP_EXPR* origexpr; 1317 1318 /* check whether to-be-visited child needs to be replaced by a new expression (representing the auxvar) 1319 * if child has children, then that is not the case 1320 * if child has no children, but also corresponding origexpr has no chilren, then this is also not the case 1321 */ 1322 child = SCIPexpriterGetChildExprDFS(it); 1323 if( SCIPexprGetNChildren(child) > 0 ) 1324 continue; 1325 1326 origexpr = (SCIP_EXPR*)SCIPhashmapGetImage(nlexpr2origexpr, (void*)child); 1327 assert(origexpr != NULL); 1328 1329 /* if child had children in original but not in copy means that we could not achieve the desired curvature 1330 * thus, we will later replace by a new child that points to the auxvar of the original expression 1331 * as we do not have the auxvar now, we will only register that we will need the auxvar later (if origexpr isn't a variable or constant) 1332 * if we are working for the concave nlhdlr, then we also indicate interest on the exprs activity for estimate (distinguish below or above) 1333 */ 1334 SCIP_CALL( SCIPregisterExprUsageNonlinear(scip, origexpr, 1335 SCIPexprGetNChildren(origexpr) > 0, FALSE, 1336 !nlhdlrdata->isnlhdlrconvex && (participating & SCIP_NLHDLR_METHOD_SEPABELOW), 1337 !nlhdlrdata->isnlhdlrconvex && (participating & SCIP_NLHDLR_METHOD_SEPAABOVE)) ); 1338 1339 /* remember that we use an auxvar */ 1340 if( SCIPexprGetNChildren(origexpr) > 0 ) 1341 usingaux = TRUE; 1342 } 1343 1344 SCIPfreeExpriter(&it); 1345 1346 #ifdef SCIP_DEBUG 1347 SCIPprintExpr(scip, nlexpr, NULL); 1348 SCIPinfoMessage(scip, NULL, " (%p) is handled as %s\n", SCIPhashmapGetImage(nlexpr2origexpr, (void*)nlexpr), 1349 SCIPexprcurvGetName(SCIPexprGetCurvature(nlexpr))); 1350 #endif 1351 1352 /* If we don't work on the extended formulation, then set curvature also in original expression 1353 * (in case someone wants to pick this up; this might be removed again). 1354 * This doesn't ensure that every convex or concave original expression is actually marked here. 1355 * Not only because our tests are incomprehensive, but also because we may not detect on sums, 1356 * prefer extended formulations (in nlhdlr_convex), or introduce auxvars for linear subexpressions 1357 * on purpose (in nlhdlr_concave). 1358 */ 1359 if( !usingaux ) 1360 SCIPexprSetCurvature(expr, SCIPexprGetCurvature(nlexpr)); 1361 1362 return SCIP_OKAY; 1363 } 1364 1365 /** adds an estimator for a vertex-polyhedral (e.g., concave) function to a given rowprep 1366 * 1367 * Calls \ref SCIPcomputeFacetVertexPolyhedralNonlinear() for given function and 1368 * box set to local bounds of auxiliary variables. 1369 */ 1370 static 1371 SCIP_RETCODE estimateVertexPolyhedral( 1372 SCIP* scip, /**< SCIP data structure */ 1373 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */ 1374 SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */ 1375 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 1376 SCIP_SOL* sol, /**< solution to use, unless usemidpoint is TRUE */ 1377 SCIP_Bool usemidpoint, /**< whether to use the midpoint of the domain instead of sol */ 1378 SCIP_Bool overestimate, /**< whether over- or underestimating */ 1379 SCIP_Real targetvalue, /**< a target value to achieve; if not reachable, then can give up early */ 1380 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 1381 SCIP_Bool* success /**< buffer to store whether successful */ 1382 ) 1383 { 1384 SCIP_NLHDLRDATA* nlhdlrdata; 1385 VERTEXPOLYFUN_EVALDATA evaldata; 1386 SCIP_Real* xstar; 1387 SCIP_Real* box; 1388 SCIP_Real facetconstant; 1389 SCIP_VAR* var; 1390 int i; 1391 SCIP_Bool allfixed; 1392 1393 assert(scip != NULL); 1394 assert(nlhdlr != NULL); 1395 assert(nlhdlrexprdata != NULL); 1396 assert(rowprep != NULL); 1397 assert(success != NULL); 1398 1399 *success = FALSE; 1400 1401 /* caller is responsible to have checked whether we can estimate, i.e., expression curvature and overestimate flag match */ 1402 assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); /* if underestimate, then must be concave */ 1403 assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX); /* if overestimate, then must be convex */ 1404 1405 #ifdef SCIP_DEBUG 1406 SCIPinfoMessage(scip, NULL, "%sestimate expression ", overestimate ? "over" : "under"); 1407 SCIPprintExpr(scip, nlhdlrexprdata->nlexpr, NULL); 1408 SCIPinfoMessage(scip, NULL, " at point\n"); 1409 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1410 { 1411 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]); 1412 assert(var != NULL); 1413 1414 SCIPinfoMessage(scip, NULL, " <%s> = %g [%g,%g]\n", SCIPvarGetName(var), 1415 usemidpoint ? 0.5 * (SCIPvarGetLbLocal(var) + SCIPvarGetUbLocal(var)) : SCIPgetSolVal(scip, sol, var), 1416 SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var)); 1417 } 1418 #endif 1419 1420 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 1421 assert(nlhdlrdata != NULL); 1422 1423 if( nlhdlrdata->evalsol == NULL ) 1424 { 1425 SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) ); 1426 } 1427 1428 evaldata.nlhdlrexprdata = nlhdlrexprdata; 1429 evaldata.evalsol = nlhdlrdata->evalsol; 1430 evaldata.scip = scip; 1431 1432 SCIP_CALL( SCIPallocBufferArray(scip, &xstar, nlhdlrexprdata->nleafs) ); 1433 SCIP_CALL( SCIPallocBufferArray(scip, &box, 2*nlhdlrexprdata->nleafs) ); 1434 1435 allfixed = TRUE; 1436 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1437 { 1438 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]); 1439 assert(var != NULL); 1440 1441 box[2*i] = SCIPvarGetLbLocal(var); 1442 if( SCIPisInfinity(scip, -box[2*i]) ) 1443 { 1444 SCIPdebugMsg(scip, "lower bound at -infinity, no estimate possible\n"); 1445 goto TERMINATE; 1446 } 1447 1448 box[2*i+1] = SCIPvarGetUbLocal(var); 1449 if( SCIPisInfinity(scip, box[2*i+1]) ) 1450 { 1451 SCIPdebugMsg(scip, "upper bound at +infinity, no estimate possible\n"); 1452 goto TERMINATE; 1453 } 1454 1455 if( !SCIPisRelEQ(scip, box[2*i], box[2*i+1]) ) 1456 allfixed = FALSE; 1457 1458 if( usemidpoint ) 1459 xstar[i] = 0.5 * (box[2*i] + box[2*i+1]); 1460 else 1461 xstar[i] = SCIPgetSolVal(scip, sol, var); 1462 assert(xstar[i] != SCIP_INVALID); 1463 } 1464 1465 if( allfixed ) 1466 { 1467 /* SCIPcomputeFacetVertexPolyhedralNonlinear prints a warning and does not succeed if all is fixed */ 1468 SCIPdebugMsg(scip, "all variables fixed, skip estimate\n"); 1469 goto TERMINATE; 1470 } 1471 1472 SCIP_CALL( SCIPensureRowprepSize(scip, rowprep, nlhdlrexprdata->nleafs + 1) ); 1473 1474 SCIP_CALL( SCIPcomputeFacetVertexPolyhedralNonlinear(scip, conshdlr, overestimate, nlhdlrExprEvalConcave, (void*)&evaldata, 1475 xstar, box, nlhdlrexprdata->nleafs, targetvalue, success, SCIProwprepGetCoefs(rowprep), &facetconstant) ); 1476 1477 if( !*success ) 1478 { 1479 SCIPdebugMsg(scip, "failed to compute facet of convex hull\n"); 1480 goto TERMINATE; 1481 } 1482 1483 SCIProwprepSetLocal(rowprep, TRUE); 1484 SCIProwprepAddConstant(rowprep, facetconstant); 1485 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1486 { 1487 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]), SCIProwprepGetCoefs(rowprep)[i]) ); 1488 } 1489 1490 #ifdef SCIP_DEBUG 1491 SCIPinfoMessage(scip, NULL, "computed estimator: "); 1492 SCIPprintRowprep(scip, rowprep, NULL); 1493 #endif 1494 1495 TERMINATE: 1496 SCIPfreeBufferArray(scip, &box); 1497 SCIPfreeBufferArray(scip, &xstar); 1498 1499 return SCIP_OKAY; 1500 } 1501 1502 /** adds an estimator computed via a gradient to a given rowprep */ 1503 static 1504 SCIP_RETCODE estimateGradientInner( 1505 SCIP* scip, /**< SCIP data structure */ 1506 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 1507 SCIP_SOL* sol, /**< solution to use */ 1508 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 1509 SCIP_Bool* success /**< buffer to store whether successful */ 1510 ) 1511 { 1512 SCIP_EXPR* nlexpr; 1513 SCIP_Real QUAD(constant); 1514 int i; 1515 1516 assert(scip != NULL); 1517 assert(nlhdlrexprdata != NULL); 1518 assert(rowprep != NULL); 1519 assert(success != NULL); 1520 1521 nlexpr = nlhdlrexprdata->nlexpr; 1522 assert(nlexpr != NULL); 1523 1524 /* compute gradient (TODO: this also re-evaluates (soltag=0), which shouldn't be necessary unless we tried ConvexSecant before or are called from Sollinearize callback) */ 1525 SCIP_CALL( SCIPevalExprGradient(scip, nlexpr, sol, 0L) ); 1526 1527 /* if gradient evaluation error, then return */ 1528 if( SCIPexprGetDerivative(nlexpr) == SCIP_INVALID ) 1529 { 1530 SCIPdebugMsg(scip, "gradient evaluation error for %p\n", (void*)nlexpr); 1531 return SCIP_OKAY; 1532 } 1533 1534 /* add gradient underestimator to rowprep: f(sol) + (x - sol) \nabla f(sol) 1535 * constant will store f(sol) - sol * \nabla f(sol) 1536 * to avoid some cancellation errors when linear variables take huge values (like 1e20), 1537 * we use double-double arithmetic here 1538 */ 1539 QUAD_ASSIGN(constant, SCIPexprGetEvalValue(nlexpr)); /* f(sol) */ 1540 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1541 { 1542 SCIP_VAR* var; 1543 SCIP_Real deriv; 1544 SCIP_Real varval; 1545 1546 assert(SCIPexprGetDiffTag(nlhdlrexprdata->leafexprs[i]) == SCIPexprGetDiffTag(nlexpr)); 1547 deriv = SCIPexprGetDerivative(nlhdlrexprdata->leafexprs[i]); 1548 if( deriv == SCIP_INVALID ) 1549 { 1550 SCIPdebugMsg(scip, "gradient evaluation error for component %d of %p\n", i, (void*)nlexpr); 1551 return SCIP_OKAY; 1552 } 1553 1554 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]); 1555 assert(var != NULL); 1556 1557 varval = SCIPgetSolVal(scip, sol, var); 1558 1559 SCIPdebugMsg(scip, "add %g * (<%s> - %g) to rowprep\n", deriv, SCIPvarGetName(var), varval); 1560 1561 /* add deriv * var to rowprep and deriv * (-varval) to constant */ 1562 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, deriv) ); 1563 SCIPquadprecSumQD(constant, constant, -deriv * varval); 1564 } 1565 1566 SCIProwprepAddConstant(rowprep, QUAD_TO_DBL(constant)); 1567 SCIProwprepSetLocal(rowprep, FALSE); 1568 1569 *success = TRUE; 1570 1571 return SCIP_OKAY; 1572 } 1573 1574 /** adds an estimator computed via a gradient to a given rowprep, possibly perturbing solution */ 1575 static 1576 SCIP_RETCODE estimateGradient( 1577 SCIP* scip, /**< SCIP data structure */ 1578 SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */ 1579 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 1580 SCIP_SOL* sol, /**< solution to use */ 1581 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 1582 SCIP_Bool* success /**< buffer to store whether successful */ 1583 ) 1584 { 1585 SCIP_NLHDLRDATA* nlhdlrdata; 1586 int i; 1587 1588 assert(nlhdlrexprdata != NULL); 1589 assert(rowprep != NULL); 1590 assert(success != NULL); 1591 1592 #ifdef SCIP_DEBUG 1593 SCIPinfoMessage(scip, NULL, "estimate expression "); 1594 SCIPprintExpr(scip, nlhdlrexprdata->nlexpr, NULL); 1595 SCIPinfoMessage(scip, NULL, " by gradient\n"); 1596 #endif 1597 1598 *success = FALSE; 1599 1600 SCIP_CALL( estimateGradientInner(scip, nlhdlrexprdata, sol, rowprep, success) ); 1601 1602 /* if succeeded, then there was no gradient evaluation error, so we are done */ 1603 if( *success ) 1604 return SCIP_OKAY; 1605 1606 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 1607 assert(nlhdlrdata != NULL); 1608 1609 /* we take maxperturb == 0 as signal to not try perturbation */ 1610 if( nlhdlrdata->maxperturb == 0.0 ) 1611 { 1612 SCIPdebugMsg(scip, "gradient evaluation error, perturbation disabled\n"); 1613 return SCIP_OKAY; 1614 } 1615 1616 SCIPdebugMsg(scip, "gradient evaluation error, try perturbed point\n"); 1617 1618 if( nlhdlrdata->evalsol == NULL ) 1619 { 1620 SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) ); 1621 } 1622 1623 if( nlhdlrdata->randnumgen == NULL ) 1624 { 1625 SCIP_CALL( SCIPcreateRandom(scip, &nlhdlrdata->randnumgen, RANDNUMINITSEED, TRUE) ); 1626 } 1627 1628 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1629 { 1630 SCIP_VAR* var; 1631 SCIP_Real lb; 1632 SCIP_Real ub; 1633 SCIP_Real val; 1634 SCIP_Real p; 1635 1636 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]); 1637 assert(var != NULL); 1638 1639 lb = SCIPvarGetLbGlobal(var); 1640 ub = SCIPvarGetUbGlobal(var); 1641 val = SCIPgetSolVal(scip, sol, var); 1642 val = MIN(ub, MAX(lb, val)); 1643 1644 p = SCIPrandomGetReal(nlhdlrdata->randnumgen, -nlhdlrdata->maxperturb, nlhdlrdata->maxperturb); 1645 if( !SCIPisZero(scip, val) ) 1646 p *= REALABS(val); 1647 1648 /* if perturbation to left underruns lower bound, then try perturb to right 1649 * if perturbation to right exceeds lower bound, then try perturb to left 1650 */ 1651 if( val + p <= lb ) 1652 val -= p; 1653 else if( val + p >= ub ) 1654 val -= p; 1655 else 1656 val += p; 1657 1658 /* if still exceeding bound, then |ub-lb| < 2*maxperturb and we pick a random value between bounds 1659 * (if ub-lb < 2eps, then SCIPrandomGetReal() still gives a reasonable number in [lb,ub]) 1660 */ 1661 if( val <= lb || val >= ub ) 1662 val = SCIPrandomGetReal(nlhdlrdata->randnumgen, lb + SCIPepsilon(scip), ub - SCIPepsilon(scip)); 1663 1664 SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, val) ); 1665 } 1666 1667 /* try again with perturbed point */ 1668 SCIP_CALL( estimateGradientInner(scip, nlhdlrexprdata, nlhdlrdata->evalsol, rowprep, success) ); 1669 1670 return SCIP_OKAY; 1671 } 1672 1673 /** adds an estimator generated by putting a secant through the coordinates given by the two closest integer points */ 1674 static 1675 SCIP_RETCODE estimateConvexSecant( 1676 SCIP* scip, /**< SCIP data structure */ 1677 SCIP_NLHDLR* nlhdlr, /**< nonlinear handler */ 1678 SCIP_NLHDLREXPRDATA* nlhdlrexprdata, /**< nonlinear handler expression data */ 1679 SCIP_SOL* sol, /**< solution to use, unless usemidpoint is TRUE */ 1680 SCIP_ROWPREP* rowprep, /**< rowprep where to store estimator */ 1681 SCIP_Bool* success /**< buffer to store whether successful */ 1682 ) 1683 { 1684 SCIP_NLHDLRDATA* nlhdlrdata; 1685 SCIP_EXPR* nlexpr; 1686 SCIP_VAR* var; 1687 SCIP_Real x; 1688 SCIP_Real left, right; 1689 SCIP_Real fleft, fright; 1690 1691 assert(nlhdlrexprdata != NULL); 1692 assert(nlhdlrexprdata->nleafs == 1); 1693 assert(rowprep != NULL); 1694 assert(success != NULL); 1695 1696 nlexpr = nlhdlrexprdata->nlexpr; 1697 assert(nlexpr != NULL); 1698 1699 *success = FALSE; 1700 1701 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 1702 assert(nlhdlrdata != NULL); 1703 1704 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[0]); 1705 assert(var != NULL); 1706 1707 x = SCIPgetSolVal(scip, sol, var); 1708 1709 #ifdef SCIP_DEBUG 1710 SCIPinfoMessage(scip, NULL, "estimate expression "); 1711 SCIPprintExpr(scip, nlexpr, NULL); 1712 SCIPinfoMessage(scip, NULL, " by secant\n"); 1713 SCIPinfoMessage(scip, NULL, "integral variable <%s> = %g [%g,%g]\n", SCIPvarGetName(var), 1714 x, SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var)); 1715 #endif 1716 1717 /* find out coordinates of var left and right to sol */ 1718 if( SCIPisIntegral(scip, x) ) 1719 { 1720 x = SCIPround(scip, x); 1721 if( SCIPisEQ(scip, x, SCIPvarGetLbGlobal(var)) ) 1722 { 1723 left = x; 1724 right = left + 1.0; 1725 } 1726 else 1727 { 1728 right = x; 1729 left = right - 1.0; 1730 } 1731 } 1732 else 1733 { 1734 left = SCIPfloor(scip, x); 1735 right = SCIPceil(scip, x); 1736 } 1737 assert(left != right); 1738 1739 /* now evaluate at left and right */ 1740 if( nlhdlrdata->evalsol == NULL ) 1741 { 1742 SCIP_CALL( SCIPcreateSol(scip, &nlhdlrdata->evalsol, NULL) ); 1743 } 1744 1745 SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, left) ); 1746 SCIP_CALL( SCIPevalExpr(scip, nlexpr, nlhdlrdata->evalsol, 0L) ); 1747 1748 /* evaluation error or a too large constant -> skip */ 1749 fleft = SCIPexprGetEvalValue(nlexpr); 1750 if( SCIPisInfinity(scip, REALABS(fleft)) ) 1751 { 1752 SCIPdebugMsg(scip, "evaluation error / too large value (%g) for %p\n", SCIPexprGetEvalValue(nlexpr), (void*)nlexpr); 1753 return SCIP_OKAY; 1754 } 1755 1756 SCIP_CALL( SCIPsetSolVal(scip, nlhdlrdata->evalsol, var, right) ); 1757 SCIP_CALL( SCIPevalExpr(scip, nlexpr, nlhdlrdata->evalsol, 0L) ); 1758 1759 /* evaluation error or a too large constant -> skip */ 1760 fright = SCIPexprGetEvalValue(nlexpr); 1761 if( SCIPisInfinity(scip, REALABS(fright)) ) 1762 { 1763 SCIPdebugMsg(scip, "evaluation error / too large value (%g) for %p\n", SCIPexprGetEvalValue(nlexpr), (void*)nlexpr); 1764 return SCIP_OKAY; 1765 } 1766 1767 SCIPdebugMsg(scip, "f(%g)=%g, f(%g)=%g\n", left, fleft, right, fright); 1768 1769 /* skip if too steep 1770 * for clay0204h, this resulted in a wrong cut from f(0)=1e12 f(1)=0.99998, 1771 * since due to limited precision, this was handled as if f(1)=1 1772 */ 1773 if( (!SCIPisZero(scip, fleft) && REALABS(fright/fleft)*SCIPepsilon(scip) > 1.0) || 1774 (!SCIPisZero(scip, fright) && REALABS(fleft/fright)*SCIPepsilon(scip) > 1.0) ) 1775 { 1776 SCIPdebugMsg(scip, "function is too steep, abandoning\n"); 1777 return SCIP_OKAY; 1778 } 1779 1780 /* now add f(left) + (f(right) - f(left)) * (x - left) as estimator to rowprep */ 1781 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, var, fright - fleft) ); 1782 SCIProwprepAddConstant(rowprep, fleft - (fright - fleft) * left); 1783 SCIProwprepSetLocal(rowprep, FALSE); 1784 1785 *success = TRUE; 1786 1787 return SCIP_OKAY; 1788 } 1789 1790 /* 1791 * Callback methods of convex nonlinear handler 1792 */ 1793 1794 /** free handler data of convex or concave nlhdlr */ 1795 static 1796 SCIP_DECL_NLHDLRFREEHDLRDATA(nlhdlrfreeHdlrDataConvexConcave) 1797 { /*lint --e{715}*/ 1798 assert(scip != NULL); 1799 assert(nlhdlrdata != NULL); 1800 assert(*nlhdlrdata != NULL); 1801 assert((*nlhdlrdata)->evalsol == NULL); 1802 assert((*nlhdlrdata)->randnumgen == NULL); 1803 1804 SCIPfreeBlockMemory(scip, nlhdlrdata); 1805 1806 return SCIP_OKAY; 1807 } 1808 1809 /** callback to free expression specific data */ 1810 static 1811 SCIP_DECL_NLHDLRFREEEXPRDATA(nlhdlrfreeExprDataConvexConcave) 1812 { /*lint --e{715}*/ 1813 assert(scip != NULL); 1814 assert(nlhdlrexprdata != NULL); 1815 assert(*nlhdlrexprdata != NULL); 1816 1817 SCIPfreeBlockMemoryArrayNull(scip, &(*nlhdlrexprdata)->leafexprs, (*nlhdlrexprdata)->nleafs); 1818 SCIP_CALL( SCIPreleaseExpr(scip, &(*nlhdlrexprdata)->nlexpr) ); 1819 SCIPhashmapFree(&(*nlhdlrexprdata)->nlexpr2origexpr); 1820 1821 SCIPfreeBlockMemory(scip, nlhdlrexprdata); 1822 1823 return SCIP_OKAY; 1824 } 1825 1826 /** deinitialization of problem-specific data */ 1827 static 1828 SCIP_DECL_NLHDLREXIT(nlhdlrExitConvex) 1829 { 1830 SCIP_NLHDLRDATA* nlhdlrdata; 1831 1832 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 1833 assert(nlhdlrdata != NULL); 1834 1835 if( nlhdlrdata->evalsol != NULL ) 1836 { 1837 SCIP_CALL( SCIPfreeSol(scip, &nlhdlrdata->evalsol) ); 1838 } 1839 1840 if( nlhdlrdata->randnumgen != NULL ) 1841 SCIPfreeRandom(scip, &nlhdlrdata->randnumgen); 1842 1843 return SCIP_OKAY; 1844 } 1845 1846 /** checks whether expression (or -expression) is convex, possibly after introducing auxiliary variables */ 1847 static 1848 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectConvex) 1849 { /*lint --e{715}*/ 1850 SCIP_NLHDLRDATA* nlhdlrdata; 1851 SCIP_EXPR* nlexpr = NULL; 1852 SCIP_HASHMAP* nlexpr2origexpr; 1853 int nleafs = 0; 1854 1855 assert(scip != NULL); 1856 assert(nlhdlr != NULL); 1857 assert(expr != NULL); 1858 assert(enforcing != NULL); 1859 assert(participating != NULL); 1860 assert(nlhdlrexprdata != NULL); 1861 1862 /* we currently do not participate if only activity computation is required */ 1863 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH ) 1864 return SCIP_OKAY; 1865 1866 /* ignore pure constants and variables */ 1867 if( SCIPexprGetNChildren(expr) == 0 ) 1868 return SCIP_OKAY; 1869 1870 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 1871 assert(nlhdlrdata != NULL); 1872 assert(nlhdlrdata->isnlhdlrconvex); 1873 1874 SCIPdebugMsg(scip, "nlhdlr_convex detect for expr %p\n", (void*)expr); 1875 1876 /* initialize mapping from copied expression to original one 1877 * 20 is not a bad estimate for the size of convex subexpressions that we can usually discover 1878 * when expressions will be allowed to store "user"data, we could get rid of this hashmap (TODO) 1879 */ 1880 SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) ); 1881 1882 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 ) /* if no separation below yet */ 1883 { 1884 SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr, 1885 SCIP_EXPRCURV_CONVEX, NULL, SCIPassumeConvexNonlinear(conshdlr), NULL) ); 1886 if( nlexpr != NULL ) 1887 { 1888 assert(SCIPexprGetNChildren(nlexpr) > 0); /* should not be trivial */ 1889 1890 *participating |= SCIP_NLHDLR_METHOD_SEPABELOW; 1891 1892 SCIPdebugMsg(scip, "detected expr %p to be convex -> can enforce expr <= auxvar\n", (void*)expr); 1893 } 1894 else 1895 { 1896 SCIP_CALL( SCIPhashmapRemoveAll(nlexpr2origexpr) ); 1897 } 1898 } 1899 1900 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 && nlexpr == NULL ) /* if no separation above and not convex */ 1901 { 1902 SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr, 1903 SCIP_EXPRCURV_CONCAVE, NULL, SCIPassumeConvexNonlinear(conshdlr), NULL) ); 1904 if( nlexpr != NULL ) 1905 { 1906 assert(SCIPexprGetNChildren(nlexpr) > 0); /* should not be trivial */ 1907 1908 *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE; 1909 1910 SCIPdebugMsg(scip, "detected expr %p to be concave -> can enforce expr >= auxvar\n", (void*)expr); 1911 } 1912 } 1913 1914 /* everything we participate in we also enforce */ 1915 *enforcing |= *participating; 1916 1917 assert(*participating || nlexpr == NULL); 1918 if( !*participating ) 1919 { 1920 SCIPhashmapFree(&nlexpr2origexpr); 1921 return SCIP_OKAY; 1922 } 1923 1924 /* create the expression data of the nonlinear handler 1925 * notify conshdlr about expr for which we will require auxiliary variables 1926 */ 1927 SCIP_CALL( createNlhdlrExprData(scip, nlhdlrdata, nlhdlrexprdata, expr, nlexpr, nlexpr2origexpr, nleafs, *participating) ); 1928 1929 return SCIP_OKAY; 1930 } 1931 1932 /** auxiliary evaluation callback */ 1933 static 1934 SCIP_DECL_NLHDLREVALAUX(nlhdlrEvalAuxConvexConcave) 1935 { /*lint --e{715}*/ 1936 assert(nlhdlrexprdata != NULL); 1937 assert(nlhdlrexprdata->nlexpr != NULL); 1938 assert(auxvalue != NULL); 1939 1940 SCIP_CALL( SCIPevalExpr(scip, nlhdlrexprdata->nlexpr, sol, 0L) ); 1941 *auxvalue = SCIPexprGetEvalValue(nlhdlrexprdata->nlexpr); 1942 1943 return SCIP_OKAY; 1944 } 1945 1946 /** init sepa callback that initializes LP */ 1947 static 1948 SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaConvex) 1949 { /*lint --e{715}*/ 1950 SCIP_EXPR* nlexpr; 1951 SCIP_EXPRCURV curvature; 1952 SCIP_Bool success; 1953 SCIP_ROWPREP* rowprep = NULL; 1954 SCIP_ROW* row; 1955 SCIP_Real lb; 1956 SCIP_Real ub; 1957 SCIP_Real lambda; 1958 SCIP_SOL* sol; 1959 int k; 1960 1961 assert(scip != NULL); 1962 assert(expr != NULL); 1963 assert(nlhdlrexprdata != NULL); 1964 1965 /* setup nlhdlrexprdata->leafexprs */ 1966 SCIP_CALL( collectLeafs(scip, nlhdlrexprdata) ); 1967 1968 nlexpr = nlhdlrexprdata->nlexpr; 1969 assert(nlexpr != NULL); 1970 assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlexpr) == expr); 1971 1972 curvature = SCIPexprGetCurvature(nlexpr); 1973 assert(curvature == SCIP_EXPRCURV_CONVEX || curvature == SCIP_EXPRCURV_CONCAVE); 1974 1975 /* we can only be estimating on the convex side */ 1976 if( curvature == SCIP_EXPRCURV_CONVEX ) 1977 overestimate = FALSE; 1978 else if( curvature == SCIP_EXPRCURV_CONCAVE ) 1979 underestimate = FALSE; 1980 if( !overestimate && !underestimate ) 1981 return SCIP_OKAY; 1982 1983 /* linearizes at 5 different points obtained as convex combination of the lower and upper bound of the variables 1984 * present in the convex expression; whether more weight is given to the lower or upper bound of a variable depends 1985 * on whether the fixing of the variable to that value is better for the objective function 1986 */ 1987 SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) ); 1988 1989 *infeasible = FALSE; 1990 1991 for( k = 0; k < 5; ++k ) 1992 { 1993 int i; 1994 lambda = 0.1 * (k+1); /* lambda = 0.1, 0.2, 0.3, 0.4, 0.5 */ 1995 1996 for( i = 0; i < nlhdlrexprdata->nleafs; ++i ) 1997 { 1998 SCIP_VAR* var; 1999 2000 var = SCIPgetVarExprVar(nlhdlrexprdata->leafexprs[i]); 2001 2002 lb = SCIPvarGetLbGlobal(var); 2003 ub = SCIPvarGetUbGlobal(var); 2004 2005 /* make sure the absolute values of bounds are not too large */ 2006 if( ub > -INITLPMAXVARVAL ) 2007 lb = MAX(lb, -INITLPMAXVARVAL); 2008 if( lb < INITLPMAXVARVAL ) 2009 ub = MIN(ub, INITLPMAXVARVAL); 2010 2011 /* in the case when ub < -maxabsbnd or lb > maxabsbnd, we still want to at least make bounds finite */ 2012 if( SCIPisInfinity(scip, -lb) ) 2013 lb = MIN(-10.0, ub - 0.1*REALABS(ub)); 2014 if( SCIPisInfinity(scip, ub) ) 2015 ub = MAX( 10.0, lb + 0.1*REALABS(lb)); 2016 2017 if( SCIPvarGetBestBoundType(var) == SCIP_BOUNDTYPE_LOWER ) 2018 SCIP_CALL( SCIPsetSolVal(scip, sol, var, lambda * ub + (1.0 - lambda) * lb) ); 2019 else 2020 SCIP_CALL( SCIPsetSolVal(scip, sol, var, lambda * lb + (1.0 - lambda) * ub) ); 2021 } 2022 2023 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 2024 SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) ); 2025 if( !success ) 2026 { 2027 SCIPdebugMsg(scip, "failed to linearize for k = %d\n", k); 2028 SCIPfreeRowprep(scip, &rowprep); 2029 continue; 2030 } 2031 2032 /* add auxiliary variable */ 2033 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) ); 2034 2035 /* straighten out numerics */ 2036 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2037 if( !success ) 2038 { 2039 SCIPdebugMsg(scip, "failed to cleanup rowprep numerics for k = %d\n", k); 2040 SCIPfreeRowprep(scip, &rowprep); 2041 continue; 2042 } 2043 2044 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_gradient%p_initsepa_%d", 2045 overestimate ? "over" : "under", (void*)expr, k); 2046 SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) ); 2047 SCIPfreeRowprep(scip, &rowprep); 2048 2049 #ifdef SCIP_DEBUG 2050 SCIPinfoMessage(scip, NULL, "initsepa computed row: "); 2051 SCIPprintRow(scip, row, NULL); 2052 #endif 2053 2054 SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) ); 2055 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 2056 2057 if( *infeasible ) 2058 break; 2059 } 2060 2061 SCIP_CALL( SCIPfreeSol(scip, &sol) ); 2062 2063 return SCIP_OKAY; 2064 } 2065 2066 /** estimator callback */ 2067 static 2068 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateConvex) 2069 { /*lint --e{715}*/ 2070 SCIP_ROWPREP* rowprep; 2071 2072 assert(scip != NULL); 2073 assert(expr != NULL); 2074 assert(nlhdlrexprdata != NULL); 2075 assert(nlhdlrexprdata->nlexpr != NULL); 2076 assert(rowpreps != NULL); 2077 assert(success != NULL); 2078 2079 assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr); 2080 2081 /* we must be called only for the side that we indicated to participate in during DETECT */ 2082 assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX 2083 || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2084 assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2085 assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX); 2086 2087 *success = FALSE; 2088 *addedbranchscores = FALSE; 2089 2090 /* we can skip eval as nlhdlrEvalAux should have been called for same solution before */ 2091 /* SCIP_CALL( nlhdlrExprEval(scip, nlexpr, sol) ); */ 2092 assert(auxvalue == SCIPexprGetEvalValue(nlhdlrexprdata->nlexpr)); /* given value (originally from 2093 nlhdlrEvalAuxConvexConcave) should coincide with the one stored in nlexpr */ 2094 2095 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 2096 2097 if( nlhdlrexprdata->nleafs == 1 && SCIPexprIsIntegral(nlhdlrexprdata->leafexprs[0]) ) 2098 { 2099 SCIP_CALL( estimateConvexSecant(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, success) ); 2100 2101 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexsecant%p_%s%" SCIP_LONGINT_FORMAT, 2102 overestimate ? "over" : "under", 2103 (void*)expr, 2104 sol != NULL ? "sol" : "lp", 2105 sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip)); 2106 } 2107 2108 /* if secant method was not used or failed, then try with gradient (unless we had an evaluation error in sol before) */ 2109 if( !*success && auxvalue != SCIP_INVALID ) 2110 { 2111 SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, success) ); 2112 2113 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexgradient%p_%s%" SCIP_LONGINT_FORMAT, 2114 overestimate ? "over" : "under", 2115 (void*)expr, 2116 sol != NULL ? "sol" : "lp", 2117 sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip)); 2118 } 2119 2120 if( *success ) 2121 { 2122 SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) ); 2123 } 2124 else 2125 { 2126 SCIPfreeRowprep(scip, &rowprep); 2127 } 2128 2129 return SCIP_OKAY; 2130 } 2131 2132 /** solution notification callback */ 2133 static 2134 SCIP_DECL_NLHDLRSOLLINEARIZE(nlhdlrSollinearizeConvex) 2135 { /*lint --e{715}*/ 2136 SCIP_ROWPREP* rowprep; 2137 SCIP_Bool success = FALSE; 2138 2139 assert(scip != NULL); 2140 assert(expr != NULL); 2141 assert(nlhdlrexprdata != NULL); 2142 assert(nlhdlrexprdata->nlexpr != NULL); 2143 assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr); 2144 assert(sol != NULL); 2145 2146 /* we must be called only for the side that we indicated to participate in during DETECT */ 2147 assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX 2148 || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2149 assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2150 assert(!underestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX); 2151 assert(underestimate == !overestimate); /* should be exactly one of under- and overestimate */ 2152 2153 /* evaluate nlexpr in solution */ 2154 SCIP_CALL( SCIPevalExpr(scip, nlhdlrexprdata->nlexpr, sol, 0L) ); 2155 2156 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 2157 2158 if( nlhdlrexprdata->nleafs == 1 && SCIPexprIsIntegral(nlhdlrexprdata->leafexprs[0]) ) 2159 { 2160 SCIP_CALL( estimateConvexSecant(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) ); 2161 2162 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexsecant%p_sol%dnotify", 2163 overestimate ? "over" : "under", (void*)expr, SCIPsolGetIndex(sol)); 2164 } 2165 2166 /* if secant method was not used or failed, then try with gradient */ 2167 if( !success ) 2168 { 2169 SCIP_CALL( estimateGradient(scip, nlhdlr, nlhdlrexprdata, sol, rowprep, &success) ); 2170 2171 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_convexgradient%p_sol%dnotify", 2172 overestimate ? "over" : "under", (void*)expr, SCIPsolGetIndex(sol)); 2173 } 2174 2175 if( success ) 2176 { 2177 /* complete estimator to cut and clean it up */ 2178 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) ); 2179 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, sol, SCIPgetHugeValue(scip), &success) ); 2180 } 2181 2182 /* if cleanup succeeded and rowprep is still global, add to cutpool */ 2183 if( success && !SCIProwprepIsLocal(rowprep) ) 2184 { 2185 SCIP_ROW* row; 2186 2187 SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) ); 2188 SCIP_CALL( SCIPaddPoolCut(scip, row) ); 2189 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 2190 } 2191 2192 SCIPfreeRowprep(scip, &rowprep); 2193 2194 return SCIP_OKAY; 2195 } 2196 2197 /** include nlhdlr in another scip instance */ 2198 static 2199 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrConvex) 2200 { /*lint --e{715}*/ 2201 assert(targetscip != NULL); 2202 assert(sourcenlhdlr != NULL); 2203 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), CONVEX_NLHDLR_NAME) == 0); 2204 2205 SCIP_CALL( SCIPincludeNlhdlrConvex(targetscip) ); 2206 2207 return SCIP_OKAY; 2208 } 2209 2210 /** includes convex nonlinear handler in nonlinear constraint handler */ 2211 SCIP_RETCODE SCIPincludeNlhdlrConvex( 2212 SCIP* scip /**< SCIP data structure */ 2213 ) 2214 { 2215 SCIP_NLHDLR* nlhdlr; 2216 SCIP_NLHDLRDATA* nlhdlrdata; 2217 2218 assert(scip != NULL); 2219 2220 SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) ); 2221 nlhdlrdata->isnlhdlrconvex = TRUE; 2222 nlhdlrdata->evalsol = NULL; 2223 nlhdlrdata->randnumgen = NULL; 2224 2225 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, CONVEX_NLHDLR_NAME, CONVEX_NLHDLR_DESC, 2226 CONVEX_NLHDLR_DETECTPRIORITY, CONVEX_NLHDLR_ENFOPRIORITY, nlhdlrDetectConvex, nlhdlrEvalAuxConvexConcave, nlhdlrdata) ); 2227 assert(nlhdlr != NULL); 2228 2229 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/detectsum", 2230 "whether to run convexity detection when the root of an expression is a non-quadratic sum", 2231 &nlhdlrdata->detectsum, FALSE, DEFAULT_DETECTSUM, NULL, NULL) ); 2232 2233 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/extendedform", 2234 "whether to create extended formulations instead of looking for maximal convex expressions", 2235 &nlhdlrdata->extendedform, FALSE, DEFAULT_EXTENDEDFORM, NULL, NULL) ); 2236 2237 SCIP_CALL( SCIPaddRealParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/maxperturb", 2238 "maximal relative perturbation of non-differentiable reference point", 2239 &nlhdlrdata->maxperturb, FALSE, DEFAULT_MAXPERTURB, 0.0, 1.0, NULL, NULL) ); 2240 2241 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxquadratic", 2242 "whether to use convexity check on quadratics", 2243 &nlhdlrdata->cvxquadratic, TRUE, DEFAULT_CVXQUADRATIC_CONVEX, NULL, NULL) ); 2244 2245 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxsignomial", 2246 "whether to use convexity check on signomials", 2247 &nlhdlrdata->cvxsignomial, TRUE, DEFAULT_CVXSIGNOMIAL, NULL, NULL) ); 2248 2249 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/cvxprodcomp", 2250 "whether to use convexity check on product composition f(h)*h", 2251 &nlhdlrdata->cvxprodcomp, TRUE, DEFAULT_CVXPRODCOMP, NULL, NULL) ); 2252 2253 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONVEX_NLHDLR_NAME "/handletrivial", 2254 "whether to also handle trivial convex expressions", 2255 &nlhdlrdata->handletrivial, TRUE, DEFAULT_HANDLETRIVIAL, NULL, NULL) ); 2256 2257 SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrfreeHdlrDataConvexConcave); 2258 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrConvex); 2259 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrfreeExprDataConvexConcave); 2260 SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaConvex, NULL, nlhdlrEstimateConvex, NULL); 2261 SCIPnlhdlrSetSollinearize(nlhdlr, nlhdlrSollinearizeConvex); 2262 SCIPnlhdlrSetInitExit(nlhdlr, NULL, nlhdlrExitConvex); 2263 2264 return SCIP_OKAY; 2265 } 2266 2267 /* 2268 * Callback methods of concave nonlinear handler 2269 */ 2270 2271 /** deinitialization of problem-specific data */ 2272 static 2273 SCIP_DECL_NLHDLREXIT(nlhdlrExitConcave) 2274 { 2275 SCIP_NLHDLRDATA* nlhdlrdata; 2276 2277 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 2278 assert(nlhdlrdata != NULL); 2279 assert(nlhdlrdata->randnumgen == NULL); /* not used for concave nlhdlr so far */ 2280 2281 if( nlhdlrdata->evalsol != NULL ) 2282 { 2283 SCIP_CALL( SCIPfreeSol(scip, &nlhdlrdata->evalsol) ); 2284 } 2285 2286 return SCIP_OKAY; 2287 } 2288 2289 /** checks whether expression (or -expression) is concave, possibly after introducing auxiliary variables */ 2290 static 2291 SCIP_DECL_NLHDLRDETECT(nlhdlrDetectConcave) 2292 { /*lint --e{715}*/ 2293 SCIP_NLHDLRDATA* nlhdlrdata; 2294 SCIP_EXPR* nlexpr = NULL; 2295 SCIP_HASHMAP* nlexpr2origexpr; 2296 int nleafs = 0; 2297 2298 assert(scip != NULL); 2299 assert(nlhdlr != NULL); 2300 assert(expr != NULL); 2301 assert(enforcing != NULL); 2302 assert(participating != NULL); 2303 assert(nlhdlrexprdata != NULL); 2304 2305 /* we currently do not participate if only activity computation is required */ 2306 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABOTH) == SCIP_NLHDLR_METHOD_SEPABOTH ) 2307 return SCIP_OKAY; 2308 2309 /* ignore pure constants and variables */ 2310 if( SCIPexprGetNChildren(expr) == 0 ) 2311 return SCIP_OKAY; 2312 2313 nlhdlrdata = SCIPnlhdlrGetData(nlhdlr); 2314 assert(nlhdlrdata != NULL); 2315 assert(!nlhdlrdata->isnlhdlrconvex); 2316 2317 SCIPdebugMsg(scip, "nlhdlr_concave detect for expr %p\n", (void*)expr); 2318 2319 /* initialize mapping from copied expression to original one 2320 * 20 is not a bad estimate for the size of concave subexpressions that we can usually discover 2321 * when expressions will be allowed to store "user"data, we could get rid of this hashmap (TODO) 2322 */ 2323 SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) ); 2324 2325 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPABELOW) == 0 ) /* if no separation below yet */ 2326 { 2327 SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr, 2328 SCIP_EXPRCURV_CONCAVE, NULL, FALSE, NULL) ); 2329 2330 if( nlexpr != NULL && nleafs > SCIP_MAXVERTEXPOLYDIM ) 2331 { 2332 SCIPdebugMsg(scip, "Too many variables (%d) in constructed expression. Will not be able to estimate. Rejecting.\n", nleafs); 2333 SCIP_CALL( SCIPreleaseExpr(scip, &nlexpr) ); 2334 } 2335 2336 if( nlexpr != NULL ) 2337 { 2338 assert(SCIPexprGetNChildren(nlexpr) > 0); /* should not be trivial */ 2339 2340 *participating |= SCIP_NLHDLR_METHOD_SEPABELOW; 2341 2342 SCIPdebugMsg(scip, "detected expr %p to be concave -> can enforce expr <= auxvar\n", (void*)expr); 2343 } 2344 else 2345 { 2346 SCIP_CALL( SCIPhashmapRemoveAll(nlexpr2origexpr) ); 2347 } 2348 } 2349 2350 if( (*enforcing & SCIP_NLHDLR_METHOD_SEPAABOVE) == 0 && nlexpr == NULL ) /* if no separation above and not concave */ 2351 { 2352 SCIP_CALL( constructExpr(scip, nlhdlrdata, &nlexpr, nlexpr2origexpr, &nleafs, expr, 2353 SCIP_EXPRCURV_CONVEX, NULL, FALSE, NULL) ); 2354 2355 if( nlexpr != NULL && nleafs > SCIP_MAXVERTEXPOLYDIM ) 2356 { 2357 SCIPdebugMsg(scip, "Too many variables (%d) in constructed expression. Will not be able to estimate. Rejecting.\n", nleafs); 2358 SCIP_CALL( SCIPreleaseExpr(scip, &nlexpr) ); 2359 } 2360 2361 if( nlexpr != NULL ) 2362 { 2363 assert(SCIPexprGetNChildren(nlexpr) > 0); /* should not be trivial */ 2364 2365 *participating |= SCIP_NLHDLR_METHOD_SEPAABOVE; 2366 2367 SCIPdebugMsg(scip, "detected expr %p to be convex -> can enforce expr >= auxvar\n", (void*)expr); 2368 } 2369 } 2370 2371 /* everything we participate in we also enforce (at the moment) */ 2372 *enforcing |= *participating; 2373 2374 assert(*participating || nlexpr == NULL); 2375 if( !*participating ) 2376 { 2377 SCIPhashmapFree(&nlexpr2origexpr); 2378 return SCIP_OKAY; 2379 } 2380 2381 /* create the expression data of the nonlinear handler 2382 * notify conshdlr about expr for which we will require auxiliary variables and use activity 2383 */ 2384 SCIP_CALL( createNlhdlrExprData(scip, nlhdlrdata, nlhdlrexprdata, expr, nlexpr, nlexpr2origexpr, nleafs, *participating) ); 2385 2386 return SCIP_OKAY; 2387 } 2388 2389 /** init sepa callback that initializes LP */ 2390 static 2391 SCIP_DECL_NLHDLRINITSEPA(nlhdlrInitSepaConcave) 2392 { 2393 SCIP_EXPR* nlexpr; 2394 SCIP_EXPRCURV curvature; 2395 SCIP_Bool success; 2396 SCIP_ROWPREP* rowprep = NULL; 2397 SCIP_ROW* row; 2398 2399 assert(scip != NULL); 2400 assert(expr != NULL); 2401 assert(nlhdlrexprdata != NULL); 2402 2403 nlexpr = nlhdlrexprdata->nlexpr; 2404 assert(nlexpr != NULL); 2405 assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlexpr) == expr); 2406 2407 /* setup nlhdlrexprdata->leafexprs */ 2408 SCIP_CALL( collectLeafs(scip, nlhdlrexprdata) ); 2409 2410 curvature = SCIPexprGetCurvature(nlexpr); 2411 assert(curvature == SCIP_EXPRCURV_CONVEX || curvature == SCIP_EXPRCURV_CONCAVE); 2412 /* we can only be estimating on non-convex side */ 2413 if( curvature == SCIP_EXPRCURV_CONCAVE ) 2414 overestimate = FALSE; 2415 else if( curvature == SCIP_EXPRCURV_CONVEX ) 2416 underestimate = FALSE; 2417 if( !overestimate && !underestimate ) 2418 return SCIP_OKAY; 2419 2420 /* compute estimator and store in rowprep */ 2421 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 2422 SCIP_CALL( estimateVertexPolyhedral(scip, conshdlr, nlhdlr, nlhdlrexprdata, NULL, TRUE, overestimate, 2423 overestimate ? SCIPinfinity(scip) : -SCIPinfinity(scip), rowprep, &success) ); 2424 if( !success ) 2425 { 2426 SCIPdebugMsg(scip, "failed to compute facet of convex hull\n"); 2427 goto TERMINATE; 2428 } 2429 2430 /* add auxiliary variable */ 2431 SCIP_CALL( SCIPaddRowprepTerm(scip, rowprep, SCIPgetExprAuxVarNonlinear(expr), -1.0) ); 2432 2433 /* straighten out numerics */ 2434 SCIP_CALL( SCIPcleanupRowprep2(scip, rowprep, NULL, SCIPgetHugeValue(scip), &success) ); 2435 if( !success ) 2436 { 2437 SCIPdebugMsg(scip, "failed to cleanup rowprep numerics\n"); 2438 goto TERMINATE; 2439 } 2440 2441 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_concave%p_initsepa", 2442 overestimate ? "over" : "under", (void*)expr); 2443 SCIP_CALL( SCIPgetRowprepRowCons(scip, &row, rowprep, cons) ); 2444 2445 #ifdef SCIP_DEBUG 2446 SCIPinfoMessage(scip, NULL, "initsepa computed row: "); 2447 SCIPprintRow(scip, row, NULL); 2448 #endif 2449 2450 SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) ); 2451 SCIP_CALL( SCIPreleaseRow(scip, &row) ); 2452 2453 TERMINATE: 2454 if( rowprep != NULL ) 2455 SCIPfreeRowprep(scip, &rowprep); 2456 2457 return SCIP_OKAY; 2458 } 2459 2460 /** estimator callback */ 2461 static 2462 SCIP_DECL_NLHDLRESTIMATE(nlhdlrEstimateConcave) 2463 { /*lint --e{715}*/ 2464 SCIP_ROWPREP* rowprep; 2465 2466 assert(scip != NULL); 2467 assert(expr != NULL); 2468 assert(nlhdlrexprdata != NULL); 2469 assert(nlhdlrexprdata->nlexpr != NULL); 2470 assert(rowpreps != NULL); 2471 assert(success != NULL); 2472 2473 assert(SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, (void*)nlhdlrexprdata->nlexpr) == expr); 2474 2475 /* we must be called only for the side that we indicated to participate in during DETECT */ 2476 assert(SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX 2477 || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2478 assert(!overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONVEX); 2479 assert( overestimate || SCIPexprGetCurvature(nlhdlrexprdata->nlexpr) == SCIP_EXPRCURV_CONCAVE); 2480 2481 *success = FALSE; 2482 *addedbranchscores = FALSE; 2483 2484 SCIP_CALL( SCIPcreateRowprep(scip, &rowprep, overestimate ? SCIP_SIDETYPE_LEFT : SCIP_SIDETYPE_RIGHT, TRUE) ); 2485 2486 SCIP_CALL( estimateVertexPolyhedral(scip, conshdlr, nlhdlr, nlhdlrexprdata, sol, FALSE, overestimate, targetvalue, rowprep, success) ); 2487 2488 if( *success ) 2489 { 2490 SCIP_CALL( SCIPsetPtrarrayVal(scip, rowpreps, 0, rowprep) ); 2491 2492 (void) SCIPsnprintf(SCIProwprepGetName(rowprep), SCIP_MAXSTRLEN, "%sestimate_concave%p_%s%" SCIP_LONGINT_FORMAT, 2493 overestimate ? "over" : "under", 2494 (void*)expr, 2495 sol != NULL ? "sol" : "lp", 2496 sol != NULL ? (SCIP_Longint) SCIPsolGetIndex(sol) : SCIPgetNLPs(scip)); 2497 } 2498 else 2499 { 2500 SCIPfreeRowprep(scip, &rowprep); 2501 } 2502 2503 if( addbranchscores ) 2504 { 2505 SCIP_Real violation; 2506 2507 /* check how much is the violation on the side that we estimate */ 2508 if( auxvalue == SCIP_INVALID ) 2509 { 2510 /* if cannot evaluate, then always branch */ 2511 violation = SCIPinfinity(scip); 2512 } 2513 else 2514 { 2515 SCIP_Real auxval; 2516 2517 /* get value of auxiliary variable of this expression */ 2518 assert(SCIPgetExprAuxVarNonlinear(expr) != NULL); 2519 auxval = SCIPgetSolVal(scip, sol, SCIPgetExprAuxVarNonlinear(expr)); 2520 2521 /* compute the violation 2522 * if we underestimate, then we enforce expr <= auxval, so violation is (positive part of) auxvalue - auxval 2523 * if we overestimate, then we enforce expr >= auxval, so violation is (positive part of) auxval - auxvalue 2524 */ 2525 if( !overestimate ) 2526 violation = MAX(0.0, auxvalue - auxval); 2527 else 2528 violation = MAX(0.0, auxval - auxvalue); 2529 } 2530 assert(violation >= 0.0); 2531 2532 /* add violation as branching-score to expressions; the core will take care distributing this onto variables */ 2533 if( nlhdlrexprdata->nleafs == 1 ) 2534 { 2535 SCIP_EXPR* e; 2536 e = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, nlhdlrexprdata->leafexprs[0]); 2537 SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, &e, 1, violation, sol, addedbranchscores) ); 2538 } 2539 else 2540 { 2541 SCIP_EXPR** exprs; 2542 int c; 2543 2544 /* map leaf expressions back to original expressions 2545 * TODO do this once at end of detect and store in nlhdlrexprdata 2546 */ 2547 SCIP_CALL( SCIPallocBufferArray(scip, &exprs, nlhdlrexprdata->nleafs) ); 2548 for( c = 0; c < nlhdlrexprdata->nleafs; ++c ) 2549 exprs[c] = (SCIP_EXPR*)SCIPhashmapGetImage(nlhdlrexprdata->nlexpr2origexpr, nlhdlrexprdata->leafexprs[c]); 2550 2551 SCIP_CALL( SCIPaddExprsViolScoreNonlinear(scip, exprs, nlhdlrexprdata->nleafs, violation, sol, addedbranchscores) ); 2552 2553 SCIPfreeBufferArray(scip, &exprs); 2554 } 2555 } 2556 2557 return SCIP_OKAY; 2558 } 2559 2560 /** includes nonlinear handler in another scip instance */ 2561 static 2562 SCIP_DECL_NLHDLRCOPYHDLR(nlhdlrCopyhdlrConcave) 2563 { /*lint --e{715}*/ 2564 assert(targetscip != NULL); 2565 assert(sourcenlhdlr != NULL); 2566 assert(strcmp(SCIPnlhdlrGetName(sourcenlhdlr), CONCAVE_NLHDLR_NAME) == 0); 2567 2568 SCIP_CALL( SCIPincludeNlhdlrConcave(targetscip) ); 2569 2570 return SCIP_OKAY; 2571 } 2572 2573 /** includes concave nonlinear handler in nonlinear constraint handler */ 2574 SCIP_RETCODE SCIPincludeNlhdlrConcave( 2575 SCIP* scip /**< SCIP data structure */ 2576 ) 2577 { 2578 SCIP_NLHDLR* nlhdlr; 2579 SCIP_NLHDLRDATA* nlhdlrdata; 2580 2581 assert(scip != NULL); 2582 2583 SCIP_CALL( SCIPallocBlockMemory(scip, &nlhdlrdata) ); 2584 nlhdlrdata->isnlhdlrconvex = FALSE; 2585 nlhdlrdata->evalsol = NULL; 2586 nlhdlrdata->randnumgen = NULL; 2587 2588 SCIP_CALL( SCIPincludeNlhdlrNonlinear(scip, &nlhdlr, CONCAVE_NLHDLR_NAME, CONCAVE_NLHDLR_DESC, 2589 CONCAVE_NLHDLR_DETECTPRIORITY, CONCAVE_NLHDLR_ENFOPRIORITY, nlhdlrDetectConcave, nlhdlrEvalAuxConvexConcave, nlhdlrdata) ); 2590 assert(nlhdlr != NULL); 2591 2592 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/detectsum", 2593 "whether to run convexity detection when the root of an expression is a sum", 2594 &nlhdlrdata->detectsum, FALSE, DEFAULT_DETECTSUM, NULL, NULL) ); 2595 2596 /* "extended" formulations of a concave expressions can give worse estimators */ 2597 nlhdlrdata->extendedform = FALSE; 2598 2599 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxquadratic", 2600 "whether to use convexity check on quadratics", 2601 &nlhdlrdata->cvxquadratic, TRUE, DEFAULT_CVXQUADRATIC_CONCAVE, NULL, NULL) ); 2602 2603 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxsignomial", 2604 "whether to use convexity check on signomials", 2605 &nlhdlrdata->cvxsignomial, TRUE, DEFAULT_CVXSIGNOMIAL, NULL, NULL) ); 2606 2607 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/cvxprodcomp", 2608 "whether to use convexity check on product composition f(h)*h", 2609 &nlhdlrdata->cvxprodcomp, TRUE, DEFAULT_CVXPRODCOMP, NULL, NULL) ); 2610 2611 SCIP_CALL( SCIPaddBoolParam(scip, "nlhdlr/" CONCAVE_NLHDLR_NAME "/handletrivial", 2612 "whether to also handle trivial convex expressions", 2613 &nlhdlrdata->handletrivial, TRUE, DEFAULT_HANDLETRIVIAL, NULL, NULL) ); 2614 2615 SCIPnlhdlrSetFreeHdlrData(nlhdlr, nlhdlrfreeHdlrDataConvexConcave); 2616 SCIPnlhdlrSetCopyHdlr(nlhdlr, nlhdlrCopyhdlrConcave); 2617 SCIPnlhdlrSetFreeExprData(nlhdlr, nlhdlrfreeExprDataConvexConcave); 2618 SCIPnlhdlrSetSepa(nlhdlr, nlhdlrInitSepaConcave, NULL, nlhdlrEstimateConcave, NULL); 2619 SCIPnlhdlrSetInitExit(nlhdlr, NULL, nlhdlrExitConcave); 2620 2621 return SCIP_OKAY; 2622 } 2623 2624 /** checks whether a given expression is convex or concave w.r.t. the original variables 2625 * 2626 * This function uses the methods that are used in the detection algorithm of the convex nonlinear handler. 2627 */ 2628 SCIP_RETCODE SCIPhasExprCurvature( 2629 SCIP* scip, /**< SCIP data structure */ 2630 SCIP_EXPR* expr, /**< expression */ 2631 SCIP_EXPRCURV curv, /**< curvature to check for */ 2632 SCIP_Bool* success, /**< buffer to store whether expression has curvature curv (w.r.t. original variables) */ 2633 SCIP_HASHMAP* assumevarfixed /**< hashmap containing variables that should be assumed to be fixed, or NULL */ 2634 ) 2635 { 2636 SCIP_NLHDLRDATA nlhdlrdata; 2637 SCIP_EXPR* rootnlexpr; 2638 SCIP_HASHMAP* nlexpr2origexpr; 2639 int nleafs; 2640 2641 assert(expr != NULL); 2642 assert(curv != SCIP_EXPRCURV_UNKNOWN); 2643 assert(success != NULL); 2644 2645 /* create temporary hashmap */ 2646 SCIP_CALL( SCIPhashmapCreate(&nlexpr2origexpr, SCIPblkmem(scip), 20) ); 2647 2648 /* prepare nonlinear handler data */ 2649 nlhdlrdata.isnlhdlrconvex = TRUE; 2650 nlhdlrdata.evalsol = NULL; 2651 nlhdlrdata.detectsum = TRUE; 2652 nlhdlrdata.extendedform = FALSE; 2653 nlhdlrdata.cvxquadratic = TRUE; 2654 nlhdlrdata.cvxsignomial = TRUE; 2655 nlhdlrdata.cvxprodcomp = TRUE; 2656 nlhdlrdata.handletrivial = TRUE; 2657 2658 SCIP_CALL( constructExpr(scip, &nlhdlrdata, &rootnlexpr, nlexpr2origexpr, &nleafs, expr, curv, assumevarfixed, FALSE, success) ); 2659 2660 /* free created expression */ 2661 if( rootnlexpr != NULL ) 2662 { 2663 SCIP_CALL( SCIPreleaseExpr(scip, &rootnlexpr) ); 2664 } 2665 2666 /* free hashmap */ 2667 SCIPhashmapFree(&nlexpr2origexpr); 2668 2669 return SCIP_OKAY; 2670 } 2671