Mantid
Loading...
Searching...
No Matches
PeakFunctionIntegrator.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 +
8
10#include "gsl/gsl_errno.h"
11
12namespace Mantid::API {
13
20PeakFunctionIntegrator::PeakFunctionIntegrator(double requiredRelativePrecision)
21 : m_integrationWorkspace(gsl_integration_workspace_alloc(1000)), m_relativePrecision(requiredRelativePrecision) {
22 /* Error handling is disabled, so error-codes are associated to the
23 * integration result
24 * and have to be checked.
25 */
26 gsl_set_error_handler_off();
27}
28
30
37
41
50 IntegrationResult result;
51
52 gsl_function f = getGSLFunction(peakFunction);
53
54 result.errorCode =
55 gsl_integration_qagi(&f, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
56 result.success = (result.errorCode == GSL_SUCCESS);
57 result.intervals = m_integrationWorkspace->size;
58
59 return result;
60}
61
69 double lowerLimit) const {
70 IntegrationResult result;
71
72 gsl_function f = getGSLFunction(peakFunction);
73
74 result.errorCode = gsl_integration_qagiu(&f, lowerLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace,
75 &result.result, &result.error);
76 result.success = (result.errorCode == GSL_SUCCESS);
77 result.intervals = m_integrationWorkspace->size;
78
79 return result;
80}
81
89 double upperLimit) const {
90 IntegrationResult result;
91
92 gsl_function f = getGSLFunction(peakFunction);
93
94 result.errorCode = gsl_integration_qagil(&f, upperLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace,
95 &result.result, &result.error);
96 result.success = (result.errorCode == GSL_SUCCESS);
97 result.intervals = m_integrationWorkspace->size;
98
99 return result;
100}
101
110 double upperLimit) const {
111 IntegrationResult result;
112
113 gsl_function f = getGSLFunction(peakFunction);
114
115 result.errorCode = gsl_integration_qags(&f, lowerLimit, upperLimit, 0, m_relativePrecision, 1000,
116 m_integrationWorkspace, &result.result, &result.error);
117 result.success = (result.errorCode == GSL_SUCCESS);
118 result.intervals = m_integrationWorkspace->size;
119 return result;
120}
121
132double PeakFunctionIntegrator::integrateError(const IPeakFunction &peakFunction, double lowerLimit,
133 double upperLimit) const {
134 const double DBL_EPS = std::numeric_limits<double>::epsilon();
135 size_t nParams = peakFunction.nParams();
136
137 std::vector<double> parameterErrors(nParams);
138 bool hasErrors = false;
139 for (size_t i = 0; i < nParams; ++i) {
140 parameterErrors[i] = peakFunction.getError(i);
141 if (parameterErrors[i] > DBL_EPS)
142 hasErrors = true;
143 }
144 if (!hasErrors)
145 return std::nan("");
146
147 // estimate the partial derivatives of the integrated intensity with respect to the fit parameters
148 std::vector<double> gradients(nParams);
149 auto peakFunctionClone = std::dynamic_pointer_cast<IPeakFunction>(peakFunction.clone());
150 for (size_t i = 0; i < nParams; i++) {
151 double parameterValue = peakFunction.getParameter(i);
152 double parameterError = parameterErrors[i];
153 if (parameterError > DBL_EPS) {
154 double step(parameterError / 2.0);
155 // evaluate the integrated intensity after increasing parameter `i`
156 peakFunctionClone->setParameter(i, parameterValue + step);
157 std::shared_ptr<const IPeakFunction> functionPlus(
158 std::const_pointer_cast<const IPeakFunction>(peakFunctionClone));
159 auto resultPlus = integrate(*functionPlus, lowerLimit, upperLimit);
160 // evaluate the integrated intensity after decreasing parameter `i`
161 peakFunctionClone->setParameter(i, parameterValue - step);
162 auto functionMinus(std::const_pointer_cast<const IPeakFunction>(peakFunctionClone));
163 auto resultMinus = integrate(*functionMinus, lowerLimit, upperLimit);
164 // restore the original value and compute the partial derivative
165 peakFunctionClone->setParameter(i, parameterValue);
166 gradients[i] = (resultPlus.result - resultMinus.result) / parameterError;
167 } else
168 gradients[i] = 0.0; // don't bother to calculate since the parameter's error is negligible
169 }
170
171 double error(0.0);
172 for (size_t i = 0; i < nParams; i++)
173 error += pow(gradients[i] * parameterErrors[i], 2);
174 return sqrt(error);
175}
176
181gsl_function PeakFunctionIntegrator::getGSLFunction(const IPeakFunction &peakFunction) const {
182 gsl_function f;
183 f.function = &Mantid::API::gsl_peak_wrapper;
184 f.params = reinterpret_cast<void *>(&const_cast<IPeakFunction &>(peakFunction));
185
186 return f;
187}
188
189double gsl_peak_wrapper(double x, void *parameters) {
190 auto *peakFunction = reinterpret_cast<IPeakFunction *>(parameters);
191
192 if (!peakFunction) {
193 throw std::runtime_error("Cannot process NULL-pointer in gsl_peak_wrapper.");
194 }
195
196 double y;
197
198 /* For the integration to work properly, functionLocal has to be used instead
199 * of the more general function-method due to the fact that the overriden
200 * function-method
201 * in IPeakFunction cuts off at some point. For slowly decaying peak functions
202 * such as Lorentzians, this is introduces large deviations for integrations
203 * from -Inf to +Inf.
204 */
205 peakFunction->functionLocal(&y, &x, 1);
206
207 return y;
208}
209} // namespace Mantid::API
double error
Definition: IndexPeaks.cpp:133
virtual std::shared_ptr< IFunction > clone() const
Virtual copy constructor.
Definition: IFunction.cpp:109
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
double getError(size_t i) const override
Get the fitting error for a parameter.
size_t nParams() const override
Total number of parameters.
Definition: ParamFunction.h:53
double getParameter(size_t i) const override
Get i-th parameter.
double requiredRelativePrecision() const
Returns the currently set precision.
void setRequiredRelativePrecision(double newPrecision)
This method sets the desired numerical relative precision that's passed on to the GSL integration-rou...
double integrateError(const IPeakFunction &peakFunction, double lowerLimit, double upperLimit) const
Error in the integrated intensity within an continuous interval due to uncertainties in the values of...
gsl_integration_workspace * m_integrationWorkspace
IntegrationResult integrateNegativeInfinity(const IPeakFunction &peakFunction, double upperLimit) const
Integration of peak function on the interval [-Inf, b].
gsl_function getGSLFunction(const IPeakFunction &peakFunction) const
Method that wraps an IPeakFunction for use with GSL functions.
IntegrationResult integratePositiveInfinity(const IPeakFunction &peakFunction, double lowerLimit) const
Integration of peak function on the interval [a, +Inf].
IntegrationResult integrate(const IPeakFunction &peakFunction, double lowerLimit, double upperLimit) const
Integration of peak function on the interval [a, b].
PeakFunctionIntegrator(double requiredRelativePrecision=1e-8)
Constructor with required relative precision argument.
IntegrationResult integrateInfinity(const IPeakFunction &peakFunction) const
Integration of peak function on the interval [-Inf, +Inf].
double MANTID_API_DLL gsl_peak_wrapper(double x, void *parameters)