Mantid
Loading...
Searching...
No Matches
Fit.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
19#include "MantidAPI/TableRow.h"
21
25
26#include <memory>
27
29
30// Register the class into the algorithm factory
32
33
34Fit::Fit() : IFittingAlgorithm(), m_maxIterations() {}
35
38void Fit::initConcrete() {
39
41 getPointerToProperty("Ties")->setDocumentation("Math expressions defining ties between parameters of "
42 "the fitting function.");
43 declareProperty("Constraints", "", Kernel::Direction::Input);
44 getPointerToProperty("Constraints")->setDocumentation("List of constraints");
45 auto mustBePositive = std::make_shared<Kernel::BoundedValidator<int>>();
46 mustBePositive->setLower(0);
47 declareProperty("MaxIterations", 500, mustBePositive->clone(),
48 "Stop after this number of iterations if a good fit is not found");
49 declareProperty("OutputStatus", "", Kernel::Direction::Output);
50 getPointerToProperty("OutputStatus")->setDocumentation("Whether the fit was successful");
51 declareProperty("OutputChi2overDoF", 0.0, "Returns the goodness of the fit", Kernel::Direction::Output);
52
53 std::vector<std::string> minimizerOptions = API::FuncMinimizerFactory::Instance().getKeys();
54 Kernel::IValidator_sptr minimizerValidator = std::make_shared<Kernel::StartsWithValidator>(minimizerOptions);
55
56 declareProperty("Minimizer", "Levenberg-Marquardt", minimizerValidator, "Minimizer to use for fitting.");
57
58 std::vector<std::string> costFuncOptions = API::CostFunctionFactory::Instance().getKeys();
59 // select only CostFuncFitting variety
60 for (auto &costFuncOption : costFuncOptions) {
61 auto costFunc = std::dynamic_pointer_cast<CostFunctions::CostFuncFitting>(
63 if (!costFunc) {
64 costFuncOption = "";
65 }
66 }
67 Kernel::IValidator_sptr costFuncValidator = std::make_shared<Kernel::ListValidator<std::string>>(costFuncOptions);
68 declareProperty("CostFunction", "Least squares", costFuncValidator,
69 "The cost function to be used for the fit, default is Least squares", Kernel::Direction::InOut);
70 declareProperty("CreateOutput", false,
71 "Set to true to create output workspaces with the results of the fit"
72 "(default is false).");
73 declareProperty("Output", "",
74 "A base name for the output workspaces (if not "
75 "given default names will be created). The "
76 "default is to use the name of the original data workspace as prefix "
77 "followed by suffixes _Workspace, _Parameters, etc.");
78 declareProperty("CalcErrors", false,
79 "Set to true to calcuate errors when output isn't created "
80 "(default is false).");
81 declareProperty("OutputCompositeMembers", false,
82 "If true and CreateOutput is true then the value of each "
83 "member of a Composite Function is also output.");
84 declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("ConvolveMembers", false),
85 "If true and OutputCompositeMembers is true members of any "
86 "Convolution are output convolved\n"
87 "with corresponding resolution");
88 declareProperty("OutputParametersOnly", false,
89 "Set to true to output only the parameters and not "
90 "workspace(s) with the calculated values\n"
91 "(default is false, ignored if CreateOutput is false and "
92 "Output is an empty string).");
93}
94
95std::map<std::string, std::string> Fit::validateInputs() {
96 std::map<std::string, std::string> issues;
97
98 const auto &possibleOperators = Mantid::API::Expression::DEFAULT_OPS_STR;
99 std::string constraints = getPropertyValue("Constraints");
100 if (constraints.size() > 0) {
101 auto operatorPresent = false;
102 for (const auto &op : possibleOperators) {
103 const auto it = constraints.find_first_of(op);
104 if (it <= constraints.size()) {
105 operatorPresent = true;
106 break;
107 }
108 }
109 if (!operatorPresent) {
110 issues["Constraints"] = "No operator is present in the constraint.";
111 }
112 }
113
114 return issues;
115}
116
118void Fit::readProperties() {
119 std::string ties = getPropertyValue("Ties");
120 if (!ties.empty()) {
121 m_function->addTies(ties);
122 }
123 std::string constraints = getPropertyValue("Constraints");
124 if (!constraints.empty()) {
125 m_function->addConstraints(constraints);
126 }
127 m_function->registerFunctionUsage(isChild());
128
129 // Try to retrieve optional properties
130 int intMaxIterations = getProperty("MaxIterations");
131 m_maxIterations = static_cast<size_t>(intMaxIterations);
132}
133
136void Fit::initializeMinimizer(size_t maxIterations) {
137 const bool unrollComposites = getProperty("OutputCompositeMembers");
138 bool convolveMembers = existsProperty("ConvolveMembers");
139 if (convolveMembers) {
140 convolveMembers = getProperty("ConvolveMembers");
141 }
142 m_domainCreator->separateCompositeMembersInOutput(unrollComposites, convolveMembers);
144 std::string minimizerName = getPropertyValue("Minimizer");
145 m_minimizer = API::FuncMinimizerFactory::Instance().createMinimizer(minimizerName);
146 m_minimizer->initialize(m_costFunction, maxIterations);
147}
148
153void Fit::copyMinimizerOutput(const API::IFuncMinimizer &minimizer) {
154 const auto &properties = minimizer.getProperties();
155 for (auto property : properties) {
156 if ((*property).direction() == Kernel::Direction::Output && (*property).isValid().empty()) {
157 auto clonedProperty = std::unique_ptr<Kernel::Property>((*property).clone());
158 declareProperty(std::move(clonedProperty));
159 }
160 }
161}
162
165size_t Fit::runMinimizer() {
166 const int64_t nsteps = m_maxIterations * m_function->estimateNoProgressCalls();
167 auto prog = std::make_shared<API::Progress>(this, 0.0, 1.0, nsteps);
168 m_function->setProgressReporter(prog);
169
170 // do the fitting until success or iteration limit is reached
171 size_t iter = 0;
172 bool isFinished = false;
173 g_log.debug("Starting minimizer iteration\n");
174 while (iter < m_maxIterations) {
175 g_log.debug() << "Starting iteration " << iter << "\n";
176 try {
177 // Perform a single iteration. isFinished is set when minimizer wants to
178 // quit.
179 m_function->iterationStarting();
180 isFinished = !m_minimizer->iterate(iter);
181 m_function->iterationFinished();
183 // This is an attempt to recover after the function changes its number of
184 // parameters or ties during the iteration.
185 if (auto cf = dynamic_cast<API::CompositeFunction *>(m_function.get())) {
186 // Make sure the composite function is valid.
187 cf->checkFunction();
188 }
189 // Re-create the cost function and minimizer.
191 }
192
193 prog->report();
194 ++iter;
195 if (isFinished) {
196 // It was the last iteration. Break out of the loop and return the number
197 // of finished iterations.
198 break;
199 }
200 }
201 g_log.debug() << "Number of minimizer iterations=" << iter << "\n";
202 return iter;
203}
204
207void Fit::finalizeMinimizer(size_t nIterations) {
208 m_minimizer->finalize();
209
210 auto errorString = m_minimizer->getError();
211 g_log.debug() << "Iteration stopped. Minimizer status string=" << errorString << "\n";
212
213 if (nIterations >= m_maxIterations) {
214 if (!errorString.empty()) {
215 errorString += '\n';
216 }
217 errorString += "Failed to converge after " + std::to_string(m_maxIterations) + " iterations.";
218 }
219
220 if (errorString.empty()) {
221 errorString = "success";
222 }
223
224 // return the status flag
225 setPropertyValue("OutputStatus", errorString);
226 if (!this->isChild()) {
227 auto &logStream = errorString == "success" ? g_log.notice() : g_log.warning();
228 logStream << "Fit status: " << errorString << '\n';
229 logStream << "Stopped after " << nIterations << " iterations" << '\n';
230 }
231}
232
234void Fit::createOutput() {
235
236 // degrees of freedom
237 size_t dof = m_costFunction->getDomain()->size() - m_costFunction->nParams();
238 if (dof == 0)
239 dof = 1;
240 double rawcostfuncval = m_minimizer->costFunctionVal();
241 double finalCostFuncVal = rawcostfuncval / double(dof);
242
243 setProperty("OutputChi2overDoF", finalCostFuncVal);
244
245 bool doCreateOutput = getProperty("CreateOutput");
246 std::string baseName = getPropertyValue("Output");
247 if (!baseName.empty()) {
248 doCreateOutput = true;
249 }
250 bool doCalcErrors = getProperty("CalcErrors");
251 if (doCreateOutput) {
252 doCalcErrors = true;
253 }
254 if (m_costFunction->nParams() == 0) {
255 doCalcErrors = false;
256 }
257
258 EigenMatrix covar;
259 if (doCalcErrors) {
260 // Calculate the covariance matrix and the errors.
261 m_costFunction->calCovarianceMatrix(covar);
262 m_costFunction->calFittingErrors(covar, rawcostfuncval);
263 }
264
265 if (doCreateOutput) {
267
268 // get the workspace
269 API::Workspace_const_sptr ws = getProperty("InputWorkspace");
270
271 if (baseName.empty()) {
272 baseName = ws->getName();
273 if (baseName.empty()) {
274 baseName = "Output";
275 }
276 }
277 baseName += "_";
278
279 declareProperty(std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>("OutputNormalisedCovarianceMatrix",
281 "The name of the TableWorkspace in which to store the final covariance "
282 "matrix");
283 setPropertyValue("OutputNormalisedCovarianceMatrix", baseName + "NormalisedCovarianceMatrix");
284
286 Mantid::API::WorkspaceFactory::Instance().createTable("TableWorkspace");
287 covariance->addColumn("str", "Name");
288 // set plot type to Label = 6
289 covariance->getColumn(covariance->columnCount() - 1)->setPlotType(6);
290 for (size_t i = 0; i < m_function->nParams(); i++) {
291 if (m_function->isActive(i)) {
292 covariance->addColumn("double", m_function->parameterName(i));
293 }
294 }
295
296 size_t np = m_function->nParams();
297 size_t ia = 0;
298 for (size_t i = 0; i < np; i++) {
299 if (!m_function->isActive(i))
300 continue;
301 Mantid::API::TableRow row = covariance->appendRow();
302 row << m_function->parameterName(i);
303 size_t ja = 0;
304 for (size_t j = 0; j < np; j++) {
305 if (!m_function->isActive(j))
306 continue;
307 if (j == i)
308 row << 100.0;
309 else {
310 if (!covar.inspector().data()) {
311 throw std::runtime_error("There was an error while allocating the covariance "
312 "matrix "
313 "which is needed to produce fitting error results.");
314 }
315 row << 100.0 * covar.get(ia, ja) / sqrt(covar.get(ia, ia) * covar.get(ja, ja));
316 }
317 ++ja;
318
319 if (ja >= covar.size2())
320 break;
321 }
322 ++ia;
323
324 if (ia >= covar.size1())
325 break;
326 }
327
328 setProperty("OutputNormalisedCovarianceMatrix", covariance);
329
330 // create output parameter table workspace to store final fit parameters
331 // including error estimates if derivative of fitting function defined
332
333 declareProperty(std::make_unique<API::WorkspaceProperty<API::ITableWorkspace>>("OutputParameters", "",
335 "The name of the TableWorkspace in which to store the "
336 "final fit parameters");
337
338 setPropertyValue("OutputParameters", baseName + "Parameters");
339
340 Mantid::API::ITableWorkspace_sptr result = Mantid::API::WorkspaceFactory::Instance().createTable("TableWorkspace");
341 result->addColumn("str", "Name");
342 // set plot type to Label = 6
343 result->getColumn(result->columnCount() - 1)->setPlotType(6);
344 result->addColumn("double", "Value");
345 result->addColumn("double", "Error");
346 // yErr = 5
347 result->getColumn(result->columnCount() - 1)->setPlotType(5);
348
349 for (size_t i = 0; i < m_function->nParams(); i++) {
350 Mantid::API::TableRow row = result->appendRow();
351 row << m_function->parameterName(i) << m_function->getParameter(i) << m_function->getError(i);
352 }
353 // Add chi-squared value at the end of parameter table
354 Mantid::API::TableRow row = result->appendRow();
355
356 std::string costfuncname = getPropertyValue("CostFunction");
357 if (costfuncname == "Rwp")
358 row << "Cost function value" << rawcostfuncval;
359 else
360 row << "Cost function value" << finalCostFuncVal;
361
362 setProperty("OutputParameters", result);
363 bool outputParametersOnly = getProperty("OutputParametersOnly");
364
365 if (!outputParametersOnly) {
366 m_domainCreator->createOutputWorkspace(baseName, m_function, m_costFunction->getDomain(),
367 m_costFunction->getValues());
368 }
369 }
370}
371
376void Fit::execConcrete() {
377
378 // Read Fit's own properties
380
381 // Get the minimizer
383
384 // Run the minimizer
385 auto nIterations = runMinimizer();
386
387 // Finilize the minimizer.
388 finalizeMinimizer(nIterations);
389
390 // fit ended, creating output
391 createOutput();
392
393 progress(1.0);
394}
395
396} // namespace Mantid::CurveFitting::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
bool existsProperty(const std::string &name) const override
Checks whether the named property is already in the list of managed property.
Definition: Algorithm.cpp:2008
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
bool isChild() const override
To query whether algorithm is a child.
Definition: Algorithm.cpp:160
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
Definition: Algorithm.cpp:1975
A composite function is a function containing other functions.
static const std::vector< std::string > DEFAULT_OPS_STR
Definition: Expression.h:119
An interface for function minimizers.
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
A generic fitting algorithm.
Definition: Fit.h:78
void finalizeMinimizer(size_t nIterations)
Finalize the minimizer.
Definition: Fit.cpp:207
std::shared_ptr< CostFunctions::CostFuncFitting > m_costFunction
The cost function.
Definition: Fit.h:103
void copyMinimizerOutput(const API::IFuncMinimizer &minimizer)
Copy all output workspace properties from the minimizer to Fit algorithm.
Definition: Fit.cpp:153
size_t m_maxIterations
Max number of iterations.
Definition: Fit.h:107
size_t runMinimizer()
Run the minimizer's iteration loop.
Definition: Fit.cpp:165
void createOutput()
Create algorithm output workspaces.
Definition: Fit.cpp:234
std::shared_ptr< API::IFuncMinimizer > m_minimizer
The minimizer.
Definition: Fit.h:105
void readProperties()
Read in the properties specific to Fit.
Definition: Fit.cpp:118
void initializeMinimizer(size_t maxIterations)
Initialize the minimizer for this fit.
Definition: Fit.cpp:136
A wrapper around Eigen::Matrix.
Definition: EigenMatrix.h:33
double get(size_t i, size_t j) const
Get an element.
size_t size1() const
First size of the matrix.
size_t size2() const
Second size of the matrix.
const map_type inspector() const
Get a const copy of the Eigen matrix.
Definition: EigenMatrix.h:58
A base class for fitting algorithms.
std::shared_ptr< API::IDomainCreator > m_domainCreator
Pointer to a domain creator.
std::shared_ptr< CostFunctions::CostFuncFitting > getCostFunctionInitialized() const
Create a cost function from the "CostFunction" property and make it ready for evaluation.
std::shared_ptr< API::IFunction > m_function
Pointer to the fitting function.
Exception thrown when a fitting function changes number of parameters during fit.
Definition: Exception.h:336
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
const std::vector< Property * > & getProperties() const override
Get the list of managed properties.
The concrete, templated class for properties.
void setDocumentation(const std::string &documentation)
Sets the user level description of the property.
Definition: Property.cpp:134
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< const Workspace > Workspace_const_sptr
shared pointer to Mantid::API::Workspace (const version)
Definition: Workspace_fwd.h:22
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::shared_ptr< IValidator > IValidator_sptr
A shared_ptr to an IValidator.
Definition: IValidator.h:26
std::string to_string(const wide_integer< Bits, Signed > &n)
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54