Mantid
Loading...
Searching...
No Matches
FitPeak.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//----------------------------------------------------------------------
10#include <algorithm>
11#include <iterator>
12#include <utility>
13
21#include "MantidAPI/TableRow.h"
27#include "MantidHistogramData/Histogram.h"
28#include "MantidHistogramData/HistogramBuilder.h"
34
35#include "boost/algorithm/string.hpp"
36#include "boost/algorithm/string/trim.hpp"
37
38using namespace Mantid;
39using namespace Mantid::API;
40using namespace Mantid::DataObjects;
41using namespace Mantid::HistogramData;
42using namespace Mantid::Kernel;
43using Mantid::HistogramData::HistogramX;
44
45using namespace std;
46
47const double MAGICNUMBER = 2.0;
48
49namespace Mantid::Algorithms {
50//----------------------------------------------------------------------------------------------
54 : API::Algorithm(), m_fitMethodSet(false), m_peakRangeSet(false), m_peakWidthSet(false), m_peakWindowSet(false),
55 m_usePeakPositionTolerance(false), m_peakFunc(), m_bkgdFunc(), m_dataWS(), m_wsIndex(0), m_minFitX(0.),
56 m_maxFitX(0.), i_minFitX(0), i_maxFitX(0), m_minPeakX(0.), m_maxPeakX(0.), i_minPeakX(0), i_maxPeakX(0),
57 m_bestPeakFunc(), m_bestBkgdFunc(), m_bkupPeakFunc(), m_bkupBkgdFunc(), m_fitErrorPeakFunc(),
58 m_fitErrorBkgdFunc(), m_minimizer("Levenberg-MarquardtMD"), m_costFunction(), m_vecFWHM(),
59 m_peakPositionTolerance(0.), m_userPeakCentre(0.), m_bestRwp(0.), m_finalGoodnessValue(0.), m_numFitCalls(0),
60 m_sstream("") {}
61
62//----------------------------------------------------------------------------------------------
65void FitOneSinglePeak::setWorskpace(const API::MatrixWorkspace_sptr &dataws, size_t wsindex) {
66 if (dataws) {
67 m_dataWS = dataws;
68 } else {
69 throw runtime_error("Input dataws is null. ");
70 }
71
72 if (wsindex < m_dataWS->getNumberHistograms()) {
73 m_wsIndex = wsindex;
74 } else {
75 throw runtime_error("Input workspace index is out of range.");
76 }
77}
78
79//----------------------------------------------------------------------------------------------
83 if (peakfunc)
84 m_peakFunc = peakfunc;
85
86 if (bkgdfunc)
87 m_bkgdFunc = bkgdfunc;
88}
89
90//----------------------------------------------------------------------------------------------
93void FitOneSinglePeak::setFitWindow(double leftwindow, double rightwindow) {
94 m_minFitX = leftwindow;
95 m_maxFitX = rightwindow;
96
97 const auto &vecX = m_dataWS->x(m_wsIndex);
98
101
102 m_peakWindowSet = true;
103}
104
105//----------------------------------------------------------------------------------------------
112void FitOneSinglePeak::setPeakRange(double xpeakleft, double xpeakright) {
113 m_minPeakX = xpeakleft;
114 m_maxPeakX = xpeakright;
115
116 const auto &vecX = m_dataWS->x(m_wsIndex);
117
120
121 m_peakRangeSet = true;
122}
123
124//----------------------------------------------------------------------------------------------
129void FitOneSinglePeak::setFittingMethod(std::string minimizer, const std::string &costfunction) {
130 m_minimizer = std::move(minimizer);
131 if (costfunction == "Chi-Square") {
132 m_costFunction = "Least squares";
133 } else if (costfunction == "Rwp") {
134 m_costFunction = "Rwp";
135 } else if (costfunction == "Least squares") {
136 m_costFunction = costfunction;
137 } else {
138 stringstream errss;
139 errss << "FitOneSinglePeak: cost function " << costfunction << " is not supported. ";
140 throw runtime_error(errss.str());
141 }
142
143 m_fitMethodSet = true;
144}
145
146//----------------------------------------------------------------------------------------------
156void FitOneSinglePeak::setupGuessedFWHM(double usrwidth, int minfwhm, int maxfwhm, int stepsize,
157 bool fitwithsteppedfwhm) {
158 m_vecFWHM.clear();
159
160 // From user specified guess value
161 if (usrwidth <= 0) {
162 // Set up default FWHM if user does not give reasonable peak width
163 m_sstream << "Client inputs user-defined peak width = " << usrwidth << "; Automatically reset to 4 as default."
164 << "\n";
165
166 if (!fitwithsteppedfwhm) {
167 fitwithsteppedfwhm = true;
168 minfwhm = 4;
169 maxfwhm = 4;
170 stepsize = 1;
171 } else {
172 if (minfwhm > 4) {
173 minfwhm = 4;
174 }
175 if (maxfwhm < minfwhm)
176 maxfwhm = 4;
177 }
178 } else {
179 m_vecFWHM.emplace_back(usrwidth);
180 m_sstream << "Add user defined FWHM = " << usrwidth << "\n";
181 }
182
183 m_peakWidthSet = true;
184
185 // From user specified minimum value to maximim value
186 if (!fitwithsteppedfwhm) {
187 if (m_vecFWHM.empty())
188 throw runtime_error("Logic error in setup guessed FWHM. ");
189 m_sstream << "No FWHM is not guessed by stepped FWHM. "
190 << "\n";
191 return;
192 }
193
194 auto &vecX = m_dataWS->x(m_wsIndex);
195
196 auto i_centre = static_cast<int>(getIndex(vecX, m_peakFunc->centre()));
197 int i_maxindex = static_cast<int>(vecX.size()) - 1;
198
199 m_sstream << "FWHM to guess. Range = " << minfwhm << ", " << maxfwhm << "; Step = " << stepsize << "\n";
200 if (stepsize == 0 || maxfwhm < minfwhm)
201 throw runtime_error("FWHM is not given right.");
202
203 for (int iwidth = minfwhm; iwidth <= maxfwhm; iwidth += stepsize) {
204 // There are 3 possible situation: peak at left edge, peak in proper range,
205 // peak at righ edge
206 int ileftside = i_centre - iwidth / 2;
207 if (ileftside < 0)
208 ileftside = 0;
209
210 int irightside = i_centre + iwidth / 2;
211 if (irightside > i_maxindex)
212 irightside = i_maxindex;
213
214 double in_fwhm = vecX[irightside] - vecX[ileftside];
215
216 if (in_fwhm < 1.0E-20) {
217 m_sstream << "It is impossible to have zero peak width as iCentre = " << i_centre << ", iWidth = " << iwidth
218 << "\n"
219 << "More information: Spectrum = " << m_wsIndex << "; Range of X is " << vecX.front() << ", "
220 << vecX.back() << "; Peak centre = " << vecX[i_centre] << "\n";
221 } else {
222 m_sstream << "Setup: i_width = " << iwidth << ", i_left = " << ileftside << ", i_right = " << irightside
223 << ", FWHM = " << in_fwhm << ", i_centre = " << i_centre << ".\n";
224 }
225
226 m_vecFWHM.emplace_back(in_fwhm);
227 }
228}
229
230//----------------------------------------------------------------------------------------------
237void FitOneSinglePeak::setFitPeakCriteria(bool usepeakpostol, double peakpostol) {
238 m_usePeakPositionTolerance = usepeakpostol;
239 if (usepeakpostol) {
240 m_peakPositionTolerance = fabs(peakpostol);
241 if (peakpostol < 1.0E-13)
242 g_log.warning("Peak position tolerance is very tight. ");
243 }
244}
245
246//----------------------------------------------------------------------------------------------
249bool FitOneSinglePeak::hasSetupToFitPeak(std::string &errmsg) {
250 errmsg = "";
251
252 if (!m_fitMethodSet)
253 errmsg += "Fitting method ";
254 if (!m_peakRangeSet)
255 errmsg += "Peak range ";
256 if (!m_peakWidthSet)
257 errmsg += "Peak width ";
258 if (!m_peakFunc)
259 errmsg += "Peak function ";
260 if (!m_bkgdFunc)
261 errmsg += "Background function ";
262 if (!m_dataWS)
263 errmsg += "Data workspace ";
264
265 if (!errmsg.empty()) {
266 errmsg = "These parameters have not been set for fitting peak: " + errmsg;
267 return false;
268 }
269
270 return true;
271}
272
273//----------------------------------------------------------------------------------------------
276std::string FitOneSinglePeak::getDebugMessage() { return m_sstream.str(); }
277
278//----------------------------------------------------------------------------------------------
282 m_numFitCalls = 0;
283 string errmsg;
284 if (!hasSetupToFitPeak(errmsg)) {
285 g_log.error(errmsg);
286 throw runtime_error("Object has not been set up completely to fit peak.");
287 }
288
289 // Initialize refinement state parameters
290 m_bestRwp = DBL_MAX;
291
292 // Set up a composite function
293 CompositeFunction_sptr compfunc = std::make_shared<CompositeFunction>();
294 compfunc->addFunction(m_peakFunc);
295 compfunc->addFunction(m_bkgdFunc);
296
297 m_sstream << "One-Step-Fit Function: " << compfunc->asString() << "\n";
298
299 // Store starting setup
302
303 // Fit with different starting values of peak width
304 size_t numfits = m_vecFWHM.size();
305
306 Progress progress(this, 0.0, 1.0, numfits);
307
308 for (size_t i = 0; i < numfits; ++i) {
309 // set FWHM
310 m_sstream << "[SingleStepFit] FWHM = " << m_vecFWHM[i] << "\n";
311 m_peakFunc->setFwhm(m_vecFWHM[i]);
312
313 // fit and process result
314 double goodndess = fitFunctionSD(compfunc, m_dataWS, m_wsIndex, m_minFitX, m_maxFitX);
315 processNStoreFitResult(goodndess, true);
316
317 // restore the function parameters
318 if (i != numfits - 1) {
321 }
322
323 progress.report();
324 }
325
326 // Retrieve the best result stored
329
331
332 m_sstream << "One-Step-Fit Best (Chi^2 = " << m_bestRwp << ") Fitted Function: " << compfunc->asString() << "\n"
333 << "Number of calls of Fit = " << m_numFitCalls << "\n";
334
335 return false;
336}
337
338//----------------------------------------------------------------------------------------------
342
343 auto &vecY = m_dataWS->y(m_wsIndex);
344
345 size_t size = i_maxFitX - i_minFitX + 1;
346 size_t ysize = size;
347 size_t ishift = i_maxFitX + 1;
348 if (ishift >= vecY.size())
349 ysize = vecY.size() - i_minFitX;
350
351 HistogramBuilder builder;
352 builder.setX(size);
353 builder.setY(ysize);
354 MatrixWorkspace_sptr purePeakWS = create<Workspace2D>(1, builder.build());
355
356 auto &vecX = m_dataWS->x(m_wsIndex);
357 auto &vecE = m_dataWS->e(m_wsIndex);
358 auto &dataX = purePeakWS->mutableX(0);
359 auto &dataY = purePeakWS->mutableY(0);
360 auto &dataE = purePeakWS->mutableE(0);
361
362 dataX.assign(vecX.cbegin() + i_minFitX, vecX.cbegin() + i_maxFitX + 1);
363 if (ishift < vecY.size()) {
364 dataY.assign(vecY.cbegin() + i_minFitX, vecY.cbegin() + i_maxFitX + 1);
365 dataE.assign(vecE.cbegin() + i_minFitX, vecE.cbegin() + i_maxFitX + 1);
366 } else {
367 dataY.assign(vecY.cbegin() + i_minFitX, vecY.cend());
368 dataE.assign(vecE.cbegin() + i_minFitX, vecE.cend());
369 }
370
371 return purePeakWS;
372}
373
374//----------------------------------------------------------------------------------------------
378 const MatrixWorkspace_sptr &dataws, size_t wsindex, size_t ixmin,
379 size_t ixmax) {
380 // Get current peak height: from current peak centre (previously setup)
381 double peakcentre = peakfunc->centre();
382 vector<double> svvec(1, peakcentre);
383 FunctionDomain1DVector svdomain(svvec);
384 FunctionValues svvalues(svdomain);
385 peakfunc->function(svdomain, svvalues);
386 double curpeakheight = svvalues[0];
387
388 const auto &vecX = dataws->x(wsindex);
389 const auto &vecY = dataws->y(wsindex);
390 double ymax = vecY[ixmin + 1];
391 size_t iymax = ixmin + 1;
392 for (size_t i = ixmin + 2; i < ixmax; ++i) {
393 double tempy = vecY[i];
394 if (tempy > ymax) {
395 ymax = tempy;
396 iymax = i;
397 }
398 }
399
400 m_sstream << "Estimate-Peak-Height: Current peak height = " << curpeakheight
401 << ". Estimate-Peak-Height: Maximum Y value between " << vecX[ixmin] << " and " << vecX[ixmax] << " is "
402 << ymax << " at X = " << vecX[iymax] << ".\n";
403
404 // Compute peak height (not the maximum peak intensity)
405 double estheight = ymax / curpeakheight * peakfunc->height();
406
407 return estheight;
408}
409
410//----------------------------------------------------------------------------------------------
415 // Calculate background
416 // FIXME - This can be costly to use FunctionDomain and FunctionValue
417 auto &vecX = purePeakWS->x(0);
418 FunctionDomain1DVector domain(MantidVec(vecX.begin(), vecX.end()));
419 FunctionValues bkgdvalues(domain);
420 m_bkgdFunc->function(domain, bkgdvalues);
421
422 // Calculate pure background and put weight on peak if using Rwp
423 purePeakWS->mutableE(0).assign(purePeakWS->y(0).size(), 1.0);
424 size_t i = 0;
425 std::transform(purePeakWS->y(0).cbegin(), purePeakWS->y(0).cend(), purePeakWS->mutableY(0).begin(),
426 [=](const double &y) mutable {
427 double newY = y - bkgdvalues[i++];
428 return std::max(0.0, newY);
429 });
430}
431
432//----------------------------------------------------------------------------------------------
440 size_t wsindex, double startx, double endx) {
441 // Check validity and debug output
442 if (!peakfunc)
443 throw std::runtime_error("fitPeakFunction's input peakfunc has not been initialized.");
444
445 m_sstream << "Function (to fit): " << peakfunc->asString() << " From " << startx << " to " << endx << ".\n";
446
447 double goodness = fitFunctionSD(peakfunc, dataws, wsindex, startx, endx);
448
449 return goodness;
450}
451
452//-----------------------------------------------------------------------
453//----------------------
460 m_numFitCalls = 0;
461
462 // Check and initialization
463 string errmsg;
464 if (!hasSetupToFitPeak(errmsg)) {
465 g_log.error(errmsg);
466 throw runtime_error("Object has not been set up completely to fit peak.");
467 } else {
468 m_sstream << "F1158: Well-setup and good to go!\n";
469 }
470
471 m_bestRwp = DBL_MAX;
472
473 // Fit background
475 // User's input peak range cannot be trusted. Data might be noisy
476 stringstream outss;
477 outss << "User specified peak range cannot be trusted! Because peak range "
478 "overlap fit window. "
479 << "Number of data points in fitting window = " << i_maxFitX - i_minFitX
480 << ". A UNRELIABLE algorithm is used to guess peak range. ";
481 g_log.warning(outss.str());
482
483 size_t numpts = i_maxFitX - i_minFitX;
484 auto shift = static_cast<size_t>(static_cast<double>(numpts) / 6.);
485 i_minPeakX += shift;
486 const auto &Xdata = m_dataWS->x(m_wsIndex);
487 if (i_minPeakX >= Xdata.size())
488 i_minPeakX = Xdata.size() - 1;
489 m_minPeakX = Xdata[i_minPeakX];
490
491 if (i_maxPeakX < shift) {
492 i_maxPeakX = 0;
493 } else {
494 i_maxPeakX -= shift;
495 }
496 m_maxPeakX = Xdata[i_maxPeakX];
497 }
498
500
501 // Generate partial workspace within given fit window
503
504 // Remove background to make a pure peak
505 removeBackground(purePeakWS);
506
507 // Estimate the peak height
508 double est_peakheight = estimatePeakHeight(m_peakFunc, purePeakWS, 0, 0, purePeakWS->x(0).size() - 1);
509 m_peakFunc->setHeight(est_peakheight);
510
511 // Store starting setup
513
514 Progress progress(this, 0.0, 1.0, m_vecFWHM.size());
515
516 // Fit with different starting values of peak width
517 for (size_t i = 0; i < m_vecFWHM.size(); ++i) {
518 // Restore
519 if (i > 0)
521
522 // Set FWHM
523 m_peakFunc->setFwhm(m_vecFWHM[i]);
524 m_sstream << "Round " << i << " of " << m_vecFWHM.size() << ". Using proposed FWHM = " << m_vecFWHM[i] << "\n";
525
526 // Fit
527 double rwp = fitPeakFunction(m_peakFunc, purePeakWS, 0, m_minFitX, m_maxFitX);
528
529 m_sstream << "Fit peak function cost = " << rwp << "\n";
530
531 // Store result
532 processNStoreFitResult(rwp, false);
533
534 progress.report();
535 }
536
537 // Get best fitting peak function and Make a combo fit
539
540 // Fit the composite function as final
542 m_bestRwp = compcost;
543
544 m_sstream << "MultStep-Fit: Best Fitted Peak: " << m_peakFunc->asString() << ". Final " << m_costFunction << " = "
545 << compcost << "\n"
546 << "Number of calls on Fit = " << m_numFitCalls << "\n";
547}
548
549//----------------------------------------------------------------------------------------------
554std::map<std::string, double> FitOneSinglePeak::backup(const IFunction_const_sptr &func) {
555 std::map<std::string, double> funcparammap;
556
557 // Set up
558 vector<string> funcparnames = func->getParameterNames();
559 size_t nParam = funcparnames.size();
560 for (size_t i = 0; i < nParam; ++i) {
561 double parvalue = func->getParameter(i);
562 funcparammap.emplace(funcparnames[i], parvalue);
563 }
564
565 return funcparammap;
566}
567
568//----------------------------------------------------------------------------------------------
574std::map<std::string, double> FitOneSinglePeak::storeFunctionError(const IFunction_const_sptr &func) {
575 // output map
576 std::map<std::string, double> paramerrormap;
577
578 // Get function error and store in output map
579 vector<string> funcparnames = func->getParameterNames();
580 size_t nParam = funcparnames.size();
581 for (size_t i = 0; i < nParam; ++i) {
582 double parerror = func->getError(i);
583 paramerrormap.emplace(funcparnames[i], parerror);
584 }
585
586 return paramerrormap;
587}
588
589//----------------------------------------------------------------------------------------------
592void FitOneSinglePeak::pop(const std::map<std::string, double> &funcparammap, const API::IFunction_sptr &func) {
593 std::map<std::string, double>::const_iterator miter;
594 for (miter = funcparammap.begin(); miter != funcparammap.end(); ++miter) {
595 string parname = miter->first;
596 double parvalue = miter->second;
597 func->setParameter(parname, parvalue);
598 }
599}
600
601//----------------------------------------------------------------------------------------------
612 size_t wsindex, double xmin, double xmax) {
613 // Set up sub algorithm fit
614 IAlgorithm_sptr fit;
615 try {
616 fit = createChildAlgorithm("CalculateChiSquared", -1, -1, false);
617 } catch (Exception::NotFoundError &) {
618 std::stringstream errss;
619 errss << "The FitPeak algorithm requires the CurveFitting library";
620 g_log.error(errss.str());
621 throw std::runtime_error(errss.str());
622 }
623
624 // Set the properties
625 fit->setProperty("Function", fitfunc);
626 fit->setProperty("InputWorkspace", dataws);
627 fit->setProperty("WorkspaceIndex", static_cast<int>(wsindex));
628 fit->setProperty("StartX", xmin);
629 fit->setProperty("EndX", xmax);
630
631 fit->executeAsChildAlg();
632 if (!fit->isExecuted()) {
633 g_log.error("Fit for background is not executed. ");
634 throw std::runtime_error("Fit for background is not executed. ");
635 }
636
637 // Retrieve result
638 const double chi2 = fit->getProperty("ChiSquaredWeightedDividedByNData");
639
640 return chi2;
641}
642
643//----------------------------------------------------------------------------------------------
650double FitOneSinglePeak::fitFunctionSD(IFunction_sptr fitfunc, const MatrixWorkspace_sptr &dataws, size_t wsindex,
651 double xmin, double xmax) {
652 // Set up sub algorithm fit
653 IAlgorithm_sptr fit;
654 try {
655 fit = createChildAlgorithm("Fit", -1, -1, false);
656 } catch (Exception::NotFoundError &) {
657 std::stringstream errss;
658 errss << "The FitPeak algorithm requires the CurveFitting library";
659 g_log.error(errss.str());
660 throw std::runtime_error(errss.str());
661 }
662
663 // Set the properties
664 fit->setProperty("Function", fitfunc);
665 fit->setProperty("InputWorkspace", dataws);
666 fit->setProperty("WorkspaceIndex", static_cast<int>(wsindex));
667 fit->setProperty("MaxIterations", 50); // magic number
668 fit->setProperty("StartX", xmin);
669 fit->setProperty("EndX", xmax);
670 fit->setProperty("Minimizer", m_minimizer);
671 fit->setProperty("CostFunction", m_costFunction);
672 fit->setProperty("CalcErrors", true);
673
674 // Execute fit and get result of fitting background
675 m_sstream << "FitSingleDomain: " << fit->asString() << ".\n";
676
677 fit->executeAsChildAlg();
678 if (!fit->isExecuted()) {
679 g_log.error("Fit for background is not executed. ");
680 throw std::runtime_error("Fit for background is not executed. ");
681 }
683
684 // Retrieve result
685 std::string fitStatus = fit->getProperty("OutputStatus");
686 double chi2 = EMPTY_DBL();
687 if (fitStatus == "success") {
688 chi2 = fit->getProperty("OutputChi2overDoF");
689 fitfunc = fit->getProperty("Function");
690 }
691
692 // Debug information
693 m_sstream << "[F1201] FitSingleDomain Fitted-Function " << fitfunc->asString() << ": Fit-status = " << fitStatus
694 << ", chi^2 = " << chi2 << ".\n";
695
696 return chi2;
697}
698
699//----------------------------------------------------------------------------------------------
708 size_t wsindex, vector<double> vec_xmin, vector<double> vec_xmax) {
709 // Validate
710 if (vec_xmin.size() != vec_xmax.size())
711 throw runtime_error("Sizes of xmin and xmax (vectors) are not equal. ");
712
713 // Set up sub algorithm fit
714 IAlgorithm_sptr fit;
715 try {
716 fit = createChildAlgorithm("Fit", -1, -1, true);
717 } catch (Exception::NotFoundError &) {
718 std::stringstream errss;
719 errss << "The FitPeak algorithm requires the CurveFitting library";
720 g_log.error(errss.str());
721 throw std::runtime_error(errss.str());
722 }
723
724 // This use multi-domain; but does not know how to set up
725 std::shared_ptr<MultiDomainFunction> funcmd = std::make_shared<MultiDomainFunction>();
726
727 // After a change of the default value of NumDeriv in MultiDomainFunction
728 // this needs to be set to false to preserve the original behaviour.
729 // Results of this algorithm as well as algorithms that use it
730 // seem to be very sensitive to the accuracy of the derivatives.
731 funcmd->setAttributeValue("NumDeriv", false);
732
733 // Set function first
734 funcmd->addFunction(fitfunc);
735
736 // set domain for function with index 0 covering both sides
737 funcmd->clearDomainIndices();
738 std::vector<size_t> ii(2);
739 ii[0] = 0;
740 ii[1] = 1;
741 funcmd->setDomainIndices(0, ii);
742
743 // Set the properties
744 fit->setProperty("Function", std::dynamic_pointer_cast<IFunction>(funcmd));
745 fit->setProperty("InputWorkspace", dataws);
746 fit->setProperty("WorkspaceIndex", static_cast<int>(wsindex));
747 fit->setProperty("StartX", vec_xmin[0]);
748 fit->setProperty("EndX", vec_xmax[0]);
749 fit->setProperty("InputWorkspace_1", dataws);
750 fit->setProperty("WorkspaceIndex_1", static_cast<int>(wsindex));
751 fit->setProperty("StartX_1", vec_xmin[1]);
752 fit->setProperty("EndX_1", vec_xmax[1]);
753 fit->setProperty("MaxIterations", 50);
754 fit->setProperty("Minimizer", m_minimizer);
755 fit->setProperty("CostFunction", "Least squares");
756
757 m_sstream << "FitMultiDomain: Funcion " << funcmd->name() << ": "
758 << "Range: (" << vec_xmin[0] << ", " << vec_xmax[0] << ") and (" << vec_xmin[1] << ", " << vec_xmax[1]
759 << "); " << funcmd->asString() << "\n";
760
761 // Execute
762 fit->execute();
763 if (!fit->isExecuted()) {
764 throw runtime_error("Fit is not executed on multi-domain function/data. ");
765 }
767
768 // Retrieve result
769 std::string fitStatus = fit->getProperty("OutputStatus");
770 m_sstream << "[DB] Multi-domain fit status: " << fitStatus << ".\n";
771
772 double chi2 = EMPTY_DBL();
773 if (fitStatus == "success") {
774 chi2 = fit->getProperty("OutputChi2overDoF");
775 m_sstream << "FitMultidomain: Successfully-Fitted Function " << fitfunc->asString() << ", Chi^2 = " << chi2 << "\n";
776 }
777
778 return chi2;
779}
780
781//----------------------------------------------------------------------------------------------
792 const API::IBackgroundFunction_sptr &bkgdfunc,
793 const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double startx,
794 double endx) {
795 // Construct composit function
796 std::shared_ptr<CompositeFunction> compfunc = std::make_shared<CompositeFunction>();
797 compfunc->addFunction(peakfunc);
798 compfunc->addFunction(bkgdfunc);
799
800 // Do calculation for starting chi^2/Rwp: as the assumption that the input the
801 // so far the best Rwp
802 // FIXME - This is not a good practise...
803 double backRwp = calChiSquareSD(bkgdfunc, dataws, wsindex, startx, endx);
804 m_sstream << "Background: Pre-fit Goodness = " << backRwp << "\n";
805 m_bestRwp = calChiSquareSD(compfunc, dataws, wsindex, startx, endx);
806 m_sstream << "Peak+Background: Pre-fit Goodness = " << m_bestRwp << "\n";
807
808 auto bkuppeakmap = backup(peakfunc);
809 auto bkupbkgdmap = backup(bkgdfunc);
812
813 // Fit
814 double goodness = fitFunctionSD(compfunc, dataws, wsindex, startx, endx);
815 string errorreason;
816
817 // Check fit result
818 goodness = checkFittedPeak(peakfunc, goodness, errorreason);
819
820 if (!errorreason.empty())
821 m_sstream << "Error reason of fit peak+background composite: " << errorreason << "\n";
822
823 double goodness_final = DBL_MAX;
824 if (goodness <= m_bestRwp && goodness <= backRwp) {
825 // Fit for composite function renders a better result
826 goodness_final = goodness;
827 processNStoreFitResult(goodness_final, true);
828 } else if (goodness > m_bestRwp && m_bestRwp < DBL_MAX && m_bestRwp <= backRwp) {
829 // A worse result is got. Revert to original function parameters
830 m_sstream << "Fit peak/background composite function FAILS to render a "
831 "better solution. "
832 << "Input cost function value = " << m_bestRwp << ", output cost function value = " << goodness << "\n";
833
834 pop(bkuppeakmap, peakfunc);
835 pop(bkupbkgdmap, bkgdfunc);
836 goodness_final = m_bestRwp;
837 } else {
838 m_sstream << "Fit peak-background function fails in all approaches! \n";
839 }
840
841 return goodness_final;
842}
843
844//----------------------------------------------------------------------------------------------
848double FitOneSinglePeak::checkFittedPeak(const IPeakFunction_sptr &peakfunc, double costfuncvalue,
849 std::string &errorreason) {
850 if (costfuncvalue < DBL_MAX) {
851 // Fit is successful. Check whether the fit result is physical
852 stringstream errorss;
853 double peakcentre = peakfunc->centre();
854 if (peakcentre < m_minPeakX || peakcentre > m_maxPeakX) {
855 errorss << "Peak centre (at " << peakcentre << " ) is out of specified range )" << m_minPeakX << ", "
856 << m_maxPeakX << "). ";
857 costfuncvalue = DBL_MAX;
858 }
859
860 double peakheight = peakfunc->height();
861 if (peakheight < 0) {
862 errorss << "Peak height (" << peakheight << ") is negative. ";
863 costfuncvalue = DBL_MAX;
864 }
865 double peakfwhm = peakfunc->fwhm();
866 if (peakfwhm > (m_maxFitX - m_minFitX) * MAGICNUMBER) {
867 errorss << "Peak width is unreasonably wide. ";
868 costfuncvalue = DBL_MAX;
869 }
870 errorreason = errorss.str();
871 } else {
872 // Fit is not successful
873 errorreason = "Fit() on peak function is NOT successful.";
874 }
875
876 return costfuncvalue;
877}
878
879//----------------------------------------------------------------------------------------------
885 // Back up background function
886 m_bkupBkgdFunc = backup(bkgdfunc);
887
888 // Fit in multiple domain
889 vector<double> vec_xmin(2);
890 vector<double> vec_xmax(2);
891 vec_xmin[0] = m_minFitX;
892 vec_xmin[1] = m_maxPeakX;
893 vec_xmax[0] = m_minPeakX;
894 vec_xmax[1] = m_maxFitX;
895 double chi2 = fitFunctionMD(std::dynamic_pointer_cast<IFunction>(bkgdfunc), m_dataWS, m_wsIndex, vec_xmin, vec_xmax);
896
897 // Process fit result
898 if (chi2 < DBL_MAX - 1) {
899 // Store fitting result
900 m_bestBkgdFunc = backup(bkgdfunc);
902 } else {
903 // Restore background function
904 pop(m_bkupBkgdFunc, bkgdfunc);
905 }
906
907 return bkgdfunc;
908}
909
910//----------------------------------------------------------------------------------------------
915void FitOneSinglePeak::processNStoreFitResult(double rwp, bool storebkgd) {
916 bool fitsuccess = true;
917 string failreason;
918
919 if (rwp < DBL_MAX) {
920 // A valid Rwp returned from Fit
921
922 // Check non-negative height
923 double f_height = m_peakFunc->height();
924 if (f_height <= 0.) {
925 rwp = DBL_MAX;
926 failreason += "Negative peak height. ";
927 fitsuccess = false;
928 }
929
930 // Check peak position
931 double f_centre = m_peakFunc->centre();
933 // Peak position criteria is on position tolerance
935 rwp = DBL_MAX;
936 failreason = "Peak centre out of tolerance. ";
937 fitsuccess = false;
938 }
939 } else if (f_centre < m_minPeakX || f_centre > m_maxPeakX) {
940 rwp = DBL_MAX;
941 failreason += "Peak centre out of input peak range ";
942 m_sstream << "Peak centre " << f_centre << " is out of peak range: " << m_minPeakX << ", " << m_maxPeakX << "\n";
943 fitsuccess = false;
944 }
945
946 } // RWP fine
947 else {
948 failreason = "(Single-step) Fit returns a DBL_MAX.";
949 fitsuccess = false;
950 }
951
952 m_sstream << "Process fit result: "
953 << "Rwp = " << rwp << ", best Rwp = " << m_bestRwp << ", Fit success = " << fitsuccess << ". ";
954
955 // Store result if
956 if (rwp < m_bestRwp && fitsuccess) {
959 if (storebkgd) {
962 }
963 m_bestRwp = rwp;
964
965 m_sstream << "Store result and new Best RWP = " << m_bestRwp << ".\n";
966 } else if (!fitsuccess) {
967 m_sstream << "Reason of fit's failure: " << failreason << "\n";
968 }
969}
970
971//----------------------------------------------------------------------------------------------
975
976//----------------------------------------------------------------------------------------------
980std::map<std::string, double> FitOneSinglePeak::getPeakError() { return m_fitErrorPeakFunc; }
981
985std::map<std::string, double> FitOneSinglePeak::getBackgroundError() { return m_fitErrorBkgdFunc; }
986
987//----------------------------------------------------------------------------------------------
988void FitOneSinglePeak::exec() { throw runtime_error("Not used."); }
989
990//----------------------------------------------------------------------------------------------
991void FitOneSinglePeak::init() { throw runtime_error("Not used."); }
992
993//----------------------------------------------------------------------------------------------
994// FitPeak Methods
995//----------------------------------------------------------------------------------------------
996
998
999//----------------------------------------------------------------------------------------------
1003 : API::Algorithm(), m_dataWS(), m_wsIndex(0), m_peakFunc(), m_bkgdFunc(), m_minFitX(0.), m_maxFitX(0.),
1004 m_minPeakX(0.), m_maxPeakX(0.), i_minFitX(0), i_maxFitX(0), i_minPeakX(0), i_maxPeakX(0), m_fitBkgdFirst(false),
1005 m_outputRawParams(false), m_userGuessedFWHM(0.), m_userPeakCentre(0.), m_minGuessedPeakWidth(0),
1006 m_maxGuessedPeakWidth(0), m_fwhmFitStep(0), m_fitWithStepPeakWidth(false), m_usePeakPositionTolerance(false),
1007 m_peakPositionTolerance(0.), m_peakParameterTableWS(), m_bkgdParameterTableWS(), m_peakParameterNames(),
1008 m_bkgdParameterNames(), m_minimizer("Levenberg-MarquardtMD"), m_bkupBkgdFunc(), m_bkupPeakFunc(),
1009 m_bestPeakFunc(), m_bestBkgdFunc(), m_bestRwp(DBL_MAX), m_finalGoodnessValue(0.), m_vecybkup(), m_vecebkup(),
1010 m_costFunction(), m_lightWeightOutput(false) {}
1011
1012//----------------------------------------------------------------------------------------------
1016 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
1017 "Name of the input workspace for peak fitting.");
1018
1019 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
1020 "Name of the output workspace containing fitted peak.");
1021
1022 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("ParameterTableWorkspace", "", Direction::Output),
1023 "Name of the table workspace containing the fitted parameters. ");
1024
1025 std::shared_ptr<BoundedValidator<int>> mustBeNonNegative = std::make_shared<BoundedValidator<int>>();
1026 mustBeNonNegative->setLower(0);
1027 declareProperty("WorkspaceIndex", 0, mustBeNonNegative, "Workspace index ");
1028
1029 std::vector<std::string> peakNames = FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
1030 vector<string> peakFullNames = addFunctionParameterNames(peakNames);
1031 declareProperty("PeakFunctionType", "", std::make_shared<StringListValidator>(peakFullNames), "Peak function type. ");
1032
1033 declareProperty(std::make_unique<ArrayProperty<string>>("PeakParameterNames"), "List of peak parameter names. ");
1034
1035 declareProperty(std::make_unique<ArrayProperty<double>>("PeakParameterValues"),
1036 "List of peak parameter values. They must have a 1-to-1 "
1037 "mapping to PeakParameterNames list. ");
1038
1039 declareProperty(std::make_unique<ArrayProperty<double>>("FittedPeakParameterValues", Direction::Output),
1040 "Fitted peak parameter values. ");
1041
1042 vector<string> bkgdtypes{"Flat", "Flat (A0)", "Linear", "Linear (A0, A1)", "Quadratic", "Quadratic (A0, A1, A2)"};
1043 declareProperty("BackgroundType", "Linear", std::make_shared<StringListValidator>(bkgdtypes), "Type of Background.");
1044
1045 declareProperty(std::make_unique<ArrayProperty<string>>("BackgroundParameterNames"),
1046 "List of background parameter names. ");
1047
1048 declareProperty(std::make_unique<ArrayProperty<double>>("BackgroundParameterValues"),
1049 "List of background parameter values. "
1050 "They must have a 1-to-1 mapping to BackgroundParameterNames list. ");
1051
1052 declareProperty(std::make_unique<ArrayProperty<double>>("FittedBackgroundParameterValues", Direction::Output),
1053 "Fitted background parameter values. ");
1054
1055 declareProperty(std::make_unique<ArrayProperty<double>>("FitWindow"),
1056 "Enter a comma-separated list of the expected X-position of "
1057 "windows to fit. "
1058 "The number of values must be 2.");
1059
1060 declareProperty(std::make_unique<ArrayProperty<double>>("PeakRange"),
1061 "Enter a comma-separated list of expected x-position as peak range. "
1062 "The number of values must be 2.");
1063
1064 declareProperty("FitBackgroundFirst", true,
1065 "If true, then the algorithm will fit background first. "
1066 "And then the peak. ");
1067
1068 declareProperty("RawParams", true,
1069 "If true, then the output table workspace contains the raw "
1070 "profile parameter. "
1071 "Otherwise, the effective parameters will be written. ");
1072
1073 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
1074 mustBePositive->setLower(1);
1075 declareProperty("MinGuessedPeakWidth", 2, mustBePositive,
1076 "Minimum guessed peak width for fit. It is in unit of number of pixels.");
1077
1078 declareProperty("MaxGuessedPeakWidth", 10, mustBePositive,
1079 "Maximum guessed peak width for fit. It is in unit of number of pixels.");
1080
1081 declareProperty("GuessedPeakWidthStep", EMPTY_INT(), mustBePositive,
1082 "Step of guessed peak width. It is in unit of number of pixels.");
1083
1084 auto mustBePostiveDbl = std::make_shared<BoundedValidator<double>>();
1085 mustBePostiveDbl->setLower(DBL_MIN);
1086 declareProperty("PeakPositionTolerance", EMPTY_DBL(), mustBePostiveDbl,
1087 "Peak position tolerance. If fitted peak's position differs "
1088 "from proposed value more than "
1089 "the given value, fit is treated as failure. ");
1090
1091 std::array<string, 2> costFuncOptions = {{"Chi-Square", "Rwp"}};
1092 declareProperty("CostFunction", "Chi-Square",
1093 Kernel::IValidator_sptr(new Kernel::ListValidator<std::string>(costFuncOptions)), "Cost functions");
1094
1095 std::vector<std::string> minimizerOptions = API::FuncMinimizerFactory::Instance().getKeys();
1096
1097 declareProperty("Minimizer", "Levenberg-Marquardt",
1099 "Minimizer to use for fitting. Minimizers available are "
1100 "\"Levenberg-Marquardt\", \"Simplex\","
1101 "\"Conjugate gradient (Fletcher-Reeves imp.)\", \"Conjugate "
1102 "gradient (Polak-Ribiere imp.)\", \"BFGS\", and "
1103 "\"Levenberg-MarquardtMD\"");
1104
1105 declareProperty("CostFunctionValue", DBL_MAX, "Value of cost function of the fitted peak. ",
1107}
1108
1109//----------------------------------------------------------------------------------------------
1113 // Get input properties
1115
1116 // Create functions
1118
1119 // Check input function, guessed value, and etc.
1121
1122 // Set parameters to fit
1123 FitOneSinglePeak fit1peakalg;
1124
1125 fit1peakalg.setFunctions(m_peakFunc, m_bkgdFunc);
1126 fit1peakalg.setWorskpace(m_dataWS, m_wsIndex);
1127
1129 fit1peakalg.setFitWindow(m_minFitX, m_maxFitX);
1130 fit1peakalg.setPeakRange(m_minPeakX, m_maxPeakX);
1133
1135
1136 if (m_fitBkgdFirst) {
1137 fit1peakalg.highBkgdFit();
1138 } else {
1139 fit1peakalg.simpleFit();
1140 }
1141 string dbmsg = fit1peakalg.getDebugMessage();
1142 g_log.information(dbmsg);
1143
1145
1146 // Output
1147 setupOutput(fit1peakalg.getPeakError(), fit1peakalg.getBackgroundError());
1148}
1149
1150//----------------------------------------------------------------------------------------------
1153std::vector<std::string> FitPeak::addFunctionParameterNames(const std::vector<std::string> &funcnames) {
1154 vector<string> vec_funcparnames;
1155
1156 for (auto &funcname : funcnames) {
1157 // Add original name in
1158 vec_funcparnames.emplace_back(funcname);
1159
1160 // Add a full function name and parameter names in
1161 IFunction_sptr tempfunc = FunctionFactory::Instance().createFunction(funcname);
1162
1163 stringstream parnamess;
1164 parnamess << funcname << " (";
1165 vector<string> funcpars = tempfunc->getParameterNames();
1166 for (size_t j = 0; j < funcpars.size(); ++j) {
1167 parnamess << funcpars[j];
1168 if (j != funcpars.size() - 1)
1169 parnamess << ", ";
1170 }
1171 parnamess << ")";
1172
1173 vec_funcparnames.emplace_back(parnamess.str());
1174 }
1175
1176 return vec_funcparnames;
1177}
1178
1179//----------------------------------------------------------------------------------------------
1184 // Data workspace (input)
1185 m_dataWS = getProperty("InputWorkspace");
1186 int tempint = getProperty("WorkspaceIndex");
1187 m_wsIndex = static_cast<size_t>(tempint);
1188
1189 // Fit window
1190 auto &vecX = m_dataWS->x(m_wsIndex);
1191
1192 vector<double> fitwindow = getProperty("FitWindow");
1193 if (fitwindow.size() != 2) {
1194 throw runtime_error("Must enter 2 and only 2 items in fit window. ");
1195 }
1196 m_minFitX = fitwindow[0];
1197 if (m_minFitX < vecX.front())
1198 m_minFitX = vecX.front();
1199 m_maxFitX = fitwindow[1];
1200 if (m_maxFitX > vecX.back())
1201 m_maxFitX = vecX.back();
1202
1203 if (m_maxFitX <= m_minFitX) {
1204 stringstream errss;
1205 errss << "Minimum X (" << m_minFitX << ") is larger and equal to maximum X (" << m_maxFitX
1206 << ") to fit. It is not allowed. ";
1207 g_log.error(errss.str());
1208 throw runtime_error(errss.str());
1209 }
1210
1211 // Peak range
1212 vector<double> peakrange = getProperty("PeakRange");
1213 if (peakrange.size() != 2) {
1214 throw runtime_error("Must enter 2 and only 2 items for PeakRange in fit window. ");
1215 }
1216 m_minPeakX = peakrange[0];
1217 m_maxPeakX = peakrange[1];
1218 if (m_maxFitX <= m_minFitX) {
1219 stringstream errss;
1220 errss << "Minimum peak range (" << m_minPeakX << ") is larger and equal to maximum X (" << m_maxPeakX
1221 << ") of the range of peak. It is not allowed. ";
1222 g_log.error(errss.str());
1223 throw runtime_error(errss.str());
1224 }
1225
1226 if (m_minPeakX < m_minFitX) {
1228 g_log.warning() << "Minimum peak range is out side of the lower boundary "
1229 "of fit window. ";
1230 }
1231 if (m_maxPeakX > m_maxFitX) {
1233 g_log.warning() << "Maximum peak range is out side of the upper boundary "
1234 "of fit window. ";
1235 }
1236
1237 // Fit strategy
1238 m_fitBkgdFirst = getProperty("FitBackgroundFirst");
1239
1240 // Trying FWHM in a certain range
1241 m_minGuessedPeakWidth = getProperty("MinGuessedPeakWidth");
1242 m_maxGuessedPeakWidth = getProperty("MaxGuessedPeakWidth");
1243 m_fwhmFitStep = getProperty("GuessedPeakWidthStep");
1245 m_fitWithStepPeakWidth = false;
1246 else {
1249 std::stringstream errss;
1250 errss << "User specified wrong guessed peak width parameters (must be "
1251 "postive and make sense). "
1252 << "User inputs are min = " << m_minGuessedPeakWidth << ", max = " << m_maxGuessedPeakWidth
1253 << ", step = " << m_fwhmFitStep;
1254 g_log.error(errss.str());
1255 throw std::runtime_error(errss.str());
1256 }
1257 }
1258
1259 // Tolerance
1260 m_peakPositionTolerance = getProperty("PeakPositionTolerance");
1263 else
1265
1266 // Cost function
1267 string costfunname = getProperty("CostFunction");
1268 if (costfunname == "Chi-Square") {
1269 m_costFunction = "Least squares";
1270 } else if (costfunname == "Rwp") {
1271 m_costFunction = "Rwp";
1272 } else {
1273 g_log.error() << "Cost function " << costfunname << " is not supported. "
1274 << "\n";
1275 throw runtime_error("Cost function is not supported. ");
1276 }
1277
1278 // Minimizer
1279 m_minimizer = getPropertyValue("Minimizer");
1280
1281 // Output option
1282 m_outputRawParams = getProperty("RawParams");
1283}
1284
1285//----------------------------------------------------------------------------------------------
1289 //=========================================================================
1290 // Generate background function
1291 //=========================================================================
1292 // Get background type
1293 string bkgdtyperaw = getPropertyValue("BackgroundType");
1294 bool usedefaultbkgdparorder = false;
1295 string bkgdtype = parseFunctionTypeFull(bkgdtyperaw, usedefaultbkgdparorder);
1296
1297 // FIXME - Fix the inconsistency in nameing the background
1298 if (bkgdtype == "Flat" || bkgdtype == "Linear")
1299 bkgdtype += "Background";
1300
1301 // Generate background function
1302 m_bkgdFunc = std::dynamic_pointer_cast<IBackgroundFunction>(FunctionFactory::Instance().createFunction(bkgdtype));
1303
1304 // Set background function parameter values
1305 m_bkgdParameterNames = getProperty("BackgroundParameterNames");
1306 if (usedefaultbkgdparorder && m_bkgdParameterNames.empty()) {
1307 m_bkgdParameterNames = m_bkgdFunc->getParameterNames();
1308 } else if (m_bkgdParameterNames.empty()) {
1309 throw runtime_error("In the non-default background parameter name mode, "
1310 "user must give out parameter names. ");
1311 }
1312
1313 vector<double> vec_bkgdparvalues = getProperty("BackgroundParameterValues");
1314 if (m_bkgdParameterNames.size() != vec_bkgdparvalues.size()) {
1315 stringstream errss;
1316 errss << "Input background properties' arrays are incorrect: # of "
1317 "parameter names = "
1318 << m_bkgdParameterNames.size() << ", # of parameter values = " << vec_bkgdparvalues.size() << "\n";
1319 g_log.error(errss.str());
1320 throw runtime_error(errss.str());
1321 }
1322
1323 // Set parameter values
1324 for (size_t i = 0; i < m_bkgdParameterNames.size(); ++i) {
1325 m_bkgdFunc->setParameter(m_bkgdParameterNames[i], vec_bkgdparvalues[i]);
1326 }
1327
1328 //=========================================================================
1329 // Generate peak function
1330 //=========================================================================
1331 string peaktypeprev = getPropertyValue("PeakFunctionType");
1332 bool defaultparorder = true;
1333 string peaktype = parseFunctionTypeFull(peaktypeprev, defaultparorder);
1334 m_peakFunc = std::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction(peaktype));
1335
1336 // Peak parameters' names
1337 m_peakParameterNames = getProperty("PeakParameterNames");
1338 if (m_peakParameterNames.empty()) {
1339 if (defaultparorder) {
1340 // Use default peak parameter names' order
1341 m_peakParameterNames = m_peakFunc->getParameterNames();
1342 } else {
1343 throw runtime_error("Peak parameter names' input is not in default mode. "
1344 "It cannot be left empty. ");
1345 }
1346 }
1347
1348 // Peak parameters' value
1349 vector<double> vec_peakparvalues = getProperty("PeakParameterValues");
1350 if (m_peakParameterNames.size() != vec_peakparvalues.size()) {
1351 stringstream errss;
1352 errss << "Input peak properties' arrays are incorrect: # of parameter "
1353 "names = "
1354 << m_peakParameterNames.size() << ", # of parameter values = " << vec_peakparvalues.size() << "\n";
1355 throw runtime_error(errss.str());
1356 }
1357
1358 // Set peak parameter values
1359 for (size_t i = 0; i < m_peakParameterNames.size(); ++i) {
1360 m_peakFunc->setParameter(m_peakParameterNames[i], vec_peakparvalues[i]);
1361 }
1362}
1363
1364//----------------------------------------------------------------------------------------------
1367std::string FitPeak::parseFunctionTypeFull(const std::string &fullstring, bool &defaultparorder) {
1368 string peaktype;
1369
1370 size_t n = std::count(fullstring.begin(), fullstring.end(), '(');
1371 if (n > 0) {
1372 peaktype = fullstring.substr(0, fullstring.find('('));
1373 boost::algorithm::trim(peaktype);
1374 defaultparorder = true;
1375 } else {
1376 peaktype = fullstring;
1377 defaultparorder = false;
1378 }
1379
1380 return peaktype;
1381}
1382
1383//----------------------------------------------------------------------------------------------
1387 // Check functions
1388 if (!m_peakFunc || !m_bkgdFunc)
1389 throw runtime_error("Either peak function or background function has not been set up.");
1390
1391 // Check validity on peak centre
1392 double centre_guess = m_peakFunc->centre();
1393 if (m_minFitX >= centre_guess || m_maxFitX <= centre_guess)
1394 throw runtime_error("Peak centre is out side of fit window. ");
1395
1396 // Peak width and centre: from user input
1397 m_userGuessedFWHM = m_peakFunc->fwhm();
1398 m_userPeakCentre = m_peakFunc->centre();
1399}
1400
1401//----------------------------------------------------------------------------------------------
1405void FitPeak::setupOutput(const std::map<std::string, double> &m_fitErrorPeakFunc,
1406 const std::map<std::string, double> &m_fitErrorBkgdFunc) {
1407 // TODO - Need to retrieve useful information from FitOneSinglePeak object
1408 // (think of how)
1409 const auto &vecX = m_dataWS->x(m_wsIndex);
1410 const size_t i_minFitX = getIndex(vecX, m_minFitX);
1411 const size_t i_maxFitX = getIndex(vecX, m_maxFitX);
1412
1413 // Data workspace
1414 const size_t nspec = 3;
1415 // Get a vector for fit window
1416
1417 vector<double> vecoutx(i_maxFitX - i_minFitX + 1);
1418 for (size_t i = i_minFitX; i <= i_maxFitX; ++i)
1419 vecoutx[i - i_minFitX] = vecX[i];
1420
1421 // Create workspace
1422 const size_t sizex = vecoutx.size();
1423 const auto sizey = sizex;
1424 HistogramBuilder builder;
1425 builder.setX(sizex);
1426 builder.setY(sizey);
1427 MatrixWorkspace_sptr outws = create<Workspace2D>(nspec, builder.build());
1428 // Calculate again
1429 FunctionDomain1DVector domain(vecoutx);
1430 FunctionValues values(domain);
1431
1432 CompositeFunction_sptr compfunc = std::make_shared<CompositeFunction>();
1433 compfunc->addFunction(m_peakFunc);
1434 compfunc->addFunction(m_bkgdFunc);
1435 compfunc->function(domain, values);
1436
1437 const auto domainVec = domain.toVector();
1438 outws->mutableX(0).assign(domainVec.cbegin(), domainVec.cend());
1439 outws->setSharedX(1, outws->sharedX(0));
1440 outws->setSharedX(2, outws->sharedX(0));
1441
1442 auto &vecY = m_dataWS->y(m_wsIndex);
1443 const auto valvec = values.toVector();
1444 outws->mutableY(0).assign(vecY.cbegin() + i_minFitX, vecY.cbegin() + i_minFitX + sizey);
1445 outws->mutableY(1).assign(valvec.cbegin(), valvec.cbegin() + sizey);
1446 std::transform(outws->y(0).cbegin(), outws->y(0).cbegin() + sizey, outws->y(1).cbegin(), outws->mutableY(2).begin(),
1447 std::minus<double>());
1448 // Set property
1449 setProperty("OutputWorkspace", outws);
1450
1451 // Function parameter table workspaces
1452 TableWorkspace_sptr peaktablews = genOutputTableWS(m_peakFunc, m_fitErrorPeakFunc, m_bkgdFunc, m_fitErrorBkgdFunc);
1453 setProperty("ParameterTableWorkspace", peaktablews);
1454
1455 // Parameter vector
1456 vector<double> vec_fitpeak;
1457 vec_fitpeak.reserve(m_peakParameterNames.size());
1458 std::transform(m_peakParameterNames.cbegin(), m_peakParameterNames.cend(), std::back_inserter(vec_fitpeak),
1459 [this](const auto &peakParameterName) { return m_peakFunc->getParameter(peakParameterName); });
1460
1461 setProperty("FittedPeakParameterValues", vec_fitpeak);
1462
1463 // Background
1464 vector<double> vec_fitbkgd;
1465 vec_fitbkgd.reserve(m_bkgdParameterNames.size());
1466 std::transform(m_bkgdParameterNames.cbegin(), m_bkgdParameterNames.cend(), std::back_inserter(vec_fitbkgd),
1467 [this](const auto &bkgdParameterName) { return m_bkgdFunc->getParameter(bkgdParameterName); });
1468
1469 setProperty("FittedBackgroundParameterValues", vec_fitbkgd);
1470
1471 // Output chi^2 or Rwp
1472 setProperty("CostFunctionValue", m_finalGoodnessValue);
1473}
1474
1475//----------------------------------------------------------------------------------------------
1479size_t getIndex(const HistogramX &vecx, double x) {
1480 size_t index;
1481 if (x <= vecx.front()) {
1482 index = 0;
1483 } else if (x >= vecx.back()) {
1484 index = vecx.size() - 1;
1485 } else {
1486 vector<double>::const_iterator fiter;
1487 fiter = lower_bound(vecx.begin(), vecx.end(), x);
1488 index = static_cast<size_t>(fiter - vecx.begin());
1489 if (index == 0)
1490 throw runtime_error("It seems impossible to have this value. ");
1491 if (x - vecx[index - 1] < vecx[index] - x)
1492 --index;
1493 }
1494
1495 return index;
1496}
1497
1498//----------------------------------------------------------------------------------------------
1501TableWorkspace_sptr FitPeak::genOutputTableWS(const IPeakFunction_sptr &peakfunc, map<string, double> peakerrormap,
1502 const IBackgroundFunction_sptr &bkgdfunc,
1503 map<string, double> bkgderrormap) {
1504 // Empty table
1505 TableWorkspace_sptr outtablews = std::make_shared<TableWorkspace>();
1506 outtablews->addColumn("str", "Name");
1507 outtablews->addColumn("double", "Value");
1508 outtablews->addColumn("double", "Error");
1509
1510 // Set chi^2
1511 TableRow newrow = outtablews->appendRow();
1512 newrow << "ChiSquare" << m_finalGoodnessValue;
1513
1514 // Set peak paraemters
1515 newrow = outtablews->appendRow();
1516 newrow << peakfunc->name();
1517
1518 if (m_outputRawParams) {
1519 vector<string> peakparnames = peakfunc->getParameterNames();
1520 for (auto &parname : peakparnames) {
1521 double parvalue = peakfunc->getParameter(parname);
1522 double error = peakerrormap[parname];
1523 newrow = outtablews->appendRow();
1524 newrow << parname << parvalue << error;
1525 }
1526 } else {
1527 newrow = outtablews->appendRow();
1528 newrow << "centre" << peakfunc->centre();
1529
1530 newrow = outtablews->appendRow();
1531 newrow << "width" << peakfunc->fwhm();
1532
1533 newrow = outtablews->appendRow();
1534 newrow << "height" << peakfunc->height();
1535 }
1536
1537 // Set background paraemters
1538 newrow = outtablews->appendRow();
1539 newrow << bkgdfunc->name();
1540
1541 if (m_outputRawParams) {
1542 vector<string> bkgdparnames = bkgdfunc->getParameterNames();
1543 for (auto &parname : bkgdparnames) {
1544 double parvalue = bkgdfunc->getParameter(parname);
1545 double error = bkgderrormap[parname];
1546 newrow = outtablews->appendRow();
1547 newrow << parname << parvalue << error;
1548 }
1549 } else {
1550 string bkgdtype = getProperty("BackgroundType");
1551
1552 newrow = outtablews->appendRow();
1553 newrow << "backgroundintercept" << bkgdfunc->getParameter("A0");
1554 if (bkgdtype != "Flat") {
1555 newrow = outtablews->appendRow();
1556 newrow << "backgroundintercept" << bkgdfunc->getParameter("A1");
1557 }
1558 if (bkgdtype == "Quadratic") {
1559 newrow = outtablews->appendRow();
1560 newrow << "A2" << bkgdfunc->getParameter("A2");
1561 }
1562 }
1563
1564 return outtablews;
1565}
1566
1567} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
CostFunctions::CostFuncFitting & m_costFunction
The cost function.
const double MAGICNUMBER
Definition: FitPeak.cpp:47
double error
Definition: IndexPeaks.cpp:133
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
Implements FunctionDomain1D with its own storage in form of a std::vector.
std::vector< double > toVector() const
Convert to a vector.
A class to store values calculated by a function.
std::vector< double > toVector() const
Return the calculated values as a vector.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
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.
FitOneSinglePeak: a class to perform peak fitting on a single peak.
Definition: FitPeak.h:30
std::map< std::string, double > m_bestBkgdFunc
Best background parameters.
Definition: FitPeak.h:178
void setPeakRange(double xpeakleft, double xpeakright)
Set peak range.
Definition: FitPeak.cpp:112
double m_finalGoodnessValue
Final goodness value (Rwp/Chi-square)
Definition: FitPeak.h:207
void setupGuessedFWHM(double usrwidth, int minfwhm, int maxfwhm, int stepsize, bool fitwithsteppedfwhm)
Set peak width to guess.
Definition: FitPeak.cpp:156
void processNStoreFitResult(double rwp, bool storebkgd)
Process and store fit result.
Definition: FitPeak.cpp:915
double fitPeakFunction(const API::IPeakFunction_sptr &peakfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double startx, double endx)
Fit peak function (flexible)
Definition: FitPeak.cpp:439
size_t m_wsIndex
Input worskpace index.
Definition: FitPeak.h:155
std::string m_minimizer
Minimzer.
Definition: FitPeak.h:191
void pop(const std::map< std::string, double > &funcparammap, const API::IFunction_sptr &func)
Restore the parameters value to a function from a string/double map.
Definition: FitPeak.cpp:592
std::string getDebugMessage()
Get debug message.
Definition: FitPeak.cpp:276
double estimatePeakHeight(const API::IPeakFunction_const_sptr &peakfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, size_t ixmin, size_t ixmax)
Estimate the peak height from a set of data containing pure peaks.
Definition: FitPeak.cpp:377
double checkFittedPeak(const API::IPeakFunction_sptr &peakfunc, double costfuncvalue, std::string &errorreason)
Check a peak function whether it is valid comparing to user specified criteria.
Definition: FitPeak.cpp:848
void setFunctions(const API::IPeakFunction_sptr &peakfunc, const API::IBackgroundFunction_sptr &bkgdfunc)
Set functions.
Definition: FitPeak.cpp:82
double m_maxPeakX
peak right boundary (client-defined)
Definition: FitPeak.h:169
API::IPeakFunction_sptr m_peakFunc
Peak function.
Definition: FitPeak.h:148
std::map< std::string, double > getPeakError()
Get fitting error for peak function.
Definition: FitPeak.cpp:980
API::IBackgroundFunction_sptr fitBackground(API::IBackgroundFunction_sptr bkgdfunc)
Fit background of a given peak in a given range.
Definition: FitPeak.cpp:884
bool m_peakRangeSet
Flag whether the peak range is set.
Definition: FitPeak.h:139
std::map< std::string, double > getBackgroundError()
Get fitting error for background function.
Definition: FitPeak.cpp:985
std::map< std::string, double > m_bkupBkgdFunc
Backed up background function parameters.
Definition: FitPeak.h:183
double fitFunctionSD(API::IFunction_sptr fitfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double xmin, double xmax)
Fit function in single domain.
Definition: FitPeak.cpp:650
double calChiSquareSD(const API::IFunction_sptr &fitfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double xmin, double xmax)
Calculate chi-square of a single domain function.
Definition: FitPeak.cpp:611
API::MatrixWorkspace_sptr m_dataWS
Input data workspace.
Definition: FitPeak.h:153
std::map< std::string, double > m_bkupPeakFunc
Backed up peak function parameters.
Definition: FitPeak.h:181
std::map< std::string, double > storeFunctionError(const API::IFunction_const_sptr &func)
Store function fitting error.
Definition: FitPeak.cpp:574
double m_peakPositionTolerance
Peak position tolerance.
Definition: FitPeak.h:198
std::map< std::string, double > m_fitErrorBkgdFunc
Fit error of background function.
Definition: FitPeak.h:188
bool hasSetupToFitPeak(std::string &errmsg)
Check whether it is ready to fit.
Definition: FitPeak.cpp:249
std::map< std::string, double > backup(const API::IFunction_const_sptr &func)
Back up fit result.
Definition: FitPeak.cpp:554
std::map< std::string, double > m_fitErrorPeakFunc
Fit error of peak function.
Definition: FitPeak.h:186
bool m_peakWindowSet
Peak widnow is set up.
Definition: FitPeak.h:143
API::IBackgroundFunction_sptr m_bkgdFunc
Background function.
Definition: FitPeak.h:150
size_t i_maxFitX
index of m_maxFitX
Definition: FitPeak.h:164
std::map< std::string, double > m_bestPeakFunc
Best peak parameters.
Definition: FitPeak.h:176
std::stringstream m_sstream
String stream.
Definition: FitPeak.h:213
double m_userPeakCentre
Peak centre provided by user.
Definition: FitPeak.h:201
void highBkgdFit()
Fit peak first considering high background.
Definition: FitPeak.cpp:459
bool simpleFit()
Fit peak and background together.
Definition: FitPeak.cpp:281
bool m_usePeakPositionTolerance
Flag to apply peak position tolerance.
Definition: FitPeak.h:145
void setFittingMethod(std::string minimizer, const std::string &costfunction)
Set fitting method.
Definition: FitPeak.cpp:129
void setFitPeakCriteria(bool usepeakpostol, double peakpostol)
Set fitted peak parameters' criterial including (a) peak position tolerance to the given one,...
Definition: FitPeak.cpp:237
void setFitWindow(double leftwindow, double rightwindow)
Set fit range.
Definition: FitPeak.cpp:93
double m_minFitX
Lower boundary of fitting range.
Definition: FitPeak.h:158
bool m_peakWidthSet
Flag whether the peak width is set.
Definition: FitPeak.h:141
std::string m_costFunction
Cost function.
Definition: FitPeak.h:193
double fitFunctionMD(const API::IFunction_sptr &fitfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, std::vector< double > vec_xmin, std::vector< double > vec_xmax)
Fit function in multiple-domain.
Definition: FitPeak.cpp:707
double m_maxFitX
Upper boundary of fitting range.
Definition: FitPeak.h:160
size_t i_maxPeakX
index of m_maxPeakX
Definition: FitPeak.h:173
API::MatrixWorkspace_sptr genFitWindowWS()
Generate a partial workspace at fit window.
Definition: FitPeak.cpp:341
double m_minPeakX
peak left boundary (client-defined)
Definition: FitPeak.h:167
std::vector< double > m_vecFWHM
Definition: FitPeak.h:195
bool m_fitMethodSet
Flag to show whether fitting parameters are set.
Definition: FitPeak.h:137
double fitCompositeFunction(const API::IPeakFunction_sptr &peakfunc, const API::IBackgroundFunction_sptr &bkgdfunc, const API::MatrixWorkspace_sptr &dataws, size_t wsindex, double startx, double endx)
Fit peak and background composite function.
Definition: FitPeak.cpp:791
size_t i_minFitX
index of m_minFitX
Definition: FitPeak.h:162
size_t i_minPeakX
index of m_minPeakX
Definition: FitPeak.h:171
void setWorskpace(const API::MatrixWorkspace_sptr &dataws, size_t wsindex)
Set workspaces.
Definition: FitPeak.cpp:65
void removeBackground(const API::MatrixWorkspace_sptr &purePeakWS)
remove background
Definition: FitPeak.cpp:414
double getFitCostFunctionValue()
Get cost function value from fitting.
Definition: FitPeak.cpp:974
FitPeak : Fit a single peak.
Definition: FitPeak.h:218
size_t i_maxFitX
Vector index of m_maxFitX.
Definition: FitPeak.h:321
std::string m_minimizer
Minimizer.
Definition: FitPeak.h:361
bool m_fitBkgdFirst
fitting strategy
Definition: FitPeak.h:328
double m_minFitX
Minimum fit position.
Definition: FitPeak.h:310
double m_maxFitX
Maximum fit position.
Definition: FitPeak.h:312
DataObjects::TableWorkspace_sptr genOutputTableWS(const API::IPeakFunction_sptr &peakfunc, std::map< std::string, double > peakerrormap, const API::IBackgroundFunction_sptr &bkgdfunc, std::map< std::string, double > bkgderrormap)
Generate table workspace.
Definition: FitPeak.cpp:1501
bool m_usePeakPositionTolerance
Use peak position tolerance as a criterial for peak fit.
Definition: FitPeak.h:348
bool m_fitWithStepPeakWidth
Flag about guessed FWHM (pixels)
Definition: FitPeak.h:345
void processProperties()
Process input propeties.
Definition: FitPeak.cpp:1183
void createFunctions()
Create functions.
Definition: FitPeak.cpp:1288
API::MatrixWorkspace_sptr m_dataWS
Input data workspace.
Definition: FitPeak.h:302
double m_maxPeakX
Maximum peak position.
Definition: FitPeak.h:316
std::vector< std::string > m_bkgdParameterNames
Background.
Definition: FitPeak.h:358
std::string parseFunctionTypeFull(const std::string &fullstring, bool &defaultparorder)
Parse peak type from full peak type/parameter names string.
Definition: FitPeak.cpp:1367
std::vector< std::string > addFunctionParameterNames(const std::vector< std::string > &funcnames)
Add function's parameter names after peak function name.
Definition: FitPeak.cpp:1153
void init() override
Declare properties.
Definition: FitPeak.cpp:1015
int m_fwhmFitStep
Step width of tried FWHM.
Definition: FitPeak.h:343
void exec() override
Declare properties.
Definition: FitPeak.cpp:1112
int m_maxGuessedPeakWidth
Maximum guessed peak width (pixels)
Definition: FitPeak.h:341
API::IBackgroundFunction_sptr m_bkgdFunc
Background function.
Definition: FitPeak.h:308
double m_minPeakX
Minimum peak position.
Definition: FitPeak.h:314
double m_peakPositionTolerance
Tolerance on peak positions as criteria.
Definition: FitPeak.h:350
double m_finalGoodnessValue
Final.
Definition: FitPeak.h:375
API::IPeakFunction_sptr m_peakFunc
Peak function.
Definition: FitPeak.h:306
std::string m_costFunction
Definition: FitPeak.h:381
bool m_outputRawParams
output option
Definition: FitPeak.h:331
double m_userGuessedFWHM
User guessed FWHM.
Definition: FitPeak.h:334
int m_minGuessedPeakWidth
Minimum guessed peak width (pixels)
Definition: FitPeak.h:339
size_t i_minFitX
Vector index of m_minFitX.
Definition: FitPeak.h:319
void setupOutput(const std::map< std::string, double > &m_fitErrorPeakFunc, const std::map< std::string, double > &m_fitErrorBkgdFunc)
Set up the output workspaces.
Definition: FitPeak.cpp:1405
std::vector< std::string > m_peakParameterNames
Peak.
Definition: FitPeak.h:356
void prescreenInputData()
Check the input properties and functions.
Definition: FitPeak.cpp:1386
double m_userPeakCentre
User guessed peak centre.
Definition: FitPeak.h:336
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
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
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
StartsWithValidator is a validator that requires the value of a property to start with one of the str...
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< IBackgroundFunction > IBackgroundFunction_sptr
std::shared_ptr< IPeakFunction > IPeakFunction_sptr
std::shared_ptr< const IPeakFunction > IPeakFunction_const_sptr
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< const IFunction > IFunction_const_sptr
shared pointer to the function base class (const version)
Definition: IFunction.h:734
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< CompositeFunction > CompositeFunction_sptr
shared pointer to the composite function base class
size_t getIndex(const HistogramData::HistogramX &vecx, double x)
Get an index of a value in a sorted vector.
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< IValidator > IValidator_sptr
A shared_ptr to an IValidator.
Definition: IValidator.h:26
Helper class which provides the Collimation Length for SANS instruments.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
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
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54