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_sum.c 26 * @ingroup DEFPLUGINS_EXPR 27 * @brief sum expression handler 28 * @author Stefan Vigerske 29 * @author Benjamin Mueller 30 * @author Felipe Serrano 31 */ 32 33 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 34 35 #include <string.h> 36 #include <stddef.h> 37 38 #include "scip/expr_sum.h" 39 #include "scip/expr_value.h" 40 #include "scip/expr_product.h" 41 #include "scip/expr_exp.h" 42 #include "scip/expr_pow.h" 43 #include "symmetry/struct_symmetry.h" 44 45 #define EXPRHDLR_NAME "sum" 46 #define EXPRHDLR_DESC "summation with coefficients and a constant" 47 #define EXPRHDLR_PRECEDENCE 40000 48 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(47161.0) 49 50 /** macro to activate/deactivate debugging information of simplify method */ 51 /*lint -emacro(774,debugSimplify) */ 52 #ifdef SIMPLIFY_DEBUG 53 #define debugSimplify printf 54 #else 55 #define debugSimplify while( FALSE ) printf 56 #endif 57 58 /* 59 * Data structures 60 */ 61 62 /** expression data */ 63 struct SCIP_ExprData 64 { 65 SCIP_Real constant; /**< constant coefficient */ 66 SCIP_Real* coefficients; /**< coefficients of children */ 67 int coefssize; /**< size of the coefficients array */ 68 }; 69 70 /* 71 * Local methods 72 */ 73 74 /** creates expression data */ 75 static 76 SCIP_RETCODE createData( 77 SCIP* scip, /**< SCIP data structure */ 78 SCIP_EXPRDATA** exprdata, /**< pointer where to store expression data */ 79 int ncoefficients, /**< number of coefficients (i.e., number of children) */ 80 SCIP_Real* coefficients, /**< array with coefficients for all children (or NULL if all 1.0) */ 81 SCIP_Real constant /**< constant term of sum */ 82 ) 83 { 84 assert(exprdata != NULL); 85 assert(ncoefficients >= 0); 86 87 SCIP_CALL( SCIPallocBlockMemory(scip, exprdata) ); 88 89 if( coefficients != NULL ) 90 { 91 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*exprdata)->coefficients, coefficients, ncoefficients) ); 92 } 93 else 94 { 95 int i; 96 97 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*exprdata)->coefficients, ncoefficients) ); 98 for( i = 0; i < ncoefficients; ++i ) 99 (*exprdata)->coefficients[i] = 1.0; 100 } 101 102 (*exprdata)->coefssize = ncoefficients; 103 (*exprdata)->constant = constant; 104 105 return SCIP_OKAY; 106 } 107 108 /** simplifies the `idx`-th child of the sum expression `duplicate` in order for it to be able to be a child of a simplified sum 109 * 110 * for example, this means that the `idx`-th child cannot be itself a sum 111 * if it is, we have to flatten it, i.e., take all its children and make them children of `duplicate` 112 */ 113 static 114 SCIP_RETCODE simplifyTerm( 115 SCIP* scip, /**< SCIP data structure */ 116 SCIP_EXPR* duplicate, /**< expression to be simplified */ 117 int idx, /**< idx of children to be simplified */ 118 SCIP_Bool* changed, /**< pointer to store if some term actually got simplified */ 119 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 120 void* ownercreatedata /**< data to pass to ownercreate */ 121 ) 122 { 123 SCIP_EXPR** children; 124 SCIP_EXPR* expr; 125 SCIP_Real* coefs; 126 SCIP_Real constant ; 127 SCIP_Real coef; 128 129 assert(duplicate != NULL); 130 assert(idx >= 0); 131 assert(idx < SCIPexprGetNChildren(duplicate)); 132 assert(changed != NULL); 133 134 children = SCIPexprGetChildren(duplicate); 135 coefs = SCIPgetCoefsExprSum(duplicate); 136 constant = SCIPgetConstantExprSum(duplicate); 137 138 coef = coefs[idx]; 139 expr = children[idx]; 140 assert(expr != NULL); 141 142 /* enforces SS3 */ 143 if( SCIPisExprValue(scip, expr) ) 144 { 145 *changed = TRUE; 146 constant += coef * SCIPgetValueExprValue(expr); 147 SCIPsetConstantExprSum(duplicate, constant); 148 149 /* TODO: remove child? */ 150 coefs[idx] = 0.0; 151 152 return SCIP_OKAY; 153 } 154 155 /* enforces SS2 */ 156 if( SCIPisExprSum(scip, expr) ) 157 { 158 *changed = TRUE; 159 160 /* pass constant to parent */ 161 constant += coef * SCIPgetConstantExprSum(expr); 162 SCIPsetConstantExprSum(duplicate, constant); 163 164 /* append all children of expr on parent except the first one */ 165 if( SCIPexprGetNChildren(expr) > 1 ) 166 { 167 int i; 168 169 for( i = 1; i < SCIPexprGetNChildren(expr); ++i ) 170 { 171 assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[i])); 172 SCIP_CALL( SCIPappendExprSumExpr(scip, duplicate, SCIPexprGetChildren(expr)[i], 173 coef * SCIPgetCoefsExprSum(expr)[i]) ); 174 } 175 } 176 177 /* replace expr with first child; need to get data again since it might be re-allocated */ 178 assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[0])); 179 180 coefs = SCIPgetCoefsExprSum(duplicate); 181 182 coefs[idx] = coef * SCIPgetCoefsExprSum(expr)[0]; 183 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, SCIPexprGetChildren(expr)[0]) ); 184 185 return SCIP_OKAY; 186 } 187 188 /* enforce SS9 */ 189 if( REALABS(coef) != 1.0 && SCIPisExprProduct(scip, expr) ) 190 { 191 SCIP_EXPR* expchild = NULL; 192 int i; 193 194 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 195 { 196 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i]; 197 assert(child != NULL); 198 199 if( SCIPisExprExp(scip, child) ) 200 { 201 expchild = child; 202 break; 203 } 204 } 205 206 /* coef != +- 1, term is product and one factor is an exponential -> enforce SS9 */ 207 if( expchild != NULL ) 208 { 209 SCIP_EXPR* sum; 210 SCIP_EXPR* prod; 211 SCIP_EXPR* simplifiedprod; 212 SCIP_EXPR* simplifiedsum; 213 SCIP_EXPR* exponential; 214 SCIP_EXPR* simplifiedexp; 215 SCIP_Real expconstant; 216 217 /* inform that expression will change */ 218 *changed = TRUE; 219 220 /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */ 221 if( coef > 0.0 ) 222 { 223 expconstant = log(coef); 224 coefs[idx] = 1.0; 225 } 226 else 227 { 228 expconstant = log(-coef); 229 coefs[idx] = -1.0; 230 } 231 232 /* add constant to exponential's child */ 233 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expchild), NULL, expconstant, ownercreate, 234 ownercreatedata) ); 235 236 /* simplify sum */ 237 SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) ); 238 SCIP_CALL( SCIPreleaseExpr(scip, &sum) ); 239 240 /* create exponential with new child */ 241 SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) ); 242 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) ); 243 244 /* simplify exponential */ 245 SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) ); 246 SCIP_CALL( SCIPreleaseExpr(scip, &exponential) ); 247 248 /* create product with new child */ 249 SCIP_CALL( SCIPcreateExprProduct(scip, &prod, 0, NULL, 1.0, ownercreate, ownercreatedata) ); 250 251 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 252 { 253 if( SCIPexprGetChildren(expr)[i] == expchild ) 254 { 255 SCIP_CALL( SCIPappendExprChild(scip, prod, simplifiedexp) ); 256 } 257 else 258 { 259 SCIP_CALL( SCIPappendExprChild(scip, prod, SCIPexprGetChildren(expr)[i]) ); 260 } 261 } 262 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) ); 263 264 /* simplify product */ 265 SCIP_CALL( SCIPcallExprSimplify(scip, prod, &simplifiedprod, ownercreate, ownercreatedata) ); 266 SCIP_CALL( SCIPreleaseExpr(scip, &prod) ); 267 268 /* replace current child with simplified product */ 269 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedprod) ); 270 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedprod) ); 271 272 /* since the simplified product can be a sum ( exp(-1)*exp(log(x+y)+1) -> x+y ), 273 * we call the function we are in again 274 * this is no endless recursion, since the coef is now +- 1 275 */ 276 SCIP_CALL( simplifyTerm(scip, duplicate, idx, changed, ownercreate, ownercreatedata) ); 277 278 return SCIP_OKAY; 279 } 280 } 281 282 /* enforce SS10 */ 283 if( REALABS(coef) != 1.0 && SCIPisExprExp(scip, expr) ) 284 { 285 /* coef != +- 1, term is exponential -> enforce SS10 by moving |coef| into argument of exponential */ 286 287 SCIP_EXPR* sum; 288 SCIP_EXPR* simplifiedsum; 289 SCIP_EXPR* exponential; 290 SCIP_EXPR* simplifiedexp; 291 SCIP_Real expconstant; 292 293 /* inform that expression will change */ 294 *changed = TRUE; 295 296 /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */ 297 if( coef > 0.0 ) 298 { 299 expconstant = log(coef); 300 coefs[idx] = 1.0; 301 } 302 else 303 { 304 expconstant = log(-coef); 305 coefs[idx] = -1.0; 306 } 307 308 /* add constant to exponential's child */ 309 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expr), NULL, expconstant, ownercreate, 310 ownercreatedata) ); /* expconstant+expchild */ 311 312 /* simplify sum */ 313 SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) ); 314 SCIP_CALL( SCIPreleaseExpr(scip, &sum) ); 315 316 /* create exponential with new child */ 317 SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) ); 318 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) ); 319 320 /* simplify exponential */ 321 SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) ); 322 SCIP_CALL( SCIPreleaseExpr(scip, &exponential) ); 323 324 /* replace current child with simplified exponential */ 325 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedexp) ); 326 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) ); 327 328 return SCIP_OKAY; 329 } 330 331 /* other types of (simplified) expressions can be a child of a simplified sum */ 332 assert(!SCIPisExprSum(scip, expr)); 333 assert(!SCIPisExprValue(scip, expr)); 334 335 return SCIP_OKAY; 336 } 337 338 /** helper struct for expressions sort */ 339 typedef struct 340 { 341 SCIP* scip; /**< SCIP data structure */ 342 SCIP_EXPR** exprs; /**< expressions */ 343 } SORTEXPRDATA; 344 345 static 346 SCIP_DECL_SORTINDCOMP(sortExprComp) 347 { 348 SORTEXPRDATA* data = (SORTEXPRDATA*)dataptr; 349 350 return SCIPcompareExpr(data->scip, data->exprs[ind1], data->exprs[ind2]); 351 } 352 353 /* 354 * Callback methods of expression handler 355 */ 356 357 /** simplifies a sum expression 358 * 359 * goes through each child and simplifies it; then sorts the simplified children; then sum the children that are equal; 360 * finally creates a sum expression with all the children that do not have a 0 coefficient and post-process so that SS6 361 * and SS7 are satisfied 362 */ 363 static 364 SCIP_DECL_EXPRSIMPLIFY(simplifySum) 365 { /*lint --e{715}*/ 366 SCIP_EXPR** children; 367 SCIP_EXPR* duplicate = NULL; 368 SCIP_EXPR** newchildren = NULL; 369 SCIP_Real* newcoefs = NULL; 370 int nnewchildren; 371 SCIP_Real newconstant; 372 SCIP_Real* coefs; 373 int i; 374 int nchildren; 375 SCIP_Bool changed; 376 SORTEXPRDATA sortdata; 377 int* order = NULL; 378 379 assert(expr != NULL); 380 assert(simplifiedexpr != NULL); 381 assert(SCIPexprGetHdlr(expr) == SCIPgetExprhdlrSum(scip)); 382 383 changed = FALSE; 384 385 /* TODO: maybe have a flag to know if it is simplified ? */ 386 /* TODO: can we do this with a shallow duplicate + copy of children pointer? currently simplifyTerm may modify children, 387 * so one would need to be careful 388 */ 389 SCIP_CALL( SCIPduplicateExpr(scip, expr, &duplicate, NULL, NULL, ownercreate, ownercreatedata) ); 390 assert(duplicate != NULL); 391 392 nchildren = SCIPexprGetNChildren(duplicate); 393 for( i = 0; i < nchildren; i++ ) 394 { 395 /* enforces SS8 TODO: remove child? */ 396 /* we have to ask for the coefs everytime, since it might get realloced in simpifyTerm */ 397 if( SCIPgetCoefsExprSum(duplicate)[i] == 0.0 ) 398 { 399 changed = TRUE; 400 continue; 401 } 402 403 /* enforces SS2, SS3, SS9, and SS10 */ 404 SCIP_CALL( simplifyTerm(scip, duplicate, i, &changed, ownercreate, ownercreatedata) ); 405 } 406 407 /* simplifyTerm can add new children to duplicate and realloc them; so get them again */ 408 nchildren = SCIPexprGetNChildren(duplicate); 409 children = SCIPexprGetChildren(duplicate); 410 coefs = SCIPgetCoefsExprSum(duplicate); 411 412 /* treat zero term case */ 413 if( nchildren == 0 ) 414 { 415 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) ); 416 goto CLEANUP; 417 } 418 419 /* treat one term case */ 420 if( nchildren == 1 ) 421 { 422 if( coefs[0] == 0.0 ) 423 { 424 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) ); 425 goto CLEANUP; 426 } 427 428 if( coefs[0] == 1.0 && SCIPgetConstantExprSum(duplicate) == 0.0 ) 429 *simplifiedexpr = children[0]; /* SS7 */ 430 else 431 *simplifiedexpr = changed ? duplicate : expr; 432 433 SCIPcaptureExpr(*simplifiedexpr); 434 435 goto CLEANUP; 436 } 437 438 /* enforces SS5: sort children */ 439 SCIP_CALL( SCIPallocBufferArray(scip, &order, nchildren) ); 440 for( i = 0; i < nchildren; i++ ) 441 order[i] = i; 442 sortdata.scip = scip; 443 sortdata.exprs = children; 444 SCIPsortInd(order, sortExprComp, (void*)&sortdata, nchildren); 445 446 /* create sorted variant of children and coefs */ 447 SCIP_CALL( SCIPallocBufferArray(scip, &newchildren, nchildren) ); 448 SCIP_CALL( SCIPallocBufferArray(scip, &newcoefs, nchildren) ); 449 for( i = 0; i < nchildren; ++i ) 450 { 451 newchildren[i] = children[order[i]]; 452 newcoefs[i] = coefs[order[i]]; 453 if( order[i] != i ) 454 changed = TRUE; 455 } 456 457 /* post-process */ 458 459 /* enforces SS4 */ 460 nnewchildren = 0; 461 for( i = 0; i < nchildren; i++ ) 462 { 463 /* eliminate zero-coefficients */ 464 if( newcoefs[i] == 0.0 ) 465 { 466 changed = TRUE; 467 continue; 468 } 469 470 /* sum equal expressions */ 471 if( i < nchildren-1 && SCIPcompareExpr(scip, newchildren[i], newchildren[i+1]) == 0 ) 472 { 473 changed = TRUE; 474 /* if we substract two almost equal not-so-small numbers, then set new coefficient to 0.0 475 * instead of some tiny value that is likely the result of some random round-off error 476 * E.g., on instance ex1221, we have x1^2 + b3 = 1.25. 477 * Probing finds an aggregation x1 = 1.11803 - 0.618034 b3. 478 * Simplify would then produce 1.25 + 1e-16 x1 = 1.25. 479 */ 480 if( SCIPisEQ(scip, newcoefs[i], -newcoefs[i+1]) && REALABS(newcoefs[i]) >= 1.0 ) 481 newcoefs[i+1] = 0.0; 482 else 483 newcoefs[i+1] += newcoefs[i]; 484 continue; 485 } 486 487 /* move i-th child to new position */ 488 newchildren[nnewchildren] = newchildren[i]; 489 newcoefs[nnewchildren] = newcoefs[i]; 490 nnewchildren++; 491 } 492 493 /* build sum expression from finalchildren and post-simplify */ 494 newconstant = SCIPgetConstantExprSum(duplicate); 495 496 debugSimplify("what to do? finalchildren has length %d\n", nnewchildren); /*lint !e506 !e681*/ 497 498 /* enforces SS6: if they are no children, return value */ 499 if( nnewchildren == 0 ) 500 { 501 debugSimplify("[sum] got empty list, return value %g\n", newconstant); /*lint !e506 !e681*/ 502 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, newconstant, ownercreate, ownercreatedata) ); 503 504 goto CLEANUP; 505 } 506 507 /* enforces SS7: if list consists of one expr with coef 1.0 and constant is 0, return that expr */ 508 if( nnewchildren == 1 && newcoefs[0] == 1.0 && newconstant == 0.0 ) 509 { 510 *simplifiedexpr = newchildren[0]; 511 SCIPcaptureExpr(*simplifiedexpr); 512 513 goto CLEANUP; 514 } 515 516 /* build sum expression from children */ 517 if( changed ) 518 { 519 SCIP_CALL( SCIPcreateExprSum(scip, simplifiedexpr, nnewchildren, newchildren, newcoefs, newconstant, 520 ownercreate, ownercreatedata) ); 521 522 goto CLEANUP; 523 } 524 525 *simplifiedexpr = expr; 526 527 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */ 528 SCIPcaptureExpr(*simplifiedexpr); 529 530 /* free memory */ 531 CLEANUP: 532 SCIPfreeBufferArrayNull(scip, &newcoefs); 533 SCIPfreeBufferArrayNull(scip, &newchildren); 534 SCIPfreeBufferArrayNull(scip, &order); 535 SCIP_CALL( SCIPreleaseExpr(scip, &duplicate) ); 536 537 return SCIP_OKAY; 538 } 539 540 /** expression callback to get information for symmetry detection */ 541 static 542 SCIP_DECL_EXPRGETSYMDATA(getSymDataSum) 543 { /*lint --e{715}*/ 544 SCIP_EXPRDATA* exprdata; 545 int i; 546 547 assert(scip != NULL); 548 assert(expr != NULL); 549 550 exprdata = SCIPexprGetData(expr); 551 assert(exprdata != NULL); 552 553 SCIP_CALL( SCIPallocBlockMemory(scip, symdata) ); 554 555 (*symdata)->nconstants = 1; 556 (*symdata)->ncoefficients = exprdata->coefssize; 557 558 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->constants, 1) ); 559 (*symdata)->constants[0] = exprdata->constant; 560 561 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->coefficients, exprdata->coefssize) ); 562 for( i = 0; i < exprdata->coefssize; ++i ) 563 (*symdata)->coefficients[i] = exprdata->coefficients[i]; 564 565 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*symdata)->children, exprdata->coefssize) ); 566 for( i = 0; i < exprdata->coefssize; ++i ) 567 (*symdata)->children[i] = SCIPexprGetChildren(expr)[i]; 568 569 return SCIP_OKAY; 570 } 571 572 /** compares two sum expressions 573 * 574 * The order of two sum expressions is a lexicographical order on the terms. 575 * 576 * Starting from the *last*, we find the first child where they differ, say, the i-th. 577 * Then u < v <=> u_i < v_i. 578 * If there are no such children and they have different number of children, then u < v <=> nchildren(u) < nchildren(v). 579 * If there are no such children and they have the same number of children, then u < v <=> const(u) < const(v). 580 * Otherwise, they are the same. 581 * 582 * Note: we are assuming expression are simplified, so within u, we have u_1 < u_2, etc 583 * 584 * Example: y + z < x + y + z, 2*x + 3*y < 3*x + 3*y 585 */ 586 static 587 SCIP_DECL_EXPRCOMPARE(compareSum) 588 { /*lint --e{715}*/ 589 SCIP_Real const1; 590 SCIP_Real* coefs1; 591 SCIP_EXPR** children1; 592 int nchildren1; 593 SCIP_Real const2; 594 SCIP_Real* coefs2; 595 SCIP_EXPR** children2; 596 int nchildren2; 597 int compareresult; 598 int i; 599 int j; 600 601 nchildren1 = SCIPexprGetNChildren(expr1); 602 nchildren2 = SCIPexprGetNChildren(expr2); 603 children1 = SCIPexprGetChildren(expr1); 604 children2 = SCIPexprGetChildren(expr2); 605 coefs1 = SCIPgetCoefsExprSum(expr1); 606 coefs2 = SCIPgetCoefsExprSum(expr2); 607 const1 = SCIPgetConstantExprSum(expr1); 608 const2 = SCIPgetConstantExprSum(expr2); 609 610 for( i = nchildren1 - 1, j = nchildren2 - 1; i >= 0 && j >= 0; --i, --j ) 611 { 612 compareresult = SCIPcompareExpr(scip, children1[i], children2[j]); 613 if( compareresult != 0 ) 614 return compareresult; 615 else 616 { 617 /* expressions are equal, compare coefficient */ 618 if( (coefs1 ? coefs1[i] : 1.0) < (coefs2 ? coefs2[j] : 1.0) ) 619 return -1; 620 if( (coefs1 ? coefs1[i] : 1.0) > (coefs2 ? coefs2[j] : 1.0) ) 621 return 1; 622 623 /* coefficients are equal, continue */ 624 } 625 } 626 627 /* all children of one expression are children of the other expression, use number of children as a tie-breaker */ 628 if( i < j ) 629 { 630 assert(i == -1); 631 /* expr1 has less elements, hence expr1 < expr2 */ 632 return -1; 633 } 634 if( i > j ) 635 { 636 assert(j == -1); 637 /* expr1 has more elements, hence expr1 > expr2 */ 638 return 1; 639 } 640 641 /* everything is equal, use constant/coefficient as tie-breaker */ 642 assert(i == -1 && j == -1); 643 if( const1 < const2 ) 644 return -1; 645 if( const1 > const2 ) 646 return 1; 647 648 /* they are equal */ 649 return 0; 650 } 651 652 /** expression handler copy callback */ 653 static 654 SCIP_DECL_EXPRCOPYHDLR(copyhdlrSum) 655 { /*lint --e{715}*/ 656 SCIP_CALL( SCIPincludeExprhdlrSum(scip) ); 657 658 return SCIP_OKAY; 659 } 660 661 /** expression data copy callback */ 662 static 663 SCIP_DECL_EXPRCOPYDATA(copydataSum) 664 { /*lint --e{715}*/ 665 SCIP_EXPRDATA* sourceexprdata; 666 667 assert(targetexprdata != NULL); 668 assert(sourceexpr != NULL); 669 670 sourceexprdata = SCIPexprGetData(sourceexpr); 671 assert(sourceexprdata != NULL); 672 673 SCIP_CALL( createData(targetscip, targetexprdata, SCIPexprGetNChildren(sourceexpr), 674 sourceexprdata->coefficients, sourceexprdata->constant) ); 675 676 return SCIP_OKAY; 677 } 678 679 /** expression data free callback */ 680 static 681 SCIP_DECL_EXPRFREEDATA(freedataSum) 682 { /*lint --e{715}*/ 683 SCIP_EXPRDATA* exprdata; 684 685 assert(expr != NULL); 686 687 exprdata = SCIPexprGetData(expr); 688 assert(exprdata != NULL); 689 690 SCIPfreeBlockMemoryArray(scip, &(exprdata->coefficients), exprdata->coefssize); 691 SCIPfreeBlockMemory(scip, &exprdata); 692 693 SCIPexprSetData(expr, NULL); 694 695 return SCIP_OKAY; 696 } 697 698 /** expression print callback */ 699 static 700 SCIP_DECL_EXPRPRINT(printSum) 701 { /*lint --e{715}*/ 702 SCIP_EXPRDATA* exprdata; 703 704 assert(expr != NULL); 705 706 exprdata = SCIPexprGetData(expr); 707 assert(exprdata != NULL); 708 709 /**! [SnippetExprPrintSum] */ 710 switch( stage ) 711 { 712 case SCIP_EXPRITER_ENTEREXPR : 713 { 714 /* print opening parenthesis, if necessary */ 715 if( EXPRHDLR_PRECEDENCE <= parentprecedence ) 716 { 717 SCIPinfoMessage(scip, file, "("); 718 } 719 720 /* print constant, if nonzero */ 721 if( exprdata->constant != 0.0 ) 722 { 723 SCIPinfoMessage(scip, file, "%g", exprdata->constant); 724 } 725 break; 726 } 727 728 case SCIP_EXPRITER_VISITINGCHILD : 729 { 730 SCIP_Real coef; 731 732 coef = exprdata->coefficients[currentchild]; 733 734 /* print coefficient, if necessary */ 735 if( coef == 1.0 ) 736 { 737 /* if coefficient is 1.0, then print only "+" if not the first term */ 738 if( exprdata->constant != 0.0 || currentchild > 0 ) 739 { 740 SCIPinfoMessage(scip, file, "+"); 741 } 742 } 743 else if( coef == -1.0 ) 744 { 745 /* if coefficient is -1.0, then print only "-" */ 746 SCIPinfoMessage(scip, file, "-"); 747 } 748 else 749 { 750 /* force "+" sign on positive coefficient if not the first term */ 751 SCIPinfoMessage(scip, file, (exprdata->constant != 0.0 || currentchild > 0) ? "%+g*" : "%g*", coef); 752 } 753 754 break; 755 } 756 757 case SCIP_EXPRITER_LEAVEEXPR : 758 { 759 /* print closing parenthesis, if necessary */ 760 if( EXPRHDLR_PRECEDENCE <= parentprecedence ) 761 { 762 SCIPinfoMessage(scip, file, ")"); 763 } 764 break; 765 } 766 767 case SCIP_EXPRITER_VISITEDCHILD : 768 default: ; 769 } 770 /**! [SnippetExprPrintSum] */ 771 772 return SCIP_OKAY; 773 } 774 775 /** expression point evaluation callback */ 776 static 777 SCIP_DECL_EXPREVAL(evalSum) 778 { /*lint --e{715}*/ 779 SCIP_EXPRDATA* exprdata; 780 int c; 781 782 assert(expr != NULL); 783 784 exprdata = SCIPexprGetData(expr); 785 assert(exprdata != NULL); 786 787 /**! [SnippetExprEvalSum] */ 788 *val = exprdata->constant; 789 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 790 { 791 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/ 792 793 *val += exprdata->coefficients[c] * SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]); 794 } 795 /**! [SnippetExprEvalSum] */ 796 797 return SCIP_OKAY; 798 } 799 800 /** expression forward derivative evaluation callback */ 801 static 802 SCIP_DECL_EXPRFWDIFF(fwdiffSum) 803 { /*lint --e{715}*/ 804 SCIP_EXPRDATA* exprdata; 805 int c; 806 807 assert(expr != NULL); 808 assert(dot != NULL); 809 810 exprdata = SCIPexprGetData(expr); 811 assert(exprdata != NULL); 812 813 *dot = 0.0; 814 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 815 { 816 assert(SCIPexprGetDot(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/ 817 818 *dot += exprdata->coefficients[c] * SCIPexprGetDot(SCIPexprGetChildren(expr)[c]); 819 } 820 821 return SCIP_OKAY; 822 } 823 824 /** expression derivative evaluation callback */ 825 static 826 SCIP_DECL_EXPRBWDIFF(bwdiffSum) 827 { /*lint --e{715}*/ 828 assert(expr != NULL); 829 assert(SCIPexprGetData(expr) != NULL); 830 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr)); 831 assert(SCIPexprGetChildren(expr)[childidx] != NULL); 832 assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[childidx])); 833 834 *val = SCIPgetCoefsExprSum(expr)[childidx]; 835 836 return SCIP_OKAY; 837 } 838 839 /** expression backward forward derivative evaluation callback */ 840 static 841 SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSum) 842 { /*lint --e{715}*/ 843 assert(bardot != NULL); 844 845 *bardot = 0.0; 846 847 return SCIP_OKAY; 848 } 849 850 /** expression interval evaluation callback */ 851 static 852 SCIP_DECL_EXPRINTEVAL(intevalSum) 853 { /*lint --e{715}*/ 854 SCIP_EXPRDATA* exprdata; 855 SCIP_INTERVAL suminterval; 856 int c; 857 858 assert(expr != NULL); 859 860 exprdata = SCIPexprGetData(expr); 861 assert(exprdata != NULL); 862 863 SCIPintervalSet(interval, exprdata->constant); 864 865 SCIPdebugMsg(scip, "inteval %p with %d children: %.20g", (void*)expr, SCIPexprGetNChildren(expr), exprdata->constant); 866 867 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 868 { 869 SCIP_INTERVAL childinterval; 870 871 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[c]); 872 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) ) 873 { 874 SCIPintervalSetEmpty(interval); 875 break; 876 } 877 878 /* compute coefficients[c] * childinterval and add the result to the so far computed interval */ 879 if( exprdata->coefficients[c] == 1.0 ) 880 { 881 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, childinterval); 882 } 883 else 884 { 885 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &suminterval, childinterval, exprdata->coefficients[c]); 886 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, suminterval); 887 } 888 889 SCIPdebugMsgPrint(scip, " %+.20g*[%.20g,%.20g]", exprdata->coefficients[c], childinterval.inf, childinterval.sup); 890 } 891 SCIPdebugMsgPrint(scip, " = [%.20g,%.20g]\n", interval->inf, interval->sup); 892 893 return SCIP_OKAY; 894 } 895 896 /** initial estimators callback */ 897 static 898 SCIP_DECL_EXPRINITESTIMATES(initEstimatesSum) 899 { /*lint --e{715}*/ 900 SCIP_EXPRDATA* exprdata; 901 902 #ifdef SCIP_DEBUG 903 SCIPinfoMessage(scip, NULL, "initEstimatesSum %d children: ", SCIPexprGetNChildren(expr)); 904 SCIPprintExpr(scip, expr, NULL); 905 SCIPinfoMessage(scip, NULL, "\n"); 906 #endif 907 assert(scip != NULL); 908 assert(expr != NULL); 909 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 910 911 assert(coefs[0] != NULL); 912 assert(constant != NULL); 913 assert(nreturned != NULL); 914 915 exprdata = SCIPexprGetData(expr); 916 assert(exprdata != NULL); 917 918 BMScopyMemoryArray(coefs[0], exprdata->coefficients, SCIPexprGetNChildren(expr)); 919 *constant = exprdata->constant; 920 *nreturned = 1; 921 922 return SCIP_OKAY; 923 } 924 925 /** expression estimate callback */ 926 static 927 SCIP_DECL_EXPRESTIMATE(estimateSum) 928 { /*lint --e{715}*/ 929 SCIP_EXPRDATA* exprdata; 930 931 assert(scip != NULL); 932 assert(expr != NULL); 933 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0); 934 assert(islocal != NULL); 935 assert(success != NULL); 936 assert(branchcand != NULL); 937 938 exprdata = SCIPexprGetData(expr); 939 assert(exprdata != NULL); 940 941 /* NOTE: nlhdlr_default assumes in nlhdlrInitSepaDefault that this estimator can be used for both under- and overestimation */ 942 943 BMScopyMemoryArray(coefs, exprdata->coefficients, SCIPexprGetNChildren(expr)); 944 *constant = exprdata->constant; 945 *islocal = FALSE; 946 *success = TRUE; 947 948 /* for none of our children, branching would improve the underestimator, so set branchcand[i]=FALSE everywhere 949 * if we branch for numerical reasons, then cons-expr-core should figure out what the candidates are 950 */ 951 BMSclearMemoryArray(branchcand, SCIPexprGetNChildren(expr)); 952 953 return SCIP_OKAY; 954 } 955 956 /** expression reverse propagation callback */ 957 static 958 SCIP_DECL_EXPRREVERSEPROP(reversepropSum) 959 { /*lint --e{715}*/ 960 SCIP_EXPRDATA* exprdata; 961 SCIP_INTERVAL* newbounds; 962 int nchildren; 963 int nreductions; 964 965 assert(scip != NULL); 966 assert(expr != NULL); 967 assert(infeasible != NULL); 968 969 nchildren = SCIPexprGetNChildren(expr); 970 assert(nchildren > 0); 971 972 exprdata = SCIPexprGetData(expr); 973 assert(exprdata != NULL); 974 975 SCIP_CALL( SCIPallocBufferArray(scip, &newbounds, nchildren) ); 976 977 nreductions = SCIPintervalPropagateWeightedSum(SCIP_INTERVAL_INFINITY, nchildren, childrenbounds, 978 exprdata->coefficients, exprdata->constant, bounds, newbounds, infeasible); 979 980 if( !*infeasible && nreductions > 0 ) 981 BMScopyMemoryArray(childrenbounds, newbounds, nchildren); 982 983 SCIPfreeBufferArray(scip, &newbounds); 984 985 return SCIP_OKAY; 986 } 987 988 /** sum hash callback */ 989 static 990 SCIP_DECL_EXPRHASH(hashSum) 991 { /*lint --e{715}*/ 992 SCIP_EXPRDATA* exprdata; 993 int c; 994 995 assert(scip != NULL); 996 assert(expr != NULL); 997 assert(hashkey != NULL); 998 assert(childrenhashes != NULL); 999 1000 exprdata = SCIPexprGetData(expr); 1001 assert(exprdata != NULL); 1002 1003 /**! [SnippetExprHashSum] */ 1004 *hashkey = EXPRHDLR_HASHKEY; 1005 *hashkey ^= SCIPcalcFibHash(exprdata->constant); 1006 1007 for( c = 0; c < SCIPexprGetNChildren(expr); ++c ) 1008 *hashkey ^= SCIPcalcFibHash(exprdata->coefficients[c]) ^ childrenhashes[c]; 1009 /**! [SnippetExprHashSum] */ 1010 1011 return SCIP_OKAY; 1012 } 1013 1014 /** expression curvature detection callback */ 1015 static 1016 SCIP_DECL_EXPRCURVATURE(curvatureSum) 1017 { /*lint --e{715}*/ 1018 SCIP_EXPRDATA* exprdata; 1019 int i; 1020 1021 assert(scip != NULL); 1022 assert(expr != NULL); 1023 assert(childcurv != NULL); 1024 assert(success != NULL); 1025 1026 exprdata = SCIPexprGetData(expr); 1027 assert(exprdata != NULL); 1028 1029 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 1030 childcurv[i] = SCIPexprcurvMultiply(exprdata->coefficients[i], exprcurvature); 1031 1032 *success = TRUE; 1033 1034 return SCIP_OKAY; 1035 } 1036 1037 /** expression monotonicity detection callback */ 1038 static 1039 SCIP_DECL_EXPRMONOTONICITY(monotonicitySum) 1040 { /*lint --e{715}*/ 1041 SCIP_EXPRDATA* exprdata; 1042 1043 assert(scip != NULL); 1044 assert(expr != NULL); 1045 assert(result != NULL); 1046 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr)); 1047 1048 exprdata = SCIPexprGetData(expr); 1049 assert(exprdata != NULL); 1050 1051 *result = exprdata->coefficients[childidx] >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC; 1052 1053 return SCIP_OKAY; 1054 } 1055 1056 /** expression integrality detection callback */ 1057 static 1058 SCIP_DECL_EXPRINTEGRALITY(integralitySum) 1059 { /*lint --e{715}*/ 1060 SCIP_EXPRDATA* exprdata; 1061 int i; 1062 1063 assert(scip != NULL); 1064 assert(expr != NULL); 1065 assert(isintegral != NULL); 1066 1067 exprdata = SCIPexprGetData(expr); 1068 assert(exprdata != NULL); 1069 1070 /**! [SnippetExprIntegralitySum] */ 1071 *isintegral = EPSISINT(exprdata->constant, 0.0); /*lint !e835*/ 1072 1073 for( i = 0; i < SCIPexprGetNChildren(expr) && *isintegral; ++i ) 1074 { 1075 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i]; 1076 assert(child != NULL); 1077 1078 *isintegral = EPSISINT(exprdata->coefficients[i], 0.0) && SCIPexprIsIntegral(child); /*lint !e835*/ 1079 } 1080 /**! [SnippetExprIntegralitySum] */ 1081 1082 return SCIP_OKAY; 1083 } 1084 1085 /** creates the handler for sum expressions and includes it into SCIP */ 1086 SCIP_RETCODE SCIPincludeExprhdlrSum( 1087 SCIP* scip /**< SCIP data structure */ 1088 ) 1089 { 1090 SCIP_EXPRHDLR* exprhdlr; 1091 1092 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalSum, NULL) ); 1093 assert(exprhdlr != NULL); 1094 1095 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrSum, NULL); 1096 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataSum, freedataSum); 1097 SCIPexprhdlrSetSimplify(exprhdlr, simplifySum); 1098 SCIPexprhdlrSetCompare(exprhdlr, compareSum); 1099 SCIPexprhdlrSetPrint(exprhdlr, printSum); 1100 SCIPexprhdlrSetIntEval(exprhdlr, intevalSum); 1101 SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesSum, estimateSum); 1102 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropSum); 1103 SCIPexprhdlrSetHash(exprhdlr, hashSum); 1104 SCIPexprhdlrSetDiff(exprhdlr, bwdiffSum, fwdiffSum, bwfwdiffSum); 1105 SCIPexprhdlrSetCurvature(exprhdlr, curvatureSum); 1106 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicitySum); 1107 SCIPexprhdlrSetIntegrality(exprhdlr, integralitySum); 1108 SCIPexprhdlrSetGetSymdata(exprhdlr, getSymDataSum); 1109 1110 return SCIP_OKAY; 1111 } 1112 1113 /** creates a sum expression */ 1114 SCIP_RETCODE SCIPcreateExprSum( 1115 SCIP* scip, /**< SCIP data structure */ 1116 SCIP_EXPR** expr, /**< pointer where to store expression */ 1117 int nchildren, /**< number of children */ 1118 SCIP_EXPR** children, /**< children */ 1119 SCIP_Real* coefficients, /**< array with coefficients for all children (or NULL if all 1.0) */ 1120 SCIP_Real constant, /**< constant term of sum */ 1121 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 1122 void* ownercreatedata /**< data to pass to ownercreate */ 1123 ) 1124 { 1125 SCIP_EXPRDATA* exprdata; 1126 1127 SCIP_CALL( createData(scip, &exprdata, nchildren, coefficients, constant) ); 1128 1129 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrSum(scip), exprdata, nchildren, children, ownercreate, ownercreatedata) ); 1130 1131 return SCIP_OKAY; 1132 } 1133 1134 /** sets the constant of a summation expression */ 1135 void SCIPsetConstantExprSum( 1136 SCIP_EXPR* expr, /**< sum expression */ 1137 SCIP_Real constant /**< constant */ 1138 ) 1139 { 1140 SCIP_EXPRDATA* exprdata; 1141 1142 assert(expr != NULL); 1143 1144 exprdata = SCIPexprGetData(expr); 1145 assert(exprdata != NULL); 1146 1147 exprdata->constant = constant; 1148 } 1149 1150 /** appends an expression to a sum expression */ 1151 SCIP_RETCODE SCIPappendExprSumExpr( 1152 SCIP* scip, /**< SCIP data structure */ 1153 SCIP_EXPR* expr, /**< sum expression */ 1154 SCIP_EXPR* child, /**< expression to be appended */ 1155 SCIP_Real childcoef /**< child's coefficient */ 1156 ) 1157 { 1158 SCIP_EXPRDATA* exprdata; 1159 int nchildren; 1160 1161 assert(expr != NULL); 1162 assert(SCIPisExprSum(scip, expr)); 1163 1164 exprdata = SCIPexprGetData(expr); 1165 assert(exprdata != NULL); 1166 1167 nchildren = SCIPexprGetNChildren(expr); 1168 1169 SCIP_CALL( SCIPensureBlockMemoryArray(scip, &exprdata->coefficients, &exprdata->coefssize, nchildren + 1) ); 1170 1171 assert(exprdata->coefssize > nchildren); 1172 exprdata->coefficients[nchildren] = childcoef; 1173 1174 SCIP_CALL( SCIPappendExprChild(scip, expr, child) ); 1175 1176 return SCIP_OKAY; 1177 } 1178 1179 /** multiplies given sum expression by a constant */ 1180 void SCIPmultiplyByConstantExprSum( 1181 SCIP_EXPR* expr, /**< sum expression */ 1182 SCIP_Real constant /**< constant that multiplies sum expression */ 1183 ) 1184 { 1185 int i; 1186 SCIP_EXPRDATA* exprdata; 1187 1188 assert(expr != NULL); 1189 1190 exprdata = SCIPexprGetData(expr); 1191 assert(exprdata != NULL); 1192 1193 for( i = 0; i < SCIPexprGetNChildren(expr); ++i ) 1194 exprdata->coefficients[i] *= constant; 1195 exprdata->constant *= constant; 1196 } 1197 1198 /** constructs the expanded product of two sum expressions */ 1199 SCIP_RETCODE SCIPmultiplyBySumExprSum( 1200 SCIP* scip, /**< SCIP data structure */ 1201 SCIP_EXPR** product, /**< buffer where to store multiplied sums (expanded as sum) */ 1202 SCIP_EXPR* factor1, /**< first sum */ 1203 SCIP_EXPR* factor2, /**< second sum */ 1204 SCIP_Bool simplify, /**< whether to simplify created terms and sum */ 1205 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 1206 void* ownercreatedata /**< data to pass to ownercreate */ 1207 ) 1208 { 1209 SCIP_Real constant1; 1210 SCIP_Real constant2; 1211 int nchildren1; 1212 int nchildren2; 1213 int i1; 1214 int i2; 1215 1216 assert(scip != NULL); 1217 assert(product != NULL); 1218 assert(factor1 != NULL); 1219 assert(SCIPisExprSum(scip, factor1)); 1220 assert(factor2 != NULL); 1221 assert(SCIPisExprSum(scip, factor2)); 1222 1223 constant1 = SCIPgetConstantExprSum(factor1); 1224 constant2 = SCIPgetConstantExprSum(factor2); 1225 nchildren1 = SCIPexprGetNChildren(factor1); 1226 nchildren2 = SCIPexprGetNChildren(factor2); 1227 1228 /* TODO might be nice to integrate more with simplify and construct a simplified sum right away */ 1229 1230 SCIP_CALL( SCIPcreateExprSum(scip, product, 0, NULL, NULL, constant1 * constant2, ownercreate, ownercreatedata) ); 1231 1232 /* first add constant1 * factor2 1233 * constant * constant2 already added above 1234 */ 1235 if( constant1 != 0.0 ) 1236 { 1237 for( i2 = 0; i2 < nchildren2; ++i2 ) 1238 { 1239 SCIP_CALL( SCIPappendExprSumExpr(scip, *product, SCIPexprGetChildren(factor2)[i2], constant1 * SCIPgetCoefsExprSum(factor2)[i2]) ); 1240 } 1241 } 1242 1243 for( i1 = 0; i1 < nchildren1; ++i1 ) 1244 { 1245 SCIP_EXPR* child1; 1246 SCIP_EXPR* child2; 1247 SCIP_Real coef1; 1248 SCIP_Real coef2; 1249 1250 coef1 = SCIPgetCoefsExprSum(factor1)[i1]; 1251 child1 = SCIPexprGetChildren(factor1)[i1]; 1252 1253 if( constant2 != 0.0 ) 1254 { 1255 /* add coef1 * child1 * constant2 */ 1256 SCIP_CALL( SCIPappendExprSumExpr(scip, *product, child1, coef1 * constant2) ); 1257 } 1258 1259 for( i2 = 0; i2 < nchildren2; ++i2 ) 1260 { 1261 /* add coef1 * child1 * coef2 * child2 */ 1262 SCIP_EXPR* termprod; 1263 SCIP_EXPR* termprodsimplified; 1264 SCIP_EXPR* termfactors[2]; 1265 1266 coef2 = SCIPgetCoefsExprSum(factor2)[i2]; 1267 child2 = SCIPexprGetChildren(factor2)[i2]; 1268 1269 /* create child1 * child2 expr */ 1270 termfactors[0] = child1; 1271 termfactors[1] = child2; 1272 SCIP_CALL( SCIPcreateExprProduct(scip, &termprod, 2, termfactors, 1.0, NULL, NULL) ); 1273 1274 if( simplify ) 1275 { 1276 SCIP_Bool changed; 1277 SCIP_Bool infeasible; 1278 1279 /* simplify child1*child2 */ 1280 SCIP_CALL( SCIPsimplifyExpr(scip, termprod, &termprodsimplified, &changed, &infeasible, ownercreate, ownercreatedata) ); 1281 assert(!infeasible); /* simplify of products should never be infeasible */ 1282 SCIP_CALL( SCIPreleaseExpr(scip, &termprod) ); 1283 } 1284 else 1285 termprodsimplified = termprod; 1286 1287 /* append to product */ 1288 SCIP_CALL( SCIPappendExprSumExpr(scip, *product, termprodsimplified, coef1 * coef2) ); 1289 SCIP_CALL( SCIPreleaseExpr(scip, &termprodsimplified) ); 1290 } 1291 } 1292 1293 if( simplify ) 1294 { 1295 SCIP_EXPR* prodsimplified; 1296 SCIP_Bool changed; 1297 SCIP_Bool infeasible; 1298 1299 SCIP_CALL( SCIPsimplifyExpr(scip, *product, &prodsimplified, &changed, &infeasible, ownercreate, ownercreatedata) ); 1300 assert(!infeasible); 1301 SCIP_CALL( SCIPreleaseExpr(scip, product) ); 1302 *product = prodsimplified; 1303 } 1304 1305 return SCIP_OKAY; 1306 } 1307 1308 /** constructs the expanded power of a sum expression 1309 * 1310 * @attention The number of terms in the expansion grows exponential with the exponent. Be aware of what you wish for. 1311 */ 1312 SCIP_EXPORT 1313 SCIP_RETCODE SCIPpowerExprSum( 1314 SCIP* scip, /**< SCIP data structure */ 1315 SCIP_EXPR** result, /**< buffer where to store expanded power of sum */ 1316 SCIP_EXPR* base, /**< sum */ 1317 int exponent, /**< exponent > 1 */ 1318 SCIP_Bool simplify, /**< whether to simplify created terms and sum */ 1319 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 1320 void* ownercreatedata /**< data to pass to ownercreate */ 1321 ) 1322 { 1323 /* for sum_i alpha_i expr_i and exponent p, this constructs 1324 * sum_{beta} (p over beta) prod_i (alpha_i expr_i)^beta_i over all multiindex beta such that sum_i beta_i = p 1325 * See also https://en.wikipedia.org/wiki/Multinomial_theorem 1326 * Calculation of multinomials adapted from https://github.com/m-j-w/MultinomialSeries.jl 1327 * 1328 * TODO might be nice to integrate more with simplify and construct a simplified sum right away 1329 */ 1330 1331 SCIP_EXPR** children; 1332 SCIP_EXPR*** childrenpower; 1333 SCIP_Real* coefs; 1334 SCIP_Real constant; 1335 SCIP_Bool haveconst; 1336 int nchildren; 1337 int nterms; 1338 int* factorials; 1339 int* beta; 1340 int betapos; 1341 int restsum; 1342 int multinomialcoef; 1343 int i; 1344 1345 SCIP_EXPR* newterm; 1346 SCIP_Real newtermcoef; 1347 1348 assert(scip != NULL); 1349 assert(result != NULL); 1350 assert(base != NULL); 1351 assert(exponent > 1); 1352 assert(SCIPisExprSum(scip, base)); 1353 assert(exponent > 1.0); 1354 1355 children = SCIPexprGetChildren(base); 1356 nchildren = SCIPexprGetNChildren(base); 1357 coefs = SCIPgetCoefsExprSum(base); 1358 constant = SCIPgetConstantExprSum(base); 1359 haveconst = constant != 0.0; 1360 nterms = nchildren + (haveconst ? 1 : 0); 1361 1362 #ifdef SCIP_DEBUG 1363 SCIPinfoMessage(scip, NULL, "expanding ("); 1364 SCIPprintExpr(scip, base, NULL); 1365 SCIPinfoMessage(scip, NULL, ")^%d\n", exponent); 1366 #endif 1367 1368 SCIP_CALL( SCIPcreateExprSum(scip, result, 0, NULL, NULL, 0.0, ownercreate, ownercreatedata) ); 1369 1370 SCIP_CALL( SCIPallocClearBufferArray(scip, &beta, nterms) ); 1371 SCIP_CALL( SCIPallocBufferArray(scip, &factorials, exponent+1) ); 1372 1373 /* precompute factorials k!, k = 0...exponent */ 1374 factorials[0] = 1; 1375 for( i = 1; i <= exponent; ++i ) 1376 factorials[i] = factorials[i-1] * i; 1377 1378 /* precompute children^k, k=1..exponent */ 1379 SCIP_CALL( SCIPallocBufferArray(scip, &childrenpower, nchildren) ); 1380 for( i = 0; i < nchildren; ++i ) 1381 { 1382 int k; 1383 SCIP_CALL( SCIPallocBufferArray(scip, &childrenpower[i], exponent+1) ); 1384 childrenpower[i][1] = children[i]; 1385 for( k = 2; k <= exponent; ++k ) 1386 { 1387 SCIP_CALL( SCIPcreateExprPow(scip, &childrenpower[i][k], children[i], (SCIP_Real)k, NULL, NULL) ); 1388 if( simplify ) 1389 { 1390 SCIP_Bool changed; 1391 SCIP_Bool infeasible; 1392 SCIP_EXPR* simplified; 1393 1394 SCIP_CALL( SCIPsimplifyExpr(scip, childrenpower[i][k], &simplified, &changed, &infeasible, ownercreate, ownercreatedata) ); 1395 assert(!infeasible); 1396 SCIP_CALL( SCIPreleaseExpr(scip, &childrenpower[i][k]) ); 1397 childrenpower[i][k] = simplified; 1398 } 1399 } 1400 } 1401 1402 /* first multinomial is (exponent,0,0,...) */ 1403 beta[0] = exponent; 1404 betapos = 0; 1405 do 1406 { 1407 /* compute multinomial coef exponent! / (beta[0]! * ... * beta[nterms-1]!) */ 1408 multinomialcoef = factorials[exponent]; 1409 for( i = 0; i < nterms; ++i ) 1410 { 1411 assert(beta[i] >= 0); 1412 assert(beta[i] <= exponent); 1413 multinomialcoef /= factorials[beta[i]]; 1414 } 1415 1416 SCIPdebugMsg(scip, "multinomial ("); 1417 for( i = 0; i < nterms; ++i ) 1418 SCIPdebugPrintf("%d ", beta[i]); 1419 SCIPdebugPrintf(") with coef %d\n", multinomialcoef); 1420 1421 /* construct new term for expanded sum */ 1422 SCIP_CALL( SCIPcreateExprProduct(scip, &newterm, 0, NULL, 1.0, ownercreate, ownercreatedata) ); 1423 newtermcoef = multinomialcoef; 1424 for( i = 0; i < nterms; ++i ) 1425 { 1426 if( beta[i] == 0 ) 1427 continue; 1428 1429 if( i == nterms-1 && haveconst ) 1430 { 1431 /* if constant term, then update newtermcoef */ 1432 newtermcoef *= pow(constant, (double)beta[i]); 1433 continue; 1434 } 1435 1436 /* alpha_i^beta_i */ 1437 newtermcoef *= pow(coefs[i], (double)beta[i]); 1438 1439 /* expr_i^beta_i*/ 1440 SCIP_CALL( SCIPappendExprChild(scip, newterm, childrenpower[i][beta[i]]) ); 1441 } 1442 1443 /* append newterm to sum */ 1444 switch( SCIPexprGetNChildren(newterm) ) 1445 { 1446 case 0: 1447 { 1448 /* no factor in product, so it is a constant -> update constant in sum */ 1449 SCIPsetConstantExprSum(*result, SCIPgetConstantExprSum(*result) + newtermcoef); 1450 break; 1451 } 1452 1453 case 1: 1454 { 1455 /* only one factor in product -> add this factor itself to sum */ 1456 SCIP_CALL( SCIPappendExprSumExpr(scip, *result, SCIPexprGetChildren(newterm)[0], newtermcoef) ); 1457 break; 1458 } 1459 1460 default: 1461 { 1462 if( simplify ) 1463 { 1464 /* simplify product */ 1465 SCIP_Bool changed; 1466 SCIP_Bool infeasible; 1467 SCIP_EXPR* simplified; 1468 1469 SCIP_CALL( SCIPsimplifyExpr(scip, newterm, &simplified, &changed, &infeasible, ownercreate, ownercreatedata) ); 1470 assert(!infeasible); 1471 SCIP_CALL( SCIPreleaseExpr(scip, &newterm) ); 1472 newterm = simplified; 1473 } 1474 1475 /* append new term to sum */ 1476 SCIP_CALL( SCIPappendExprSumExpr(scip, *result, newterm, newtermcoef) ); 1477 break; 1478 } 1479 } 1480 SCIP_CALL( SCIPreleaseExpr(scip, &newterm) ); 1481 1482 /* determine next beta */ 1483 while( beta[betapos] == 0 ) 1484 --betapos; 1485 1486 if( betapos == nterms-1 ) 1487 { 1488 do 1489 { 1490 if( betapos == 0 ) 1491 goto TERMINATE; 1492 --betapos; 1493 } 1494 while( beta[betapos] == 0 ); 1495 1496 restsum = 0; 1497 for( i = betapos+1; i < nterms; ++i ) 1498 { 1499 restsum += beta[i]; 1500 beta[i] = 0; 1501 } 1502 beta[betapos+1] = restsum; 1503 } 1504 1505 if( beta[betapos] > 0 ) 1506 { 1507 --beta[betapos]; 1508 ++beta[++betapos]; 1509 } 1510 } 1511 while( TRUE ); /*lint !e506*/ 1512 1513 TERMINATE: 1514 1515 if( simplify ) 1516 { 1517 SCIP_Bool changed; 1518 SCIP_Bool infeasible; 1519 SCIP_EXPR* simplified; 1520 1521 SCIP_CALL( SCIPsimplifyExpr(scip, *result, &simplified, &changed, &infeasible, ownercreate, ownercreatedata) ); 1522 assert(!infeasible); 1523 SCIP_CALL( SCIPreleaseExpr(scip, result) ); 1524 *result = simplified; 1525 } 1526 1527 for( i = nchildren-1; i >= 0; --i ) 1528 { 1529 int k; 1530 for( k = exponent; k >= 2; --k ) 1531 { 1532 SCIP_CALL( SCIPreleaseExpr(scip, &childrenpower[i][k]) ); 1533 } 1534 SCIPfreeBufferArray(scip, &childrenpower[i]); 1535 } 1536 SCIPfreeBufferArray(scip, &childrenpower); 1537 SCIPfreeBufferArray(scip, &factorials); 1538 SCIPfreeBufferArray(scip, &beta); 1539 1540 #ifdef SCIP_DEBUG 1541 SCIPinfoMessage(scip, NULL, "-> "); 1542 SCIPprintExpr(scip, *result, NULL); 1543 SCIPinfoMessage(scip, NULL, "\n"); 1544 #endif 1545 1546 return SCIP_OKAY; 1547 } 1548 1549 /* from pub_expr.h */ 1550 1551 /** gets the coefficients of a summation expression */ 1552 SCIP_Real* SCIPgetCoefsExprSum( 1553 SCIP_EXPR* expr /**< sum expression */ 1554 ) 1555 { 1556 SCIP_EXPRDATA* exprdata; 1557 1558 assert(expr != NULL); 1559 1560 exprdata = SCIPexprGetData(expr); 1561 assert(exprdata != NULL); 1562 1563 return exprdata->coefficients; 1564 } 1565 1566 /** gets the constant of a summation expression */ 1567 SCIP_Real SCIPgetConstantExprSum( 1568 SCIP_EXPR* expr /**< sum expression */ 1569 ) 1570 { 1571 SCIP_EXPRDATA* exprdata; 1572 1573 assert(expr != NULL); 1574 1575 exprdata = SCIPexprGetData(expr); 1576 assert(exprdata != NULL); 1577 1578 return exprdata->constant; 1579 } 1580