Scippy

SCIP

Solving Constraint Integer Programs

expr_trig.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (c) 2002-2024 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
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
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
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
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
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
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
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
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,
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
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,
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
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,
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
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
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
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
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
789 { /*lint --e{715}*/
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
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
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
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
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
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
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
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
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
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
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
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
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
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);
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
1099 { /*lint --e{715}*/
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
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
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
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
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
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
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
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
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 
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
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
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
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);
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 */
1381  SCIP* scip /**< SCIP data structure */
1382  )
1383 {
1384  SCIP_EXPRHDLR* exprhdlr;
1385 
1386  /* include expression handler */
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 */
1406  SCIP* scip /**< SCIP data structure */
1407  )
1408 {
1409  SCIP_EXPRHDLR* exprhdlr;
1410 
1411  /* include expression handler */
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 */
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 */
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}*/
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}*/
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 }
#define NULL
Definition: def.h:267
void SCIPintervalCos(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
#define COSEXPRHDLR_HASHKEY
Definition: expr_trig.c:67
#define COSEXPRHDLR_DESC
Definition: expr_trig.c:65
SCIP_Bool SCIPisRelEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1717
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3854
void SCIPexprhdlrSetDiff(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRBWDIFF((*bwdiff)), SCIP_DECL_EXPRFWDIFF((*fwdiff)), SCIP_DECL_EXPRBWFWDIFF((*bwfwdiff)))
Definition: expr.c:473
const char * SCIPexprhdlrGetName(SCIP_EXPRHDLR *exprhdlr)
Definition: expr.c:545
static SCIP_RETCODE computeInitialCutsTrig(SCIP *scip, SCIP_EXPR *expr, SCIP_Real childlb, SCIP_Real childub, SCIP_Bool underestimate, SCIP_Real **coefs, SCIP_Real *constant, int *nreturned)
Definition: expr_trig.c:648
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_DECL_EXPRHASH(hashSin)
Definition: expr_trig.c:1011
void SCIPexprhdlrSetParse(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRPARSE((*parse)))
Definition: expr.c:407
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define FALSE
Definition: def.h:94
static SCIP_Bool computeRightSecantSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:390
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
#define TRUE
Definition: def.h:93
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
#define SINEXPRHDLR_PRECEDENCE
Definition: expr_trig.c:61
static SCIP_Bool computeSolTangentSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub, SCIP_Real solpoint)
Definition: expr_trig.c:228
void SCIPexprhdlrSetHash(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRHASH((*hash)))
Definition: expr.c:451
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition: expr.c:4010
static SCIP_RETCODE computeRevPropIntervalSin(SCIP *scip, SCIP_INTERVAL parentbounds, SCIP_INTERVAL childbounds, SCIP_INTERVAL *newbounds)
Definition: expr_trig.c:472
static SCIP_Bool computeLeftTangentSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb)
Definition: expr_trig.c:166
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void SCIPcaptureExpr(SCIP_EXPR *expr)
Definition: scip_expr.c:1409
#define SCIPdebugMsg
Definition: scip_message.h:78
SCIP_EXPRHDLR * SCIPfindExprhdlr(SCIP *scip, const char *name)
Definition: scip_expr.c:868
SCIP_Real SCIPcalcRootNewton(SCIP_DECL_NEWTONEVAL((*function)), SCIP_DECL_NEWTONEVAL((*derivative)), SCIP_Real *params, int nparams, SCIP_Real x, SCIP_Real eps, int k)
Definition: misc.c:9865
#define COSEXPRHDLR_PRECEDENCE
Definition: expr_trig.c:66
#define COSEXPRHDLR_NAME
Definition: expr_trig.c:64
static SCIP_DECL_EXPRCURVATURE(curvatureSin)
Definition: expr_trig.c:1027
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3864
SCIP_Real inf
Definition: intervalarith.h:55
static SCIP_DECL_EXPRESTIMATE(estimateSin)
Definition: expr_trig.c:973
SCIP_Real SCIPexprGetEvalValue(SCIP_EXPR *expr)
Definition: expr.c:3928
static SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSin)
Definition: expr_trig.c:914
SCIP_RETCODE SCIPcreateExpr(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPRHDLR *exprhdlr, SCIP_EXPRDATA *exprdata, int nchildren, SCIP_EXPR **children, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: scip_expr.c:974
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
#define SINEXPRHDLR_NAME
Definition: expr_trig.c:59
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
static SCIP_DECL_EXPRREVERSEPROP(reversepropSin)
Definition: expr_trig.c:995
SCIP_Bool SCIPisExprValue(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1442
SCIP_RETCODE SCIPincludeExprhdlrCos(SCIP *scip)
Definition: expr_trig.c:1405
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
#define SINEXPRHDLR_DESC
Definition: expr_trig.c:60
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
SCIP_Real SCIPexprGetDot(SCIP_EXPR *expr)
Definition: expr.c:3968
static SCIP_DECL_EXPRMONOTONICITY(monotonicitySin)
Definition: expr_trig.c:1059
#define SCIP_CALL(x)
Definition: def.h:380
SCIP_Real sup
Definition: intervalarith.h:56
SCIP_RETCODE SCIPcreateExprValue(SCIP *scip, SCIP_EXPR **expr, SCIP_Real value, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_value.c:270
SCIP_Bool SCIPisFeasLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void SCIPexprhdlrSetCopyFreeHdlr(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCOPYHDLR((*copyhdlr)), SCIP_DECL_EXPRFREEHDLR((*freehdlr)))
Definition: expr.c:370
static SCIP_EXPRCURV computeCurvatureSin(SCIP_EXPRCURV childcurvature, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:733
#define SCIP_INTERVAL_INFINITY
Definition: def.h:195
static SCIP_DECL_NEWTONEVAL(function1)
Definition: expr_trig.c:82
static SCIP_Bool computeRightTangentSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real ub)
Definition: expr_trig.c:199
static SCIP_Bool computeLeftSecantSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:303
static SCIP_DECL_EXPRINITESTIMATES(initEstimatesSin)
Definition: expr_trig.c:953
#define MAXCHILDABSVAL
Definition: expr_trig.c:69
#define SCIP_Bool
Definition: def.h:91
#define SCIP_DECL_EXPR_OWNERCREATE(x)
Definition: type_expr.h:143
SCIP_EXPRCURV
Definition: type_expr.h:60
void SCIPexprhdlrSetMonotonicity(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRMONOTONICITY((*monotonicity)))
Definition: expr.c:429
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1417
handler for sin expressions
static SCIP_DECL_EXPRSIMPLIFY(simplifySin)
Definition: expr_trig.c:802
void SCIPexprhdlrSetReverseProp(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRREVERSEPROP((*reverseprop)))
Definition: expr.c:510
SCIP_RETCODE SCIPparseExpr(SCIP *scip, SCIP_EXPR **expr, const char *exprstr, const char **finalpos, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: scip_expr.c:1380
#define SINEXPRHDLR_HASHKEY
Definition: expr_trig.c:62
SCIP_EXPRHDLR * SCIPexprGetHdlr(SCIP_EXPR *expr)
Definition: expr.c:3877
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
constant value expression handler
#define M_PI
Definition: pricer_rpa.c:97
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void SCIPexprhdlrSetCurvature(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCURVATURE((*curvature)))
Definition: expr.c:418
static SCIP_DECL_EXPRFWDIFF(fwdiffSin)
Definition: expr_trig.c:892
SCIP_RETCODE SCIPcreateExprSin(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_trig.c:1430
SCIP_RETCODE SCIPcreateExprCos(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_trig.c:1450
SCIP_VAR * a
Definition: circlepacking.c:66
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
static SCIP_Bool computeSecantSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:134
static SCIP_Bool computeEstimatorsTrig(SCIP *scip, SCIP_EXPR *expr, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real refpoint, SCIP_Real childlb, SCIP_Real childub, SCIP_Bool underestimate)
Definition: expr_trig.c:568
#define SCIP_Real
Definition: def.h:173
static SCIP_DECL_EXPREVAL(evalSin)
Definition: expr_trig.c:857
#define NEWTON_PRECISION
Definition: expr_trig.c:71
#define SCIP_INVALID
Definition: def.h:193
SCIP_Real SCIPgetValueExprValue(SCIP_EXPR *expr)
Definition: expr_value.c:294
SCIP_RETCODE SCIPincludeExprhdlrSin(SCIP *scip)
Definition: expr_trig.c:1380
SCIP_Bool SCIPisExprSin(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_trig.c:1469
#define NEWTON_NITERATIONS
Definition: expr_trig.c:70
void SCIPexprhdlrSetEstimate(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRINITESTIMATES((*initestimates)), SCIP_DECL_EXPRESTIMATE((*estimate)))
Definition: expr.c:532
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_DECL_EXPRINTEVAL(intevalSin)
Definition: expr_trig.c:934
static SCIP_DECL_EXPRCOPYHDLR(copyhdlrSin)
Definition: expr_trig.c:788
SCIP_Bool SCIPisExprCos(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_trig.c:1480
static SCIP_DECL_EXPRPARSE(parseSin)
Definition: expr_trig.c:833
SCIP_RETCODE SCIPincludeExprhdlr(SCIP *scip, SCIP_EXPRHDLR **exprhdlr, const char *name, const char *desc, unsigned int precedence, SCIP_DECL_EXPREVAL((*eval)), SCIP_EXPRHDLRDATA *data)
Definition: scip_expr.c:823
void SCIPexprhdlrSetSimplify(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRSIMPLIFY((*simplify)))
Definition: expr.c:499
void SCIPintervalSin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPexprhdlrSetIntEval(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRINTEVAL((*inteval)))
Definition: expr.c:488
#define ABS(x)
Definition: def.h:235
static SCIP_DECL_EXPRBWDIFF(bwdiffSin)
Definition: expr_trig.c:870