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