Alexandria  2.27.0
SDC-CH common library for the Euclid project
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 
30 namespace Euclid {
31 namespace MathUtils {
32 
33 class LinearInterpolator final : public PiecewiseBase {
34 public:
36  : PiecewiseBase(std::move(knots)), m_coef0(std::move(coef0)), m_coef1(std::move(coef1)){};
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
54  auto first_x = std::lower_bound(xs.begin(), xs.end(), m_knots.front());
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
72  auto current_knot = std::lower_bound(m_knots.begin(), m_knots.end(), *x);
73  while (o != out.end() && current_knot < m_knots.end()) {
74  current_knot = std::lower_bound(current_knot, m_knots.end(), *x);
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 
86  }
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;
101  auto knotIter = std::upper_bound(m_knots.begin(), m_knots.end(), min);
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;
114  result += antiderivative(i, up) - antiderivative(i, down);
115  }
116  ++i;
117  }
118  return direction * result;
119  }
120 
121 private:
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) {
131  std::vector<double> coef0(x.size()), coef1(x.size()), knots(x);
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 
146  return std::unique_ptr<Function>(new LinearInterpolator(std::move(knots), std::move(coef0), std::move(coef1)));
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
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
std::unique_ptr< NAryFunction > clone() const override
Definition: linear.cpp:84
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::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)