Alexandria  2.19
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 
29 #include <limits>
30 
31 namespace Euclid {
32 namespace MathUtils {
33 
35 
36  // Number of intervals
37  int n = x.size() - 1;
38 
39  // Differences between knot points
40  std::vector<double> h(n, 0.);
41  for (int i = 0; i < n; i++)
42  h[i] = x[i + 1] - x[i];
43 
44  std::vector<double> mu(n, 0.);
45  std::vector<double> z(n + 1, 0.);
46  for (int i = 1; i < n; ++i) {
47  double g = 2. * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
48  mu[i] = h[i] / g;
49  z[i] =
50  (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]) - h[i - 1] * z[i - 1]) / g;
51  }
52 
53  // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
54  std::vector<double> a(n, 0.);
55  std::vector<double> b(n, 0.);
56  std::vector<double> c(n + 1, 0.);
57  std::vector<double> d(n, 0.);
58 
59  z[n] = 0.;
60  c[n] = 0.;
61 
62  for (int j = n - 1; j >= 0; j--) {
63  a[j] = y[j];
64  c[j] = z[j] - mu[j] * c[j + 1];
65  b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2. * c[j]) / 3.;
66  d[j] = (c[j + 1] - c[j]) / (3. * h[j]);
67  }
68 
69  // The above were taken from SplineInterpolator from Apache commons math. These
70  // polynomials need to be shifted by -x[i] in our case.
71  for (int i = 0; i < n; i++) {
72  double x_1 = -x[i];
73  double x_2 = x_1 * x_1;
74  double x_3 = x_1 * x_2;
75  a[i] = a[i] + b[i] * x_1 + c[i] * x_2 + d[i] * x_3;
76  b[i] = b[i] + 2. * c[i] * x_1 + 3. * d[i] * x_2;
77  c[i] = c[i] + 3. * d[i] * x_1;
78  // d[i] keeps the same value
79  }
80 
82  functions.reserve(n);
83  for (int i = 0; i < n; i++) {
84  functions.emplace_back(std::unique_ptr<Function>(new Polynomial{{a[i], b[i], c[i], d[i]}}));
85  }
86 
87  if (extrapolate) {
88  std::vector<double> x_copy(x);
91  return std::unique_ptr<Function>(new Piecewise{x_copy, std::move(functions)});
92  }
93 
94  return std::unique_ptr<Function>(new Piecewise{x, std::move(functions)});
95 }
96 
97 } // namespace MathUtils
98 } // end of namespace Euclid
std::move
T move(T... args)
std::vector::reserve
T reserve(T... args)
interpolation.h
std::vector< double >
std::vector::size
T size(T... args)
Euclid::MathUtils::Piecewise
Represents a piecewise function.
Definition: Piecewise.h:48
std::vector::back
T back(T... args)
g
constexpr double g
std::vector::front
T front(T... args)
Piecewise.h
std::numeric_limits::lowest
T lowest(T... args)
Euclid::MathUtils::Polynomial
Represents a polynomial function.
Definition: Polynomial.h:43
Exception.h
Polynomial.h
std::vector::emplace_back
T emplace_back(T... args)
Euclid::MathUtils::splineInterpolation
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:34
std::numeric_limits::max
T max(T... args)
std::unique_ptr
STL class.
Euclid
Definition: InstOrRefHolder.h:29