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_trig.c
26   	 * @ingroup DEFPLUGINS_EXPR
27   	 * @brief  handler for sine and cosine expressions
28   	 * @author Fabian Wegscheider
29   	 *
30   	 * The estimator/separator code always computes underestimators for sin(x).
31   	 * For overestimators of cos(x), we first reduce to underestimators of sin(x).
32   	 *
33   	 * Overestimator for sin(x):
34   	 *   Assume that a*y+b <= sin(y) for y in [-ub,-lb].
35   	 *   Then we have a*(-y)-b >= -sin(y) = sin(-y) for y in [-ub,-lb].
36   	 *   Thus, a*x-b >= sin(x) for x in [lb,ub].
37   	 *
38   	 * Underestimator for cos(x):
39   	 *   Assume that a*y+b <= sin(y) for y in [lb+pi/2,ub+pi/2].
40   	 *   Then we have a*(x+pi/2) + b <= sin(x+pi/2) = cos(x) for x in [lb,ub].
41   	 *   Thus, a*x + (b+a*pi/2) <= cos(x) for x in [lb,ub].
42   	 *
43   	 * Overestimator for cos(x):
44   	 *   Assume that a*z+b <= sin(z) for z in [-(ub+pi/2),-(lb+pi/2)].
45   	 *   Then, a*y-b >= sin(y) for y in [lb+pi/2,ub+pi/2].
46   	 *   Then, a*x-b+a*pi/2 >= cos(x) for x in [lb,ub].
47   	 */
48   	
49   	/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
50   	
51   	#define _USE_MATH_DEFINES   /* to get M_PI on Windows */  /*lint !750 */
52   	
53   	#include <string.h>
54   	#include <math.h>
55   	#include "scip/expr_trig.h"
56   	#include "scip/expr_value.h"
57   	
58   	/* fundamental expression handler properties */
59   	#define SINEXPRHDLR_NAME         "sin"
60   	#define SINEXPRHDLR_DESC         "sine expression"
61   	#define SINEXPRHDLR_PRECEDENCE   91000
62   	#define SINEXPRHDLR_HASHKEY      SCIPcalcFibHash(82457.0)
63   	
64   	#define COSEXPRHDLR_NAME         "cos"
65   	#define COSEXPRHDLR_DESC         "cosine expression"
66   	#define COSEXPRHDLR_PRECEDENCE   92000
67   	#define COSEXPRHDLR_HASHKEY      SCIPcalcFibHash(82463.0)
68   	
69   	#define MAXCHILDABSVAL        1e+6                       /**< maximum absolute value that is accepted for propagation */
70   	#define NEWTON_NITERATIONS    100
71   	#define NEWTON_PRECISION      1e-12
72   	
73   	/*
74   	 * Local methods
75   	 */
76   	
77   	/** evaluates the function a*x + b - sin(x) for some coefficient a and constant b at a given point p
78   	 *
79   	 *  the constants a and b are expected to be stored in that order in params
80   	 */
81   	static
82   	SCIP_DECL_NEWTONEVAL(function1)
83   	{  /*lint --e{715}*/
84   	   assert(params != NULL);
85   	   assert(nparams == 2);
86   	
87   	   return params[0]*point + params[1] - sin(point);
88   	}
89   	
90   	/** evaluates the derivative of a*x + b - sin(x) for some coefficient a and constant b at a given point p
91   	 *
92   	 *  the constants a and b are expected to be stored in that order in params
93   	 */
94   	static
95   	SCIP_DECL_NEWTONEVAL(derivative1)
96   	{  /*lint --e{715}*/
97   	   assert(params != NULL);
98   	   assert(nparams == 2);
99   	
100  	   return params[0] - cos(point);
101  	}
102  	
103  	/** evaluates the function sin(x) + (alpha - x)*cos(x) - sin(alpha) for some constant alpha at a given point p
104  	 *
105  	 *  the constant alpha is expected to be stored in params
106  	 */
107  	static
108  	SCIP_DECL_NEWTONEVAL(function2)
109  	{  /*lint --e{715}*/
110  	   assert(params != NULL);
111  	   assert(nparams == 1);
112  	
113  	   return sin(point) + (params[0] - point) * cos(point) - sin(params[0]);
114  	}
115  	
116  	/** evaluates the derivative of sin(x) + (alpha - x)*cos(x) - sin(alpha) for some constant alpha at a given point p
117  	 *
118  	 *  the constant alpha is expected to be stored in params
119  	 */
120  	static
121  	SCIP_DECL_NEWTONEVAL(derivative2)
122  	{  /*lint --e{715}*/
123  	   assert(params != NULL);
124  	   assert(nparams == 1);
125  	
126  	   return (point - params[0]) * sin(point);
127  	}
128  	
129  	/** helper function to compute the secant if it is a valid underestimator
130  	 *
131  	 *  returns true if the estimator was computed successfully
132  	 */
133  	static
134  	SCIP_Bool computeSecantSin(
135  	   SCIP*                 scip,               /**< SCIP data structure */
136  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of secant */
137  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of secant */
138  	   SCIP_Real             lb,                 /**< lower bound of argument variable */
139  	   SCIP_Real             ub                  /**< upper bound of argument variable */
140  	   )
141  	{
142  	   assert(scip != NULL);
143  	   assert(lincoef != NULL);
144  	   assert(linconst != NULL);
145  	   assert(lb < ub);
146  	
147  	   /* if range is too big, secant is not underestimating */
148  	   if( ub - lb >= M_PI )
149  	      return FALSE;
150  	
151  	   /* if bounds are not within positive bay, secant is not underestimating */
152  	   if( sin(lb) < 0.0 || sin(ub) < 0.0  || (sin(lb) == 0.0 && cos(lb) < 0.0) )
153  	      return FALSE;
154  	
155  	   *lincoef = (sin(ub) - sin(lb)) / (ub - lb);
156  	   *linconst = sin(ub) - (*lincoef) * ub;
157  	
158  	   return TRUE;
159  	}
160  	
161  	/** helper function to compute the tangent at lower bound if it is underestimating
162  	 *
163  	 *  returns true if the underestimator was computed successfully
164  	 */
165  	static
166  	SCIP_Bool computeLeftTangentSin(
167  	   SCIP*                 scip,               /**< SCIP data structure */
168  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of tangent */
169  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of tangent */
170  	   SCIP_Real             lb                  /**< lower bound of argument variable */
171  	   )
172  	{
173  	   assert(scip != NULL);
174  	   assert(lincoef != NULL);
175  	   assert(linconst != NULL);
176  	
177  	   if( SCIPisInfinity(scip, -lb) )
178  	      return FALSE;
179  	
180  	   /* left tangent is only underestimating in [pi, 1.5*pi) *2kpi */
181  	   if( sin(lb) > 0.0 || cos(lb) >= 0.0 )
182  	      return FALSE;
183  	
184  	   *lincoef = cos(lb);
185  	   *linconst = sin(lb) - (*lincoef) * lb;
186  	
187  	   return TRUE;
188  	}
189  	
190  	/* TODO: fix this, more cases can be considered, see at unit test
191  	 * the underestimating of the tangents depends not only on the ub but also on the lower bound.
192  	 * right now, this function is only checking whether the tangent underestimates independently of the lower bound!
193  	 */
194  	/** helper function to compute the tangent at upper bound if it is an underestimator
195  	 *
196  	 *  returns true if the underestimator was computed successfully
197  	 */
198  	static
199  	SCIP_Bool computeRightTangentSin(
200  	   SCIP*                 scip,               /**< SCIP data structure */
201  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of tangent */
202  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of tangent */
203  	   SCIP_Real             ub                  /**< upper bound of argument variable */
204  	   )
205  	{
206  	   assert(scip != NULL);
207  	   assert(lincoef != NULL);
208  	   assert(linconst != NULL);
209  	
210  	   if( SCIPisInfinity(scip, ub) )
211  	      return FALSE;
212  	
213  	   /* right tangent is only underestimating in (1.5*pi, 2*pi] *2kpi */
214  	   if( sin(ub) > 0.0 || cos(ub) <= 0.0 )
215  	      return FALSE;
216  	
217  	   *lincoef = cos(ub);
218  	   *linconst = sin(ub) - (*lincoef) * ub;
219  	
220  	   return TRUE;
221  	}
222  	
223  	/** helper function to compute the tangent at solution point if it is an underestimator
224  	 *
225  	 *  returns true if the underestimator was computed successfully
226  	 */
227  	static
228  	SCIP_Bool computeSolTangentSin(
229  	   SCIP*                 scip,               /**< SCIP data structure */
230  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of tangent */
231  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of tangent */
232  	   SCIP_Real             lb,                 /**< lower bound of argument variable */
233  	   SCIP_Real             ub,                 /**< upper bound of argument variable */
234  	   SCIP_Real             solpoint            /**< solution point to be separated */
235  	   )
236  	{
237  	   SCIP_Real params[2];
238  	   SCIP_Real startingpoints[3];
239  	   SCIP_Real solpointmodpi;
240  	   SCIP_Real intersection;
241  	   int i;
242  	
243  	   assert(scip != NULL);
244  	   assert(lincoef != NULL);
245  	   assert(linconst != NULL);
246  	
247  	   /* tangent is only underestimating in negative bay */
248  	   if( sin(solpoint) > 0.0 )
249  	      return FALSE;
250  	
251  	   /* compute solution point mod pi */
252  	   solpointmodpi = fmod(solpoint, M_PI);
253  	   if( solpoint < 0.0 )
254  	      solpointmodpi += M_PI;
255  	
256  	   /* if the point is too far away from the bounds or is at a multiple of pi, then tangent is not underestimating */
257  	   if( SCIPisGE(scip, solpoint - lb, 2*M_PI) || SCIPisGE(scip, ub - solpoint, 2*M_PI)
258  	      || SCIPisZero(scip, solpointmodpi) )
259  	      return FALSE;
260  	
261  	   params[0] = cos(solpoint);
262  	   params[1] = sin(solpoint) - params[0] * solpoint;
263  	
264  	   /* choose starting points for Newton procedure */
265  	   if( SCIPisGT(scip, solpointmodpi, M_PI_2) )
266  	   {
267  	      startingpoints[0] = solpoint + (M_PI - solpointmodpi) + M_PI_2;
268  	      startingpoints[1] = startingpoints[0] + M_PI_2;
269  	      startingpoints[2] = startingpoints[1] + M_PI_2;
270  	   }
271  	   else
272  	   {
273  	      startingpoints[0] = solpoint - solpointmodpi - M_PI_2;
274  	      startingpoints[1] = startingpoints[0] - M_PI_2;
275  	      startingpoints[2] = startingpoints[1] - M_PI_2;
276  	   }
277  	
278  	   /* use Newton procedure to test if cut is valid */
279  	   for( i = 0; i < 3; ++i )
280  	   {
281  	      intersection = SCIPcalcRootNewton(function1, derivative1, params, 2, startingpoints[i], NEWTON_PRECISION,
282  	         NEWTON_NITERATIONS);
283  	
284  	      if( intersection != SCIP_INVALID && !SCIPisEQ(scip, intersection, solpoint) ) /*lint !e777*/
285  	         break;
286  	   }
287  	
288  	   /* if Newton failed or intersection point lies within bounds, underestimator is not valid */
289  	   if( intersection == SCIP_INVALID || (intersection >= lb && intersection <= ub) ) /*lint !e777*/
290  	      return FALSE;
291  	
292  	   *lincoef = params[0];
293  	   *linconst = params[1];
294  	
295  	   return TRUE;
296  	}
297  	
298  	/** helper function to compute the secant between lower bound and some point of the graph such that it underestimates
299  	 *
300  	 *  returns true if the underestimator was computed successfully
301  	 */
302  	static
303  	SCIP_Bool computeLeftSecantSin(
304  	   SCIP*                 scip,               /**< SCIP data structure */
305  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of tangent */
306  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of tangent */
307  	   SCIP_Real             lb,                 /**< lower bound of argument variable */
308  	   SCIP_Real             ub                  /**< upper bound of argument variable */
309  	   )
310  	{
311  	   SCIP_Real lbmodpi;
312  	   SCIP_Real tangentpoint;
313  	   SCIP_Real startingpoint;
314  	
315  	   assert(scip != NULL);
316  	   assert(lincoef != NULL);
317  	   assert(linconst != NULL);
318  	   assert(lb < ub);
319  	
320  	   if( SCIPisInfinity(scip, -lb) )
321  	      return FALSE;
322  	
323  	   /* compute shifted bounds for case evaluation */
324  	   lbmodpi = fmod(lb, M_PI);
325  	   if( lb < 0.0 )
326  	      lbmodpi += M_PI;
327  	
328  	   /* choose starting point for Newton procedure */
329  	   if( cos(lb) < 0.0 )
330  	   {
331  	      /* in [pi/2,pi] underestimating doesn't work; otherwise, take the midpoint of possible area */
332  	      if( SCIPisLE(scip, sin(lb), 0.0) )
333  	         return FALSE;
334  	      else
335  	         startingpoint = lb + 1.25*M_PI - lbmodpi;
336  	   }
337  	   else
338  	   {
339  	      /* in ascending area, take the midpoint of the possible area in descending part */
340  	      /* for lb < 0 but close to zero, we may have sin(lb) = 0 but lbmodpi = pi, which gives a starting point too close to lb
341  	       * but for sin(lb) around 0 we know that the tangent point needs to be in [lb+pi,lb+pi+pi/2]
342  	       */
343  	      if( SCIPisZero(scip, sin(lb)) )
344  	         startingpoint = lb + 1.25*M_PI;
345  	      else if( sin(lb) < 0.0 )
346  	         startingpoint = lb + 2.25*M_PI - lbmodpi;
347  	      else
348  	         startingpoint = lb + 1.25*M_PI - lbmodpi;
349  	   }
350  	
351  	   /* use Newton procedure to find the point where the tangent intersects sine at lower bound */
352  	   tangentpoint = SCIPcalcRootNewton(function2, derivative2, &lb, 1, startingpoint, NEWTON_PRECISION,
353  	      NEWTON_NITERATIONS);
354  	
355  	   /* if Newton procedure failed, no cut is added */
356  	   if( tangentpoint == SCIP_INVALID ) /*lint !e777*/
357  	      return FALSE;
358  	
359  	   /* if the computed point lies outside the bounds, it is shifted to upper bound */
360  	   if( SCIPisGE(scip, tangentpoint, ub) )
361  	   {
362  	      tangentpoint = ub;
363  	
364  	      /* check whether affine function is still underestimating */
365  	      if( SCIPisLE(scip, sin(0.5 * (ub + lb)), sin(lb) + 0.5*(sin(ub) - sin(lb))) )
366  	         return FALSE;
367  	   }
368  	
369  	   if( SCIPisEQ(scip, tangentpoint, lb) )  /*lint !e777 */
370  	      return FALSE;
371  	
372  	   /* compute secant between lower bound and connection point */
373  	   *lincoef = (sin(tangentpoint) - sin(lb)) / (tangentpoint - lb);
374  	   *linconst = sin(lb) - (*lincoef) * lb;
375  	
376  	   /* if the bounds are too close to each other, it's possible that the underestimator is not valid */
377  	   if( *lincoef >= cos(lb) )
378  	      return FALSE;
379  	
380  	   SCIPdebugMsg(scip, "left secant: %g + %g*x <= sin(x) on [%g,%g]\n", *linconst, *lincoef, lb, ub);
381  	
382  	   return TRUE;
383  	}
384  	
385  	/** helper function to compute the secant between upper bound and some point of the graph such that it underestimates
386  	 *
387  	 *  returns true if the underestimator was computed successfully
388  	 */
389  	static
390  	SCIP_Bool computeRightSecantSin(
391  	   SCIP*                 scip,               /**< SCIP data structure */
392  	   SCIP_Real*            lincoef,            /**< buffer to store linear coefficient of tangent */
393  	   SCIP_Real*            linconst,           /**< buffer to store linear constant of tangent */
394  	   SCIP_Real             lb,                 /**< lower bound of argument variable */
395  	   SCIP_Real             ub                  /**< upper bound of argument variable */
396  	   )
397  	{
398  	   SCIP_Real ubmodpi;
399  	   SCIP_Real tangentpoint;
400  	   SCIP_Real startingpoint;
401  	
402  	   assert(scip != NULL);
403  	   assert(lincoef != NULL);
404  	   assert(linconst != NULL);
405  	   assert(lb < ub);
406  	
407  	   if( SCIPisInfinity(scip, ub) )
408  	      return FALSE;
409  	
410  	   /* compute shifted bounds for case evaluation */
411  	   ubmodpi = fmod(ub, M_PI);
412  	   if( ub < 0.0 )
413  	      ubmodpi += M_PI;
414  	
415  	   /* choose starting point for Newton procedure */
416  	   if( cos(ub) > 0.0 )
417  	   {
418  	      /* in [3*pi/2,2*pi] underestimating doesn't work; otherwise, take the midpoint of possible area */
419  	      if( SCIPisLE(scip, sin(ub), 0.0) )
420  	         return FALSE;
421  	      else
422  	         startingpoint = ub - M_PI_4 - ubmodpi;
423  	   }
424  	   else
425  	   {
426  	      /* in descending area, take the midpoint of the possible area in ascending part */
427  	      /* for ub < 0 but close to zero, we may have sin(ub) = 0 but ubmodpi = pi, which gives a starting point too close to ub
428  	       * but for sin(ub) around 0 we know that the tangent point needs to be in [ub-(pi+pi/2),ub-pi]
429  	       */
430  	      if( SCIPisZero(scip, sin(ub)) )
431  	         startingpoint = ub - 1.25*M_PI;
432  	      else if( sin(ub) < 0.0 )
433  	         startingpoint = ub - 1.25*M_PI - ubmodpi;
434  	      else
435  	         startingpoint = ub - M_PI_4 - ubmodpi;
436  	   }
437  	
438  	   /* use Newton procedure to find the point where the tangent intersects sine at lower bound */
439  	   tangentpoint = SCIPcalcRootNewton(function2, derivative2, &ub, 1, startingpoint, NEWTON_PRECISION,
440  	      NEWTON_NITERATIONS);
441  	
442  	   /* if Newton procedure failed, no underestimator is found */
443  	   if( tangentpoint == SCIP_INVALID ) /*lint !e777*/
444  	      return FALSE;
445  	
446  	   /* if the computed point lies outside the bounds, it is shifted to upper bound */
447  	   if( SCIPisLE(scip, tangentpoint, lb) )
448  	   {
449  	      tangentpoint = lb;
450  	
451  	      /* check whether affine function is still underestimating */
452  	      if( SCIPisLE(scip, sin(0.5 * (ub + lb)), sin(lb) + 0.5*(sin(ub) - sin(lb))) )
453  	         return FALSE;
454  	   }
455  	
456  	   if( SCIPisEQ(scip, tangentpoint, ub) )  /*lint !e777 */
457  	      return FALSE;
458  	
459  	   /* compute secant between lower bound and connection point */
460  	   *lincoef = (sin(tangentpoint) - sin(ub)) / (tangentpoint - ub);
461  	   *linconst = sin(ub) - (*lincoef) * ub;
462  	
463  	   /* if the bounds are to close to each other, it's possible that the underestimator is not valid */
464  	   if( *lincoef <= cos(lb) )
465  	      return FALSE;
466  	
467  	   return TRUE;
468  	}
469  	
470  	/** helper function to compute the new interval for child in reverse propagation */
471  	static
472  	SCIP_RETCODE computeRevPropIntervalSin(
473  	   SCIP*                 scip,               /**< SCIP data structure */
474  	   SCIP_INTERVAL         parentbounds,       /**< bounds for sine expression */
475  	   SCIP_INTERVAL         childbounds,        /**< bounds for child expression */
476  	   SCIP_INTERVAL*        newbounds           /**< buffer to store new child bounds */
477  	   )
478  	{
479  	   SCIP_Real newinf = childbounds.inf;
480  	   SCIP_Real newsup = childbounds.sup;
481  	
482  	   /* if the absolute values of the bounds are too large, skip reverse propagation
483  	    * TODO: if bounds are close but too large, shift them to [0,2pi] and do the computation there
484  	    */
485  	   if( ABS(newinf) > MAXCHILDABSVAL || ABS(newsup) > MAXCHILDABSVAL )
486  	   {
487  	      SCIPintervalSetBounds(newbounds, newinf, newsup);
488  	      return SCIP_OKAY;
489  	   }
490  	
491  	   if( !SCIPisInfinity(scip, -newinf) )
492  	   {
493  	      /* l(x) and u(x) are lower/upper bound of child, l(s) and u(s) are lower/upper bound of sin expr
494  	       *
495  	       * if sin(l(x)) < l(s), we are looking for k minimal s.t. a + 2k*pi > l(x) where a = asin(l(s))
496  	       * then the new lower bound is a + 2k*pi
497  	       */
498  	      if( SCIPisLT(scip, sin(newinf), parentbounds.inf) )
499  	      {
500  	         SCIP_Real a = asin(parentbounds.inf);
501  	         int k = (int) ceil((newinf - a) / (2.0*M_PI));
502  	         newinf = a + 2.0*M_PI * k;
503  	      }
504  	
505  	      /* if sin(l(x)) > u(s), we are looking for k minimal s.t. pi - a + 2k*pi > l(x) where a = asin(u(s))
506  	       * then the new lower bound is pi - a + 2k*pi
507  	       */
508  	      else if( SCIPisGT(scip, sin(newinf), parentbounds.sup) )
509  	      {
510  	         SCIP_Real a = asin(parentbounds.sup);
511  	         int k = (int) ceil((newinf + a) / (2.0*M_PI) - 0.5);
512  	         newinf = M_PI * (2.0*k + 1.0) - a;
513  	      }
514  	
515  	      assert(newinf >= childbounds.inf);
516  	      assert(SCIPisFeasGE(scip, sin(newinf), parentbounds.inf));
517  	      assert(SCIPisFeasLE(scip, sin(newinf), parentbounds.sup));
518  	   }
519  	
520  	   if( !SCIPisInfinity(scip, newsup) )
521  	   {
522  	      /* if sin(u(x)) > u(s), we are looking for k minimal s.t. a + 2k*pi > u(x) - 2*pi where a = asin(u(s))
523  	       * then the new upper bound is a + 2k*pi
524  	       */
525  	      if ( SCIPisGT(scip, sin(newsup), parentbounds.sup) )
526  	      {
527  	         SCIP_Real a = asin(parentbounds.sup);
528  	         int k = (int) ceil((newsup - a ) / (2.0*M_PI)) - 1;
529  	         newsup = a + 2.0*M_PI * k;
530  	      }
531  	
532  	      /* if sin(u(x)) < l(s), we are looking for k minimal s.t. pi - a + 2k*pi > l(x) - 2*pi where a = asin(l(s))
533  	       * then the new upper bound is pi - a + 2k*pi
534  	       */
535  	      if( SCIPisLT(scip, sin(newsup), parentbounds.inf) )
536  	      {
537  	         SCIP_Real a = asin(parentbounds.inf);
538  	         int k = (int) ceil((newsup + a) / (2.0*M_PI) - 0.5) - 1;
539  	         newsup = M_PI * (2.0*k + 1.0) - a;
540  	      }
541  	
542  	      assert(newsup <= childbounds.sup);
543  	      assert(SCIPisFeasGE(scip, sin(newsup), parentbounds.inf));
544  	      assert(SCIPisFeasLE(scip, sin(newsup), parentbounds.sup));
545  	   }
546  	
547  	   /* if the new interval is invalid, the old one was already invalid */
548  	   if( newinf <= newsup )
549  	      SCIPintervalSetBounds(newbounds, newinf, newsup);
550  	   else
551  	      SCIPintervalSetEmpty(newbounds);
552  	
553  	   return SCIP_OKAY;
554  	}
555  	
556  	/** helper function to compute coefficients and constant term of a linear estimator at a given point
557  	 *
558  	 *  The function will try to compute the following estimators in that order:
559  	 *  - soltangent: tangent at specified refpoint
560  	 *  - secant: secant between the points (lb,sin(lb)) and (ub,sin(ub))
561  	 *  - left secant: secant between lower bound and some point of the graph
562  	 *  - right secant: secant between upper bound and some point of the graph
563  	 *
564  	 *  They are ordered such that a successful computation for one of them cannot be improved by following ones in terms
565  	 *  of value at the reference point.
566  	 */
567  	static
568  	SCIP_Bool computeEstimatorsTrig(
569  	   SCIP*                 scip,               /**< SCIP data structure */
570  	   SCIP_EXPR*            expr,               /**< sin or cos expression */
571  	   SCIP_Real*            lincoef,            /**< buffer to store the linear coefficient */
572  	   SCIP_Real*            linconst,           /**< buffer to store the constant term */
573  	   SCIP_Real             refpoint,           /**< point at which to underestimate (can be SCIP_INVALID) */
574  	   SCIP_Real             childlb,            /**< lower bound of child variable */
575  	   SCIP_Real             childub,            /**< upper bound of child variable */
576  	   SCIP_Bool             underestimate       /**< whether the estimator should be underestimating */
577  	   )
578  	{
579  	   SCIP_Bool success;
580  	   SCIP_Bool iscos;
581  	
582  	   assert(scip != NULL);
583  	   assert(expr != NULL);
584  	   assert(SCIPexprGetNChildren(expr) == 1);
585  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "sin") == 0
586  	         || strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") == 0);
587  	   assert(SCIPisLE(scip, childlb, childub));
588  	
589  	   /* if child is essentially constant, then there should be no point in estimation */
590  	   if( SCIPisEQ(scip, childlb, childub) ) /* @todo maybe return a constant estimator? */
591  	      return FALSE;
592  	
593  	   iscos = strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") == 0;
594  	
595  	   /* for cos expressions, the bounds have to be shifted before and after computation */
596  	   if( iscos )
597  	   {
598  	      childlb += M_PI_2;
599  	      childub += M_PI_2;
600  	      refpoint += M_PI_2;
601  	   }
602  	
603  	   if( !underestimate )
604  	   {
605  	      SCIP_Real tmp = childlb;
606  	      childlb = -childub;
607  	      childub = -tmp;
608  	      refpoint *= -1;
609  	   }
610  	
611  	   /* try out tangent at solution point */
612  	   success = computeSolTangentSin(scip, lincoef, linconst, childlb, childub, refpoint);
613  	
614  	   /* otherwise, try out secant */
615  	   if( !success )
616  	      success = computeSecantSin(scip, lincoef, linconst, childlb, childub);
617  	
618  	   /* otherwise, try left secant */
619  	   if( !success )
620  	      success = computeLeftSecantSin(scip, lincoef, linconst, childlb, childub);
621  	
622  	   /* otherwise, try right secant */
623  	   if( !success )
624  	      success = computeRightSecantSin(scip, lincoef, linconst, childlb, childub);
625  	
626  	   if( !success )
627  	      return FALSE;
628  	
629  	   /* for overestimators, mirror back */
630  	   if( !underestimate )
631  	      (*linconst) *= -1.0;
632  	
633  	   /* for cos expressions, shift back */
634  	   if( iscos )
635  	      (*linconst) += (*lincoef) * M_PI_2;
636  	
637  	   return TRUE;
638  	}
639  	
640  	/** helper function to create initial cuts for sine and cosine separation
641  	 *
642  	 *  The following 5 cuts can be generated:
643  	 *  - secant: secant between the bounds (lb,sin(lb)) and (ub,sin(ub))
644  	 *  - left/right secant: secant between lower/upper bound and some point of the graph
645  	 *  - left/right tangent: tangents at the lower/upper bounds
646  	 */
647  	static
648  	SCIP_RETCODE computeInitialCutsTrig(
649  	   SCIP*                 scip,               /**< SCIP data structure */
650  	   SCIP_EXPR*            expr,               /**< sin or cos expression */
651  	   SCIP_Real             childlb,            /**< lower bound of child variable */
652  	   SCIP_Real             childub,            /**< upper bound of child variable */
653  	   SCIP_Bool             underestimate,      /**< whether the cuts should be underestimating */
654  	   SCIP_Real**           coefs,              /**< buffer to store coefficients of computed estimators */
655  	   SCIP_Real*            constant,           /**< buffer to store constant of computed estimators */
656  	   int*                  nreturned           /**< buffer to store number of estimators that have been computed */
657  	   )
658  	{
659  	   SCIP_Bool iscos;
660  	   int i;
661  	
662  	   assert(scip != NULL);
663  	   assert(expr != NULL);
664  	   assert(SCIPexprGetNChildren(expr) == 1);
665  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "sin") == 0 || strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") == 0);
666  	   assert(SCIPisLE(scip, childlb, childub));
667  	
668  	   /* caller must ensure that variable is not already fixed */
669  	   assert(!SCIPisEQ(scip, childlb, childub));
670  	
671  	   *nreturned = 0;
672  	
673  	   /* for cos expressions, the bounds have to be shifted before and after computation */
674  	   iscos = strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") == 0;
675  	   if( iscos )
676  	   {
677  	      childlb += M_PI_2;
678  	      childub += M_PI_2;
679  	   }
680  	
681  	   /*
682  	    * Compute all initial cuts
683  	    * For each linear equation z = a*x + b with bounds [lb,ub] the parameters can be computed by:
684  	    *
685  	    * a = cos(x^)    and     b = sin(x^) - a * x^        where x^ is any known point in [lb,ub]
686  	    *
687  	    * and the resulting cut is       a*x + b <=/>= z           depending on over-/underestimation
688  	    */
689  	
690  	   if( ! underestimate )
691  	   {
692  	      SCIP_Real aux;
693  	      aux = childlb;
694  	      childlb = -childub;
695  	      childub = -aux;
696  	   }
697  	
698  	   /* if we can generate a secant between the bounds, then we have convex (concave) hull */
699  	   if( computeSecantSin(scip, coefs[*nreturned], &constant[*nreturned], childlb, childub) )
700  	      (*nreturned)++;
701  	   else
702  	   {
703  	      /* try generating a secant between lb (ub) and some point < ub (> lb); otherwise try with tangent at lb (ub)*/
704  	      if( computeLeftSecantSin(scip, coefs[*nreturned], &constant[*nreturned], childlb, childub) )
705  	         (*nreturned)++;
706  	      else if( computeLeftTangentSin(scip, coefs[*nreturned], &constant[*nreturned], childlb) )
707  	         (*nreturned)++;
708  	
709  	      /* try generating a secant between ub (lb) and some point > lb (< ub); otherwise try with tangent at ub (lb)*/
710  	      if( computeRightSecantSin(scip, coefs[*nreturned], &constant[*nreturned], childlb, childub) )
711  	         (*nreturned)++;
712  	      else if( computeRightTangentSin(scip, coefs[*nreturned], &constant[*nreturned], childub) )
713  	         (*nreturned)++;
714  	   }
715  	
716  	   /* for cos expressions, the estimator needs to be shifted back to match original bounds */
717  	   for( i = 0; i < *nreturned; ++i )
718  	   {
719  	      if( ! underestimate )
720  	         constant[i] *= -1.0;
721  	
722  	      if( iscos)
723  	      {
724  	         constant[i] += coefs[i][0] * M_PI_2;
725  	      }
726  	   }
727  	
728  	   return SCIP_OKAY;
729  	}
730  	
731  	/* helper function that computes the curvature of a sine expression for given bounds and curvature of child */
732  	static
733  	SCIP_EXPRCURV computeCurvatureSin(
734  	   SCIP_EXPRCURV         childcurvature,     /**< curvature of child */
735  	   SCIP_Real             lb,                 /**< lower bound of child */
736  	   SCIP_Real             ub                  /**< upper bound of child */
737  	   )
738  	{
739  	   SCIP_Real lbsin = sin(lb);
740  	   SCIP_Real ubsin = sin(ub);
741  	   SCIP_Real lbcos = cos(lb);
742  	   SCIP_Real ubcos = cos(ub);
743  	
744  	   /* curvature can only be determined if bounds lie within one bay*/
745  	   if( (ub - lb <= M_PI) && (lbsin * ubsin >= 0.0) )
746  	   {
747  	      /* special case that both sin(ub) and sin(lb) are 0 (i.e. ub - lb = pi) */
748  	      if( lbsin == 0.0 && ubsin == 0.0 )
749  	      {
750  	         if( childcurvature == SCIP_EXPRCURV_LINEAR )
751  	            return (fmod(lb, 2.0*M_PI) == 0.0) ? SCIP_EXPRCURV_CONCAVE : SCIP_EXPRCURV_CONVEX;
752  	      }
753  	
754  	      /* if sine is monotone on the interval, the curvature depends on the child curvature and on the segment */
755  	      else if( lbcos * ubcos >= 0.0 )
756  	      {
757  	         /* on [0, pi/2], sine is concave iff child is concave */
758  	         if( lbsin >= 0.0 && lbcos >= 0.0 && ((int)(childcurvature & SCIP_EXPRCURV_CONCAVE) != 0))
759  	            return SCIP_EXPRCURV_CONCAVE;
760  	
761  	         /* on [pi/2, pi], sine is concave iff child is convex */
762  	         if( lbsin >= 0.0 && lbcos <= 0.0 && ((int)(childcurvature & SCIP_EXPRCURV_CONVEX) != 0))
763  	            return SCIP_EXPRCURV_CONCAVE;
764  	
765  	         /* on [pi, 3pi/2], sine is convex iff child is concave */
766  	         if( lbsin <= 0.0 && lbcos <= 0.0 && ((int)(childcurvature & SCIP_EXPRCURV_CONCAVE) != 0))
767  	            return SCIP_EXPRCURV_CONVEX;
768  	
769  	         /* on [3pi/2, 2pi], sine is convex iff child is convex */
770  	         if( lbsin <= 0.0 && lbcos >= 0.0 && ((int)(childcurvature & SCIP_EXPRCURV_CONVEX) != 0))
771  	            return SCIP_EXPRCURV_CONVEX;
772  	      }
773  	
774  	      /* otherwise, we can only say something if the child is linear */
775  	      else if( childcurvature == SCIP_EXPRCURV_LINEAR )
776  	         return (lbsin >= 0.0 && ubsin >= 0.0) ? SCIP_EXPRCURV_CONCAVE : SCIP_EXPRCURV_CONVEX;
777  	   }
778  	
779  	   return SCIP_EXPRCURV_UNKNOWN;
780  	}
781  	
782  	/*
783  	 * Callback methods of expression handler
784  	 */
785  	
786  	/** expression handler copy callback */
787  	static
788  	SCIP_DECL_EXPRCOPYHDLR(copyhdlrSin)
789  	{  /*lint --e{715}*/
790  	   SCIP_CALL( SCIPincludeExprhdlrSin(scip) );
791  	
792  	   return SCIP_OKAY;
793  	}
794  	
795  	/** simplifies a sine expression
796  	 *
797  	 * Evaluates the sine value function when its child is a value expression.
798  	 *
799  	 * TODO: add further simplifications
800  	 */
801  	static
802  	SCIP_DECL_EXPRSIMPLIFY(simplifySin)
803  	{  /*lint --e{715}*/
804  	   SCIP_EXPR* child;
805  	
806  	   assert(scip != NULL);
807  	   assert(expr != NULL);
808  	   assert(simplifiedexpr != NULL);
809  	   assert(SCIPexprGetNChildren(expr) == 1);
810  	
811  	   child = SCIPexprGetChildren(expr)[0];
812  	   assert(child != NULL);
813  	
814  	   /* check for value expression */
815  	   if( SCIPisExprValue(scip, child) )
816  	   {
817  	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, sin(SCIPgetValueExprValue(child)), ownercreate,
818  	            ownercreatedata) );
819  	   }
820  	   else
821  	   {
822  	      *simplifiedexpr = expr;
823  	
824  	      /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
825  	      SCIPcaptureExpr(*simplifiedexpr);
826  	   }
827  	
828  	   return SCIP_OKAY;
829  	}
830  	
831  	/** expression parse callback */
832  	static
833  	SCIP_DECL_EXPRPARSE(parseSin)
834  	{  /*lint --e{715}*/
835  	   SCIP_EXPR* childexpr;
836  	
837  	   assert(expr != NULL);
838  	
839  	   /* parse child expression from remaining string */
840  	   SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
841  	   assert(childexpr != NULL);
842  	
843  	   /* create sine expression */
844  	   SCIP_CALL( SCIPcreateExprSin(scip, expr, childexpr, ownercreate, ownercreatedata) );
845  	   assert(*expr != NULL);
846  	
847  	   /* release child expression since it has been captured by the sine expression */
848  	   SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
849  	
850  	   *success = TRUE;
851  	
852  	   return SCIP_OKAY;
853  	}
854  	
855  	/** expression (point-) evaluation callback */
856  	static
857  	SCIP_DECL_EXPREVAL(evalSin)
858  	{  /*lint --e{715}*/
859  	   assert(expr != NULL);
860  	   assert(SCIPexprGetNChildren(expr) == 1);
861  	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
862  	
863  	   *val = sin(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
864  	
865  	   return SCIP_OKAY;
866  	}
867  	
868  	/** expression derivative evaluation callback */
869  	static
870  	SCIP_DECL_EXPRBWDIFF(bwdiffSin)
871  	{  /*lint --e{715}*/
872  	   SCIP_EXPR* child;
873  	
874  	   assert(expr != NULL);
875  	   assert(childidx == 0);
876  	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
877  	
878  	   child = SCIPexprGetChildren(expr)[0];
879  	   assert(child != NULL);
880  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
881  	
882  	   *val = cos(SCIPexprGetEvalValue(child));
883  	
884  	   return SCIP_OKAY;
885  	}
886  	
887  	/** derivative evaluation callback
888  	 *
889  	 * Computes <gradient, children.dot>, that is, cos(child) dot(child).
890  	 */
891  	static
892  	SCIP_DECL_EXPRFWDIFF(fwdiffSin)
893  	{  /*lint --e{715}*/
894  	   SCIP_EXPR* child;
895  	
896  	   assert(expr != NULL);
897  	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
898  	
899  	   child = SCIPexprGetChildren(expr)[0];
900  	   assert(child != NULL);
901  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
902  	   assert(SCIPexprGetDot(child) != SCIP_INVALID); /*lint !e777*/
903  	
904  	   *dot = cos(SCIPexprGetEvalValue(child)) * SCIPexprGetDot(child);
905  	
906  	   return SCIP_OKAY;
907  	}
908  	
909  	/** expression backward forward derivative evaluation callback
910  	 *
911  	 * Computes partial/partial child ( <gradient, children.dot> ), that is, -sin(child) dot(child).
912  	 */
913  	static
914  	SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSin)
915  	{  /*lint --e{715}*/
916  	   SCIP_EXPR* child;
917  	
918  	   assert(expr != NULL);
919  	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
920  	   assert(childidx == 0);
921  	
922  	   child = SCIPexprGetChildren(expr)[0];
923  	   assert(child != NULL);
924  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
925  	   assert(SCIPexprGetDot(child) != SCIP_INVALID); /*lint !e777*/
926  	
927  	   *bardot = -sin(SCIPexprGetEvalValue(child)) * SCIPexprGetDot(child);
928  	
929  	   return SCIP_OKAY;
930  	}
931  	
932  	/** expression interval evaluation callback */
933  	static
934  	SCIP_DECL_EXPRINTEVAL(intevalSin)
935  	{  /*lint --e{715}*/
936  	   SCIP_INTERVAL childinterval;
937  	
938  	   assert(expr != NULL);
939  	   assert(SCIPexprGetNChildren(expr) == 1);
940  	
941  	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
942  	
943  	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
944  	      SCIPintervalSetEmpty(interval);
945  	   else
946  	      SCIPintervalSin(SCIP_INTERVAL_INFINITY, interval, childinterval);
947  	
948  	   return SCIP_OKAY;
949  	}
950  	
951  	/** separation initialization callback */
952  	static
953  	SCIP_DECL_EXPRINITESTIMATES(initEstimatesSin)
954  	{  /*lint --e{715}*/
955  	   SCIP_Real childlb;
956  	   SCIP_Real childub;
957  	
958  	   childlb = bounds[0].inf;
959  	   childub = bounds[0].sup;
960  	
961  	   /* no need for cut if child is fixed */
962  	   if( SCIPisRelEQ(scip, childlb, childub) )
963  	      return SCIP_OKAY;
964  	
965  	   /* compute cuts */
966  	   SCIP_CALL( computeInitialCutsTrig(scip, expr, childlb, childub, ! overestimate, coefs, constant, nreturned) );
967  	
968  	   return SCIP_OKAY;
969  	}
970  	
971  	/** expression estimator callback */
972  	static
973  	SCIP_DECL_EXPRESTIMATE(estimateSin)
974  	{  /*lint --e{715}*/
975  	   assert(scip != NULL);
976  	   assert(expr != NULL);
977  	   assert(SCIPexprGetNChildren(expr) == 1);
978  	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), SINEXPRHDLR_NAME) == 0);
979  	   assert(coefs != NULL);
980  	   assert(constant != NULL);
981  	   assert(islocal != NULL);
982  	   assert(branchcand != NULL);
983  	   assert(*branchcand == TRUE);
984  	   assert(success != NULL);
985  	
986  	   *success = computeEstimatorsTrig(scip, expr, coefs, constant, refpoint[0], localbounds[0].inf,
987  	         localbounds[0].sup, ! overestimate);
988  	   *islocal = TRUE;  /* TODO there are cases where cuts would be globally valid */
989  	
990  	   return SCIP_OKAY;
991  	}
992  	
993  	/** expression reverse propagation callback */
994  	static
995  	SCIP_DECL_EXPRREVERSEPROP(reversepropSin)
996  	{  /*lint --e{715}*/
997  	   assert(scip != NULL);
998  	   assert(expr != NULL);
999  	   assert(SCIPexprGetNChildren(expr) == 1);
1000 	   assert(SCIPintervalGetInf(bounds) >= -1.0);
1001 	   assert(SCIPintervalGetSup(bounds) <= 1.0);
1002 	
1003 	   /* compute the new child interval */
1004 	   SCIP_CALL( computeRevPropIntervalSin(scip, bounds, childrenbounds[0], childrenbounds) );
1005 	
1006 	   return SCIP_OKAY;
1007 	}
1008 	
1009 	/** sine hash callback */
1010 	static
1011 	SCIP_DECL_EXPRHASH(hashSin)
1012 	{  /*lint --e{715}*/
1013 	   assert(scip != NULL);
1014 	   assert(expr != NULL);
1015 	   assert(SCIPexprGetNChildren(expr) == 1);
1016 	   assert(hashkey != NULL);
1017 	   assert(childrenhashes != NULL);
1018 	
1019 	   *hashkey = SINEXPRHDLR_HASHKEY;
1020 	   *hashkey ^= childrenhashes[0];
1021 	
1022 	   return SCIP_OKAY;
1023 	}
1024 	
1025 	/** expression curvature detection callback */
1026 	static
1027 	SCIP_DECL_EXPRCURVATURE(curvatureSin)
1028 	{  /*lint --e{715}*/
1029 	   SCIP_EXPR* child;
1030 	   SCIP_INTERVAL childinterval;
1031 	
1032 	   assert(scip != NULL);
1033 	   assert(expr != NULL);
1034 	   assert(childcurv != NULL);
1035 	   assert(success != NULL);
1036 	   assert(SCIPexprGetNChildren(expr) == 1);
1037 	
1038 	   child = SCIPexprGetChildren(expr)[0];
1039 	   assert(child != NULL);
1040 	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
1041 	   childinterval = SCIPexprGetActivity(child);
1042 	
1043 	   /* TODO rewrite SCIPcomputeCurvatureSin so it provides the reverse operation */
1044 	   *success = TRUE;
1045 	   if( computeCurvatureSin(SCIP_EXPRCURV_CONVEX, childinterval.inf, childinterval.sup) == exprcurvature )
1046 	      *childcurv = SCIP_EXPRCURV_CONVEX;
1047 	   else if( computeCurvatureSin(SCIP_EXPRCURV_CONCAVE, childinterval.inf, childinterval.sup) == exprcurvature )
1048 	      *childcurv = SCIP_EXPRCURV_CONCAVE;
1049 	   if( computeCurvatureSin(SCIP_EXPRCURV_LINEAR, childinterval.inf, childinterval.sup) == exprcurvature )
1050 	      *childcurv = SCIP_EXPRCURV_LINEAR;
1051 	   else
1052 	      *success = FALSE;
1053 	
1054 	   return SCIP_OKAY;
1055 	}
1056 	
1057 	/** expression monotonicity detection callback */
1058 	static
1059 	SCIP_DECL_EXPRMONOTONICITY(monotonicitySin)
1060 	{  /*lint --e{715}*/
1061 	   SCIP_INTERVAL interval;
1062 	   SCIP_Real inf;
1063 	   SCIP_Real sup;
1064 	   int k;
1065 	
1066 	   assert(scip != NULL);
1067 	   assert(expr != NULL);
1068 	   assert(result != NULL);
1069 	   assert(childidx == 0);
1070 	
1071 	   assert(SCIPexprGetChildren(expr)[0] != NULL);
1072 	   SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(expr)[0]) );
1073 	   interval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
1074 	
1075 	   *result = SCIP_MONOTONE_UNKNOWN;
1076 	   inf = SCIPintervalGetInf(interval);
1077 	   sup = SCIPintervalGetSup(interval);
1078 	
1079 	   /* expression is not monotone because the interval is too large */
1080 	   if( SCIPisGT(scip, sup - inf, M_PI) )
1081 	      return SCIP_OKAY;
1082 	
1083 	   /* compute k s.t. PI * (2k+1) / 2 <= interval.inf <= PI * (2k+3) / 2 */
1084 	   k = (int)floor(inf/M_PI - 0.5);
1085 	   assert(SCIPisLE(scip, M_PI * (2.0*k + 1.0) / 2.0, inf));
1086 	   assert(SCIPisGE(scip, M_PI * (2.0*k + 3.0) / 2.0, inf));
1087 	
1088 	   /* check whether [inf,sup] are contained in an interval for which the sine function is monotone */
1089 	   if( SCIPisLE(scip, sup, M_PI * (2.0*k + 3.0) / 2.0) )
1090 	      *result = ((k % 2 + 2) % 2) == 1 ? SCIP_MONOTONE_INC : SCIP_MONOTONE_DEC;
1091 	
1092 	   return SCIP_OKAY;
1093 	}
1094 	
1095 	
1096 	/** expression handler copy callback */
1097 	static
1098 	SCIP_DECL_EXPRCOPYHDLR(copyhdlrCos)
1099 	{  /*lint --e{715}*/
1100 	   SCIP_CALL( SCIPincludeExprhdlrCos(scip) );
1101 	
1102 	   return SCIP_OKAY;
1103 	}
1104 	
1105 	/** simplifies a cosine expression
1106 	 *
1107 	 * Evaluates the cosine value function when its child is a value expression.
1108 	 *
1109 	 * TODO: add further simplifications
1110 	 */
1111 	static
1112 	SCIP_DECL_EXPRSIMPLIFY(simplifyCos)
1113 	{  /*lint --e{715}*/
1114 	   SCIP_EXPR* child;
1115 	
1116 	   assert(scip != NULL);
1117 	   assert(expr != NULL);
1118 	   assert(simplifiedexpr != NULL);
1119 	   assert(SCIPexprGetNChildren(expr) == 1);
1120 	
1121 	   child = SCIPexprGetChildren(expr)[0];
1122 	   assert(child != NULL);
1123 	
1124 	   /* check for value expression */
1125 	   if( SCIPisExprValue(scip, child) )
1126 	   {
1127 	      SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, cos(SCIPgetValueExprValue(child)), ownercreate,
1128 	               ownercreatedata) );
1129 	   }
1130 	   else
1131 	   {
1132 	      *simplifiedexpr = expr;
1133 	
1134 	      /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */
1135 	      SCIPcaptureExpr(*simplifiedexpr);
1136 	   }
1137 	
1138 	   return SCIP_OKAY;
1139 	}
1140 	
1141 	/** expression parse callback */
1142 	static
1143 	SCIP_DECL_EXPRPARSE(parseCos)
1144 	{  /*lint --e{715}*/
1145 	   SCIP_EXPR* childexpr;
1146 	
1147 	   assert(expr != NULL);
1148 	
1149 	   /* parse child expression from remaining string */
1150 	   SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) );
1151 	   assert(childexpr != NULL);
1152 	
1153 	   /* create cosine expression */
1154 	   SCIP_CALL( SCIPcreateExprCos(scip, expr, childexpr, ownercreate, ownercreatedata) );
1155 	   assert(*expr != NULL);
1156 	
1157 	   /* release child expression since it has been captured by the cosine expression */
1158 	   SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) );
1159 	
1160 	   *success = TRUE;
1161 	
1162 	   return SCIP_OKAY;
1163 	}
1164 	
1165 	/** expression (point-) evaluation callback */
1166 	static
1167 	SCIP_DECL_EXPREVAL(evalCos)
1168 	{  /*lint --e{715}*/
1169 	   assert(expr != NULL);
1170 	   assert(SCIPexprGetNChildren(expr) == 1);
1171 	   assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/
1172 	
1173 	   *val = cos(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]));
1174 	
1175 	   return SCIP_OKAY;
1176 	}
1177 	
1178 	/** expression derivative evaluation callback */
1179 	static
1180 	SCIP_DECL_EXPRBWDIFF(bwdiffCos)
1181 	{  /*lint --e{715}*/
1182 	   SCIP_EXPR* child;
1183 	
1184 	   assert(expr != NULL);
1185 	   assert(childidx == 0);
1186 	   assert(SCIPexprGetEvalValue(expr) != SCIP_INVALID); /*lint !e777*/
1187 	
1188 	   child = SCIPexprGetChildren(expr)[0];
1189 	   assert(child != NULL);
1190 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(child)), "val") != 0);
1191 	
1192 	   *val = -sin(SCIPexprGetEvalValue(child));
1193 	
1194 	   return SCIP_OKAY;
1195 	}
1196 	
1197 	/** expression interval evaluation callback */
1198 	static
1199 	SCIP_DECL_EXPRINTEVAL(intevalCos)
1200 	{  /*lint --e{715}*/
1201 	   SCIP_INTERVAL childinterval;
1202 	
1203 	   assert(expr != NULL);
1204 	   assert(SCIPexprGetNChildren(expr) == 1);
1205 	
1206 	   childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
1207 	
1208 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) )
1209 	      SCIPintervalSetEmpty(interval);
1210 	   else
1211 	      SCIPintervalCos(SCIP_INTERVAL_INFINITY, interval, childinterval);
1212 	
1213 	   return SCIP_OKAY;
1214 	}
1215 	
1216 	/** separation initialization callback */
1217 	static
1218 	SCIP_DECL_EXPRINITESTIMATES(initEstimatesCos)
1219 	{
1220 	   SCIP_Real childlb;
1221 	   SCIP_Real childub;
1222 	
1223 	   childlb = bounds[0].inf;
1224 	   childub = bounds[0].sup;
1225 	
1226 	   /* no need for cut if child is fixed */
1227 	   if( SCIPisRelEQ(scip, childlb, childub) )
1228 	      return SCIP_OKAY;
1229 	
1230 	   /* compute cuts */
1231 	   SCIP_CALL( computeInitialCutsTrig(scip, expr, childlb, childub, ! overestimate, coefs, constant, nreturned) );
1232 	
1233 	   return SCIP_OKAY;
1234 	}
1235 	
1236 	/** expression estimator callback */
1237 	static
1238 	SCIP_DECL_EXPRESTIMATE(estimateCos)
1239 	{  /*lint --e{715}*/
1240 	   assert(scip != NULL);
1241 	   assert(expr != NULL);
1242 	   assert(SCIPexprGetNChildren(expr) == 1);
1243 	   assert(strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), COSEXPRHDLR_NAME) == 0);
1244 	   assert(coefs != NULL);
1245 	   assert(constant != NULL);
1246 	   assert(islocal != NULL);
1247 	   assert(branchcand != NULL);
1248 	   assert(*branchcand == TRUE);
1249 	   assert(success != NULL);
1250 	
1251 	   *success = computeEstimatorsTrig(scip, expr, coefs, constant, refpoint[0], localbounds[0].inf,
1252 	         localbounds[0].sup, ! overestimate);
1253 	   *islocal = TRUE;  /* TODO there are cases where cuts would be globally valid */
1254 	
1255 	   return SCIP_OKAY;
1256 	}
1257 	
1258 	/** expression reverse propagation callback */
1259 	static
1260 	SCIP_DECL_EXPRREVERSEPROP(reversepropCos)
1261 	{  /*lint --e{715}*/
1262 	   SCIP_INTERVAL newbounds;
1263 	
1264 	   assert(scip != NULL);
1265 	   assert(expr != NULL);
1266 	   assert(SCIPexprGetNChildren(expr) == 1);
1267 	   /* bounds should have been intersected with activity, which is within [-1,1] */
1268 	   assert(SCIPintervalGetInf(bounds) >= -1.0);
1269 	   assert(SCIPintervalGetSup(bounds) <= 1.0);
1270 	
1271 	   /* get the child interval */
1272 	   newbounds = childrenbounds[0];
1273 	
1274 	   /* shift child interval to match sine */
1275 	   SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &newbounds, newbounds, M_PI_2);  /* TODO use bounds on Pi/2 instead of approximation of Pi/2 */
1276 	
1277 	   /* compute the new child interval */
1278 	   SCIP_CALL( computeRevPropIntervalSin(scip, bounds, newbounds, &newbounds) );
1279 	
1280 	   if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, newbounds) )
1281 	   {
1282 	      *infeasible = TRUE;
1283 	      return SCIP_OKAY;
1284 	   }
1285 	
1286 	   /* shift the new interval back */
1287 	   SCIPintervalAddScalar(SCIP_INTERVAL_INFINITY, &childrenbounds[0], newbounds, -M_PI_2);  /* TODO use bounds on Pi/2 instead of approximation of Pi/2 */
1288 	
1289 	   return SCIP_OKAY;
1290 	}
1291 	
1292 	/** cosine hash callback */
1293 	static
1294 	SCIP_DECL_EXPRHASH(hashCos)
1295 	{  /*lint --e{715}*/
1296 	   assert(scip != NULL);
1297 	   assert(expr != NULL);
1298 	   assert(SCIPexprGetNChildren(expr) == 1);
1299 	   assert(hashkey != NULL);
1300 	   assert(childrenhashes != NULL);
1301 	
1302 	   *hashkey = COSEXPRHDLR_HASHKEY;
1303 	   *hashkey ^= childrenhashes[0];
1304 	
1305 	   return SCIP_OKAY;
1306 	}
1307 	
1308 	/** expression curvature detection callback */
1309 	static
1310 	SCIP_DECL_EXPRCURVATURE(curvatureCos)
1311 	{  /*lint --e{715}*/
1312 	   SCIP_EXPR* child;
1313 	   SCIP_INTERVAL childinterval;
1314 	
1315 	   assert(scip != NULL);
1316 	   assert(expr != NULL);
1317 	   assert(exprcurvature != SCIP_EXPRCURV_UNKNOWN);
1318 	   assert(childcurv != NULL);
1319 	   assert(success != NULL);
1320 	   assert(SCIPexprGetNChildren(expr) == 1);
1321 	
1322 	   child = SCIPexprGetChildren(expr)[0];
1323 	   assert(child != NULL);
1324 	   SCIP_CALL( SCIPevalExprActivity(scip, child) );
1325 	   childinterval = SCIPexprGetActivity(child);
1326 	
1327 	   /* TODO rewrite SCIPcomputeCurvatureSin so it provides the reverse operation */
1328 	   *success = TRUE;
1329 	   if( computeCurvatureSin(SCIP_EXPRCURV_CONCAVE, childinterval.inf + M_PI_2, childinterval.sup + M_PI_2) == exprcurvature )
1330 	      *childcurv = SCIP_EXPRCURV_CONCAVE;
1331 	   else if( computeCurvatureSin(SCIP_EXPRCURV_CONVEX, childinterval.inf + M_PI_2, childinterval.sup + M_PI_2) == exprcurvature )
1332 	      *childcurv = SCIP_EXPRCURV_CONVEX;
1333 	   else if( computeCurvatureSin(SCIP_EXPRCURV_LINEAR, childinterval.inf + M_PI_2, childinterval.sup + M_PI_2) == exprcurvature )
1334 	      *childcurv = SCIP_EXPRCURV_LINEAR;
1335 	   else
1336 	      *success = FALSE;
1337 	
1338 	   return SCIP_OKAY;
1339 	}
1340 	
1341 	/** expression monotonicity detection callback */
1342 	static
1343 	SCIP_DECL_EXPRMONOTONICITY(monotonicityCos)
1344 	{  /*lint --e{715}*/
1345 	   SCIP_INTERVAL interval;
1346 	   SCIP_Real inf;
1347 	   SCIP_Real sup;
1348 	   int k;
1349 	
1350 	   assert(scip != NULL);
1351 	   assert(expr != NULL);
1352 	   assert(result != NULL);
1353 	   assert(childidx == 0);
1354 	
1355 	   assert(SCIPexprGetChildren(expr)[0] != NULL);
1356 	   SCIP_CALL( SCIPevalExprActivity(scip, SCIPexprGetChildren(expr)[0]) );
1357 	   interval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]);
1358 	
1359 	   *result = SCIP_MONOTONE_UNKNOWN;
1360 	   inf = SCIPintervalGetInf(interval);
1361 	   sup = SCIPintervalGetSup(interval);
1362 	
1363 	   /* expression is not monotone because the interval is too large */
1364 	   if( SCIPisGT(scip, sup - inf, M_PI) )
1365 	      return SCIP_OKAY;
1366 	
1367 	   /* compute k s.t. PI * k <= interval.inf <= PI * (k+1) */
1368 	   k = (int)floor(inf/M_PI);
1369 	   assert(SCIPisLE(scip, M_PI * k, inf));
1370 	   assert(SCIPisGE(scip, M_PI * (k+1), inf));
1371 	
1372 	   /* check whether [inf,sup] are contained in an interval for which the cosine function is monotone */
1373 	   if( SCIPisLE(scip, sup, M_PI * (k+1)) )
1374 	      *result = ((k % 2 + 2) % 2) == 0 ? SCIP_MONOTONE_DEC : SCIP_MONOTONE_INC;
1375 	
1376 	   return SCIP_OKAY;
1377 	}
1378 	
1379 	/** creates the handler for sin expressions and includes it into SCIP */
1380 	SCIP_RETCODE SCIPincludeExprhdlrSin(
1381 	   SCIP*                 scip                /**< SCIP data structure */
1382 	   )
1383 	{
1384 	   SCIP_EXPRHDLR* exprhdlr;
1385 	
1386 	   /* include expression handler */
1387 	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, SINEXPRHDLR_NAME, SINEXPRHDLR_DESC, SINEXPRHDLR_PRECEDENCE, evalSin, NULL) );
1388 	   assert(exprhdlr != NULL);
1389 	
1390 	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrSin, NULL);
1391 	   SCIPexprhdlrSetSimplify(exprhdlr, simplifySin);
1392 	   SCIPexprhdlrSetParse(exprhdlr, parseSin);
1393 	   SCIPexprhdlrSetIntEval(exprhdlr, intevalSin);
1394 	   SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesSin, estimateSin);
1395 	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropSin);
1396 	   SCIPexprhdlrSetHash(exprhdlr, hashSin);
1397 	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffSin, fwdiffSin, bwfwdiffSin);
1398 	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureSin);
1399 	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicitySin);
1400 	
1401 	   return SCIP_OKAY;
1402 	}
1403 	
1404 	/** creates the handler for cos expressions and includes it SCIP */
1405 	SCIP_RETCODE SCIPincludeExprhdlrCos(
1406 	   SCIP*                 scip                /**< SCIP data structure */
1407 	   )
1408 	{
1409 	   SCIP_EXPRHDLR* exprhdlr;
1410 	
1411 	   /* include expression handler */
1412 	   SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, COSEXPRHDLR_NAME, COSEXPRHDLR_DESC, COSEXPRHDLR_PRECEDENCE, evalCos, NULL) );
1413 	   assert(exprhdlr != NULL);
1414 	
1415 	   SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrCos, NULL);
1416 	   SCIPexprhdlrSetSimplify(exprhdlr, simplifyCos);
1417 	   SCIPexprhdlrSetParse(exprhdlr, parseCos);
1418 	   SCIPexprhdlrSetIntEval(exprhdlr, intevalCos);
1419 	   SCIPexprhdlrSetEstimate(exprhdlr, initEstimatesCos, estimateCos);
1420 	   SCIPexprhdlrSetReverseProp(exprhdlr, reversepropCos);
1421 	   SCIPexprhdlrSetHash(exprhdlr, hashCos);
1422 	   SCIPexprhdlrSetDiff(exprhdlr, bwdiffCos, NULL, NULL);
1423 	   SCIPexprhdlrSetCurvature(exprhdlr, curvatureCos);
1424 	   SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityCos);
1425 	
1426 	   return SCIP_OKAY;
1427 	}
1428 	
1429 	/** creates a sin expression */
1430 	SCIP_RETCODE SCIPcreateExprSin(
1431 	   SCIP*                 scip,               /**< SCIP data structure */
1432 	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
1433 	   SCIP_EXPR*            child,              /**< single child */
1434 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1435 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
1436 	   )
1437 	{
1438 	   assert(expr != NULL);
1439 	   assert(child != NULL);
1440 	   assert(SCIPfindExprhdlr(scip, SINEXPRHDLR_NAME) != NULL);
1441 	
1442 	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, SINEXPRHDLR_NAME), NULL, 1, &child, ownercreate,
1443 	            ownercreatedata) );
1444 	
1445 	   return SCIP_OKAY;
1446 	}
1447 	
1448 	
1449 	/** creates a cos expression */
1450 	SCIP_RETCODE SCIPcreateExprCos(
1451 	   SCIP*                 scip,               /**< SCIP data structure */
1452 	   SCIP_EXPR**           expr,               /**< pointer where to store expression */
1453 	   SCIP_EXPR*            child,              /**< single child */
1454 	   SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */
1455 	   void*                 ownercreatedata     /**< data to pass to ownercreate */
1456 	   )
1457 	{
1458 	   assert(expr != NULL);
1459 	   assert(child != NULL);
1460 	   assert(SCIPfindExprhdlr(scip, COSEXPRHDLR_NAME) != NULL);
1461 	
1462 	   SCIP_CALL( SCIPcreateExpr(scip, expr, SCIPfindExprhdlr(scip, COSEXPRHDLR_NAME), NULL, 1, &child, ownercreate,
1463 	            ownercreatedata) );
1464 	
1465 	   return SCIP_OKAY;
1466 	}
1467 	
1468 	/** indicates whether expression is of sine-type */  /*lint -e{715}*/
1469 	SCIP_Bool SCIPisExprSin(
1470 	   SCIP*                 scip,               /**< SCIP data structure */
1471 	   SCIP_EXPR*            expr                /**< expression */
1472 	   )
1473 	{  /*lint --e{715}*/
1474 	   assert(expr != NULL);
1475 	
1476 	   return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), SINEXPRHDLR_NAME) == 0;
1477 	}
1478 	
1479 	/** indicates whether expression is of cosine-type */  /*lint -e{715}*/
1480 	SCIP_Bool SCIPisExprCos(
1481 	   SCIP*                 scip,               /**< SCIP data structure */
1482 	   SCIP_EXPR*            expr                /**< expression */
1483 	   )
1484 	{  /*lint --e{715}*/
1485 	   assert(expr != NULL);
1486 	
1487 	   return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), COSEXPRHDLR_NAME) == 0;
1488 	}
1489