1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2016 Klaus Spanderen
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 #ifndef quantlib_lagrange_interpolation_hpp
21 #define quantlib_lagrange_interpolation_hpp
22 
23 #include <ql/math/array.hpp>
24 #include <ql/math/interpolation.hpp>
25 #if defined(QL_EXTRA_SAFETY_CHECKS)
26 #include <set>
27 #endif
28 
29 namespace QuantLib {
30     /*! References: J-P. Berrut and L.N. Trefethen,
31                     Barycentric Lagrange interpolation,
32                     SIAM Review, 46(3):501–517, 2004.
33         https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
34     */
35 
36     namespace detail {
37         class UpdatedYInterpolation {
38           public:
~UpdatedYInterpolation()39             virtual ~UpdatedYInterpolation() {}
40             virtual Real value(const Array& yValues, Real x) const = 0;
41         };
42 
43         template <class I1, class I2>
44         class LagrangeInterpolationImpl
45             : public Interpolation::templateImpl<I1,I2>,
46               public UpdatedYInterpolation {
47 
48           public:
LagrangeInterpolationImpl(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin)49             LagrangeInterpolationImpl(const I1& xBegin, const I1& xEnd,
50                                       const I2& yBegin)
51             : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
52               n_(std::distance(xBegin, xEnd)),
53               lambda_(n_) {
54                 #if defined(QL_EXTRA_SAFETY_CHECKS)
55                 QL_REQUIRE(std::set<Real>(xBegin, xEnd).size() == n_,
56                         "x values must not contain duplicates");
57                 #endif
58             }
59 
update()60             void update() {
61                 const Real cM1 = 4.0/(*(this->xEnd_-1) - *(this->xBegin_));
62 
63                 for (Size i=0; i < n_; ++i) {
64                     lambda_[i] = 1.0;
65 
66                     const Real x_i = this->xBegin_[i];
67                     for (Size j=0; j < n_; ++j) {
68                         if (i != j)
69                             lambda_[i] *= cM1*(x_i-this->xBegin_[j]);
70                     }
71                     lambda_[i] = 1.0/lambda_[i];
72                 }
73             }
74 
value(Real x) const75             Real value(Real x) const {
76                 return _value(this->yBegin_, x);
77             }
78 
derivative(Real x) const79             Real derivative(Real x) const {
80                 Real n=0.0, d=0.0, nd=0.0, dd=0.0;
81                 for (Size i=0; i < n_; ++i) {
82                     const Real x_i = this->xBegin_[i];
83 
84                     if (close_enough(x, x_i)) {
85                         Real p=0.0;
86                         for (Size j=0; j < n_; ++j)
87                             if (i != j) {
88                                 p+=lambda_[j]/(x-this->xBegin_[j])
89                                     *(this->yBegin_[j] - this->yBegin_[i]);
90                             }
91                         return p/lambda_[i];
92                     }
93 
94                     const Real alpha = lambda_[i]/(x-x_i);
95                     const Real alphad = -alpha/(x-x_i);
96                     n += alpha * this->yBegin_[i];
97                     d += alpha;
98                     nd += alphad * this->yBegin_[i];
99                     dd += alphad;
100                 }
101                 return (nd * d - n * dd)/(d*d);
102             }
103 
primitive(Real) const104             Real primitive(Real) const {
105                 QL_FAIL("LagrangeInterpolation primitive is not implemented");
106             }
107 
secondDerivative(Real) const108             Real secondDerivative(Real) const {
109                 QL_FAIL("LagrangeInterpolation secondDerivative "
110                         "is not implemented");
111             }
112 
value(const Array & y,Real x) const113             Real value(const Array& y, Real x) const {
114                 return _value(y.begin(), x);
115             }
116 
117             private:
118               template <class Iterator>
_value(const Iterator & yBegin,Real x) const119               Real _value(const Iterator& yBegin, Real x) const {
120                   Real n=0.0, d=0.0;
121                   for (Size i=0; i < n_; ++i) {
122                       const Real x_i = this->xBegin_[i];
123                       if (close_enough(x, x_i))
124                           return yBegin[i];
125 
126                       const Real alpha = lambda_[i]/(x-x_i);
127                       n += alpha * yBegin[i];
128                       d += alpha;
129                   }
130                   return n/d;
131               }
132 
133               const Size n_;
134               Array lambda_;
135         };
136     }
137 
138     /*! \ingroup interpolations
139         \warning See the Interpolation class for information about the
140                  required lifetime of the underlying data.
141     */
142     class LagrangeInterpolation : public Interpolation {
143       public:
144         template <class I1, class I2>
LagrangeInterpolation(const I1 & xBegin,const I1 & xEnd,const I2 & yBegin)145         LagrangeInterpolation(const I1& xBegin, const I1& xEnd,
146                               const I2& yBegin) {
147             impl_ = ext::make_shared<detail::LagrangeInterpolationImpl<I1,I2> >(
148                 xBegin, xEnd, yBegin);
149             impl_->update();
150         }
151 
152         // interpolate with new set of y values for a new x value
value(const Array & y,Real x) const153         Real value(const Array& y, Real x) const {
154             return ext::dynamic_pointer_cast<detail::UpdatedYInterpolation>
155                 (impl_)->value(y, x);
156         }
157     };
158 
159 }
160 
161 #endif
162