Mantid
Loading...
Searching...
No Matches
BSpline.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//----------------------------------------------------------------------
14
15#include <algorithm>
16
18
19using namespace CurveFitting;
20using namespace Kernel;
21using namespace API;
22
23DECLARE_FUNCTION(BSpline)
24
25
29 const size_t nbreak = 10;
30
31 auto orderValidator = BoundedValidator<int>();
32 orderValidator.setLower(1);
33 declareAttribute("Order", Attribute(3), orderValidator);
34
35 auto NBreakValidator = BoundedValidator<int>();
36 NBreakValidator.setLower(2);
37 declareAttribute("NBreak", Attribute(static_cast<int>(nbreak)), NBreakValidator);
38
39 auto startXValidator = BoundedValidator<double>();
40 startXValidator.setUpper(1.0);
41 startXValidator.setUpperExclusive(true);
42 declareAttribute("StartX", Attribute(0.0), startXValidator);
43
44 auto endXValidator = BoundedValidator<double>();
45 endXValidator.setLower(0.0);
46 endXValidator.setLowerExclusive(true);
47 declareAttribute("EndX", Attribute(1.0), endXValidator);
48
49 auto breakPointsValidator = ArrayBoundedValidator<double>(0.0, 1.0);
50 breakPointsValidator.setError(1e-8);
51 declareAttribute("BreakPoints", Attribute(std::vector<double>(nbreak)), breakPointsValidator);
52
53 declareAttribute("Uniform", Attribute(true));
54
55 m_spline = Spline1D();
56 resetParameters();
57 resetKnots();
58}
59
61 auto attStartX = getAttribute("StartX");
62 auto attEndX = getAttribute("EndX");
63
64 auto startXValidator = dynamic_cast<BoundedValidator<double> *>(attStartX.getValidator().get());
65 startXValidator->setUpper(attEndX.asDouble());
66
67 auto endXValidator = dynamic_cast<BoundedValidator<double> *>(attEndX.getValidator().get());
68 endXValidator->setLower(attStartX.asDouble());
69
70 auto breakPointsValidator =
71 dynamic_cast<ArrayBoundedValidator<double> *>(getAttribute("BreakPoints").getValidator().get());
72 breakPointsValidator->setLower(attStartX.asDouble());
73 breakPointsValidator->setUpper(attEndX.asDouble());
74}
75
82void BSpline::function1D(double *out, const double *xValues, const size_t nData) const {
83 size_t np = nParams();
84 EigenVector B(np);
85 size_t currentBBase = 0;
86 double startX = getAttribute("StartX").asDouble();
87 double endX = getAttribute("EndX").asDouble();
88
89 for (size_t i = 0; i < nData; ++i) {
90 double x = xValues[i];
91 if (x < startX || x > endX) {
92 out[i] = 0.0;
93 } else {
94 currentBBase = evaluateBasisFunctions(B, x, currentBBase);
95
96 double val = 0.0;
97 for (size_t j = 0; j < np; ++j) {
98 val += getParameter(j) * B.get(j);
99 }
100 out[i] = val;
101 }
102 }
103}
104
109void BSpline::initialiseSpline(const std::vector<double> &breakPoints) {
110 m_spline = Spline1D(EigenVector(m_knots).mutator(), EigenVector(breakPoints).mutator());
111}
112
121size_t BSpline::evaluateBasisFunctions(EigenVector &B, const double x, size_t currentBBase) const {
122 int degree = getDegree();
123 auto res = m_spline.BasisFunctions(x, degree, m_spline.knots());
124
125 currentBBase = getSpanIndex(x, currentBBase);
126 B.zero();
127 for (int i = 0; i < res.size(); i++) { // Populate B
128 B[currentBBase + i] = res.data()[i];
129 }
130 return currentBBase;
131}
132
142size_t BSpline::getSpanIndex(const double x, const size_t currentBBase, const bool clamped) const {
143 const size_t clampedKnots = clamped ? static_cast<size_t>(getClampedKnots()) : 1u;
144 size_t nKnots = m_knots.size();
145 for (size_t i = currentBBase + clampedKnots; i < nKnots; i++) {
146 if (x < m_knots[i]) {
147 return i - clampedKnots;
148 }
149 }
150 return nKnots - clampedKnots * 2;
151}
152
160void BSpline::derivative1D(double *out, const double *xValues, size_t nData, const size_t order) const {
161 double startX = getAttribute("StartX").asDouble();
162 double endX = getAttribute("EndX").asDouble();
163 size_t jstart = 0;
164 size_t splineOrder = static_cast<size_t>(getAttribute("Order").asInt());
165
166 for (size_t i = 0; i < nData; ++i) {
167 double x = xValues[i];
168 if (x < startX || x > endX) {
169 out[i] = 0.0;
170 } else {
171 jstart = getSpanIndex(x, jstart);
172 size_t jend = jstart + splineOrder;
173
174 auto res = evaluateBasisFnDerivatives(x, order);
175 double val = 0.0;
176 for (size_t j = jstart; j < jend; ++j) {
177 val += getParameter(j) * res(order, j - jstart);
178 }
179 out[i] = val;
180 }
181 }
182}
183
190EigenMatrix BSpline::evaluateBasisFnDerivatives(const double x, const size_t derivOrder) const {
191 int degree = getDegree();
192 EigenMatrix res(derivOrder + 1u, degree + 1);
193 res.mutator() = m_spline.BasisFunctionDerivatives(x, derivOrder, degree, m_spline.knots());
194 return res;
195}
196
202void BSpline::setAttribute(const std::string &attName, const API::IFunction::Attribute &att) {
203 bool isUniform = attName == "Uniform" && att.asBool();
204 storeAttributeValue(attName, att);
205
206 if (attName == "BreakPoints") {
207 storeAttributeValue("NBreak", Attribute(static_cast<int>(att.asVector().size())));
209 resetKnots();
210 } else if (isUniform || attName == "StartX" || attName == "EndX") {
212 resetKnots();
213 } else if (attName == "NBreak" || attName == "Order") {
215 resetKnots();
216 }
217}
218
222std::vector<std::string> BSpline::getAttributeNames() const {
223 return {"Uniform", "Order", "NBreak", "EndX", "StartX", "BreakPoints"};
224}
225
230 if (nParams() > 0) {
232 }
233 size_t np = getNBSplineCoefficients();
234 for (size_t i = 0; i < np; ++i) {
235 std::string pname = "A" + std::to_string(i);
236 declareParameter(pname);
237 }
238}
239
244 bool isUniform = getAttribute("Uniform").asBool();
245 std::vector<double> breakPoints;
246 if (isUniform) {
247 // create uniform knots in the interval [StartX, EndX]
248 double startX = getAttribute("StartX").asDouble();
249 double endX = getAttribute("EndX").asDouble();
250 // calc uniform break points
251 breakPoints = calcUniformBreakPoints(startX, endX);
252 storeAttributeValue("BreakPoints", Attribute(breakPoints));
253 // calc uniform knots
254 resetKnotVector(breakPoints);
255 } else {
256 // set the break points from BreakPoints vector attribute, update other attributes
257 breakPoints = getAttribute("BreakPoints").asVector();
258 // check that points are in ascending order
259 double prev = breakPoints[0];
260 for (size_t i = 1; i < breakPoints.size(); ++i) {
261 double next = breakPoints[i];
262 if (next < prev) {
263 throw std::invalid_argument("BreakPoints must be in ascending order.");
264 }
265 prev = next;
266 }
267 int nbreaks = getNBreakPoints();
268 // if number of break points change do necessary updates
269 if (static_cast<size_t>(nbreaks) != breakPoints.size()) {
270 storeAttributeValue("NBreak", Attribute(static_cast<int>(breakPoints.size())));
272 }
273 resetKnotVector(breakPoints);
274 }
275 initialiseSpline(breakPoints);
276}
277
284std::vector<double> BSpline::calcUniformBreakPoints(const double startX, const double endX) {
285 const int nBreak = getNBreakPoints();
286 std::vector<double> breakPoints(nBreak);
287 const double interval = (endX - startX) / (nBreak - 1.0);
288 std::generate(breakPoints.begin(), breakPoints.end(),
289 [n = 0, &interval, &startX]() mutable { return n++ * interval + startX; });
290 return breakPoints;
291}
292
297void BSpline::resetKnotVector(const std::vector<double> &breakPoints) {
298 const int nKnots = getNKnots();
299 const int clampedKnots = getClampedKnots();
300 m_knots.clear();
301 m_knots.resize(nKnots);
302
303 for (int i = 0; i < nKnots; i++) {
304 if (i < clampedKnots) {
305 m_knots[i] = breakPoints[0];
306 } else if (i >= nKnots - clampedKnots) {
307 m_knots[i] = breakPoints[breakPoints.size() - 1];
308 } else {
309 m_knots[i] = breakPoints[i - clampedKnots + 1];
310 }
311 }
312}
313
314} // namespace Mantid::CurveFitting::Functions
#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
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
std::vector< double > asVector() const
Returns vector<double> if attribute is vector<double>, throws exception otherwise.
Definition: IFunction.cpp:765
int asInt() const
Returns int value if attribute is a int, throws exception otherwise.
Definition: IFunction.cpp:726
Kernel::IValidator_sptr getValidator()
Return a clone of the attribute validator;.
Definition: IFunction.h:317
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
bool asBool() const
Returns bool value if attribute is a bool, throws exception otherwise.
Definition: IFunction.cpp:752
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition: IFunction.cpp:1394
void storeAttributeValue(const std::string &name, const API::IFunction::Attribute &value)
Store an attribute's value.
Definition: IFunction.cpp:1464
void clearAllParameters()
Nonvirtual member which removes all declared parameters.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
size_t nParams() const override
Total number of parameters.
Definition: ParamFunction.h:53
double getParameter(size_t i) const override
Get i-th parameter.
A wrapper around Eigen::Matrix.
Definition: EigenMatrix.h:33
map_type & mutator()
Get the map to Eigen matrix.
Definition: EigenMatrix.h:56
A wrapper around Eigen::Vector.
Definition: EigenVector.h:27
double get(const size_t i) const
Get an element.
A wrapper around Eigen functions implementing a B-spline.
Definition: BSpline.h:28
void resetKnotVector(const std::vector< double > &breakPoints)
Reset knot vector based upon break points.
Definition: BSpline.cpp:297
size_t evaluateBasisFunctions(EigenVector &B, const double x, size_t currentBBase) const
evaluate non-zero basis functions, return which index to use as the base of the results vector.
Definition: BSpline.cpp:121
std::vector< double > calcUniformBreakPoints(const double startX, const double endX)
Populate a vector with a uniform set of break points.
Definition: BSpline.cpp:284
void initialiseSpline(const std::vector< double > &breakPoints)
initialise the m_spline variable based on the knot vector and given a vector of breakpoints
Definition: BSpline.cpp:109
void function1D(double *out, const double *xValues, const size_t nData) const override
Execute the function.
Definition: BSpline.cpp:82
int getNBreakPoints() const
Get number of Break Points.
Definition: BSpline.h:65
void setAttribute(const std::string &attName, const Attribute &) override
Set a value to attribute attName.
Definition: BSpline.cpp:202
int getClampedKnots() const
Get number of knots required to clamp the spline closed.
Definition: BSpline.h:71
void resetKnots()
Reset b-spline knots.
Definition: BSpline.cpp:243
int getDegree() const
Get the degree of constituent polynomial functions.
Definition: BSpline.h:73
void resetValidators()
Reset Function Attribute Validators.
Definition: BSpline.cpp:60
void resetParameters()
Reset fitting parameters.
Definition: BSpline.cpp:229
int getNBSplineCoefficients()
Get number of B-Spline coefficients.
Definition: BSpline.h:63
Spline1D m_spline
Member variable for spline.
Definition: BSpline.h:76
int getNKnots() const
Get number of Knots.
Definition: BSpline.h:67
size_t getSpanIndex(const double x, const size_t currentBBase, const bool clamped=true) const
get the index of the span/interval which x is in
Definition: BSpline.cpp:142
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.
Definition: BSpline.cpp:160
std::vector< std::string > getAttributeNames() const override
Returns a list of attribute names.
Definition: BSpline.cpp:222
EigenMatrix evaluateBasisFnDerivatives(const double x, const size_t derivOrder) const
Evaluate derivatives up to a specified order for each non-zero basis function.
Definition: BSpline.cpp:190
std::vector< double > m_knots
Definition: BSpline.h:77
Kernel/ArrayBoundedValidator.h.
void setLower(const TYPE &value) noexcept
Set lower bound value.
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
void setLower(const TYPE &value) noexcept
Set lower bound value.
void setUpper(const TYPE &value) noexcept
Set upper bound value.
Eigen::Spline< double, 1, Eigen::Dynamic > Spline1D
Definition: BSpline.h:22
std::string to_string(const wide_integer< Bits, Signed > &n)