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