Alexandria  2.22.0
Please provide a description of the project.
linear.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 
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 == m_knots.front()) {
64  *o = m_coef1.front() * *x + m_coef0.front();
65  ++o, ++x;
66  }
67 
68  // Interpolate values within range
69  auto current_knot = std::lower_bound(m_knots.begin(), m_knots.end(), *x);
70  while (o != out.end() && current_knot < m_knots.end()) {
71  auto i = current_knot - m_knots.begin();
72  *o = m_coef1[i - 1] * *x + m_coef0[i - 1];
73  ++o, ++x;
74  current_knot = std::find_if(current_knot, m_knots.end(), [x](double k) { return k >= *x; });
75  }
76 
77  // Trailing 0s
78  std::fill(o, out.end(), 0.);
79  }
80 
83  }
84 
85  double integrate(const double a, const double b) const override {
86  if (a == b) {
87  return 0;
88  }
89  int direction = 1;
90  double min = a;
91  double max = b;
92  if (min > max) {
93  direction = -1;
94  min = b;
95  max = a;
96  }
97  double result = 0;
98  auto knotIter = std::upper_bound(m_knots.begin(), m_knots.end(), min);
99  if (knotIter != m_knots.begin()) {
100  --knotIter;
101  }
102  auto i = knotIter - m_knots.begin();
103  while (++knotIter != m_knots.end()) {
104  auto prevKnotIter = knotIter - 1;
105  if (max <= *prevKnotIter) {
106  break;
107  }
108  if (min < *knotIter) {
109  double down = (min > *prevKnotIter) ? min : *prevKnotIter;
110  double up = (max < *knotIter) ? max : *knotIter;
111  result += antiderivative(i, up) - antiderivative(i, down);
112  }
113  ++i;
114  }
115  return direction * result;
116  }
117 
118 private:
120 
121  double antiderivative(int i, double x) const {
122  return m_coef0[i] * x + m_coef1[i] * x * x / 2;
123  }
124 };
125 
127  bool extrapolate) {
128  std::vector<double> coef0(x.size()), coef1(x.size()), knots(x);
129 
130  for (size_t i = 0; i < x.size() - 1; i++) {
131  coef1[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
132  coef0[i] = y[i] - coef1[i] * x[i];
133  }
134 
135  if (extrapolate) {
138  }
139 
140  return std::unique_ptr<Function>(new LinearInterpolator(std::move(knots), std::move(coef0), std::move(coef1)));
141 }
142 
143 } // namespace MathUtils
144 } // 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:119
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:85
double antiderivative(int i, double x) const
Definition: linear.cpp:121
const std::vector< double > m_coef1
Definition: linear.cpp:119
double operator()(double x) const override
Definition: linear.cpp:40
std::unique_ptr< NAryFunction > clone() const override
Definition: linear.cpp:81
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 find_if(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:126
STL namespace.
T resize(T... args)
T size(T... args)
T upper_bound(T... args)