Mantid
Loading...
Searching...
No Matches
Gaussian.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
12
13#include <cmath>
14#include <numeric>
15
17
18using namespace CurveFitting;
19using namespace Kernel;
20using namespace API;
21
22DECLARE_FUNCTION(Gaussian)
23
24Gaussian::Gaussian() : IPeakFunction(), m_intensityCache(0.0) {}
25
27 declareParameter("Height", 0.0, "Height of peak");
28 declareParameter("PeakCentre", 0.0, "Centre of peak");
29 declareParameter("Sigma", 0.0, "Width parameter");
30}
31
32void Gaussian::functionLocal(double *out, const double *xValues, const size_t nData) const {
33 const double height = getParameter("Height");
34 const double peakCentre = getParameter("PeakCentre");
35 const double weight = pow(1 / getParameter("Sigma"), 2);
36
37 for (size_t i = 0; i < nData; i++) {
38 double diff = xValues[i] - peakCentre;
39 out[i] = height * exp(-0.5 * diff * diff * weight);
40 }
41}
42
43void Gaussian::functionDerivLocal(Jacobian *out, const double *xValues, const size_t nData) {
44 const double height = getParameter("Height");
45 const double peakCentre = getParameter("PeakCentre");
46 const double weight = pow(1 / getParameter("Sigma"), 2);
47
48 for (size_t i = 0; i < nData; i++) {
49 double diff = xValues[i] - peakCentre;
50 double e = exp(-0.5 * diff * diff * weight);
51 out->set(i, 0, e);
52 out->set(i, 1, diff * height * e * weight);
53 out->set(i, 2,
54 -0.5 * diff * diff * height * e); // derivative with respect to weight not sigma
55 }
56}
57
58void Gaussian::setActiveParameter(size_t i, double value) {
59 if (!isActive(i)) {
60 throw std::runtime_error("Attempt to use an inactive parameter");
61 }
62 if (parameterName(i) == "Sigma")
63 setParameter(i, sqrt(fabs(1. / value)), false);
64 else
65 setParameter(i, value, false);
66}
67
68double Gaussian::activeParameter(size_t i) const {
69 if (!isActive(i)) {
70 throw std::runtime_error("Attempt to use an inactive parameter");
71 }
72 if (parameterName(i) == "Sigma")
73 return 1. / pow(getParameter(i), 2);
74 else
75 return getParameter(i);
76}
77
78double Gaussian::centre() const { return getParameter("PeakCentre"); }
79double Gaussian::height() const { return getParameter("Height"); }
80double Gaussian::fwhm() const { return 2.0 * sqrt(2.0 * M_LN2) * getParameter("Sigma"); }
81double Gaussian::intensity() const {
82 auto sigma = getParameter("Sigma");
83 if (sigma == 0.0) {
84 auto height = getParameter("Height");
85 if (std::isfinite(height)) {
87 }
88 } else {
89 m_intensityCache = getParameter("Height") * getParameter("Sigma") * sqrt(2.0 * M_PI);
90 }
91 return m_intensityCache;
92}
94 const double heightError = getError("Height");
95 const double sigmaError = getError("Sigma");
96
97 return intensity() * sqrt(pow(heightError / getParameter("Height"), 2) + pow(sigmaError / getParameter("Sigma"), 2));
98}
99
100void Gaussian::setCentre(const double c) { setParameter("PeakCentre", c); }
101void Gaussian::setHeight(const double h) { setParameter("Height", h); }
102void Gaussian::setFwhm(const double w) { setParameter("Sigma", w / (2.0 * sqrt(2.0 * M_LN2))); }
103void Gaussian::setIntensity(const double i) {
105 auto sigma = getParameter("Sigma");
106 if (sigma == 0.0) {
107 setParameter("Height", i);
108 } else {
109 setParameter("Height", i / (sigma * sqrt(2.0 * M_PI)));
110 }
111}
112
113void Gaussian::fixCentre(bool isDefault) { fixParameter("PeakCentre", isDefault); }
114
115void Gaussian::unfixCentre() { unfixParameter("PeakCentre"); }
116
117void Gaussian::fixIntensity(bool isDefault) {
118 std::string formula = std::to_string(intensity() / sqrt(2.0 * M_PI)) + "/Sigma";
119 tie("Height", formula, isDefault);
120}
121
123
131void Gaussian::histogram1D(double *out, double left, const double *right, const size_t nBins) const {
132
133 double amplitude = intensity();
134 const double peakCentre = getParameter("PeakCentre");
135 const double sigma2 = getParameter("Sigma") * sqrt(2.0);
136
137 auto cumulFun = [sigma2, peakCentre](double x) { return 0.5 * erf((x - peakCentre) / sigma2); };
138 double cLeft = cumulFun(left);
139 for (size_t i = 0; i < nBins; ++i) {
140 double cRight = cumulFun(right[i]);
141 out[i] = amplitude * (cRight - cLeft);
142 cLeft = cRight;
143 }
144}
145
152void Gaussian::histogramDerivative1D(Jacobian *jacobian, double left, const double *right, const size_t nBins) const {
153 const double h = getParameter("Height");
154 const double c = getParameter("PeakCentre");
155 const double s = getParameter("Sigma");
156 const double w = pow(1 / s, 2);
157 const double sw = sqrt(w);
158
159 auto cumulFun = [sw, c](double x) { return sqrt(M_PI / 2) / sw * erf(sw / sqrt(2.0) * (x - c)); };
160 auto fun = [w, c](double x) { return exp(-w / 2 * pow(x - c, 2)); };
161
162 double xl = left;
163 double fLeft = fun(left);
164 double cLeft = cumulFun(left);
165 const double h_over_2w = h / (2 * w);
166 for (size_t i = 0; i < nBins; ++i) {
167 double xr = right[i];
168 double fRight = fun(xr);
169 double cRight = cumulFun(xr);
170 jacobian->set(i, 0, cRight - cLeft); // height
171 jacobian->set(i, 1, -h * (fRight - fLeft)); // centre
172 jacobian->set(i, 2,
173 h_over_2w * ((xr - c) * fRight - (xl - c) * fLeft + cLeft - cRight)); // weight
174 fLeft = fRight;
175 cLeft = cRight;
176 xl = xr;
177 }
178}
179
180} // namespace Mantid::CurveFitting::Functions
double value
The value of the point.
Definition: FitMW.cpp:51
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
double sigma
Definition: GetAllEi.cpp:156
double left
Definition: LineProfile.cpp:80
double right
Definition: LineProfile.cpp:81
#define fabs(x)
Definition: Matrix.cpp:22
bool isActive(size_t i) const
Check if an active parameter i is actually active.
Definition: IFunction.cpp:160
void unfixParameter(const std::string &name)
Free a parameter.
Definition: IFunction.cpp:1522
virtual void tie(const std::string &parName, const std::string &expr, bool isDefault=false)
Tie a parameter to other parameters (or a constant)
Definition: IFunction.cpp:214
void fixParameter(const std::string &name, bool isDefault=false)
Fix a parameter.
Definition: IFunction.cpp:1515
virtual void removeTie(const std::string &parName)
Removes the tie off a parameter.
Definition: IFunction.cpp:254
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
void setParameter(const std::string &name, const double &value, bool explicitlySet=true) override
Override parent so that we may bust the cache when a parameter is set.
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
virtual void set(size_t iY, size_t iP, double value)=0
Set a value to a Jacobian matrix element.
double getError(size_t i) const override
Get the fitting error for a parameter.
std::string parameterName(size_t i) const override
Returns the name of parameter i.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
double getParameter(size_t i) const override
Get i-th parameter.
Provide gaussian peak shape function interface to IPeakFunction.
Definition: Gaussian.h:36
void fixIntensity(bool isDefault=false) override
Fix a parameter or set up a tie such that value returned by intensity() is constant during fitting.
Definition: Gaussian.cpp:117
void setHeight(const double h) override
Sets the parameters such that height == h.
Definition: Gaussian.cpp:101
double centre() const override
overwrite IPeakFunction base class methods
Definition: Gaussian.cpp:78
void histogramDerivative1D(API::Jacobian *jacobian, double left, const double *right, const size_t nBins) const override
Devivatives of the histogram.
Definition: Gaussian.cpp:152
void setFwhm(const double w) override
Sets the parameters such that FWHM = w.
Definition: Gaussian.cpp:102
void setIntensity(const double i) override
Sets the integral intensity of the peak.
Definition: Gaussian.cpp:103
void setCentre(const double c) override
Sets the parameters such that centre == c.
Definition: Gaussian.cpp:100
double intensity() const override
Returns the integral intensity of the peak.
Definition: Gaussian.cpp:81
double fwhm() const override
Returns the peak FWHM.
Definition: Gaussian.cpp:80
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Function evaluation method to be implemented in the inherited classes.
Definition: Gaussian.cpp:32
double activeParameter(size_t i) const override
Value of i-th active parameter.
Definition: Gaussian.cpp:68
double intensityError() const override
Error in the integrated intensity of the peak due to uncertainties in the values of the fit parameter...
Definition: Gaussian.cpp:93
void unfixIntensity() override
Free the intensity parameter.
Definition: Gaussian.cpp:122
void init() override
overwrite IFunction base class method, which declare function parameters
Definition: Gaussian.cpp:26
void unfixCentre() override
Free the centre parameter.
Definition: Gaussian.cpp:115
void histogram1D(double *out, double left, const double *right, const size_t nBins) const override
Calculate histogram data.
Definition: Gaussian.cpp:131
double height() const override
Returns the height of the function.
Definition: Gaussian.cpp:79
double m_intensityCache
Intensity cache to help recover form Sigma==0 situation.
Definition: Gaussian.h:72
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Derivative evaluation method. Default is to calculate numerically.
Definition: Gaussian.cpp:43
void fixCentre(bool isDefault=false) override
Fix a parameter or set up a tie such that value returned by centre() is constant during fitting.
Definition: Gaussian.cpp:113
void setActiveParameter(size_t i, double value) override
Set new value of i-th active parameter.
Definition: Gaussian.cpp:58
std::string to_string(const wide_integer< Bits, Signed > &n)