21#include <gsl/gsl_blas.h>
22#include <gsl/gsl_multifit_nlin.h>
23#include <gsl/gsl_multimin.h>
24#include <gsl/gsl_statistics.h>
25#include <gsl/gsl_version.h>
33using namespace Kernel;
35using API::MatrixWorkspace;
38using API::WorkspaceProperty;
54 void set(
size_t iY,
size_t iP,
double value)
override {
55 int j =
m_map[
static_cast<int>(iP)];
64 double get(
size_t iY,
size_t iP)
override {
65 int j =
m_map[
static_cast<int>(iP)];
67 return gsl_matrix_get(
m_J, iY, j);
72 void zero()
override { gsl_matrix_set_zero(
m_J); }
116static int gsl_f(
const gsl_vector *
x,
void *params, gsl_vector *f) {
117 auto fitParams =
static_cast<FitData *
>(params);
118 for (
size_t i = 0, j = 0; i < fitParams->active.size(); i++)
119 if (fitParams->active[i])
122 fitParams->fit1D->function(fitParams->parameters, f->data, fitParams->X, fitParams->n);
127 for (
size_t i = 0; i < fitParams->n; i++)
128 f->data[i] = (f->data[i] - fitParams->Y[i]) / fitParams->sigmaData[i];
139static int gsl_df(
const gsl_vector *
x,
void *params, gsl_matrix *J) {
140 auto fitParams =
static_cast<FitData *
>(params);
141 for (
size_t i = 0, j = 0; i < fitParams->active.size(); i++) {
142 if (fitParams->active[i])
146 fitParams->J.setJ(J);
148 fitParams->fit1D->functionDeriv(fitParams->parameters, &fitParams->J, fitParams->X, fitParams->n);
154 for (
size_t iY = 0; iY < fitParams->n; iY++)
155 for (
size_t iP = 0; iP < fitParams->p; iP++)
156 J->data[iY * fitParams->p + iP] /= fitParams->sigmaData[iY];
168static int gsl_fdf(
const gsl_vector *
x,
void *params, gsl_vector *f, gsl_matrix *J) {
182 auto fitParams =
static_cast<FitData *
>(params);
185 for (
size_t i = 0, j = 0; i < fitParams->active.size(); i++)
186 if (fitParams->active[i])
187 fitParams->parameters[i] =
x->data[j++];
189 fitParams->fit1D->function(fitParams->parameters, l_forSimplexLSwrap, fitParams->X, fitParams->n);
193 for (
size_t i = 0; i < fitParams->n; i++)
194 l_forSimplexLSwrap[i] = (l_forSimplexLSwrap[i] - fitParams->Y[i]) / fitParams->sigmaData[i];
198 for (
unsigned int i = 0; i < fitParams->n; i++)
199 retVal += l_forSimplexLSwrap[i] * l_forSimplexLSwrap[i];
218void Fit1D::functionDeriv(
const double *in,
Jacobian *out,
const double *xValues,
const size_t nData) {
235void Fit1D::modifyInitialFittedParameters(std::vector<double> &fittedParameter) {
236 (void)fittedParameter;
246void Fit1D::modifyFinalFittedParameters(std::vector<double> &fittedParameter) {
247 (void)fittedParameter;
254 "Name of the input Workspace");
256 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
257 mustBePositive->setLower(0);
259 "The Workspace to fit, uses the workspace numbering of the "
260 "spectra (default 0)");
262 "A value of x in, or on the low x "
263 "boundary of, the first bin to "
265 "the fit (default lowest value of x)");
267 "A value in, or on the high x boundary "
268 "of, the last bin the fitting range\n"
269 "(default the highest value of x)");
278 for (
size_t i = i0; i < props.size(); i++) {
283 "A list of comma separated parameter names which "
284 "should be fixed in the fit");
286 "Stop after this number of iterations if a good fit is not found");
291 gsl_set_error_handler_off();
296 "If not empty OutputParameters TableWorksace "
297 "and OutputWorkspace will be created.");
310 bool isDerivDefined =
true;
311 gsl_matrix *M =
nullptr;
315 const double xValuesTest = 0;
323 isDerivDefined =
false;
329 const int maxInterations =
getProperty(
"MaxIterations");
335 const size_t numberOfSpectra = localworkspace->getNumberHistograms();
337 if (histNumber >=
static_cast<int>(numberOfSpectra)) {
338 g_log.
warning(
"Invalid Workspace index given, using first Workspace");
343 const MantidVec &XValues = localworkspace->readX(histNumber);
344 const MantidVec &YValues = localworkspace->readY(histNumber);
345 const MantidVec &YErrors = localworkspace->readE(histNumber);
352 startX = XValues.front();
357 endX = XValues.back();
366 if (startX < XValues.front()) {
367 g_log.
warning(
"StartX out of range! Set to start of frame.");
368 startX = XValues.front();
372 for (m_minX = 0; XValues[m_minX + 1] < startX; ++m_minX) {
377 if (endX >= XValues.back() || endX < startX) {
378 g_log.
warning(
"EndX out of range! Set to end of frame");
379 endX = XValues.back();
380 m_maxX =
static_cast<int>(YValues.size());
382 for (m_maxX = m_minX; XValues[m_maxX] < endX; ++m_maxX) {
395 l_data.
n = m_maxX - m_minX;
398 throw std::runtime_error(
"The data set is empty.");
400 if (l_data.
n < l_data.
p) {
401 g_log.
error(
"Number of data points less than number of parameters to be fitted.");
402 throw std::runtime_error(
"Number of data points less than number of parameters to be fitted.");
404 l_data.
X =
new double[l_data.
n];
411 const bool isHistogram = localworkspace->isHistogramData();
412 for (
unsigned int i = 0; i < l_data.
n; ++i) {
414 l_data.
X[i] = 0.5 * (XValues[m_minX + i] + XValues[m_minX + i + 1]);
416 l_data.
X[i] = XValues[m_minX + i];
419 l_data.
Y = &YValues[m_minX];
422 for (
unsigned int i = 0; i < l_data.
n; ++i) {
423 if (YErrors[m_minX + i] <= 0.0) {
426 l_data.
sigmaData[i] = YErrors[m_minX + i];
442 for (
size_t i = 0; i <
nParams(); i++) {
446 for (
size_t i = 0; i <
nParams(); i++) {
452 gsl_vector *initFuncArg;
453 initFuncArg = gsl_vector_alloc(l_data.
p);
455 for (
size_t i = 0, j = 0; i <
nParams(); i++) {
462 gsl_multimin_function gslSimplexContainer;
463 gslSimplexContainer.n = l_data.
p;
465 gslSimplexContainer.params = &l_data;
469 gsl_multifit_function_fdf f;
479 const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
480 gsl_multifit_fdfsolver *s =
nullptr;
481 if (isDerivDefined) {
482 s = gsl_multifit_fdfsolver_alloc(T, l_data.
n, l_data.
p);
483 gsl_multifit_fdfsolver_set(s, &f, initFuncArg);
488 const gsl_multimin_fminimizer_type *simplexType = gsl_multimin_fminimizer_nmsimplex;
489 gsl_multimin_fminimizer *simplexMinimizer =
nullptr;
490 gsl_vector *simplexStepSize =
nullptr;
491 if (!isDerivDefined) {
492 simplexMinimizer = gsl_multimin_fminimizer_alloc(simplexType, l_data.
p);
493 simplexStepSize = gsl_vector_alloc(l_data.
p);
494 gsl_vector_set_all(simplexStepSize,
496 gsl_multimin_fminimizer_set(simplexMinimizer, &gslSimplexContainer, initFuncArg, simplexStepSize);
503 double finalCostFuncVal;
504 auto dof =
static_cast<double>(l_data.
n - l_data.
p);
508 Progress prog(
this, 0.0, 1.0, maxInterations);
509 if (isDerivDefined) {
513 status = gsl_multifit_fdfsolver_iterate(s);
518 status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
520 }
while (status == GSL_CONTINUE && iter < maxInterations);
522 double chi = gsl_blas_dnrm2(s->f);
523 finalCostFuncVal = chi * chi / dof;
526 for (
size_t i = 0, j = 0; i <
nParams(); i++)
532 status = gsl_multimin_fminimizer_iterate(simplexMinimizer);
537 double size = gsl_multimin_fminimizer_size(simplexMinimizer);
538 status = gsl_multimin_test_size(size, 1e-2);
540 }
while (status == GSL_CONTINUE && iter < maxInterations);
542 finalCostFuncVal = simplexMinimizer->fval / dof;
554 std::string reportOfFit = gsl_strerror(status);
557 <<
"Status = " << reportOfFit <<
"\n"
558 <<
"Chi^2/DoF = " << finalCostFuncVal <<
"\n";
565 setProperty(
"OutputChi2overDoF", finalCostFuncVal);
570 if (!output.empty()) {
573 gsl_matrix *covar(
nullptr);
574 std::vector<double> standardDeviations;
575 std::vector<double> sdExtended;
576 if (isDerivDefined) {
577 covar = gsl_matrix_alloc(l_data.
p, l_data.
p);
578#if GSL_MAJOR_VERSION < 2
579 gsl_multifit_covar(s->J, 0.0, covar);
581 gsl_matrix *J = gsl_matrix_alloc(l_data.
n, l_data.
p);
582 gsl_multifit_fdfsolver_jac(s, J);
583 gsl_multifit_covar(J, 0.0, covar);
588 for (
size_t i = 0; i <
nParams(); i++) {
589 sdExtended.emplace_back(1.0);
591 sdExtended[i] = sqrt(gsl_matrix_get(covar, iPNotFixed, iPNotFixed));
596 for (
size_t i = 0; i <
nParams(); i++)
598 standardDeviations.emplace_back(sdExtended[i]);
602 "The name of the TableWorkspace in which to store the final "
603 "covariance matrix");
604 setPropertyValue(
"OutputNormalisedCovarianceMatrix", output +
"_NormalisedCovarianceMatrix");
608 m_covariance->addColumn(
"str",
"Name");
609 std::vector<std::string> paramThatAreFitted;
610 for (
size_t i = 0; i <
nParams(); i++) {
617 for (
size_t i = 0; i < l_data.
p; i++) {
620 row << paramThatAreFitted[i];
621 for (
size_t j = 0; j < l_data.
p; j++) {
625 row << 100.0 * gsl_matrix_get(covar, i, j) /
626 sqrt(gsl_matrix_get(covar, i, i) * gsl_matrix_get(covar, j, j));
631 setProperty(
"OutputNormalisedCovarianceMatrix", m_covariance);
636 "The name of the TableWorkspace in which to store the "
637 "final fit parameters");
639 "Name of the output Workspace holding resulting simlated spectrum");
647 m_result->addColumn(
"str",
"Name");
648 m_result->addColumn(
"double",
"Value");
650 m_result->addColumn(
"double",
"Error");
652 firstRow <<
"Chi^2/DoF" << finalCostFuncVal;
654 for (
size_t i = 0; i <
nParams(); i++) {
657 if (isDerivDefined && l_data.
active[i]) {
659 row << sdExtended[i];
667 const MantidVec &inputX = inputWorkspace->readX(iSpec);
668 const MantidVec &inputY = inputWorkspace->readY(iSpec);
670 int histN = isHistogram ? 1 : 0;
674 ws->getAxis(0)->unit() = inputWorkspace->getAxis(0)->unit();
676 for (
int i = 0; i < 3; i++)
677 ws->dataX(i).assign(inputX.begin() + m_minX, inputX.begin() + m_maxX + histN);
679 ws->dataY(0).assign(inputY.begin() + m_minX, inputY.begin() + m_maxX);
684 auto lOut =
new double[l_data.
n];
692 for (
unsigned int i = 0; i < l_data.
n; i++) {
694 E[i] = l_data.
Y[i] -
Y[i];
699 setProperty(
"OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(ws));
702 gsl_matrix_free(covar);
708 gsl_multifit_fdfsolver_free(s);
710 gsl_vector_free(simplexStepSize);
711 gsl_multimin_fminimizer_free(simplexMinimizer);
718 gsl_vector_free(initFuncArg);
725FitData::FitData(
Fit1D *fit,
const std::string &fixed)
726 :
n(0),
X(nullptr),
Y(nullptr), sigmaData(nullptr), fit1D(fit), forSimplexLSwrap(nullptr), parameters(nullptr) {
729 for (
const auto &name : namesStrTok) {
730 std::vector<std::string>::const_iterator i =
735 throw std::invalid_argument(
"Attempt to fix non-existing parameter " + name);
738 for (
int i = 0; i < int(
active.size()); i++)
740 J.
m_map[i] =
static_cast<int>(
p++);
double value
The value of the point.
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
const std::vector< Kernel::Property * > & getProperties() const override
Get the list of managed properties.
Represents the Jacobian in IFitFunction::functionDeriv.
Helper class for reporting progress from algorithms.
TableRow represents a row in a TableWorkspace.
A property class for workspaces.
Deprecation notice: instead of using this algorithm please use the Fit algorithm instead.
virtual void declareAdditionalProperties()
Declare additional properties other than fitting parameters.
std::vector< std::string > m_parameterNames
Holds a copy of the names of the fitting parameters.
std::vector< double > m_fittedParameter
Holds a copy of the value of the parameters that are actually least-squared fitted.
virtual void afterDataRangedDetermined(const int &m_minX, const int &m_maxX)
Called after the data ranged has been determined but before the fitting starts.
virtual void function(const double *in, double *out, const double *xValues, const size_t nData)=0
Function you want to fit to.
virtual void functionDeriv(const double *in, API::Jacobian *out, const double *xValues, const size_t nData)
Derivatives of function with respect to parameters you are trying to fit.
size_t nParams() const
Number of parameters (incuding fixed).
virtual void prepare()
Called in the beginning of exec(). Custom initialization.
virtual void modifyInitialFittedParameters(std::vector< double > &fittedParameter)
Overload this function if the actual fitted parameters are different from those the user specifies.
virtual void modifyStartOfRange(double &startX)
Option for providing intelligent range starting value based e.g.
virtual void modifyFinalFittedParameters(std::vector< double > &fittedParameter)
If modifyInitialFittedParameters is overloaded this method must also be overloaded to reverse the eff...
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
virtual void modifyEndOfRange(double &endX)
Option for providing intelligent range finishing value based e.g.
virtual void declareParameters()=0
Declare parameters specific to fitting function.
The implementation of Jacobian.
std::map< int, int > m_map
The index map.
void setJ(gsl_matrix *J)
Set the pointer to the GSL's jacobian.
JacobianImpl()
Default constructor.
void set(size_t iY, size_t iP, double value) override
Set a value to a Jacobian matrix element.
void zero() override
Zero all matrix elements.
double get(size_t iY, size_t iP) override
Get a value to a Jacobian matrix element.
gsl_matrix * m_J
The pointer to the GSL's internal jacobian matrix.
Marks code as not implemented yet.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
@ TOK_TRIM
remove leading and trailing whitespace from tokens
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
static double gsl_costFunction(const gsl_vector *v, void *params)
The gsl_costFunction is optimized by GSL simplex.
int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
Fit derivatives and function GSL wrapper.
int gsl_f(const gsl_vector *x, void *params, gsl_vector *f)
Fit GSL function wrapper.
int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J)
Fit GSL derivative function wrapper.
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
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::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
double gsl_costFunction(const gsl_vector *x, void *params)
Structure to contain least squares data and used by GSL.
JacobianImpl J
Jacobi matrix interface.
std::vector< bool > active
Holds a boolean for each parameter, true if it's active or false if it's fixed.
Fit1D * fit1D
pointer to instance of Fit1D
size_t n
number of points to be fitted (size of X, Y and sigmaData arrays)
double * X
the data to be fitted (abscissae)
size_t p
number of (active) fit parameters
double * forSimplexLSwrap
Needed to use the simplex algorithm within the gsl least-squared framework.
const double * Y
the data to be fitted (ordinates)
double * parameters
A copy of the parameters.
double * sigmaData
the standard deviations of the Y data points
@ Input
An input workspace.
@ Output
An output workspace.