1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/*
Copyright (C) 2002 Ferdinando Ametrano
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 ferdinando@ametrano.net
The license is also available online at http://quantlib.org/html/license.html
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.
*/
/*! \file bilinearinterpolation.hpp
\brief bilinear interpolation between discrete points
\fullpath
ql/Math/%bilinearinterpolation.hpp
*/
// $Id: bilinearinterpolation.hpp,v 1.7 2002/03/12 10:55:29 nando Exp $
#ifndef quantlib_bilinear_interpolation_h
#define quantlib_bilinear_interpolation_h
#include <ql/Math/interpolation2D.hpp>
#include <algorithm>
namespace QuantLib {
namespace Math {
//! bilinear interpolation between discrete points
template <class RandomAccessIteratorX,
class RandomAccessIteratorY,
class MatricialData>
class BilinearInterpolation
: public Interpolation2D<RandomAccessIteratorX,
RandomAccessIteratorY,
MatricialData> {
public:
typedef
typename QL_ITERATOR_TRAITS<RandomAccessIteratorX>::value_type
first_argument_type;
typedef
typename QL_ITERATOR_TRAITS<RandomAccessIteratorY>::value_type
second_argument_type;
typedef double result_type;
BilinearInterpolation(
const RandomAccessIteratorX& xBegin,
const RandomAccessIteratorX& xEnd,
const RandomAccessIteratorY& yBegin,
const RandomAccessIteratorY& yEnd,
const MatricialData& data,
bool allowExtrapolation)
: Interpolation2D<RandomAccessIteratorX,
RandomAccessIteratorY,
MatricialData>(
xBegin,xEnd,yBegin,yEnd,data,allowExtrapolation) {}
double operator()(
const first_argument_type& x,
const second_argument_type& y) const;
};
// inline definitions
template <class I1, class I2, class M>
double BilinearInterpolation<I1,I2,M>::operator()(
const BilinearInterpolation<I1,I2,M>::first_argument_type& x,
const BilinearInterpolation<I1,I2,M>::second_argument_type& y)
const {
I1 i; // column
if (x < *xBegin_) {
QL_REQUIRE(allowExtrapolation_,
"BilinearInterpolation::operator() : "
"extrapolation not allowed "
"[x<xMin]");
i = xBegin_;
} else if (x > *(xEnd_-1)) {
QL_REQUIRE(allowExtrapolation_,
"BilinearInterpolation::operator() : "
"extrapolation not allowed "
"[x>xMax]");
i = xEnd_-2;
} else
i = std::upper_bound(xBegin_,xEnd_-1,x)-1;
I2 j; // row
if (y < *yBegin_) {
QL_REQUIRE(allowExtrapolation_,
"BilinearInterpolation::operator() : "
"extrapolation not allowed "
"[y<yMin]");
j = yBegin_;
} else if (y > *(yEnd_-1)) {
QL_REQUIRE(allowExtrapolation_,
"BilinearInterpolation::operator() : "
"extrapolation not allowed "
"[y>yMax]");
j = yEnd_-2;
} else
j = std::upper_bound(yBegin_,yEnd_-1,y)-1;
double z1=data_[j-yBegin_] [i-xBegin_];
double z2=data_[j-yBegin_] [i-xBegin_+1];
double z3=data_[j-yBegin_+1] [i-xBegin_];
double z4=data_[j-yBegin_+1] [i-xBegin_+1];
double t=(x-*i)/(*(i+1)-*i);
double u=(y-*j)/(*(j+1)-*j);
return (1.0-t) * (1.0-u) * z1+
t * (1.0-u) * z2+
(1.0-t) * u * z3+
t * u * z4;
}
}
}
#endif
|
