1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 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_exp.c
17 * @ingroup DEFPLUGINS_EXPR
18 * @brief exponential expression handler
19 * @author Stefan Vigerske
20 * @author Benjamin Mueller
21 * @author Ksenia Bestuzheva
22 */
23
24 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
25
26 #define _USE_MATH_DEFINES /* to get M_E on Windows */ /*lint !750 */
27
28 #include <string.h>
29 #include <math.h>
30
31 #include "scip/expr_exp.h"
32 #include "scip/expr_value.h"
33
34 #define EXPRHDLR_NAME "exp"
35 #define EXPRHDLR_DESC "exponential expression"
36 #define EXPRHDLR_PRECEDENCE 85000
37 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(10181.0)
38
39 /*
40 * Data structures
41 */
42
43 /*
44 * Local methods
45 */
46
47 /** computes coefficients of secant of an exponential term */
48 static
49 void addExpSecant(
50 SCIP* scip, /**< SCIP data structure */
51 SCIP_Real lb, /**< lower bound on variable */
52 SCIP_Real ub, /**< upper bound on variable */
53 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */
54 SCIP_Real* linconstant, /**< buffer to add constant of secant */
55 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
56 )
57 {
58 SCIP_Real coef;
59 SCIP_Real constant;
60
61 assert(scip != NULL);
62 assert(!SCIPisInfinity(scip, lb));
63 assert(!SCIPisInfinity(scip, -ub));
64 assert(SCIPisLE(scip, lb, ub));
65 assert(lincoef != NULL);
66 assert(linconstant != NULL);
67 assert(success != NULL);
68
69 if( SCIPisInfinity(scip, -lb) || SCIPisInfinity(scip, ub) )
70 {
71 /* unboundedness */
72 *success = FALSE;
73 return;
74 }
75
76 /* if lb and ub are too close use a safe secant */
77 if( SCIPisEQ(scip, lb, ub) )
78 {
79 coef = 0.0;
80 constant = exp(ub);
81 }
82 else
83 {
84 coef = (exp(ub) - exp(lb)) / (ub - lb);
85 constant = exp(ub) - coef * ub;
86 }
87
88 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
89 {
90 *success = FALSE;
91 return;
92 }
93
94 *lincoef += coef;
95 *linconstant += constant;
96 }
97
98 /** computes coefficients of linearization of an exponential term in a reference point */
99 static
100 void addExpLinearization(
101 SCIP* scip, /**< SCIP data structure */
102 SCIP_Real refpoint, /**< point for which to compute value of linearization */
103 SCIP_Bool isint, /**< whether corresponding variable is a discrete variable, and thus linearization could be moved */
104 SCIP_Real* lincoef, /**< buffer to add coefficient of secant */
105 SCIP_Real* linconstant, /**< buffer to add constant of secant */
106 SCIP_Bool* success /**< buffer to set to FALSE if secant has failed due to large numbers or unboundedness */
107 )
108 {
109 SCIP_Real constant;
110 SCIP_Real coef;
111
112 assert(scip != NULL);
113 assert(lincoef != NULL);
114 assert(linconstant != NULL);
115 assert(success != NULL);
116
117 if( SCIPisInfinity(scip, REALABS(refpoint)) )
118 {
119 *success = FALSE;
120 return;
121 }
122
123 if( !isint || SCIPisIntegral(scip, refpoint) )
124 {
125 coef = exp(refpoint);
126 constant = exp(refpoint) * (1.0 - refpoint);
127 }
128 else
129 {
130 /* exp(x) -> secant between f=floor(refpoint) and f+1 = ((e-1)*e^f) * x + e^f - f * ((e-1)*e^f) */
131 SCIP_Real f;
132
133 f = SCIPfloor(scip, refpoint);
134
135 coef = (M_E - 1.0) * exp(f);
136 constant = exp(f) - f * coef;
137 }
138
139 if( SCIPisInfinity(scip, REALABS(coef)) || SCIPisInfinity(scip, REALABS(constant)) )
140 {
141 *success = FALSE;
142 return;
143 }
144
145 *lincoef += coef;
146 *linconstant += constant;
147 }
148
149 /*
150 * Callback methods of expression handler
151 */
152
153 /** simplifies an exp expression
154 *
155 * Evaluates the exponential function when its child is a value expression.
156 *
157 * TODO: exp(log(*)) = *
158 */
159 static
160 SCIP_DECL_EXPRSIMPLIFY(simplifyExp)
161 {
162 SCIP_EXPR* child;
163
164 assert(scip != NULL);
165 assert(expr != NULL);
166 assert(simplifiedexpr != NULL);
167 assert(SCIPexprGetNChildren(expr) == 1);
168
169 child = SCIPexprGetChildren(expr)[0];
170 assert(child != NULL);
171
172 /**! [SnippetExprSimplifyExp] */
173 /* check for value expression */
174 if( SCIPisExprValue(scip, child) )
175 {
176 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, exp(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) );
177 }
178 else
179 {
180 *simplifiedexpr = expr;
181
182 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
183 SCIPcaptureExpr(*simplifiedexpr);
184 }
185 /**! [SnippetExprSimplifyExp] */
186
187 return SCIP_OKAY;
188 }
189
190 /** expression handler copy callback */
191 static
192 SCIP_DECL_EXPRCOPYHDLR(copyhdlrExp)
193 { /*lint --e{715}*/
194 SCIP_CALL( SCIPincludeExprhdlrExp(scip) );
195
196 return SCIP_OKAY;
197 }
198
199 /** expression data copy callback */
200 static
201 SCIP_DECL_EXPRCOPYDATA(copydataExp)
202 { /*lint --e{715}*/
203 assert(targetexprdata != NULL);
204 assert(sourceexpr != NULL);
205 assert(SCIPexprGetData(sourceexpr) == NULL);
206
207 *targetexprdata = NULL;
208
209 return SCIP_OKAY;
210 }
211
212 /** expression data free callback */
213 static
214 SCIP_DECL_EXPRFREEDATA(freedataExp)
215 { /*lint --e{715}*/
216 assert(expr != NULL);
217
218 SCIPexprSetData(expr, NULL);
219
220 return SCIP_OKAY;
221 }
222
223 /** expression parse callback */
224 static
225 SCIP_DECL_EXPRPARSE(parseExp)
226 { /*lint --e{715}*/
227 SCIP_EXPR* childexpr;
228
229 assert(expr != NULL);
230
231 /**! [SnippetExprParseExp] */
232 /* parse child expression from remaining string */
233 SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
234 assert(childexpr != NULL);
235
236 /* create exponential expression */
237 SCIP_CALL( SCIPcreateExprExp(scip, expr, childexpr, ownercreate, ownercreatedata) );
238 assert(*expr != NULL);
239
240 /* release child expression since it has been captured by the exponential expression */
241 SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
242
243 *success = TRUE;
244 /**! [SnippetExprParseExp] */
245
246 return SCIP_OKAY;
247 }
248
249 /** expression point evaluation callback */
250 static
251 SCIP_DECL_EXPREVAL(evalExp)
252 { /*lint --e{715}*/
253 assert(expr != NULL);
254 assert(SCIPexprGetData(expr) == NULL);
255 assert(SCIPexprGetNChildren(expr) == 1);
256 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
257
258 *val = exp(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
259
260 return SCIP_OKAY;
261 }
262
263 /** expression derivative evaluation callback */
264 static
265 SCIP_DECL_EXPRBWDIFF(bwdiffExp)
266 { /*lint --e{715}*/
267 assert(expr != NULL);
268 assert(childidx == 0);
269 assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[0]));
270 assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
271
272 *val = SCIPexprGetEvalValue(expr);
273
274 return SCIP_OKAY;
275 }
276
277 /** expression interval evaluation callback */
278 static
279 SCIP_DECL_EXPRINTEVAL(intevalExp)
280 { /*lint --e{715}*/
281 SCIP_INTERVAL childinterval;
282
283 assert(expr != NULL);
284 assert(SCIPexprGetData(expr) == NULL);
285 assert(SCIPexprGetNChildren(expr) == 1);
286
287 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
288
289 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
290 SCIPintervalSetEmpty(interval);
291 else
292 SCIPintervalExp(SCIP_INTERVAL_INFINITY, interval, childinterval);
293
294 return SCIP_OKAY;
295 }
296
297 /** expression estimator callback */
298 static
299 SCIP_DECL_EXPRESTIMATE(estimateExp)
300 { /*lint --e{715}*/
301 assert(scip != NULL);
302 assert(expr != NULL);
303 assert(SCIPexprGetNChildren(expr) == 1);
304 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
305 assert(coefs != NULL);
306 assert(constant != NULL);
307 assert(islocal != NULL);
308 assert(branchcand != NULL);
309 assert(*branchcand == TRUE);
310 assert(success != NULL);
311
312 *success = TRUE;
313 *coefs = 0.0;
314 *constant = 0.0;
315
316 if( overestimate )
317 {
318 addExpSecant(scip, localbounds[0].inf, localbounds[0].sup, coefs, constant, success);
319 *islocal = TRUE; /* secants are only valid locally */
320 }
321 else
322 {
323 addExpLinearization(scip, refpoint[0], SCIPexprIsIntegral(SCIPexprGetChildren(expr)[0]), coefs, constant, success);
324 *islocal = FALSE; /* linearization are globally valid */
325 *branchcand = FALSE;
326 }
327
328 return SCIP_OKAY;
329 }
330
331 /** initital estimates callback for an exponential expression */
332 static
333 SCIP_DECL_EXPRINITESTIMATES(initestimatesExp)
334 {
335 SCIP_Real refpointsunder[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
336 SCIP_Bool overest[4] = {FALSE, FALSE, FALSE, TRUE};
337 SCIP_EXPR* child;
338 SCIP_Real lb;
339 SCIP_Real ub;
340 SCIP_Bool success;
341 int i;
342
343 assert(scip != NULL);
344 assert(expr != NULL);
345 assert(SCIPexprGetNChildren(expr) == 1);
346 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
347
348 /* get expression data */
349 child = SCIPexprGetChildren(expr)[0];
350 assert(child != NULL);
351
352 lb = bounds[0].inf;
353 ub = bounds[0].sup;
354
|
(1) Event cond_true: |
Condition "!overestimate", taking true branch. |
355 if( !overestimate )
356 {
357 SCIP_Real lbfinite;
358 SCIP_Real ubfinite;
359
360 /* make bounds finite */
|
(2) Event cond_true: |
Condition "-lb >= scip->set->num_infinity", taking true branch. |
|
(3) Event cond_true: |
Condition "-5. <= ub - 0.1 * fabs(ub)", taking true branch. |
361 lbfinite = SCIPisInfinity(scip, -lb) ? MIN(-5.0, ub - 0.1 * REALABS(ub)) : lb; /*lint !e666*/
|
(4) Event cond_true: |
Condition "ub >= scip->set->num_infinity", taking true branch. |
|
(5) Event cond_true: |
Condition "3. >= lb + 0.1 * fabs(lb)", taking true branch. |
362 ubfinite = SCIPisInfinity(scip, ub) ? MAX( 3.0, lb + 0.1 * REALABS(lb)) : ub; /*lint !e666*/
363
364 refpointsunder[0] = (7.0 * lbfinite + ubfinite) / 8.0;
365 refpointsunder[1] = (lbfinite + ubfinite) / 2.0;
366 refpointsunder[2] = (lbfinite + 7.0 * ubfinite) / 8.0;
367 }
368
369 *nreturned = 0;
|
(6) Event cond_true: |
Condition "i < 4", taking true branch. |
|
(17) Event loop_begin: |
Jumped back to beginning of loop. |
|
(18) Event cond_true: |
Condition "i < 4", taking true branch. |
|
(19) Event cond_at_most: |
Checking "i < 4" implies that "i" may be up to 3 on the true branch. |
| Also see events: |
[overrun-local] |
370 for( i = 0; i < 4; ++i )
371 {
|
(7) Event cond_true: |
Condition "!overest[i]", taking true branch. |
|
(8) Event cond_false: |
Condition "overestimate", taking false branch. |
|
(20) Event cond_true: |
Condition "!overest[i]", taking true branch. |
|
(21) Event cond_false: |
Condition "overestimate", taking false branch. |
372 if( !overest[i] && overestimate )
|
(9) Event if_end: |
End of if statement. |
|
(22) Event if_end: |
End of if statement. |
373 continue;
374
|
(10) Event cond_false: |
Condition "overest[i]", taking false branch. |
|
(23) Event cond_false: |
Condition "overest[i]", taking false branch. |
375 if( overest[i] && (!overestimate || SCIPisInfinity(scip, ub) || SCIPisInfinity(scip, -lb)) )
|
(11) Event if_end: |
End of if statement. |
|
(24) Event if_end: |
End of if statement. |
376 continue;
377
378 assert(overest[i] || (SCIPisLE(scip, refpointsunder[i], ub) && SCIPisGE(scip, refpointsunder[i], lb))); /*lint !e661*/
379
380 coefs[*nreturned][0] = 0.0;
381 constant[*nreturned] = 0.0;
382
383 success = TRUE;
384
|
(12) Event cond_true: |
Condition "!overest[i]", taking true branch. |
|
(25) Event cond_true: |
Condition "!overest[i]", taking true branch. |
385 if( !overest[i] )
386 {
387 assert(i < 3);
|
(26) Event overrun-local: |
Overrunning array "refpointsunder" of 3 8-byte elements at element index 3 (byte offset 31) using index "i" (which evaluates to 3). |
| Also see events: |
[cond_at_most] |
388 addExpLinearization(scip, refpointsunder[i], SCIPexprIsIntegral(child), coefs[*nreturned], &constant[*nreturned], &success); /*lint !e661*/
|
(13) Event if_fallthrough: |
Falling through to end of if statement. |
389 }
390 else
|
(14) Event if_end: |
End of if statement. |
391 addExpSecant(scip, lb, ub, coefs[*nreturned], &constant[*nreturned], &success);
392
|
(15) Event cond_true: |
Condition "success", taking true branch. |
393 if( success )
394 ++*nreturned;
|
(16) Event loop: |
Jumping back to the beginning of the loop. |
395 }
396
397 return SCIP_OKAY;
398 }
399
400 /** expression reverse propagation callback */
401 static
402 SCIP_DECL_EXPRREVERSEPROP(reversepropExp)
403 { /*lint --e{715}*/
404 assert(scip != NULL);
405 assert(expr != NULL);
406 assert(SCIPexprGetNChildren(expr) == 1);
407 assert(SCIPintervalGetInf(bounds) >= 0.0);
408
409 if( SCIPintervalGetSup(bounds) <= 0.0 )
410 {
411 *infeasible = TRUE;
412 return SCIP_OKAY;
413 }
414
415 /* f = exp(c0) -> c0 = log(f) */
416 SCIPintervalLog(SCIP_INTERVAL_INFINITY, childrenbounds, bounds);
417
418 return SCIP_OKAY;
419 }
420
421 /** expression hash callback */
422 static
423 SCIP_DECL_EXPRHASH(hashExp)
424 { /*lint --e{715}*/
425 assert(scip != NULL);
426 assert(expr != NULL);
427 assert(SCIPexprGetNChildren(expr) == 1);
428 assert(hashkey != NULL);
429 assert(childrenhashes != NULL);
430
431 *hashkey = EXPRHDLR_HASHKEY;
432 *hashkey ^= childrenhashes[0];
433
434 return SCIP_OKAY;
435 }
436
437 /** expression curvature detection callback */
438 static
439 SCIP_DECL_EXPRCURVATURE(curvatureExp)
440 { /*lint --e{715}*/
441 assert(scip != NULL);
442 assert(expr != NULL);
443 assert(childcurv != NULL);
444 assert(success != NULL);
445 assert(SCIPexprGetNChildren(expr) == 1);
446
447 /* expression is convex if child is convex; expression cannot be concave or linear */
448 if( exprcurvature == SCIP_EXPRCURV_CONVEX )
449 {
450 *success = TRUE;
451 *childcurv = SCIP_EXPRCURV_CONVEX;
452 }
453 else
454 *success = FALSE;
455
456 return SCIP_OKAY;
457 }
458
459 /** expression monotonicity detection callback */
460 static
461 SCIP_DECL_EXPRMONOTONICITY(monotonicityExp)
462 { /*lint --e{715}*/
463 assert(scip != NULL);
464 assert(expr != NULL);
465 assert(result != NULL);
466 assert(childidx == 0);
467
468 *result = SCIP_MONOTONE_INC;
469
470 return SCIP_OKAY;
471 }
472
473 /** creates the handler for exponential expressions and includes it into SCIP */
474 SCIP_RETCODE SCIPincludeExprhdlrExp(
475 SCIP* scip /**< SCIP data structure */
476 )
477 {
478 SCIP_EXPRHDLR* exprhdlr;
479
480 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC,
481 EXPRHDLR_PRECEDENCE, evalExp, NULL) );
482 assert(exprhdlr != NULL);
483
484 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrExp, NULL);
485 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataExp, freedataExp);
486 SCIPexprhdlrSetSimplify(exprhdlr, simplifyExp);
487 SCIPexprhdlrSetParse(exprhdlr, parseExp);
488 SCIPexprhdlrSetIntEval(exprhdlr, intevalExp);
489 SCIPexprhdlrSetEstimate(exprhdlr, initestimatesExp, estimateExp);
490 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropExp);
491 SCIPexprhdlrSetHash(exprhdlr, hashExp);
492 SCIPexprhdlrSetDiff(exprhdlr, bwdiffExp, NULL, NULL);
493 SCIPexprhdlrSetCurvature(exprhdlr, curvatureExp);
494 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityExp);
495
496 return SCIP_OKAY;
497 }
498
499 /** creates an exponential expression */
500 SCIP_RETCODE SCIPcreateExprExp(
501 SCIP* scip, /**< SCIP data structure */
502 SCIP_EXPR** expr, /**< pointer where to store expression */
503 SCIP_EXPR* child, /**< single child */
504 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
505 void* ownercreatedata /**< data to pass to ownercreate */
506 )
507 {
508 assert(expr != NULL);
509 assert(child != NULL);
510 assert(SCIPfindExprhdlr(scip, EXPRHDLR_NAME) != NULL);
511
512 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, EXPRHDLR_NAME), NULL, 1, &child, ownercreate, ownercreatedata) );
513
514 return SCIP_OKAY;
515 }
516
517 /** indicates whether expression is of exp-type */ /*lint -e{715}*/
518 SCIP_Bool SCIPisExprExp(
519 SCIP* scip, /**< SCIP data structure */
520 SCIP_EXPR* expr /**< expression */
521 )
522 { /*lint --e{715}*/
523 assert(expr != NULL);
524
525 return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0;
526 }
527