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