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
|
// 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
// 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 routines here all allow 2-5, 3-5, and/or 4-5 to be repeated as needed.
// 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.
//
// One could also imagine doing these steps in the order 12435, i.e.
// interpolating many different sets of date to the same position(s).
// The Cactus interpolator interface supports this nicely, but alas
// the interface here doesn't support this at all. :( :(
//
//
// 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.
//
//******************************************************************************
//
// This class provides 1D Lagrange global polynomial interpolation
// for uniformly spaced data.
//
// ... const qualifiers refer to results of interpolate()
//
namespace interpolation
{
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
void setup_data_y(const fpt y_array_in[]);
// the actual interpolation
// ... needs setup_data_x() and setup_data_y()
// to have already been called
fpt interpolate(fpt x, char end_action = 'o') const;
// 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 (i_center in source code)
// as function value
// ... needs setup_data_x() to have already been called,
// but does *not* need setup_data_y() to have already been called
int bdod(fpt x, char end_action = 'o', int& i_lo, int& i_hi) const;
// Jacobian d interpolate(x,end_action) / d y[i]
// ... needs setup_data_x() to have already been called,
// but does *not* need setup_data_y() to have already been called
fpt Jacobian(fpt x, char end_action = 'o', int i) const;
// constructor, destructor
Lagrange_uniform_1d(int order_in);
~Lagrange_uniform_1d() { }
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
// 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
};
};
|
