1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file expr_sum.c
17 * @ingroup DEFPLUGINS_EXPR
18 * @brief sum expression handler
19 * @author Stefan Vigerske
20 * @author Benjamin Mueller
21 * @author Felipe Serrano
22 */
23
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25
26 #include <string.h>
27 #include <stddef.h>
28
29 #include "scip/expr_sum.h"
30 #include "scip/expr_value.h"
31 #include "scip/expr_product.h"
32 #include "scip/expr_exp.h"
33
34 #define EXPRHDLR_NAME "sum"
35 #define EXPRHDLR_DESC "summation with coefficients and a constant"
36 #define EXPRHDLR_PRECEDENCE 40000
37 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(47161.0)
38
39 /** macro to activate/deactivate debugging information of simplify method */
40 /*lint -emacro(774,debugSimplify) */
41 #ifdef SIMPLIFY_DEBUG
42 #define debugSimplify printf
43 #else
44 #define debugSimplify while( FALSE ) printf
45 #endif
46
47 /*
48 * Data structures
49 */
50
51 /** expression data */
52 struct SCIP_ExprData
53 {
54 SCIP_Real constant; /**< constant coefficient */
55 SCIP_Real* coefficients; /**< coefficients of children */
56 int coefssize; /**< size of the coefficients array */
57 };
58
59 /*
60 * Local methods
61 */
62
63 /** creates expression data */
64 static
65 SCIP_RETCODE createData(
66 SCIP* scip, /**< SCIP data structure */
67 SCIP_EXPRDATA** exprdata, /**< pointer where to store expression data */
68 int ncoefficients, /**< number of coefficients (i.e., number of children) */
69 SCIP_Real* coefficients, /**< array with coefficients for all children (or NULL if all 1.0) */
70 SCIP_Real constant /**< constant term of sum */
71 )
72 {
73 assert(exprdata != NULL);
74 assert(ncoefficients >= 0);
75
76 SCIP_CALL( SCIPallocBlockMemory(scip, exprdata) );
77
78 if( coefficients != NULL )
79 {
80 SCIP_CALL( SCIPduplicateBlockMemoryArray(scip, &(*exprdata)->coefficients, coefficients, ncoefficients) );
81 }
82 else
83 {
84 int i;
85
86 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(*exprdata)->coefficients, ncoefficients) );
87 for( i = 0; i < ncoefficients; ++i )
88 (*exprdata)->coefficients[i] = 1.0;
89 }
90
91 (*exprdata)->coefssize = ncoefficients;
92 (*exprdata)->constant = constant;
93
94 return SCIP_OKAY;
95 }
96
97 /** 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
98 *
99 * for example, this means that the `idx`-th child cannot be itself a sum
100 * if it is, we have to flatten it, i.e., take all its children and make them children of `duplicate`
101 */
102 static
103 SCIP_RETCODE simplifyTerm(
104 SCIP* scip, /**< SCIP data structure */
105 SCIP_EXPR* duplicate, /**< expression to be simplified */
106 int idx, /**< idx of children to be simplified */
107 SCIP_Bool* changed, /**< pointer to store if some term actually got simplified */
108 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
109 void* ownercreatedata /**< data to pass to ownercreate */
110 )
111 {
112 SCIP_EXPR** children;
113 SCIP_EXPR* expr;
114 SCIP_Real* coefs;
115 SCIP_Real constant ;
116 SCIP_Real coef;
117
118 assert(duplicate != NULL);
119 assert(idx >= 0);
120 assert(idx < SCIPexprGetNChildren(duplicate));
121 assert(changed != NULL);
122
123 children = SCIPexprGetChildren(duplicate);
124 coefs = SCIPgetCoefsExprSum(duplicate);
125 constant = SCIPgetConstantExprSum(duplicate);
126
127 coef = coefs[idx];
128 expr = children[idx];
129 assert(expr != NULL);
130
131 /* enforces SS3 */
132 if( SCIPisExprValue(scip, expr) )
133 {
134 *changed = TRUE;
135 constant += coef * SCIPgetValueExprValue(expr);
136 SCIPsetConstantExprSum(duplicate, constant);
137
138 /* TODO: remove child? */
139 coefs[idx] = 0.0;
140
141 return SCIP_OKAY;
142 }
143
144 /* enforces SS2 */
145 if( SCIPisExprSum(scip, expr) )
146 {
147 *changed = TRUE;
148
149 /* pass constant to parent */
150 constant += coef * SCIPgetConstantExprSum(expr);
151 SCIPsetConstantExprSum(duplicate, constant);
152
153 /* append all children of expr on parent except the first one */
154 if( SCIPexprGetNChildren(expr) > 1 )
155 {
156 int i;
157
158 for( i = 1; i < SCIPexprGetNChildren(expr); ++i )
159 {
160 assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[i]));
161 SCIP_CALL( SCIPappendExprSumExpr(scip, duplicate, SCIPexprGetChildren(expr)[i],
162 coef * SCIPgetCoefsExprSum(expr)[i]) );
163 }
164 }
165
166 /* replace expr with first child; need to get data again since it might be re-allocated */
167 assert(!SCIPisExprSum(scip, SCIPexprGetChildren(expr)[0]));
168
169 coefs = SCIPgetCoefsExprSum(duplicate);
170
171 coefs[idx] = coef * SCIPgetCoefsExprSum(expr)[0];
172 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, SCIPexprGetChildren(expr)[0]) );
173
174 return SCIP_OKAY;
175 }
176
177 /* enforce SS9 */
178 if( REALABS(coef) != 1.0 && SCIPisExprProduct(scip, expr) )
179 {
180 SCIP_EXPR* expchild = NULL;
181 int i;
182
183 for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
184 {
185 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
186 assert(child != NULL);
187
188 if( SCIPisExprExp(scip, child) )
189 {
190 expchild = child;
191 break;
192 }
193 }
194
195 /* coef != +- 1, term is product and one factor is an exponential -> enforce SS9 */
196 if( expchild != NULL )
197 {
198 SCIP_EXPR* sum;
199 SCIP_EXPR* prod;
200 SCIP_EXPR* simplifiedprod;
201 SCIP_EXPR* simplifiedsum;
202 SCIP_EXPR* exponential;
203 SCIP_EXPR* simplifiedexp;
204 SCIP_Real expconstant;
205
206 /* inform that expression will change */
207 *changed = TRUE;
208
209 /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */
210 if( coef > 0.0 )
211 {
212 expconstant = log(coef);
213 coefs[idx] = 1.0;
214 }
215 else
216 {
217 expconstant = log(-coef);
218 coefs[idx] = -1.0;
219 }
220
221 /* add constant to exponential's child */
222 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expchild), NULL, expconstant, ownercreate,
223 ownercreatedata) );
224
225 /* simplify sum */
226 SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) );
227 SCIP_CALL( SCIPreleaseExpr(scip, &sum) );
228
229 /* create exponential with new child */
230 SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) );
231 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) );
232
233 /* simplify exponential */
234 SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) );
235 SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
236
237 /* create product with new child */
238 SCIP_CALL( SCIPcreateExprProduct(scip, &prod, 0, NULL, 1.0, ownercreate, ownercreatedata) );
239
240 for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
241 {
242 if( SCIPexprGetChildren(expr)[i] == expchild )
243 {
244 SCIP_CALL( SCIPappendExprChild(scip, prod, simplifiedexp) );
245 }
246 else
247 {
248 SCIP_CALL( SCIPappendExprChild(scip, prod, SCIPexprGetChildren(expr)[i]) );
249 }
250 }
251 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) );
252
253 /* simplify product */
254 SCIP_CALL( SCIPcallExprSimplify(scip, prod, &simplifiedprod, ownercreate, ownercreatedata) );
255 SCIP_CALL( SCIPreleaseExpr(scip, &prod) );
256
257 /* replace current child with simplified product */
258 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedprod) );
259 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedprod) );
260
261 /* since the simplified product can be a sum ( exp(-1)*exp(log(x+y)+1) -> x+y ),
262 * we call the function we are in again
263 * this is no endless recursion, since the coef is now +- 1
264 */
265 SCIP_CALL( simplifyTerm(scip, duplicate, idx, changed, ownercreate, ownercreatedata) );
266
267 return SCIP_OKAY;
268 }
269 }
270
271 /* enforce SS10 */
272 if( REALABS(coef) != 1.0 && SCIPisExprExp(scip, expr) )
273 {
274 /* coef != +- 1, term is exponential -> enforce SS10 by moving |coef| into argument of exponential */
275
276 SCIP_EXPR* sum;
277 SCIP_EXPR* simplifiedsum;
278 SCIP_EXPR* exponential;
279 SCIP_EXPR* simplifiedexp;
280 SCIP_Real expconstant;
281
282 /* inform that expression will change */
283 *changed = TRUE;
284
285 /* compute expchild's coefficient as +- 1.0 * exp(log(abs(coef))) */
286 if( coef > 0.0 )
287 {
288 expconstant = log(coef);
289 coefs[idx] = 1.0;
290 }
291 else
292 {
293 expconstant = log(-coef);
294 coefs[idx] = -1.0;
295 }
296
297 /* add constant to exponential's child */
298 SCIP_CALL( SCIPcreateExprSum(scip, &sum, 1, SCIPexprGetChildren(expr), NULL, expconstant, ownercreate,
299 ownercreatedata) ); /* expconstant+expchild */
300
301 /* simplify sum */
302 SCIP_CALL( SCIPcallExprSimplify(scip, sum, &simplifiedsum, ownercreate, ownercreatedata) );
303 SCIP_CALL( SCIPreleaseExpr(scip, &sum) );
304
305 /* create exponential with new child */
306 SCIP_CALL( SCIPcreateExprExp(scip, &exponential, simplifiedsum, ownercreate, ownercreatedata) );
307 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedsum) );
308
309 /* simplify exponential */
310 SCIP_CALL( SCIPcallExprSimplify(scip, exponential, &simplifiedexp, ownercreate, ownercreatedata) );
311 SCIP_CALL( SCIPreleaseExpr(scip, &exponential) );
312
313 /* replace current child with simplified exponential */
314 SCIP_CALL( SCIPreplaceExprChild(scip, duplicate, idx, simplifiedexp) );
315 SCIP_CALL( SCIPreleaseExpr(scip, &simplifiedexp) );
316
317 return SCIP_OKAY;
318 }
319
320 /* other types of (simplified) expressions can be a child of a simplified sum */
321 assert(!SCIPisExprSum(scip, expr));
322 assert(!SCIPisExprValue(scip, expr));
323
324 return SCIP_OKAY;
325 }
326
327 /** helper struct for expressions sort */
328 typedef struct
329 {
330 SCIP* scip; /**< SCIP data structure */
331 SCIP_EXPR** exprs; /**< expressions */
332 } SORTEXPRDATA;
333
334 static
335 SCIP_DECL_SORTINDCOMP(sortExprComp)
336 {
337 SORTEXPRDATA* data = (SORTEXPRDATA*)dataptr;
338
339 return SCIPcompareExpr(data->scip, data->exprs[ind1], data->exprs[ind2]);
340 }
341
342 /*
343 * Callback methods of expression handler
344 */
345
346 /** simplifies a sum expression
347 *
348 * goes through each child and simplifies it; then sorts the simplified children; then sum the children that are equal;
349 * finally creates a sum expression with all the children that do not have a 0 coefficient and post-process so that SS6
350 * and SS7 are satisfied
351 */
352 static
353 SCIP_DECL_EXPRSIMPLIFY(simplifySum)
354 { /*lint --e{715}*/
355 SCIP_EXPR** children;
356 SCIP_EXPR* duplicate = NULL;
357 SCIP_EXPR** newchildren = NULL;
358 SCIP_Real* newcoefs = NULL;
359 int nnewchildren;
360 SCIP_Real newconstant;
361 SCIP_Real* coefs;
362 int i;
363 int nchildren;
364 SCIP_Bool changed;
365 SORTEXPRDATA sortdata;
366 int* order = NULL;
367
368 assert(expr != NULL);
369 assert(simplifiedexpr != NULL);
370 assert(SCIPexprGetHdlr(expr) == SCIPgetExprhdlrSum(scip));
371
372 changed = FALSE;
373
374 /* TODO: maybe have a flag to know if it is simplified ? */
375 /* TODO: can we do this with a shallow duplicate + copy of children pointer? currently simplifyTerm may modify children,
376 * so one would need to be careful
377 */
378 SCIP_CALL( SCIPduplicateExpr(scip, expr, &duplicate, NULL, NULL, ownercreate, ownercreatedata) );
379 assert(duplicate != NULL);
380
381 nchildren = SCIPexprGetNChildren(duplicate);
382 for( i = 0; i < nchildren; i++ )
383 {
384 /* enforces SS8 TODO: remove child? */
385 /* we have to ask for the coefs everytime, since it might get realloced in simpifyTerm */
386 if( SCIPgetCoefsExprSum(duplicate)[i] == 0.0 )
387 {
388 changed = TRUE;
389 continue;
390 }
391
392 /* enforces SS2, SS3, SS9, and SS10 */
393 SCIP_CALL( simplifyTerm(scip, duplicate, i, &changed, ownercreate, ownercreatedata) );
394 }
395
396 /* simplifyTerm can add new children to duplicate and realloc them; so get them again */
397 nchildren = SCIPexprGetNChildren(duplicate);
398 children = SCIPexprGetChildren(duplicate);
399 coefs = SCIPgetCoefsExprSum(duplicate);
400
401 /* treat zero term case */
402 if( nchildren == 0 )
403 {
404 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) );
405 goto CLEANUP;
406 }
407
408 /* treat one term case */
409 if( nchildren == 1 )
410 {
411 if( coefs[0] == 0.0 )
412 {
413 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, SCIPgetConstantExprSum(duplicate), ownercreate, ownercreatedata) );
414 goto CLEANUP;
415 }
416
417 if( coefs[0] == 1.0 && SCIPgetConstantExprSum(duplicate) == 0.0 )
418 *simplifiedexpr = children[0]; /* SS7 */
419 else
420 *simplifiedexpr = changed ? duplicate : expr;
421
422 SCIPcaptureExpr(*simplifiedexpr);
423
424 goto CLEANUP;
425 }
426
427 /* enforces SS5: sort children */
428 SCIP_CALL( SCIPallocBufferArray(scip, &order, nchildren) );
429 for( i = 0; i < nchildren; i++ )
430 order[i] = i;
431 sortdata.scip = scip;
432 sortdata.exprs = children;
433 SCIPsortInd(order, sortExprComp, (void*)&sortdata, nchildren);
434
435 /* create sorted variant of children and coefs */
436 SCIP_CALL( SCIPallocBufferArray(scip, &newchildren, nchildren) );
437 SCIP_CALL( SCIPallocBufferArray(scip, &newcoefs, nchildren) );
438 for( i = 0; i < nchildren; ++i )
439 {
440 newchildren[i] = children[order[i]];
441 newcoefs[i] = coefs[order[i]];
442 if( order[i] != i )
443 changed = TRUE;
444 }
445
446 /* post-process */
447
448 /* enforces SS4 */
449 nnewchildren = 0;
450 for( i = 0; i < nchildren; i++ )
451 {
452 /* eliminate zero-coefficients */
453 if( newcoefs[i] == 0.0 )
454 {
455 changed = TRUE;
456 continue;
457 }
458
459 /* sum equal expressions */
460 if( i < nchildren-1 && SCIPcompareExpr(scip, newchildren[i], newchildren[i+1]) == 0 )
461 {
462 changed = TRUE;
463 /* if we substract two almost equal not-so-small numbers, then set new coefficient to 0.0
464 * instead of some tiny value that is likely the result of some random round-off error
465 * E.g., on instance ex1221, we have x1^2 + b3 = 1.25.
466 * Probing finds an aggregation x1 = 1.11803 - 0.618034 b3.
467 * Simplify would then produce 1.25 + 1e-16 x1 = 1.25.
468 */
469 if( SCIPisEQ(scip, newcoefs[i], -newcoefs[i+1]) && REALABS(newcoefs[i]) >= 1.0 )
470 newcoefs[i+1] = 0.0;
471 else
472 newcoefs[i+1] += newcoefs[i];
473 continue;
474 }
475
476 /* move i-th child to new position */
477 newchildren[nnewchildren] = newchildren[i];
478 newcoefs[nnewchildren] = newcoefs[i];
479 nnewchildren++;
480 }
481
482 /* build sum expression from finalchildren and post-simplify */
483 newconstant = SCIPgetConstantExprSum(duplicate);
484
485 debugSimplify("what to do? finalchildren has length %d\n", nnewchildren); /*lint !e506 !e681*/
486
487 /* enforces SS6: if they are no children, return value */
488 if( nnewchildren == 0 )
489 {
490 debugSimplify("[sum] got empty list, return value %g\n", newconstant); /*lint !e506 !e681*/
491 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, newconstant, ownercreate, ownercreatedata) );
492
493 goto CLEANUP;
494 }
495
496 /* enforces SS7: if list consists of one expr with coef 1.0 and constant is 0, return that expr */
497 if( nnewchildren == 1 && newcoefs[0] == 1.0 && newconstant == 0.0 )
498 {
499 *simplifiedexpr = newchildren[0];
500 SCIPcaptureExpr(*simplifiedexpr);
501
502 goto CLEANUP;
503 }
504
505 /* build sum expression from children */
506 if( changed )
507 {
508 SCIP_CALL( SCIPcreateExprSum(scip, simplifiedexpr, nnewchildren, newchildren, newcoefs, newconstant,
509 ownercreate, ownercreatedata) );
510
511 goto CLEANUP;
512 }
513
514 *simplifiedexpr = expr;
515
516 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
517 SCIPcaptureExpr(*simplifiedexpr);
518
519 /* free memory */
520 CLEANUP:
521 SCIPfreeBufferArrayNull(scip, &newcoefs);
522 SCIPfreeBufferArrayNull(scip, &newchildren);
523 SCIPfreeBufferArrayNull(scip, &order);
524 SCIP_CALL( SCIPreleaseExpr(scip, &duplicate) );
525
526 return SCIP_OKAY;
527 }
528
529 /** compares two sum expressions.
530 *
531 * The order of two sum expressions is a lexicographical order on the terms.
532 *
533 * Starting from the *last*, we find the first child where they differ, say, the i-th.
534 * Then u < v <=> u_i < v_i.
535 * If there are no such children and they have different number of children, then u < v <=> nchildren(u) < nchildren(v).
536 * If there are no such children and they have the same number of children, then u < v <=> const(u) < const(v).
537 * Otherwise, they are the same.
538 *
539 * Note: we are assuming expression are simplified, so within u, we have u_1 < u_2, etc
540 *
541 * Example: y + z < x + y + z, 2*x + 3*y < 3*x + 3*y
542 */
543 static
544 SCIP_DECL_EXPRCOMPARE(compareSum)
545 { /*lint --e{715}*/
546 SCIP_Real const1;
547 SCIP_Real* coefs1;
548 SCIP_EXPR** children1;
549 int nchildren1;
550 SCIP_Real const2;
551 SCIP_Real* coefs2;
552 SCIP_EXPR** children2;
553 int nchildren2;
554 int compareresult;
555 int i;
556 int j;
557
558 nchildren1 = SCIPexprGetNChildren(expr1);
559 nchildren2 = SCIPexprGetNChildren(expr2);
560 children1 = SCIPexprGetChildren(expr1);
561 children2 = SCIPexprGetChildren(expr2);
562 coefs1 = SCIPgetCoefsExprSum(expr1);
563 coefs2 = SCIPgetCoefsExprSum(expr2);
564 const1 = SCIPgetConstantExprSum(expr1);
565 const2 = SCIPgetConstantExprSum(expr2);
566
567 for( i = nchildren1 - 1, j = nchildren2 - 1; i >= 0 && j >= 0; --i, --j )
568 {
569 compareresult = SCIPcompareExpr(scip, children1[i], children2[j]);
570 if( compareresult != 0 )
571 return compareresult;
572 else
573 {
574 /* expressions are equal, compare coefficient */
575 if( (coefs1 ? coefs1[i] : 1.0) < (coefs2 ? coefs2[j] : 1.0) )
576 return -1;
577 if( (coefs1 ? coefs1[i] : 1.0) > (coefs2 ? coefs2[j] : 1.0) )
578 return 1;
579
580 /* coefficients are equal, continue */
581 }
582 }
583
584 /* all children of one expression are children of the other expression, use number of children as a tie-breaker */
585 if( i < j )
586 {
587 assert(i == -1);
588 /* expr1 has less elements, hence expr1 < expr2 */
589 return -1;
590 }
591 if( i > j )
592 {
593 assert(j == -1);
594 /* expr1 has more elements, hence expr1 > expr2 */
595 return 1;
596 }
597
598 /* everything is equal, use constant/coefficient as tie-breaker */
599 assert(i == -1 && j == -1);
600 if( const1 < const2 )
601 return -1;
602 if( const1 > const2 )
603 return 1;
604
605 /* they are equal */
606 return 0;
607 }
608
609 /** expression handler copy callback */
610 static
611 SCIP_DECL_EXPRCOPYHDLR(copyhdlrSum)
612 { /*lint --e{715}*/
613 SCIP_CALL( SCIPincludeExprhdlrSum(scip) );
614
615 return SCIP_OKAY;
616 }
617
618 /** expression data copy callback */
619 static
620 SCIP_DECL_EXPRCOPYDATA(copydataSum)
621 { /*lint --e{715}*/
622 SCIP_EXPRDATA* sourceexprdata;
623
624 assert(targetexprdata != NULL);
625 assert(sourceexpr != NULL);
626
627 sourceexprdata = SCIPexprGetData(sourceexpr);
628 assert(sourceexprdata != NULL);
629
630 SCIP_CALL( createData(targetscip, targetexprdata, SCIPexprGetNChildren(sourceexpr),
631 sourceexprdata->coefficients, sourceexprdata->constant) );
632
633 return SCIP_OKAY;
634 }
635
636 /** expression data free callback */
637 static
638 SCIP_DECL_EXPRFREEDATA(freedataSum)
639 { /*lint --e{715}*/
640 SCIP_EXPRDATA* exprdata;
641
642 assert(expr != NULL);
643
644 exprdata = SCIPexprGetData(expr);
645 assert(exprdata != NULL);
646
647 SCIPfreeBlockMemoryArray(scip, &(exprdata->coefficients), exprdata->coefssize);
648 SCIPfreeBlockMemory(scip, &exprdata);
649
650 SCIPexprSetData(expr, NULL);
651
652 return SCIP_OKAY;
653 }
654
655 /** expression print callback */
656 static
657 SCIP_DECL_EXPRPRINT(printSum)
658 { /*lint --e{715}*/
659 SCIP_EXPRDATA* exprdata;
660
661 assert(expr != NULL);
662
663 exprdata = SCIPexprGetData(expr);
664 assert(exprdata != NULL);
665
666 switch( stage )
667 {
668 case SCIP_EXPRITER_ENTEREXPR :
669 {
670 /* print opening parenthesis, if necessary */
671 if( EXPRHDLR_PRECEDENCE <= parentprecedence )
672 {
673 SCIPinfoMessage(scip, file, "(");
674 }
675
676 /* print constant, if nonzero */
677 if( exprdata->constant != 0.0 )
678 {
679 SCIPinfoMessage(scip, file, "%g", exprdata->constant);
680 }
681 break;
682 }
683
684 case SCIP_EXPRITER_VISITINGCHILD :
685 {
686 SCIP_Real coef;
687
688 coef = exprdata->coefficients[currentchild];
689
690 /* print coefficient, if necessary */
691 if( coef == 1.0 )
692 {
693 /* if coefficient is 1.0, then print only "+" if not the first term */
694 if( exprdata->constant != 0.0 || currentchild > 0 )
695 {
696 SCIPinfoMessage(scip, file, "+");
697 }
698 }
699 else if( coef == -1.0 )
700 {
701 /* if coefficient is -1.0, then print only "-" */
702 SCIPinfoMessage(scip, file, "-");
703 }
704 else
705 {
706 /* force "+" sign on positive coefficient if not the first term */
707 SCIPinfoMessage(scip, file, (exprdata->constant != 0.0 || currentchild > 0) ? "%+g*" : "%g*", coef);
708 }
709
710 break;
711 }
712
713 case SCIP_EXPRITER_LEAVEEXPR :
714 {
715 /* print closing parenthesis, if necessary */
716 if( EXPRHDLR_PRECEDENCE <= parentprecedence )
717 {
718 SCIPinfoMessage(scip, file, ")");
719 }
720 break;
721 }
722
723 case SCIP_EXPRITER_VISITEDCHILD :
724 default: ;
725 }
726
727 return SCIP_OKAY;
728 }
729
730 /** expression point evaluation callback */
731 static
732 SCIP_DECL_EXPREVAL(evalSum)
733 { /*lint --e{715}*/
734 SCIP_EXPRDATA* exprdata;
735 int c;
736
737 assert(expr != NULL);
738
739 exprdata = SCIPexprGetData(expr);
740 assert(exprdata != NULL);
741
742 *val = exprdata->constant;
743 for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
744 {
745 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/
746
747 *val += exprdata->coefficients[c] * SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[c]);
748 }
749
750 return SCIP_OKAY;
751 }
752
753 /** expression forward derivative evaluation callback */
754 static
755 SCIP_DECL_EXPRFWDIFF(fwdiffSum)
756 { /*lint --e{715}*/
757 SCIP_EXPRDATA* exprdata;
758 int c;
759
760 assert(expr != NULL);
761 assert(dot != NULL);
762
763 exprdata = SCIPexprGetData(expr);
764 assert(exprdata != NULL);
765
766 *dot = 0.0;
767 for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
768 {
769 assert(SCIPexprGetDot(SCIPexprGetChildren(expr)[c]) != SCIP_INVALID); /*lint !e777*/
770
771 *dot += exprdata->coefficients[c] * SCIPexprGetDot(SCIPexprGetChildren(expr)[c]);
772 }
773
774 return SCIP_OKAY;
775 }
776
777 /** expression derivative evaluation callback */
778 static
779 SCIP_DECL_EXPRBWDIFF(bwdiffSum)
780 { /*lint --e{715}*/
781 assert(expr != NULL);
782 assert(SCIPexprGetData(expr) != NULL);
783 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr));
784 assert(SCIPexprGetChildren(expr)[childidx] != NULL);
785 assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[childidx]));
786
787 *val = SCIPgetCoefsExprSum(expr)[childidx];
788
789 return SCIP_OKAY;
790 }
791
792 /** expression backward forward derivative evaluation callback */
793 static
794 SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSum)
795 { /*lint --e{715}*/
796 assert(bardot != NULL);
797
798 *bardot = 0.0;
799
800 return SCIP_OKAY;
801 }
802
803 /** expression interval evaluation callback */
804 static
805 SCIP_DECL_EXPRINTEVAL(intevalSum)
806 { /*lint --e{715}*/
807 SCIP_EXPRDATA* exprdata;
808 SCIP_INTERVAL suminterval;
809 int c;
810
811 assert(expr != NULL);
812
813 exprdata = SCIPexprGetData(expr);
814 assert(exprdata != NULL);
815
816 SCIPintervalSet(interval, exprdata->constant);
817
818 SCIPdebugMsg(scip, "inteval %p with %d children: %.20g", (void*)expr, SCIPexprGetNChildren(expr), exprdata->constant);
819
820 for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
821 {
822 SCIP_INTERVAL childinterval;
823
824 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[c]);
825 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
826 {
827 SCIPintervalSetEmpty(interval);
828 break;
829 }
830
831 /* compute coefficients[c] * childinterval and add the result to the so far computed interval */
832 if( exprdata->coefficients[c] == 1.0 )
833 {
834 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, childinterval);
835 }
836 else
837 {
838 SCIPintervalMulScalar(SCIP_INTERVAL_INFINITY, &suminterval, childinterval, exprdata->coefficients[c]);
839 SCIPintervalAdd(SCIP_INTERVAL_INFINITY, interval, *interval, suminterval);
840 }
841
842 SCIPdebugMsgPrint(scip, " %+.20g*[%.20g,%.20g]", exprdata->coefficients[c], childinterval.inf, childinterval.sup);
843 }
844 SCIPdebugMsgPrint(scip, " = [%.20g,%.20g]\n", interval->inf, interval->sup);
845
846 return SCIP_OKAY;
847 }
848
849 /** initial estimators callback */
850 static
851 SCIP_DECL_EXPRINITESTIMATES(initEstimatesSum)
852 { /*lint --e{715}*/
853 SCIP_EXPRDATA* exprdata;
854
855 #ifdef SCIP_DEBUG
856 SCIPinfoMessage(scip, NULL, "initEstimatesSum %d children: ", SCIPexprGetNChildren(expr));
857 SCIPprintExpr(scip, expr, NULL);
858 SCIPinfoMessage(scip, NULL, "\n");
859 #endif
860 assert(scip != NULL);
861 assert(expr != NULL);
862 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
863
864 assert(coefs[0] != NULL);
865 assert(constant != NULL);
866 assert(nreturned != NULL);
867
868 exprdata = SCIPexprGetData(expr);
869 assert(exprdata != NULL);
870
871 BMScopyMemoryArray(coefs[0], exprdata->coefficients, SCIPexprGetNChildren(expr));
872 *constant = exprdata->constant;
873 *nreturned = 1;
874
875 return SCIP_OKAY;
876 }
877
878 /** expression estimate callback */
879 static
880 SCIP_DECL_EXPRESTIMATE(estimateSum)
881 { /*lint --e{715}*/
882 SCIP_EXPRDATA* exprdata;
883
884 assert(scip != NULL);
885 assert(expr != NULL);
886 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
887 assert(islocal != NULL);
888 assert(success != NULL);
889 assert(branchcand != NULL);
890
891 exprdata = SCIPexprGetData(expr);
892 assert(exprdata != NULL);
893
894 /* NOTE: nlhdlr_default assumes in nlhdlrInitSepaDefault that this estimator can be used for both under- and overestimation */
895
896 BMScopyMemoryArray(coefs, exprdata->coefficients, SCIPexprGetNChildren(expr));
897 *constant = exprdata->constant;
898 *islocal = FALSE;
899 *success = TRUE;
900
901 /* for none of our children, branching would improve the underestimator, so set branchcand[i]=FALSE everywhere
902 * if we branch for numerical reasons, then cons-expr-core should figure out what the candidates are
903 */
904 BMSclearMemoryArray(branchcand, SCIPexprGetNChildren(expr));
905
906 return SCIP_OKAY;
907 }
908
909 /** expression reverse propagation callback */
910 static
911 SCIP_DECL_EXPRREVERSEPROP(reversepropSum)
912 { /*lint --e{715}*/
913 SCIP_EXPRDATA* exprdata;
914 SCIP_INTERVAL* newbounds;
915 int nchildren;
916 int nreductions;
917
918 assert(scip != NULL);
919 assert(expr != NULL);
920 assert(infeasible != NULL);
921
922 nchildren = SCIPexprGetNChildren(expr);
923 assert(nchildren > 0);
924
925 exprdata = SCIPexprGetData(expr);
926 assert(exprdata != NULL);
927
928 SCIP_CALL( SCIPallocBufferArray(scip, &newbounds, nchildren) );
929
930 nreductions = SCIPintervalPropagateWeightedSum(SCIP_INTERVAL_INFINITY, nchildren, childrenbounds,
931 exprdata->coefficients, exprdata->constant, bounds, newbounds, infeasible);
932
933 if( !*infeasible && nreductions > 0 )
934 BMScopyMemoryArray(childrenbounds, newbounds, nchildren);
935
936 SCIPfreeBufferArray(scip, &newbounds);
937
938 return SCIP_OKAY;
939 }
940
941 /** sum hash callback */
942 static
943 SCIP_DECL_EXPRHASH(hashSum)
944 { /*lint --e{715}*/
945 SCIP_EXPRDATA* exprdata;
946 int c;
947
948 assert(scip != NULL);
949 assert(expr != NULL);
950 assert(hashkey != NULL);
951 assert(childrenhashes != NULL);
952
953 exprdata = SCIPexprGetData(expr);
954 assert(exprdata != NULL);
955
956 *hashkey = EXPRHDLR_HASHKEY;
957 *hashkey ^= SCIPcalcFibHash(exprdata->constant);
958
959 for( c = 0; c < SCIPexprGetNChildren(expr); ++c )
960 *hashkey ^= SCIPcalcFibHash(exprdata->coefficients[c]) ^ childrenhashes[c];
961
962 return SCIP_OKAY;
963 }
964
965 /** expression curvature detection callback */
966 static
967 SCIP_DECL_EXPRCURVATURE(curvatureSum)
968 { /*lint --e{715}*/
969 SCIP_EXPRDATA* exprdata;
970 int i;
971
972 assert(scip != NULL);
973 assert(expr != NULL);
974 assert(childcurv != NULL);
975 assert(success != NULL);
976
977 exprdata = SCIPexprGetData(expr);
978 assert(exprdata != NULL);
979
|
(1) Event cond_true: |
Condition "i < SCIPexprGetNChildren(expr)", taking true branch. |
980 for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
|
(2) Event ptr_arith: |
Performing pointer arithmetic on "childcurv" in expression "childcurv + i". |
981 childcurv[i] = SCIPexprcurvMultiply(exprdata->coefficients[i], exprcurvature);
982
983 *success = TRUE;
984
985 return SCIP_OKAY;
986 }
987
988 /** expression monotonicity detection callback */
989 static
990 SCIP_DECL_EXPRMONOTONICITY(monotonicitySum)
991 { /*lint --e{715}*/
992 SCIP_EXPRDATA* exprdata;
993
994 assert(scip != NULL);
995 assert(expr != NULL);
996 assert(result != NULL);
997 assert(childidx >= 0 && childidx < SCIPexprGetNChildren(expr));
998
999 exprdata = SCIPexprGetData(expr);
1000 assert(exprdata != NULL);
1001
1002 *result = exprdata->coefficients[childidx] >= 0.0 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
1003
1004 return SCIP_OKAY;
1005 }
1006
1007 /** expression integrality detection callback */
1008 static
1009 SCIP_DECL_EXPRINTEGRALITY(integralitySum)
1010 { /*lint --e{715}*/
1011 SCIP_EXPRDATA* exprdata;
1012 int i;
1013
1014 assert(scip != NULL);
1015 assert(expr != NULL);
1016 assert(isintegral != NULL);
1017
1018 exprdata = SCIPexprGetData(expr);
1019 assert(exprdata != NULL);
1020
1021 *isintegral = EPSISINT(exprdata->constant, 0.0); /*lint !e835*/
1022
1023 for( i = 0; i < SCIPexprGetNChildren(expr) && *isintegral; ++i )
1024 {
1025 SCIP_EXPR* child = SCIPexprGetChildren(expr)[i];
1026 assert(child != NULL);
1027
1028 *isintegral = EPSISINT(exprdata->coefficients[i], 0.0) && SCIPexprIsIntegral(child); /*lint !e835*/
1029 }
1030
1031 return SCIP_OKAY;
1032 }
1033
1034 /** creates the handler for sum expressions and includes it into SCIP */
1035 SCIP_RETCODE SCIPincludeExprhdlrSum(
1036 SCIP* scip /**< SCIP data structure */
1037 )
1038 {
1039 SCIP_EXPRHDLR* exprhdlr;
1040
1041 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalSum, NULL) );
1042 assert(exprhdlr != NULL);
1043
1044 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrSum, NULL);
1045 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataSum, freedataSum);
1046 SCIPexprhdlrSetSimplify(exprhdlr, simplifySum);
1047 SCIPexprhdlrSetCompare(exprhdlr, compareSum);
1048 SCIPexprhdlrSetPrint(exprhdlr, printSum);
1049 SCIPexprhdlrSetIntEval(exprhdlr, intevalSum);
1050 SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesSum, estimateSum);
1051 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropSum);
1052 SCIPexprhdlrSetHash(exprhdlr, hashSum);
1053 SCIPexprhdlrSetDiff(exprhdlr, bwdiffSum, fwdiffSum, bwfwdiffSum);
1054 SCIPexprhdlrSetCurvature(exprhdlr, curvatureSum);
1055 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicitySum);
1056 SCIPexprhdlrSetIntegrality(exprhdlr, integralitySum);
1057
1058 return SCIP_OKAY;
1059 }
1060
1061 /** creates a sum expression */
1062 SCIP_RETCODE SCIPcreateExprSum(
1063 SCIP* scip, /**< SCIP data structure */
1064 SCIP_EXPR** expr, /**< pointer where to store expression */
1065 int nchildren, /**< number of children */
1066 SCIP_EXPR** children, /**< children */
1067 SCIP_Real* coefficients, /**< array with coefficients for all children (or NULL if all 1.0) */
1068 SCIP_Real constant, /**< constant term of sum */
1069 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1070 void* ownercreatedata /**< data to pass to ownercreate */
1071 )
1072 {
1073 SCIP_EXPRDATA* exprdata;
1074
1075 SCIP_CALL( createData(scip, &exprdata, nchildren, coefficients, constant) );
1076
1077 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPgetExprhdlrSum(scip), exprdata, nchildren, children, ownercreate, ownercreatedata) );
1078
1079 return SCIP_OKAY;
1080 }
1081
1082 /** sets the constant of a summation expression */
1083 void SCIPsetConstantExprSum(
1084 SCIP_EXPR* expr, /**< sum expression */
1085 SCIP_Real constant /**< constant */
1086 )
1087 {
1088 SCIP_EXPRDATA* exprdata;
1089
1090 assert(expr != NULL);
1091
1092 exprdata = SCIPexprGetData(expr);
1093 assert(exprdata != NULL);
1094
1095 exprdata->constant = constant;
1096 }
1097
1098 /** appends an expression to a sum expression */
1099 SCIP_RETCODE SCIPappendExprSumExpr(
1100 SCIP* scip, /**< SCIP data structure */
1101 SCIP_EXPR* expr, /**< sum expression */
1102 SCIP_EXPR* child, /**< expression to be appended */
1103 SCIP_Real childcoef /**< child's coefficient */
1104 )
1105 {
1106 SCIP_EXPRDATA* exprdata;
1107 int nchildren;
1108
1109 assert(expr != NULL);
1110 assert(SCIPisExprSum(scip, expr));
1111
1112 exprdata = SCIPexprGetData(expr);
1113 assert(exprdata != NULL);
1114
1115 nchildren = SCIPexprGetNChildren(expr);
1116
1117 SCIP_CALL( SCIPensureBlockMemoryArray(scip, &exprdata->coefficients, &exprdata->coefssize, nchildren + 1) );
1118
1119 assert(exprdata->coefssize > nchildren);
1120 exprdata->coefficients[nchildren] = childcoef;
1121
1122 SCIP_CALL( SCIPappendExprChild(scip, expr, child) );
1123
1124 return SCIP_OKAY;
1125 }
1126
1127 /** multiplies given sum expression by a constant */
1128 void SCIPmultiplyByConstantExprSum(
1129 SCIP_EXPR* expr, /**< sum expression */
1130 SCIP_Real constant /**< constant that multiplies sum expression */
1131 )
1132 {
1133 int i;
1134 SCIP_EXPRDATA* exprdata;
1135
1136 assert(expr != NULL);
1137
1138 exprdata = SCIPexprGetData(expr);
1139 assert(exprdata != NULL);
1140
1141 for( i = 0; i < SCIPexprGetNChildren(expr); ++i )
1142 exprdata->coefficients[i] *= constant;
1143 exprdata->constant *= constant;
1144 }
1145
1146 /* from pub_expr.h */
1147
1148 /** gets the coefficients of a summation expression */
1149 SCIP_Real* SCIPgetCoefsExprSum(
1150 SCIP_EXPR* expr /**< sum expression */
1151 )
1152 {
1153 SCIP_EXPRDATA* exprdata;
1154
1155 assert(expr != NULL);
1156
1157 exprdata = SCIPexprGetData(expr);
1158 assert(exprdata != NULL);
1159
1160 return exprdata->coefficients;
1161 }
1162
1163 /** gets the constant of a summation expression */
1164 SCIP_Real SCIPgetConstantExprSum(
1165 SCIP_EXPR* expr /**< sum expression */
1166 )
1167 {
1168 SCIP_EXPRDATA* exprdata;
1169
1170 assert(expr != NULL);
1171
1172 exprdata = SCIPexprGetData(expr);
1173 assert(exprdata != NULL);
1174
1175 return exprdata->constant;
1176 }
1177