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_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 /* check for value expression */
173 if( SCIPisExprValue(scip, child) )
174 {
175 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, exp(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) );
176 }
177 else
178 {
179 *simplifiedexpr = expr;
180
181 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
182 SCIPcaptureExpr(*simplifiedexpr);
183 }
184
185 return SCIP_OKAY;
186 }
187
188 /** expression handler copy callback */
189 static
190 SCIP_DECL_EXPRCOPYHDLR(copyhdlrExp)
191 { /*lint --e{715}*/
192 SCIP_CALL( SCIPincludeExprhdlrExp(scip) );
193
194 return SCIP_OKAY;
195 }
196
197 /** expression data copy callback */
198 static
199 SCIP_DECL_EXPRCOPYDATA(copydataExp)
200 { /*lint --e{715}*/
201 assert(targetexprdata != NULL);
202 assert(sourceexpr != NULL);
203 assert(SCIPexprGetData(sourceexpr) == NULL);
204
205 *targetexprdata = NULL;
206
207 return SCIP_OKAY;
208 }
209
210 /** expression data free callback */
211 static
212 SCIP_DECL_EXPRFREEDATA(freedataExp)
213 { /*lint --e{715}*/
214 assert(expr != NULL);
215
216 SCIPexprSetData(expr, NULL);
217
218 return SCIP_OKAY;
219 }
220
221 /** expression parse callback */
222 static
223 SCIP_DECL_EXPRPARSE(parseExp)
224 { /*lint --e{715}*/
225 SCIP_EXPR* childexpr;
226
227 assert(expr != NULL);
228
229 /* parse child expression from remaining string */
230 SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
231 assert(childexpr != NULL);
232
233 /* create exponential expression */
234 SCIP_CALL( SCIPcreateExprExp(scip, expr, childexpr, ownercreate, ownercreatedata) );
235 assert(*expr != NULL);
236
237 /* release child expression since it has been captured by the exponential expression */
238 SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
239
240 *success = TRUE;
241
242 return SCIP_OKAY;
243 }
244
245 /** expression point evaluation callback */
246 static
247 SCIP_DECL_EXPREVAL(evalExp)
248 { /*lint --e{715}*/
249 assert(expr != NULL);
250 assert(SCIPexprGetData(expr) == NULL);
251 assert(SCIPexprGetNChildren(expr) == 1);
252 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
253
254 *val = exp(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
255
256 return SCIP_OKAY;
257 }
258
259 /** expression derivative evaluation callback */
260 static
261 SCIP_DECL_EXPRBWDIFF(bwdiffExp)
262 { /*lint --e{715}*/
263 assert(expr != NULL);
264 assert(childidx == 0);
265 assert(!SCIPisExprValue(scip, SCIPexprGetChildren(expr)[0]));
266 assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
267
268 *val = SCIPexprGetEvalValue(expr);
269
270 return SCIP_OKAY;
271 }
272
273 /** expression interval evaluation callback */
274 static
275 SCIP_DECL_EXPRINTEVAL(intevalExp)
276 { /*lint --e{715}*/
277 SCIP_INTERVAL childinterval;
278
279 assert(expr != NULL);
280 assert(SCIPexprGetData(expr) == NULL);
281 assert(SCIPexprGetNChildren(expr) == 1);
282
283 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
284
285 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
286 SCIPintervalSetEmpty(interval);
287 else
288 SCIPintervalExp(SCIP_INTERVAL_INFINITY, interval, childinterval);
289
290 return SCIP_OKAY;
291 }
292
293 /** expression estimator callback */
294 static
295 SCIP_DECL_EXPRESTIMATE(estimateExp)
296 { /*lint --e{715}*/
297 assert(scip != NULL);
298 assert(expr != NULL);
299 assert(SCIPexprGetNChildren(expr) == 1);
300 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
301 assert(coefs != NULL);
302 assert(constant != NULL);
303 assert(islocal != NULL);
304 assert(branchcand != NULL);
305 assert(*branchcand == TRUE);
306 assert(success != NULL);
307
308 *success = TRUE;
309 *coefs = 0.0;
310 *constant = 0.0;
311
312 if( overestimate )
313 {
314 addExpSecant(scip, localbounds[0].inf, localbounds[0].sup, coefs, constant, success);
315 *islocal = TRUE; /* secants are only valid locally */
316 }
317 else
318 {
319 addExpLinearization(scip, refpoint[0], SCIPexprIsIntegral(SCIPexprGetChildren(expr)[0]), coefs, constant, success);
320 *islocal = FALSE; /* linearization are globally valid */
321 *branchcand = FALSE;
322 }
323
324 return SCIP_OKAY;
325 }
326
327 /** initital estimates callback for an exponential expression */
328 static
329 SCIP_DECL_EXPRINITESTIMATES(initestimatesExp)
330 {
331 SCIP_Real refpointsunder[3] = {SCIP_INVALID, SCIP_INVALID, SCIP_INVALID};
332 SCIP_Bool overest[4] = {FALSE, FALSE, FALSE, TRUE};
333 SCIP_EXPR* child;
334 SCIP_Real lb;
335 SCIP_Real ub;
336 SCIP_Bool success;
337 int i;
338
339 assert(scip != NULL);
340 assert(expr != NULL);
341 assert(SCIPexprGetNChildren(expr) == 1);
342 assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
343
344 /* get expression data */
345 child = SCIPexprGetChildren(expr)[0];
346 assert(child != NULL);
347
348 lb = bounds[0].inf;
349 ub = bounds[0].sup;
350
|
(1) Event cond_true: |
Condition "!overestimate", taking true branch. |
351 if( !overestimate )
352 {
353 SCIP_Real lbfinite;
354 SCIP_Real ubfinite;
355
356 /* 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. |
357 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. |
358 ubfinite = SCIPisInfinity(scip, ub) ? MAX( 3.0, lb + 0.1 * REALABS(lb)) : ub; /*lint !e666*/
359
360 refpointsunder[0] = (7.0 * lbfinite + ubfinite) / 8.0;
361 refpointsunder[1] = (lbfinite + ubfinite) / 2.0;
362 refpointsunder[2] = (lbfinite + 7.0 * ubfinite) / 8.0;
363 }
364
365 *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] |
366 for( i = 0; i < 4; ++i )
367 {
|
(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. |
368 if( !overest[i] && overestimate )
|
(9) Event if_end: |
End of if statement. |
|
(22) Event if_end: |
End of if statement. |
369 continue;
370
|
(10) Event cond_false: |
Condition "overest[i]", taking false branch. |
|
(23) Event cond_false: |
Condition "overest[i]", taking false branch. |
371 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. |
372 continue;
373
374 assert(overest[i] || (SCIPisLE(scip, refpointsunder[i], ub) && SCIPisGE(scip, refpointsunder[i], lb))); /*lint !e661*/
375
376 coefs[*nreturned][0] = 0.0;
377 constant[*nreturned] = 0.0;
378
379 success = TRUE;
380
|
(12) Event cond_true: |
Condition "!overest[i]", taking true branch. |
|
(25) Event cond_true: |
Condition "!overest[i]", taking true branch. |
381 if( !overest[i] )
382 {
383 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] |
384 addExpLinearization(scip, refpointsunder[i], SCIPexprIsIntegral(child), coefs[*nreturned], &constant[*nreturned], &success); /*lint !e661*/
|
(13) Event if_fallthrough: |
Falling through to end of if statement. |
385 }
386 else
|
(14) Event if_end: |
End of if statement. |
387 addExpSecant(scip, lb, ub, coefs[*nreturned], &constant[*nreturned], &success);
388
|
(15) Event cond_true: |
Condition "success", taking true branch. |
389 if( success )
390 ++*nreturned;
|
(16) Event loop: |
Jumping back to the beginning of the loop. |
391 }
392
393 return SCIP_OKAY;
394 }
395
396 /** expression reverse propagation callback */
397 static
398 SCIP_DECL_EXPRREVERSEPROP(reversepropExp)
399 { /*lint --e{715}*/
400 assert(scip != NULL);
401 assert(expr != NULL);
402 assert(SCIPexprGetNChildren(expr) == 1);
403 assert(SCIPintervalGetInf(bounds) >= 0.0);
404
405 if( SCIPintervalGetSup(bounds) <= 0.0 )
406 {
407 *infeasible = TRUE;
408 return SCIP_OKAY;
409 }
410
411 /* f = exp(c0) -> c0 = log(f) */
412 SCIPintervalLog(SCIP_INTERVAL_INFINITY, childrenbounds, bounds);
413
414 return SCIP_OKAY;
415 }
416
417 /** expression hash callback */
418 static
419 SCIP_DECL_EXPRHASH(hashExp)
420 { /*lint --e{715}*/
421 assert(scip != NULL);
422 assert(expr != NULL);
423 assert(SCIPexprGetNChildren(expr) == 1);
424 assert(hashkey != NULL);
425 assert(childrenhashes != NULL);
426
427 *hashkey = EXPRHDLR_HASHKEY;
428 *hashkey ^= childrenhashes[0];
429
430 return SCIP_OKAY;
431 }
432
433 /** expression curvature detection callback */
434 static
435 SCIP_DECL_EXPRCURVATURE(curvatureExp)
436 { /*lint --e{715}*/
437 assert(scip != NULL);
438 assert(expr != NULL);
439 assert(childcurv != NULL);
440 assert(success != NULL);
441 assert(SCIPexprGetNChildren(expr) == 1);
442
443 /* expression is convex if child is convex; expression cannot be concave or linear */
444 if( exprcurvature == SCIP_EXPRCURV_CONVEX )
445 {
446 *success = TRUE;
447 *childcurv = SCIP_EXPRCURV_CONVEX;
448 }
449 else
450 *success = FALSE;
451
452 return SCIP_OKAY;
453 }
454
455 /** expression monotonicity detection callback */
456 static
457 SCIP_DECL_EXPRMONOTONICITY(monotonicityExp)
458 { /*lint --e{715}*/
459 assert(scip != NULL);
460 assert(expr != NULL);
461 assert(result != NULL);
462 assert(childidx == 0);
463
464 *result = SCIP_MONOTONE_INC;
465
466 return SCIP_OKAY;
467 }
468
469 /** creates the handler for exponential expressions and includes it into SCIP */
470 SCIP_RETCODE SCIPincludeExprhdlrExp(
471 SCIP* scip /**< SCIP data structure */
472 )
473 {
474 SCIP_EXPRHDLR* exprhdlr;
475
476 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC,
477 EXPRHDLR_PRECEDENCE, evalExp, NULL) );
478 assert(exprhdlr != NULL);
479
480 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrExp, NULL);
481 SCIPexprhdlrSetCopyFreeData(exprhdlr, copydataExp, freedataExp);
482 SCIPexprhdlrSetSimplify(exprhdlr, simplifyExp);
483 SCIPexprhdlrSetParse(exprhdlr, parseExp);
484 SCIPexprhdlrSetIntEval(exprhdlr, intevalExp);
485 SCIPexprhdlrSetEstimate(exprhdlr, initestimatesExp, estimateExp);
486 SCIPexprhdlrSetReverseProp(exprhdlr, reversepropExp);
487 SCIPexprhdlrSetHash(exprhdlr, hashExp);
488 SCIPexprhdlrSetDiff(exprhdlr, bwdiffExp, NULL, NULL);
489 SCIPexprhdlrSetCurvature(exprhdlr, curvatureExp);
490 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityExp);
491
492 return SCIP_OKAY;
493 }
494
495 /** creates an exponential expression */
496 SCIP_RETCODE SCIPcreateExprExp(
497 SCIP* scip, /**< SCIP data structure */
498 SCIP_EXPR** expr, /**< pointer where to store expression */
499 SCIP_EXPR* child, /**< single child */
500 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
501 void* ownercreatedata /**< data to pass to ownercreate */
502 )
503 {
504 assert(expr != NULL);
505 assert(child != NULL);
506 assert(SCIPfindExprhdlr(scip, EXPRHDLR_NAME) != NULL);
507
508 SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, EXPRHDLR_NAME), NULL, 1, &child, ownercreate, ownercreatedata) );
509
510 return SCIP_OKAY;
511 }
512
513 /** indicates whether expression is of exp-type */ /*lint -e{715}*/
514 SCIP_Bool SCIPisExprExp(
515 SCIP* scip, /**< SCIP data structure */
516 SCIP_EXPR* expr /**< expression */
517 )
518 { /*lint --e{715}*/
519 assert(expr != NULL);
520
521 return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0;
522 }
523