Alexandria  2.22.0
Please provide a description of the project.
interpolation.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2021 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]) {
69  if (y[i] == y[i - 1]) {
70  continue;
71  }
72  }
73  if (x[i] < x[i - 1]) {
74  throw InterpolationException() << "Only increasing X values are supported "
75  << "but found " << x[i] << " after " << x[i - 1];
76  }
77  final_x.emplace_back(x[i]);
78  final_y.emplace_back(y[i]);
79  }
80 
81  switch (type) {
83  return linearInterpolation(final_x, final_y, extrapolate);
85  return splineInterpolation(final_x, final_y, extrapolate);
86  }
87  return nullptr;
88 }
89 
91  bool extrapolate) {
94  x.reserve(dataset.size());
95  y.reserve(dataset.size());
96  for (auto& pair : dataset) {
97  x.emplace_back(pair.first);
98  y.emplace_back(pair.second);
99  }
100  return interpolate(x, y, type, extrapolate);
101 }
102 
103 double simple_interpolation(double x, const std::vector<double>& xp, const std::vector<double>& yp, bool extrapolate) {
104  if (xp.empty()) {
105  throw InterpolationException() << "Can not interpolate an empty list of values";
106  }
107  if (xp.size() != yp.size()) {
108  throw InterpolationException() << "Number of knots and number of values do not match";
109  }
110  if (xp.size() == 1) {
111  return (extrapolate || xp.front() == x) ? yp.front() : 0.;
112  }
113  auto gt = std::upper_bound(xp.begin(), xp.end(), x);
114  std::size_t i = gt - xp.begin();
115  if (i == 0) {
116  if (!extrapolate && x < xp.front())
117  return 0.;
118  ++i;
119  }
120  if (i == xp.size()) {
121  if (!extrapolate && x > xp.back())
122  return 0.;
123  --i;
124  }
125  double x1 = xp[i], x0 = xp[i - 1];
126  double y1 = yp[i], y0 = yp[i - 1];
127 
128  if (x1 == x0) {
129  return (x1 + y1) / 2.;
130  }
131 
132  double coef1 = (y1 - y0) / (x1 - x0);
133  double coef0 = y0 - coef1 * x0;
134  return coef0 + coef1 * x;
135 }
136 
137 } // namespace MathUtils
138 } // 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
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:136
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:126
T reserve(T... args)
T size(T... args)
T upper_bound(T... args)