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 */
81static
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 */
94static
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 */
107static
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 */
120static
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 */
133static
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 */
165static
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 */
198static
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 */
227static
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 */
302static
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 */
389static
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 */
471static
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 */
567static
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 */
647static
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 */
732static
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))
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))
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))
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))
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
780}
781
782/*
783 * Callback methods of expression handler
784 */
785
786/** expression handler copy callback */
787static
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 */
801static
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 */
832static
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 */
856static
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 */
869static
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 */
891static
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 */
913static
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 */
933static
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 */
952static
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 */
972static
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 */
994static
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 */
1010static
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 */
1026static
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);
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 */
1058static
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 */
1097static
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 */
1111static
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 */
1142static
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 */
1166static
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 */
1179static
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 */
1198static
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 */
1217static
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 */
1237static
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 */
1259static
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 */
1293static
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 */
1309static
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);
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 */
1342static
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);
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);
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}
SCIP_VAR * a
Definition: circlepacking.c:66
#define NULL
Definition: def.h:267
#define SCIP_INVALID
Definition: def.h:193
#define SCIP_INTERVAL_INFINITY
Definition: def.h:195
#define SCIP_Bool
Definition: def.h:91
#define SCIP_Real
Definition: def.h:173
#define ABS(x)
Definition: def.h:235
#define TRUE
Definition: def.h:93
#define FALSE
Definition: def.h:94
#define SCIP_CALL(x)
Definition: def.h:374
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 COSEXPRHDLR_NAME
Definition: expr_trig.c:64
static SCIP_DECL_EXPRINITESTIMATES(initEstimatesSin)
Definition: expr_trig.c:953
static SCIP_Bool computeRightSecantSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:390
#define SINEXPRHDLR_PRECEDENCE
Definition: expr_trig.c:61
static SCIP_DECL_EXPRBWFWDIFF(bwfwdiffSin)
Definition: expr_trig.c:914
static SCIP_RETCODE computeRevPropIntervalSin(SCIP *scip, SCIP_INTERVAL parentbounds, SCIP_INTERVAL childbounds, SCIP_INTERVAL *newbounds)
Definition: expr_trig.c:472
#define SINEXPRHDLR_HASHKEY
Definition: expr_trig.c:62
static SCIP_DECL_EXPRBWDIFF(bwdiffSin)
Definition: expr_trig.c:870
static SCIP_EXPRCURV computeCurvatureSin(SCIP_EXPRCURV childcurvature, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:733
#define COSEXPRHDLR_HASHKEY
Definition: expr_trig.c:67
static SCIP_DECL_EXPRHASH(hashSin)
Definition: expr_trig.c:1011
#define SINEXPRHDLR_DESC
Definition: expr_trig.c:60
static SCIP_DECL_EXPREVAL(evalSin)
Definition: expr_trig.c:857
static SCIP_DECL_EXPRESTIMATE(estimateSin)
Definition: expr_trig.c:973
#define COSEXPRHDLR_DESC
Definition: expr_trig.c:65
static SCIP_DECL_EXPRCOPYHDLR(copyhdlrSin)
Definition: expr_trig.c:788
static SCIP_DECL_NEWTONEVAL(function1)
Definition: expr_trig.c:82
#define NEWTON_NITERATIONS
Definition: expr_trig.c:70
static SCIP_DECL_EXPRSIMPLIFY(simplifySin)
Definition: expr_trig.c:802
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_EXPRMONOTONICITY(monotonicitySin)
Definition: expr_trig.c:1059
#define MAXCHILDABSVAL
Definition: expr_trig.c:69
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
static SCIP_Bool computeSecantSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb, SCIP_Real ub)
Definition: expr_trig.c:134
#define COSEXPRHDLR_PRECEDENCE
Definition: expr_trig.c:66
#define SINEXPRHDLR_NAME
Definition: expr_trig.c:59
#define NEWTON_PRECISION
Definition: expr_trig.c:71
static SCIP_Bool computeRightTangentSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real ub)
Definition: expr_trig.c:199
static SCIP_DECL_EXPRFWDIFF(fwdiffSin)
Definition: expr_trig.c:892
static SCIP_Bool computeLeftTangentSin(SCIP *scip, SCIP_Real *lincoef, SCIP_Real *linconst, SCIP_Real lb)
Definition: expr_trig.c:166
static SCIP_DECL_EXPRPARSE(parseSin)
Definition: expr_trig.c:833
static SCIP_DECL_EXPRREVERSEPROP(reversepropSin)
Definition: expr_trig.c:995
static SCIP_DECL_EXPRCURVATURE(curvatureSin)
Definition: expr_trig.c:1027
static SCIP_DECL_EXPRINTEVAL(intevalSin)
Definition: expr_trig.c:934
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
handler for sin expressions
constant value expression handler
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_Bool SCIPisExprCos(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_trig.c:1480
SCIP_Bool SCIPisExprSin(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_trig.c:1469
SCIP_RETCODE SCIPcreateExprValue(SCIP *scip, SCIP_EXPR **expr, SCIP_Real value, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_value.c:270
SCIP_RETCODE SCIPincludeExprhdlrCos(SCIP *scip)
Definition: expr_trig.c:1405
SCIP_RETCODE SCIPincludeExprhdlrSin(SCIP *scip)
Definition: expr_trig.c:1380
#define SCIPdebugMsg
Definition: scip_message.h:78
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
const char * SCIPexprhdlrGetName(SCIP_EXPRHDLR *exprhdlr)
Definition: expr.c:545
void SCIPexprhdlrSetHash(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRHASH((*hash)))
Definition: expr.c:451
void SCIPexprhdlrSetCopyFreeHdlr(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCOPYHDLR((*copyhdlr)), SCIP_DECL_EXPRFREEHDLR((*freehdlr)))
Definition: expr.c:370
void SCIPexprhdlrSetDiff(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRBWDIFF((*bwdiff)), SCIP_DECL_EXPRFWDIFF((*fwdiff)), SCIP_DECL_EXPRBWFWDIFF((*bwfwdiff)))
Definition: expr.c:473
void SCIPexprhdlrSetReverseProp(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRREVERSEPROP((*reverseprop)))
Definition: expr.c:510
void SCIPexprhdlrSetParse(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRPARSE((*parse)))
Definition: expr.c:407
void SCIPexprhdlrSetEstimate(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRINITESTIMATES((*initestimates)), SCIP_DECL_EXPRESTIMATE((*estimate)))
Definition: expr.c:532
void SCIPexprhdlrSetMonotonicity(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRMONOTONICITY((*monotonicity)))
Definition: expr.c:429
void SCIPexprhdlrSetIntEval(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRINTEVAL((*inteval)))
Definition: expr.c:488
void SCIPexprhdlrSetCurvature(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRCURVATURE((*curvature)))
Definition: expr.c:418
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
SCIP_EXPRHDLR * SCIPfindExprhdlr(SCIP *scip, const char *name)
Definition: scip_expr.c:868
void SCIPexprhdlrSetSimplify(SCIP_EXPRHDLR *exprhdlr, SCIP_DECL_EXPRSIMPLIFY((*simplify)))
Definition: expr.c:499
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
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3860
SCIP_Bool SCIPisExprValue(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1442
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1417
SCIP_Real SCIPexprGetDot(SCIP_EXPR *expr)
Definition: expr.c:3974
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
SCIP_Real SCIPgetValueExprValue(SCIP_EXPR *expr)
Definition: expr_value.c:294
SCIP_Real SCIPexprGetEvalValue(SCIP_EXPR *expr)
Definition: expr.c:3934
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3870
SCIP_INTERVAL SCIPexprGetActivity(SCIP_EXPR *expr)
Definition: expr.c:4016
void SCIPcaptureExpr(SCIP_EXPR *expr)
Definition: scip_expr.c:1409
SCIP_RETCODE SCIPevalExprActivity(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1717
SCIP_EXPRHDLR * SCIPexprGetHdlr(SCIP_EXPR *expr)
Definition: expr.c:3883
SCIP_Real SCIPintervalGetInf(SCIP_INTERVAL interval)
void SCIPintervalCos(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
SCIP_Bool SCIPintervalIsEmpty(SCIP_Real infinity, SCIP_INTERVAL operand)
void SCIPintervalSin(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand)
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
void SCIPintervalAddScalar(SCIP_Real infinity, SCIP_INTERVAL *resultant, SCIP_INTERVAL operand1, SCIP_Real operand2)
SCIP_Real SCIPintervalGetSup(SCIP_INTERVAL interval)
void SCIPintervalSetEmpty(SCIP_INTERVAL *resultant)
SCIP_Bool SCIPisRelEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisInfinity(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisFeasLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define M_PI
Definition: pricer_rpa.c:97
SCIP_Real sup
Definition: intervalarith.h:56
SCIP_Real inf
Definition: intervalarith.h:55
#define SCIP_DECL_EXPR_OWNERCREATE(x)
Definition: type_expr.h:143
SCIP_EXPRCURV
Definition: type_expr.h:61
@ SCIP_EXPRCURV_CONVEX
Definition: type_expr.h:63
@ SCIP_EXPRCURV_LINEAR
Definition: type_expr.h:65
@ SCIP_EXPRCURV_UNKNOWN
Definition: type_expr.h:62
@ SCIP_EXPRCURV_CONCAVE
Definition: type_expr.h:64
@ SCIP_MONOTONE_UNKNOWN
Definition: type_expr.h:71
@ SCIP_MONOTONE_INC
Definition: type_expr.h:72
@ SCIP_MONOTONE_DEC
Definition: type_expr.h:73
@ SCIP_OKAY
Definition: type_retcode.h:42
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63