Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
linear.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 <limits>
29
30namespace Euclid {
31namespace MathUtils {
32
34public:
37
38 virtual ~LinearInterpolator() = default;
39
40 double operator()(double x) const override {
41 auto i = findKnot(x);
42 if (i < 0 || i >= static_cast<ssize_t>(m_knots.size())) {
43 return 0;
44 }
45 if (i == 0) {
46 return m_coef1.front() * x + m_coef0.front();
47 }
48 return m_coef1[i - 1] * x + m_coef0[i - 1];
49 }
50
51 void operator()(const std::vector<double>& xs, std::vector<double>& out) const override {
52 out.resize(xs.size());
53 // Find the first X that is within the range
55 auto n_less = first_x - xs.begin();
56
57 // Opening 0s
58 auto o = out.begin() + n_less;
59 std::fill(out.begin(), o, 0.);
60
61 // To avoid if within the loop, treat values exactly equal to the first knot here
62 auto x = xs.begin() + n_less;
63 while (x < xs.end() && *x == m_knots.front()) {
64 *o = m_coef1.front() * *x + m_coef0.front();
65 ++o, ++x;
66 }
67 if (x == xs.end()) {
68 return;
69 }
70
71 // Interpolate values within range
73 while (o != out.end() && current_knot < m_knots.end()) {
75 auto i = current_knot - m_knots.begin();
76 *o = m_coef1[i - 1] * *x + m_coef0[i - 1];
77 ++o, ++x;
78 }
79
80 // Trailing 0s
81 std::fill(o, out.end(), 0.);
82 }
83
87
88 double integrate(const double a, const double b) const override {
89 if (a == b) {
90 return 0;
91 }
92 int direction = 1;
93 double min = a;
94 double max = b;
95 if (min > max) {
96 direction = -1;
97 min = b;
98 max = a;
99 }
100 double result = 0;
102 if (knotIter != m_knots.begin()) {
103 --knotIter;
104 }
105 auto i = knotIter - m_knots.begin();
106 while (++knotIter != m_knots.end()) {
107 auto prevKnotIter = knotIter - 1;
108 if (max <= *prevKnotIter) {
109 break;
110 }
111 if (min < *knotIter) {
112 double down = (min > *prevKnotIter) ? min : *prevKnotIter;
113 double up = (max < *knotIter) ? max : *knotIter;
115 }
116 ++i;
117 }
118 return direction * result;
119 }
120
121private:
123
124 double antiderivative(int i, double x) const {
125 return m_coef0[i] * x + m_coef1[i] * x * x / 2;
126 }
127};
128
130 bool extrapolate) {
132
133 for (size_t i = 0; i < x.size() - 1; i++) {
134 double dx = (x[i + 1] - x[i]);
135 if (dx > 0) {
136 coef1[i] = (y[i + 1] - y[i]) / dx;
137 coef0[i] = y[i] - coef1[i] * x[i];
138 }
139 }
140
141 if (extrapolate) {
144 }
145
147}
148
149} // namespace MathUtils
150} // end of namespace Euclid
T back(T... args)
T begin(T... args)
LinearInterpolator(std::vector< double > knots, std::vector< double > coef0, std::vector< double > coef1)
Definition linear.cpp:35
const std::vector< double > m_coef0
Definition linear.cpp:122
std::unique_ptr< NAryFunction > clone() const override
Definition linear.cpp:84
void operator()(const std::vector< double > &xs, std::vector< double > &out) const override
Definition linear.cpp:51
double integrate(const double a, const double b) const override
Definition linear.cpp:88
double antiderivative(int i, double x) const
Definition linear.cpp:124
const std::vector< double > m_coef1
Definition linear.cpp:122
double operator()(double x) const override
Definition linear.cpp:40
Represents a piecewise function.
Definition Piecewise.h:48
ssize_t findKnot(double x) const
Definition Piecewise.h:60
std::vector< double > m_knots
A vector where the knots are kept.
Definition Piecewise.h:74
T end(T... args)
T fill(T... args)
T front(T... args)
T lower_bound(T... args)
T lowest(T... args)
T max(T... args)
T move(T... args)
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
STL namespace.
T resize(T... args)
T size(T... args)
T upper_bound(T... args)