diff --git a/configure.ac b/configure.ac index b02e023..7cb157f 100644 --- a/configure.ac +++ b/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]); diff --git a/examples/Makefile.am b/examples/Makefile.am index 545aa32..ee27c4c 100644 --- a/examples/Makefile.am +++ b/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 diff --git a/examples/exIntegrator.cpp b/examples/exIntegrator.cpp new file mode 100644 index 0000000..730d93a --- /dev/null +++ b/examples/exIntegrator.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +using namespace std; +using namespace Latan; + +int main(int argc, char* argv[]) +{ + string source; + double xMin, xMax; + + if (argc != 4) + { + cerr << "usage: " << argv[0] << " " << endl; + + return EXIT_FAILURE; + } + source = argv[1]; + xMin = strTo(argv[2]); + xMax = strTo(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() <. + */ + +#include +#include + +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(fPt))(&x); + }; + + gsl_function gslF; + double result; + + gslF.function = fWrap; + gslF.params = reinterpret_cast(&const_cast(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_; +} diff --git a/lib/GslQagsIntegrator.hpp b/lib/GslQagsIntegrator.hpp new file mode 100644 index 0000000..2559bb4 --- /dev/null +++ b/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 . + */ + +#ifndef Latan_GslQagsIntegrator_hpp_ +#define Latan_GslQagsIntegrator_hpp_ + +#include +#include +#include +#include + +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_ diff --git a/lib/Integrator.hpp b/lib/Integrator.hpp new file mode 100644 index 0000000..1e2cfe9 --- /dev/null +++ b/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 . + */ + +#ifndef Latan_Integrator_hpp_ +#define Latan_Integrator_hpp_ + +#include +#include + +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_ diff --git a/lib/Makefile.am b/lib/Makefile.am index b54679f..4e60d79 100644 --- a/lib/Makefile.am +++ b/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 \