1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file expr_product.c 26 * @ingroup DEFPLUGINS_EXPR 27 * @brief product expression handler 28 * @author Stefan Vigerske 29 * @author Benjamin Mueller 30 * @author Felipe Serrano 31 * @author Ksenia Bestuzheva 32 */ 33 34 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 35 36 #include <string.h> 37 38 #include "scip/pub_expr.h" 39 #include "scip/expr_product.h" 40 #include "scip/expr_sum.h" 41 #include "scip/expr_pow.h" 42 #include "scip/expr_value.h" 43 #include "scip/expr_exp.h" 44 #include "scip/expr_abs.h" 45 #include "scip/expr_entropy.h" 46 #include "scip/cons_nonlinear.h" 47 #include "scip/pub_misc.h" 48 #include "scip/nlhdlr_bilinear.h" 49 #include "symmetry/struct_symmetry.h" 50 51 #define EXPRHDLR_NAME "prod" 52 #define EXPRHDLR_DESC "product expression" 53 #define EXPRHDLR_PRECEDENCE 50000 54 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(54949.0) 55 56 /** macro to activate/deactivate debugging information of simplify method */ 57 /*lint -emacro(681,debugSimplify) */ 58 /*lint -emacro(506,debugSimplify) */ 59 /*lint -emacro(774,debugSimplify) */ 60 #ifdef SIMPLIFY_DEBUG 61 #define debugSimplify printf 62 #else 63 #define debugSimplify while( FALSE ) printf 64 #endif 65 66 67 /*lint -e777*/ 68 69 /* 70 * Data structures 71 */ 72 73 /** expression data */ 74 struct SCIP_ExprData 75 { 76 SCIP_Real coefficient; /**< coefficient */ 77 }; 78 79 struct SCIP_ExprhdlrData 80 { 81 SCIP_CONSHDLR* conshdlr; /**< nonlinear constraint handler (to compute estimates for > 2-dim products) */ 82 83 SCIP_Bool expandalways; /**< whether to expand products of a sum and several factors in simplify (SP12b) */ 84 }; 85 86 /** node for linked list of expressions */ 87 struct exprnode 88 { 89 SCIP_EXPR* expr; /**< expression in node */ 90 struct exprnode* next; /**< next node */ 91 }; 92 93 typedef struct exprnode EXPRNODE; 94 95 /* 96 * Local methods 97 */ 98 99 /** evaluation callback for (vertex-polyhedral) functions used as input for facet computation of its envelopes */ 100 static 101 SCIP_DECL_VERTEXPOLYFUN(prodfunction) 102 { 103 /* funcdata is a pointer to the double holding the coefficient */ 104 SCIP_Real ret = *(SCIP_Real*)funcdata; 105 int i; 106 107 for( i = 0; i < nargs; ++i ) 108 ret *= args[i]; 109 110 return ret; 111 } 112 113 static 114 SCIP_RETCODE buildSimplifiedProduct( 115 SCIP* scip, /**< SCIP data structure */ 116 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * PI simplifiedfactors */ 117 EXPRNODE** simplifiedfactors, /**< factors of simplified product */ 118 SCIP_Bool expandalways, /**< whether to expand products of a sum and several factors in simplify (SP12b) */ 119 SCIP_Bool changed, /**< indicates whether some of the simplified factors was changed */ 120 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 121 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 122 void* ownercreatedata /**< data to pass to ownercreate */ 123 ); 124 125 /* methods for handling linked list of expressions */ 126 127 /** inserts newnode at beginning of list */ 128 static 129 void insertFirstList( 130 EXPRNODE* newnode, /**< node to insert */ 131 EXPRNODE** list /**< list */ 132 ) 133 { 134 assert(list != NULL); 135 assert(newnode != NULL); 136 137 newnode->next = *list; 138 *list = newnode; 139 } 140 141 /** removes first element of list and returns it */ 142 static 143 EXPRNODE* listPopFirst( 144 EXPRNODE** list /**< list */ 145 ) 146 { 147 EXPRNODE* first; 148 149 assert(list != NULL); 150 151 if( *list == NULL ) 152 return NULL; 153 154 first = *list; 155 *list = (*list)->next; 156 first->next = NULL; 157 158 return first; 159 } 160 161 /** returns length of list */ 162 static 163 int listLength( 164 EXPRNODE* list /**< list */ 165 ) 166 { 167 int length; 168 169 if( list == NULL ) 170 return 0; 171 172 length = 1; 173 while( (list=list->next) != NULL ) 174 ++length; 175 176 return length; 177 } 178 179 /** creates expression node and captures expression */ 180 static 181 SCIP_RETCODE createExprNode( 182 SCIP* scip, /**< SCIP data structure */ 183 SCIP_EXPR* expr, /**< expression stored at node */ 184 EXPRNODE** newnode /**< pointer to store node */ 185 ) 186 { 187 SCIP_CALL( SCIPallocBlockMemory(scip, newnode) ); 188 189 (*newnode)->expr = expr; 190 (*newnode)->next = NULL; 191 SCIPcaptureExpr(expr); 192 193 return SCIP_OKAY; 194 } 195 196 /** creates expression list from expressions */ 197 static 198 SCIP_RETCODE createExprlistFromExprs( 199 SCIP* scip, /**< SCIP data structure */ 200 SCIP_EXPR** exprs, /**< expressions stored in list */ 201 int nexprs, /**< number of expressions */ 202 EXPRNODE** list /**< pointer to store list */ 203 ) 204 { 205 int i; 206 207 assert(*list == NULL); 208 assert(nexprs > 0); 209 210 debugSimplify("building expr list from %d expressions\n", nexprs); 211 for( i = nexprs - 1; i >= 0; --i ) 212 { 213 EXPRNODE* newnode; 214 215 SCIP_CALL( createExprNode(scip, exprs[i], &newnode) ); 216 insertFirstList(newnode, list); 217 } 218 assert(nexprs > 1 || (*list)->next == NULL); 219 220 return SCIP_OKAY; 221 } 222 223 /** frees expression node and releases expressions */ 224 static 225 SCIP_RETCODE freeExprNode( 226 SCIP* scip, /**< SCIP data structure */ 227 EXPRNODE** node /**< node to be freed */ 228 ) 229 { 230 assert(node != NULL && *node != NULL); 231 232 SCIP_CALL( SCIPreleaseExpr(scip, &(*node)->expr) ); 233 SCIPfreeBlockMemory(scip, node); 234 235 return SCIP_OKAY; 236 } 237 238 /** frees an expression list */ 239 static 240 SCIP_RETCODE freeExprlist( 241 SCIP* scip, /**< SCIP data structure */ 242 EXPRNODE** exprlist /**< list */ 243 ) 244 { 245 EXPRNODE* current; 246 247 if( *exprlist == NULL ) 248 return SCIP_OKAY; 249 250 current = *exprlist; 251 while( current != NULL ) 252 { 253 EXPRNODE* tofree; 254 255 tofree = current; 256 current = current->next; 257 SCIP_CALL( freeExprNode(scip, &tofree) ); 258 } 259 assert(current == NULL); 260 *exprlist = NULL; 261 262 return SCIP_OKAY; 263 } 264 265 /* helper functions for simplifying expressions */ 266 267 /** creates a product expression with the elements of exprlist as its children */ 268 static 269 SCIP_RETCODE createExprProductFromExprlist( 270 SCIP* scip, /**< SCIP data structure */ 271 EXPRNODE* exprlist, /**< list containing the children of expr */ 272 SCIP_Real coef, /**< coef of expr */ 273 SCIP_EXPR** expr, /**< pointer to store the product expression */ 274 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 275 void* ownercreatedata /**< data to pass to ownercreate */ 276 ) 277 { 278 int i; 279 int nchildren; 280 SCIP_EXPR** children; 281 282 /* asserts SP8 */ 283 assert(coef == 1.0); 284 nchildren = listLength(exprlist); 285 286 SCIP_CALL( SCIPallocBufferArray(scip, &children, nchildren) ); 287 288 for( i = 0; i < nchildren; ++i ) 289 { 290 children[i] = exprlist->expr; 291 exprlist = exprlist->next; 292 } 293 294 assert(exprlist == NULL); 295 296 SCIP_CALL( SCIPcreateExprProduct(scip, expr, nchildren, children, coef, ownercreate, ownercreatedata) ); 297 298 SCIPfreeBufferArray(scip, &children); 299 300 return SCIP_OKAY; 301 } 302 303 /** simplifies a factor of a product expression: base, so that it is a valid children of a simplified product expr 304 * 305 * @note In contrast to other simplify methods, this does *not* return a simplified expression. 306 * Instead, the method is intended to be called only when simplifying a product expression. 307 * Since in general, base is not a simplified child of a product expression, this method returns 308 * a list of expressions L, such that (prod L) = baset *and* each expression in L 309 * is a valid child of a simplified product expression. 310 */ 311 static 312 SCIP_RETCODE simplifyFactor( 313 SCIP* scip, /**< SCIP data structure */ 314 SCIP_EXPR* factor, /**< expression to be simplified */ 315 SCIP_Real* simplifiedcoef, /**< coefficient of parent product expression */ 316 EXPRNODE** simplifiedfactor, /**< pointer to store the resulting expression node/list of nodes */ 317 SCIP_Bool* changed /**< pointer to store if some term actually got simplified */ 318 ) 319 { 320 assert(simplifiedfactor != NULL); 321 assert(*simplifiedfactor == NULL); 322 assert(factor != NULL); 323 assert(changed != NULL); 324 325 /* enforces SP7 */ 326 if( SCIPisExprValue(scip, factor) ) 327 { 328 *changed = TRUE; 329 *simplifiedcoef *= SCIPgetValueExprValue(factor); 330 return SCIP_OKAY; 331 } 332 333 /* enforces SP2 */ 334 if( SCIPisExprProduct(scip, factor) ) 335 { 336 *changed = TRUE; 337 338 /* assert SP8 */ 339 assert(SCIPgetCoefExprProduct(factor) == 1.0); 340 debugSimplify("[simplifyFactor] seeing a product: include its children\n"); 341 342 SCIP_CALL( createExprlistFromExprs(scip, SCIPexprGetChildren(factor), SCIPexprGetNChildren(factor), simplifiedfactor) ); 343 344 return SCIP_OKAY; 345 } 346 347 /* enforces SP13: a sum with a unique child and no constant -> take the coefficient and use its child as factor */ 348 if( SCIPisExprSum(scip, factor) && SCIPexprGetNChildren(factor) == 1 && SCIPgetConstantExprSum(factor) == 0.0 ) 349 { 350 *changed = TRUE; 351 352 /* assert SS8 and SS7 */ 353 assert(SCIPgetCoefsExprSum(factor)[0] != 0.0 && SCIPgetCoefsExprSum(factor)[0] != 1.0); 354 debugSimplify("[simplifyFactor] seeing a sum of the form coef * child : take coef and child apart\n"); 355 356 if( SCIPisExprProduct(scip, SCIPexprGetChildren(factor)[0]) ) 357 { 358 /* if child is a product, then add its children to exprlist */ 359 SCIP_CALL( createExprlistFromExprs(scip, SCIPexprGetChildren(SCIPexprGetChildren(factor)[0]), SCIPexprGetNChildren(SCIPexprGetChildren(factor)[0]), simplifiedfactor) ); 360 *simplifiedcoef *= SCIPgetCoefExprProduct(SCIPexprGetChildren(factor)[0]); 361 } 362 else 363 { 364 SCIP_CALL( createExprlistFromExprs(scip, SCIPexprGetChildren(factor), 1, simplifiedfactor) ); 365 } 366 *simplifiedcoef *= SCIPgetCoefsExprSum(factor)[0]; 367 368 return SCIP_OKAY; 369 } 370 371 /* the given (simplified) expression `factor`, can be a child of a simplified product */ 372 assert(!SCIPisExprProduct(scip, factor)); 373 assert(!SCIPisExprValue(scip, factor)); 374 375 SCIP_CALL( createExprNode(scip, factor, simplifiedfactor) ); 376 377 return SCIP_OKAY; 378 } 379 380 /** merges tomerge into finalchildren 381 * 382 * Both, tomerge and finalchildren contain expressions that could be the children of a simplified product 383 * (except for SP8 and SP10 which are enforced later). 384 * However, the concatenation of both lists will not in general yield a simplified product expression, 385 * because SP4, SP5 and SP14 could be violated. So the purpose of this method is to enforce SP4, SP5 and SP14. 386 * In the process of enforcing SP4, it could happen that SP2 is violated. Since enforcing SP2 387 * could generate further violations, we remove the affected children from finalchildren 388 * and include them in unsimplifiedchildren for further processing. 389 * @note if tomerge has more than one element, then they are the children of a simplified product expression 390 */ 391 static 392 SCIP_RETCODE mergeProductExprlist( 393 SCIP* scip, /**< SCIP data structure */ 394 EXPRNODE* tomerge, /**< list to merge */ 395 EXPRNODE** finalchildren, /**< pointer to store the result of merge between tomerge and *finalchildren */ 396 EXPRNODE** unsimplifiedchildren,/**< the list of children that should go to the product expression; 397 * they are unsimplified when seen as children of a simplified product */ 398 SCIP_Bool* changed, /**< pointer to store if some term actually got simplified */ 399 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 400 void* ownercreatedata /**< data to pass to ownercreate */ 401 ) 402 { 403 EXPRNODE* tomergenode; 404 EXPRNODE* current; 405 EXPRNODE* previous; 406 407 if( tomerge == NULL ) 408 return SCIP_OKAY; 409 410 if( *finalchildren == NULL ) 411 { 412 *finalchildren = tomerge; 413 return SCIP_OKAY; 414 } 415 416 tomergenode = tomerge; 417 current = *finalchildren; 418 previous = NULL; 419 420 while( tomergenode != NULL && current != NULL ) 421 { 422 int compareres; 423 EXPRNODE* aux; 424 SCIP_EXPR* base1; 425 SCIP_EXPR* base2; 426 SCIP_Real expo1; 427 SCIP_Real expo2; 428 SCIP_Bool issignpower1; 429 SCIP_Bool issignpower2; 430 431 /* assert invariants */ 432 assert(!SCIPisExprValue(scip, tomergenode->expr)); 433 assert(!SCIPisExprValue(scip, current->expr)); 434 assert(previous == NULL || previous->next == current); 435 436 /* we are going to multiply the two exprs: current and tomergenode 437 * we first check if they are both exponentials 438 * if so, we multiply them 439 * otherwise, we interpret them as base^exponent 440 * the base of an expr is the expr itself 441 * if type(expr) != pow, otherwise it is the child of pow 442 */ 443 444 /* if both are exponentials, create a new exponential with the sum of their children */ 445 if( SCIPisExprExp(scip, current->expr) && SCIPisExprExp(scip, tomergenode->expr) ) 446 { 447 SCIP_EXPR* sum; 448 SCIP_EXPR* simplifiedsum; 449 SCIP_EXPR* expexpr; 450 SCIP_EXPR* simplifiedexp; 451 452 /* inform that expressions changed */ 453 *changed = TRUE; 454 455 /* create sum */ 456 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(current->expr), NULL, 0.0, ownercreate, ownercreatedata) ); 457 SCIP_CALL( SCIPappendExprSumExpr(scip, sum, SCIPexprGetChildren(tomergenode->expr)[0], 1.0) ); 458 459 /* simplify sum */ 460 SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) ); 461 SCIP_CALL( SCIPreleaseExpr(scip, &sum) ); 462 463 /* create exponential */ 464 SCIP_CALL( SCIPcreateExprExp(scip, &expexpr, simplifiedsum, ownercreate, ownercreatedata) ); 465 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) ); 466 467 /* simplify exponential */ 468 SCIP_CALL( SCIPcallExprSimplify(scip, expexpr, &simplifiedexp, ownercreate, ownercreatedata) ); 469 SCIP_CALL( SCIPreleaseExpr(scip, &expexpr) ); 470 471 /* note that simplified exponential might be a product exp(x) * exp(-x + log(y*z)) -> y*z and so it is not a 472 * valid child of a simplified product; therefore we add it to the unsimplifiedchildren's list 473 */ 474 475 /* replace tomergenode's expression with simplifiedexp */ 476 /* TODO: this code repeats below; add new function to avoid duplication */ 477 SCIP_CALL( SCIPreleaseExpr(scip, &tomergenode->expr) ); 478 tomergenode->expr = simplifiedexp; 479 480 /* move tomergenode to unsimplifiedchildren */ 481 aux = tomergenode; 482 tomergenode = tomergenode->next; 483 insertFirstList(aux, unsimplifiedchildren); 484 485 /* remove current */ 486 if( current == *finalchildren ) 487 { 488 assert(previous == NULL); 489 aux = listPopFirst(finalchildren); 490 assert(aux == current); 491 current = *finalchildren; 492 } 493 else 494 { 495 assert(previous != NULL); 496 aux = current; 497 current = current->next; 498 previous->next = current; 499 } 500 SCIP_CALL( freeExprNode(scip, &aux) ); 501 502 continue; 503 } 504 505 /* they were not exponentials, so collect bases and exponents */ 506 if( SCIPisExprPower(scip, current->expr) ) 507 { 508 base1 = SCIPexprGetChildren(current->expr)[0]; 509 expo1 = SCIPgetExponentExprPow(current->expr); 510 issignpower1 = FALSE; 511 } 512 else if( SCIPisExprSignpower(scip, current->expr) ) 513 { 514 base1 = SCIPexprGetChildren(current->expr)[0]; 515 expo1 = SCIPgetExponentExprPow(current->expr); 516 issignpower1 = TRUE; 517 } 518 else 519 { 520 base1 = current->expr; 521 expo1 = 1.0; 522 issignpower1 = FALSE; 523 } 524 if( SCIPisExprPower(scip, tomergenode->expr) ) 525 { 526 base2 = SCIPexprGetChildren(tomergenode->expr)[0]; 527 expo2 = SCIPgetExponentExprPow(tomergenode->expr); 528 issignpower2 = FALSE; 529 } 530 else if( SCIPisExprSignpower(scip, tomergenode->expr) ) 531 { 532 base2 = SCIPexprGetChildren(tomergenode->expr)[0]; 533 expo2 = SCIPgetExponentExprPow(tomergenode->expr); 534 issignpower2 = TRUE; 535 } 536 else 537 { 538 base2 = tomergenode->expr; 539 expo2 = 1.0; 540 issignpower2 = FALSE; 541 } 542 543 if( SCIPcompareExpr(scip, base1, base2) == 0 ) 544 { 545 /* the bases are the same, so we should try to merge the multiplication of the powers */ 546 SCIP_EXPR* power = NULL; 547 548 if( !issignpower1 && !issignpower2 ) 549 { 550 /* and both are normal power, then add to unsimplifiedchildren the resulting expr of simplify(base^(expo1 + expo2)) */ 551 #if 0 /* TODO we should not loose the implicit base >= 0 constraint, if there is one, but then we should look at bounds on base; simplify currently doesn't */ 552 /* 553 * unless expo1 or expo2 are fractional but expo1+expo2 is not fractional, then we better keep the original 554 * the reason for that is that x^fractional implies a constraint x >= 0 555 */ 556 if( (EPSISINT(expo1, 0.0) && EPSISINT(expo2, 0.0)) || !EPSISINT(expo1+expo2, 0.0) ) /*lint !e835*/ 557 #endif 558 { 559 SCIP_CALL( SCIPcreateExprPow(scip, &power, base1, expo1 + expo2, ownercreate, ownercreatedata) ); 560 } 561 } 562 else if( issignpower1 ^ issignpower2 ) 563 { 564 /* exactly one is signpower: sign(x) |x|^expo1 x^expo2 = sign(x)^(1+expo2) |x|^(expo1+expo2), with x = base */ 565 if( EPSISINT(expo2, 0.0) ) /*lint !e835*/ 566 { 567 if( (int)expo2 % 2 == 0 ) 568 { 569 /* if expo2 is even, then sign(x)^(1+expo2) = sign(x), so we have signpower: sign(x) |x|^(expo1+expo2) 570 * TODO: we can remove this case distinction once the simplification of power expressions tranform 571 * |expr|^even -> expr^even, since the call to SCIPcallExprSimplify(scip, conshdlr, power, 572 * &simplifiedpower) below will take care of this. 573 */ 574 SCIP_CALL( SCIPcreateExprSignpower(scip, &power, base1, expo1 + expo2, ownercreate, ownercreatedata) ); 575 } 576 else 577 { 578 /* if expo2 is odd, then sign(x)^(1+expo2) = 1, so we have |x|^(expo1+expo2) */ 579 SCIP_EXPR* absbase; 580 581 SCIP_CALL( SCIPcreateExprAbs(scip, &absbase, base1, ownercreate, ownercreatedata) ); 582 SCIP_CALL( SCIPcreateExprPow(scip, &power, absbase, expo1 + expo2, ownercreate, ownercreatedata) ); 583 SCIP_CALL( SCIPreleaseExpr(scip, &absbase) ); 584 } 585 } 586 else if( !EPSISINT(expo1+expo2, 0.0) ) /*lint !e835*/ 587 { 588 /* if expo2 is fractional and expo1+expo2 is fractional, then we need x >= 0, so we can use x^(expo1+expo2) */ 589 SCIP_CALL( SCIPcreateExprPow(scip, &power, base1, expo1 + expo2, ownercreate, ownercreatedata) ); 590 } 591 /* else: expo2 is fractional but expo1+expo2 is integral, then we better do not do anything for now 592 * (leave power at NULL) 593 */ 594 } 595 else 596 { 597 /* if both are signpower, then we have |base|^(expo1+expo2) 598 * if expo1+expo2 is even, then we can change this to base^(expo1+expo2) 599 */ 600 if( EPSISINT(expo1+expo2, 0.0) && (int)(expo1+expo2)%2 == 0 ) /*lint !e835*/ 601 { 602 SCIP_CALL( SCIPcreateExprPow(scip, &power, base1, expo1 + expo2, ownercreate, ownercreatedata) ); 603 } 604 else 605 { 606 SCIP_EXPR* absbase; 607 608 SCIP_CALL( SCIPcreateExprAbs(scip, &absbase, base1, ownercreate, ownercreatedata) ); 609 SCIP_CALL( SCIPcreateExprPow(scip, &power, absbase, expo1 + expo2, ownercreate, ownercreatedata) ); 610 SCIP_CALL( SCIPreleaseExpr(scip, &absbase) ); 611 } 612 } 613 614 if( power != NULL ) 615 { 616 /* we have created a new power: simplify again and continue */ 617 SCIP_EXPR* simplifiedpower; 618 619 /* call simplifyPow or simplifySignpower */ 620 SCIP_CALL( SCIPcallExprSimplify(scip, power, &simplifiedpower, ownercreate, ownercreatedata) ); 621 SCIP_CALL( SCIPreleaseExpr(scip, &power) ); 622 623 /* replace tomergenode's expression with simplifiedpower */ 624 SCIP_CALL( SCIPreleaseExpr(scip, &tomergenode->expr) ); 625 tomergenode->expr = simplifiedpower; 626 627 *changed = TRUE; 628 629 /* move tomergenode to unsimplifiedchildren */ 630 aux = tomergenode; 631 tomergenode = tomergenode->next; 632 insertFirstList(aux, unsimplifiedchildren); 633 634 /* remove current */ 635 if( current == *finalchildren ) 636 { 637 assert(previous == NULL); 638 aux = listPopFirst(finalchildren); 639 assert(aux == current); 640 current = *finalchildren; 641 } 642 else 643 { 644 assert(previous != NULL); 645 aux = current; 646 current = current->next; 647 previous->next = current; 648 } 649 SCIP_CALL( freeExprNode(scip, &aux) ); 650 651 continue; 652 } 653 } 654 655 /* bases are not the same, or we do not want to merge them 656 * then expressions cannot be the same 657 * therefore we need to insert tomergenode in finalchildren 658 * for this, we need to take care of the order 659 */ 660 compareres = SCIPcompareExpr(scip, current->expr, tomergenode->expr); 661 if( compareres == -1 ) 662 { 663 /* current < tomergenode => move current */ 664 previous = current; 665 current = current->next; 666 } 667 else 668 { 669 *changed = TRUE; 670 assert(compareres == 1); 671 672 /* insert: if current is the first node, then insert at beginning; otherwise, insert between previous and current */ 673 if( current == *finalchildren ) 674 { 675 assert(previous == NULL); 676 aux = tomergenode; 677 tomergenode = tomergenode->next; 678 insertFirstList(aux, finalchildren); 679 previous = *finalchildren; 680 } 681 else 682 { 683 assert(previous != NULL); 684 /* extract */ 685 aux = tomergenode; 686 tomergenode = tomergenode->next; 687 /* insert */ 688 previous->next = aux; 689 aux->next = current; 690 previous = aux; 691 } 692 } 693 } 694 695 /* if all nodes of tomerge were merged, we are done */ 696 if( tomergenode == NULL ) 697 return SCIP_OKAY; 698 699 assert(current == NULL); 700 701 /* if all nodes of finalchildren were cancelled by nodes of tomerge (i.e., transfered to unsimplifiedchildren), 702 * then the rest of tomerge is finalchildren 703 */ 704 if( *finalchildren == NULL ) 705 { 706 assert(previous == NULL); 707 *finalchildren = tomergenode; 708 return SCIP_OKAY; 709 } 710 711 /* there are still nodes of tomerge unmerged; these nodes are larger than finalchildren, so append at end */ 712 assert(previous != NULL && previous->next == NULL); 713 previous->next = tomergenode; 714 715 return SCIP_OKAY; 716 } 717 718 /** simplifies the given (simplified) exprs so that they can be factors of a simplified product 719 * 720 * in particular, it will sort and multiply factors whose product leads to new expressions 721 */ 722 static 723 SCIP_RETCODE simplifyMultiplyChildren( 724 SCIP* scip, /**< SCIP data structure */ 725 SCIP_EXPR** exprs, /**< factors to be simplified */ 726 int nexprs, /**< number of factors */ 727 SCIP_Real* simplifiedcoef, /**< buffer to store coefficient of PI exprs; needs to be initialized */ 728 EXPRNODE** finalchildren, /**< expr node list to store the simplified factors */ 729 SCIP_Bool* changed, /**< buffer to store whether some factor changed */ 730 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 731 void* ownercreatedata /**< data to pass to ownercreate */ 732 ) 733 { 734 EXPRNODE* unsimplifiedchildren; 735 736 /* set up list of current children (when looking at each of them individually, they are simplified, but as 737 * children of a product expression they might be unsimplified) 738 */ 739 unsimplifiedchildren = NULL; 740 SCIP_CALL( createExprlistFromExprs(scip, exprs, nexprs, &unsimplifiedchildren) ); 741 742 *changed = FALSE; 743 744 /* while there are still children to process */ 745 *finalchildren = NULL; 746 while( unsimplifiedchildren != NULL ) 747 { 748 EXPRNODE* tomerge; 749 EXPRNODE* first; 750 751 first = listPopFirst(&unsimplifiedchildren); 752 assert(first != NULL); 753 754 #ifdef SIMPLIFY_DEBUG 755 debugSimplify("simplifying factor:\n"); 756 SCIP_CALL( SCIPprintExpr(scip, first->expr, NULL) ); 757 SCIPinfoMessage(scip, NULL, "\n"); 758 #endif 759 760 /* enforces SP2, SP7 and SP13 */ 761 tomerge = NULL; 762 SCIP_CALL( simplifyFactor(scip, first->expr, simplifiedcoef, &tomerge, changed) ); 763 764 /* enforces SP4 and SP5 note: merge frees (or uses) the nodes of the tomerge list */ 765 SCIP_CALL( mergeProductExprlist(scip, tomerge, finalchildren, &unsimplifiedchildren, changed, ownercreate, ownercreatedata) ); 766 767 /* free first */ 768 SCIP_CALL( freeExprlist(scip, &first) ); 769 770 /* if the simplified coefficient is 0, we can return value 0 */ 771 if( *simplifiedcoef == 0.0 ) 772 { 773 *changed = TRUE; 774 SCIP_CALL( freeExprlist(scip, finalchildren) ); 775 SCIP_CALL( freeExprlist(scip, &unsimplifiedchildren) ); 776 assert(*finalchildren == NULL); 777 break; 778 } 779 } 780 return SCIP_OKAY; 781 } 782 783 /* make sure product has at least two children 784 * - if it is empty; return value 785 * - if it has one child and coef = 1; return child 786 * - if it has one child and coef != 1; return (sum 0 coef expr) 787 */ 788 static 789 SCIP_RETCODE enforceSP10( 790 SCIP* scip, /**< SCIP data structure */ 791 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * PI simplifiedfactors */ 792 EXPRNODE* finalchildren, /**< factors of simplified product */ 793 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 794 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 795 void* ownercreatedata /**< data to pass to ownercreate */ 796 ) 797 { 798 /* empty list --> return value */ 799 if( finalchildren == NULL ) 800 { 801 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, simplifiedcoef, ownercreate, ownercreatedata) ); 802 return SCIP_OKAY; 803 } 804 805 /* one child and coef equal to 1 --> return child */ 806 if( finalchildren->next == NULL && simplifiedcoef == 1.0 ) 807 { 808 *simplifiedexpr = finalchildren->expr; 809 SCIPcaptureExpr(*simplifiedexpr); 810 return SCIP_OKAY; 811 } 812 813 /* one child and coef different from 1 --> return (sum 0 coef child) */ 814 if( finalchildren->next == NULL ) 815 { 816 SCIP_EXPR* sum; 817 818 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, &(finalchildren->expr), &simplifiedcoef, 0.0, ownercreate, ownercreatedata) ); 819 820 /* simplifying here is necessary, the product could have sums as children e.g., (prod 2 (sum 1 <x>)) 821 * -> (sum 0 2 (sum 1 <x>)) and that needs to be simplified to (sum 0 2 <x>) 822 */ 823 SCIP_CALL( SCIPcallExprSimplify(scip, sum, simplifiedexpr, ownercreate, ownercreatedata) ); 824 SCIP_CALL( SCIPreleaseExpr(scip, &sum) ); 825 return SCIP_OKAY; 826 } 827 828 return SCIP_OKAY; 829 } 830 831 /** checks if it is entropy expression */ 832 static 833 SCIP_RETCODE enforceSP11( 834 SCIP* scip, /**< SCIP data structure */ 835 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * PI simplifiedfactors */ 836 EXPRNODE* finalchildren, /**< factors of simplified product */ 837 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 838 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 839 void* ownercreatedata /**< data to pass to ownercreate */ 840 ) 841 { 842 SCIP_EXPR* entropicchild = NULL; 843 844 if( !(finalchildren != NULL && finalchildren->next != NULL && finalchildren->next->next == NULL) ) 845 return SCIP_OKAY; 846 847 /* could be log(expr) * expr, e.g., log(sin(x)) * sin(x) (OR11) */ 848 if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(finalchildren->expr)), "log") == 0 ) 849 { 850 assert(SCIPexprGetNChildren(finalchildren->expr) == 1); 851 if( 0 == SCIPcompareExpr(scip, SCIPexprGetChildren(finalchildren->expr)[0], finalchildren->next->expr) ) 852 entropicchild = finalchildren->next->expr; 853 } 854 /* could be expr * log(expr), e.g., (1 + abs(x)) log(1 + abs(x)) (OR11) */ 855 else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(finalchildren->next->expr)), "log") == 0 ) 856 { 857 assert(SCIPexprGetNChildren(finalchildren->next->expr) == 1); 858 if( 0 == SCIPcompareExpr(scip, SCIPexprGetChildren(finalchildren->next->expr)[0], finalchildren->expr) ) 859 entropicchild = finalchildren->expr; 860 } 861 862 /* success --> replace finalchildren by entropy expression */ 863 if( entropicchild != NULL ) 864 { 865 SCIP_EXPR* entropy; 866 867 simplifiedcoef *= -1.0; 868 869 SCIP_CALL( SCIPcreateExprEntropy(scip, &entropy, entropicchild, ownercreate, ownercreatedata) ); 870 871 /* enforces SP8: if simplifiedcoef != 1.0, transform it into a sum with the (simplified) entropy as child */ 872 if( simplifiedcoef != 1.0 ) 873 { 874 SCIP_CALL( SCIPcreateExprSum(scip, simplifiedexpr, 1, &entropy, &simplifiedcoef, 0.0, ownercreate, ownercreatedata) ); 875 SCIP_CALL( SCIPreleaseExpr(scip, &entropy) ); 876 } 877 else 878 *simplifiedexpr = entropy; 879 } 880 881 return SCIP_OKAY; 882 } 883 884 /* expands product of two sums or one sum and another expression 885 * -) two sums: (prod (sum c1 s1 ... sn) (sum c2 t1 ... tm) 886 * Builds a sum representing the expansion, where all of its children are simplified, and then simplify the sum 887 * - constant != 0 --> c1 ti or c2 * sj is simplified (ti, sj are not sums, because they are children of a simplified sum) 888 * - sj * ti may be not be simplified, so put them in a product list and simplify them from there 889 * -) one sum: (prod factor (sum c s1 ... sn)) 890 * - c != 0 --> c * factor is simplified (i.e. factor is not sum!) 891 * - factor * si may be not be simplified, so put them in a product list and simplify them from there 892 */ 893 static 894 SCIP_RETCODE enforceSP12( 895 SCIP* scip, /**< SCIP data structure */ 896 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * PI simplifiedfactors */ 897 EXPRNODE* finalchildren, /**< factors of simplified product */ 898 SCIP_Bool expandalways, /**< whether to expand products of a sum and several factors in simplify (SP12b) */ 899 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 900 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 901 void* ownercreatedata /**< data to pass to ownercreate */ 902 ) 903 { 904 /* we need only two children */ 905 if( ! (finalchildren != NULL && finalchildren->next != NULL && finalchildren->next->next == NULL) ) 906 return SCIP_OKAY; 907 908 /* handle both sums case */ 909 if( SCIPisExprSum(scip, finalchildren->expr) && SCIPisExprSum(scip, finalchildren->next->expr) ) 910 { 911 SCIP_EXPR* expanded = NULL; 912 SCIP_Real c1 = SCIPgetConstantExprSum(finalchildren->expr); 913 SCIP_Real c2 = SCIPgetConstantExprSum(finalchildren->next->expr); 914 int nchildren1 = SCIPexprGetNChildren(finalchildren->expr); 915 int nchildren2 = SCIPexprGetNChildren(finalchildren->next->expr); 916 int j; 917 int k; 918 919 #ifdef SIMPLIFY_DEBUG 920 debugSimplify("Multiplying sum1 * sum2\n"); 921 debugSimplify("sum1: \n"); 922 SCIP_CALL( SCIPprintExpr(scip, finalchildren->expr, NULL) ); 923 SCIPinfoMessage(scip, NULL, "\n"); 924 debugSimplify("sum2: \n"); 925 SCIP_CALL( SCIPprintExpr(scip, finalchildren->next->expr, NULL) ); 926 SCIPinfoMessage(scip, NULL, "\n"); 927 #endif 928 SCIP_CALL( SCIPcreateExprSum(scip, &expanded, 0, NULL, NULL, c1 * c2 * simplifiedcoef, ownercreate, ownercreatedata) ); 929 930 /* multiply c1 * sum2 */ 931 if( c1 != 0.0 ) 932 { 933 int i; 934 935 for( i = 0; i < nchildren2; ++i ) 936 { 937 SCIP_EXPR* term; 938 939 term = SCIPexprGetChildren(finalchildren->next->expr)[i]; 940 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, SCIPgetCoefsExprSum(finalchildren->next->expr)[i] * c1 * simplifiedcoef) ); 941 /* we are just re-using a child here, so do not release term! */ 942 #ifdef SIMPLIFY_DEBUG 943 debugSimplify("Multiplying %f * summand2_i\n", c1); 944 debugSimplify("summand2_i: \n"); 945 SCIP_CALL( SCIPprintExpr(scip, term, NULL) ); 946 SCIPinfoMessage(scip, NULL, "\n"); 947 #endif 948 } 949 } 950 /* multiply c2 * sum1 */ 951 if( c2 != 0.0 ) 952 { 953 int i; 954 955 for( i = 0; i < nchildren1; ++i ) 956 { 957 SCIP_EXPR* term; 958 959 term = SCIPexprGetChildren(finalchildren->expr)[i]; 960 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, SCIPgetCoefsExprSum(finalchildren->expr)[i] * c2 * simplifiedcoef) ); 961 /* we are just re-using a child here, so do not release term! */ 962 #ifdef SIMPLIFY_DEBUG 963 debugSimplify("Multiplying summand1_i * %f\n", c2); 964 debugSimplify("summand1_i: \n"); 965 SCIP_CALL( SCIPprintExpr(scip, term, NULL) ); 966 SCIPinfoMessage(scip, NULL, "\n"); 967 #endif 968 } 969 } 970 /* multiply sum1 * sum2 without constants */ 971 for( j = 0; j < nchildren1; ++j ) 972 { 973 SCIP_EXPR* factors[2]; 974 SCIP_Real coef1; 975 976 coef1 = SCIPgetCoefsExprSum(finalchildren->expr)[j]; 977 factors[0] = SCIPexprGetChildren(finalchildren->expr)[j]; 978 for( k = 0; k < nchildren2; ++k ) 979 { 980 EXPRNODE* finalfactors; 981 SCIP_Real factorscoef; 982 SCIP_Real coef2; 983 SCIP_EXPR* term = NULL; 984 SCIP_Bool dummy; 985 986 coef2 = SCIPgetCoefsExprSum(finalchildren->next->expr)[k]; 987 factors[1] = SCIPexprGetChildren(finalchildren->next->expr)[k]; 988 989 #ifdef SIMPLIFY_DEBUG 990 debugSimplify("multiplying %g expr1 * %g expr2\n", coef1, coef2); 991 debugSimplify("expr1:\n"); 992 SCIP_CALL( SCIPprintExpr(scip, factors[0], NULL) ); 993 SCIPinfoMessage(scip, NULL, "\n"); 994 debugSimplify("expr2\n"); 995 SCIP_CALL( SCIPprintExpr(scip, factors[1], NULL) ); 996 SCIPinfoMessage(scip, NULL, "\n"); 997 #endif 998 999 factorscoef = coef1 * coef2; 1000 SCIP_CALL( simplifyMultiplyChildren(scip, factors, 2, &factorscoef, &finalfactors, &dummy, ownercreate, ownercreatedata) ); 1001 assert(factorscoef != 0.0); 1002 1003 #ifdef SIMPLIFY_DEBUG 1004 { 1005 EXPRNODE* node; 1006 int i; 1007 1008 debugSimplify("Building product from simplified factors\n"); 1009 node = finalfactors; 1010 i = 0; 1011 while( node != NULL ) 1012 { 1013 debugSimplify("factor %d (nuses %d):\n", i, SCIPexprGetNUses(node->expr)); 1014 SCIP_CALL( SCIPprintExpr(scip, node->expr, NULL) ); 1015 SCIPinfoMessage(scip, NULL, "\n"); 1016 node = node->next; 1017 i++; 1018 } 1019 } 1020 #endif 1021 1022 SCIP_CALL( buildSimplifiedProduct(scip, 1.0, &finalfactors, expandalways, TRUE, &term, ownercreate, ownercreatedata) ); 1023 assert(finalfactors == NULL); 1024 assert(term != NULL); 1025 1026 #ifdef SIMPLIFY_DEBUG 1027 debugSimplify("%g expr1 * %g expr2 = %g * product\n", coef1, coef2, coef1 * coef2); 1028 debugSimplify("product: (nused %d)\n", SCIPexprGetNUses(term)); 1029 SCIP_CALL( SCIPprintExpr(scip, term, NULL) ); 1030 SCIPinfoMessage(scip, NULL, "\n"); 1031 #endif 1032 1033 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, factorscoef * simplifiedcoef) ); 1034 1035 SCIP_CALL( SCIPreleaseExpr(scip, &term) ); 1036 } 1037 } 1038 1039 /* simplify the sum */ 1040 SCIP_CALL( SCIPcallExprSimplify(scip, expanded, simplifiedexpr, ownercreate, ownercreatedata) ); 1041 SCIP_CALL( SCIPreleaseExpr(scip, &expanded) ); 1042 1043 return SCIP_OKAY; 1044 } 1045 1046 /* handle one sum case */ 1047 if( SCIPisExprSum(scip, finalchildren->expr) || SCIPisExprSum(scip, finalchildren->next->expr) ) 1048 { 1049 SCIP_EXPR* expanded = NULL; 1050 SCIP_EXPR* factors[2]; 1051 SCIP_EXPR* sum = NULL; 1052 SCIP_Real constant; 1053 int nchildren; 1054 int j; 1055 1056 if( SCIPisExprSum(scip, finalchildren->expr) ) 1057 { 1058 assert(!SCIPisExprSum(scip, finalchildren->next->expr)); 1059 sum = finalchildren->expr; 1060 factors[0] = finalchildren->next->expr; 1061 } 1062 else 1063 { 1064 assert(!SCIPisExprSum(scip, finalchildren->expr)); 1065 sum = finalchildren->next->expr; 1066 factors[0] = finalchildren->expr; 1067 } 1068 constant = simplifiedcoef * SCIPgetConstantExprSum(sum); 1069 nchildren = SCIPexprGetNChildren(sum); 1070 1071 SCIP_CALL( SCIPcreateExprSum(scip, &expanded, 1, &factors[0], &constant, 0.0, ownercreate, ownercreatedata) ); 1072 /* we are just re-using a child here, so do not release factor! */ 1073 1074 for( j = 0; j < nchildren; ++j ) 1075 { 1076 SCIP_Real coef; 1077 SCIP_Real termcoef; 1078 SCIP_Bool dummy; 1079 EXPRNODE* finalfactors; 1080 SCIP_EXPR* term = NULL; 1081 1082 coef = SCIPgetCoefsExprSum(sum)[j]; 1083 factors[1] = SCIPexprGetChildren(sum)[j]; 1084 1085 termcoef = coef; 1086 SCIP_CALL( simplifyMultiplyChildren(scip, factors, 2, &termcoef, &finalfactors, &dummy, ownercreate, ownercreatedata) ); 1087 assert(termcoef != 0.0); 1088 1089 SCIP_CALL( buildSimplifiedProduct(scip, 1.0, &finalfactors, expandalways, TRUE, &term, ownercreate, ownercreatedata) ); 1090 assert(finalfactors == NULL); 1091 assert(term != NULL); 1092 1093 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, termcoef * simplifiedcoef) ); 1094 SCIP_CALL( SCIPreleaseExpr(scip, &term) ); 1095 } 1096 1097 /* simplify the sum */ 1098 SCIP_CALL( SCIPcallExprSimplify(scip, expanded, simplifiedexpr, ownercreate, ownercreatedata) ); 1099 SCIP_CALL( SCIPreleaseExpr(scip, &expanded) ); 1100 } 1101 1102 return SCIP_OKAY; 1103 } 1104 1105 /* expands product of one sum and other expressions 1106 * -) (prod factor ... factor (sum c s1 ... sn) factor ... factor ) 1107 * - c != 0 --> c * factor is simplified (i.e. factor is not sum!) 1108 * - factor ... factor * si may be not be simplified, so put them in a product list and simplify them from there 1109 */ 1110 static 1111 SCIP_RETCODE enforceSP12b( 1112 SCIP* scip, /**< SCIP data structure */ 1113 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * prod simplifiedfactors */ 1114 EXPRNODE* finalchildren, /**< factors of simplified product */ 1115 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 1116 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 1117 void* ownercreatedata /**< data to pass to ownercreate */ 1118 ) 1119 { 1120 EXPRNODE* sum_node = NULL; 1121 EXPRNODE* n; 1122 int nfactors = 0; 1123 SCIP_EXPR** factors; 1124 int nchildren; 1125 SCIP_EXPR* expanded; 1126 int j; 1127 1128 /* check whether there is exactly one sum, calc number of factors */ 1129 for( n = finalchildren; n != NULL; n = n->next ) 1130 { 1131 if( SCIPisExprSum(scip, n->expr) ) 1132 { 1133 if( sum_node == NULL ) 1134 sum_node = n; 1135 else 1136 return SCIP_OKAY; /* more than one sum */ 1137 } 1138 else 1139 { 1140 ++nfactors; 1141 } 1142 } 1143 if( sum_node == NULL || nfactors == 0 ) /* no sum or no other factors */ 1144 return SCIP_OKAY; 1145 1146 /* collect exprs of all factors other than the sum */ 1147 SCIP_CALL( SCIPallocBufferArray(scip, &factors, nfactors + 1) ); 1148 for( n = finalchildren, j = 0; n != NULL; n = n->next ) 1149 if( n != sum_node ) 1150 factors[j++] = n->expr; 1151 1152 /* build new sum expression */ 1153 nchildren = SCIPexprGetNChildren(sum_node->expr); 1154 SCIP_CALL( SCIPcreateExprSum(scip, &expanded, 0, NULL, NULL, 0.0, ownercreate, ownercreatedata) ); 1155 1156 /* handle constant-from-sum * factors */ 1157 if( SCIPgetConstantExprSum(sum_node->expr) != 0.0 ) 1158 { 1159 if( nfactors == 1 ) 1160 { 1161 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, factors[0], simplifiedcoef * SCIPgetConstantExprSum(sum_node->expr)) ); 1162 } 1163 else 1164 { 1165 SCIP_Real termcoef = 1.0; 1166 SCIP_Bool dummy; 1167 EXPRNODE* finalfactors; 1168 SCIP_EXPR* term = NULL; 1169 1170 SCIP_CALL( simplifyMultiplyChildren(scip, factors, nfactors, &termcoef, &finalfactors, &dummy, ownercreate, ownercreatedata) ); 1171 assert(termcoef != 0.0); 1172 1173 SCIP_CALL( buildSimplifiedProduct(scip, 1.0, &finalfactors, TRUE, TRUE, &term, ownercreate, ownercreatedata) ); 1174 assert(finalfactors == NULL); 1175 assert(term != NULL); 1176 1177 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, termcoef * simplifiedcoef * SCIPgetConstantExprSum(sum_node->expr)) ); 1178 SCIP_CALL( SCIPreleaseExpr(scip, &term) ); 1179 } 1180 } 1181 1182 for( j = 0; j < nchildren; ++j ) 1183 { 1184 SCIP_Real coef; 1185 SCIP_Real termcoef; 1186 SCIP_Bool dummy; 1187 EXPRNODE* finalfactors; 1188 SCIP_EXPR* term = NULL; 1189 1190 coef = SCIPgetCoefsExprSum(sum_node->expr)[j]; 1191 factors[nfactors] = SCIPexprGetChildren(sum_node->expr)[j]; 1192 1193 termcoef = coef; 1194 SCIP_CALL( simplifyMultiplyChildren(scip, factors, nfactors + 1, &termcoef, &finalfactors, &dummy, ownercreate, ownercreatedata) ); 1195 assert(termcoef != 0.0); 1196 1197 SCIP_CALL( buildSimplifiedProduct(scip, 1.0, &finalfactors, TRUE, TRUE, &term, ownercreate, ownercreatedata) ); 1198 assert(finalfactors == NULL); 1199 assert(term != NULL); 1200 1201 SCIP_CALL( SCIPappendExprSumExpr(scip, expanded, term, termcoef * simplifiedcoef) ); 1202 SCIP_CALL( SCIPreleaseExpr(scip, &term) ); 1203 } 1204 1205 /* simplify the sum */ 1206 SCIP_CALL( SCIPcallExprSimplify(scip, expanded, simplifiedexpr, ownercreate, ownercreatedata) ); 1207 SCIP_CALL( SCIPreleaseExpr(scip, &expanded) ); 1208 1209 SCIPfreeBufferArray(scip, &factors); 1210 1211 return SCIP_OKAY; 1212 } 1213 1214 /** builds a simplified product from simplifiedfactors 1215 * 1216 * @note this function also releases simplifiedfactors 1217 */ 1218 static 1219 SCIP_RETCODE buildSimplifiedProduct( 1220 SCIP* scip, /**< SCIP data structure */ 1221 SCIP_Real simplifiedcoef, /**< simplified product should be simplifiedcoef * PI simplifiedfactors */ 1222 EXPRNODE** simplifiedfactors, /**< factors of simplified product */ 1223 SCIP_Bool expandalways, /**< whether to expand products of a sum and several factors in simplify (SP12b) */ 1224 SCIP_Bool changed, /**< indicates whether some of the simplified factors was changed */ 1225 SCIP_EXPR** simplifiedexpr, /**< buffer to store the simplified expression */ 1226 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 1227 void* ownercreatedata /**< data to pass to ownercreate */ 1228 ) 1229 { 1230 EXPRNODE* finalchildren = *simplifiedfactors; 1231 1232 /* build product expression from finalchildren and post-simplify */ 1233 debugSimplify("[simplifyProduct] finalchildren has length %d\n", listLength(finalchildren)); 1234 1235 *simplifiedexpr = NULL; 1236 1237 SCIP_CALL( enforceSP11(scip, simplifiedcoef, *simplifiedfactors, simplifiedexpr, ownercreate, ownercreatedata) ); 1238 if( *simplifiedexpr != NULL ) 1239 goto CLEANUP; 1240 1241 SCIP_CALL( enforceSP12(scip, simplifiedcoef, *simplifiedfactors, expandalways, simplifiedexpr, ownercreate, ownercreatedata) ); 1242 if( *simplifiedexpr != NULL ) 1243 goto CLEANUP; 1244 1245 if( expandalways ) 1246 { 1247 SCIP_CALL( enforceSP12b(scip, simplifiedcoef, *simplifiedfactors, simplifiedexpr, ownercreate, ownercreatedata) ); 1248 if( *simplifiedexpr != NULL ) 1249 goto CLEANUP; 1250 } 1251 1252 SCIP_CALL( enforceSP10(scip, simplifiedcoef, *simplifiedfactors, simplifiedexpr, ownercreate, ownercreatedata) ); 1253 if( *simplifiedexpr != NULL ) 1254 goto CLEANUP; 1255 1256 /* enforces SP8: if simplifiedcoef != 1.0, transform it into a sum with the (simplified) product as child */ 1257 if( simplifiedcoef != 1.0 ) 1258 { 1259 SCIP_EXPR* aux; 1260 SCIP_EXPR* sum; 1261 1262 /* create sum */ 1263 SCIP_CALL( createExprProductFromExprlist(scip, finalchildren, 1.0, &aux, ownercreate, ownercreatedata) ); 1264 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, &aux, &simplifiedcoef, 0.0, ownercreate, ownercreatedata) ); 1265 SCIP_CALL( SCIPreleaseExpr(scip, &aux) ); 1266 1267 /* simplify sum */ 1268 SCIP_CALL( SCIPcallExprSimplify(scip, sum, simplifiedexpr, ownercreate, ownercreatedata) ); 1269 SCIP_CALL( SCIPreleaseExpr(scip, &sum) ); 1270 1271 goto CLEANUP; 1272 } 1273 1274 /* build product expression from list */ 1275 if( changed ) 1276 { 1277 SCIP_CALL( createExprProductFromExprlist(scip, finalchildren, simplifiedcoef, simplifiedexpr, ownercreate, ownercreatedata) ); 1278 goto CLEANUP; 1279 } 1280 1281 CLEANUP: 1282 1283 SCIP_CALL( freeExprlist(scip, simplifiedfactors) ); 1284 return SCIP_OKAY; 1285 } 1286 1287 /** computes an estimator for a product as a vertex polyhedral function 1288 * 1289 * Since the product is multilinear, its convex and concave envelopes are piecewise linear. 1290 */ 1291 static 1292 SCIP_RETCODE estimateVertexPolyhedralProduct( 1293 SCIP* scip, /**< SCIP data structure */ 1294 SCIP_CONSHDLR* conshdlr, /**< nonlinear constraint handler */ 1295 int nfactors, /**< number of factors */ 1296 SCIP_INTERVAL* bounds, /**< bound for each factor */ 1297 SCIP_Real constantfactor, /**< another constant factor */ 1298 SCIP_Real* refpoint, /**< reference point where to estimate, or NULL if called from initestimates */ 1299 SCIP_Bool overestimate, /**< should estimator overestimate expr (TRUE) or underestimate (FALSE) */ 1300 SCIP_Real targetvalue, /**< no need to compute facet if value in xstar would be worse than target value */ 1301 SCIP_Real* coefs, /**< array to store cut coefficients */ 1302 SCIP_Real* constant, /**< pointer to store cut constant */ 1303 SCIP_Bool* success /**< pointer to store whether estimation was successful */ 1304 ) 1305 { 1306 SCIP_Real* box; 1307 SCIP_Real* xstar; 1308 int nfixed; 1309 int i; 1310 1311 assert(conshdlr != NULL); 1312 assert(nfactors > 0); 1313 assert(bounds != NULL); 1314 assert(constantfactor != 0.0); 1315 assert(coefs != NULL); 1316 assert(constant != NULL); 1317 assert(success != NULL); 1318 1319 *success = FALSE; 1320 1321 /* assemble box, check for unbounded variables, assemble xstar */ 1322 SCIP_CALL( SCIPallocBufferArray(scip, &box, 2*nfactors) ); 1323 SCIP_CALL( SCIPallocBufferArray(scip, &xstar, nfactors) ); 1324 for( i = 0, nfixed = 0; i < nfactors; ++i ) 1325 { 1326 assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bounds[i])); 1327 1328 if( SCIPisInfinity(scip, -bounds[i].inf) || SCIPisInfinity(scip, bounds[i].sup) ) 1329 { 1330 SCIPdebugMsg(scip, "a factor is unbounded, no cut is possible\n"); 1331 goto CLEANUP; 1332 } 1333 1334 box[2*i] = bounds[i].inf; 1335 box[2*i+1] = bounds[i].sup; 1336 1337 xstar[i] = refpoint != NULL ? refpoint[i] : 0.5 * (box[2*i] + box[2*i+1]); 1338 1339 if( SCIPisRelEQ(scip, box[2*i], box[2*i+1]) ) 1340 ++nfixed; 1341 } 1342 1343 if( nfixed < nfactors && nfactors - nfixed <= SCIP_MAXVERTEXPOLYDIM ) 1344 { 1345 SCIP_CALL( SCIPcomputeFacetVertexPolyhedralNonlinear(scip, conshdlr, 1346 overestimate, prodfunction, &constantfactor, xstar, box, nfactors, targetvalue, success, coefs, constant) ); 1347 } 1348 1349 CLEANUP: 1350 SCIPfreeBufferArray(scip, &xstar); 1351 SCIPfreeBufferArray(scip, &box); 1352 1353 return SCIP_OKAY; 1354 } 1355 1356 /* 1357 * Callback methods of expression handler 1358 */ 1359 1360 /** simplifies a product expression 1361 * 1362 * Summary: we first build a list of expressions (called finalchildren) which will be the children of the simplified product 1363 * and then we process this list in order to enforce SP8 and SP10. 1364 * 1365 * Description: In order to build finalchildren, we first build a list of unsimplified children (called unsimplifiedchildren) 1366 * with the children of the product. Each node of the list is manipulated (see simplifyFactor) in order to satisfy 1367 * SP2 and SP7 as follows: 1368 * - SP7: if the node's expression is a value, multiply the value to the products's coef 1369 * - SP2: if the node's expression is a product, then build a list with the child's children 1370 * 1371 * Then, we merge the built list (or the simplified node) into finalchildren. While merging, nodes from finalchildren 1372 * can go back to unsimplifiedchildren for further processing (see mergeProductExprlist() for more details). 1373 * After building finalchildren, we create the simplified product out of it, taking care that SP8 and SP10 are satisfied 1374 */ 1375 static 1376 SCIP_DECL_EXPRSIMPLIFY(simplifyProduct) 1377 { /*lint --e{715}*/ 1378 EXPRNODE* finalchildren; 1379 SCIP_Real simplifiedcoef; 1380 SCIP_Bool changed; 1381 1382 assert(expr != NULL); 1383 assert(simplifiedexpr != NULL); 1384 1385 simplifiedcoef = SCIPgetCoefExprProduct(expr); 1386 1387 #ifdef SIMPLIFY_DEBUG 1388 debugSimplify("Simplifying expr:\n"); 1389 SCIP_CALL( SCIPprintExpr(scip, expr, NULL) ); 1390 SCIPinfoMessage(scip, NULL, "\n"); 1391 debugSimplify("First multiplying children\n"); 1392 #endif 1393 1394 /* simplify and multiply factors */ 1395 SCIP_CALL( simplifyMultiplyChildren(scip, SCIPexprGetChildren(expr), SCIPexprGetNChildren(expr), &simplifiedcoef, 1396 &finalchildren, &changed, ownercreate, ownercreatedata) ); 1397 1398 #ifdef SIMPLIFY_DEBUG 1399 { 1400 EXPRNODE* node; 1401 int i; 1402 1403 debugSimplify("Building product from simplified factors\n"); 1404 node = finalchildren; 1405 i = 0; 1406 while( node != NULL ) 1407 { 1408 debugSimplify("factor %d:\n", i); 1409 SCIP_CALL( SCIPprintExpr(scip, node->expr, NULL) ); 1410 SCIPinfoMessage(scip, NULL, "\n"); 1411 node = node->next; 1412 i++; 1413 } 1414 } 1415 #endif 1416 1417 /* get simplified product from simplified factors in finalchildren */ 1418 SCIP_CALL( buildSimplifiedProduct(scip, simplifiedcoef, &finalchildren, SCIPexprhdlrGetData(SCIPexprGetHdlr(expr))->expandalways, changed, simplifiedexpr, ownercreate, 1419 ownercreatedata) ); 1420 assert(finalchildren == NULL); 1421 1422 if( *simplifiedexpr == NULL ) 1423 { 1424 *simplifiedexpr = expr; 1425 1426 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */ 1427 SCIPcaptureExpr(*simplifiedexpr); 1428 } 1429 assert(*simplifiedexpr != NULL); 1430 1431 return SCIP_OKAY; 1432 } 1433 1434 /** compare two product expressions 1435 * 1436 * The order of two product expressions, u and v, is a lexicographical order on the factors. 1437 * 1438 * Starting from the *last*, we find the first child where they differ, say, the i-th. 1439 * Then u < v <=> u_i < v_i. 1440 * If there is no such children and they have different number of children, then u < v <=> nchildren(u) < nchildren(v). 1441 * If all children are the same and they have the same number of children, then u < v <=> coeff(u) < coeff(v). 1442 * Otherwise, they are the same. 1443 * 1444 * Note: we are assuming expression are simplified, so within u, we have u_1 < u_2, etc. 1445 * 1446 * Example: y * z < x * y * z 1447 */ 1448 static 1449 SCIP_DECL_EXPRCOMPARE(compareProduct) 1450 { /*lint --e{715}*/ 1451 int compareresult; 1452 int i; 1453 int j; 1454 int nchildren1; 1455 int nchildren2; 1456 SCIP_EXPR** children1; 1457 SCIP_EXPR** children2; 1458 1459 nchildren1 = SCIPexprGetNChildren(expr1); 1460 nchildren2 = SCIPexprGetNChildren(expr2); 1461 children1 = SCIPexprGetChildren(expr1); 1462 children2 = SCIPexprGetChildren(expr2); 1463 1464 for( i = nchildren1 - 1, j = nchildren2 - 1; i >= 0 && j >= 0; --i, --j ) 1465 { 1466 compareresult = SCIPcompareExpr(scip, children1[i], children2[j]); 1467 if( compareresult != 0 ) 1468 return compareresult; 1469 /* expressions are equal, continue */ 1470 } 1471 1472 /* all children of one expression are children of the other expression, use number of children as a tie-breaker */ 1473 if( i < j ) 1474 { 1475 assert(i == -1); 1476 /* expr1 has less elements, hence expr1 < expr2 */ 1477 return -1; 1478 } 1479 if( i > j ) 1480 { 1481 assert(j == -1); 1482 /* expr1 has more elements, hence expr1 > expr2 */ 1483 return 1; 1484 } 1485 1486 /* everything is equal, use coefficient as tie-breaker */ 1487 assert(i == -1 && j == -1); 1488 if( SCIPgetCoefExprProduct(expr1) < SCIPgetCoefExprProduct(expr2) ) 1489 return -1; 1490 if( SCIPgetCoefExprProduct(expr1) > SCIPgetCoefExprProduct(expr2) ) 1491 return 1; 1492 1493 /* they are equal */ 1494 return 0; 1495 } 1496 1497 /** expression handler copy callback */ 1498 static 1499 SCIP_DECL_EXPRCOPYHDLR(copyhdlrProduct) 1500 { /*lint --e{715}*/ 1501 SCIP_CALL( SCIPincludeExprhdlrProduct(scip) ); 1502 1503 return SCIP_OKAY; 1504 } 1505 1506 /** expression handler free callback */ 1507 static 1508 SCIP_DECL_EXPRFREEHDLR(freehdlrProduct) 1509 { /*lint --e{715}*/ 1510 assert(scip != NULL); 1511 assert(exprhdlr != NULL); 1512 assert(exprhdlrdata != NULL); 1513 assert(*exprhdlrdata != NULL); 1514 1515 SCIPfreeBlockMemory(scip, exprhdlrdata); 1516 assert(*exprhdlrdata == NULL); 1517 1518 return SCIP_OKAY; 1519 } 1520 1521 /** expression data copy callback */ 1522 static 1523 SCIP_DECL_EXPRCOPYDATA(copydataProduct) 1524 { /*lint --e{715}*/ 1525 SCIP_EXPRDATA* sourceexprdata; 1526 1527 assert(targetexprdata != NULL); 1528 assert(sourceexpr != NULL); 1529 1530 sourceexprdata = SCIPexprGetData(sourceexpr); 1531 assert(sourceexprdata != NULL); 1532 1533 SCIP_CALL( SCIPduplicateBlockMemory(targetscip, targetexprdata, sourceexprdata) ); 1534 1535 return SCIP_OKAY; 1536 } 1537 1538 /** expression data free callback */ 1539 static 1540 SCIP_DECL_EXPRFREEDATA(freedataProduct) 1541 { /*lint --e{715}*/ 1542 SCIP_EXPRDATA* exprdata; 1543 1544 assert(expr != NULL); 1545 1546 exprdata = SCIPexprGetData(expr); 1547 assert(exprdata != NULL); 1548 1549 SCIPfreeBlockMemory(scip, &exprdata); 1550 1551 SCIPexprSetData(expr, NULL); 1552 1553 return SCIP_OKAY; 1554 } 1555 1556 /** expression print callback */ 1557 static 1558 SCIP_DECL_EXPRPRINT(printProduct) 1559 { /*lint --e{715}*/ 1560 SCIP_EXPRDATA* exprdata; 1561 1562 assert(expr != NULL); 1563 1564 exprdata = SCIPexprGetData(expr); 1565 assert(exprdata != NULL); 1566 1567 switch( stage ) 1568 { 1569 case SCIP_EXPRITER_ENTEREXPR : 1570 { 1571 /* print opening parenthesis, if necessary */ 1572 if( EXPRHDLR_PRECEDENCE <= parentprecedence ) 1573 { 1574 SCIPinfoMessage(scip, file, "("); 1575 } 1576 1577 /* print coefficient, if not one */ 1578 if( exprdata->coefficient != 1.0 ) 1579 { 1580 if( exprdata->coefficient < 0.0 && EXPRHDLR_PRECEDENCE > parentprecedence ) 1581 { 1582 SCIPinfoMessage(scip, file, "(%g)", exprdata->coefficient); 1583 } 1584 else 1585 { 1586 SCIPinfoMessage(scip, file, "%g", exprdata->coefficient); 1587 } 1588 } 1589 break; 1590 } 1591 1592 case SCIP_EXPRITER_VISITINGCHILD : 1593 { 1594 /* print multiplication sign, if not first factor */ 1595 if( exprdata->coefficient != 1.0 || currentchild > 0 ) 1596 { 1597 SCIPinfoMessage(scip, file, "*"); 1598 } 1599 break; 1600 } 1601 1602 case SCIP_EXPRITER_VISITEDCHILD : 1603 { 1604 break; 1605 } 1606 1607 case SCIP_EXPRITER_LEAVEEXPR : 1608 { 1609 /* print closing parenthesis, if necessary */ 1610 if( EXPRHDLR_PRECEDENCE <= parentprecedence ) 1611 { 1612 SCIPinfoMessage(scip, file, ")"); 1613 } 1614 break; 1615 } 1616 1617 default: 1618 /* all stages should have been covered above */ 1619 SCIPABORT(); 1620 } 1621 1622 return SCIP_OKAY; 1623 } 1624 1625 /** product hash callback */ 1626 static 1627 SCIP_DECL_EXPRHASH(hashProduct) 1628 { /*lint --e{715}*/ 1629 SCIP_EXPRDATA* exprdata; 1630 int c; 1631 1632 assert(scip != NULL); 1633 assert(expr != NULL); 1634 assert(hashkey != NULL); 1635 assert(childrenhashes != NULL); 1636 1637 exprdata = SCIPexprGetData(expr); 1638 assert(exprdata != NULL); 1639 1640 *hashkey = EXPRHDLR_HASHKEY; 1641 *hashkey ^= SCIPcalcFibHash(exprdata->coefficient); 1642 1643 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1644 *hashkey ^= childrenhashes[c]; 1645 1646 return SCIP_OKAY; 1647 } 1648 1649 /** expression point evaluation callback */ 1650 static 1651 SCIP_DECL_EXPREVAL(evalProduct) 1652 { /*lint --e{715}*/ 1653 SCIP_EXPRDATA* exprdata; 1654 SCIP_Real childval; 1655 int c; 1656 1657 assert(expr != NULL); 1658 1659 exprdata = SCIPexprGetData(expr); 1660 assert(exprdata != NULL); 1661 1662 *val = exprdata->coefficient; 1663 for( c = 0; c < SCIPexprGetNChildren(expr) && (*val != 0.0); ++c ) 1664 { 1665 childval = SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]); 1666 assert(childval != SCIP_INVALID); 1667 1668 *val *= childval; 1669 } 1670 1671 return SCIP_OKAY; 1672 } 1673 1674 /** derivative evaluation callback computing <gradient, children.dot> 1675 * 1676 * If expr is \f$\prod_i x_i\f$, then computes \f$\sum_j \prod_{i\neq j} x_i x^{\text{dot}}_j\f$. 1677 */ 1678 static 1679 SCIP_DECL_EXPRFWDIFF(fwdiffProduct) 1680 { /*lint --e{715}*/ 1681 int c; 1682 1683 assert(expr != NULL); 1684 assert(dot != NULL); 1685 1686 assert(SCIPexprGetData(expr) != NULL); 1687 1688 /* TODO add special handling for nchildren == 2 */ 1689 1690 /**! [SnippetExprFwdiffProduct] */ 1691 *dot = 0.0; 1692 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1693 { 1694 SCIP_EXPR* child; 1695 1696 child = SCIPexprGetChildren(expr)[c]; 1697 1698 assert(SCIPexprGetEvalValue(child) != SCIP_INVALID); 1699 assert(SCIPexprGetDot(child) != SCIP_INVALID); 1700 1701 if( SCIPexprGetDot(child) == 0.0 ) 1702 continue; 1703 1704 if( SCIPexprGetEvalValue(child) != 0.0 ) 1705 *dot += SCIPexprGetEvalValue(expr) / SCIPexprGetEvalValue(child) * SCIPexprGetDot(child); 1706 else 1707 { 1708 SCIP_Real partial; 1709 int i; 1710 1711 partial = SCIPexprGetData(expr)->coefficient; 1712 for( i = 0; i < SCIPexprGetNChildren(expr) && (partial != 0.0); ++i ) 1713 { 1714 if( i == c ) 1715 continue; 1716 1717 partial *= SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[i]); 1718 } 1719 *dot += partial * SCIPexprGetDot(child); 1720 } 1721 } 1722 /**! [SnippetExprFwdiffProduct] */ 1723 1724 return SCIP_OKAY; 1725 } 1726 1727 /** expression backward forward derivative evaluation callback 1728 * 1729 * Computes \f$\frac{\partial}{\partial \text{childidx}} ( \langle \text{gradient}, \text{children.dot}\rangle )\f$. 1730 * 1731 * If expr is \f$\prod_i x_i\f$, and childidx is \f$k\f$ then computes 1732 * \f$\partial_k \sum_j \prod_{i \neq j} x_i x^{\text{dot}}_j 1733 * = \sum_{j \neq k} \prod_{i \neq j, k} x_i x^{\text{dot}}_j\f$ 1734 */ 1735 static 1736 SCIP_DECL_EXPRBWFWDIFF(bwfwdiffProduct) 1737 { /*lint --e{715}*/ 1738 SCIP_EXPR* partialchild; 1739 int c; 1740 1741 assert(expr != NULL); 1742 assert(bardot != NULL); 1743 assert(SCIPexprGetData(expr) != NULL); 1744 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr)); 1745 1746 partialchild = SCIPexprGetChildren(expr)[childidx]; 1747 assert(partialchild != NULL); 1748 assert(!SCIPisExprValue(scip, partialchild)); 1749 assert(SCIPexprGetEvalValue(partialchild) != SCIP_INVALID); 1750 1751 /* TODO add special handling for nchildren == 2 */ 1752 1753 /**! [SnippetExprBwfwdiffProduct] */ 1754 *bardot = 0.0; 1755 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1756 { 1757 SCIP_EXPR* child; 1758 1759 if( c == childidx ) 1760 continue; 1761 1762 child = SCIPexprGetChildren(expr)[c]; 1763 1764 assert(SCIPexprGetEvalValue(child) != SCIP_INVALID); 1765 assert(SCIPexprGetDot(child) != SCIP_INVALID); 1766 1767 if( SCIPexprGetDot(child) == 0.0 ) 1768 continue; 1769 1770 if( SCIPexprGetEvalValue(child) != 0.0 && SCIPexprGetEvalValue(partialchild) != 0.0 ) 1771 *bardot += SCIPexprGetEvalValue(expr) / (SCIPexprGetEvalValue(child) * SCIPexprGetEvalValue(partialchild)) * SCIPexprGetDot(child); 1772 else 1773 { 1774 SCIP_Real partial; 1775 int i; 1776 1777 partial = SCIPexprGetData(expr)->coefficient; 1778 for( i = 0; i < SCIPexprGetNChildren(expr) && (partial != 0.0); ++i ) 1779 { 1780 if( i == c || i == childidx ) 1781 continue; 1782 1783 partial *= SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[i]); 1784 } 1785 *bardot += partial * SCIPexprGetDot(child); 1786 } 1787 } 1788 /**! [SnippetExprBwfwdiffProduct] */ 1789 1790 return SCIP_OKAY; 1791 } 1792 1793 /** expression derivative evaluation callback */ 1794 static 1795 SCIP_DECL_EXPRBWDIFF(bwdiffProduct) 1796 { /*lint --e{715}*/ 1797 SCIP_EXPR* child; 1798 1799 assert(expr != NULL); 1800 assert(SCIPexprGetData(expr) != NULL); 1801 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr)); 1802 1803 child = SCIPexprGetChildren(expr)[childidx]; 1804 assert(child != NULL); 1805 assert(!SCIPisExprValue(scip, child)); 1806 assert(SCIPexprGetEvalValue(child) != SCIP_INVALID); 1807 1808 /* TODO add special handling for nchildren == 2 */ 1809 1810 /**! [SnippetExprBwdiffProduct] */ 1811 if( !SCIPisZero(scip, SCIPexprGetEvalValue(child)) ) 1812 { 1813 *val = SCIPexprGetEvalValue(expr) / SCIPexprGetEvalValue(child); 1814 } 1815 else 1816 { 1817 int i; 1818 1819 *val = SCIPexprGetData(expr)->coefficient; 1820 for( i = 0; i < SCIPexprGetNChildren(expr) && (*val != 0.0); ++i ) 1821 { 1822 if( i == childidx ) 1823 continue; 1824 1825 *val *= SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[i]); 1826 } 1827 } 1828 /**! [SnippetExprBwdiffProduct] */ 1829 1830 return SCIP_OKAY; 1831 } 1832 1833 /** expression interval evaluation callback */ 1834 static 1835 SCIP_DECL_EXPRINTEVAL(intevalProduct) 1836 { /*lint --e{715}*/ 1837 SCIP_EXPRDATA* exprdata; 1838 int c; 1839 1840 assert(expr != NULL); 1841 1842 exprdata = SCIPexprGetData(expr); 1843 assert(exprdata != NULL); 1844 1845 /**! [SnippetExprIntevalProduct] */ 1846 SCIPintervalSet(interval, exprdata->coefficient); 1847 1848 SCIPdebugMsg(scip, "inteval %p with %d children: %.20g", (void*)expr, SCIPexprGetNChildren(expr), exprdata->coefficient); 1849 1850 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1851 { 1852 SCIP_INTERVAL childinterval; 1853 1854 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[c]); 1855 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) ) 1856 { 1857 SCIPintervalSetEmpty(interval); 1858 break; 1859 } 1860 1861 /* multiply childinterval with the so far computed interval */ 1862 SCIPintervalMul(SCIP_INTERVAL_INFINITY, interval, *interval, childinterval); 1863 1864 SCIPdebugMsgPrint(scip, " *[%.20g,%.20g]", childinterval.inf, childinterval.sup); 1865 } 1866 SCIPdebugMsgPrint(scip, " = [%.20g,%.20g]\n", interval->inf, interval->sup); 1867 /**! [SnippetExprIntevalProduct] */ 1868 1869 return SCIP_OKAY; 1870 } 1871 1872 /** estimates a multilinear function of the form \f$ f(x) := a \prod_{i = 1}^n x_i \f$ 1873 * 1874 * \f$ x_i \f$ are the auxiliary variables of the children. 1875 * If !overestimate, then we look for an affine underestimator of \f$ f(x) \f$ which has a value above targetvalue at \f$ x^* \f$, 1876 * i.e., \f$ g(x) := \alpha^T x + \beta \le f(x)\f$ for all \f$ x \f$ in the domain, such that \f$ \alpha x^* + \beta > \text{targetvalue}\f$. 1877 * 1878 * Since \f$ f(x) \f$ is componentwise linear, its convex envelope is piecewise linear and its value can be computed by 1879 * finding the largest affine underestimator. 1880 * This is done either explicitly (if n=2) or by solving an LP, see SCIPcomputeFacetVertexPolyhedralNonlinear(). 1881 */ 1882 static 1883 SCIP_DECL_EXPRESTIMATE(estimateProduct) 1884 { /*lint --e{715}*/ 1885 SCIP_EXPRDATA* exprdata; 1886 int nchildren; 1887 1888 assert(scip != NULL); 1889 assert(expr != NULL); 1890 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 1891 assert(refpoint != NULL); 1892 assert(coefs != NULL); 1893 assert(constant != NULL); 1894 assert(islocal != NULL); 1895 assert(branchcand != NULL); 1896 assert(*branchcand == TRUE); 1897 assert(success != NULL); 1898 1899 exprdata = SCIPexprGetData(expr); 1900 assert(exprdata != NULL); 1901 1902 *success = FALSE; 1903 *islocal = TRUE; 1904 1905 nchildren = SCIPexprGetNChildren(expr); 1906 1907 /* debug output: prints expression we are trying to estimate, bounds of variables and point */ 1908 #ifdef SCIP_DEBUG 1909 { 1910 int c; 1911 1912 SCIPdebugMsg(scip, "%sestimating product with %d variables\n", overestimate ? "over": "under", SCIPexprGetNChildren(expr)); 1913 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1914 { 1915 SCIPdebugMsg(scip, "child %d = %g in [%g, %g]\n", c, refpoint[c], localbounds[c].inf, localbounds[c].sup); 1916 1917 if( SCIPisInfinity(scip, localbounds[c].sup) || SCIPisInfinity(scip, -localbounds[c].inf) ) 1918 { 1919 SCIPdebugMsg(scip, "unbounded factor related to\n"); 1920 SCIP_CALL( SCIPdismantleExpr(scip, NULL, SCIPexprGetChildren(expr)[0]) ); 1921 } 1922 } 1923 } 1924 #endif 1925 1926 /* bilinear term */ 1927 if( nchildren == 2 ) 1928 { 1929 SCIP_Real refpointx; 1930 SCIP_Real refpointy; 1931 SCIP_INTERVAL bndx; 1932 SCIP_INTERVAL bndy; 1933 1934 /* collect first variable */ 1935 refpointx = refpoint[0]; 1936 bndx = localbounds[0]; 1937 assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bndx)); 1938 1939 /* collect second variable */ 1940 refpointy = refpoint[1]; 1941 bndy = localbounds[1]; 1942 assert(!SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, bndy)); 1943 1944 /* adjust the reference points */ 1945 refpointx = MIN(MAX(refpointx, bndx.inf), bndx.sup); 1946 refpointy = MIN(MAX(refpointy, bndy.inf), bndy.sup); 1947 1948 coefs[0] = 0.0; 1949 coefs[1] = 0.0; 1950 *constant = 0.0; 1951 *success = TRUE; 1952 1953 SCIPaddBilinMcCormick(scip, exprdata->coefficient, bndx.inf, bndx.sup, refpointx, 1954 bndy.inf, bndy.sup, refpointy, overestimate, &coefs[0], &coefs[1], constant, 1955 success); 1956 } 1957 else 1958 { 1959 SCIP_EXPRHDLRDATA* exprhdlrdata; 1960 1961 exprhdlrdata = SCIPexprhdlrGetData(SCIPexprGetHdlr(expr)); 1962 assert(exprhdlrdata != NULL); 1963 1964 if( exprhdlrdata->conshdlr != NULL ) 1965 { 1966 SCIP_CALL( estimateVertexPolyhedralProduct(scip, exprhdlrdata->conshdlr, nchildren, localbounds, exprdata->coefficient, refpoint, overestimate, 1967 targetvalue, coefs, constant, success) ); 1968 } 1969 else 1970 { 1971 SCIPdebugMsg(scip, "no cons_nonlinear included in SCIP, cannot estimate vertex-polyhedral product function\n"); 1972 } 1973 } 1974 1975 return SCIP_OKAY; 1976 } 1977 1978 /** initial estimators callback */ 1979 static 1980 SCIP_DECL_EXPRINITESTIMATES(initestimatesProduct) 1981 { 1982 SCIP_EXPRDATA* exprdata; 1983 SCIP_Bool success = TRUE; 1984 int nchildren; 1985 1986 assert(scip != NULL); 1987 assert(expr != NULL); 1988 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 1989 assert(nreturned != NULL); 1990 1991 nchildren = SCIPexprGetNChildren(expr); 1992 1993 exprdata = SCIPexprGetData(expr); 1994 assert(exprdata != NULL); 1995 1996 if( nchildren == 2 ) 1997 { 1998 SCIP_INTERVAL bndx = bounds[0]; 1999 SCIP_INTERVAL bndy = bounds[1]; 2000 2001 constant[0] = 0.0; 2002 coefs[0][0] = 0.0; 2003 coefs[0][1] = 0.0; 2004 2005 /* build estimator */ 2006 SCIPaddBilinMcCormick(scip, exprdata->coefficient, bndx.inf, bndx.sup, (bndx.inf + bndx.sup) / 2.0, 2007 bndy.inf, bndy.sup, (bndy.inf + bndy.sup ) / 2.0, overestimate, &coefs[0][0], &coefs[0][1], 2008 constant, &success); 2009 } 2010 else 2011 { 2012 SCIP_EXPRHDLRDATA* exprhdlrdata; 2013 2014 exprhdlrdata = SCIPexprhdlrGetData(SCIPexprGetHdlr(expr)); 2015 assert(exprhdlrdata != NULL); 2016 2017 if( exprhdlrdata->conshdlr != NULL ) 2018 { 2019 SCIP_CALL( estimateVertexPolyhedralProduct(scip, exprhdlrdata->conshdlr, nchildren, bounds, exprdata->coefficient, NULL, overestimate, 2020 overestimate ? SCIPinfinity(scip) : -SCIPinfinity(scip), coefs[0], constant, &success) ); 2021 } 2022 else 2023 { 2024 SCIPdebugMsg(scip, "no cons_nonlinear included in SCIP, cannot estimate vertex-polyhedral product function\n"); 2025 } 2026 } 2027 2028 if( success ) 2029 *nreturned = 1; 2030 2031 return SCIP_OKAY; 2032 } 2033 2034 /** expression reverse propagation callback */ 2035 static 2036 SCIP_DECL_EXPRREVERSEPROP(reversepropProduct) 2037 { /*lint --e{715}*/ 2038 SCIP_EXPRDATA* exprdata; 2039 SCIP_INTERVAL childbounds; 2040 SCIP_INTERVAL otherfactor; 2041 SCIP_INTERVAL zero; 2042 int i; 2043 int j; 2044 2045 assert(scip != NULL); 2046 assert(expr != NULL); 2047 assert(SCIPexprGetNChildren(expr) > 0); 2048 assert(infeasible != NULL); 2049 assert(childrenbounds != NULL); 2050 2051 *infeasible = FALSE; 2052 2053 /* too expensive (runtime here is quadratic in number of children) 2054 * TODO implement something faster for larger numbers of factors, e.g., split product into smaller products 2055 */ 2056 if( SCIPexprGetNChildren(expr) > 10 ) 2057 return SCIP_OKAY; 2058 2059 /* not possible to learn bounds on children if expression bounds are unbounded in both directions */ 2060 if( SCIPintervalIsEntire(SCIP_INTERVAL_INFINITY, bounds) ) 2061 return SCIP_OKAY; 2062 2063 exprdata = SCIPexprGetData(expr); 2064 assert(exprdata != NULL); 2065 2066 /**! [SnippetExprReversepropProduct] */ 2067 SCIPintervalSet(&zero, 0.0); 2068 2069 /* f = const * prod_k c_k => c_i solves c_i * (const * prod_{j:j!=i} c_j) = f */ 2070 for( i = 0; i < SCIPexprGetNChildren(expr) && !(*infeasible); ++i ) 2071 { 2072 SCIPintervalSet(&otherfactor, exprdata->coefficient); 2073 2074 /* compute prod_{j:j!=i} c_j */ 2075 for( j = 0; j < SCIPexprGetNChildren(expr); ++j ) 2076 { 2077 if( i == j ) 2078 continue; 2079 2080 /* TODO we should compute these only one time instead of repeating this for almost every i */ 2081 childbounds = childrenbounds[j]; 2082 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childbounds) ) 2083 { 2084 *infeasible = TRUE; 2085 return SCIP_OKAY; 2086 } 2087 2088 SCIPintervalMul(SCIP_INTERVAL_INFINITY, &otherfactor, otherfactor, childbounds); 2089 } 2090 2091 childbounds = childrenbounds[i]; 2092 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childbounds) ) 2093 { 2094 *infeasible = TRUE; 2095 return SCIP_OKAY; 2096 } 2097 2098 /* solve x*otherfactor = f for x in c_i */ 2099 SCIPintervalSolveUnivariateQuadExpression(SCIP_INTERVAL_INFINITY, &childbounds, zero, otherfactor, bounds, childbounds); 2100 2101 SCIPdebugMsg(scip, "child %d: solved [%g,%g]*x = [%g,%g] with x in [%g,%g] -> x = [%g,%g]\n", i, otherfactor.inf, otherfactor.sup, 2102 bounds.inf, bounds.sup, 2103 childrenbounds[i].inf, childrenbounds[i].sup, 2104 childbounds.inf, childbounds.sup); 2105 2106 /* store computed bounds of the expression */ 2107 SCIPintervalIntersect(&childrenbounds[i], childrenbounds[i], childbounds); 2108 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childrenbounds[i]) ) 2109 { 2110 *infeasible = TRUE; 2111 return SCIP_OKAY; 2112 } 2113 } 2114 /**! [SnippetExprReversepropProduct] */ 2115 2116 return SCIP_OKAY; 2117 } 2118 2119 /** expression curvature detection callback */ 2120 static 2121 SCIP_DECL_EXPRCURVATURE(curvatureProduct) 2122 { /*lint --e{715}*/ 2123 assert(scip != NULL); 2124 assert(expr != NULL); 2125 assert(success != NULL); 2126 assert(SCIPexprGetNChildren(expr) > 0); 2127 2128 if( SCIPexprGetNChildren(expr) == 1 ) 2129 { 2130 *childcurv = SCIPexprcurvMultiply(SCIPgetCoefExprProduct(expr), exprcurvature); 2131 *success = TRUE; 2132 } 2133 else 2134 { 2135 *success = FALSE; 2136 } 2137 2138 return SCIP_OKAY; 2139 } 2140 2141 /** expression monotonicity detection callback */ 2142 static 2143 SCIP_DECL_EXPRMONOTONICITY(monotonicityProduct) 2144 { /*lint --e{715}*/ 2145 SCIP_Real coef; 2146 int i; 2147 int nneg; 2148 2149 assert(scip != NULL); 2150 assert(expr != NULL); 2151 assert(result != NULL); 2152 assert(SCIPexprGetNChildren(expr) >= 1); 2153 assert(childidx >= 0); 2154 assert(childidx < SCIPexprGetNChildren(expr)); 2155 2156 coef = SCIPgetCoefExprProduct(expr); 2157 2158 /* count the number of negative children (except for childidx); if some children changes sign 2159 * -> monotonicity unknown 2160 */ 2161 nneg = 0; 2162 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 2163 { 2164 SCIP_INTERVAL interval; 2165 2166 if( i == childidx ) 2167 continue; 2168 2169 assert(SCIPexprGetChildren(expr)[i] != NULL); 2170 SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(expr)[i]) ); 2171 interval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[i]); 2172 2173 if( SCIPintervalGetSup(interval) <= 0.0 ) 2174 nneg++; 2175 else if( SCIPintervalGetInf(interval) < 0.0 ) 2176 { 2177 *result = SCIP_MONOTONE_UNKNOWN; 2178 return SCIP_OKAY; 2179 } 2180 } 2181 2182 /* note that the monotonicity depends on the sign of the coefficient */ 2183 if( nneg % 2 == 0 ) 2184 *result = (coef >= 0.0) ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC; 2185 else 2186 *result = (coef >= 0.0) ? SCIP_MONOTONE_DEC : SCIP_MONOTONE_INC; 2187 2188 return SCIP_OKAY; 2189 } 2190 2191 /** expression integrality detection callback */ 2192 static 2193 SCIP_DECL_EXPRINTEGRALITY(integralityProduct) 2194 { /*lint --e{715}*/ 2195 SCIP_EXPRDATA* exprdata; 2196 int i; 2197 2198 assert(scip != NULL); 2199 assert(expr != NULL); 2200 assert(isintegral != NULL); 2201 2202 exprdata = SCIPexprGetData(expr); 2203 assert(exprdata != NULL); 2204 2205 *isintegral = EPSISINT(exprdata->coefficient, 0.0); /*lint !e835*/ 2206 2207 for( i = 0; i < SCIPexprGetNChildren(expr) && *isintegral; ++i ) 2208 { 2209 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i]; 2210 assert(child != NULL); 2211 2212 *isintegral = SCIPexprIsIntegral(child); 2213 } 2214 2215 return SCIP_OKAY; 2216 } 2217 2218 /** expression callback to get information for symmetry detection */ 2219 static 2220 SCIP_DECL_EXPRGETSYMDATA(getSymDataProduct) 2221 { /*lint --e{715}*/ 2222 assert(scip != NULL); 2223 assert(expr != NULL); 2224 assert(symdata != NULL); 2225 2226 SCIP_CALL( SCIPallocBlockMemory(scip, symdata) ); 2227 2228 (*symdata)->nconstants = 1; 2229 (*symdata)->ncoefficients = 0; 2230 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->constants, 1) ); 2231 (*symdata)->constants[0] = SCIPgetCoefExprProduct(expr); 2232 2233 return SCIP_OKAY; 2234 } 2235 2236 /** creates the handler for product expressions and includes it into SCIP */ 2237 SCIP_RETCODE SCIPincludeExprhdlrProduct( 2238 SCIP* scip /**< SCIP data structure */ 2239 ) 2240 { 2241 SCIP_EXPRHDLRDATA* exprhdlrdata; 2242 SCIP_EXPRHDLR* exprhdlr; 2243 2244 /* allocate expression handler data */ 2245 SCIP_CALL( SCIPallocClearBlockMemory(scip, &exprhdlrdata) ); 2246 exprhdlrdata->conshdlr = SCIPfindConshdlr(scip, "nonlinear"); 2247 2248 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalProduct, 2249 exprhdlrdata) ); 2250 assert(exprhdlr != NULL); 2251 2252 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrProduct, freehdlrProduct); 2253 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataProduct, freedataProduct); 2254 SCIPexprhdlrSetSimplify(exprhdlr, simplifyProduct); 2255 SCIPexprhdlrSetCompare(exprhdlr, compareProduct); 2256 SCIPexprhdlrSetPrint(exprhdlr, printProduct); 2257 SCIPexprhdlrSetIntEval(exprhdlr, intevalProduct); 2258 SCIPexprhdlrSetEstimate(exprhdlr, initestimatesProduct, estimateProduct); 2259 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropProduct); 2260 SCIPexprhdlrSetHash(exprhdlr, hashProduct); 2261 SCIPexprhdlrSetDiff(exprhdlr, bwdiffProduct, fwdiffProduct, bwfwdiffProduct); 2262 SCIPexprhdlrSetCurvature(exprhdlr, curvatureProduct); 2263 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityProduct); 2264 SCIPexprhdlrSetIntegrality(exprhdlr, integralityProduct); 2265 SCIPexprhdlrSetGetSymdata(exprhdlr, getSymDataProduct); 2266 2267 SCIP_CALL( SCIPaddBoolParam(scip, "expr/" EXPRHDLR_NAME "/expandalways", 2268 "whether to expand products of a sum and several factors in simplify", 2269 &exprhdlrdata->expandalways, FALSE, FALSE, NULL, NULL) ); 2270 2271 return SCIP_OKAY; 2272 } 2273 2274 /** creates a product expression */ 2275 SCIP_RETCODE SCIPcreateExprProduct( 2276 SCIP* scip, /**< SCIP data structure */ 2277 SCIP_EXPR** expr, /**< pointer where to store expression */ 2278 int nchildren, /**< number of children */ 2279 SCIP_EXPR** children, /**< children */ 2280 SCIP_Real coefficient, /**< constant coefficient of product */ 2281 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 2282 void* ownercreatedata /**< data to pass to ownercreate */ 2283 ) 2284 { 2285 SCIP_EXPRDATA* exprdata; 2286 2287 /**! [SnippetCreateExprProduct] */ 2288 SCIP_CALL( SCIPallocBlockMemory(scip, &exprdata) ); 2289 exprdata->coefficient = coefficient; 2290 2291 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrProduct(scip), exprdata, nchildren, children, ownercreate, ownercreatedata) ); 2292 /**! [SnippetCreateExprProduct] */ 2293 2294 return SCIP_OKAY; 2295 } 2296 2297 /* from pub_expr.h */ 2298 2299 /** gets the constant coefficient of a product expression */ 2300 SCIP_Real SCIPgetCoefExprProduct( 2301 SCIP_EXPR* expr /**< product expression */ 2302 ) 2303 { 2304 SCIP_EXPRDATA* exprdata; 2305 2306 assert(expr != NULL); 2307 2308 exprdata = SCIPexprGetData(expr); 2309 assert(exprdata != NULL); 2310 2311 return exprdata->coefficient; 2312 } 2313