/src/quantlib/ql/math/interpolations/lagrangeinterpolation.hpp

Source (jump to first uncovered line)
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2016 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#ifndef quantlib_lagrange_interpolation_hpp
#define quantlib_lagrange_interpolation_hpp

#include <ql/math/array.hpp>
#include <ql/math/interpolation.hpp>
#if defined(QL_EXTRA_SAFETY_CHECKS)
#include <set>
#endif

namespace QuantLib {
    /*! References: J-P. Berrut and L.N. Trefethen,
                    Barycentric Lagrange interpolation,
                    SIAM Review, 46(3):501–517, 2004.
        https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
    */

    namespace detail {
        class UpdatedYInterpolation {
          public:
            virtual ~UpdatedYInterpolation() = default;
            virtual Real value(const Array& yValues, Real x) const = 0;
        };

        template <class I1, class I2>
        class LagrangeInterpolationImpl
            : public Interpolation::templateImpl<I1,I2>,
              public UpdatedYInterpolation {

          public:
            LagrangeInterpolationImpl(const I1& xBegin, const I1& xEnd,
                                      const I2& yBegin)
            : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
              n_(std::distance(xBegin, xEnd)),
              lambda_(n_) {
                #if defined(QL_EXTRA_SAFETY_CHECKS)
                QL_REQUIRE(std::set<Real>(xBegin, xEnd).size() == n_,
                        "x values must not contain duplicates");
                #endif
            }

            void update() override {
                const Real cM1 = 4.0/(*(this->xEnd_-1) - *(this->xBegin_));

                for (Size i=0; i < n_; ++i) {
                    lambda_[i] = 1.0;

                    const Real x_i = this->xBegin_[i];
                    for (Size j=0; j < n_; ++j) {
                        if (i != j)
                            lambda_[i] *= cM1*(x_i-this->xBegin_[j]);
                    }
                    lambda_[i] = 1.0/lambda_[i];
                }
            }

            Real value(Real x) const override { return _value(this->yBegin_, x); }

            Real derivative(Real x) const override {
                Real n=0.0, d=0.0, nd=0.0, dd=0.0;
                for (Size i=0; i < n_; ++i) {
                    const Real x_i = this->xBegin_[i];

                    if (close_enough(x, x_i)) {
                        Real p=0.0;
                        for (Size j=0; j < n_; ++j)
                            if (i != j) {
                                p+=lambda_[j]/(x-this->xBegin_[j])
                                    *(this->yBegin_[j] - this->yBegin_[i]);
                            }
                        return p/lambda_[i];
                    }

                    const Real alpha = lambda_[i]/(x-x_i);
                    const Real alphad = -alpha/(x-x_i);
                    n += alpha * this->yBegin_[i];
                    d += alpha;
                    nd += alphad * this->yBegin_[i];
                    dd += alphad;
                }
                return (nd * d - n * dd)/(d*d);
            }

            Real primitive(Real) const override {
                QL_FAIL("LagrangeInterpolation primitive is not implemented");
            }

            Real secondDerivative(Real) const override {
                QL_FAIL("LagrangeInterpolation secondDerivative "
                        "is not implemented");
            }

            Real value(const Array& y, Real x) const override { return _value(y.begin(), x); }

          private:
            template <class Iterator>
            Real _value(const Iterator& yBegin, Real x) const {

                const Real eps = 10*QL_EPSILON*std::abs(x);
                const auto iter = std::lower_bound(
                    this->xBegin_, this->xEnd_, x - eps);
                if (iter != this->xEnd_ && *iter - x < eps) {
                    return yBegin[std::distance(this->xBegin_, iter)];
                }

                Real n = 0.0, d = 0.0;
                for (Size i = 0; i < n_; ++i) {
                    const Real alpha = lambda_[i] / (x - this->xBegin_[i]);
                    n += alpha * yBegin[i];
                    d += alpha;
                }
                return n / d;
              }

              const Size n_;
              Array lambda_;
        };
    }

    /*! \ingroup interpolations
        \warning See the Interpolation class for information about the
                 required lifetime of the underlying data.
    */
    class LagrangeInterpolation : public Interpolation {
      public:
        template <class I1, class I2>
        LagrangeInterpolation(const I1& xBegin, const I1& xEnd,
                              const I2& yBegin) {
            impl_ = ext::make_shared<detail::LagrangeInterpolationImpl<I1,I2> >(
                xBegin, xEnd, yBegin);
            impl_->update();
        }

        // interpolate with new set of y values for a new x value
        Real value(const Array& y, Real x) const {
            return ext::dynamic_pointer_cast<detail::UpdatedYInterpolation>
                (impl_)->value(y, x);
        }
    };

}

#endif

