Mantid
Loading...
Searching...
No Matches
AsymmetricPearsonVII.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2022 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//----------------------------------------------------------------------
13
14#include <cmath>
15#include <iostream>
16
18
19using namespace CurveFitting;
20using namespace Constraints;
21using namespace std;
22using namespace Kernel;
23using namespace API;
24
25DECLARE_FUNCTION(AsymmetricPearsonVII)
26
28
30 declareParameter("PeakHeight", 1.0, "Hight of the peak");
31 declareParameter("PeakCentre", 0.0, "Location of the peak");
32 declareParameter("Width", 0.1, "Full width at half maximum");
33 declareParameter("LeftShape", 1.0, "Left shape parameter");
34 declareParameter("RightShape", 1.0, "Right shape parameter");
35
36 // constrain shape parameters to the suggested range of values
37 auto leftShapeConstraint = std::make_unique<BoundaryConstraint>(this, "LeftShape", 0.0, true);
38 addConstraint(std::move(leftShapeConstraint));
39 auto rightShapeConstraint = std::make_unique<BoundaryConstraint>(this, "RightShape", 0.0, true);
40 addConstraint(std::move(rightShapeConstraint));
41}
42
43void AsymmetricPearsonVII::functionLocal(double *out, const double *xValues, const size_t nData) const {
44 const double peak_height = getParameter("PeakHeight");
45 const double peak_centre = getParameter("PeakCentre");
46 const double weight = 1.0 / getParameter("Width");
47 const double ml = getParameter("LeftShape");
48 const double mr = getParameter("RightShape");
49
50 for (size_t i = 0; i < nData; ++i) {
51 double offset = xValues[i] - peak_centre;
52
53 if (xValues[i] <= peak_centre) {
54 // left_part
55 if (ml == 0)
56 out[i] = getPearsonVIILimitmEq0(peak_height);
57 else
58 out[i] = getPearsonVII(peak_height, offset, weight, ml);
59 } else {
60 // right_part
61 if (mr == 0)
62 out[i] = getPearsonVIILimitmEq0(peak_height);
63 else
64 out[i] = getPearsonVII(peak_height, offset, weight, mr);
65 }
66 }
67}
68
69void AsymmetricPearsonVII::functionDerivLocal(Jacobian *out, const double *xValues, const size_t nData) {
70 const double peak_height = getParameter("PeakHeight");
71 const double peak_centre = getParameter("PeakCentre");
72 const double weight = 1.0 / getParameter("Width");
73 const double ml = getParameter("LeftShape");
74 const double mr = getParameter("RightShape");
75
76 // derivatives
77 for (size_t i = 0; i < nData; ++i) {
78 double offset = xValues[i] - peak_centre;
79
80 if (xValues[i] <= peak_centre) {
81 // left_part
82 if (ml == 0) {
86 out->set(i, 3, getPearsonVIIDerivWRTmLimitmEq0(peak_height, offset, weight));
87 } else {
88 out->set(i, 0, getPearsonVIIDerivWRTh(offset, weight, ml));
89 out->set(i, 1, getPearsonVIIDerivWRTc(peak_height, offset, weight, ml));
90 out->set(i, 2, getPearsonVIIDerivWRTw(peak_height, offset, weight, ml));
91 out->set(i, 3, getPearsonVIIDerivWRTm(peak_height, offset, weight, ml));
92 }
93 out->set(i, 4, 0.0);
94 } else {
95 // right_part
96 if (mr == 0) {
100 out->set(i, 4, getPearsonVIIDerivWRTmLimitmEq0(peak_height, offset, weight));
101 } else {
102 out->set(i, 0, getPearsonVIIDerivWRTh(offset, weight, mr));
103 out->set(i, 1, getPearsonVIIDerivWRTc(peak_height, offset, weight, mr));
104 out->set(i, 2, getPearsonVIIDerivWRTw(peak_height, offset, weight, mr));
105 out->set(i, 4, getPearsonVIIDerivWRTm(peak_height, offset, weight, mr));
106 }
107 out->set(i, 3, 0.0);
108 }
109 }
110}
111
113 if (!isActive(i)) {
114 throw std::runtime_error("Attempt to use an inactive parameter");
115 }
116 if (parameterName(i) == "Width")
117 setParameter(i, 1.0 / value, false);
118 else
119 setParameter(i, value, false);
120}
121
123 if (!isActive(i)) {
124 throw std::runtime_error("Attempt to use an inactive parameter");
125 }
126 if (parameterName(i) == "Width")
127 return 1.0 / getParameter(i);
128 else
129 return getParameter(i);
130}
131
132double AsymmetricPearsonVII::height() const { return getParameter("PeakHeight"); }
133double AsymmetricPearsonVII::centre() const { return getParameter("PeakCentre"); }
134double AsymmetricPearsonVII::fwhm() const { return getParameter("Width"); }
135double AsymmetricPearsonVII::leftShape() const { return getParameter("LeftShape"); }
136double AsymmetricPearsonVII::rightShape() const { return getParameter("RightShape"); }
137
138void AsymmetricPearsonVII::setCentre(const double c) { setParameter("PeakCentre", c); }
139void AsymmetricPearsonVII::setHeight(const double h) { setParameter("PeakHeight", h); }
140void AsymmetricPearsonVII::setFwhm(const double w) { setParameter("Width", w); }
141void AsymmetricPearsonVII::setLeftShape(const double ml) { setParameter("LeftShape", ml); }
142void AsymmetricPearsonVII::setRightShape(const double mr) { setParameter("RightShape", mr); }
143
144double AsymmetricPearsonVII::getPearsonVII(double peak_height, double offset, double weight, double m) const {
145 return peak_height / pow(denominator_function(offset, weight, m), m);
146}
147double AsymmetricPearsonVII::getPearsonVIILimitmEq0(double peak_height) const { return peak_height / 2.0; }
148
149double AsymmetricPearsonVII::getPearsonVIIDerivWRTh(double offset, double weight, double m) const {
150 return 1.0 / pow(denominator_function(offset, weight, m), m);
151}
153
154double AsymmetricPearsonVII::getPearsonVIIDerivWRTc(double peak_height, double offset, double weight, double m) const {
155 return derivative_function(peak_height, offset, weight, m) * weight;
156}
158
159double AsymmetricPearsonVII::getPearsonVIIDerivWRTw(double peak_height, double offset, double weight, double m) const {
160 return -derivative_function(peak_height, offset, weight, m) * offset;
161}
163
164double AsymmetricPearsonVII::getPearsonVIIDerivWRTm(double peak_height, double offset, double weight, double m) const {
165 return m_derivative_function(peak_height, offset, weight, m);
166}
167double AsymmetricPearsonVII::getPearsonVIIDerivWRTmLimitmEq0(double peak_height, double offset, double weight) const {
168 double offset_sq = offset * offset;
169 double weight_sq = weight * weight;
170 return -peak_height / 2.0 * log(4.0 * offset_sq * weight_sq);
171}
172
173double denominator_function(double offset, double weight, double m) {
174 double offset_sq = offset * offset;
175 double weight_sq = weight * weight;
176 return 1.0 + 4.0 * offset_sq * (pow(2, 1.0 / m) - 1.0) * weight_sq;
177}
178
179double derivative_function(double peak_height, double offset, double weight, double m) {
180 return 8.0 * peak_height * m * offset * (pow(2, 1.0 / m) - 1.0) * weight /
181 pow(denominator_function(offset, weight, m), 1.0 + m);
182}
183
184double m_derivative_function(double peak_height, double offset, double weight, double m) {
185 double offset_sq = offset * offset;
186 double weight_sq = weight * weight;
187 return peak_height *
188 (4.0 * offset_sq * weight_sq * pow(2, 1.0 / m) * log(2.0) / m / denominator_function(offset, weight, m) -
189 log(denominator_function(offset, weight, m))) /
190 pow(denominator_function(offset, weight, m), m);
191}
192
193} // 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.
bool isActive(size_t i) const
Check if an active parameter i is actually active.
Definition: IFunction.cpp:160
virtual void addConstraint(std::unique_ptr< IConstraint > ic)
Add a constraint to function.
Definition: IFunction.cpp:376
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.
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.
Provides an implementation of the asymmetric PearsonVII function (sometimes it is also referred to as...
double fwhm() const override
Returns the peak FWHM.
void setCentre(const double newCentre) override
Sets the parameters such that centre == c.
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Derivative evaluation method. Default is to calculate numerically.
void setActiveParameter(size_t i, double value) override
Set new value of i-th active parameter.
double getPearsonVIIDerivWRTh(double offset, double weight, double m) const
double getPearsonVII(double peak_height, double offset, double weight, double m) const
double getPearsonVIIDerivWRTmLimitmEq0(double peak_height, double offset, double weight) const
void setFwhm(const double newFwhm) override
Sets the parameters such that FWHM = w.
double getPearsonVIIDerivWRTm(double peak_height, double offset, double weight, double m) const
double activeParameter(size_t i) const override
Value of i-th active parameter.
void setHeight(const double newHight) override
Sets the parameters such that height == h.
void init() override
Override IFunction base class method, which declares function parameters.
double height() const override
Override IPeakFunction base class methods.
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Function evaluation method to be implemented in the inherited classes.
double getPearsonVIIDerivWRTw(double peak_height, double offset, double weight, double m) const
double getPearsonVIIDerivWRTc(double peak_height, double offset, double weight, double m) const
double centre() const override
Returns the centre of the function, which may be something as simple as the centre of the fitting ran...
double m_derivative_function(double peak_height, double offset_sq, double weight_sq, double m)
double denominator_function(double offset_sq, double weight_sq, double m)
double derivative_function(double peak_height, double offset, double weight, double m)
STL namespace.