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_abs.c
26   	 * @ingroup DEFPLUGINS_EXPR
27   	 * @brief  absolute expression handler
28   	 * @author Stefan Vigerske
29   	 * @author Benjamin Mueller
30   	 */
31   	
32   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
33   	
34   	#include <string.h>
35   	
36   	#include "scip/expr_value.h"
37   	#include "scip/expr_abs.h"
38   	#include "scip/expr.h"
39   	
40   	#define EXPRHDLR_NAME         "abs"
41   	#define EXPRHDLR_DESC         "absolute value expression"
42   	#define EXPRHDLR_PRECEDENCE   70000
43   	#define EXPRHDLR_HASHKEY      SCIPcalcFibHash(7187.0)
44   	
45   	/*
46   	 * Data structures
47   	 */
48   	
49   	/*
50   	 * Local methods
51   	 */
52   	
53   	/** computes both tangent underestimates and secant */
54   	static
55   	SCIP_RETCODE computeCutsAbs(
56   	   SCIP*                 scip,               /**< SCIP data structure */
57   	   SCIP_INTERVAL         bounds,             /**< bounds of child */
58   	   SCIP_Bool             overestimate,       /**< whether the expression shall be overestimated or underestimated */
59   	   SCIP_Real**           coefs,              /**< buffer to store coefficients of computed estimators */
60   	   SCIP_Real*            constant,           /**< buffer to store constant of computed estimators */
61   	   int*                  nreturned           /**< buffer to store number of estimators that have been computed */
62   	   )
63   	{
64   	   assert(scip != NULL);
65   	
66   	   *nreturned = 0;
67   	
68   	   /**! [SnippetExprInitestimatesAbs] */
69   	   if( !overestimate )
70   	   {
71   	      /* compute left tangent -x <= z */
72   	      coefs[*nreturned][0] = -1.0;
73   	      constant[*nreturned] = 0.0;
74   	      (*nreturned)++;
75   	
76   	      /* compute right tangent x <= z */
77   	      coefs[*nreturned][0] = 1.0;
78   	      constant[*nreturned] = 0.0;
79   	      (*nreturned)++;
80   	   }
81   	
82   	   /* compute secant */
83   	   if( overestimate )
84   	   {
85   	      SCIP_Real lb;
86   	      SCIP_Real ub;
87   	
88   	      lb = bounds.inf;
89   	      ub = bounds.sup;
90   	
91   	      /* it does not make sense to add a cut if child variable is unbounded or fixed */
92   	      if( !SCIPisInfinity(scip, -lb) && !SCIPisInfinity(scip, ub) && !SCIPisEQ(scip, lb, ub) )
93   	      {
94   	         if( !SCIPisPositive(scip, ub) )
95   	         {
96   	            /* z = -x, so add z <= -x here (-x <= z is the underestimator that is added above) */
97   	            coefs[*nreturned][0] = -1.0;
98   	            constant[*nreturned] = 0.0;
99   	            (*nreturned)++;
100  	         }
101  	         else if( !SCIPisNegative(scip, lb) )
102  	         {
103  	            /* z =  x, so add z <= x here (x <= z is the underestimator that is added above) */
104  	            coefs[*nreturned][0] = 1.0;
105  	            constant[*nreturned] = 0.0;
106  	            (*nreturned)++;
107  	         }
108  	         else
109  	         {
110  	            /* z = abs(x), x still has mixed sign */
111  	            SCIP_Real alpha;
112  	
113  	            /* let alpha = (|ub|-|lb|) / (ub-lb) then the resulting secant looks like
114  	             *
115  	             * z - |ub| <= alpha * (x - ub)  <=> z <= alpha * x + |ub| - alpha * ub
116  	             */
117  	            alpha = (REALABS(ub) - REALABS(lb)) / (ub - lb);
118  	
119  	            coefs[*nreturned][0] = alpha;
120  	            constant[*nreturned] = REALABS(ub) - alpha * ub;
121  	            (*nreturned)++;
122  	         }
123  	      }
124  	   }
125  	   /**! [SnippetExprInitestimatesAbs] */
126  	
127  	   return SCIP_OKAY;
128  	}
129  	
130  	
131  	/*
132  	 * Callback methods of expression handler
133  	 */
134  	
135  	/** simplifies an abs expression
136  	 *
137  	 * Evaluates the absolute value function when its child is a value expression.
138  	 *
139  	 * TODO: abs(*) = * if * >= 0 or - * if * < 0
140  	 */
141  	static
142  	SCIP_DECL_EXPRSIMPLIFY(simplifyAbs)
143  	{  /*lint --e{715}*/
144  	   SCIP_EXPR* child;
145  	
146  	   assert(scip != NULL);
147  	   assert(expr != NULL);
148  	   assert(simplifiedexpr != NULL);
149  	   assert(SCIPexprGetNChildren(expr) == 1);
150  	
151  	   child = SCIPexprGetChildren(expr)[0];
152  	   assert(child != NULL);
153  	
154  	   /* check for value expression */
155  	   if( SCIPisExprValue(scip, child) )
156  	   {
157  	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, REALABS(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) );
158  	   }
159  	   else
160  	   {
161  	      *simplifiedexpr = expr;
162  	
163  	      /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
164  	      SCIPcaptureExpr(*simplifiedexpr);
165  	   }
166  	
167  	   return SCIP_OKAY;
168  	}
169  	
170  	static
171  	SCIP_DECL_EXPRCOPYHDLR(copyhdlrAbs)
172  	{  /*lint --e{715}*/
173  	   SCIP_CALL( SCIPincludeExprhdlrAbs(scip) );
174  	
175  	   return SCIP_OKAY;
176  	}
177  	
178  	static
179  	SCIP_DECL_EXPRPARSE(parseAbs)
180  	{  /*lint --e{715}*/
181  	   SCIP_EXPR* childexpr;
182  	
183  	   assert(expr != NULL);
184  	
185  	   /* parse child expression from remaining string */
186  	   SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
187  	   assert(childexpr != NULL);
188  	
189  	   /* create absolute expression */
190  	   SCIP_CALL( SCIPcreateExprAbs(scip, expr, childexpr, ownercreate, ownercreatedata) );
191  	   assert(*expr != NULL);
192  	
193  	   /* release child expression since it has been captured by the absolute expression */
194  	   SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
195  	
196  	   *success = TRUE;
197  	
198  	   return SCIP_OKAY;
199  	}
200  	
201  	/** expression point evaluation callback */
202  	static
203  	SCIP_DECL_EXPREVAL(evalAbs)
204  	{  /*lint --e{715}*/
205  	   assert(expr != NULL);
206  	   assert(SCIPexprGetNChildren(expr) == 1);
207  	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
208  	
209  	   *val = REALABS(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
210  	
211  	   return SCIP_OKAY;
212  	}
213  	
214  	
215  	/** expression derivative evaluation callback */
216  	static
217  	SCIP_DECL_EXPRBWDIFF(bwdiffAbs)
218  	{  /*lint --e{715}*/
219  	   SCIP_EXPR* child;
220  	
221  	   assert(expr != NULL);
222  	   assert(childidx == 0);
223  	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
224  	
225  	   child = SCIPexprGetChildren(expr)[0];
226  	   assert(child != NULL);
227  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
228  	
229  	   *val = (SCIPexprGetEvalValue(child) >= 0.0) ? 1.0 : -1.0;
230  	
231  	   return SCIP_OKAY;
232  	}
233  	
234  	/** expression interval evaluation callback */
235  	static
236  	SCIP_DECL_EXPRINTEVAL(intevalAbs)
237  	{  /*lint --e{715}*/
238  	   SCIP_INTERVAL childinterval;
239  	
240  	   assert(expr != NULL);
241  	   assert(SCIPexprGetNChildren(expr) == 1);
242  	
243  	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
244  	
245  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
246  	      SCIPintervalSetEmpty(interval);
247  	   else
248  	      SCIPintervalAbs(SCIP_INTERVAL_INFINITY, interval, childinterval);
249  	
250  	   return SCIP_OKAY;
251  	}
252  	
253  	/** expression estimator callback */
254  	static
255  	SCIP_DECL_EXPRESTIMATE(estimateAbs)
256  	{  /*lint --e{715}*/
257  	   assert(scip != NULL);
258  	   assert(expr != NULL);
259  	   assert(SCIPexprGetNChildren(expr) == 1);
260  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
261  	   assert(coefs != NULL);
262  	   assert(constant != NULL);
263  	   assert(islocal != NULL);
264  	   assert(branchcand != NULL);
265  	   assert(*branchcand == TRUE);
266  	   assert(success != NULL);
267  	
268  	   SCIPdebugMsg(scip, "%sestimate |child| over locdom=[%g,%g] glbdom=[%g,%g]\n", overestimate ? "over" : "under",
269  	      localbounds[0].inf, localbounds[0].sup, globalbounds[0].inf, globalbounds[0].sup);
270  	
271  	   /**! [SnippetExprEstimateAbs] */
272  	   if( !overestimate )
273  	   {
274  	      *constant = 0.0;
275  	
276  	      if( refpoint[0] <= 0.0 )
277  	         *coefs = -1.0;
278  	      else
279  	         *coefs = 1.0;
280  	
281  	      *islocal = FALSE;
282  	      *branchcand = FALSE;
283  	   }
284  	   else
285  	   {
286  	      /* overestimator */
287  	      SCIP_Real lb;
288  	      SCIP_Real ub;
289  	
290  	      lb = localbounds[0].inf;
291  	      ub = localbounds[0].sup;
292  	
293  	      if( !SCIPisPositive(scip, ub) )
294  	      {
295  	         /* |x| = -x */
296  	         *coefs = -1.0;
297  	         *constant = 0.0;
298  	         *islocal = SCIPisPositive(scip, globalbounds[0].sup);
299  	         *branchcand = FALSE;
300  	      }
301  	      else if( !SCIPisNegative(scip, lb) )
302  	      {
303  	         /* |x| = x */
304  	         *coefs =  1.0;
305  	         *constant = 0.0;
306  	         *islocal = SCIPisNegative(scip, globalbounds[0].inf);
307  	         *branchcand = FALSE;
308  	      }
309  	      else if( !SCIPisRelEQ(scip, lb, -ub) )
310  	      {
311  	         /* |x| with x having mixed sign and ub+lb does not cancel out -> secant */
312  	         SCIP_Real alpha;
313  	
314  	         assert(lb < 0.0);
315  	         assert(ub > 0.0);
316  	
317  	         /* let alpha = (|ub|-|lb|) / (ub-lb) = (ub+lb)/(ub-lb)
318  	          * then the resulting secant is -lb + alpha * (x - lb) = -lb - alpha*lb + alpha*x
319  	          */
320  	         alpha = (ub + lb) / (ub - lb);
321  	
322  	         *coefs = alpha;
323  	         *constant = -lb - alpha * lb;
324  	         *islocal = TRUE;
325  	      }
326  	      else if( lb == -ub ) /*lint !e777*/
327  	      {
328  	         /* alpha = 0 */
329  	         *coefs = 0.0;
330  	         *constant = -lb;
331  	         *islocal = TRUE;
332  	      }
333  	      else
334  	      {
335  	         *success = FALSE;
336  	         return SCIP_OKAY;
337  	      }
338  	   }
339  	   /**! [SnippetExprEstimateAbs] */
340  	
341  	   SCIPdebugMsg(scip, "-> %g * <child> %+g, local=%u branchcand=%u\n", *coefs, *constant, *islocal, *branchcand);
342  	
343  	   *success = TRUE;
344  	
345  	   return SCIP_OKAY;
346  	}
347  	
348  	/** expression estimate initialization callback */
349  	static
350  	SCIP_DECL_EXPRINITESTIMATES(initEstimatesAbs)
351  	{  /*lint --e{715}*/
352  	   assert(expr != NULL);
353  	   assert(SCIPexprGetNChildren(expr) == 1);
354  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0);
355  	
356  	   /* compute initial cuts */
357  	   SCIP_CALL( computeCutsAbs(scip, bounds[0], overestimate, coefs, constant, nreturned) );
358  	
359  	   return SCIP_OKAY;
360  	}
361  	
362  	/** expression reverse propagation callback */
363  	static
364  	SCIP_DECL_EXPRREVERSEPROP(reversepropAbs)
365  	{  /*lint --e{715}*/
366  	   SCIP_INTERVAL childbounds;
367  	   SCIP_INTERVAL left;
368  	   SCIP_INTERVAL right;
369  	
370  	   assert(scip != NULL);
371  	   assert(expr != NULL);
372  	   assert(SCIPexprGetNChildren(expr) == 1);
373  	   assert(bounds.inf >= 0.0);  /* bounds should have been intersected with activity, which is >= 0 */
374  	
375  	   /**! [SnippetExprReversepropAbs] */
376  	   /* abs(x) in I -> x \in (-I \cup I) \cap bounds(x) */
377  	   right = bounds;  /* I */
378  	   SCIPintervalSetBounds(&left, -right.sup, -right.inf); /* -I */
379  	
380  	   childbounds = childrenbounds[0];
381  	   /* childbounds can be empty here already, but that should work fine here */
382  	
383  	   SCIPintervalIntersect(&left, left, childbounds);    /* -I \cap bounds(x), could become empty */
384  	   SCIPintervalIntersect(&right, right, childbounds);  /*  I \cap bounds(x), could become empty */
385  	
386  	   /* compute smallest interval containing (-I \cap bounds(x)) \cup (I \cap bounds(x)) = (-I \cup I) \cap bounds(x)
387  	    * this works also if left or right is empty
388  	    */
389  	   SCIPintervalUnify(&childrenbounds[0], left, right);
390  	   /**! [SnippetExprReversepropAbs] */
391  	
392  	   return SCIP_OKAY;
393  	}
394  	
395  	/** expression hash callback */
396  	static
397  	SCIP_DECL_EXPRHASH(hashAbs)
398  	{  /*lint --e{715}*/
399  	   assert(scip != NULL);
400  	   assert(expr != NULL);
401  	   assert(SCIPexprGetNChildren(expr) == 1);
402  	   assert(hashkey != NULL);
403  	   assert(childrenhashes != NULL);
404  	
405  	   *hashkey = EXPRHDLR_HASHKEY;
406  	   *hashkey ^= childrenhashes[0];
407  	
408  	   return SCIP_OKAY;
409  	}
410  	
411  	/** expression curvature detection callback */
412  	static
413  	SCIP_DECL_EXPRCURVATURE(curvatureAbs)
414  	{  /*lint --e{715}*/
415  	   SCIP_EXPR* child;
416  	   SCIP_INTERVAL childbounds;
417  	   SCIP_Real childinf;
418  	   SCIP_Real childsup;
419  	
420  	   assert(scip != NULL);
421  	   assert(expr != NULL);
422  	   assert(exprcurvature != SCIP_EXPRCURV_UNKNOWN);
423  	   assert(success != NULL);
424  	   assert(childcurv != NULL);
425  	   assert(SCIPexprGetNChildren(expr) == 1);
426  	
427  	   child = SCIPexprGetChildren(expr)[0];
428  	   assert(child != NULL);
429  	
430  	   /**! [SnippetExprCurvatureAbs] */
431  	   /* expression is |child|, get domain of child */
432  	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
433  	   childbounds = SCIPexprGetActivity(child);
434  	   childinf = SCIPintervalGetInf(childbounds);
435  	   childsup = SCIPintervalGetSup(childbounds);
436  	
437  	   *success = TRUE;
438  	   if( childinf >= 0.0 )  /* |f(x)| = f(x) */
439  	      childcurv[0] = exprcurvature;
440  	   else if( childsup <= 0.0 ) /* |f(x)| = -f(x) */
441  	      childcurv[0] = SCIPexprcurvNegate(exprcurvature);
442  	   else if( exprcurvature == SCIP_EXPRCURV_CONVEX )   /* |f(x)|, f mixed sign, is convex if f is linear */
443  	      childcurv[0] = SCIP_EXPRCURV_LINEAR;
444  	   else /* |f(x)|, f mixed sign, is never concave nor linear */
445  	      *success = FALSE;
446  	   /**! [SnippetExprCurvatureAbs] */
447  	
448  	   return SCIP_OKAY;
449  	}
450  	
451  	/** expression monotonicity detection callback */
452  	static
453  	SCIP_DECL_EXPRMONOTONICITY(monotonicityAbs)
454  	{  /*lint --e{715}*/
455  	   SCIP_EXPR* child;
456  	   SCIP_INTERVAL childbounds;
457  	
458  	   assert(scip != NULL);
459  	   assert(expr != NULL);
460  	   assert(result != NULL);
461  	   assert(childidx == 0);
462  	
463  	   child = SCIPexprGetChildren(expr)[0];
464  	   assert(child != NULL);
465  	
466  	   /**! [SnippetExprMonotonicityAbs] */
467  	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
468  	   childbounds = SCIPexprGetActivity(child);
469  	
470  	   if( childbounds.sup <= 0.0 )
471  	      *result = SCIP_MONOTONE_DEC;
472  	   else if( childbounds.inf >= 0.0 )
473  	      *result = SCIP_MONOTONE_INC;
474  	   else
475  	      *result = SCIP_MONOTONE_UNKNOWN;
476  	   /**! [SnippetExprMonotonicityAbs] */
477  	
478  	   return SCIP_OKAY;
479  	}
480  	
481  	/** expression integrality detection callback */
482  	static
483  	SCIP_DECL_EXPRINTEGRALITY(integralityAbs)
484  	{  /*lint --e{715}*/
485  	   SCIP_EXPR* child;
486  	
487  	   assert(scip != NULL);
488  	   assert(expr != NULL);
489  	   assert(isintegral != NULL);
490  	   assert(SCIPexprGetNChildren(expr) == 1);
491  	
492  	   child = SCIPexprGetChildren(expr)[0];
493  	   assert(child != NULL);
494  	
495  	   *isintegral = SCIPexprIsIntegral(child);
496  	
497  	   return SCIP_OKAY;
498  	}
499  	
500  	
501  	/** creates the handler for absolute expression and includes it into SCIP */
502  	SCIP_RETCODE SCIPincludeExprhdlrAbs(
503  	   SCIP*                 scip                /**< SCIP data structure */
504  	   )
505  	{
506  	   SCIP_EXPRHDLR* exprhdlr;
507  	
508  	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC,
509  	         EXPRHDLR_PRECEDENCE, evalAbs, NULL) );
510  	   assert(exprhdlr != NULL);
511  	
512  	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrAbs, NULL);
513  	   SCIPexprhdlrSetSimplify(exprhdlr, simplifyAbs);
514  	   SCIPexprhdlrSetParse(exprhdlr, parseAbs);
515  	   SCIPexprhdlrSetIntEval(exprhdlr, intevalAbs);
516  	   SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesAbs, estimateAbs);
517  	   SCIPexprhdlrSetHash(exprhdlr, hashAbs);
518  	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropAbs);
519  	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffAbs, NULL, NULL);
520  	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureAbs);
521  	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityAbs);
522  	   SCIPexprhdlrSetIntegrality(exprhdlr, integralityAbs);
523  	
524  	   return SCIP_OKAY;
525  	}
526  	
527  	/** creates an absolute value expression */
528  	SCIP_RETCODE SCIPcreateExprAbs(
529  	   SCIP*                 scip,               /**< SCIP data structure */
530  	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
531  	   SCIP_EXPR*            child,              /**< single child */
532  	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
533  	   void*                 ownercreatedata     /**< data to pass to ownercreate */
534  	   )
535  	{
536  	   assert(expr != NULL);
537  	   assert(child != NULL);
538  	   assert(SCIPfindExprhdlr(scip, EXPRHDLR_NAME) != NULL);
539  	
540  	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, EXPRHDLR_NAME), NULL, 1, &child, ownercreate, ownercreatedata) );
541  	
542  	   return SCIP_OKAY;
543  	}
544  	
545  	/** indicates whether expression is of abs-type */  /*lint -e{715}*/
546  	SCIP_Bool SCIPisExprAbs(
547  	   SCIP*                 scip,               /**< SCIP data structure */
548  	   SCIP_EXPR*            expr                /**< expression */
549  	   )
550  	{  /*lint --e{715}*/
551  	   assert(expr != NULL);
552  	
553  	   return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0;
554  	}
555