QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
multicubicspline.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) 2003, 2004 Roman Gitlin
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/*! \file multicubicspline.hpp
21 \brief N-dimensional cubic spline interpolation between discrete points
22*/
23
24#ifndef quantlib_multi_cubic_spline_hpp
25#define quantlib_multi_cubic_spline_hpp
26
27#include <ql/errors.hpp>
28#include <ql/types.hpp>
29#include <algorithm>
30#include <functional>
31#include <utility>
32#include <vector>
33
34namespace QuantLib {
35
36 namespace detail {
37
38 // data structures
39
40 typedef std::vector<std::vector<Real> > SplineGrid;
41
42 // different termination markers are necessary
43 // to maintain separation of possibly overlapping types
44 struct EmptyArg {}; // arg_type termination marker
45 struct EmptyRes {}; // res_type termination marker
46 struct EmptyDim {}; // size_t termination marker
47
48 template<class X> struct DataTable {
49 DataTable(const std::vector<Size>::const_iterator &i) {
50 std::vector<X> temp(*i, X(i + 1));
51 data_table_.swap(temp);
52 }
53 DataTable(const SplineGrid::const_iterator &i) {
54 std::vector<X> temp(i->size(), X(i + 1));
55 data_table_.swap(temp);
56 }
57 template<class U> DataTable(const std::vector<U> &v) {
58 DataTable temp(v.begin());
59 data_table_.swap(temp.data_table_);
60 }
61 Size size() const {
62 return data_table_.size();
63 }
64 const X &operator[](Size n) const {return data_table_[n];}
65 X &operator[](Size n) {return data_table_[n];}
66 std::vector<X> data_table_;
67 };
68
69 template<> struct DataTable<Real> {
71 DataTable(const std::vector<Size>::const_iterator& i)
72 : data_table_(*i) {}
73 DataTable(const SplineGrid::const_iterator &i)
74 : data_table_(i->size()) {}
75 template<class U> DataTable(const std::vector<U> &v) {
76 DataTable temp(v.begin());
77 data_table_.swap(temp.data_table_);
78 }
79 Size size() const {
80 return data_table_.size();
81 }
82 Real operator[](Size n) const {return data_table_[n];}
84 std::vector<Real> data_table_;
85 };
86
88
89 template<class X, class Y> struct Data {
91 : first(X()), second(Y()) {}
92 Data(const SplineGrid::const_iterator &i)
93 : first(*i), second(i + 1) {}
95 : first(v[0]), second(v.begin()+1) {}
96 void swap(Data<X, Y> &d) noexcept {
97 first.swap(d.first);
98 second.swap(d.second);
99 }
102 };
103
104 template<> struct Data<std::vector<Real>, EmptyArg> {
105 Data() = default;
106 Data(const SplineGrid::const_iterator &i)
107 : first(*i) {}
109 : first(v[0]) {}
110 Data(std::vector<Real> v) : first(std::move(v)) {}
111 void swap(Data<std::vector<Real>, EmptyArg> &d) noexcept {
112 first.swap(d.first);
113 }
114 Real operator[](Size n) const {return first[n];}
116 std::vector<Real> first;
118 };
119
121
122 template<class X, class Y> struct Point {
123 typedef X data_type;
125 : first(data_type()), second(Y()) {}
126 Point(const std::vector<Real>::const_iterator &i)
127 : first(*i), second(i + 1) {}
128 Point(const std::vector<Real> &v)
129 : first(v[0]), second(v.begin()+1) {}
130 Point(const SplineGrid::const_iterator &i)
131 : first(i->size()), second(i + 1) {}
132 Point(const SplineGrid &grid)
133 : first(grid[0].size()), second(grid.begin()+1) {}
134 operator data_type() const {
135 return first;
136 }
137 data_type operator[](Size n) const { return n != 0U ? second[n - 1] : first; }
138 data_type& operator[](Size n) { return n != 0U ? second[n - 1] : first; }
141 };
142
143 template<> struct Point<Real, EmptyArg> {
146 : first(s) {}
147 Point(const std::vector<Real>::const_iterator &i)
148 : first(*i) {}
149 Point(const std::vector<Real> &v)
150 : first(v[0]) {}
151 operator data_type() const {return first;}
153 QL_REQUIRE(n == 0, "operator[] : access violation");
154 return first;
155 }
157 QL_REQUIRE(n == 0, "operator[] : access violation");
158 return first;
159 }
162 };
163
165
166 template<> struct Point<Real, EmptyRes> {
169 : first(data_type()) {}
171 : first(s) {}
172 operator data_type() const {return first;}
173 const data_type &operator[](Size n) const {
174 QL_REQUIRE(n == 0, "operator[] : access violation");
175 return first;
176 }
178 QL_REQUIRE(n == 0, "operator[] : access violation");
179 return first;
180 }
183 };
184
186
187 template<> struct Point<Size, EmptyDim> {
190 : first(data_type()) {}
192 : first(s) {}
193 operator data_type() const {return first;}
195 QL_REQUIRE(n == 0, "operator[] : access violation");
196 return first;
197 }
199 QL_REQUIRE(n == 0, "operator[] : access violation");
200 return first;
201 }
204 };
205
207
208 template<> struct Point<base_data_table, EmptyRes> {
210 Point(data_type s) : first(std::move(s)) {}
211 Point(const SplineGrid::const_iterator &i)
212 : first(i->size()) {}
213 Point(const SplineGrid &grid)
214 : first(grid[0].size()) {}
215 Real operator[](Size n) const {return first[n];}
219 };
220
222
223
224 // cubic spline iplementations
225
226 // no heap memory is allocated
227 // in any of the recursive calls
229 public:
235 base_cubic_spline(const data &d, const data &d2,
236 const data_table& y, data_table &y2,
237 output_data &v) {
238 Size dim = d.first.size();
239 Size j = 1, k = 2, l = 3;
240 result_type &w = ((y2[dim] = y[1]) -= y[0]) /= d[0],
241 &u = ((y2[0] = y[2]) -= y[1]) /= d[1], &t = v[dim];
242 y2[1] = -d[1] / d2[0], v[1] = 6.0 * (u - w) / d2[0];
243 for(; k < dim; u = w, j = k, k = l, ++l) {
244 w = (y[l]-y[k])/d[k];
245 u = (u-w)*6.0;
246 (y2[k] = d[k]) /= ((t = -y2[j]) *= d[j]) -= d2[j];
247 (v[k] = (u += d[j] * v[j])) /= t;
248 }
249 y2[0] = y2[dim] = 0.0;
250 while (k != 0U) {
251 (y2[k-1] *= y2[l-1]) += v[k-1];
252 --k; --l;
253 }
254 }
255 };
256
257 template<class X>
259 public:
263 n_cubic_spline(const data &d, const data &d2,
264 const data_table &y, data_table &y2, output_data &v)
265 : d_(d), d2_(d2), y_(y), y2_(y2), v_(v) {
266 for(Size j = 0, dim = y_.size(); j < dim; ++j)
267 X(d_.second, d2_.second, y_[j], y2_[j], v_.second);
268 }
269 ~n_cubic_spline() = default;
270
271 private:
272 const data &d_, &d2_;
276 };
277
279 public:
288 const return_type &a2, const return_type &b2,
289 const dimensions &i,
290 const data&, const data&,
291 const data_table &y, data_table &y2,
294 result_type &res) {
295 res = a * y[i] + b * y[i + 1] + a2 * y2[i] + b2 * y2[i + 1];
296 }
297 };
298
299 template<class X>
301 public:
308 typedef Point<result_type,
309 typename X::return_type> return_type;
311 const return_type &a2, const return_type &b2,
312 const dimensions &i, const data &d, const data &d2,
313 const data_table &y, data_table &y2, output_data &v,
314 output_data &v1, output_data &v2,
315 result_type& r)
316 : a_(a), b_(b), a2_(a2), b2_(b2), i_(i), d_(d), d2_(d2),
317 y_(y), y2_(y2), v_(v), v1_(v1), v2_(v2) {
318 for(Size j = 0, dim = y_.size(); j < dim; ++j)
320 d_.second, d2_.second, y_[j], y2_[j], v_.second,
323 v1_.first.first, v2_.first.first, v_.first);
326 v1_.first.first, v2_.first.first,
327 v_.first, v1_.first, v2_.first, r);
328 }
329 ~n_cubic_splint() = default;
330
331 private:
332 const return_type &a_, &b_, &a2_, &b2_;
334 const data &d_, &d2_;
338 };
339
355
371
372 template<Size i> struct Int2Type {
375 };
376
377 template<> struct Int2Type<2> {
380 };
381
382 template<> struct Int2Type<3> {
385 };
386
387 template<> struct Int2Type<4> {
390 };
391
392 template<> struct Int2Type<5> {
395 };
396
397 template<> struct Int2Type<6> {
400 };
401
402 template<> struct Int2Type<7> {
405 };
406
407 template<> struct Int2Type<8> {
410 };
411
412 template<> struct Int2Type<9> {
415 };
416
417 template<> struct Int2Type<10> {
420 };
421
422 template<> struct Int2Type<11> {
425 };
426
427 template<> struct Int2Type<12> {
430 };
431
432 template<> struct Int2Type<13> {
435 };
436
437 template<> struct Int2Type<14> {
440 };
441
442 template<> struct Int2Type<15> {
445 };
446
447 }
448
449
450 // Multi-cubic spline
451
453
454 //! N-dimensional cubic spline interpolation between discrete points
455 /*! \test interpolated values are checked against the original
456 function.
457
458 \todo
459 - allow extrapolation as for the other interpolations
460 - investigate if and how to implement Hyman filters and
461 different boundary conditions
462
463 \bug cannot interpolate at the grid points on the boundary
464 surface of the N-dimensional region
465 */
466 template <Size i> class MultiCubicSpline {
469 public:
476 typedef typename c_splint::data data;
478 const std::vector<bool>& ae =
479 std::vector<bool>(20, false))
480 : grid_(grid), y_(y), ae_(ae), v_(grid), v1_(grid),
481 v2_(grid), y2_(grid) {
483 c_spline(d_, d2_, y_, y2_, v_);
484 }
487 c_splint(a_, b_, a2_, b2_, i_, d_, d2_, y_, y2_,
488 v_, v1_, v2_, res_);
489 return res_;
490 }
491 void set_shared_increments() const;
493 private:
496 const std::vector<bool> &ae_;
497 mutable return_type a_, b_, a2_, b2_;
500 mutable dimensions i_;
501 mutable data d_, d2_;
503 };
504
505 // the data is checked and, in case of insufficient number of points,
506 // exception is thrown BEFORE the main body of interpolation begins
507 template <Size i>
509 SplineGrid x(i), y(i);
510 Size k = 0, dim = 0;
511 for(Size j = 0; j < i; k = 0, ++j) {
512 const std::vector<Real> &v = grid_[j];
513 if((dim = v.size() - 1) > 2) {
514 std::vector<Real> tmp1(dim);
515 x[j].swap(tmp1);
516 std::vector<Real> tmp2(dim - 1);
517 y[j].swap(tmp2);
518 for(; k < dim; ++k) {
519 if((x[j][k] = v[k + 1] - v[k]) <= 0.0) break;
520 if (k != 0U)
521 y[j][k - 1] = 2.0 * (v[k + 1] - v[k - 1]);
522 }
523 }
524 QL_REQUIRE(dim >= 3,
525 "Dimension " << j
526 << " : not enough points for interpolation");
527 QL_REQUIRE(k >= dim,
528 "Dimension " << j << " : invalid data");
529 }
530
531 typename c_splint::data tmp1(x), tmp2(y);
532 d_.swap(tmp1);
533 d2_.swap(tmp2);
534 }
535
536 #ifndef __DOXYGEN__
537 // the argument value is checked and, in out of boundaries case,
538 // exception is thrown BEFORE the main body of interpolation begins
539 template <Size i>
541 const typename MultiCubicSpline<i>::argument_type &x) const {
542 for(Size j = 0; j < i; ++j) {
543 Size &k = i_[j], sz = grid_[j].size() - 1;
544 const std::vector<Real> &v = grid_[j];
545 if(x[j] < v[0] || x[j] >= v[sz]) {
546 QL_REQUIRE(ae_[j],
547 "Dimension " << j
548 << ": extrapolation is not allowed.");
549 a_[j] = 1.0, a2_[j] = b_[j] = b2_[j] = 0.0;
550 k = x[j] < v[0] ? 0 : sz;
551 }
552 else {
553 k = v[k] <= x[j] && x[j] < v[k + 1] ? k :
554 std::upper_bound(v.begin(),v.end(),x[j])-v.begin()-1;
555 Real h = v[k + 1] - v[k];
556 a_[j] = (v[k + 1] - x[j]) / h, b_[j] = (x[j] - v[k]) / h;
557 a2_[j] = (a_[j] * a_[j] * a_[j] - a_[j]) * h * h / 6.0,
558 b2_[j] = (b_[j] * b_[j] * b_[j] - b_[j]) * h * h / 6.0;
559 }
560 }
561 }
562 #endif
563
564}
565
566
567#endif
568
N-dimensional cubic spline interpolation between discrete points.
c_splint::argument_type argument_type
c_splint::result_type result_type
c_splint::dimensions dimensions
c_splint::output_data output_data
const std::vector< bool > & ae_
result_type operator()(const argument_type &x) const
void set_shared_coefficients(const argument_type &x) const
c_splint::return_type return_type
MultiCubicSpline(const SplineGrid &grid, const data_table &y, const std::vector< bool > &ae=std::vector< bool >(20, false))
c_splint::data_table data_table
detail::Int2Type< i >::c_splint c_splint
detail::Int2Type< i >::c_spline c_spline
base_cubic_spline(const data &d, const data &d2, const data_table &y, data_table &y2, output_data &v)
base_cubic_splint(const return_type &a, const return_type &b, const return_type &a2, const return_type &b2, const dimensions &i, const data &, const data &, const data_table &y, data_table &y2, output_data &, output_data &, output_data &, result_type &res)
DataTable< typename X::data_table > data_table
n_cubic_spline(const data &d, const data &d2, const data_table &y, data_table &y2, output_data &v)
Data< base_data, typename X::data > data
Point< base_output_data, typename X::output_data > output_data
DataTable< typename X::data_table > data_table
Point< Size, typename X::dimensions > dimensions
Point< Real, typename X::argument_type > argument_type
Data< base_data, typename X::data > data
Point< base_output_data, typename X::output_data > output_data
n_cubic_splint(const return_type &a, const return_type &b, const return_type &a2, const return_type &b2, const dimensions &i, const data &d, const data &d2, const data_table &y, data_table &y2, output_data &v, output_data &v1, output_data &v2, result_type &r)
Point< result_type, typename X::return_type > return_type
Real a_
const DefaultType & t
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Date d
ext::function< Real(Real)> b
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
const VF_R b_
n_cubic_splint< cubic_splint_13 > cubic_splint_14
n_cubic_splint< cubic_splint_05 > cubic_splint_06
Point< Size, EmptyDim > base_dimensions
base_cubic_spline cubic_spline_01
Data< std::vector< Real >, EmptyArg > base_data
n_cubic_splint< cubic_splint_12 > cubic_splint_13
n_cubic_spline< cubic_spline_12 > cubic_spline_13
n_cubic_splint< cubic_splint_02 > cubic_splint_03
n_cubic_splint< cubic_splint_01 > cubic_splint_02
Point< Real, EmptyRes > base_return_type
n_cubic_spline< cubic_spline_02 > cubic_spline_03
n_cubic_splint< cubic_splint_11 > cubic_splint_12
n_cubic_spline< cubic_spline_07 > cubic_spline_08
n_cubic_spline< cubic_spline_10 > cubic_spline_11
n_cubic_splint< cubic_splint_03 > cubic_splint_04
base_cubic_splint cubic_splint_01
n_cubic_spline< cubic_spline_04 > cubic_spline_05
n_cubic_spline< cubic_spline_05 > cubic_spline_06
Point< base_data_table, EmptyRes > base_output_data
n_cubic_splint< cubic_splint_10 > cubic_splint_11
n_cubic_splint< cubic_splint_08 > cubic_splint_09
n_cubic_splint< cubic_splint_14 > cubic_splint_15
Point< Real, EmptyArg > base_arg_type
DataTable< Real > base_data_table
n_cubic_splint< cubic_splint_06 > cubic_splint_07
n_cubic_splint< cubic_splint_09 > cubic_splint_10
n_cubic_splint< cubic_splint_04 > cubic_splint_05
n_cubic_spline< cubic_spline_14 > cubic_spline_15
std::vector< std::vector< Real > > SplineGrid
n_cubic_spline< cubic_spline_03 > cubic_spline_04
n_cubic_splint< cubic_splint_07 > cubic_splint_08
n_cubic_spline< cubic_spline_11 > cubic_spline_12
n_cubic_spline< cubic_spline_09 > cubic_spline_10
n_cubic_spline< cubic_spline_01 > cubic_spline_02
n_cubic_spline< cubic_spline_06 > cubic_spline_07
n_cubic_spline< cubic_spline_08 > cubic_spline_09
n_cubic_spline< cubic_spline_13 > cubic_spline_14
Definition: any.hpp:35
detail::SplineGrid SplineGrid
STL namespace.
ext::shared_ptr< YieldTermStructure > r
ext::shared_ptr< BlackVolTermStructure > v
void swap(Data< std::vector< Real >, EmptyArg > &d) noexcept
void swap(Data< X, Y > &d) noexcept
Data(const SplineGrid &v)
Data(const SplineGrid::const_iterator &i)
DataTable(const std::vector< Size >::const_iterator &i)
DataTable(const std::vector< U > &v)
DataTable(const SplineGrid::const_iterator &i)
DataTable(const std::vector< Size >::const_iterator &i)
const X & operator[](Size n) const
DataTable(const std::vector< U > &v)
DataTable(const SplineGrid::const_iterator &i)
Point(const std::vector< Real >::const_iterator &i)
const data_type & operator[](Size n) const
data_type & operator[](Size n)
Point(const std::vector< Real >::const_iterator &i)
Point(const std::vector< Real > &v)
data_type operator[](Size n) const
Point(const SplineGrid::const_iterator &i)
Point(const SplineGrid &grid)
Custom types.