QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
lagrangeinterpolation.hpp
Go to the documentation of this file.
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>
25#if defined(QL_EXTRA_SAFETY_CHECKS)
26#include <set>
27#endif
28
29namespace 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 {
38 public:
39 virtual ~UpdatedYInterpolation() = default;
40 virtual Real value(const Array& yValues, Real x) const = 0;
41 };
42
43 template <class I1, class I2>
45 : public Interpolation::templateImpl<I1,I2>,
47
48 public:
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
60 void update() override {
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
75 Real value(Real x) const override { return _value(this->yBegin_, x); }
76
77 Real derivative(Real x) const override {
78 Real n=0.0, d=0.0, nd=0.0, dd=0.0;
79 for (Size i=0; i < n_; ++i) {
80 const Real x_i = this->xBegin_[i];
81
82 if (close_enough(x, x_i)) {
83 Real p=0.0;
84 for (Size j=0; j < n_; ++j)
85 if (i != j) {
86 p+=lambda_[j]/(x-this->xBegin_[j])
87 *(this->yBegin_[j] - this->yBegin_[i]);
88 }
89 return p/lambda_[i];
90 }
91
92 const Real alpha = lambda_[i]/(x-x_i);
93 const Real alphad = -alpha/(x-x_i);
94 n += alpha * this->yBegin_[i];
95 d += alpha;
96 nd += alphad * this->yBegin_[i];
97 dd += alphad;
98 }
99 return (nd * d - n * dd)/(d*d);
100 }
101
102 Real primitive(Real) const override {
103 QL_FAIL("LagrangeInterpolation primitive is not implemented");
104 }
105
106 Real secondDerivative(Real) const override {
107 QL_FAIL("LagrangeInterpolation secondDerivative "
108 "is not implemented");
109 }
110
111 Real value(const Array& y, Real x) const override { return _value(y.begin(), x); }
112
113 private:
114 template <class Iterator>
115 inline Real _value(const Iterator& yBegin, Real x) const {
116
117 const Real eps = 10*QL_EPSILON*std::abs(x);
118 const auto iter = std::lower_bound(
119 this->xBegin_, this->xEnd_, x - eps);
120 if (iter != this->xEnd_ && *iter - x < eps) {
121 return yBegin[std::distance(this->xBegin_, iter)];
122 }
123
124 Real n = 0.0, d = 0.0;
125 for (Size i = 0; i < n_; ++i) {
126 const Real alpha = lambda_[i] / (x - this->xBegin_[i]);
127 n += alpha * yBegin[i];
128 d += alpha;
129 }
130 return n / d;
131 }
132
133 const Size n_;
135 };
136 }
137
138 /*! \ingroup interpolations
139 \warning See the Interpolation class for information about the
140 required lifetime of the underlying data.
141 */
143 public:
144 template <class I1, class I2>
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
153 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
1-D array used in linear algebra.
1-D array used in linear algebra.
Definition: array.hpp:52
basic template implementation
templateImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const int requiredPoints=2)
base class for 1-D interpolations.
ext::shared_ptr< Impl > impl_
Real value(const Array &y, Real x) const
LagrangeInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Real _value(const Iterator &yBegin, Real x) const
Real value(const Array &y, Real x) const override
LagrangeInterpolationImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
virtual Real value(const Array &yValues, Real x) const =0
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
base class for 1-D interpolations
Definition: any.hpp:35
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
STL namespace.
Real alpha
Definition: sabr.cpp:200