Mantid
Loading...
Searching...
No Matches
TabulatedFunction.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//----------------------------------------------------------------------
11#include "MantidAPI/Algorithm.h"
16
17#include <algorithm>
18#include <cmath>
19#include <fstream>
20#include <sstream>
21#include <utility>
22
24
25using namespace CurveFitting;
26
27using namespace Kernel;
28
29using namespace API;
30
31namespace {
33Logger g_log("TabulatedFunction");
34} // namespace
35
36DECLARE_FUNCTION(TabulatedFunction)
37
39
41TabulatedFunction::TabulatedFunction() : m_setupFinished(false), m_explicitXY(false) {
42 declareParameter("Scaling", 1.0, "A scaling factor");
43 declareParameter("Shift", 0.0, "Shift in the abscissa");
44 declareParameter("XScaling", 1.0, "Scaling factor in X");
45 declareAttribute("FileName", Attribute("", true));
46 declareAttribute("Workspace", Attribute(""));
48 declareAttribute("X", Attribute(std::vector<double>()));
49 declareAttribute("Y", Attribute(std::vector<double>()));
50}
51
53void TabulatedFunction::eval(double scaling, double xshift, double xscale, double *out, const double *xValues,
54 const size_t nData) const {
55 if (nData == 0)
56 return;
57 setupData();
58
59 if (size() == 0)
60 return;
61
62 // shift and scale the domain over which the function is defined
63 std::vector<double> xData(m_xData);
64 for (double &value : xData) {
65 value *= xscale;
66 value += xshift;
67 }
68 const double xStart = xData.front();
69 const double xEnd = xData.back();
70
71 if (xStart >= xValues[nData - 1] || xEnd <= xValues[0])
72 return;
73 size_t i = 0;
74 while (i < nData - 1 && xValues[i] < xStart) {
75 out[i] = 0;
76 i++;
77 }
78 size_t j = 0;
79 for (; i < nData; i++) {
80 double xi = xValues[i];
81 while (j < size() - 1 && xi > xData[j])
82 j++;
83 if (j > size() - 1) {
84 out[i] = 0;
85 } else {
86 if (xi == xData[j]) {
87 out[i] = m_yData[j] * scaling;
88 } else if (xi > xData[j]) {
89 out[i] = 0;
90 } else if (j > 0) {
91 double x0 = xData[j - 1];
92 double x1 = xData[j];
93 double y0 = m_yData[j - 1];
94 double y1 = m_yData[j];
95 out[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
96 out[i] *= scaling;
97 } else {
98 out[i] = 0;
99 }
100 }
101 }
102}
103
110void TabulatedFunction::function1D(double *out, const double *xValues, const size_t nData) const {
111 const double scaling = getParameter("Scaling");
112 const double xshift = getParameter("Shift");
113 const double xscale = getParameter("XScaling");
114 eval(scaling, xshift, xscale, out, xValues, nData);
115}
116
124void TabulatedFunction::functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) {
125 const double scaling = getParameter("Scaling");
126 const double xshift = getParameter("Shift");
127 const double xscale = getParameter("XScaling");
128 std::vector<double> tmp(nData);
129 // derivative with respect to Scaling parameter
130 eval(1.0, xshift, xscale, tmp.data(), xValues, nData);
131 for (size_t i = 0; i < nData; ++i) {
132 out->set(i, 0, tmp[i]);
133 }
134
135 const double dx = (xValues[nData - 1] - xValues[0]) / static_cast<double>(nData);
136 std::vector<double> tmpplus(nData);
137 std::vector<double> tmpminus(nData);
138
139 // There is no unique definition for the partial derivative with respect
140 // to the Shift parameter. Here we take the central difference,
141 eval(scaling, xshift + dx, xscale, tmpplus.data(), xValues, nData);
142 eval(scaling, xshift - dx, xscale, tmpminus.data(), xValues, nData);
143 for (size_t i = 0; i < nData; ++i) {
144 out->set(i, 1, (tmpplus[i] - tmpminus[i]) / (2 * dx));
145 }
146
147 eval(scaling, xshift, xscale + dx, tmpplus.data(), xValues, nData);
148 eval(scaling, xshift, xscale - dx, tmpminus.data(), xValues, nData);
149 for (size_t i = 0; i < nData; ++i) {
150 out->set(i, 2, (tmpplus[i] - tmpminus[i]) / (2 * dx));
151 }
152}
153
156 m_xData.clear();
157 m_yData.clear();
158 m_setupFinished = false;
159}
160
165void TabulatedFunction::setAttribute(const std::string &attName, const IFunction::Attribute &value) {
166 if (attName == "FileName") {
167 std::string fileName = value.asUnquotedString();
168 if (fileName.empty()) {
169 storeAttributeValue("FileName", Attribute("", true));
170 return;
171 }
172 FileValidator fval;
173 std::string error = fval.isValid(fileName);
174 if (error.empty()) {
175 storeAttributeValue(attName, Attribute(fileName, true));
176 storeAttributeValue("Workspace", Attribute(""));
177 } else {
178 // file not found
179 throw Kernel::Exception::FileError(error, fileName);
180 }
181 load(fileName);
182 m_setupFinished = false;
183 m_explicitXY = false;
184 } else if (attName == "Workspace") {
185 std::string wsName = value.asString();
186 if (!wsName.empty()) {
187 storeAttributeValue(attName, value);
188 storeAttributeValue("FileName", Attribute("", true));
189 loadWorkspace(wsName);
190 m_setupFinished = false;
191 m_explicitXY = false;
192 }
193 } else if (attName == "X") {
194 m_xData = value.asVector();
195 if (m_xData.empty()) {
196 m_setupFinished = false;
197 m_explicitXY = false;
198 if (!m_yData.empty()) {
199 m_yData.clear();
200 }
201 return;
202 }
203 if (m_xData.size() != m_yData.size()) {
204 m_yData.resize(m_xData.size());
205 }
206 storeAttributeValue("FileName", Attribute("", true));
207 storeAttributeValue("Workspace", Attribute(""));
208 m_setupFinished = true;
209 m_explicitXY = true;
210 } else if (attName == "Y") {
211 m_yData = value.asVector();
212 if (m_yData.empty()) {
213 m_setupFinished = false;
214 m_explicitXY = false;
215 if (!m_xData.empty()) {
216 m_xData.clear();
217 }
218 return;
219 }
220 if (m_xData.size() != m_yData.size()) {
221 m_xData.resize(m_yData.size());
222 }
223 storeAttributeValue("FileName", Attribute("", true));
224 storeAttributeValue("Workspace", Attribute(""));
225 m_setupFinished = true;
226 m_explicitXY = true;
227 } else {
229 m_setupFinished = false;
230 }
231}
232
235
237std::vector<std::string> TabulatedFunction::getAttributeNames() const { return IFunction::getAttributeNames(); }
238
241IFunction::Attribute TabulatedFunction::getAttribute(const std::string &attName) const {
242 if (attName == "X") {
243 return m_explicitXY ? Attribute(m_xData) : Attribute(std::vector<double>());
244 } else if (attName == "Y") {
245 return m_explicitXY ? Attribute(m_yData) : Attribute(std::vector<double>());
246 }
247 return IFunction::getAttribute(attName);
248}
253void TabulatedFunction::load(const std::string &fname) {
254 auto loadAlg = Mantid::API::AlgorithmFactory::Instance().create("Load", -1);
255 loadAlg->initialize();
256 loadAlg->setChild(true);
257 loadAlg->setLogging(false);
258 try {
259 loadAlg->setPropertyValue("Filename", fname);
260 loadAlg->setPropertyValue("OutputWorkspace", "_TabulatedFunction_fit_data_");
261 loadAlg->execute();
262 } catch (std::runtime_error &) {
263 throw std::runtime_error("Unable to load Nexus file for TabulatedFunction function.");
264 }
265
266 Workspace_sptr ws = loadAlg->getProperty("OutputWorkspace");
267 MatrixWorkspace_sptr resData = std::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>(ws);
268 loadWorkspace(resData);
269}
270
275void TabulatedFunction::loadWorkspace(const std::string &wsName) const {
276 auto ws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(wsName);
277 loadWorkspace(ws);
278 if (!m_workspace) {
279 throw std::runtime_error("Unable to set " + wsName + " as workspace attribute. Expected a MatrixWorkspace.");
280 }
281}
282
287void TabulatedFunction::loadWorkspace(std::shared_ptr<API::MatrixWorkspace> ws) const {
288 m_workspace = std::move(ws);
289 m_setupFinished = false;
290}
291
296 if (m_setupFinished) {
297 if (m_xData.size() != m_yData.size()) {
298 throw std::invalid_argument(this->name() + ": X and Y vectors have different sizes.");
299 }
300 g_log.debug() << "Re-setting isn't required.";
301 return;
302 }
303
304 if (!m_workspace) {
305 std::string wsName = getAttribute("Workspace").asString();
306 if (wsName.empty())
307 throw std::invalid_argument("Data not set for function " + this->name());
308 else
309 loadWorkspace(wsName);
310 }
311
312 size_t index = static_cast<size_t>(getAttribute("WorkspaceIndex").asInt());
313
314 g_log.debug() << "Setting up " << m_workspace->getName() << " index " << index << '\n';
315
316 const auto &xData = m_workspace->points(index);
317 const auto &yData = m_workspace->y(index);
318 m_xData.assign(xData.begin(), xData.end());
319 m_yData.assign(yData.begin(), yData.end());
320
321 m_workspace.reset();
322 m_setupFinished = true;
323}
324
325} // namespace Mantid::CurveFitting::Functions
gsl_vector * tmp
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.
double error
Definition: IndexPeaks.cpp:133
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
static Kernel::Logger g_log
Logger instance.
Definition: IFunction1D.h:74
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
std::string asString() const
Returns string value if attribute is a string, throws exception otherwise.
Definition: IFunction.cpp:660
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
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
Definition: IFunction.cpp:1409
virtual std::vector< std::string > getAttributeNames() const
Returns a list of attribute names.
Definition: IFunction.cpp:1368
void storeAttributeValue(const std::string &name, const API::IFunction::Attribute &value)
Store an attribute's value.
Definition: IFunction.cpp:1464
virtual size_t nAttributes() const
Returns the number of attributes associated with the function.
Definition: IFunction.cpp:1336
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
virtual void set(size_t iY, size_t iP, double value)=0
Set a value to a Jacobian matrix element.
Base MatrixWorkspace Abstract Class.
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.
static const int defaultIndexValue
The default value for the workspace index.
void load(const std::string &fname)
Call the appropriate load function.
size_t nAttributes() const override
Returns the number of attributes associated with the function.
std::shared_ptr< API::MatrixWorkspace > m_workspace
Temporary workspace holder.
std::vector< double > m_yData
Stores y-values.
void loadWorkspace(const std::string &wsName) const
Load the points from a MatrixWorkspace.
bool m_explicitXY
Flag of explicit x-y data setup.
bool m_setupFinished
Flag of completing data setup.
void eval(double scaling, double xshift, double xscale, double *out, const double *xValues, const size_t nData) const
Evaluate the function for a list of arguments and given scaling factor.
void function1D(double *out, const double *xValues, const size_t nData) const override
Calculate the function values.
std::vector< std::string > getAttributeNames() const override
Returns a list of attribute names.
std::string name() const override
overwrite IFunction base class methods
void setAttribute(const std::string &attName, const IFunction::Attribute &value) override
Set a value to attribute attName.
std::vector< double > m_xData
Stores x-values.
void setupData() const
Fill in the x and y value containers (m_xData and m_yData)
void functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) override
function derivatives
Attribute getAttribute(const std::string &attName) const override
Return a value of attribute attName.
Records the filename and the description of failure.
Definition: Exception.h:98
FileValidator is a validator that checks that a filepath is valid.
Definition: FileValidator.h:25
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class