Mantid
Loading...
Searching...
No Matches
PseudoVoigt.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 +
10
11#include <iostream>
12
13#include <cmath>
14
15namespace {
16
17const double LN2_DIV_PI = 2. * std::sqrt(M_LN2 / M_PI);
18constexpr double LN2_M4 = 4 * M_LN2;
19
24inline double cal_ag(const double gamma) { return LN2_DIV_PI / gamma; }
25
30inline constexpr double cal_bg(const double gamma) { return LN2_M4 / (gamma * gamma); }
31
38inline double cal_gaussian(const double ag, const double bg, const double xdiffsq) {
39 return ag * std::exp(-bg * xdiffsq);
40}
41
48inline constexpr double cal_lorentzian(const double gamma_div_2, const double gammasq_div_4, const double xdiffsq) {
49 return gamma_div_2 / (xdiffsq + gammasq_div_4) / M_PI;
50}
51
52} // namespace
53
55
56using namespace CurveFitting;
57using namespace Constraints;
58using namespace std;
59using namespace API;
60
61namespace {
63Kernel::Logger g_log("PseudoVoigt");
64} // namespace
65
66DECLARE_FUNCTION(PseudoVoigt)
67
68
72void PseudoVoigt::init() {
73 declareParameter("Mixing", 0.5);
74 declareParameter("Intensity", 1.);
75 declareParameter("PeakCentre", 0.);
76 declareParameter("FWHM", 1.);
77
78 auto mixingConstraint = std::make_unique<BoundaryConstraint>(this, "Mixing", 0.0, 1.0, true);
79 mixingConstraint->setPenaltyFactor(1e9);
80 addConstraint(std::move(mixingConstraint));
81
82 auto fwhm_constraint = std::make_unique<BoundaryConstraint>(this, "FWHM", 1.E-20, true);
83 fwhm_constraint->setPenaltyFactor(1e9);
84 addConstraint(std::move(fwhm_constraint));
85
86 // init the peak height setup parameters
87 m_height = 1.; // peak height set by user
88
89 // parameter set history: all start from an arbitary out of boundary value
90 // (100 as easy). if set, order set to 0
91 m_set_history_distances.resize(4, 128);
92 // 0: mixing, 1: intensity, 2: fwhm, 3: height
93}
94
100void PseudoVoigt::functionLocal(double *out, const double *xValues, const size_t nData) const {
101 // read values
102 const double peak_intensity = getParameter("Intensity");
103 const double x0 = getParameter("PeakCentre");
104 const double gamma = fabs(getParameter("FWHM"));
105 const double gFraction = getParameter("Mixing");
106
107 if (gamma < 1.E-20)
108 throw std::runtime_error("Pseudo-voigt has FWHM as 0. It will generate "
109 "infinity at center in the Lorentzian part.");
110
111 const double lFraction = 1.0 - gFraction;
112
113 // calculate constants
114 const double a_g = cal_ag(gamma);
115 const double b_g = cal_bg(gamma);
116 const double gamma_div_2 = 0.5 * gamma;
117 const double gammasq_div_4 = gamma_div_2 * gamma_div_2;
118
119 for (size_t i = 0; i < nData; ++i) {
120 double xDiffSquared = (xValues[i] - x0) * (xValues[i] - x0);
121 out[i] = peak_intensity * (gFraction * cal_gaussian(a_g, b_g, xDiffSquared) +
122 lFraction * cal_lorentzian(gamma_div_2, gammasq_div_4, xDiffSquared));
123 }
124}
125
131void PseudoVoigt::functionDerivLocal(Jacobian *out, const double *xValues, const size_t nData) {
132 // get parameters
133 double peak_intensity = getParameter("Intensity");
134 double x0 = getParameter("PeakCentre");
135 double gamma = getParameter("FWHM");
136 double gFraction = getParameter("Mixing");
137 double lFraction = 1.0 - gFraction;
138
139 // calcualte constants
140 const double a_g = cal_ag(gamma);
141 const double b_g = cal_bg(gamma);
142 const double gamma_div_2 = 0.5 * gamma;
143 const double gammasq_div_4 = gamma_div_2 * gamma_div_2;
144
145 // derivatives
146 for (size_t i = 0; i < nData; ++i) {
147 double xDiff = (xValues[i] - x0);
148 double xDiffSquared = xDiff * xDiff;
149
150 // calculate normalized Gaussian and Lorentzian
151 double gaussian_term = cal_gaussian(a_g, b_g, xDiffSquared);
152 double lorentzian_term = cal_lorentzian(gamma_div_2, gammasq_div_4, xDiffSquared);
153
154 // derivative to mixing/eta (0-th)
155 out->set(i, 0, peak_intensity * (gaussian_term - lorentzian_term));
156
157 // derivative on intensity (1-st)
158 out->set(i, 1, gFraction * gaussian_term + lFraction * lorentzian_term);
159
160 // peak center: x0
161 const double derive_g_x0 = 2. * b_g * xDiff * gaussian_term;
162 const double derive_l_x0 = 4. * M_PI * xDiff / gamma * lorentzian_term * lorentzian_term;
163 const double deriv_x0 = peak_intensity * (gFraction * derive_g_x0 + lFraction * derive_l_x0);
164 out->set(i, 2, deriv_x0);
165
166 // peak width: gamma or H
167 const double t1 = -1. / gamma * gaussian_term;
168 const double t2 = 2. * b_g * xDiffSquared * gaussian_term / gamma;
169 const double t3 = lorentzian_term / gamma;
170 const double t4 = -M_PI * lorentzian_term * lorentzian_term;
171 const double derive_gamma = peak_intensity * (gFraction * (t1 + t2) + lFraction * (t3 + t4));
172 out->set(i, 3, derive_gamma);
173 }
174}
175
177void PseudoVoigt::setParameter(const std::string &name, const double &value, bool explicitlySet) {
178 g_log.debug() << "[PV] Set " << name << " as " << value << "\n";
180}
181
199void PseudoVoigt::setParameter(size_t i, const double &value, bool explicitlySet) {
200 API::IPeakFunction::setParameter(i, value, explicitlySet);
201
202 g_log.debug() << "[PV] Set " << i << "-th parameter with value " << value << "\n";
203
204 // explicitly set means that there is a chance that some parameter shall be
205 // re-calculated
206 // peak center shall be excluded as it does nothing to do with height, H,
207 // intensity and mixing
208 if (explicitlySet && i != 2) {
209 g_log.debug() << "Update set history for " << i << "-th parameter"
210 << "\n";
213 }
214
215 return;
216}
217
219 // peak center
220 size_t to_calculate_index = get_parameter_to_calculate_from_set();
221 g_log.debug() << "Time to calculate " << to_calculate_index << "-th parameter"
222 << "\n";
223
224 bool some_value_updated{true};
225 if (to_calculate_index == 0) {
226 // calculate mixings
227 double peak_intensity = getParameter(1);
228 double gamma = getParameter(3);
229 double mixing = (m_height * 0.5 * M_PI * gamma / peak_intensity - 1) / (sqrt(M_PI * M_LN2) - 1);
230 if (mixing < 0) {
231 g_log.debug() << "Calculated mixing (eta) = " << mixing << " < 0. Set to 0 instead"
232 << "\n";
233 mixing = 0;
234 } else if (mixing > 1) {
235 g_log.debug() << "Calculated mixing (eta) = " << mixing << " > 1. Set to 1 instead"
236 << "\n";
237 mixing = 1;
238 }
239
240 setParameter(0, mixing,
241 false); // no explicitly set: won't cause further re-calcualting
242 } else if (to_calculate_index == 1) {
243 // calculate intensity
244 double eta = getParameter(0);
245 double gamma = getParameter(3);
246 double intensity = m_height / 2. / (1 + (sqrt(M_PI * M_LN2) - 1) * eta) * (M_PI * gamma);
247 setParameter(1, intensity, false);
248 g_log.debug() << "Estimate peak intensity = " << intensity << "\n";
249 } else if (to_calculate_index == 3) {
250 // calculate peak width
251 double eta = getParameter(0);
252 double intensity = getParameter(1);
253
254 g_log.debug() << "Intensity = " << intensity << ", height = " << m_height << ", mixing = " << eta << "\n";
255 double gamma = intensity * 2 * (1 + (sqrt(M_PI * M_LN2) - 1) * eta) / (M_PI * m_height);
256 if (gamma < 1.E-10) {
257 g_log.debug() << "Peak width (H or Gamma) = " << gamma << ". Set 1.E-2 instead"
258 << "\n";
259 }
260 setParameter(3, gamma, false);
261 } else {
262 // no value updated
263 some_value_updated = false;
264 }
265
266 return some_value_updated;
267}
268
274void PseudoVoigt::update_set_history(size_t set_index) {
275 if (set_index > 3)
276 throw std::runtime_error("Parameter set up history index out of boundary");
277
278 size_t prev_distance = m_set_history_distances[set_index];
279 for (size_t i = 0; i < 4; ++i) {
280 // no operation on any never-set or out-dated index
281 if (m_set_history_distances[i] >= 3)
282 continue;
283
284 // only increase the historoy distance by 1 in case
285 // the other parameters with those with distance less than the previous one
286 // i.e., out of date or never been set
287 if (m_set_history_distances[i] < prev_distance) {
289 // if out of range, don't worry about it
290 }
291 }
292
293 // set the current one to 0 (most recent)
294 m_set_history_distances[set_index] = 0;
295}
296
304 // go through setting-history-distance
305 size_t index_largest_distance{0};
306 size_t largest_distance{0};
307 size_t num_been_set{0};
308
309 for (size_t i = 0; i < m_set_history_distances.size(); ++i) {
310 g_log.debug() << "distance[" << i << "] = " << m_set_history_distances[i] << "\n";
311 if (m_set_history_distances[i] >= largest_distance) {
312 // for those to set for: any number counts
313 index_largest_distance = i;
314 largest_distance = m_set_history_distances[i];
315 }
316 // only those been set and not out of date counted
317 if (m_set_history_distances[i] < 3) {
318 ++num_been_set;
319 }
320 }
321
322 // signal for no enough index been set up
323 if (num_been_set < 3)
324 index_largest_distance = 128;
325
326 return index_largest_distance;
327}
328
332double PseudoVoigt::height() const {
333 double peak_intensity = getParameter("Intensity");
334 double gamma = getParameter("FWHM");
335 double eta = getParameter("Mixing");
336 const double height = 2. * peak_intensity * (1 + (sqrt(M_PI * M_LN2) - 1) * eta) / (M_PI * gamma);
337
338 return height;
339}
340
345void PseudoVoigt::setHeight(const double h) {
346 // set height
347 double prev_height = m_height;
348 m_height = h;
349
350 // update set history
352 g_log.debug() << "PV: set height = " << m_height << "\n";
353
354 // update other parameters' value
355 bool updated_any_value = estimate_parameter_value();
356 if (!updated_any_value) {
357 // force to update peak intensity
358 double old_intensity = getParameter("Intensity");
359 // This is what Roman suggests!
360 double new_intensity{1.};
361 if (prev_height > 0)
362 new_intensity = old_intensity * m_height / prev_height;
363 else
364 new_intensity = m_height; // this is a better than nothing estimation
365 setParameter("Intensity", new_intensity,
366 false); // no need to update other parameters
367 }
368}
369
373void PseudoVoigt::setFwhm(const double w) { setParameter("FWHM", w, true); }
374
375} // 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.
#define fabs(x)
Definition: Matrix.cpp:22
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 getParameter(size_t i) const override
Get i-th parameter.
void setFwhm(const double w) override
set FWHM
void functionLocal(double *out, const double *xValues, const size_t nData) const override
calculate pseudo voigt
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
calcualte derivative analytically
void setParameter(size_t i, const double &value, bool explicitlySet=true) override
Set i-th parameter.
void setHeight(const double h) override
set non-native parameter peak height
double height() const override
calculate effective parameter: height
size_t get_parameter_to_calculate_from_set()
get the parameter (by index) to calculate according to parameter set history
void update_set_history(size_t set_index)
The purpose of this is to track which user-specified parameter is outdated and to be calculated.
double intensity() const override
Returns the integral intensity of the peak.
Definition: PseudoVoigt.h:39
std::string name() const override
Returns the function's name.
Definition: PseudoVoigt.h:49
std::vector< size_t > m_set_history_distances
historty of the order to be set
Definition: PseudoVoigt.h:66
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
Kernel::Logger g_log("ExperimentInfo")
static logger object
STL namespace.