1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2008 Simon Ibbotson 5 6 This file is part of QuantLib, a free-software/open-source library 7 for financial quantitative analysts and developers - http://quantlib.org/ 8 9 QuantLib is free software: you can redistribute it and/or modify it 10 under the terms of the QuantLib license. You should have received a 11 copy of the license along with this program; if not, please email 12 <quantlib-dev@lists.sf.net>. The license is also available online at 13 <http://quantlib.org/license.shtml>. 14 15 This program is distributed in the hope that it will be useful, but WITHOUT 16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 17 FOR A PARTICULAR PURPOSE. See the license for more details. 18 */ 19 20 /*! \file convexmonotoneinterpolation.hpp 21 \brief convex monotone interpolation method 22 */ 23 24 #ifndef quantlib_convex_monotone_interpolation_hpp 25 #define quantlib_convex_monotone_interpolation_hpp 26 27 #include <ql/math/interpolation.hpp> 28 #include <map> 29 30 namespace QuantLib { 31 32 namespace detail { 33 template<class I1, class I2> class ConvexMonotoneImpl; 34 class SectionHelper; 35 } 36 37 //! Convex monotone yield-curve interpolation method. 38 /*! Enhances implementation of the convex monotone method 39 described in "Interpolation Methods for Curve Construction" by 40 Hagan & West AMF Vol 13, No2 2006. 41 42 A setting of monotonicity = 1 and quadraticity = 0 will 43 reproduce the basic Hagan/West method. However, this can 44 produce excessive gradients which can mean P&L swings for some 45 curves. Setting monotonicity < 1 and/or quadraticity > 0 46 produces smoother curves. Extra enhancement to avoid negative 47 values (if required) is in place. 48 49 \ingroup interpolations 50 \warning See the Interpolation class for information about the 51 required lifetime of the underlying data. 52 */ 53 template <class I1, class I2> 54 class ConvexMonotoneInterpolation : public Interpolation { 55 typedef std::map<Real, ext::shared_ptr<detail::SectionHelper> > 56 helper_map; 57 public: ConvexMonotoneInterpolation(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin,Real quadraticity,Real monotonicity,bool forcePositive,bool flatFinalPeriod=false,const helper_map & preExistingHelpers=helper_map ())58 ConvexMonotoneInterpolation(const I1& xBegin, const I1& xEnd, 59 const I2& yBegin, Real quadraticity, 60 Real monotonicity, bool forcePositive, 61 bool flatFinalPeriod = false, 62 const helper_map& preExistingHelpers = 63 helper_map()) { 64 impl_ = ext::shared_ptr<Interpolation::Impl>( 65 new detail::ConvexMonotoneImpl<I1,I2>(xBegin, 66 xEnd, 67 yBegin, 68 quadraticity, 69 monotonicity, 70 forcePositive, 71 flatFinalPeriod, 72 preExistingHelpers)); 73 impl_->update(); 74 } 75 ConvexMonotoneInterpolation(Interpolation & interp)76 ConvexMonotoneInterpolation(Interpolation& interp) 77 : Interpolation(interp) {} 78 79 std::map<Real, ext::shared_ptr<detail::SectionHelper> > getExistingHelpers()80 getExistingHelpers() { 81 ext::shared_ptr<detail::ConvexMonotoneImpl<I1,I2> > derived = 82 ext::dynamic_pointer_cast<detail::ConvexMonotoneImpl<I1,I2>, 83 Interpolation::Impl>(impl_); 84 return derived->getExistingHelpers(); 85 } 86 }; 87 88 //! Convex-monotone interpolation factory and traits 89 /*! \ingroup interpolations */ 90 class ConvexMonotone { 91 public: 92 static const bool global = true; 93 static const Size requiredPoints = 2; 94 static const Size dataSizeAdjustment = 1; 95 ConvexMonotone(Real quadraticity=0.3,Real monotonicity=0.7,bool forcePositive=true)96 explicit ConvexMonotone(Real quadraticity = 0.3, 97 Real monotonicity = 0.7, 98 bool forcePositive = true) 99 : quadraticity_(quadraticity), monotonicity_(monotonicity), 100 forcePositive_(forcePositive) {} 101 102 template <class I1, class I2> interpolate(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin) const103 Interpolation interpolate(const I1& xBegin, const I1& xEnd, 104 const I2& yBegin) const { 105 return ConvexMonotoneInterpolation<I1,I2>(xBegin, xEnd, yBegin, 106 quadraticity_, 107 monotonicity_, 108 forcePositive_, 109 false); 110 } 111 112 template <class I1, class I2> localInterpolate(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin,Size localisation,Interpolation & prevInterpolation,Size finalSize) const113 Interpolation localInterpolate(const I1& xBegin, const I1& xEnd, 114 const I2& yBegin, Size localisation, 115 Interpolation& prevInterpolation, 116 Size finalSize) const { 117 Size length = std::distance(xBegin, xEnd); 118 if (length - localisation == 1) { // the first time this 119 // function is called 120 if (length == finalSize) { 121 return ConvexMonotoneInterpolation<I1,I2>(xBegin, xEnd, 122 yBegin, 123 quadraticity_, 124 monotonicity_, 125 forcePositive_, 126 false); 127 } else { 128 return ConvexMonotoneInterpolation<I1,I2>(xBegin, xEnd, 129 yBegin, 130 quadraticity_, 131 monotonicity_, 132 forcePositive_, 133 true); 134 } 135 } 136 137 ConvexMonotoneInterpolation<I1,I2> interp(prevInterpolation); 138 if (length == finalSize) { 139 return ConvexMonotoneInterpolation<I1,I2>( 140 xBegin, xEnd, yBegin, 141 quadraticity_, 142 monotonicity_, 143 forcePositive_, 144 false, 145 interp.getExistingHelpers()); 146 } else { 147 return ConvexMonotoneInterpolation<I1,I2>( 148 xBegin, xEnd, yBegin, 149 quadraticity_, 150 monotonicity_, 151 forcePositive_, 152 true, 153 interp.getExistingHelpers()); 154 } 155 } 156 private: 157 Real quadraticity_, monotonicity_; 158 bool forcePositive_; 159 }; 160 161 162 namespace detail { 163 164 class SectionHelper { 165 public: ~SectionHelper()166 virtual ~SectionHelper() {} 167 virtual Real value(Real x) const = 0; 168 virtual Real primitive(Real x) const = 0; 169 virtual Real fNext() const = 0; 170 }; 171 172 //the first value in the y-vector is ignored. 173 template <class I1, class I2> 174 class ConvexMonotoneImpl : public Interpolation::templateImpl<I1, I2> { 175 typedef std::map<Real, ext::shared_ptr<SectionHelper> > 176 helper_map; 177 public: 178 enum SectionType { 179 EverywhereConstant, 180 ConstantGradient, 181 QuadraticMinimum, 182 QuadraticMaximum 183 }; 184 ConvexMonotoneImpl(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin,Real quadraticity,Real monotonicity,bool forcePositive,bool constantLastPeriod,const helper_map & preExistingHelpers)185 ConvexMonotoneImpl(const I1& xBegin, 186 const I1& xEnd, 187 const I2& yBegin, 188 Real quadraticity, 189 Real monotonicity, 190 bool forcePositive, 191 bool constantLastPeriod, 192 const helper_map& preExistingHelpers) 193 : Interpolation::templateImpl<I1,I2>(xBegin,xEnd,yBegin, 194 ConvexMonotone::requiredPoints), 195 preSectionHelpers_(preExistingHelpers), 196 forcePositive_(forcePositive), 197 constantLastPeriod_(constantLastPeriod), 198 quadraticity_(quadraticity), monotonicity_(monotonicity), 199 length_(xEnd-xBegin) { 200 201 QL_REQUIRE(monotonicity_ >= 0 && monotonicity_ <= 1, 202 "Monotonicity must lie between 0 and 1"); 203 QL_REQUIRE(quadraticity_ >= 0 && quadraticity_ <= 1, 204 "Quadraticity must lie between 0 and 1"); 205 QL_REQUIRE(length_ >= 2, 206 "Single point provided, not supported by convex " 207 "monotone method as first point is ignored"); 208 QL_REQUIRE((length_ - preExistingHelpers.size()) > 1, 209 "Too many existing helpers have been supplied"); 210 } 211 212 void update(); 213 214 Real value(Real x) const; 215 Real primitive(Real x) const; derivative(Real) const216 Real derivative(Real) const { 217 QL_FAIL("Convex-monotone spline derivative not implemented"); 218 } secondDerivative(Real) const219 Real secondDerivative(Real) const { 220 QL_FAIL("Convex-monotone spline second derivative " 221 "not implemented"); 222 } 223 getExistingHelpers()224 helper_map getExistingHelpers() { 225 helper_map retArray(sectionHelpers_); 226 if (constantLastPeriod_) 227 retArray.erase(*(this->xEnd_-1)); 228 return retArray; 229 } 230 private: 231 helper_map sectionHelpers_; 232 helper_map preSectionHelpers_; 233 ext::shared_ptr<SectionHelper> extrapolationHelper_; 234 bool forcePositive_, constantLastPeriod_; 235 Real quadraticity_; 236 Real monotonicity_; 237 Size length_; 238 }; 239 240 241 class ComboHelper : public SectionHelper { 242 public: ComboHelper(ext::shared_ptr<SectionHelper> & quadraticHelper,ext::shared_ptr<SectionHelper> & convMonoHelper,Real quadraticity)243 ComboHelper(ext::shared_ptr<SectionHelper>& quadraticHelper, 244 ext::shared_ptr<SectionHelper>& convMonoHelper, 245 Real quadraticity) 246 : quadraticity_(quadraticity), 247 quadraticHelper_(quadraticHelper), 248 convMonoHelper_(convMonoHelper) { 249 QL_REQUIRE(quadraticity < 1.0 && quadraticity > 0.0, 250 "Quadratic value must lie between 0 and 1"); } 251 value(Real x) const252 Real value(Real x) const { 253 return( quadraticity_*quadraticHelper_->value(x) + (1.0-quadraticity_)*convMonoHelper_->value(x) ); 254 } primitive(Real x) const255 Real primitive(Real x) const { 256 return( quadraticity_*quadraticHelper_->primitive(x) + (1.0-quadraticity_)*convMonoHelper_->primitive(x) ); 257 } fNext() const258 Real fNext() const { 259 return( quadraticity_*quadraticHelper_->fNext() + (1.0-quadraticity_)*convMonoHelper_->fNext() ); 260 } 261 262 private: 263 Real quadraticity_; 264 ext::shared_ptr<SectionHelper> quadraticHelper_; 265 ext::shared_ptr<SectionHelper> convMonoHelper_; 266 }; 267 268 class EverywhereConstantHelper : public SectionHelper { 269 public: EverywhereConstantHelper(Real value,Real prevPrimitive,Real xPrev)270 EverywhereConstantHelper(Real value, Real prevPrimitive, Real xPrev) 271 : value_(value), prevPrimitive_(prevPrimitive), xPrev_(xPrev) 272 {} 273 value(Real) const274 Real value(Real) const { 275 return value_; 276 } primitive(Real x) const277 Real primitive(Real x) const { 278 return prevPrimitive_ + (x-xPrev_)*value_; 279 } fNext() const280 Real fNext() const { 281 return value_; 282 } 283 private: 284 Real value_; 285 Real prevPrimitive_; 286 Real xPrev_; 287 }; 288 289 class ConvexMonotone2Helper : public SectionHelper { 290 public: ConvexMonotone2Helper(Real xPrev,Real xNext,Real gPrev,Real gNext,Real fAverage,Real eta2,Real prevPrimitive)291 ConvexMonotone2Helper(Real xPrev, Real xNext, 292 Real gPrev, Real gNext, 293 Real fAverage, Real eta2, 294 Real prevPrimitive) 295 : xPrev_(xPrev), xScaling_(xNext-xPrev), gPrev_(gPrev), 296 gNext_(gNext), fAverage_(fAverage), eta2_(eta2), 297 prevPrimitive_(prevPrimitive) 298 {} 299 value(Real x) const300 Real value(Real x) const { 301 Real xVal = (x-xPrev_)/xScaling_; 302 if (xVal <= eta2_) { 303 return( fAverage_ + gPrev_ ); 304 } else { 305 return( fAverage_ + gPrev_ + (gNext_-gPrev_)/((1-eta2_)*(1-eta2_))*(xVal-eta2_)*(xVal-eta2_) ); 306 } 307 } 308 primitive(Real x) const309 Real primitive(Real x) const { 310 Real xVal = (x-xPrev_)/xScaling_; 311 if (xVal <= eta2_) { 312 return( prevPrimitive_ + xScaling_*(fAverage_*xVal + gPrev_*xVal) ); 313 } else { 314 return( prevPrimitive_ + xScaling_*(fAverage_*xVal + gPrev_*xVal + (gNext_-gPrev_)/((1-eta2_)*(1-eta2_)) * 315 (1.0/3.0*(xVal*xVal*xVal - eta2_*eta2_*eta2_) - eta2_*xVal*xVal + eta2_*eta2_*xVal) ) ); 316 } 317 } fNext() const318 Real fNext() const { 319 return(fAverage_+gNext_); 320 } 321 private: 322 Real xPrev_, xScaling_, gPrev_, gNext_, fAverage_, eta2_, prevPrimitive_; 323 }; 324 325 class ConvexMonotone3Helper : public SectionHelper { 326 public: ConvexMonotone3Helper(Real xPrev,Real xNext,Real gPrev,Real gNext,Real fAverage,Real eta3,Real prevPrimitive)327 ConvexMonotone3Helper(Real xPrev, Real xNext, 328 Real gPrev, Real gNext, 329 Real fAverage, Real eta3, 330 Real prevPrimitive) 331 : xPrev_(xPrev), xScaling_(xNext-xPrev), gPrev_(gPrev), 332 gNext_(gNext), fAverage_(fAverage), eta3_(eta3), prevPrimitive_(prevPrimitive) 333 {} 334 value(Real x) const335 Real value(Real x) const { 336 Real xVal = (x-xPrev_)/xScaling_; 337 if (xVal <= eta3_) { 338 return( fAverage_ + gNext_ + (gPrev_-gNext_) / (eta3_*eta3_) * (eta3_-xVal)*(eta3_-xVal) ); 339 } else { 340 return( fAverage_ + gNext_ ); 341 } 342 } 343 primitive(Real x) const344 Real primitive(Real x) const { 345 Real xVal = (x-xPrev_)/xScaling_; 346 if (xVal <= eta3_) { 347 return( prevPrimitive_ + xScaling_ * (fAverage_*xVal + gNext_*xVal + (gPrev_-gNext_)/(eta3_*eta3_) * 348 (1.0/3.0 * xVal*xVal*xVal - eta3_*xVal*xVal + eta3_*eta3_*xVal) ) ); 349 } else { 350 return( prevPrimitive_ + xScaling_ * (fAverage_*xVal + gNext_*xVal + (gPrev_-gNext_)/(eta3_*eta3_) * 351 (1.0/3.0 * eta3_*eta3_*eta3_)) ); 352 } 353 } fNext() const354 Real fNext() const { 355 return(fAverage_+gNext_); 356 } 357 private: 358 Real xPrev_, xScaling_, gPrev_, gNext_, fAverage_, eta3_, prevPrimitive_; 359 }; 360 361 class ConvexMonotone4Helper : public SectionHelper { 362 public: ConvexMonotone4Helper(Real xPrev,Real xNext,Real gPrev,Real gNext,Real fAverage,Real eta4,Real prevPrimitive)363 ConvexMonotone4Helper(Real xPrev, Real xNext, 364 Real gPrev, Real gNext, 365 Real fAverage, Real eta4, 366 Real prevPrimitive) 367 : xPrev_(xPrev), xScaling_(xNext-xPrev), gPrev_(gPrev), 368 gNext_(gNext), fAverage_(fAverage), eta4_(eta4), prevPrimitive_(prevPrimitive) { 369 A_ = -0.5*(eta4_*gPrev_ + (1-eta4_)*gNext_); 370 } 371 value(Real x) const372 Real value(Real x) const { 373 Real xVal = (x-xPrev_)/xScaling_; 374 if (xVal <= eta4_) { 375 return(fAverage_ + A_ + (gPrev_-A_)*(eta4_-xVal)*(eta4_-xVal)/(eta4_*eta4_) ); 376 } else { 377 return(fAverage_ + A_ + (gNext_-A_)*(xVal-eta4_)*(xVal-eta4_)/((1-eta4_)*(1-eta4_)) ); 378 } 379 } 380 primitive(Real x) const381 Real primitive(Real x) const { 382 Real xVal = (x-xPrev_)/xScaling_; 383 Real retVal; 384 if (xVal <= eta4_) { 385 retVal = prevPrimitive_ + xScaling_ * (fAverage_ + A_ + (gPrev_-A_)/(eta4_*eta4_) * 386 (eta4_*eta4_ - eta4_*xVal + 1.0/3.0*xVal*xVal)) * xVal; 387 } else { 388 retVal = prevPrimitive_ + xScaling_ *(fAverage_*xVal + A_*xVal + (gPrev_-A_)*(1.0/3.0*eta4_) + 389 (gNext_-A_)/((1-eta4_)*(1-eta4_)) * 390 (1.0/3.0*xVal*xVal*xVal - eta4_*xVal*xVal + eta4_*eta4_*xVal - 1.0/3.0*eta4_*eta4_*eta4_)); 391 } 392 return retVal; 393 } fNext() const394 Real fNext() const { 395 return(fAverage_+gNext_); 396 } 397 protected: 398 Real xPrev_, xScaling_, gPrev_, gNext_, fAverage_, eta4_, prevPrimitive_; 399 Real A_; 400 }; 401 402 class ConvexMonotone4MinHelper : public ConvexMonotone4Helper { 403 public: ConvexMonotone4MinHelper(Real xPrev,Real xNext,Real gPrev,Real gNext,Real fAverage,Real eta4,Real prevPrimitive)404 ConvexMonotone4MinHelper(Real xPrev, Real xNext, 405 Real gPrev, Real gNext, 406 Real fAverage, Real eta4, 407 Real prevPrimitive) 408 : ConvexMonotone4Helper(xPrev, xNext, gPrev, gNext, 409 fAverage, eta4, prevPrimitive), 410 splitRegion_(false) { 411 if ( A_+ fAverage_ <= 0.0 ) { 412 splitRegion_ = true; 413 Real fPrev = gPrev_+fAverage_; 414 Real fNext = gNext_+fAverage_; 415 Real reqdShift = (eta4_*fPrev + (1-eta4_)*fNext)/3.0 - fAverage_; 416 Real reqdPeriod = reqdShift * xScaling_ / (fAverage_+reqdShift); 417 Real xAdjust = xScaling_ - reqdPeriod; 418 xRatio_ = xAdjust/xScaling_; 419 420 fAverage_ += reqdShift; 421 gNext_ = fNext - fAverage_; 422 gPrev_ = fPrev - fAverage_; 423 A_ = -(eta4_ * gPrev_ + (1.0-eta4)*gNext_)/2.0; 424 x2_ = xPrev_ + xAdjust * eta4_; 425 x3_ = xPrev_ + xScaling_ - xAdjust*(1.0-eta4_); 426 } 427 } 428 value(Real x) const429 Real value(Real x) const { 430 if (!splitRegion_) 431 return ConvexMonotone4Helper::value(x); 432 433 Real xVal = (x-xPrev_)/xScaling_; 434 if (x <= x2_) { 435 xVal /= xRatio_; 436 return(fAverage_ + A_ + (gPrev_-A_)*(eta4_-xVal)*(eta4_-xVal)/(eta4_*eta4_)); 437 } else if (x < x3_) { 438 return 0.0; 439 } else { 440 xVal = 1.0 - (1.0 - xVal) / xRatio_; 441 return(fAverage_ + A_ + (gNext_-A_)*(xVal-eta4_)*(xVal-eta4_)/((1-eta4_)*(1-eta4_)) ); 442 } 443 } 444 primitive(Real x) const445 Real primitive(Real x) const { 446 if (!splitRegion_) 447 return ConvexMonotone4Helper::primitive(x); 448 449 Real xVal = (x-xPrev_)/xScaling_; 450 if (x <= x2_) { 451 xVal /= xRatio_; 452 return( prevPrimitive_ + xScaling_*xRatio_*(fAverage_ + A_ + (gPrev_-A_)/(eta4_*eta4_) * 453 (eta4_*eta4_ - eta4_*xVal + 1.0/3.0*xVal*xVal)) * xVal ); 454 } else if (x <= x3_) { 455 return( prevPrimitive_ + xScaling_*xRatio_*(fAverage_*eta4_ + A_*eta4_ + (gPrev_-A_)/(eta4_*eta4_) * 456 (1.0/3.0*eta4_*eta4_*eta4_)) ); 457 } else { 458 xVal = 1.0 - (1.0-xVal)/xRatio_; 459 return( prevPrimitive_ + xScaling_*xRatio_*(fAverage_*xVal + A_*xVal + (gPrev_-A_)*(1.0/3.0*eta4_) + 460 (gNext_-A_) / ((1.0-eta4_)*(1.0-eta4_)) * 461 (1.0/3.0*xVal*xVal*xVal - eta4_*xVal*xVal + eta4_*eta4_*xVal - 1.0/3.0*eta4_*eta4_*eta4_)) ); 462 } 463 } 464 private: 465 bool splitRegion_; 466 Real xRatio_, x2_, x3_; 467 }; 468 469 class ConstantGradHelper : public SectionHelper { 470 public: ConstantGradHelper(Real fPrev,Real prevPrimitive,Real xPrev,Real xNext,Real fNext)471 ConstantGradHelper(Real fPrev, Real prevPrimitive, 472 Real xPrev, Real xNext, Real fNext) 473 : fPrev_(fPrev), prevPrimitive_(prevPrimitive), 474 xPrev_(xPrev), fGrad_((fNext-fPrev)/(xNext-xPrev)),fNext_(fNext) 475 {} 476 value(Real x) const477 Real value(Real x) const { 478 return (fPrev_ + (x-xPrev_)*fGrad_); 479 } primitive(Real x) const480 Real primitive(Real x) const { 481 return (prevPrimitive_+(x-xPrev_)*(fPrev_+0.5*(x-xPrev_)*fGrad_)); 482 } fNext() const483 Real fNext() const { 484 return fNext_; 485 } 486 private: 487 Real fPrev_, prevPrimitive_, xPrev_, fGrad_, fNext_; 488 }; 489 490 class QuadraticHelper : public SectionHelper { 491 public: QuadraticHelper(Real xPrev,Real xNext,Real fPrev,Real fNext,Real fAverage,Real prevPrimitive)492 QuadraticHelper(Real xPrev, Real xNext, 493 Real fPrev, Real fNext, 494 Real fAverage, 495 Real prevPrimitive) 496 : xPrev_(xPrev), xNext_(xNext), fPrev_(fPrev), 497 fNext_(fNext), fAverage_(fAverage), 498 prevPrimitive_(prevPrimitive) { 499 a_ = 3*fPrev_ + 3*fNext_ - 6*fAverage_; 500 b_ = -(4*fPrev_ + 2*fNext_ - 6*fAverage_); 501 c_ = fPrev_; 502 xScaling_ = xNext_-xPrev_; 503 } 504 value(Real x) const505 Real value(Real x) const { 506 Real xVal = (x-xPrev_)/xScaling_; 507 return( a_*xVal*xVal + b_*xVal + c_ ); 508 } 509 primitive(Real x) const510 Real primitive(Real x) const { 511 Real xVal = (x-xPrev_)/xScaling_; 512 return( prevPrimitive_ + xScaling_ * (a_/3*xVal*xVal + b_/2*xVal + c_) * xVal ); 513 } 514 fNext() const515 Real fNext() const { 516 return fNext_; 517 } 518 private: 519 Real xPrev_, xNext_, fPrev_, fNext_, fAverage_, prevPrimitive_; 520 Real xScaling_, a_, b_, c_; 521 }; 522 523 class QuadraticMinHelper : public SectionHelper { 524 public: QuadraticMinHelper(Real xPrev,Real xNext,Real fPrev,Real fNext,Real fAverage,Real prevPrimitive)525 QuadraticMinHelper(Real xPrev, Real xNext, 526 Real fPrev, Real fNext, 527 Real fAverage, 528 Real prevPrimitive) 529 : splitRegion_(false), x1_(xPrev), x4_(xNext), 530 primitive1_(prevPrimitive), fAverage_(fAverage), 531 fPrev_(fPrev), fNext_(fNext), xRatio_(1.0) { 532 a_ = 3*fPrev_ + 3*fNext_ - 6*fAverage_; 533 b_ = -(4*fPrev_ + 2*fNext_ - 6*fAverage_); 534 c_ = fPrev_; 535 Real d = b_*b_-4*a_*c_; 536 xScaling_ = x4_-x1_; 537 if (d > 0) { 538 Real aAv = 36; 539 Real bAv = -24*(fPrev_+fNext_); 540 Real cAv = 4*(fPrev_*fPrev_ + fPrev_*fNext_ + fNext_*fNext_); 541 Real dAv = bAv*bAv - 4.0*aAv*cAv; 542 if (dAv >= 0.0) { 543 splitRegion_ = true; 544 Real avRoot = (-bAv - std::sqrt(dAv))/(2*aAv); 545 546 xRatio_ = fAverage_ / avRoot; 547 xScaling_ *= xRatio_; 548 549 a_ = 3*fPrev_ + 3*fNext_ - 6*avRoot; 550 b_ = -(4*fPrev_ + 2*fNext_ - 6*avRoot); 551 c_ = fPrev_; 552 Real xRoot = -b_/(2*a_); 553 x2_ = x1_ + xRatio_ * (x4_-x1_) * xRoot; 554 x3_ = x4_ - xRatio_ * (x4_-x1_) * (1-xRoot); 555 primitive2_ = 556 primitive1_ + xScaling_*(a_/3*xRoot*xRoot + b_/2*xRoot + c_)*xRoot; 557 } 558 } 559 } 560 value(Real x) const561 Real value(Real x) const { 562 Real xVal = (x - x1_) / (x4_-x1_); 563 if (splitRegion_) { 564 if (x <= x2_) { 565 xVal /= xRatio_; 566 } else if (x < x3_) { 567 return 0.0; 568 } else { 569 xVal = 1.0 - (1.0 - xVal) / xRatio_; 570 } 571 } 572 573 return c_ + b_*xVal + a_*xVal*xVal; 574 } 575 primitive(Real x) const576 Real primitive(Real x) const { 577 Real xVal = (x - x1_) / (x4_-x1_); 578 if (splitRegion_) { 579 if (x < x2_) { 580 xVal /= xRatio_; 581 } else if (x < x3_) { 582 return primitive2_; 583 } else { 584 xVal = 1.0 - (1.0 - xVal) / xRatio_; 585 } 586 } 587 return primitive1_ + xScaling_ * (a_/3*xVal*xVal+ b_/2*xVal+c_)*xVal; 588 } 589 fNext() const590 Real fNext() const { 591 return fNext_; 592 } 593 594 private: 595 bool splitRegion_; 596 Real x1_, x2_, x3_, x4_; 597 Real a_, b_, c_; 598 Real primitive1_, primitive2_; 599 Real fAverage_, fPrev_, fNext_, xScaling_, xRatio_; 600 }; 601 602 template <class I1, class I2> update()603 void ConvexMonotoneImpl<I1,I2>::update() { 604 sectionHelpers_.clear(); 605 if (length_ == 2) { //single period 606 ext::shared_ptr<SectionHelper> singleHelper( 607 new EverywhereConstantHelper(this->yBegin_[1], 608 0.0, 609 this->xBegin_[0])); 610 sectionHelpers_[this->xBegin_[1]] = singleHelper; 611 extrapolationHelper_ = singleHelper; 612 return; 613 } 614 615 std::vector<Real> f(length_); 616 sectionHelpers_ = preSectionHelpers_; 617 Size startPoint = sectionHelpers_.size()+1; 618 619 //first derive the boundary forwards. 620 for (Size i=startPoint; i<length_-1; ++i) { 621 Real dxPrev = this->xBegin_[i] - this->xBegin_[i-1]; 622 Real dx = this->xBegin_[i+1] - this->xBegin_[i]; 623 f[i] = dx/(dx+dxPrev) * this->yBegin_[i] 624 + dxPrev/(dx+dxPrev) * this->yBegin_[i+1]; 625 } 626 627 if (startPoint > 1) 628 f[startPoint-1] = preSectionHelpers_.rbegin()->second->fNext(); 629 if (startPoint == 1) 630 f[0] = 1.5 * this->yBegin_[1] - 0.5 * f[1]; 631 632 f[length_-1] = 1.5 * this->yBegin_[length_-1] - 0.5 * f[length_-2]; 633 634 if (forcePositive_) { 635 if (f[0] < 0) 636 f[0] = 0.0; 637 if (f[length_-1] < 0.0) 638 f[length_-1] = 0.0; 639 } 640 641 Real primitive = 0.0; 642 for (Size i = 0; i < startPoint-1; ++i) 643 primitive += 644 this->yBegin_[i+1] * (this->xBegin_[i+1]-this->xBegin_[i]); 645 646 Size endPoint = length_; 647 //constantLastPeriod_ = false; 648 if (constantLastPeriod_) 649 endPoint = endPoint-1; 650 651 for (Size i=startPoint; i< endPoint; ++i) { 652 Real gPrev = f[i-1] - this->yBegin_[i]; 653 Real gNext = f[i] - this->yBegin_[i]; 654 //first deal with the zero gradient case 655 if ( std::fabs(gPrev) < 1.0E-14 656 && std::fabs(gNext) < 1.0E-14 ) { 657 ext::shared_ptr<SectionHelper> singleHelper( 658 new ConstantGradHelper(f[i-1], primitive, 659 this->xBegin_[i-1], 660 this->xBegin_[i], 661 f[i])); 662 sectionHelpers_[this->xBegin_[i]] = singleHelper; 663 } else { 664 Real quadraticity = quadraticity_; 665 ext::shared_ptr<SectionHelper> quadraticHelper; 666 ext::shared_ptr<SectionHelper> convMonotoneHelper; 667 if (quadraticity_ > 0.0) { 668 if (gPrev >= -2.0*gNext && gPrev > -0.5*gNext && forcePositive_) { 669 quadraticHelper = 670 ext::shared_ptr<SectionHelper>( 671 new QuadraticMinHelper(this->xBegin_[i-1], 672 this->xBegin_[i], 673 f[i-1], f[i], 674 this->yBegin_[i], 675 primitive) ); 676 } else { 677 quadraticHelper = 678 ext::shared_ptr<SectionHelper>( 679 new QuadraticHelper(this->xBegin_[i-1], 680 this->xBegin_[i], 681 f[i-1], f[i], 682 this->yBegin_[i], 683 primitive) ); 684 } 685 } 686 if (quadraticity_ < 1.0) { 687 688 if ((gPrev > 0.0 && -0.5*gPrev >= gNext && gNext >= -2.0*gPrev) || 689 (gPrev < 0.0 && -0.5*gPrev <= gNext && gNext <= -2.0*gPrev)) { 690 quadraticity = 1.0; 691 if (quadraticity_ == 0) { 692 if (forcePositive_) { 693 quadraticHelper = 694 ext::shared_ptr<SectionHelper>( 695 new QuadraticMinHelper( 696 this->xBegin_[i-1], 697 this->xBegin_[i], 698 f[i-1], f[i], 699 this->yBegin_[i], 700 primitive) ); 701 } else { 702 quadraticHelper = 703 ext::shared_ptr<SectionHelper>( 704 new QuadraticHelper( 705 this->xBegin_[i-1], 706 this->xBegin_[i], 707 f[i-1], f[i], 708 this->yBegin_[i], 709 primitive) ); 710 } 711 } 712 } 713 else if ( (gPrev < 0.0 && gNext > -2.0*gPrev) || 714 (gPrev > 0.0 && gNext < -2.0*gPrev)) { 715 716 Real eta = (gNext + 2.0*gPrev)/(gNext - gPrev); 717 Real b2 = (1.0 + monotonicity_)/2.0; 718 if (eta < b2) { 719 convMonotoneHelper = 720 ext::shared_ptr<SectionHelper>( 721 new ConvexMonotone2Helper( 722 this->xBegin_[i-1], 723 this->xBegin_[i], 724 gPrev, gNext, 725 this->yBegin_[i], 726 eta, primitive)); 727 } else { 728 if (forcePositive_) { 729 convMonotoneHelper = 730 ext::shared_ptr<SectionHelper>( 731 new ConvexMonotone4MinHelper( 732 this->xBegin_[i-1], 733 this->xBegin_[i], 734 gPrev, gNext, 735 this->yBegin_[i], 736 b2, primitive)); 737 } else { 738 convMonotoneHelper = 739 ext::shared_ptr<SectionHelper>( 740 new ConvexMonotone4Helper( 741 this->xBegin_[i-1], 742 this->xBegin_[i], 743 gPrev, gNext, 744 this->yBegin_[i], 745 b2, primitive)); 746 } 747 } 748 } 749 else if ( (gPrev > 0.0 && gNext < 0.0 && gNext > -0.5*gPrev) || 750 (gPrev < 0.0 && gNext > 0.0 && gNext < -0.5*gPrev) ) { 751 Real eta = gNext/(gNext-gPrev) * 3.0; 752 Real b3 = (1.0 - monotonicity_)/2.0; 753 if (eta > b3) { 754 convMonotoneHelper = 755 ext::shared_ptr<SectionHelper>( 756 new ConvexMonotone3Helper( 757 this->xBegin_[i-1], 758 this->xBegin_[i], 759 gPrev, gNext, 760 this->yBegin_[i], 761 eta, primitive)); 762 } else { 763 if (forcePositive_) { 764 convMonotoneHelper = 765 ext::shared_ptr<SectionHelper>( 766 new ConvexMonotone4MinHelper( 767 this->xBegin_[i-1], 768 this->xBegin_[i], 769 gPrev, gNext, 770 this->yBegin_[i], 771 b3, primitive)); 772 } else { 773 convMonotoneHelper = 774 ext::shared_ptr<SectionHelper>( 775 new ConvexMonotone4Helper( 776 this->xBegin_[i-1], 777 this->xBegin_[i], 778 gPrev, gNext, 779 this->yBegin_[i], 780 b3, primitive)); 781 } 782 } 783 } else { 784 Real eta = gNext/(gPrev + gNext); 785 Real b2 = (1.0 + monotonicity_)/2.0; 786 Real b3 = (1.0 - monotonicity_)/2.0; 787 if (eta > b2) 788 eta = b2; 789 if (eta < b3) 790 eta = b3; 791 if (forcePositive_) { 792 convMonotoneHelper = 793 ext::shared_ptr<SectionHelper>( 794 new ConvexMonotone4MinHelper( 795 this->xBegin_[i-1], 796 this->xBegin_[i], 797 gPrev, gNext, 798 this->yBegin_[i], 799 eta, primitive)); 800 } else { 801 convMonotoneHelper = 802 ext::shared_ptr<SectionHelper>( 803 new ConvexMonotone4Helper( 804 this->xBegin_[i-1], 805 this->xBegin_[i], 806 gPrev, gNext, 807 this->yBegin_[i], 808 eta, primitive)); 809 } 810 } 811 } 812 813 if (quadraticity == 1.0) { 814 sectionHelpers_[this->xBegin_[i]] = quadraticHelper; 815 } else if (quadraticity == 0.0) { 816 sectionHelpers_[this->xBegin_[i]] = convMonotoneHelper; 817 } else { 818 sectionHelpers_[this->xBegin_[i]] = 819 ext::shared_ptr<SectionHelper>( 820 new ComboHelper(quadraticHelper, 821 convMonotoneHelper, 822 quadraticity)); 823 } 824 825 } 826 primitive += 827 this->yBegin_[i] * (this->xBegin_[i]-this->xBegin_[i-1]); 828 } 829 830 if (constantLastPeriod_) { 831 sectionHelpers_[this->xBegin_[length_-1]] = 832 ext::shared_ptr<SectionHelper>( 833 new EverywhereConstantHelper(this->yBegin_[length_-1], 834 primitive, 835 this->xBegin_[length_-2])); 836 extrapolationHelper_ = sectionHelpers_[this->xBegin_[length_-1]]; 837 } else { 838 extrapolationHelper_ = 839 ext::shared_ptr<SectionHelper>( 840 new EverywhereConstantHelper( 841 (sectionHelpers_.rbegin())->second->value(*(this->xEnd_-1)), 842 primitive, 843 *(this->xEnd_-1))); 844 } 845 } 846 847 template <class I1, class I2> value(Real x) const848 Real ConvexMonotoneImpl<I1,I2>::value(Real x) const { 849 if (x >= *(this->xEnd_-1)) { 850 return extrapolationHelper_->value(x); 851 } 852 853 return sectionHelpers_.upper_bound(x)->second->value(x); 854 } 855 856 template <class I1, class I2> primitive(Real x) const857 Real ConvexMonotoneImpl<I1,I2>::primitive(Real x) const { 858 if (x >= *(this->xEnd_-1)) { 859 return extrapolationHelper_->primitive(x); 860 } 861 862 return sectionHelpers_.upper_bound(x)->second->primitive(x); 863 } 864 865 } 866 867 } 868 869 #endif 870