Alexandria  2.27.0
SDC-CH common library for the Euclid project
interpolation.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2022 Euclid Science Ground Segment
3  *
4  * This library is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the Free
6  * Software Foundation; either version 3.0 of the License, or (at your option)
7  * any later version.
8  *
9  * This library is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU Lesser General Public License
15  * along with this library; if not, write to the Free Software Foundation, Inc.,
16  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
26 #include "ElementsKernel/Logging.h"
28 #include "implementations.h"
30 
31 namespace Euclid {
32 namespace MathUtils {
33 
34 namespace {
35 
37 
38 } // Anonymous namespace
39 
41  InterpolationType type, bool extrapolate) {
42 
43  if (x.size() != y.size()) {
44  throw InterpolationException() << "Interpolation using vectors of incompatible "
45  << "size: X=" << x.size() << ", Y=" << y.size();
46  }
47 
48  if (x.size() == 1) {
49  auto c = y.front();
50  if (extrapolate) {
51  return make_unique<FunctionAdapter>([c](double) { return c; });
52  }
53  auto sx = x.front();
54  return make_unique<FunctionAdapter>([c, sx](double v) { return c * (v == sx); });
55  }
56 
57  // We remove any duplicate lines and we check that we have only increasing
58  // X values and no step functions
59  std::vector<double> final_x;
60  std::vector<double> final_y;
61 
62  final_x.reserve(x.size());
63  final_y.reserve(x.size());
64 
65  final_x.emplace_back(x[0]);
66  final_y.emplace_back(y[0]);
67  for (std::size_t i = 1; i < x.size(); ++i) {
68  if (x[i] == x[i - 1] && y[i] == y[i - 1]) {
69  continue;
70  }
71  if (x[i] < x[i - 1]) {
72  throw InterpolationException() << "Only increasing X values are supported "
73  << "but found " << x[i] << " after " << x[i - 1];
74  }
75  final_x.emplace_back(x[i]);
76  final_y.emplace_back(y[i]);
77  }
78 
79  switch (type) {
81  return linearInterpolation(final_x, final_y, extrapolate);
83  return splineInterpolation(final_x, final_y, extrapolate);
84  }
85  return nullptr;
86 }
87 
89  bool extrapolate) {
90  if (dataset.size() == 1) {
91  auto c = dataset.front();
92  if (extrapolate) {
93  return make_unique<FunctionAdapter>([c](double) { return c.second; });
94  }
95  return make_unique<FunctionAdapter>([c](double v) { return c.second * (v == c.first); });
96  }
97 
100  x.reserve(dataset.size());
101  y.reserve(dataset.size());
102  for (auto& pair : dataset) {
103  if (!x.empty()) {
104  if (pair.first == x.back() && pair.second == y.back()) {
105  continue;
106  }
107  if (pair.first < x.back()) {
108  throw InterpolationException() << "Only increasing X values are supported "
109  << "but found " << pair.first << " after " << x.back();
110  }
111  }
112  x.emplace_back(pair.first);
113  y.emplace_back(pair.second);
114  }
115 
116  switch (type) {
118  return linearInterpolation(x, y, extrapolate);
120  return splineInterpolation(x, y, extrapolate);
121  }
122  return nullptr;
123 }
124 
125 double simple_interpolation(double x, const std::vector<double>& xp, const std::vector<double>& yp, bool extrapolate) {
126  if (xp.empty()) {
127  throw InterpolationException() << "Can not interpolate an empty list of values";
128  }
129  if (xp.size() != yp.size()) {
130  throw InterpolationException() << "Number of knots and number of values do not match";
131  }
132  if (xp.size() == 1) {
133  return (extrapolate || xp.front() == x) ? yp.front() : 0.;
134  }
135 
136  if ((x < xp.front() || x > xp.back()) && !extrapolate) {
137  return 0;
138  }
139 
140  auto gt = std::upper_bound(xp.begin(), xp.end(), x);
141  std::size_t i = gt - xp.begin();
142  if (i == 0) {
143  ++i;
144  }
145  if (i == xp.size()) {
146  if (!extrapolate && x > xp.back())
147  return 0.;
148  --i;
149  }
150  return simple_interpolation(x, xp[i - 1], xp[i], yp[i - 1], yp[i], extrapolate);
151 }
152 
153 double simple_interpolation(double x, double x0, double x1, double y0, double y1, bool extrapolate) noexcept {
154  // Looks convolved, but this avoids branching
155  double within = ((x >= x0) & (x <= x1)) | extrapolate;
156  double coef1 = (y1 - y0) / (x1 - x0);
157  double coef0 = y0 - coef1 * x0;
158  return within * (coef0 + coef1 * x);
159 }
160 
161 } // namespace MathUtils
162 } // end of namespace Euclid
static Elements::Logging logger
Logger.
Definition: Example.cpp:55
T back(T... args)
T begin(T... args)
static Logging getLogger(const std::string &name="")
This module provides an interface for accessing two dimensional datasets (pairs of (X,...
Definition: XYDataset.h:59
const std::pair< double, double > & front() const
Returns a reference to the first pair of the dataset.
Definition: XYDataset.cpp:44
size_t size() const
Get the size of the vector container.
Definition: XYDataset.h:150
T emplace_back(T... args)
T empty(T... args)
T end(T... args)
T front(T... args)
ELEMENTS_API std::unique_ptr< Function > interpolate(const std::vector< double > &x, const std::vector< double > &y, InterpolationType type, bool extrapolate=false)
ELEMENTS_API double simple_interpolation(double x, const std::vector< double > &xp, const std::vector< double > &yp, bool extrapolate=false)
std::unique_ptr< Function > splineInterpolation(const std::vector< double > &x, const std::vector< double > &y, bool extrapolate)
Performs cubic spline interpolation for the given set of data points.
Definition: spline.cpp:139
InterpolationType
Enumeration of the different supported interpolation types.
Definition: interpolation.h:43
std::unique_ptr< Function > linearInterpolation(const std::vector< double > &x, const std::vector< double > &y, bool extrapolate)
Performs linear interpolation for the given set of data points.
Definition: linear.cpp:129
T reserve(T... args)
T size(T... args)
T upper_bound(T... args)