Alexandria  2.22.0
Please provide a description of the project.
spline.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 CubicInterpolator final : public PiecewiseBase {
34 public:
37  : PiecewiseBase(std::move(knots))
38  , m_coef0(std::move(coef0))
39  , m_coef1(std::move(coef1))
40  , m_coef2(std::move(coef2))
41  , m_coef3(std::move(coef3)) {}
42 
43  virtual ~CubicInterpolator() = default;
44 
45  double operator()(double x) const override {
46  auto i = findKnot(x);
47  if (i < 0 || i >= static_cast<ssize_t>(m_knots.size())) {
48  return 0.;
49  }
50  if (i == 0) {
51  return m_coef3.front() * x * x * x + m_coef2.front() * x * x + m_coef1.front() * x + m_coef0.front();
52  }
53  --i;
54  return m_coef3[i] * x * x * x + m_coef2[i] * x * x + m_coef1[i] * x + m_coef0[i];
55  }
56 
57  void operator()(const std::vector<double>& xs, std::vector<double>& out) const override {
58  out.resize(xs.size());
59  // Find the first X that is within the range
60  auto first_x = std::lower_bound(xs.begin(), xs.end(), m_knots.front());
61  auto n_less = first_x - xs.begin();
62 
63  // Opening 0s
64  auto o = out.begin() + n_less;
65  std::fill(out.begin(), o, 0.);
66 
67  // To avoid if within the loop, treat values exactly equal to the first knot here
68  auto x = xs.begin() + n_less;
69  while (*x == m_knots.front()) {
70  auto x2 = *x * *x;
71  *o = m_coef3.front() * x2 * *x + m_coef2.front() * x2 + m_coef1.front() * *x + m_coef0.front();
72  ++o, ++x;
73  }
74 
75  // Interpolate values within range
76  auto current_knot = std::lower_bound(m_knots.begin(), m_knots.end(), *x);
77  while (o != out.end() && current_knot < m_knots.end()) {
78  auto i = current_knot - m_knots.begin();
79  auto x2 = *x * *x;
80  --i;
81  *o = m_coef3[i] * x2 * *x + m_coef2[i] * x2 + m_coef1[i] * *x + m_coef0[i];
82  ++o, ++x;
83  current_knot = std::find_if(current_knot, m_knots.end(), [x](double k) { return k >= *x; });
84  }
85 
86  // Trailing 0s
87  std::fill(o, out.end(), 0.);
88  }
89 
92  }
93 
94  double integrate(const double a, const double b) const override {
95  if (a == b) {
96  return 0;
97  }
98  int direction = 1;
99  double min = a;
100  double max = b;
101  if (min > max) {
102  direction = -1;
103  min = b;
104  max = a;
105  }
106  double result = 0;
107  auto knotIter = std::upper_bound(m_knots.begin(), m_knots.end(), min);
108  if (knotIter != m_knots.begin()) {
109  --knotIter;
110  }
111  auto i = knotIter - m_knots.begin();
112  while (++knotIter != m_knots.end()) {
113  auto prevKnotIter = knotIter - 1;
114  if (max <= *prevKnotIter) {
115  break;
116  }
117  if (min < *knotIter) {
118  double down = (min > *prevKnotIter) ? min : *prevKnotIter;
119  double up = (max < *knotIter) ? max : *knotIter;
120  result += antiderivative(i, up) - antiderivative(i, down);
121  }
122  ++i;
123  }
124  return direction * result;
125  }
126 
127 private:
129 
130  double antiderivative(int i, double x) const {
131  double x2 = x * x;
132  return m_coef0[i] * x + m_coef1[i] * x2 / 2. + m_coef2[i] * x2 * x / 3. + m_coef3[i] * x2 * x2 / 4.;
133  }
134 };
135 
137  bool extrapolate) {
138  std::vector<double> knots(x);
139 
140  // Number of intervals
141  int n = x.size() - 1;
142 
143  // Differences between knot points
144  std::vector<double> h(n, 0.);
145  for (int i = 0; i < n; i++)
146  h[i] = x[i + 1] - x[i];
147 
148  std::vector<double> mu(n, 0.);
149  std::vector<double> z(n + 1, 0.);
150  for (int i = 1; i < n; ++i) {
151  double g = 2. * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
152  mu[i] = h[i] / g;
153  z[i] = (3. * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1]) + y[i - 1] * h[i]) / (h[i - 1] * h[i]) -
154  h[i - 1] * z[i - 1]) /
155  g;
156  }
157 
158  // cubic spline coefficients
159  std::vector<double> coef0(n, 0.);
160  std::vector<double> coef1(n, 0.);
161  std::vector<double> coef2(n + 1, 0.);
162  std::vector<double> coef3(n, 0.);
163 
164  z[n] = 0.;
165  coef2[n] = 0.;
166 
167  for (int j = n - 1; j >= 0; j--) {
168  coef0[j] = y[j];
169  coef2[j] = z[j] - mu[j] * coef2[j + 1];
170  coef1[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (coef2[j + 1] + 2. * coef2[j]) / 3.;
171  coef3[j] = (coef2[j + 1] - coef2[j]) / (3. * h[j]);
172  }
173 
174  // The above were taken from SplineInterpolator from Apache commons math. These
175  // polynomials need to be shifted by -x[i] in our case.
176  for (int i = 0; i < n; i++) {
177  double x_1 = -x[i];
178  double x_2 = x_1 * x_1;
179  double x_3 = x_1 * x_2;
180  coef0[i] = coef0[i] + coef1[i] * x_1 + coef2[i] * x_2 + coef3[i] * x_3;
181  coef1[i] = coef1[i] + 2. * coef2[i] * x_1 + 3. * coef3[i] * x_2;
182  coef2[i] = coef2[i] + 3. * coef3[i] * x_1;
183  // d[i] keeps the same value
184  }
185 
186  if (extrapolate) {
189  }
190 
192  new CubicInterpolator(std::move(knots), std::move(coef0), std::move(coef1), std::move(coef2), std::move(coef3)));
193 }
194 
195 } // namespace MathUtils
196 } // end of namespace Euclid
T back(T... args)
T begin(T... args)
std::vector< double > m_coef2
Definition: spline.cpp:128
std::vector< double > m_coef0
Definition: spline.cpp:128
CubicInterpolator(std::vector< double > knots, std::vector< double > coef0, std::vector< double > coef1, std::vector< double > coef2, std::vector< double > coef3)
Definition: spline.cpp:35
double integrate(const double a, const double b) const override
Definition: spline.cpp:94
std::vector< double > m_coef3
Definition: spline.cpp:128
std::unique_ptr< NAryFunction > clone() const override
Definition: spline.cpp:90
std::vector< double > m_coef1
Definition: spline.cpp:128
void operator()(const std::vector< double > &xs, std::vector< double > &out) const override
Definition: spline.cpp:57
double antiderivative(int i, double x) const
Definition: spline.cpp:130
double operator()(double x) const override
Definition: spline.cpp:45
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)
constexpr double g
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
STL namespace.
T resize(T... args)
T size(T... args)
T upper_bound(T... args)