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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
|
// interpolate.hh -- classes for 1D interpolation
// $Id$
//
// ***** design notes *****
// Lagrange_interpolator_1d_uniform
//
//
// prerequisites:
// <assert.h>
// <jt/util.hh>
// <jt/linear_map.hh>
//
//******************************************************************************
//
// ***** design notes *****
//
//
// The basic problem of interpolation is this: given a set of data points
// (x[i],y[i]) sampled from a smooth function y = y(x) , we wish to
// estimate y or one if its derivatives, at arbitrary x values
// (presumably within or close to the range of the x[i] values).
//
// We classify interpolation functions in 4 different ways:
// * Is the interpolation 1-dimensional (the x value are real numbers),
// 2-dimensional (the x values are points in R^2), 3-dimensional, ...?
// The routines in this file are for the 1-dimensional case only.
// * Is the interpolation "local" (the underlying routines, where all
// data points are used), or "global" (search for the right place in
// a large array, then do a local interpolation)?
// * Is the interpolation "uniform" (the x[i] values are uniformly
// spaced, and may be implicitly specified as a linear function of i),
// or non-uniform (the x[i] values are in general non-uniformly
// spaced, and hence must be supplied explicitly to the interpolation
// function along with the y[i] values)?
// * Are doing "function" interpolation (we want to estimate y ), or
// "derivative" interpolation (we want to estimate dy/dx ? (Higher
// derivatives are also possible, but I've never found a use for them.)
//
// We may also want meta-information about the interpolation, for example
// its backward domain of dependence (for a given xout , the set of i
// such that x[i] is needed to interpolate data at xout ) or the Jacobian
// coefficients (the partial derivatives d(interpolated data)/d y[i]).
//
//
// Most interpolators have an "order" parameter, but the terminology
// for this isn't quite standardized. We use the convention that order
// always refers to the order of the interpolating function, i.e.
// linear interpolation is order=1, quadratic is order=2, etc etc.
// Thus the interpolator error is typically one power higher than this
// (eg. linear interpolation has O(dx^2) error).
//
//
// All the actual interpolation functions here take an argument
// end_action . This specifies what action we should take if the
// nominal (centered) local interpolant would need data from (fuzzily)
// outside the domain of the data points. In general there are three
// possibilities for this:
// 'e' ==> Do an error_exit() .
// 'o' ==> Off-center the local interpolant as needed. (This is the
// default.)
// In both of these cases, an error_exit() is still signalled if the
// interpolation position lies (fuzzily) outside the data range.
// '*' ==> Like 'o', but also allow the interpolation position to lie
// (fuzzily) outside the data range. In other words, this case
// accepts *any* interpolation position.
//
// For Lagrange interpolation all three possibilities are valid;
// for spline interpolation the "off-centering as needed" is intrinsic
// to the method , so we take 'e' as a synonym for 'o'.
//
// The default choice is 'o' for all interpolation methods.
//
//
// Some types of interpolatators have to do expensive "setup" computations.
// We'd like to avoid having to redo these for each new interpolation.
// Conceptually, we can divide the interpolation process into 5 steps:
// 1. Set up the type and order of interpolation scheme, and the number
// of data points, and any other similar options such as uniform vs
// nonuniform spacing of the data points' x coordinates.
// 2. Set up the x coordinates of the data points.
// 3. Set up the y coordinates of the data points.
// 4. Set up the x coordinates where interpolated data is desired.
// 5. Do the actual interpolation.
//
// The idea is that earlier steps can preallocate storage arrays and/or
// precompute coefficients to be used by later stages, so repeating the
// later stages may have lower overheads than redoing the entire process.
// The interface here requires that 1 be done first, then allows 2-4
// to be done in any order before 5.
//
//
// Note also that these routines are allowed to just save the pointers
// or references supplied by the caller, i.e. the actual data need not
// be copied. Thus any such data should remain valid throughout the
// lifetime of the interpolator object.
//
//******************************************************************************
namespace interpolation
{
//**************************************
//
// This class provides 1D Lagrange global polynomial interpolation
// for uniformly spaced data.
//
// ... const qualifiers refer to results of interpolate()
//
template <typename fpt>
class Lagrange_uniform_1d
{
public:
// this function sets up the x data values
// and the range of valid [i] and the mapping between [i] and x
void setup_data_x(const linear_map<fpt>& x_map_in);
// this function sets up the y data values
// n.b. indexing off this pointer is [i], i.e. it need not be 0-origin
void setup_data_y(const fpt y_array_in[]);
// this function sets up
// * the x coordinate where interpolation is to be done
// * what sort of off-the-end-of-the-data handling is desired
void setup_interp_x(fpt x_interp_in, char end_action_in = 'o');
// the actual interpolation
// ... requires that
// setup_data_x()
// setup_data_y()
// setup_interp_x()
// have already been called
fpt interpolate();
// this function computes the backwards domain of dependence
// of the interpolation, i.e. the interval of [i] such that
// y[i] is needed to compute interpolate(x,end_action);
// equivalently, this function defines the centering policy
// for interpolate()
//
// ... returns domain-of-dependence interval in [i_lo,i_hi]
// ... returns nominal center point of underlying local interpolant
// (i_center in source code) as function value
// ... requires that
// setup_data_x()
// setup_interp_x()
// have already been called
int bdod(int& i_lo, int& i_hi) const;
// Jacobian d (value returned by interpolate()) / d y[i]
// ... requires that
// setup_data_x()
// setup_interp_x()
// have already been called
fpt Jacobian(int i) const;
// constructor, destructor
// ... compiler-generated no-op destructor is ok
Lagrange_uniform_1d(int order_in);
private:
// we forbid copying and passing by value
// by declaring the copy constructor and assignment operator
// private, but never defining them
Lagrange_uniform_1d(const Lagrange_uniform_1d<fpt> &rhs);
Lagrange_uniform_1d<fpt>& operator=
(const Lagrange_uniform_1d<fpt>& rhs);
private:
const int order;
const int N_interp; // number of points in local
// interpolator, i.e.
// linear = 2 points,
// quadratic = 3 points, etc
bool x_interp_valid; // valid flag for x_interp, end_action
fpt x_interp;
char end_action;
// n.b. non-const pointers to const data!
const linear_map<fpt> *x_map_ptr; // --> client-owned data
const fpt *y_array; // --> client-owned data
};
//**************************************
//
// ***** underlying local interpolants of varying orders *****
//
// ... x here is in units of the grid spacing,
// relative to the center point
// ... these make no checks for data being in range
template <typename fpt>
fpt interpolate_local_order1(const fpt x, const fpt *y_center);
template <typename fpt>
fpt interpolate_local_order2(const fpt x, const fpt *y_center);
template <typename fpt>
fpt interpolate_local_order3(const fpt x, const fpt *y_center);
template <typename fpt>
fpt interpolate_local_order4(const fpt x, const fpt *y_center);
fpt interpolate_local_order5(const fpt x, const fpt *y_center);
fpt interpolate_local_order6(const fpt x, const fpt *y_center);
//**************************************
}; // close namespace interpolation::
|
