Mantid
Loading...
Searching...
No Matches
CubicSpline.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
12#include "MantidKernel/Logger.h"
13
14#include <algorithm>
15#include <ostream>
16#include <stdexcept>
17#include <vector>
18
20
21using namespace CurveFitting;
22namespace {
24Kernel::Logger g_log("CubicSpline");
25} // namespace
26
27using namespace Kernel;
28
29using namespace API;
30
31DECLARE_FUNCTION(CubicSpline)
32
33
37 : m_min_points(3), m_acc(gsl_interp_accel_alloc(), m_gslFree),
38 m_spline(gsl_spline_alloc(gsl_interp_cspline, m_min_points), m_gslFree), m_recalculateSpline(true) {
39 // setup class with a default set of attributes
40 declareAttribute("n", Attribute(m_min_points));
41
42 // declare corresponding attributes and parameters
43 declareAttribute("x0", Attribute(0.0));
44 declareAttribute("x1", Attribute(1.0));
45 declareAttribute("x2", Attribute(2.0));
46
47 declareParameter("y0", 0);
48 declareParameter("y1", 0);
49 declareParameter("y2", 0);
50}
51
58void CubicSpline::function1D(double *out, const double *xValues, const size_t nData) const {
59 // check if spline needs recalculating
60 int n = getAttribute("n").asInt();
61
62 boost::scoped_array<double> x(new double[n]);
63 boost::scoped_array<double> y(new double[n]);
64
65 // setup the reference points and calculate
67 setupInput(x, y, n);
68
69 calculateSpline(out, xValues, nData);
70}
71
78void CubicSpline::setupInput(boost::scoped_array<double> &x, boost::scoped_array<double> &y, int n) const {
79 // Populate data points from the input attributes and parameters
80 bool isSorted = true;
81
82 for (int i = 0; i < n; ++i) {
83
84 x[i] = getAttribute("x" + std::to_string(i)).asDouble();
85 y[i] = getParameter(i);
86
87 if (isSorted) {
88 // if x[i] is out of order with its neighbours
89 if (i > 1 && i < n && (x[i - 1] < x[i - 2] || x[i - 1] > x[i])) {
90 isSorted = false;
91 }
92 }
93 }
94
95 // sort the data points if necessary
96 if (!isSorted) {
97 g_log.warning() << "Spline x parameters are not in ascending order. Values "
98 "will be sorted.\n";
99
100 using point = std::pair<double, double>;
101 std::vector<point> pairs;
102 pairs.reserve(n);
103 for (int i = 0; i < n; ++i) {
104 pairs.emplace_back(x[i], y[i]);
105 }
106
107 std::sort(pairs.begin(), pairs.end(), [](const point &xy1, const point &xy2) { return xy1.first < xy2.first; });
108
109 for (int i = 0; i < n; ++i) {
110 x[i] = pairs[i].first;
111 y[i] = pairs[i].second;
112 }
113 }
114
115 // pass values to GSL objects
116 initGSLObjects(x, y, n);
117 m_recalculateSpline = false;
118}
119
127void CubicSpline::derivative1D(double *out, const double *xValues, size_t nData, const size_t order) const {
128 int n = getAttribute("n").asInt();
129
130 boost::scoped_array<double> x(new double[n]);
131 boost::scoped_array<double> y(new double[n]);
132
133 // setup the reference points and calculate
135 setupInput(x, y, n);
136 calculateDerivative(out, xValues, nData, order);
137}
138
144bool CubicSpline::checkXInRange(double x) const { return (x >= m_spline->interp->xmin && x <= m_spline->interp->xmax); }
145
152void CubicSpline::calculateSpline(double *out, const double *xValues, const size_t nData) const {
153 // calculate spline for given input set
154 bool outOfRange(false);
155 for (size_t i = 0; i < nData; ++i) {
156 if (checkXInRange(xValues[i])) {
157 // calculate the y value
158 out[i] = splineEval(xValues[i]);
159 } else {
160 // if out of range, set it to constant of fringe values
161 outOfRange = true;
162 if (xValues[i] < m_spline->interp->xmin) {
163 out[i] = splineEval(m_spline->interp->xmin);
164 } else {
165 out[i] = splineEval(m_spline->interp->xmax);
166 }
167 }
168 }
169
170 // inform user that some values weren't calculated
171 if (outOfRange) {
172 g_log.information() << "Some x values where out of range and will not be calculated.\n";
173 }
174}
175
181double CubicSpline::splineEval(const double x) const {
182 // calculate the y value
183 double y = gsl_spline_eval(m_spline.get(), x, m_acc.get());
184 int errorCode = gsl_spline_eval_e(m_spline.get(), x, m_acc.get(), &y);
185
186 // check if GSL function returned an error
187 checkGSLError(errorCode, GSL_EDOM);
188
189 return y;
190}
191
199void CubicSpline::calculateDerivative(double *out, const double *xValues, const size_t nData,
200 const size_t order) const {
201 double xDeriv = 0;
202 int errorCode = 0;
203 bool outOfRange(false);
204
205 // throw error if the order is not the 1st or 2nd derivative
206 if (order < 1)
207 throw std::invalid_argument("CubicSpline: order of derivative must be 1 or greater");
208
209 for (size_t i = 0; i < nData; ++i) {
210 if (checkXInRange(xValues[i])) {
211 // choose the order of the derivative
212 if (order == 1) {
213 xDeriv = gsl_spline_eval_deriv(m_spline.get(), xValues[i], m_acc.get());
214 errorCode = gsl_spline_eval_deriv_e(m_spline.get(), xValues[i], m_acc.get(), &xDeriv);
215 } else if (order == 2) {
216 xDeriv = gsl_spline_eval_deriv2(m_spline.get(), xValues[i], m_acc.get());
217 errorCode = gsl_spline_eval_deriv2_e(m_spline.get(), xValues[i], m_acc.get(), &xDeriv);
218 }
219 } else {
220 // if out of range, just set it to zero
221 outOfRange = true;
222 xDeriv = 0;
223 }
224
225 // check GSL functions didn't return an error
226 checkGSLError(errorCode, GSL_EDOM);
227
228 // record the value
229 out[i] = xDeriv;
230 }
231
232 // inform user that some values weren't calculated
233 if (outOfRange) {
234 g_log.information() << "Some x values where out of range and will not be calculated.\n";
235 }
236}
237
244void CubicSpline::setParameter(size_t i, const double &value, bool explicitlySet) {
245 // Call parent setParameter implementation
246 ParamFunction::setParameter(i, value, explicitlySet);
247
248 // recalculate if necessary
249 m_recalculateSpline = true;
250}
251
257void CubicSpline::setAttribute(const std::string &attName, const API::IFunction::Attribute &att) {
258
259 if (attName == "n") {
260 // get the new and old number of data points
261 int n = att.asInt();
262 int oldN = getAttribute("n").asInt();
263
264 // check that the number of data points is in a valid range
265 if (n > oldN) {
266 // get the name of the last x data point
267 std::string oldXName = "x" + std::to_string(oldN - 1);
268 double oldX = getAttribute(oldXName).asDouble();
269
270 // reallocate gsl object to new size
272
273 // create blank a number of new blank parameters and attributes
274 for (int i = oldN; i < n; ++i) {
275 std::string num = std::to_string(i);
276
277 std::string newXName = "x" + num;
278 std::string newYName = "y" + num;
279
280 declareAttribute(newXName, Attribute(oldX + static_cast<double>(i - oldN + 1)));
281 declareParameter(newYName, 0);
282 }
283
284 // flag that the spline + derivatives will now need to be recalculated
285 m_recalculateSpline = true;
286 } else if (n < oldN) {
287 throw std::invalid_argument("Cubic Spline: Can't decrease the number of attributes");
288 }
289 }
290
291 storeAttributeValue(attName, att);
292}
293
299void CubicSpline::setXAttribute(const size_t index, double x) {
300 size_t n = static_cast<size_t>(getAttribute("n").asInt());
301
302 // check that setting the x attribute is within our range
303 if (index < n) {
304 std::string xName = "x" + std::to_string(index);
305 setAttributeValue(xName, x);
306
307 // attributes updated, flag for recalculation
308 m_recalculateSpline = true;
309 } else {
310 throw std::range_error("Cubic Spline: x index out of range.");
311 }
312}
313
320void CubicSpline::checkGSLError(const int status, const int errorType) const {
321 // check GSL functions didn't return an error
322 if (status == errorType) {
323 m_recalculateSpline = true;
324
325 std::string message("CubicSpline: ");
326 message.append(gsl_strerror(errorType));
327
328 throw std::runtime_error(message);
329 }
330}
331
338void CubicSpline::initGSLObjects(boost::scoped_array<double> &x, boost::scoped_array<double> &y, int n) const {
339 int status = gsl_spline_init(m_spline.get(), x.get(), y.get(), n);
340 checkGSLError(status, GSL_EINVAL);
341}
342
348 m_spline.reset(gsl_spline_alloc(gsl_interp_cspline, n), m_gslFree);
349 gsl_interp_accel_reset(m_acc.get());
350}
351
352} // namespace Mantid::CurveFitting::Functions
double value
The value of the point.
Definition: FitMW.cpp:51
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
Mantid::API::IFunction::Attribute Attribute
Definition: IsoRotDiff.cpp:25
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
int asInt() const
Returns int value if attribute is a int, throws exception otherwise.
Definition: IFunction.cpp:726
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition: IFunction.cpp:1394
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition: IFunction.cpp:1418
void setAttributeValue(const std::string &attName, const T &value)
Set an attribute value.
Definition: IFunction.h:597
void storeAttributeValue(const std::string &name, const API::IFunction::Attribute &value)
Store an attribute's value.
Definition: IFunction.cpp:1464
void setParameter(size_t, const double &value, bool explicitlySet=true) override
Set i-th parameter.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
double getParameter(size_t i) const override
Get i-th parameter.
A wrapper around GSL functions implementing cubic spline interpolation.
Definition: CubicSpline.h:31
bool m_recalculateSpline
Flag for checking if the spline needs recalculating.
Definition: CubicSpline.h:68
void setupInput(boost::scoped_array< double > &x, boost::scoped_array< double > &y, int n) const
Method to setup the gsl function.
Definition: CubicSpline.cpp:78
void setParameter(size_t i, const double &value, bool explicitlySet=true) override
Set a parameter for the function and flags the spline for re-calculation.
bool checkXInRange(double x) const
Check if an x value falls within the range of the spline.
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
struct Mantid::CurveFitting::Functions::CubicSpline::GSLFree m_gslFree
void reallocGSLObjects(const int n)
Reallocate the spline object to use n data points.
double splineEval(const double x) const
Evaluate a point on the spline, with basic error handling.
std::shared_ptr< gsl_spline > m_spline
GSL data structure used to calculate spline.
Definition: CubicSpline.h:65
void setXAttribute(const size_t index, double x)
Set the value of a data point location to x.
void function1D(double *out, const double *xValues, const size_t nData) const override
Execute the function.
Definition: CubicSpline.cpp:58
void checkGSLError(const int status, const int errorType) const
Check if an error occurred and throw appropriate message.
void calculateSpline(double *out, const double *xValues, const size_t nData) const
Calculate the spline.
std::shared_ptr< gsl_interp_accel > m_acc
GSL interpolation accelerator object.
Definition: CubicSpline.h:62
void calculateDerivative(double *out, const double *xValues, const size_t nData, const size_t order) const
Calculate the derivative.
void derivative1D(double *out, const double *xValues, size_t nData, const size_t order) const override
Calculate the derivatives for a set of points on the spline.
void initGSLObjects(boost::scoped_array< double > &x, boost::scoped_array< double > &y, int n) const
Initialise GSL objects if required.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::string to_string(const wide_integer< Bits, Signed > &n)