10#include "gsl/gsl_errno.h"
21 : m_integrationWorkspace(gsl_integration_workspace_alloc(1000)), m_relativePrecision(requiredRelativePrecision) {
26 gsl_set_error_handler_off();
69 double lowerLimit)
const {
89 double upperLimit)
const {
110 double upperLimit)
const {
133 double upperLimit)
const {
134 const double DBL_EPS = std::numeric_limits<double>::epsilon();
135 size_t nParams = peakFunction.
nParams();
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)
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++) {
152 double parameterError = parameterErrors[i];
153 if (parameterError > DBL_EPS) {
154 double step(parameterError / 2.0);
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);
161 peakFunctionClone->setParameter(i, parameterValue - step);
162 auto functionMinus(std::const_pointer_cast<const IPeakFunction>(peakFunctionClone));
163 auto resultMinus =
integrate(*functionMinus, lowerLimit, upperLimit);
165 peakFunctionClone->setParameter(i, parameterValue);
166 gradients[i] = (resultPlus.result - resultMinus.result) / parameterError;
172 for (
size_t i = 0; i < nParams; i++)
173 error += pow(gradients[i] * parameterErrors[i], 2);
184 f.params =
reinterpret_cast<void *
>(&
const_cast<IPeakFunction &
>(peakFunction));
190 auto *peakFunction =
reinterpret_cast<IPeakFunction *
>(parameters);
193 throw std::runtime_error(
"Cannot process NULL-pointer in gsl_peak_wrapper.");
205 peakFunction->functionLocal(&
y, &
x, 1);
virtual std::shared_ptr< IFunction > clone() const
Virtual copy constructor.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
double getError(size_t i) const override
Get the fitting error for a parameter.
size_t nParams() const override
Total number of parameters.
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...
virtual ~PeakFunctionIntegrator()
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.
double m_relativePrecision
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)