LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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 <map>
31 #include "boost/format.hpp"
32 #include <memory>
33 #include "gsl/gsl_errno.h"
34 #include "gsl/gsl_interp.h"
35 #include "ndarray.h"
36 #include "lsst/pex/exceptions.h"
38 
39 namespace lsst {
40 namespace afw {
41 namespace math {
42 
43 namespace {
45  std::vector<double> const &y) {
46  if (x.size() != y.size()) {
47  throw LSST_EXCEPT(
48  pex::exceptions::InvalidParameterError,
49  str(boost::format("Dimensions of x and y must match; %ul != %ul") % x.size() % y.size()));
50  }
51  std::size_t const len = x.size();
52  if (len == 0) {
53  throw LSST_EXCEPT(pex::exceptions::OutOfRangeError, "You must provide at least 1 point");
54  } else if (len == 1) {
55  return std::make_pair(x, y);
56  }
57 
58  std::vector<double> recentered_x(len + 1);
59  std::vector<double> recentered_y(len + 1);
60 
61  recentered_x[0] = 0.5 * (3 * x[0] - x[1]);
62  recentered_y[0] = y[0];
63 
64  for (std::size_t i = 0, j = 1; i < len - 1; ++i, ++j) {
65  recentered_x[j] = 0.5 * (x[i] + x[i + 1]);
66  recentered_y[j] = 0.5 * (y[i] + y[i + 1]);
67  }
68  recentered_x[len] = 0.5 * (3 * x[len - 1] - x[len - 2]);
69  recentered_y[len] = y[len - 1];
70 
71  return std::make_pair(recentered_x, recentered_y);
72 }
73 } // namespace
74 
77  std::vector<double> const &y,
78  Interpolate::Style const style);
79 
80 public:
81  ~InterpolateConstant() override = default;
82  double interpolate(double const x) const override;
83 
84 private:
86  std::vector<double> const &y,
87  Interpolate::Style const style
88  )
89  : Interpolate(recenter(x, y)), _old(_x.begin()) {}
90  mutable std::vector<double>::const_iterator _old; // last position we found xInterp at
91 };
92 
94 double InterpolateConstant::interpolate(double const xInterp // the value we want to interpolate to
95  ) const {
96  //
97  // Look for the interval wherein lies xInterp. We could naively use std::upper_bound, but that requires a
98  // logarithmic time lookup so we'll cache the previous answer in _old -- this is a good idea if people
99  // usually call this routine repeatedly for a range of x
100  //
101  // We start by searching up from _old
102  //
103  if (xInterp < *_old) { // We're to the left of the cache
104  if (_old == _x.begin()) { // ... actually off the array
105  return _y[0];
106  }
107  _old = _x.begin(); // reset the cached point to the start of the array
108  } else { // see if we're still in the same interval
109  if (_old < _x.end() - 1 and xInterp < *(_old + 1)) { // we are, so we're done
110  return _y[_old - _x.begin()];
111  }
112  }
113  // We're to the right of the cached point and not in the same inverval, so search up from _old
115  //
116  // Did that work?
117  if (low == _old && _old != _x.begin()) {
118  // No. Sigh. Search the entire range.
119  low = std::upper_bound(_x.begin(), low + 1, xInterp);
120  }
121  //
122  // OK, we've found the right interval. Return the desired value, being careful at the ends
123  //
124  if (low == _x.end()) {
125  return _y[_y.size() - 1];
126  } else if (low == _x.begin()) {
127  return _y[0];
128  } else {
129  --low;
130  _old = low;
131  return _y[low - _x.begin()];
132  }
133 }
134 
135 namespace {
136 /*
137  * Conversion function to switch an Interpolate::Style to a gsl_interp_type.
138  */
139 ::gsl_interp_type const *styleToGslInterpType(Interpolate::Style const style) {
140  switch (style) {
143  "CONSTANT interpolation not supported.");
144  case Interpolate::LINEAR:
145  return ::gsl_interp_linear;
147  return ::gsl_interp_cspline;
149  return ::gsl_interp_cspline;
151  return ::gsl_interp_cspline_periodic;
153  return ::gsl_interp_akima;
155  return ::gsl_interp_akima_periodic;
158  "I am unable to make an interpolator of type UNKNOWN");
161  str(boost::format("You can't get here: style == %") % style));
162  }
163  // don't use default statement to let compiler check switch coverage
165  str(boost::format("You can't get here: style == %") % style));
166 }
167 } // namespace
168 
169 class InterpolateGsl : public Interpolate {
171  std::vector<double> const &y,
172  Interpolate::Style const style);
173 
174 public:
175  ~InterpolateGsl() override;
176  double interpolate(double const x) const override;
177 
178 private:
180  Interpolate::Style const style);
181 
182  ::gsl_interp_type const *_interpType;
183  ::gsl_interp_accel *_acc;
184  ::gsl_interp *_interp;
185 };
186 
187 InterpolateGsl::InterpolateGsl(std::vector<double> const &x,
188  std::vector<double> const &y,
189  Interpolate::Style const style
190  )
191  : Interpolate(x, y), _interpType(styleToGslInterpType(style)) {
192  // Turn the gsl error handler off, we want to use our own exceptions
193  ::gsl_set_error_handler_off();
194 
195  _acc = ::gsl_interp_accel_alloc();
196  if (!_acc) {
197  throw std::bad_alloc();
198  }
199 
200  _interp = ::gsl_interp_alloc(_interpType, _y.size());
201  if (!_interp) {
203  str(boost::format("Failed to initialise spline for type %s, length %d") %
204  _interpType->name % _y.size()));
205  }
206  // Note, "x" and "y" are vector<double>; gsl_inter_init requires double[].
207  // The &(x[0]) here is valid because std::vector guarantees that the values are
208  // stored contiguously in memory (for types other than bool); C++0X 23.3.6.1 for
209  // those of you reading along.
210  int const status = ::gsl_interp_init(_interp, &x[0], &y[0], _y.size());
211  if (status != 0) {
212  throw LSST_EXCEPT(
214  str(boost::format("gsl_interp_init failed: %s [%d]") % ::gsl_strerror(status) % status));
215  }
216 }
217 
219  ::gsl_interp_free(_interp);
220  ::gsl_interp_accel_free(_acc);
221 }
222 
223 double InterpolateGsl::interpolate(double const xInterp) const {
224  // New GSL versions refuse to extrapolate.
225  // gsl_interp_init() requires x to be ordered, so can just check
226  // the array endpoints for out-of-bounds.
227  if ((xInterp < _x.front() || (xInterp > _x.back()))) {
228  // do our own quadratic extrapolation.
229  // (GSL only provides first and second derivative functions)
230  /* could also just fail via:
231  throw LSST_EXCEPT(
232  pex::exceptions::InvalidParameterError,
233  (boost::format("Interpolation point %f outside range [%f, %f]")
234  % x % _x.front() % _x.back()).str()
235  );
236  */
237  double x0, y0;
238  if (xInterp < _x.front()) {
239  x0 = _x.front();
240  y0 = _y.front();
241  } else {
242  x0 = _x.back();
243  y0 = _y.back();
244  }
245  // first derivative at endpoint
246  double d = ::gsl_interp_eval_deriv(_interp, &_x[0], &_y[0], x0, _acc);
247  // second derivative at endpoint
248  double d2 = ::gsl_interp_eval_deriv2(_interp, &_x[0], &_y[0], x0, _acc);
249  return y0 + (xInterp - x0) * d + 0.5 * (xInterp - x0) * (xInterp - x0) * d2;
250  }
251  assert(xInterp >= _x.front());
252  assert(xInterp <= _x.back());
253  return ::gsl_interp_eval(_interp, &_x[0], &_y[0], xInterp, _acc);
254 }
255 
257  static std::map<std::string, Interpolate::Style> gslInterpTypeStrings;
258  if (gslInterpTypeStrings.empty()) {
259  gslInterpTypeStrings["CONSTANT"] = Interpolate::CONSTANT;
260  gslInterpTypeStrings["LINEAR"] = Interpolate::LINEAR;
261  gslInterpTypeStrings["CUBIC_SPLINE"] = Interpolate::CUBIC_SPLINE;
262  gslInterpTypeStrings["NATURAL_SPLINE"] = Interpolate::NATURAL_SPLINE;
263  gslInterpTypeStrings["CUBIC_SPLINE_PERIODIC"] = Interpolate::CUBIC_SPLINE_PERIODIC;
264  gslInterpTypeStrings["AKIMA_SPLINE"] = Interpolate::AKIMA_SPLINE;
265  gslInterpTypeStrings["AKIMA_SPLINE_PERIODIC"] = Interpolate::AKIMA_SPLINE_PERIODIC;
266  }
267 
268  if (gslInterpTypeStrings.find(style) == gslInterpTypeStrings.end()) {
269  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Interp style not found: " + style);
270  }
271  return gslInterpTypeStrings[style];
272 }
273 
275  if (n < 1) {
276  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "n must be greater than 0");
277  } else if (n > 4) {
279  } else {
280  static std::vector<Interpolate::Style> styles;
281  if (styles.empty()) {
282  styles.resize(5);
283 
284  styles[0] = Interpolate::UNKNOWN; // impossible to reach as we check for n < 1
285  styles[1] = Interpolate::CONSTANT;
286  styles[2] = Interpolate::LINEAR;
287  styles[3] = Interpolate::CUBIC_SPLINE;
288  styles[4] = Interpolate::CUBIC_SPLINE;
289  }
290  return styles[n];
291  }
292 }
293 
295  size_t const num = x.size();
296  std::vector<double> out(num);
297  for (size_t i = 0; i < num; ++i) {
298  out[i] = interpolate(x[i]);
299  }
300  return out;
301 }
302 
303 ndarray::Array<double, 1> Interpolate::interpolate(ndarray::Array<double const, 1> const &x) const {
304  int const num = x.getShape()[0];
305  ndarray::Array<double, 1> out = ndarray::allocate(ndarray::makeVector(num));
306  for (int i = 0; i < num; ++i) {
307  std::cout << "Interpolating " << x[i] << std::endl;
308  out[i] = interpolate(x[i]);
309  }
310  return out;
311 }
312 
314  static std::vector<int> minPoints;
315  if (minPoints.empty()) {
316  minPoints.resize(Interpolate::NUM_STYLES);
317  minPoints[Interpolate::CONSTANT] = 1;
318  minPoints[Interpolate::LINEAR] = 2;
319  minPoints[Interpolate::NATURAL_SPLINE] = 3;
320  minPoints[Interpolate::CUBIC_SPLINE] = 3;
321  minPoints[Interpolate::CUBIC_SPLINE_PERIODIC] = 3;
322  minPoints[Interpolate::AKIMA_SPLINE] = 5;
323  minPoints[Interpolate::AKIMA_SPLINE_PERIODIC] = 5;
324  }
325 
326  if (style >= 0 && style < Interpolate::NUM_STYLES) {
327  return minPoints[style];
328  } else {
329  throw LSST_EXCEPT(
331  str(boost::format("Style %d is out of range 0..%d") % style % (Interpolate::NUM_STYLES - 1)));
332  }
333 }
334 
336 
337  Interpolate::Style const style)
338  : _x(xy.first), _y(xy.second), _style(style) {
339  ;
340 }
341 
343  Interpolate::Style const style) {
344  switch (style) {
347  default: // use GSL
348  return std::shared_ptr<Interpolate>(new InterpolateGsl(x, y, style));
349  }
350 }
351 
352 std::shared_ptr<Interpolate> makeInterpolate(ndarray::Array<double const, 1> const &x,
353  ndarray::Array<double const, 1> const &y,
354  Interpolate::Style const style) {
355  return makeInterpolate(std::vector<double>(x.begin(), x.end()), std::vector<double>(y.begin(), y.end()),
356  style);
357 }
358 } // namespace math
359 } // namespace afw
360 } // namespace lsst
double x
#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.
Definition: Interpolate.cc:342
double interpolate(double const x) const override
Definition: Interpolate.cc:94
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.
Definition: Interpolate.cc:342
double interpolate(double const x) const override
Definition: Interpolate.cc:223
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 make_pair(T... args)
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for 'n' points.
Definition: Interpolate.cc:274
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style.
Definition: Interpolate.cc:313
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:256
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.
Definition: Interpolate.cc:342
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T resize(T... args)
T size(T... args)
T upper_bound(T... args)