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