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_erf.c
26   	 * @brief  handler for Gaussian error function expressions
27   	 * @author Benjamin Mueller
28   	 */
29   	
30   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
31   	
32   	#include "scip/expr_erf.h"
33   	#include "scip/expr_value.h"
34   	
35   	/* fundamental expression handler properties */
36   	#define EXPRHDLR_NAME         "erf"
37   	#define EXPRHDLR_DESC         "Gaussian error function"
38   	#define EXPRHDLR_PRECEDENCE   79000
39   	#define EXPRHDLR_HASHKEY      SCIPcalcFibHash(131071.0)
40   	
41   	/*
42   	 * Data structures
43   	 */
44   	
45   	/*
46   	 * Local methods
47   	 */
48   	
49   	/** evaluates the Gaussian error function at a given point */
50   	static
51   	SCIP_Real errorf(
52   	   SCIP_Real             x                   /**< point to evaluate */
53   	   )
54   	{
55   	   SCIP_Real a1 = +0.254829592;
56   	   SCIP_Real a2 = -0.284496736;
57   	   SCIP_Real a3 = +1.421413741;
58   	   SCIP_Real a4 = -1.453152027;
59   	   SCIP_Real a5 = +1.061405429;
60   	   SCIP_Real p  = +0.3275911;
61   	   int sign  = (x >= 0.0) ? 1 : -1;
62   	   SCIP_Real t = 1.0 / (1.0 + p * REALABS(x));
63   	   SCIP_Real y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * exp(-x*x);
64   	
65   	   return sign*y;
66   	}
67   	
68   	/*
69   	 * Callback methods of expression handler
70   	 */
71   	
72   	/** expression handler copy callback */
73   	static
74   	SCIP_DECL_EXPRCOPYHDLR(copyhdlrErf)
75   	{  /*lint --e{715}*/
76   	   SCIP_CALL( SCIPincludeExprhdlrErf(scip) );
77   	
78   	   return SCIP_OKAY;
79   	}
80   	
81   	/** simplifies an erf expression */
82   	static
83   	SCIP_DECL_EXPRSIMPLIFY(simplifyErf)
84   	{  /*lint --e{715}*/
85   	   SCIP_EXPR* child;
86   	
87   	   assert(scip != NULL);
88   	   assert(expr != NULL);
89   	   assert(simplifiedexpr != NULL);
90   	   assert(SCIPexprGetNChildren(expr) == 1);
91   	
92   	   child = SCIPexprGetChildren(expr)[0];
93   	   assert(child != NULL);
94   	
95   	   /* check for value expression */
96   	   if( SCIPisExprValue(scip, child) )
97   	   {
98   	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, errorf(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) );
99   	   }
100  	   else
101  	   {
102  	      *simplifiedexpr = expr;
103  	
104  	      /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
105  	      SCIPcaptureExpr(*simplifiedexpr);
106  	   }
107  	
108  	   return SCIP_OKAY;
109  	}
110  	
111  	/** expression parse callback */
112  	static
113  	SCIP_DECL_EXPRPARSE(parseErf)
114  	{  /*lint --e{715}*/
115  	   SCIP_EXPR* childexpr;
116  	
117  	   assert(expr != NULL);
118  	
119  	   /* parse child expression from remaining string */
120  	   SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
121  	   assert(childexpr != NULL);
122  	
123  	   /* create gaussian error function expression */
124  	   SCIP_CALL( SCIPcreateExprErf(scip, expr, childexpr, ownercreate, ownercreatedata) );
125  	   assert(*expr != NULL);
126  	
127  	   /* release child expression since it has been captured by the gaussian error function expression */
128  	   SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
129  	
130  	   *success = TRUE;
131  	
132  	   return SCIP_OKAY;
133  	}
134  	
135  	/** expression (point-) evaluation callback */
136  	static
137  	SCIP_DECL_EXPREVAL(evalErf)
138  	{  /*lint --e{715}*/
139  	   assert(expr != NULL);
140  	   assert(SCIPexprGetData(expr) == NULL);
141  	   assert(SCIPexprGetNChildren(expr) == 1);
142  	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
143  	
144  	   *val = errorf(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
145  	
146  	   return SCIP_OKAY;
147  	}
148  	
149  	/** expression derivative evaluation callback */
150  	static
151  	SCIP_DECL_EXPRBWDIFF(bwdiffErf)
152  	{  /*lint --e{715}*/
153  	   assert(expr != NULL);
154  	
155  	   SCIPerrorMessage("method of erf expression handler not implemented yet\n");
156  	   SCIPABORT(); /*lint --e{527}*/
157  	
158  	   return SCIP_OKAY;
159  	}
160  	
161  	/** expression interval evaluation callback */
162  	static
163  	SCIP_DECL_EXPRINTEVAL(intevalErf)
164  	{  /*lint --e{715}*/
165  	   SCIP_INTERVAL childinterval;
166  	
167  	   assert(expr != NULL);
168  	   assert(SCIPexprGetData(expr) == NULL);
169  	   assert(SCIPexprGetNChildren(expr) == 1);
170  	
171  	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
172  	
173  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
174  	      SCIPintervalSetEmpty(interval);
175  	   else
176  	   {
177  	      SCIP_Real childinf = SCIPintervalGetInf(childinterval);
178  	      SCIP_Real childsup = SCIPintervalGetSup(childinterval);
179  	      SCIP_Real inf = childinf <= -SCIP_INTERVAL_INFINITY ? -1.0 : errorf(childinf);
180  	      SCIP_Real sup = childsup >= +SCIP_INTERVAL_INFINITY ? +1.0 : errorf(childsup);
181  	      assert(inf <= sup);
182  	      SCIPintervalSetBounds(interval, inf, sup);
183  	   }
184  	
185  	   return SCIP_OKAY;
186  	}
187  	
188  	/** erf hash callback */
189  	static
190  	SCIP_DECL_EXPRHASH(hashErf)
191  	{  /*lint --e{715}*/
192  	   assert(scip != NULL);
193  	   assert(expr != NULL);
194  	   assert(SCIPexprGetNChildren(expr) == 1);
195  	   assert(hashkey != NULL);
196  	   assert(childrenhashes != NULL);
197  	
198  	   *hashkey = EXPRHDLR_HASHKEY;
199  	   *hashkey ^= childrenhashes[0];
200  	
201  	   return SCIP_OKAY;
202  	}
203  	
204  	/** expression curvature detection callback */
205  	static
206  	SCIP_DECL_EXPRCURVATURE(curvatureErf)
207  	{  /*lint --e{715}*/
208  	   assert(scip != NULL);
209  	   assert(expr != NULL);
210  	
211  	   /* expression is
212  	    *  - convex if child is convex and child <= 0
213  	    *  - concave if child is concave and child >= 0
214  	    */
215  	   if( exprcurvature == SCIP_EXPRCURV_CONVEX )
216  	   {
217  	      *success = TRUE;
218  	      *childcurv = SCIP_EXPRCURV_CONVEX;
219  	   }
220  	   else
221  	      *success = FALSE;
222  	
223  	   return SCIP_OKAY;
224  	}
225  	
226  	/** expression monotonicity detection callback */
227  	static
228  	SCIP_DECL_EXPRMONOTONICITY(monotonicityErf)
229  	{  /*lint --e{715}*/
230  	   assert(scip != NULL);
231  	   assert(expr != NULL);
232  	   assert(result != NULL);
233  	
234  	   *result = SCIP_MONOTONE_INC;
235  	
236  	   return SCIP_OKAY;
237  	}
238  	
239  	/** expression integrality detection callback */
240  	static
241  	SCIP_DECL_EXPRINTEGRALITY(integralityErf)
242  	{  /*lint --e{715}*/
243  	   assert(scip != NULL);
244  	   assert(expr != NULL);
245  	   assert(isintegral != NULL);
246  	
247  	   *isintegral = FALSE;
248  	
249  	   return SCIP_OKAY;
250  	}
251  	
252  	/** creates an erf expression
253  	 *
254  	 * @attention The implementation of `erf` expressions is incomplete.
255  	 * They are not usable for most use cases so far.
256  	 */
257  	SCIP_RETCODE SCIPcreateExprErf(
258  	   SCIP*                 scip,               /**< SCIP data structure */
259  	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
260  	   SCIP_EXPR*            child,              /**< child expression */
261  	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
262  	   void*                 ownercreatedata     /**< data to pass to ownercreate */
263  	   )
264  	{
265  	   SCIP_EXPRHDLR* exprhdlr;
266  	
267  	   assert(expr != NULL);
268  	   assert(child != NULL);
269  	
270  	   exprhdlr = SCIPfindExprhdlr(scip, EXPRHDLR_NAME);
271  	   if( exprhdlr == NULL )
272  	   {
273  	      SCIPerrorMessage("could not find %s expression handler -> abort\n", EXPRHDLR_NAME);
274  	      SCIPABORT();
275  	      return SCIP_PLUGINNOTFOUND;
276  	   }
277  	
278  	   /* create expression */
279  	   SCIP_CALL( SCIPcreateExpr(scip, expr, exprhdlr, NULL, 1, &child, ownercreate, ownercreatedata) );
280  	
281  	   return SCIP_OKAY;
282  	}
283  	
284  	/** indicates whether expression is of erf-type */  /*lint -e{715}*/
285  	SCIP_Bool SCIPisExprErf(
286  	   SCIP*                 scip,               /**< SCIP data structure */
287  	   SCIP_EXPR*            expr                /**< expression */
288  	   )
289  	{  /*lint --e{715}*/
290  	   assert(expr != NULL);
291  	
292  	   return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0;
293  	}
294  	
295  	/** creates the handler for erf expressions and includes it into SCIP
296  	 *
297  	 * @attention The implementation of this expression handler is incomplete.
298  	 * It is not usable for most use cases so far.
299  	 */
300  	SCIP_RETCODE SCIPincludeExprhdlrErf(
301  	   SCIP*                 scip                /**< SCIP data structure */
302  	   )
303  	{
304  	   SCIP_EXPRHDLR* exprhdlr;
305  	
306  	   /* include expression handler */
307  	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalErf, NULL) );
308  	   assert(exprhdlr != NULL);
309  	
310  	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrErf, NULL);
311  	   SCIPexprhdlrSetSimplify(exprhdlr, simplifyErf);
312  	   SCIPexprhdlrSetParse(exprhdlr, parseErf);
313  	   SCIPexprhdlrSetIntEval(exprhdlr, intevalErf);
314  	   SCIPexprhdlrSetHash(exprhdlr, hashErf);
315  	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffErf, NULL, NULL);
316  	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureErf);
317  	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityErf);
318  	   SCIPexprhdlrSetIntegrality(exprhdlr, integralityErf);
319  	
320  	   return SCIP_OKAY;
321  	}
322