Refs #9445. Added PeakFunctionIntegrator

PeakFunctionIntegrator takes an IPeakFunction and integrates it numerically using GSL routines.

Code/Mantid/Framework/SINQ/CMakeLists.txt

@@ -13,6 +13,7 @@ set ( SRC_FILES
 	src/PoldiPeakSearch.cpp
 	src/PoldiRemoveDeadWires.cpp
 	src/PoldiUtilities/MillerIndices.cpp
+        src/PoldiUtilities/PeakFunctionIntegrator.cpp
 	src/PoldiUtilities/PoldiAutoCorrelationCore.cpp
 	src/PoldiUtilities/PoldiBasicChopper.cpp
 	src/PoldiUtilities/PoldiChopperFactory.cpp
@@ -52,6 +53,7 @@ set ( INC_FILES
 	inc/MantidSINQ/PoldiRemoveDeadWires.h
 	inc/MantidSINQ/PoldiUtilities/MillerIndices.h
 	inc/MantidSINQ/PoldiUtilities/MillerIndicesIO.h
+        inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h
 	inc/MantidSINQ/PoldiUtilities/PoldiAbstractChopper.h
 	inc/MantidSINQ/PoldiUtilities/PoldiAbstractDetector.h
 	inc/MantidSINQ/PoldiUtilities/PoldiAutoCorrelationCore.h
@@ -83,6 +85,7 @@ set ( TEST_FILES
 	MDHistoToWorkspace2DTest.h
 	MillerIndicesIOTest.h
 	MillerIndicesTest.h
+        PeakFunctionIntegratorTest.h
 	PoldiAutoCorrelationCoreTest.h
 	PoldiBasicChopperTest.h
 	PoldiChopperFactoryTest.h

Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h

@@ -0,0 +1,79 @@
+#ifndef PEAKFUNCTIONINTEGRATOR_H
+#define PEAKFUNCTIONINTEGRATOR_H
+
+#include "MantidSINQ/DllConfig.h"
+#include "MantidAPI/IPeakFunction.h"
+#include "gsl/gsl_integration.h"
+
+namespace Mantid {
+namespace Poldi {
+
+/** PeakFunctionIntegrator :
+ *
+  General integration of peaks (in the form of IPeakFunction) by wrapping the
+  corresponding GSL-functions. Integration with infinity limits is supported.
+
+    @author Michael Wedel, Paul Scherrer Institut - SINQ
+    @date 24/04/2014
+
+    Copyright © 2014 PSI-MSS
+
+  This file is part of Mantid.
+
+  Mantid is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 3 of the License, or
+  (at your option) any later version.
+
+  Mantid is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+  File change history is stored at: <https://github.com/mantidproject/mantid>
+  Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+
+using namespace Mantid::API;
+
+struct MANTID_SINQ_DLL IntegrationResult {
+    double result;
+    double error;
+    size_t intervals;
+
+    int errorCode;
+    bool success;
+};
+
+class MANTID_SINQ_DLL PeakFunctionIntegrator {
+public:
+    PeakFunctionIntegrator(double requiredRelativePrecision = 1e-8);
+    virtual ~PeakFunctionIntegrator();
+
+    void setRequiredRelativePrecision(double newPrecision);
+    double requiredRelativePrecision() const;
+
+    IntegrationResult integrateInfinity(IPeakFunction_const_sptr peakFunction) const;
+    IntegrationResult integratePositiveInfinity(IPeakFunction_const_sptr peakFunction, double lowerLimit) const;
+    IntegrationResult integrateNegativeInfinity(IPeakFunction_const_sptr peakFunction, double upperLimit) const;
+
+    IntegrationResult integrate(IPeakFunction_const_sptr peakFunction, double lowerLimit, double upperLimit) const;
+
+protected:
+    gsl_function getGSLFunction(IPeakFunction_const_sptr peakFunction) const;
+    void throwIfInvalid(IPeakFunction_const_sptr peakFunction) const;
+
+    gsl_integration_workspace *m_integrationWorkspace;
+
+    double m_relativePrecision;
+};
+
+double gsl_peak_wrapper(double x, void *parameters);
+
+}
+}
+
+#endif // PEAKFUNCTIONINTEGRATOR_H

Code/Mantid/Framework/SINQ/src/PoldiUtilities/PeakFunctionIntegrator.cpp

@@ -0,0 +1,173 @@
+#include "MantidSINQ/PoldiUtilities/PeakFunctionIntegrator.h"
+
+#include "MantidAPI/FunctionDomain1D.h"
+#include "gsl/gsl_errno.h"
+#include <iostream>
+#include <iomanip>
+
+namespace Mantid {
+namespace Poldi {
+
+/** Constructor with required relative precision argument. The default is 1e-8.
+ *  See also PeakFunctionIntegrator::setRequiredRelativePrecision.
+ *
+ * @param requiredRelativePrecision :: Desired relative precision of the integral estimations.
+ */
+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);
+}
+
+/** This method sets the desired numerical relative precision that's passed on to the
+ *  GSL integration-routines.
+ *
+ *  @param newPrecision :: Desired relative precision for integrations.
+ */
+void PeakFunctionIntegrator::setRequiredRelativePrecision(double newPrecision)
+{
+    m_relativePrecision = newPrecision;
+}
+
+/** Returns the currently set precision
+ */
+double PeakFunctionIntegrator::requiredRelativePrecision() const
+{
+    return m_relativePrecision;
+}
+
+/** Integration of peak function on the interval [-Inf, +Inf]. Internally, gsl_integration_qagi is used
+ *  for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *  The results are returned as IntegrationResult-struct, which contains the approximation of the integral along
+ *  with other information such as an error estimate (absolute).
+ *
+ *  @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrateInfinity(IPeakFunction_const_sptr peakFunction) const
+{
+    throwIfInvalid(peakFunction);
+
+    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;
+}
+
+/** Integration of peak function on the interval [a, +Inf]. Internally, gsl_integration_qagiu is used
+ *  for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ *  @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integratePositiveInfinity(IPeakFunction_const_sptr peakFunction, double lowerLimit) const
+{
+    throwIfInvalid(peakFunction);
+
+    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;
+}
+
+/** Integration of peak function on the interval [-Inf, b]. Internally, gsl_integration_qagil is used
+ *  for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ *  @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrateNegativeInfinity(IPeakFunction_const_sptr peakFunction, double upperLimit) const
+{
+    throwIfInvalid(peakFunction);
+
+    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;
+}
+
+/** Integration of peak function on the interval [a, b]. Internally, gsl_integration_qags is used
+ *  for this. If a default constructed IPeakFunction_const_sptr is passed to the function, std::invalid_argument is thrown.
+ *
+ *  @param peakFunction :: Peak function to integrate.
+ */
+IntegrationResult PeakFunctionIntegrator::integrate(IPeakFunction_const_sptr peakFunction, double lowerLimit, double upperLimit) const
+{
+    throwIfInvalid(peakFunction);
+
+    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;
+}
+
+/** Method that wraps an IPeakFunction for use with GSL functions.
+ *
+ *  @param peakFunction :: Peak function to wrap.
+ */
+gsl_function PeakFunctionIntegrator::getGSLFunction(IPeakFunction_const_sptr peakFunction) const
+{
+    gsl_function f;
+    f.function = &Mantid::Poldi::gsl_peak_wrapper;
+    f.params = &peakFunction;
+
+    return f;
+}
+
+void PeakFunctionIntegrator::throwIfInvalid(IPeakFunction_const_sptr peakFunction) const
+{
+    if(!peakFunction) {
+        throw std::invalid_argument("Can not integrate NULL-function.");
+    }
+}
+
+double gsl_peak_wrapper(double x, void *parameters)
+{
+    IPeakFunction_const_sptr peakFunction = *(IPeakFunction_const_sptr *) parameters;
+
+    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;
+}
+
+
+
+}
+
+
+
+}