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