LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
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 "gsl/gsl_spline.h"
36 #include "ndarray.h"
37 #include "lsst/pex/exceptions.h"
39 
40 namespace lsst {
41 namespace afw {
42 namespace math {
43 
44 namespace {
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 
78  std::vector<double> const &y,
79  Interpolate::Style const style);
80 
81 public:
82  ~InterpolateConstant() override {}
83  double interpolate(double const x) const override;
84 
85 private:
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 
95 double 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
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 
136 namespace {
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.");
145  case Interpolate::LINEAR:
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 
170 class InterpolateGsl : public Interpolate {
172  std::vector<double> const &y,
173  Interpolate::Style const style);
174 
175 public:
176  ~InterpolateGsl() override;
177  double interpolate(double const x) const override;
178 
179 private:
181  Interpolate::Style const style);
182 
183  ::gsl_interp_type const *_interpType;
184  ::gsl_interp_accel *_acc;
185  ::gsl_interp *_interp;
186 };
187 
188 InterpolateGsl::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 
224 double InterpolateGsl::interpolate(double const xInterp) const {
225  // New GSL versions refuse to extrapolate.
226  // gsl_interp_init() requires x to be ordered, so can just check
227  // the array endpoints for out-of-bounds.
228  if ((xInterp < _x.front() || (xInterp > _x.back()))) {
229  // do our own quadratic extrapolation.
230  // (GSL only provides first and second derivative functions)
231  /* could also just fail via:
232  throw LSST_EXCEPT(
233  pex::exceptions::InvalidParameterError,
234  (boost::format("Interpolation point %f outside range [%f, %f]")
235  % x % _x.front() % _x.back()).str()
236  );
237  */
238  double x0, y0;
239  if (xInterp < _x.front()) {
240  x0 = _x.front();
241  y0 = _y.front();
242  } else {
243  x0 = _x.back();
244  y0 = _y.back();
245  }
246  // first derivative at endpoint
247  double d = ::gsl_interp_eval_deriv(_interp, &_x[0], &_y[0], x0, _acc);
248  // second derivative at endpoint
249  double d2 = ::gsl_interp_eval_deriv2(_interp, &_x[0], &_y[0], x0, _acc);
250  return y0 + (xInterp - x0) * d + 0.5 * (xInterp - x0) * (xInterp - x0) * d2;
251  }
252  assert(xInterp >= _x.front());
253  assert(xInterp <= _x.back());
254  return ::gsl_interp_eval(_interp, &_x[0], &_y[0], xInterp, _acc);
255 }
256 
258  static std::map<std::string, Interpolate::Style> gslInterpTypeStrings;
259  if (gslInterpTypeStrings.empty()) {
260  gslInterpTypeStrings["CONSTANT"] = Interpolate::CONSTANT;
261  gslInterpTypeStrings["LINEAR"] = Interpolate::LINEAR;
262  gslInterpTypeStrings["CUBIC_SPLINE"] = Interpolate::CUBIC_SPLINE;
263  gslInterpTypeStrings["NATURAL_SPLINE"] = Interpolate::NATURAL_SPLINE;
264  gslInterpTypeStrings["CUBIC_SPLINE_PERIODIC"] = Interpolate::CUBIC_SPLINE_PERIODIC;
265  gslInterpTypeStrings["AKIMA_SPLINE"] = Interpolate::AKIMA_SPLINE;
266  gslInterpTypeStrings["AKIMA_SPLINE_PERIODIC"] = Interpolate::AKIMA_SPLINE_PERIODIC;
267  }
268 
269  if (gslInterpTypeStrings.find(style) == gslInterpTypeStrings.end()) {
270  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Interp style not found: " + style);
271  }
272  return gslInterpTypeStrings[style];
273 }
274 
276  if (n < 1) {
277  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "n must be greater than 0");
278  } else if (n > 4) {
280  } else {
281  static std::vector<Interpolate::Style> styles;
282  if (styles.empty()) {
283  styles.resize(5);
284 
285  styles[0] = Interpolate::UNKNOWN; // impossible to reach as we check for n < 1
286  styles[1] = Interpolate::CONSTANT;
287  styles[2] = Interpolate::LINEAR;
288  styles[3] = Interpolate::CUBIC_SPLINE;
289  styles[4] = Interpolate::CUBIC_SPLINE;
290  }
291  return styles[n];
292  }
293 }
294 
296  size_t const num = x.size();
297  std::vector<double> out(num);
298  for (size_t i = 0; i < num; ++i) {
299  out[i] = interpolate(x[i]);
300  }
301  return out;
302 }
303 
304 ndarray::Array<double, 1> Interpolate::interpolate(ndarray::Array<double const, 1> const &x) const {
305  int const num = x.getShape()[0];
306  ndarray::Array<double, 1> out = ndarray::allocate(ndarray::makeVector(num));
307  for (int i = 0; i < num; ++i) {
308  std::cout << "Interpolating " << x[i] << std::endl;
309  out[i] = interpolate(x[i]);
310  }
311  return out;
312 }
313 
315  static std::vector<int> minPoints;
316  if (minPoints.empty()) {
317  minPoints.resize(Interpolate::NUM_STYLES);
318  minPoints[Interpolate::CONSTANT] = 1;
319  minPoints[Interpolate::LINEAR] = 2;
320  minPoints[Interpolate::NATURAL_SPLINE] = 3;
321  minPoints[Interpolate::CUBIC_SPLINE] = 3;
322  minPoints[Interpolate::CUBIC_SPLINE_PERIODIC] = 3;
323  minPoints[Interpolate::AKIMA_SPLINE] = 5;
324  minPoints[Interpolate::AKIMA_SPLINE_PERIODIC] = 5;
325  }
326 
327  if (style >= 0 && style < Interpolate::NUM_STYLES) {
328  return minPoints[style];
329  } else {
330  throw LSST_EXCEPT(
332  str(boost::format("Style %d is out of range 0..%d") % style % (Interpolate::NUM_STYLES - 1)));
333  }
334 }
335 
337 
338  Interpolate::Style const style)
339  : _x(xy.first), _y(xy.second), _style(style) {
340  ;
341 }
342 
344  Interpolate::Style const style) {
345  switch (style) {
347  return std::shared_ptr<Interpolate>(new InterpolateConstant(x, y, style));
348  default: // use GSL
349  return std::shared_ptr<Interpolate>(new InterpolateGsl(x, y, style));
350  }
351 }
352 
353 std::shared_ptr<Interpolate> makeInterpolate(ndarray::Array<double const, 1> const &x,
354  ndarray::Array<double const, 1> const &y,
355  Interpolate::Style const style) {
356  return makeInterpolate(std::vector<double>(x.begin(), x.end()), std::vector<double>(y.begin(), y.end()),
357  style);
358 }
359 } // namespace math
360 } // namespace afw
361 } // namespace lsst
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
Interpolate::Style stringToInterpStyle(std::string const &style)
Conversion function to switch a string to an Interpolate::Style.
Definition: Interpolate.cc:257
T empty(T... args)
Interpolate(Interpolate const &)=delete
double interpolate(double const x) const override
Definition: Interpolate.cc:95
T front(T... args)
virtual double interpolate(double const x) const =0
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:343
T endl(T... args)
T upper_bound(T... args)
int y
Definition: SpanSet.cc:49
T end(T... args)
STL class.
T resize(T... args)
STL class.
A base class for image defects.
double interpolate(double const x) const override
Definition: Interpolate.cc:224
T make_pair(T... args)
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
std::vector< double > const _x
Definition: Interpolate.h:83
double x
Interpolate::Style lookupMaxInterpStyle(int const n)
Get the highest order Interpolation::Style available for &#39;n&#39; points.
Definition: Interpolate.cc:275
T find(T... args)
T size(T... args)
int lookupMinInterpPoints(Interpolate::Style const style)
Get the minimum number of points needed to use the requested interpolation style. ...
Definition: Interpolate.cc:314
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:343
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
T begin(T... args)
STL class.
std::vector< double > const _y
Definition: Interpolate.h:84
Interpolate::Style const _style
Definition: Interpolate.h:85
T back(T... args)
Reports invalid arguments.
Definition: Runtime.h:66
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104