PeakFunctionIntegrator.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAPI/PeakFunctionIntegrator.h"
 
#include "MantidAPI/FunctionDomain1D.h"
#include "gsl/gsl_errno.h"
 
namespace Mantid::API {
 
PeakFunctionIntegrator::PeakFunctionIntegrator(double requiredRelativePrecision)
    : m_integrationWorkspace(gsl_integration_workspace_alloc(1000)), m_relativePrecision(requiredRelativePrecision) {
  /* Error handling is disabled, so error-codes are associated to the
   * integration result
   * and have to be checked.
   */
  gsl_set_error_handler_off();
}
 
PeakFunctionIntegrator::~PeakFunctionIntegrator() { gsl_integration_workspace_free(m_integrationWorkspace); }
 
void PeakFunctionIntegrator::setRequiredRelativePrecision(double newPrecision) { m_relativePrecision = newPrecision; }
 
double PeakFunctionIntegrator::requiredRelativePrecision() const { return m_relativePrecision; }
 
IntegrationResult PeakFunctionIntegrator::integrateInfinity(const IPeakFunction &peakFunction) const {
  IntegrationResult result;
 
  gsl_function f = getGSLFunction(peakFunction);
 
  result.errorCode =
      gsl_integration_qagi(&f, 0, m_relativePrecision, 1000, m_integrationWorkspace, &result.result, &result.error);
  result.success = (result.errorCode == GSL_SUCCESS);
  result.intervals = m_integrationWorkspace->size;
 
  return result;
}
 
IntegrationResult PeakFunctionIntegrator::integratePositiveInfinity(const IPeakFunction &peakFunction,
                                                                    double lowerLimit) const {
  IntegrationResult result;
 
  gsl_function f = getGSLFunction(peakFunction);
 
  result.errorCode = gsl_integration_qagiu(&f, lowerLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace,
                                           &result.result, &result.error);
  result.success = (result.errorCode == GSL_SUCCESS);
  result.intervals = m_integrationWorkspace->size;
 
  return result;
}
 
IntegrationResult PeakFunctionIntegrator::integrateNegativeInfinity(const IPeakFunction &peakFunction,
                                                                    double upperLimit) const {
  IntegrationResult result;
 
  gsl_function f = getGSLFunction(peakFunction);
 
  result.errorCode = gsl_integration_qagil(&f, upperLimit, 0, m_relativePrecision, 1000, m_integrationWorkspace,
                                           &result.result, &result.error);
  result.success = (result.errorCode == GSL_SUCCESS);
  result.intervals = m_integrationWorkspace->size;
 
  return result;
}
 
IntegrationResult PeakFunctionIntegrator::integrate(const IPeakFunction &peakFunction, double lowerLimit,
                                                    double upperLimit) const {
  IntegrationResult result;
 
  gsl_function f = getGSLFunction(peakFunction);
 
  result.errorCode = gsl_integration_qags(&f, lowerLimit, upperLimit, 0, m_relativePrecision, 1000,
                                          m_integrationWorkspace, &result.result, &result.error);
  result.success = (result.errorCode == GSL_SUCCESS);
  result.intervals = m_integrationWorkspace->size;
  return result;
}
 
double PeakFunctionIntegrator::integrateError(const IPeakFunction &peakFunction, double lowerLimit,
                                              double upperLimit) const {
  const double DBL_EPS = std::numeric_limits<double>::epsilon();
  size_t nParams = peakFunction.nParams();
 
  std::vector<double> parameterErrors(nParams);
  bool hasErrors = false;
  for (size_t i = 0; i < nParams; ++i) {
    parameterErrors[i] = peakFunction.getError(i);
    if (parameterErrors[i] > DBL_EPS)
      hasErrors = true;
  }
  if (!hasErrors)
    return std::nan("");
 
  // estimate the partial derivatives of the integrated intensity with respect to the fit parameters
  std::vector<double> gradients(nParams);
  auto peakFunctionClone = std::dynamic_pointer_cast<IPeakFunction>(peakFunction.clone());
  for (size_t i = 0; i < nParams; i++) {
    double parameterValue = peakFunction.getParameter(i);
    double parameterError = parameterErrors[i];
    if (parameterError > DBL_EPS) {
      double step(parameterError / 2.0);
      // evaluate the integrated intensity after increasing parameter `i`
      peakFunctionClone->setParameter(i, parameterValue + step);
      std::shared_ptr<const IPeakFunction> functionPlus(
          std::const_pointer_cast<const IPeakFunction>(peakFunctionClone));
      auto resultPlus = integrate(*functionPlus, lowerLimit, upperLimit);
      // evaluate the integrated intensity after decreasing parameter `i`
      peakFunctionClone->setParameter(i, parameterValue - step);
      auto functionMinus(std::const_pointer_cast<const IPeakFunction>(peakFunctionClone));
      auto resultMinus = integrate(*functionMinus, lowerLimit, upperLimit);
      // restore the original value and compute the partial derivative
      peakFunctionClone->setParameter(i, parameterValue);
      gradients[i] = (resultPlus.result - resultMinus.result) / parameterError;
    } else
      gradients[i] = 0.0; // don't bother to calculate since the parameter's error is negligible
  }
 
  double error(0.0);
  for (size_t i = 0; i < nParams; i++)
    error += pow(gradients[i] * parameterErrors[i], 2);
  return sqrt(error);
}
 
gsl_function PeakFunctionIntegrator::getGSLFunction(const IPeakFunction &peakFunction) const {
  gsl_function f;
  f.function = &Mantid::API::gsl_peak_wrapper;
  f.params = reinterpret_cast<void *>(&const_cast<IPeakFunction &>(peakFunction));
 
  return f;
}
 
double gsl_peak_wrapper(double x, void *parameters) {
  auto const *peakFunction = reinterpret_cast<IPeakFunction *>(parameters);
 
  if (!peakFunction) {
    throw std::runtime_error("Cannot process NULL-pointer in gsl_peak_wrapper.");
  }
 
  double y;
 
  /* For the integration to work properly, functionLocal has to be used instead
   * of the more general function-method due to the fact that the overriden
   * function-method
   * in IPeakFunction cuts off at some point. For slowly decaying peak functions
   * such as Lorentzians, this is introduces large deviations for integrations
   * from -Inf to +Inf.
   */
  peakFunction->functionLocal(&y, &x, 1);
 
  return y;
}
} // namespace Mantid::API