Mantid
Loading...
Searching...
No Matches
Fit1D.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/Axis.h"
13#include "MantidAPI/TableRow.h"
20
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>
26
27#include <cmath>
28#include <numeric>
29#include <sstream>
30
32
33using namespace Kernel;
34using API::Jacobian;
35using API::MatrixWorkspace;
37using API::Progress;
38using API::WorkspaceProperty;
39
41class JacobianImpl : public Jacobian {
42public:
44 JacobianImpl() : Jacobian(), m_J(nullptr){};
45
47 std::map<int, int> m_map;
54 void set(size_t iY, size_t iP, double value) override {
55 int j = m_map[static_cast<int>(iP)];
56 if (j >= 0)
57 gsl_matrix_set(m_J, iY, j, value);
58 }
64 double get(size_t iY, size_t iP) override {
65 int j = m_map[static_cast<int>(iP)];
66 if (j >= 0)
67 return gsl_matrix_get(m_J, iY, j);
68 return 0.0;
69 }
72 void zero() override { gsl_matrix_set_zero(m_J); }
74 void setJ(gsl_matrix *J) { m_J = J; }
75
76private:
78 gsl_matrix *m_J;
79};
80
82struct FitData {
84 FitData(Fit1D *fit, const std::string &fixed);
86 size_t n;
88 size_t p;
90 double *X;
92 const double *Y;
94 double *sigmaData;
102 double *parameters;
105 std::vector<bool> active;
108};
109
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])
120 fitParams->parameters[i] = x->data[j++];
121
122 fitParams->fit1D->function(fitParams->parameters, f->data, fitParams->X, fitParams->n);
123
124 // function() return calculated data values. Need to convert this values into
125 // calculated-observed devided by error values used by GSL
126
127 for (size_t i = 0; i < fitParams->n; i++)
128 f->data[i] = (f->data[i] - fitParams->Y[i]) / fitParams->sigmaData[i];
129
130 return GSL_SUCCESS;
131}
132
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])
143 fitParams->parameters[i] = x->data[j++];
144 }
145
146 fitParams->J.setJ(J);
147
148 fitParams->fit1D->functionDeriv(fitParams->parameters, &fitParams->J, fitParams->X, fitParams->n);
149
150 // functionDeriv() return derivatives of calculated data values. Need to
151 // convert this values into derivatives of calculated-observed devided
152 // by error values used by GSL
153
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];
157
158 return GSL_SUCCESS;
159}
160
168static int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J) {
169 gsl_f(x, params, f);
170 gsl_df(x, params, J);
171 return GSL_SUCCESS;
172}
173
181static double gsl_costFunction(const gsl_vector *x, void *params) {
182 auto fitParams = static_cast<FitData *>(params);
183 double *l_forSimplexLSwrap = fitParams->forSimplexLSwrap;
184
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++];
188
189 fitParams->fit1D->function(fitParams->parameters, l_forSimplexLSwrap, fitParams->X, fitParams->n);
190
191 // function() return calculated data values. Need to convert this values into
192 // calculated-observed devided by error values used by GSL
193 for (size_t i = 0; i < fitParams->n; i++)
194 l_forSimplexLSwrap[i] = (l_forSimplexLSwrap[i] - fitParams->Y[i]) / fitParams->sigmaData[i];
195
196 double retVal = 0.0;
197
198 for (unsigned int i = 0; i < fitParams->n; i++)
199 retVal += l_forSimplexLSwrap[i] * l_forSimplexLSwrap[i];
200
201 return retVal;
202}
203
218void Fit1D::functionDeriv(const double *in, Jacobian *out, const double *xValues, const size_t nData) {
219 UNUSED_ARG(in);
220 UNUSED_ARG(out);
221 UNUSED_ARG(xValues);
222 UNUSED_ARG(nData);
223 throw Exception::NotImplementedError("No derivative function provided");
224}
225
235void Fit1D::modifyInitialFittedParameters(std::vector<double> &fittedParameter) {
236 (void)fittedParameter; // Avoid compiler warning
237}
238
246void Fit1D::modifyFinalFittedParameters(std::vector<double> &fittedParameter) {
247 (void)fittedParameter; // Avoid compiler warning
248}
249
252void Fit1D::init() {
253 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
254 "Name of the input Workspace");
255
256 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
257 mustBePositive->setLower(0);
258 declareProperty("WorkspaceIndex", 0, mustBePositive,
259 "The Workspace to fit, uses the workspace numbering of the "
260 "spectra (default 0)");
261 declareProperty("StartX", EMPTY_DBL(),
262 "A value of x in, or on the low x "
263 "boundary of, the first bin to "
264 "include in\n"
265 "the fit (default lowest value of x)");
266 declareProperty("EndX", EMPTY_DBL(),
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)");
270
271 size_t i0 = getProperties().size();
272
273 // declare parameters specific to a given fitting function
275
276 // load the name of these specific parameter into a vector for later use
277 const std::vector<Property *> props = getProperties();
278 for (size_t i = i0; i < props.size(); i++) {
279 m_parameterNames.emplace_back(props[i]->name());
280 }
281
282 declareProperty("Fix", "",
283 "A list of comma separated parameter names which "
284 "should be fixed in the fit");
285 declareProperty("MaxIterations", 500, mustBePositive,
286 "Stop after this number of iterations if a good fit is not found");
287 declareProperty("OutputStatus", "", Direction::Output);
288 declareProperty("OutputChi2overDoF", 0.0, Direction::Output);
289
290 // Disable default gsl error handler (which is to call abort!)
291 gsl_set_error_handler_off();
292
294
295 declareProperty("Output", "",
296 "If not empty OutputParameters TableWorksace "
297 "and OutputWorkspace will be created.");
298}
299
304void Fit1D::exec() {
305
306 // Custom initialization
307 prepare();
308
309 // check if derivative defined in derived class
310 bool isDerivDefined = true;
311 gsl_matrix *M = nullptr;
312 try {
313 const std::vector<double> inTest(m_parameterNames.size(), 1.0);
314 std::vector<double> outTest(m_parameterNames.size());
315 const double xValuesTest = 0;
316 JacobianImpl J;
317 M = gsl_matrix_alloc(m_parameterNames.size(), 1);
318 J.setJ(M);
319 // note nData set to zero (last argument) hence this should avoid further
320 // memory problems
321 functionDeriv(&(inTest.front()), &J, &xValuesTest, 0);
323 isDerivDefined = false;
324 }
325 gsl_matrix_free(M);
326
327 // Try to retrieve optional properties
328 int histNumber = getProperty("WorkspaceIndex");
329 const int maxInterations = getProperty("MaxIterations");
330
331 // Get the input workspace
332 MatrixWorkspace_const_sptr localworkspace = getProperty("InputWorkspace");
333
334 // number of histogram is equal to the number of spectra
335 const size_t numberOfSpectra = localworkspace->getNumberHistograms();
336 // Check that the index given is valid
337 if (histNumber >= static_cast<int>(numberOfSpectra)) {
338 g_log.warning("Invalid Workspace index given, using first Workspace");
339 histNumber = 0;
340 }
341
342 // Retrieve the spectrum into a vector
343 const MantidVec &XValues = localworkspace->readX(histNumber);
344 const MantidVec &YValues = localworkspace->readY(histNumber);
345 const MantidVec &YErrors = localworkspace->readE(histNumber);
346
347 // Read in the fitting range data that we were sent
348 double startX = getProperty("StartX");
349 double endX = getProperty("EndX");
350 // check if the values had been set, otherwise use defaults
351 if (isEmpty(startX)) {
352 startX = XValues.front();
353 modifyStartOfRange(startX); // does nothing by default but derived class may
354 // provide a more intelligent value
355 }
356 if (isEmpty(endX)) {
357 endX = XValues.back();
358 modifyEndOfRange(endX); // does nothing by default but derived class may
359 // previde a more intelligent value
360 }
361
362 int m_minX;
363 int m_maxX;
364
365 // Check the validity of startX
366 if (startX < XValues.front()) {
367 g_log.warning("StartX out of range! Set to start of frame.");
368 startX = XValues.front();
369 }
370 // Get the corresponding bin boundary that comes before (or coincides with)
371 // this value
372 for (m_minX = 0; XValues[m_minX + 1] < startX; ++m_minX) {
373 }
374
375 // Check the validity of endX and get the bin boundary that come after (or
376 // coincides with) it
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());
381 } else {
382 for (m_maxX = m_minX; XValues[m_maxX] < endX; ++m_maxX) {
383 }
384 }
385
386 afterDataRangedDetermined(m_minX, m_maxX);
387
388 // create and populate GSL data container warn user if l_data.n < l_data.p
389 // since as a rule of thumb this is required as a minimum to obtained
390 // 'accurate'
391 // fitting parameter values.
392
393 FitData l_data(this, getProperty("Fix"));
394
395 l_data.n = m_maxX - m_minX; // m_minX and m_maxX are array index markers. I.e. e.g. 0 & 19.
396 if (l_data.n == 0) {
397 g_log.error("The data set is empty.");
398 throw std::runtime_error("The data set is empty.");
399 }
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.");
403 }
404 l_data.X = new double[l_data.n];
405 l_data.sigmaData = new double[l_data.n];
406 l_data.forSimplexLSwrap = new double[l_data.n];
407 l_data.parameters = new double[nParams()];
408
409 // check if histogram data in which case use mid points of histogram bins
410
411 const bool isHistogram = localworkspace->isHistogramData();
412 for (unsigned int i = 0; i < l_data.n; ++i) {
413 if (isHistogram)
414 l_data.X[i] = 0.5 * (XValues[m_minX + i] + XValues[m_minX + i + 1]); // take mid-point if histogram bin
415 else
416 l_data.X[i] = XValues[m_minX + i];
417 }
418
419 l_data.Y = &YValues[m_minX];
420
421 // check that no error is negative or zero
422 for (unsigned int i = 0; i < l_data.n; ++i) {
423 if (YErrors[m_minX + i] <= 0.0) {
424 l_data.sigmaData[i] = 1.0;
425 } else
426 l_data.sigmaData[i] = YErrors[m_minX + i];
427 }
428
429 // create array of fitted parameter. Take these to those input by the user.
430 // However, for doing the
431 // underlying fitting it might be more efficient to actually perform the
432 // fitting on some of other
433 // form of the fitted parameters. For instance, take the Gaussian sigma
434 // parameter. In practice it
435 // in fact more efficient to perform the fitting not on sigma but 1/sigma^2.
436 // The methods
437 // modifyInitialFittedParameters() and modifyFinalFittedParameters() are used
438 // to allow for this;
439 // by default these function do nothing.
440
441 m_fittedParameter.clear();
442 for (size_t i = 0; i < nParams(); i++) {
444 }
445 modifyInitialFittedParameters(m_fittedParameter); // does nothing except if overwritten by derived class
446 for (size_t i = 0; i < nParams(); i++) {
447 l_data.parameters[i] = m_fittedParameter[i];
448 }
449
450 // set-up initial guess for fit parameters
451
452 gsl_vector *initFuncArg;
453 initFuncArg = gsl_vector_alloc(l_data.p);
454
455 for (size_t i = 0, j = 0; i < nParams(); i++) {
456 if (l_data.active[i])
457 gsl_vector_set(initFuncArg, j++, m_fittedParameter[i]);
458 }
459
460 // set-up GSL container to be used with GSL simplex algorithm
461
462 gsl_multimin_function gslSimplexContainer;
463 gslSimplexContainer.n = l_data.p; // n here refers to number of parameters
464 gslSimplexContainer.f = &gsl_costFunction;
465 gslSimplexContainer.params = &l_data;
466
467 // set-up GSL least squares container
468
469 gsl_multifit_function_fdf f;
470 f.f = &gsl_f;
471 f.df = &gsl_df;
472 f.fdf = &gsl_fdf;
473 f.n = l_data.n;
474 f.p = l_data.p;
475 f.params = &l_data;
476
477 // set-up remaining GSL machinery for least squared
478
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);
484 }
485
486 // set-up remaining GSL machinery to use simplex algorithm
487
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,
495 1.0); // is this always a sensible starting step size?
496 gsl_multimin_fminimizer_set(simplexMinimizer, &gslSimplexContainer, initFuncArg, simplexStepSize);
497 }
498
499 // finally do the fitting
500
501 int iter = 0;
502 int status;
503 double finalCostFuncVal;
504 auto dof = static_cast<double>(l_data.n - l_data.p); // dof stands for degrees of freedom
505
506 // Standard least-squares used if derivative function defined otherwise
507 // simplex
508 Progress prog(this, 0.0, 1.0, maxInterations);
509 if (isDerivDefined) {
510
511 do {
512 iter++;
513 status = gsl_multifit_fdfsolver_iterate(s);
514
515 if (status) // break if error
516 break;
517
518 status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
519 prog.report();
520 } while (status == GSL_CONTINUE && iter < maxInterations);
521
522 double chi = gsl_blas_dnrm2(s->f);
523 finalCostFuncVal = chi * chi / dof;
524
525 // put final converged fitting values back into m_fittedParameter
526 for (size_t i = 0, j = 0; i < nParams(); i++)
527 if (l_data.active[i])
528 m_fittedParameter[i] = gsl_vector_get(s->x, j++);
529 } else {
530 do {
531 iter++;
532 status = gsl_multimin_fminimizer_iterate(simplexMinimizer);
533
534 if (status) // break if error
535 break;
536
537 double size = gsl_multimin_fminimizer_size(simplexMinimizer);
538 status = gsl_multimin_test_size(size, 1e-2);
539 prog.report();
540 } while (status == GSL_CONTINUE && iter < maxInterations);
541
542 finalCostFuncVal = simplexMinimizer->fval / dof;
543
544 // put final converged fitting values back into m_fittedParameter
545 for (unsigned int i = 0, j = 0; i < m_fittedParameter.size(); i++)
546 if (l_data.active[i])
547 m_fittedParameter[i] = gsl_vector_get(simplexMinimizer->x, j++);
548 }
549
550 modifyFinalFittedParameters(m_fittedParameter); // do nothing except if overwritten by derived class
551
552 // Output summary to log file
553
554 std::string reportOfFit = gsl_strerror(status);
555
556 g_log.information() << "Iteration = " << iter << "\n"
557 << "Status = " << reportOfFit << "\n"
558 << "Chi^2/DoF = " << finalCostFuncVal << "\n";
559 for (size_t i = 0; i < m_fittedParameter.size(); i++)
560 g_log.information() << m_parameterNames[i] << " = " << m_fittedParameter[i] << " \n";
561
562 // also output summary to properties
563
564 setProperty("OutputStatus", reportOfFit);
565 setProperty("OutputChi2overDoF", finalCostFuncVal);
566 for (size_t i = 0; i < m_fittedParameter.size(); i++)
568
569 std::string output = getProperty("Output");
570 if (!output.empty()) {
571 // calculate covariance matrix if derivatives available
572
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);
580#else
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);
584 gsl_matrix_free(J);
585#endif
586
587 int iPNotFixed = 0;
588 for (size_t i = 0; i < nParams(); i++) {
589 sdExtended.emplace_back(1.0);
590 if (l_data.active[i]) {
591 sdExtended[i] = sqrt(gsl_matrix_get(covar, iPNotFixed, iPNotFixed));
592 iPNotFixed++;
593 }
594 }
595 modifyFinalFittedParameters(sdExtended);
596 for (size_t i = 0; i < nParams(); i++)
597 if (l_data.active[i])
598 standardDeviations.emplace_back(sdExtended[i]);
599
600 declareProperty(std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("OutputNormalisedCovarianceMatrix", "",
602 "The name of the TableWorkspace in which to store the final "
603 "covariance matrix");
604 setPropertyValue("OutputNormalisedCovarianceMatrix", output + "_NormalisedCovarianceMatrix");
605
607 Mantid::API::WorkspaceFactory::Instance().createTable("TableWorkspace");
608 m_covariance->addColumn("str", "Name");
609 std::vector<std::string> paramThatAreFitted; // used for populating 1st "name" column
610 for (size_t i = 0; i < nParams(); i++) {
611 if (l_data.active[i]) {
612 m_covariance->addColumn("double", m_parameterNames[i]);
613 paramThatAreFitted.emplace_back(m_parameterNames[i]);
614 }
615 }
616
617 for (size_t i = 0; i < l_data.p; i++) {
618
619 Mantid::API::TableRow row = m_covariance->appendRow();
620 row << paramThatAreFitted[i];
621 for (size_t j = 0; j < l_data.p; j++) {
622 if (j == i)
623 row << 1.0;
624 else {
625 row << 100.0 * gsl_matrix_get(covar, i, j) /
626 sqrt(gsl_matrix_get(covar, i, i) * gsl_matrix_get(covar, j, j));
627 }
628 }
629 }
630
631 setProperty("OutputNormalisedCovarianceMatrix", m_covariance);
632 }
633
635 std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("OutputParameters", "", Direction::Output),
636 "The name of the TableWorkspace in which to store the "
637 "final fit parameters");
638 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
639 "Name of the output Workspace holding resulting simlated spectrum");
640
641 setPropertyValue("OutputParameters", output + "_Parameters");
642 setPropertyValue("OutputWorkspace", output + "_Workspace");
643
644 // Save the final fit parameters in the output table workspace
646 Mantid::API::WorkspaceFactory::Instance().createTable("TableWorkspace");
647 m_result->addColumn("str", "Name");
648 m_result->addColumn("double", "Value");
649 if (isDerivDefined)
650 m_result->addColumn("double", "Error");
651 Mantid::API::TableRow firstRow = m_result->appendRow();
652 firstRow << "Chi^2/DoF" << finalCostFuncVal;
653
654 for (size_t i = 0; i < nParams(); i++) {
655 Mantid::API::TableRow row = m_result->appendRow();
656 row << m_parameterNames[i] << m_fittedParameter[i];
657 if (isDerivDefined && l_data.active[i]) {
658 // perhaps want to scale standard deviations with sqrt(finalCostFuncVal)
659 row << sdExtended[i];
660 }
661 }
662 setProperty("OutputParameters", m_result);
663
664 // Save the fitted and simulated spectra in the output workspace
665 MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
666 int iSpec = getProperty("WorkspaceIndex");
667 const MantidVec &inputX = inputWorkspace->readX(iSpec);
668 const MantidVec &inputY = inputWorkspace->readY(iSpec);
669
670 int histN = isHistogram ? 1 : 0;
671 Mantid::DataObjects::Workspace2D_sptr ws = std::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>(
672 Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", 3, l_data.n + histN, l_data.n));
673 ws->setTitle("");
674 ws->getAxis(0)->unit() = inputWorkspace->getAxis(0)->unit(); // UnitFactory::Instance().create("TOF");
675
676 for (int i = 0; i < 3; i++)
677 ws->dataX(i).assign(inputX.begin() + m_minX, inputX.begin() + m_maxX + histN);
678
679 ws->dataY(0).assign(inputY.begin() + m_minX, inputY.begin() + m_maxX);
680
681 MantidVec &Y = ws->dataY(1);
682 MantidVec &E = ws->dataY(2);
683
684 auto lOut = new double[l_data.n]; // to capture output from call to function()
685 modifyInitialFittedParameters(m_fittedParameter); // does nothing except if
686 // overwritten by derived
687 // class
688 function(&m_fittedParameter[0], lOut, l_data.X, l_data.n);
689 modifyInitialFittedParameters(m_fittedParameter); // reverse the effect of
690 // modifyInitialFittedParameters - if any
691
692 for (unsigned int i = 0; i < l_data.n; i++) {
693 Y[i] = lOut[i];
694 E[i] = l_data.Y[i] - Y[i];
695 }
696
697 delete[] lOut;
698
699 setProperty("OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(ws));
700
701 if (isDerivDefined)
702 gsl_matrix_free(covar);
703 }
704
705 // clean up dynamically allocated gsl stuff
706
707 if (isDerivDefined)
708 gsl_multifit_fdfsolver_free(s);
709 else {
710 gsl_vector_free(simplexStepSize);
711 gsl_multimin_fminimizer_free(simplexMinimizer);
712 }
713
714 delete[] l_data.X;
715 delete[] l_data.sigmaData;
716 delete[] l_data.forSimplexLSwrap;
717 delete[] l_data.parameters;
718 gsl_vector_free(initFuncArg);
719}
720
725FitData::FitData(Fit1D *fit, const std::string &fixed)
726 : n(0), X(nullptr), Y(nullptr), sigmaData(nullptr), fit1D(fit), forSimplexLSwrap(nullptr), parameters(nullptr) {
728 active.insert(active.begin(), fit1D->m_parameterNames.size(), true);
729 for (const auto &name : namesStrTok) {
730 std::vector<std::string>::const_iterator i =
731 std::find(fit1D->m_parameterNames.begin(), fit1D->m_parameterNames.end(), name);
732 if (i != fit1D->m_parameterNames.end()) {
733 active[i - fit1D->m_parameterNames.begin()] = false;
734 } else
735 throw std::invalid_argument("Attempt to fix non-existing parameter " + name);
736 }
737 p = 0;
738 for (int i = 0; i < int(active.size()); i++)
739 if (active[i])
740 J.m_map[i] = static_cast<int>(p++);
741 else
742 J.m_map[i] = -1;
743}
744
745} // namespace Mantid::CurveFitting::Algorithms
double value
The value of the point.
Definition: FitMW.cpp:51
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
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
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
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
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.
Definition: Algorithm.cpp:2050
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
Deprecation notice: instead of using this algorithm please use the Fit algorithm instead.
Definition: Fit1D.h:46
virtual void declareAdditionalProperties()
Declare additional properties other than fitting parameters.
Definition: Fit1D.h:78
std::vector< std::string > m_parameterNames
Holds a copy of the names of the fitting parameters.
Definition: Fit1D.h:114
std::vector< double > m_fittedParameter
Holds a copy of the value of the parameters that are actually least-squared fitted.
Definition: Fit1D.h:111
virtual void afterDataRangedDetermined(const int &m_minX, const int &m_maxX)
Called after the data ranged has been determined but before the fitting starts.
Definition: Fit1D.h:91
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.
Definition: Fit1D.cpp:218
size_t nParams() const
Number of parameters (incuding fixed).
Definition: Fit1D.h:117
virtual void prepare()
Called in the beginning of exec(). Custom initialization.
Definition: Fit1D.h:80
virtual void modifyInitialFittedParameters(std::vector< double > &fittedParameter)
Overload this function if the actual fitted parameters are different from those the user specifies.
Definition: Fit1D.cpp:235
virtual void modifyStartOfRange(double &startX)
Option for providing intelligent range starting value based e.g.
Definition: Fit1D.h:68
virtual void modifyFinalFittedParameters(std::vector< double > &fittedParameter)
If modifyInitialFittedParameters is overloaded this method must also be overloaded to reverse the eff...
Definition: Fit1D.cpp:246
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
Definition: Fit1D.h:49
virtual void modifyEndOfRange(double &endX)
Option for providing intelligent range finishing value based e.g.
Definition: Fit1D.h:73
virtual void declareParameters()=0
Declare parameters specific to fitting function.
The implementation of Jacobian.
Definition: Fit1D.cpp:41
std::map< int, int > m_map
The index map.
Definition: Fit1D.cpp:47
void setJ(gsl_matrix *J)
Set the pointer to the GSL's jacobian.
Definition: Fit1D.cpp:74
void set(size_t iY, size_t iP, double value) override
Set a value to a Jacobian matrix element.
Definition: Fit1D.cpp:54
void zero() override
Zero all matrix elements.
Definition: Fit1D.cpp:72
double get(size_t iY, size_t iP) override
Get a value to a Jacobian matrix element.
Definition: Fit1D.cpp:64
gsl_matrix * m_J
The pointer to the GSL's internal jacobian matrix.
Definition: Fit1D.cpp:78
Marks code as not implemented yet.
Definition: Exception.h:138
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.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
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
Definition: cow_ptr.h:172
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
double gsl_costFunction(const gsl_vector *x, void *params)
Structure to contain least squares data and used by GSL.
Definition: Fit1D.cpp:82
JacobianImpl J
Jacobi matrix interface.
Definition: Fit1D.cpp:107
std::vector< bool > active
Holds a boolean for each parameter, true if it's active or false if it's fixed.
Definition: Fit1D.cpp:105
Fit1D * fit1D
pointer to instance of Fit1D
Definition: Fit1D.cpp:96
size_t n
number of points to be fitted (size of X, Y and sigmaData arrays)
Definition: Fit1D.cpp:86
double * X
the data to be fitted (abscissae)
Definition: Fit1D.cpp:90
size_t p
number of (active) fit parameters
Definition: Fit1D.cpp:88
double * forSimplexLSwrap
Needed to use the simplex algorithm within the gsl least-squared framework.
Definition: Fit1D.cpp:100
const double * Y
the data to be fitted (ordinates)
Definition: Fit1D.cpp:92
double * parameters
A copy of the parameters.
Definition: Fit1D.cpp:102
double * sigmaData
the standard deviations of the Y data points
Definition: Fit1D.cpp:94
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54