-
Notifications
You must be signed in to change notification settings - Fork 121
/
UserFunction1D.cpp
162 lines (140 loc) · 5.27 KB
/
UserFunction1D.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/Functions/UserFunction1D.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/UnitFactory.h"
namespace Mantid {
namespace CurveFitting {
namespace Functions {
using namespace CurveFitting;
// Register the class into the algorithm factory
DECLARE_ALGORITHM(UserFunction1D)
using namespace Kernel;
using namespace API;
/** Static callback function used by MuParser to initialize variables implicitly
* @param varName :: The name of a new variable
* @param palg :: Pointer to the algorithm
* @return A pointer to the initialized variable
*/
double *UserFunction1D::AddVariable(const char *varName, void *palg) {
UserFunction1D &alg = *reinterpret_cast<UserFunction1D *>(palg);
if (std::string(varName) != "x") {
alg.declareProperty(varName, 0.0);
alg.m_parameterNames.emplace_back(varName);
} else {
alg.m_x_set = true;
alg.m_x = 0.;
return &alg.m_x;
}
return &alg.m_parameters[alg.m_nPars++];
}
/** Declare properties that are not fit parameters
*/
void UserFunction1D::declareAdditionalProperties() {
declareProperty("Function", "", std::make_shared<MandatoryValidator<std::string>>(), "The fit function");
declareProperty("InitialParameters", "",
"The comma separated list of initial values of the fit "
"parameters in the form varName=value");
}
/** Declare fit parameters using muParser's implicit variable initialization.
*/
void UserFunction1D::prepare() {
m_parser.SetVarFactory(AddVariable, this);
std::string funct = getProperty("Function");
m_parser.SetExpr(funct);
// Call Eval() to implicitly initialize the variables
m_parser.Eval();
if (!m_x_set)
throw std::runtime_error("Formula does not contain the x variable");
// Set the initial values to the fit parameters
std::string initParams = getProperty("InitialParameters");
if (!initParams.empty()) {
Mantid::Kernel::StringTokenizer values(initParams, ",", Mantid::Kernel::StringTokenizer::TOK_TRIM);
for (const auto &it : values) {
size_t ieq = it.find('=');
if (ieq == std::string::npos)
throw std::invalid_argument("Property InitialParameters is malformed");
std::string varName = it.substr(0, ieq);
std::string varValue = it.substr(ieq + 1);
size_t i0 = varName.find_first_not_of(" \t");
size_t i1 = varName.find_last_not_of(" \t");
if (i0 == std::string::npos)
throw std::invalid_argument("Property InitialParameters is malformed");
varName = varName.substr(i0, i1 - i0 + 1);
if (varName.empty() || varValue.empty())
throw std::invalid_argument("Property InitialParameters is malformed");
double value = std::stod(varValue);
if (!existsProperty(varName))
throw std::invalid_argument("Fit parameter " + varName + " does not exist");
setProperty(varName, value);
}
}
}
/*double UserFunction1D::function(const double* in, const double& x)
{
for(int i=0;i<m_nPars;i++)
m_parameters[i] = in[i];
m_x = x;
return m_parser.Eval();
}*/
/** Calculate the fitting function.
* @param in :: A pointer ot the input function parameters
* @param out :: A pointer to the output fitting function buffer. The buffer
* must be large enough to receive nData double values.
* The fitting procedure will try to minimise Sum(out[i]^2)
* @param xValues :: The array of nData x-values.
* @param nData :: The size of the fitted data.
*/
void UserFunction1D::function(const double *in, double *out, const double *xValues, const size_t nData) {
for (size_t i = 0; i < static_cast<size_t>(m_nPars); i++)
m_parameters[i] = in[i];
for (size_t i = 0; i < nData; i++) {
m_x = xValues[i];
out[i] = m_parser.Eval();
}
}
/**
* @param in :: Input fitting parameter values
* @param out :: Derivatives
* @param xValues :: X values for data points
* @param nData :: Number of data points
*/
void UserFunction1D::functionDeriv(const double *in, Jacobian *out, const double *xValues, const size_t nData) {
// throw Exception::NotImplementedError("No derivative function provided");
if (nData == 0)
return;
std::vector<double> dp(m_nPars);
std::vector<double> in1(m_nPars);
for (int i = 0; i < m_nPars; i++) {
in1[i] = in[i];
m_parameters[i] = in[i];
if (m_parameters[i] != 0.0)
dp[i] = m_parameters[i] * 0.01;
else
dp[i] = 0.01;
}
if (m_tmp.empty()) {
m_tmp.resize(nData);
m_tmp1.resize(nData);
}
function(in, m_tmp.data(), xValues, nData);
for (int j = 0; j < m_nPars; j++) {
in1[j] += dp[j];
function(&in1[0], m_tmp1.data(), xValues, nData);
for (size_t i = 0; i < nData; i++) {
out->set(i, j, (m_tmp1[i] - m_tmp[i]) / dp[j]);
}
in1[j] -= dp[j];
}
}
} // namespace Functions
} // namespace CurveFitting
} // namespace Mantid