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