Line	Count	Source (jump to first uncovered line)
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:
39	0	virtual ~UpdatedYInterpolation() = default;
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:
49		LagrangeInterpolationImpl(const I1& xBegin, const I1& xEnd,
50		const I2& yBegin)
51	0	: Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
52	0	n_(std::distance(xBegin, xEnd)),
53	0	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	0	} Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::LagrangeInterpolationImpl(double const* const&, double const* const&, double const* const&) Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::LagrangeInterpolationImpl(double* const&, double* const&, double* const&)
59
60	0	void update() override {
61	0	const Real cM1 = 4.0/((this->xEnd_-1) - (this->xBegin_));
62
63	0	for (Size i=0; i < n_; ++i) {
64	0	lambda_[i] = 1.0;
65
66	0	const Real x_i = this->xBegin_[i];
67	0	for (Size j=0; j < n_; ++j) {
68	0	if (i != j)
69	0	lambda_[i] = cM1(x_i-this->xBegin_[j]);
70	0	}
71	0	lambda_[i] = 1.0/lambda_[i];
72	0	}
73	0	} Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::update() Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::update()
74
75	0	Real value(Real x) const override { return _value(this->yBegin_, x); } Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::value(double) const Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::value(double) const
76
77	0	Real derivative(Real x) const override {
78	0	Real n=0.0, d=0.0, nd=0.0, dd=0.0;
79	0	for (Size i=0; i < n_; ++i) {
80	0	const Real x_i = this->xBegin_[i];
81
82	0	if (close_enough(x, x_i)) {
83	0	Real p=0.0;
84	0	for (Size j=0; j < n_; ++j)
85	0	if (i != j) {
86	0	p+=lambda_[j]/(x-this->xBegin_[j])
87	0	*(this->yBegin_[j] - this->yBegin_[i]);
88	0	}
89	0	return p/lambda_[i];
90	0	}
91
92	0	const Real alpha = lambda_[i]/(x-x_i);
93	0	const Real alphad = -alpha/(x-x_i);
94	0	n += alpha * this->yBegin_[i];
95	0	d += alpha;
96	0	nd += alphad * this->yBegin_[i];
97	0	dd += alphad;
98	0	}
99	0	return (nd * d - n * dd)/(d*d);
100	0	} Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::derivative(double) const Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::derivative(double) const
101
102	0	Real primitive(Real) const override {
103	0	QL_FAIL("LagrangeInterpolation primitive is not implemented");
104	0	} Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::primitive(double) const Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::primitive(double) const
105
106	0	Real secondDerivative(Real) const override {
107	0	QL_FAIL("LagrangeInterpolation secondDerivative "
108	0	"is not implemented");
109	0	} Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::secondDerivative(double) const Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::secondDerivative(double) const
110
111	0	Real value(const Array& y, Real x) const override { return _value(y.begin(), x); } Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::value(QuantLib::Array const&, double) const Unexecuted instantiation: QuantLib::detail::LagrangeInterpolationImpl<double, double>::value(QuantLib::Array const&, double) const
112
113		private:
114		template <class Iterator>
115	0	Real _value(const Iterator& yBegin, Real x) const {
116
117	0	const Real eps = 10QL_EPSILONstd::abs(x);
118	0	const auto iter = std::lower_bound(
119	0	this->xBegin_, this->xEnd_, x - eps);
120	0	if (iter != this->xEnd_ && *iter - x < eps) {
121	0	return yBegin[std::distance(this->xBegin_, iter)];
122	0	}
123
124	0	Real n = 0.0, d = 0.0;
125	0	for (Size i = 0; i < n_; ++i) {
126	0	const Real alpha = lambda_[i] / (x - this->xBegin_[i]);
127	0	n += alpha * yBegin[i];
128	0	d += alpha;
129	0	}
130	0	return n / d;
131	0	} Unexecuted instantiation: double QuantLib::detail::LagrangeInterpolationImpl<double const, double const>::_value<double const>(double const const&, double) const Unexecuted instantiation: double QuantLib::detail::LagrangeInterpolationImpl<double, double>::_value<double>(double const&, double) const Unexecuted instantiation: double QuantLib::detail::LagrangeInterpolationImpl<double, double>::_value<double const>(double const const&, double) const
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>
145		LagrangeInterpolation(const I1& xBegin, const I1& xEnd,
146	0	const I2& yBegin) {
147	0	impl_ = ext::make_shared<detail::LagrangeInterpolationImpl<I1,I2> >(
148	0	xBegin, xEnd, yBegin);
149	0	impl_->update();
150	0	} Unexecuted instantiation: QuantLib::LagrangeInterpolation::LagrangeInterpolation<double const, double const>(double const* const&, double const* const&, double const* const&) Unexecuted instantiation: QuantLib::LagrangeInterpolation::LagrangeInterpolation<double, double>(double* const&, double* const&, double* const&)
151
152		// interpolate with new set of y values for a new x value
153	0	Real value(const Array& y, Real x) const {
154	0	return ext::dynamic_pointer_cast<detail::UpdatedYInterpolation>
155	0	(impl_)->value(y, x);
156	0	}
157		};
158
159		}
160
161		#endif

Coverage Report

Created: 2025-08-05 06:45