LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
Interpolate.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
25/*
26 * Interpolate values for a set of x,y vector<>s
27 */
28#include <limits>
29#include <algorithm>
30#include <cmath>
31#include <map>
32#include "boost/format.hpp"
33#include <memory>
34#include "gsl/gsl_errno.h"
35#include "gsl/gsl_interp.h"
36#include "ndarray.h"
37#include "lsst/pex/exceptions.h"
39
40namespace lsst {
41namespace afw {
42namespace math {
43
44namespace {
46 std::vector<double> const &y) {
47 if (x.size() != y.size()) {
48 throw LSST_EXCEPT(
49 pex::exceptions::InvalidParameterError,
50 str(boost::format("Dimensions of x and y must match; %ul != %ul") % x.size() % y.size()));
51 }
52 std::size_t const len = x.size();
53 if (len == 0) {
54 throw LSST_EXCEPT(pex::exceptions::OutOfRangeError, "You must provide at least 1 point");
55 } else if (len == 1) {
56 return std::make_pair(x, y);
57 }
58
59 std::vector<double> recentered_x(len + 1);
60 std::vector<double> recentered_y(len + 1);
61
62 recentered_x[0] = 0.5 * (3 * x[0] - x[1]);
63 recentered_y[0] = y[0];
64
65 for (std::size_t i = 0, j = 1; i < len - 1; ++i, ++j) {
66 recentered_x[j] = 0.5 * (x[i] + x[i + 1]);
67 recentered_y[j] = 0.5 * (y[i] + y[i + 1]);
68 }
69 recentered_x[len] = 0.5 * (3 * x[len - 1] - x[len - 2]);
70 recentered_y[len] = y[len - 1];
71
72 return std::make_pair(recentered_x, recentered_y);
73}
74} // namespace
75
79 Interpolate::Style const style);
80
81public:
82 ~InterpolateConstant() override = default;
83 double interpolate(double const x) const override;
84
85private:
87 std::vector<double> const &y,
88 Interpolate::Style const style
89 )
90 : Interpolate(recenter(x, y)), _old(_x.begin()) {}
91 mutable std::vector<double>::const_iterator _old; // last position we found xInterp at
92};
93
95double InterpolateConstant::interpolate(double const xInterp // the value we want to interpolate to
96 ) const {
97 //
98 // Look for the interval wherein lies xInterp. We could naively use std::upper_bound, but that requires a
99 // logarithmic time lookup so we'll cache the previous answer in _old -- this is a good idea if people
100 // usually call this routine repeatedly for a range of x
101 //
102 // We start by searching up from _old
103 //
104 if (xInterp < *_old) { // We're to the left of the cache
105 if (_old == _x.begin()) { // ... actually off the array
106 return _y[0];
107 }
108 _old = _x.begin(); // reset the cached point to the start of the array
109 } else { // see if we're still in the same interval
110 if (_old < _x.end() - 1 and xInterp < *(_old + 1)) { // we are, so we're done
111 return _y[_old - _x.begin()];
112 }
113 }
114 // We're to the right of the cached point and not in the same inverval, so search up from _old
115 std::vector<double>::const_iterator low = std::upper_bound(_old, _x.end(), xInterp);
116 //
117 // Did that work?
118 if (low == _old && _old != _x.begin()) {
119 // No. Sigh. Search the entire range.
120 low = std::upper_bound(_x.begin(), low + 1, xInterp);
121 }
122 //
123 // OK, we've found the right interval. Return the desired value, being careful at the ends
124 //
125 if (low == _x.end()) {
126 return _y[_y.size() - 1];
127 } else if (low == _x.begin()) {
128 return _y[0];
129 } else {
130 --low;
131 _old = low;
132 return _y[low - _x.begin()];
133 }
134}
135
136namespace {
137/*
138 * Conversion function to switch an Interpolate::Style to a gsl_interp_type.
139 */
140::gsl_interp_type const *styleToGslInterpType(Interpolate::Style const style) {
141 switch (style) {
144 "CONSTANT interpolation not supported.");
146 return ::gsl_interp_linear;
148 return ::gsl_interp_cspline;
150 return ::gsl_interp_cspline;
152 return ::gsl_interp_cspline_periodic;
154 return ::gsl_interp_akima;
156 return ::gsl_interp_akima_periodic;
159 "I am unable to make an interpolator of type UNKNOWN");
162 str(boost::format("You can't get here: style == %") % style));
163 }
164 // don't use default statement to let compiler check switch coverage
166 str(boost::format("You can't get here: style == %") % style));
167}
168} // namespace
169
172 std::vector<double> const &y,
173 Interpolate::Style const style);
174
175public:
176 ~InterpolateGsl() override;
177 double interpolate(double const x) const override;
178
179private:
181 Interpolate::Style const style);
182
183 ::gsl_interp_type const *_interpType;
184 ::gsl_interp_accel *_acc;
185 ::gsl_interp *_interp;
186};
187
188InterpolateGsl::InterpolateGsl(std::vector<double> const &x,
189 std::vector<double> const &y,
190 Interpolate::Style const style
191 )
192 : Interpolate(x, y), _interpType(styleToGslInterpType(style)) {
193 // Turn the gsl error handler off, we want to use our own exceptions
194 ::gsl_set_error_handler_off();
195
196 _acc = ::gsl_interp_accel_alloc();
197 if (!_acc) {
198 throw std::bad_alloc();
199 }
200
201 _interp = ::gsl_interp_alloc(_interpType, _y.size());
202 if (!_interp) {
204 str(boost::format("Failed to initialise spline for type %s, length %d") %
205 _interpType->name % _y.size()));
206 }
207 // Note, "x" and "y" are vector<double>; gsl_inter_init requires double[].
208 // The &(x[0]) here is valid because std::vector guarantees that the values are
209 // stored contiguously in memory (for types other than bool); C++0X 23.3.6.1 for
210 // those of you reading along.
211 int const status = ::gsl_interp_init(_interp, &x[0], &y[0], _y.size());
212 if (status != 0) {
213 throw LSST_EXCEPT(
215 str(boost::format("gsl_interp_init failed: %s [%d]") % ::gsl_strerror(status) % status));
216 }
217}
218
220 ::gsl_interp_free(_interp);
221 ::gsl_interp_accel_free(_acc);
222}
223
224double InterpolateGsl::interpolate(double const xInterp) const {
225 // NaNs interpolate as NaN.
226 if (std::isnan(xInterp)) {
227 return xInterp;
228 }
229 // New GSL versions refuse to extrapolate.
230 // gsl_interp_init() requires x to be ordered, so can just check
231 // the array endpoints for out-of-bounds.
232 if ((xInterp < _x.front() || (xInterp > _x.back()))) {
233 // do our own quadratic extrapolation.
234 // (GSL only provides first and second derivative functions)
235 /* could also just fail via:
236 throw LSST_EXCEPT(
237 pex::exceptions::InvalidParameterError,
238 (boost::format("Interpolation point %f outside range [%f, %f]")
239 % x % _x.front() % _x.back()).str()
240 );
241 */
242 double x0, y0;
243 if (xInterp < _x.front()) {
244 x0 = _x.front();
245 y0 = _y.front();
246 } else {
247 x0 = _x.back();
248 y0 = _y.back();
249 }
250 // first derivative at endpoint
251 double d = ::gsl_interp_eval_deriv(_interp, &_x[0], &_y[0], x0, _acc);
252 // second derivative at endpoint
253 double d2 = ::gsl_interp_eval_deriv2(_interp, &_x[0], &_y[0], x0, _acc);
254 return y0 + (xInterp - x0) * d + 0.5 * (xInterp - x0) * (xInterp - x0) * d2;
255 }
256 assert(xInterp >= _x.front());
257 assert(xInterp <= _x.back());
258 return ::gsl_interp_eval(_interp, &_x[0], &_y[0], xInterp, _acc);
259}
260
262 static std::map<std::string, Interpolate::Style> gslInterpTypeStrings;
263 if (gslInterpTypeStrings.empty()) {
264 gslInterpTypeStrings["CONSTANT"] = Interpolate::CONSTANT;
265 gslInterpTypeStrings["LINEAR"] = Interpolate::LINEAR;
266 gslInterpTypeStrings["CUBIC_SPLINE"] = Interpolate::CUBIC_SPLINE;
267 gslInterpTypeStrings["NATURAL_SPLINE"] = Interpolate::NATURAL_SPLINE;
268 gslInterpTypeStrings["CUBIC_SPLINE_PERIODIC"] = Interpolate::CUBIC_SPLINE_PERIODIC;
269 gslInterpTypeStrings["AKIMA_SPLINE"] = Interpolate::AKIMA_SPLINE;
270 gslInterpTypeStrings["AKIMA_SPLINE_PERIODIC"] = Interpolate::AKIMA_SPLINE_PERIODIC;
271 }
272
273 if (gslInterpTypeStrings.find(style) == gslInterpTypeStrings.end()) {
274 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Interp style not found: " + style);
275 }
276 return gslInterpTypeStrings[style];
277}
278
280 if (n < 1) {
281 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "n must be greater than 0");
282 } else if (n > 4) {
284 } else {
286 if (styles.empty()) {
287 styles.resize(5);
288
289 styles[0] = Interpolate::UNKNOWN; // impossible to reach as we check for n < 1
290 styles[1] = Interpolate::CONSTANT;
291 styles[2] = Interpolate::LINEAR;
292 styles[3] = Interpolate::CUBIC_SPLINE;
293 styles[4] = Interpolate::CUBIC_SPLINE;
294 }
295 return styles[n];
296 }
297}
298
300 size_t const num = x.size();
301 std::vector<double> out(num);
302 for (size_t i = 0; i < num; ++i) {
303 out[i] = interpolate(x[i]);
304 }
305 return out;
306}
307
308ndarray::Array<double, 1> Interpolate::interpolate(ndarray::Array<double const, 1> const &x) const {
309 int const num = x.getShape()[0];
310 ndarray::Array<double, 1> out = ndarray::allocate(ndarray::makeVector(num));
311 for (int i = 0; i < num; ++i) {
312 std::cout << "Interpolating " << x[i] << std::endl;
313 out[i] = interpolate(x[i]);
314 }
315 return out;
316}
317
319 static std::vector<int> minPoints;
320 if (minPoints.empty()) {
322 minPoints[Interpolate::CONSTANT] = 1;
323 minPoints[Interpolate::LINEAR] = 2;
324 minPoints[Interpolate::NATURAL_SPLINE] = 3;
325 minPoints[Interpolate::CUBIC_SPLINE] = 3;
327 minPoints[Interpolate::AKIMA_SPLINE] = 5;
329 }
330
331 if (style >= 0 && style < Interpolate::NUM_STYLES) {
332 return minPoints[style];
333 } else {
334 throw LSST_EXCEPT(
336 str(boost::format("Style %d is out of range 0..%d") % style % (Interpolate::NUM_STYLES - 1)));
337 }
338}
339
341
342 Interpolate::Style const style)
343 : _x(xy.first), _y(xy.second), _style(style) {
344 ;
345}
346
348 Interpolate::Style const style) {
349 switch (style) {
352 default: // use GSL
353 return std::shared_ptr<Interpolate>(new InterpolateGsl(x, y, style));
354 }
355}
356
357std::shared_ptr<Interpolate> makeInterpolate(ndarray::Array<double const, 1> const &x,
358 ndarray::Array<double const, 1> const &y,
359 Interpolate::Style const style) {
360 return makeInterpolate(std::vector<double>(x.begin(), x.end()), std::vector<double>(y.begin(), y.end()),
361 style);
362}
363} // namespace math
364} // namespace afw
365} // namespace lsst
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
int y
Definition SpanSet.cc:48
T back(T... args)
T begin(T... args)
~InterpolateConstant() override=default
friend std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style)
A factory function to make Interpolate objects.
double interpolate(double const x) const override
friend std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style)
A factory function to make Interpolate objects.
double interpolate(double const x) const override
virtual double interpolate(double const x) const =0
std::vector< double > const _x
Definition Interpolate.h:83
std::vector< double > const _y
Definition Interpolate.h:84
Interpolate(Interpolate const &)=delete
Reports invalid arguments.
Definition Runtime.h:66
Reports errors in the logical structure of the program.
Definition Runtime.h:46
Reports attempts to access elements outside a valid range of indices.
Definition Runtime.h:89
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T empty(T... args)
T end(T... args)
T endl(T... args)
T find(T... args)
T front(T... args)
T isnan(T... args)
T make_pair(T... args)
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for 'n' points.
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style.
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
std::shared_ptr< Interpolate > makeInterpolate(std::vector< double > const &x, std::vector< double > const &y, Interpolate::Style const style=Interpolate::AKIMA_SPLINE)
A factory function to make Interpolate objects.
T resize(T... args)
T size(T... args)
T upper_bound(T... args)