GSL integrator added - LatAnalyze now depends on GSL

configure.ac

@@ -35,6 +35,14 @@ AC_ARG_ENABLE([SSE],
 				[Define to 1 if your CPU support SSE instructions.])],
 	[]
 )
+AC_ARG_WITH([gsl],
+    [AS_HELP_STRING([--with-gsl=prefix],
+		[try this for a non-standard install prefix of the GSL library])],
+    [AM_CFLAGS="$AM_CFLAGS -I$with_gsl/include"]
+	[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_gsl/include"]
+	[AM_LDFLAGS="$AM_LDFLAGS -L$with_gsl/lib"],
+	[]
+)
 AC_ARG_WITH([Minuit2],
     [AS_HELP_STRING([--with-Minuit2=prefix],
 		[try this for a non-standard install prefix of the Minuit2 library])],
@@ -70,6 +78,7 @@ AC_DEFINE_UNQUOTED([GXX_VERSION],["$GXX_VERSION"],
 
 # Checks for libraries.
 AC_CHECK_LIB([m],[cos],[],[AC_MSG_ERROR([libm library not found])])
+AC_CHECK_LIB([gsl],[gsl_blas_dgemm],[],[AC_MSG_ERROR([GSL library not found])])
 SAVED_LDFLAGS=$LDFLAGS
 LDFLAGS="$LDFLAGS -lMinuit2"
 AC_MSG_CHECKING([for ROOT::Minuit2::BasicMinimumError in -lMinuit2]);

examples/Makefile.am

@@ -17,6 +17,7 @@ endif
 noinst_PROGRAMS =           \
     exCompiledDoubleFunction\
     exFit                   \
+    exIntegrator            \
     exMat                   \
     exMathInterpreter       \
     exMin                   \
@@ -31,6 +32,10 @@ exFit_SOURCES                    = exFit.cpp
 exFit_CFLAGS                     = -g -O2
 exFit_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
+exIntegrator_SOURCES             = exIntegrator.cpp
+exIntegrator_CFLAGS              = -g -O2
+exIntegrator_LDFLAGS             = -L../lib/.libs -lLatAnalyze
+
 exMat_SOURCES                    = exMat.cpp
 exMat_CFLAGS                     = -g -O2
 exMat_LDFLAGS                    = -L../lib/.libs -lLatAnalyze

examples/exIntegrator.cpp

@@ -0,0 +1,32 @@
+#include <iostream>
+#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/GslQagsIntegrator.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char* argv[])
+{
+    string source;
+    double xMin, xMax;
+    
+    if (argc != 4)
+    {
+        cerr << "usage: " << argv[0] << " <function> <xMin> <xMax>" << endl;
+        
+        return EXIT_FAILURE;
+    }
+    source = argv[1];
+    xMin   = strTo<double>(argv[2]);
+    xMax   = strTo<double>(argv[3]);
+    
+    CompiledDoubleFunction f(1, source);
+    GslQagsIntegrator      integrator;
+    double                 result;
+    
+    result = integrator(f, xMin, xMax);
+    cout << "function integral on [" << xMin << ", " << xMax << "] = ";
+    cout << result << " +/- " << integrator.getLastError() <<endl;
+    
+    return EXIT_SUCCESS;
+}

lib/GslQagsIntegrator.cpp

@@ -0,0 +1,69 @@
+/*
+ * GslQagsIntegrator.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 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.
+ *
+ * LatAnalyze 3 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 LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/GslQagsIntegrator.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                    GslQagIntegrator implementation                         *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GslQagsIntegrator::GslQagsIntegrator(const unsigned int limit,
+                                     const double precision)
+: limit_(limit)
+, precision_(precision)
+{
+    workspace_ = gsl_integration_workspace_alloc(limit);
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+GslQagsIntegrator::~GslQagsIntegrator(void)
+{
+    gsl_integration_workspace_free(workspace_);
+}
+
+// integral calculation ////////////////////////////////////////////////////////
+double GslQagsIntegrator::operator()(const DoubleFunction &f, const double xMin,
+                                     const double xMax)
+{
+    double (*fWrap)(double, void *) = [](double x, void *fPt)
+    {
+        return (*static_cast<DoubleFunction *>(fPt))(&x);
+    };
+    
+    gsl_function gslF;
+    double       result;
+    
+    gslF.function = fWrap;
+    gslF.params   = reinterpret_cast<void *>(&const_cast<DoubleFunction &>(f));
+    
+    gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_, workspace_,
+                         &result, &error_);
+    
+    return result;
+}
+
+// get last error //////////////////////////////////////////////////////////////
+double GslQagsIntegrator::getLastError(void) const
+{
+    return error_;
+}

lib/GslQagsIntegrator.hpp

@@ -0,0 +1,58 @@
+/*
+ * GslQagsIntegrator.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 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.
+ *
+ * LatAnalyze 3 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 LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_GslQagsIntegrator_hpp_
+#define Latan_GslQagsIntegrator_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Integrator.hpp>
+#include <gsl/gsl_integration.h>
+
+BEGIN_NAMESPACE
+
+/******************************************************************************
+ *       GSL general quadrature adaptive integration with singularities       *
+ ******************************************************************************/
+
+class GslQagsIntegrator: public Integrator
+{
+public:
+    static const     unsigned int defaultLimit = 1000;
+    static constexpr double       defaultPrec  = 1.0e-7;
+public:
+    // constructor
+    GslQagsIntegrator(const unsigned int limit = defaultLimit,
+                      const double precision = defaultPrec);
+    // destructor
+    virtual ~GslQagsIntegrator(void);
+    // integral calculation
+    virtual double operator()(const DoubleFunction &f, const double xMin,
+                              const double xMax);
+    // get last error
+    double getLastError(void) const;
+private:
+    unsigned int              limit_;
+    double                    precision_, error_;
+    gsl_integration_workspace *workspace_;
+};
+
+END_NAMESPACE
+
+#endif // Latan_GslQagsIntegrator_hpp_

lib/Integrator.hpp

@@ -0,0 +1,46 @@
+/*
+ * Integrator.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 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.
+ *
+ * LatAnalyze 3 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 LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_Integrator_hpp_
+#define Latan_Integrator_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Function.hpp>
+
+BEGIN_NAMESPACE
+
+/******************************************************************************
+ *                        abstract integrator class                           *
+ ******************************************************************************/
+
+class Integrator
+{
+public:
+    // constructor
+    Integrator(void) = default;
+    // destructor
+    virtual ~Integrator(void) = default;
+    // integral calculation
+    virtual double operator()(const DoubleFunction &f, const double xMin,
+                              const double xMax) = 0;
+};
+
+END_NAMESPACE
+
+#endif // Latan_Integrator_hpp_

lib/Makefile.am

@@ -32,11 +32,12 @@ libLatAnalyze_la_SOURCES =\
     CompiledFunction.cpp  \
     CompiledModel.cpp     \
     Exceptions.cpp        \
-    Function.cpp          \
-    Global.cpp            \
-    includes.hpp          \
     File.cpp              \
     FitInterface.cpp      \
+    Function.cpp          \
+    Global.cpp            \
+    GslQagsIntegrator.cpp \
+    includes.hpp          \
     Mat.cpp               \
     Math.cpp              \
     MathInterpreter.cpp   \
@@ -62,8 +63,10 @@ libLatAnalyze_la_HEADERS =\
     Exceptions.hpp        \
     FitInterface.hpp      \
     Function.hpp          \
-    Global.hpp            \
     File.hpp              \
+    Global.hpp            \
+    GslQagsIntegrator.hpp \
+    Integrator.hpp        \
     IoObject.hpp          \
     Mat.hpp               \
     Math.hpp              \