/src/quantlib/ql/math/interpolations/mixedinterpolation.hpp
Line | Count | Source |
1 | | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ |
2 | | |
3 | | /* |
4 | | Copyright (C) 2010 Ferdinando Ametrano |
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 | | <https://www.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 mixedinterpolation.hpp |
21 | | \brief mixed interpolation between discrete points |
22 | | */ |
23 | | |
24 | | #ifndef quantlib_mixed_interpolation_hpp |
25 | | #define quantlib_mixed_interpolation_hpp |
26 | | |
27 | | #include <ql/math/interpolations/linearinterpolation.hpp> |
28 | | #include <ql/math/interpolations/cubicinterpolation.hpp> |
29 | | #include <ql/utilities/dataformatters.hpp> |
30 | | |
31 | | namespace QuantLib { |
32 | | |
33 | | namespace detail { |
34 | | template <class I1, class I2, class SwitchFn> |
35 | | class MixedInterpolationImpl; |
36 | | } |
37 | | |
38 | | struct MixedInterpolation { |
39 | | enum Behavior { |
40 | | ShareRanges, /*!< Define both interpolations over the |
41 | | whole range defined by the passed |
42 | | iterators. This is the default |
43 | | behavior. */ |
44 | | SplitRanges /*!< Define the first interpolation over the |
45 | | first part of the range, and the second |
46 | | interpolation over the second part. */ |
47 | | }; |
48 | | }; |
49 | | |
50 | | //! mixed linear/cubic interpolation between discrete points |
51 | | // |
52 | | // When using SplitRanges one can set leftC to FirstDerivative and |
53 | | // leftConditionValue to Null to match the first derivatives at the |
54 | | // switch point. |
55 | | /*! \ingroup interpolations |
56 | | \warning See the Interpolation class for information about the |
57 | | required lifetime of the underlying data. |
58 | | */ |
59 | | class MixedLinearCubicInterpolation : public Interpolation { |
60 | | public: |
61 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
62 | | template <class I1, class I2> |
63 | | MixedLinearCubicInterpolation(const I1& xBegin, const I1& xEnd, |
64 | | const I2& yBegin, Size n, |
65 | | MixedInterpolation::Behavior behavior, |
66 | | CubicInterpolation::DerivativeApprox da, |
67 | | bool monotonic, |
68 | | CubicInterpolation::BoundaryCondition leftC, |
69 | | Real leftConditionValue, |
70 | | CubicInterpolation::BoundaryCondition rightC, |
71 | | Real rightConditionValue, |
72 | | bool update = true) { |
73 | | bool matchDerivatives = leftC == CubicInterpolation::FirstDerivative && |
74 | | leftConditionValue == Null<Real>(); |
75 | | QL_REQUIRE(!matchDerivatives || behavior == MixedInterpolation::SplitRanges, |
76 | | "matching derivatives is only supported with SplitRanges"); |
77 | | |
78 | | auto switchFn = [=](Interpolation& left, Interpolation& right, Real x) { |
79 | | if (!matchDerivatives) return; |
80 | | // NOLINTNEXTLINE(cppcoreguidelines-pro-type-static-cast-downcast) |
81 | | auto& cubicImpl = static_cast<detail::CubicInterpolationBaseImpl&>(*right.impl_); |
82 | | cubicImpl.leftValue_ = left.derivative(x, true); |
83 | | }; |
84 | | impl_ = ext::make_shared<detail::MixedInterpolationImpl<I1, I2, decltype(switchFn)>>( |
85 | | xBegin, xEnd, yBegin, n, behavior, |
86 | | Linear(), |
87 | | Cubic(da, monotonic, |
88 | | leftC, leftConditionValue, |
89 | | rightC, rightConditionValue), |
90 | | std::move(switchFn)); |
91 | | if (update) |
92 | | impl_->update(); |
93 | | } |
94 | | }; |
95 | | |
96 | | //! mixed linear/cubic interpolation factory and traits |
97 | | /*! \ingroup interpolations */ |
98 | | class MixedLinearCubic { |
99 | | public: |
100 | | MixedLinearCubic(Size n, |
101 | | MixedInterpolation::Behavior behavior, |
102 | | CubicInterpolation::DerivativeApprox da, |
103 | | bool monotonic = true, |
104 | | CubicInterpolation::BoundaryCondition leftCondition |
105 | | = CubicInterpolation::SecondDerivative, |
106 | | Real leftConditionValue = 0.0, |
107 | | CubicInterpolation::BoundaryCondition rightCondition |
108 | | = CubicInterpolation::SecondDerivative, |
109 | | Real rightConditionValue = 0.0) |
110 | | : n_(n), behavior_(behavior), da_(da), monotonic_(monotonic), |
111 | | leftType_(leftCondition), rightType_(rightCondition), |
112 | 0 | leftValue_(leftConditionValue), rightValue_(rightConditionValue) {} |
113 | | template <class I1, class I2> |
114 | | Interpolation interpolate(const I1& xBegin, const I1& xEnd, |
115 | | const I2& yBegin, bool update = true) const { |
116 | | return MixedLinearCubicInterpolation(xBegin, xEnd, |
117 | | yBegin, n_, behavior_, |
118 | | da_, monotonic_, |
119 | | leftType_, leftValue_, |
120 | | rightType_, rightValue_, |
121 | | update); |
122 | | } |
123 | | // fix below |
124 | | static const bool global = true; |
125 | | static const Size requiredPoints = 3; |
126 | | private: |
127 | | Size n_; |
128 | | MixedInterpolation::Behavior behavior_; |
129 | | CubicInterpolation::DerivativeApprox da_; |
130 | | bool monotonic_; |
131 | | CubicInterpolation::BoundaryCondition leftType_, rightType_; |
132 | | Real leftValue_, rightValue_; |
133 | | }; |
134 | | |
135 | | // convenience classes |
136 | | |
137 | | class MixedLinearCubicNaturalSpline : public MixedLinearCubicInterpolation { |
138 | | public: |
139 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
140 | | template <class I1, class I2> |
141 | | MixedLinearCubicNaturalSpline(const I1& xBegin, const I1& xEnd, |
142 | | const I2& yBegin, Size n, |
143 | | MixedInterpolation::Behavior behavior |
144 | | = MixedInterpolation::ShareRanges) |
145 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
146 | | CubicInterpolation::Spline, false, |
147 | | CubicInterpolation::SecondDerivative, 0.0, |
148 | | CubicInterpolation::SecondDerivative, 0.0) {} |
149 | | }; |
150 | | |
151 | | class MixedLinearMonotonicCubicNaturalSpline : public MixedLinearCubicInterpolation { |
152 | | public: |
153 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
154 | | template <class I1, class I2> |
155 | | MixedLinearMonotonicCubicNaturalSpline(const I1& xBegin, const I1& xEnd, |
156 | | const I2& yBegin, Size n, |
157 | | MixedInterpolation::Behavior behavior |
158 | | = MixedInterpolation::ShareRanges) |
159 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
160 | | CubicInterpolation::Spline, true, |
161 | | CubicInterpolation::SecondDerivative, 0.0, |
162 | | CubicInterpolation::SecondDerivative, 0.0) {} |
163 | | }; |
164 | | |
165 | | class MixedLinearKrugerCubic : public MixedLinearCubicInterpolation { |
166 | | public: |
167 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
168 | | template <class I1, class I2> |
169 | | MixedLinearKrugerCubic(const I1& xBegin, const I1& xEnd, |
170 | | const I2& yBegin, Size n, |
171 | | MixedInterpolation::Behavior behavior |
172 | | = MixedInterpolation::ShareRanges) |
173 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
174 | | CubicInterpolation::Kruger, false, |
175 | | CubicInterpolation::SecondDerivative, 0.0, |
176 | | CubicInterpolation::SecondDerivative, 0.0) {} |
177 | | }; |
178 | | |
179 | | class MixedLinearFritschButlandCubic : public MixedLinearCubicInterpolation { |
180 | | public: |
181 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
182 | | template <class I1, class I2> |
183 | | MixedLinearFritschButlandCubic(const I1& xBegin, const I1& xEnd, |
184 | | const I2& yBegin, Size n, |
185 | | MixedInterpolation::Behavior behavior |
186 | | = MixedInterpolation::ShareRanges) |
187 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
188 | | CubicInterpolation::FritschButland, false, |
189 | | CubicInterpolation::SecondDerivative, 0.0, |
190 | | CubicInterpolation::SecondDerivative, 0.0) {} |
191 | | }; |
192 | | |
193 | | class MixedLinearParabolic : public MixedLinearCubicInterpolation { |
194 | | public: |
195 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
196 | | template <class I1, class I2> |
197 | | MixedLinearParabolic(const I1& xBegin, const I1& xEnd, |
198 | | const I2& yBegin, Size n, |
199 | | MixedInterpolation::Behavior behavior |
200 | | = MixedInterpolation::ShareRanges) |
201 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
202 | | CubicInterpolation::Parabolic, false, |
203 | | CubicInterpolation::SecondDerivative, 0.0, |
204 | | CubicInterpolation::SecondDerivative, 0.0) {} |
205 | | }; |
206 | | |
207 | | class MixedLinearMonotonicParabolic : public MixedLinearCubicInterpolation { |
208 | | public: |
209 | | /*! \pre the \f$ x \f$ values must be sorted. */ |
210 | | template <class I1, class I2> |
211 | | MixedLinearMonotonicParabolic(const I1& xBegin, const I1& xEnd, |
212 | | const I2& yBegin, Size n, |
213 | | MixedInterpolation::Behavior behavior |
214 | | = MixedInterpolation::ShareRanges) |
215 | | : MixedLinearCubicInterpolation(xBegin, xEnd, yBegin, n, behavior, |
216 | | CubicInterpolation::Parabolic, true, |
217 | | CubicInterpolation::SecondDerivative, 0.0, |
218 | | CubicInterpolation::SecondDerivative, 0.0) {} |
219 | | }; |
220 | | |
221 | | namespace detail { |
222 | | |
223 | | template <class I1, class I2, class SwitchFn> |
224 | | class MixedInterpolationImpl final |
225 | | : public Interpolation::templateImpl<I1, I2> { |
226 | | public: |
227 | | template <class Interpolator1, class Interpolator2> |
228 | | MixedInterpolationImpl(const I1& xBegin, const I1& xEnd, |
229 | | const I2& yBegin, Size n, |
230 | | MixedInterpolation::Behavior behavior, |
231 | | const Interpolator1& factory1, |
232 | | const Interpolator2& factory2, |
233 | | SwitchFn switchFn) |
234 | | : Interpolation::templateImpl<I1, I2>(xBegin, xEnd, yBegin, 1), |
235 | | switchFn_(std::move(switchFn)) { |
236 | | // the switch point xBegin_ + n is dereferenced in update() and |
237 | | // value(), so it must be a valid element; for SplitRanges the |
238 | | // first segment additionally reads xBegin_ + n + 1 |
239 | | Size maxN = static_cast<Size>(xEnd - xBegin) - 1; |
240 | | // This only checks that we pass valid iterators into interpolate() |
241 | | // calls below. The calls themselves check requiredPoints for each |
242 | | // of the segments. |
243 | | QL_REQUIRE(n <= maxN, "n is too large (" << n << " > " << maxN << ")"); |
244 | | |
245 | | xBegin2_ = this->xBegin_ + n; |
246 | | |
247 | | switch (behavior) { |
248 | | case MixedInterpolation::ShareRanges: |
249 | | interpolation1_ = detail::interpolateWithoutUpdate( |
250 | | factory1, this->xBegin_, this->xEnd_, this->yBegin_); |
251 | | interpolation2_ = detail::interpolateWithoutUpdate( |
252 | | factory2, this->xBegin_, this->xEnd_, this->yBegin_); |
253 | | break; |
254 | | case MixedInterpolation::SplitRanges: |
255 | | interpolation1_ = detail::interpolateWithoutUpdate( |
256 | | factory1, this->xBegin_, this->xBegin2_ + 1, this->yBegin_); |
257 | | interpolation2_ = detail::interpolateWithoutUpdate( |
258 | | factory2, this->xBegin2_, this->xEnd_, this->yBegin_ + n); |
259 | | break; |
260 | | default: |
261 | | QL_FAIL("unknown mixed-interpolation behavior: " << behavior); |
262 | | } |
263 | | } |
264 | | |
265 | | void update() override { |
266 | | interpolation1_.update(); |
267 | | switchFn_(interpolation1_, interpolation2_, *xBegin2_); |
268 | | interpolation2_.update(); |
269 | | } |
270 | | Real value(Real x) const override { |
271 | | if (x<*xBegin2_) |
272 | | return interpolation1_(x, true); |
273 | | return interpolation2_(x, true); |
274 | | } |
275 | | Real primitive(Real x) const override { |
276 | | if (x<*xBegin2_) |
277 | | return interpolation1_.primitive(x, true); |
278 | | return interpolation2_.primitive(x, true) - |
279 | | interpolation2_.primitive(*xBegin2_, true) + |
280 | | interpolation1_.primitive(*xBegin2_, true); |
281 | | } |
282 | | Real derivative(Real x) const override { |
283 | | if (x<*xBegin2_) |
284 | | return interpolation1_.derivative(x, true); |
285 | | return interpolation2_.derivative(x, true); |
286 | | } |
287 | | Real secondDerivative(Real x) const override { |
288 | | if (x<*xBegin2_) |
289 | | return interpolation1_.secondDerivative(x, true); |
290 | | return interpolation2_.secondDerivative(x, true); |
291 | | } |
292 | | private: |
293 | | I1 xBegin2_; |
294 | | Interpolation interpolation1_, interpolation2_; |
295 | | SwitchFn switchFn_; |
296 | | }; |
297 | | |
298 | | } |
299 | | |
300 | | } |
301 | | |
302 | | #endif |