Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
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
28#include "implementations.h"
30
31namespace Euclid {
32namespace MathUtils {
33
34namespace {
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
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) {
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) {
121 }
122 return nullptr;
123}
124
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
153double 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
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.
std::array< std::vector< double >, N > Coordinates
Used to pass the grid coordinates to interpn. Internally will make a copy of the required values.
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 size(T... args)
T upper_bound(T... args)