From adc013e2895de042898e685c4cce29a387b4321d Mon Sep 17 00:00:00 2001 From: Alexander Schmidt Date: Tue, 26 Aug 2014 23:12:57 +0200 Subject: [PATCH 1/1] Initial commit. --- .gitignore | 7 ++ CMakeLists.txt | 10 +++ cubaint.cpp | 209 +++++++++++++++++++++++++++++++++++++++++++ cubaint.hpp | 100 +++++++++++++++++++++ cubaint_internal.hpp | 59 ++++++++++++ demo.cpp | 106 ++++++++++++++++++++++ 6 files changed, 491 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 cubaint.cpp create mode 100644 cubaint.hpp create mode 100644 cubaint_internal.hpp create mode 100644 demo.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..943b5a8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +CMakeFiles +CMakeCache.txt +*~ +Makefile +cmake_install.cmake +demo +libcubaint.a diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..deb7b3a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,10 @@ +cmake_minimum_required (VERSION 3.0) +project (CUBAINT) + +list (APPEND CMAKE_CXX_FLAGS "-std=c++11") + +add_library (cubaint cubaint.cpp) +target_link_libraries (cubaint cuba) + +add_executable (demo demo.cpp) +target_link_libraries (demo cubaint) diff --git a/cubaint.cpp b/cubaint.cpp new file mode 100644 index 0000000..939da69 --- /dev/null +++ b/cubaint.cpp @@ -0,0 +1,209 @@ +#include "cubaint.hpp" +#include "cubaint_internal.hpp" + +#include +#include +#include + +int cubaint::integrate (integrand_t integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, double integral[], + double error[], double prob[], + const unsigned int& verbosity, + void *userdata, const options& Options) +{ + assert (verbosity < 4); + + _integrate( integrand, IntMethod, Limits, ncomp, + nregions, neval, fail, integral, + error, prob, verbosity, userdata, Options ); +} + +int cubaint::integrate (integrand_tc integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, std::complex integral[], + std::complex error[], std::complex prob[], + const unsigned int& verbosity, + void *userdata, const options& Options) +{ + assert (verbosity < 4); + + _integrate( integrand, IntMethod, Limits, 2*ncomp, + nregions, neval, fail, integral, + error, prob, verbosity, userdata, Options ); +} + +template +int cubaint::_integrate (integrandtype integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, compreal integral[], + compreal error[], compreal prob[], + const unsigned int& verbosity, + void *userdata, const options& Options) +{ + metauserdata MetaUserData(userdata, Limits, integrand); + + switch (IntMethod) { + case SUAVE: + Suave( Limits.size(), ncomp, + &masterintegrand, &MetaUserData, Options.nvec, + Options.epsrel, Options.epsabs, + Options.flags | verbosity, Options.seed, + Options.mineval, Options.maxeval, + Options.Suave.nnew, Options.Suave.flatness, + Options.statefile, + &nregions, &neval, &fail, + reinterpret_cast(integral), + reinterpret_cast(error), + reinterpret_cast(prob) ); + break; + case DIVONNE: + Divonne( Limits.size(), ncomp, + &masterintegrand, + &MetaUserData, Options.nvec, + Options.epsrel, Options.epsabs, + Options.flags | verbosity, Options.seed, + Options.mineval, Options.maxeval, + Options.Divonne.key1, Options.Divonne.key2, + Options.Divonne.key3, Options.Divonne.maxpass, + Options.Divonne.border, Options.Divonne.maxchisq, + Options.Divonne.mindeviation, + Options.Divonne.ngiven, Options.Divonne.ldxgiven, + Options.Divonne.xgiven, + Options.Divonne.nextra, Options.Divonne.peakfinder, + Options.statefile, + &nregions, &neval, &fail, + reinterpret_cast(integral), + reinterpret_cast(error), + reinterpret_cast(prob) ); + break; + case VEGAS: + Vegas( Limits.size(), ncomp, + &masterintegrand, + &MetaUserData, Options.nvec, + Options.epsrel, Options.epsabs, + Options.flags | verbosity, Options.seed, + Options.mineval, Options.maxeval, + Options.Vegas.nstart, Options.Vegas.nincrease, Options.Vegas.nbatch, + Options.Vegas.gridno, Options.statefile, + &neval, &fail, + reinterpret_cast(integral), + reinterpret_cast(error), + reinterpret_cast(prob) ); + break; + case CUHRE: + Cuhre( Limits.size(), ncomp, + &masterintegrand, + &MetaUserData, Options.nvec, + Options.epsrel, Options.epsabs, + Options.flags, Options.mineval, Options.maxeval, + Options.Cuhre.key, + Options.statefile, + &nregions, &neval, &fail, + reinterpret_cast(integral), + reinterpret_cast(error), + reinterpret_cast(prob) ); + break; + } + + return 0; +} + +const cubaint::options cubaint::DefaultOptions; + +void cubaint::kernelmap_intvar_numnum(const double& x01, double& xab, + const std::array& range) +{ + xab = range[0].number + x01*( range[1].number - range[0].number ); +} + +void cubaint::kernelmap_intvar_numinf(const double& x01, double& xab, + const std::array& range) +{ + xab = range[0].number + (1-x01)/x01; +} + +void cubaint::kernelmap_intvar_infnum(const double& x01, double& xab, + const std::array& range) +{ + xab = range[1].number - (1-x01)/x01; +} + +void cubaint::kernelmap_intvar_infinf(const double& x01, double& xab, + const std::array& range) +{ + xab = (2*x01-1) / ( 1 - pow( 2*x01-1, 2 ) ); +} + +double cubaint::kernelmap_intmulti_numnum(const double& x01, + const std::array& range) +{ + return range[1].number - range[0].number; +} + +double cubaint::kernelmap_intmulti_numinf(const double& x01, + const std::array& range) +{ + return pow(x01,-2); +} + +double cubaint::kernelmap_intmulti_infnum(const double& x01, + const std::array& range) +{ + return pow(x01,-2); +} + +double cubaint::kernelmap_intmulti_infinf(const double& x01, + const std::array& range) +{ + const double dings = pow( 2*x01-1, 2 ); + return 2 * (1+dings) * pow( 1-dings, -2 ); +} + +const cubaint::kernelmap_intvar cubaint::kernelmap_intvars[] = + {&cubaint::kernelmap_intvar_numnum, &cubaint::kernelmap_intvar_numinf, + &cubaint::kernelmap_intvar_infnum, &cubaint::kernelmap_intvar_infinf}; + +const cubaint::kernelmap_intmulti cubaint::kernelmap_intmultis[] = + {&cubaint::kernelmap_intmulti_numnum, &cubaint::kernelmap_intmulti_numinf, + &cubaint::kernelmap_intmulti_infnum, &cubaint::kernelmap_intmulti_infinf}; + +cubaint::limit::limit(const double& _number) : + LimType(limtype::NUM), number(_number) +{ +} + +cubaint::limit::limit(const limtype& _LimType) : + LimType(_LimType), number(0) +{ +} + +template +int cubaint::masterintegrand(const int *ndim, const double x01[], + const int *ncomp, double f[], void *_MetaUserData) +{ + metauserdata *MetaUserData = (metauserdata*)_MetaUserData; + + double xab[*ndim]; + for (int idim=0; idim<*ndim; idim++) + kernelmap_intvars[ MetaUserData->Limits[idim][0].LimType*2 + + MetaUserData->Limits[idim][1].LimType ] + (x01[idim], xab[idim], MetaUserData->Limits[idim]); + + MetaUserData->integrand( ndim, xab, + ncomp, + reinterpret_cast(f), + MetaUserData->userdata ); + + for (int idim=0; idim<*ndim; idim++) + for (int icomp=0; icomp<*ncomp; icomp++) + f[icomp] *= + kernelmap_intmultis[ MetaUserData->Limits[idim][0].LimType*2 + + MetaUserData->Limits[idim][1].LimType ] + (x01[idim], MetaUserData->Limits[idim]); + + return 0; +} diff --git a/cubaint.hpp b/cubaint.hpp new file mode 100644 index 0000000..103ea86 --- /dev/null +++ b/cubaint.hpp @@ -0,0 +1,100 @@ +#ifndef CUBAINT_HPP +#define CUBAINT_HPP + +#include +#include +#include + +#include "cuba.h" + +namespace cubaint +{ + enum intmethod {SUAVE, DIVONNE, VEGAS, CUHRE}; + + enum limtype {NUM, INF}; + + struct limit { + limit(const double& _number); + limit(const limtype& _LimType); + const limtype LimType; + const double number; + }; + + typedef int (*integrand_t)(const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata); + + typedef int (*integrand_tc)(const int *ndim, const double x[], + const int *ncomp, std::complex f[], void *userdata); + + typedef std::vector> limits; + + struct options { + int seed, flags, nvec, mineval, maxeval; + double epsrel, epsabs; + const char *statefile; + int *spin; + options() : + seed(0), flags(0), nvec(1), mineval(0), maxeval(50000), + epsrel(1e-3), epsabs(1e-12), statefile(NULL), spin(NULL) + {}; + + struct suave { + int nnew; + double flatness; + suave() : + nnew(1000), flatness(25.) + {}; + }; + + struct divonne { + int key1, key2, key3, maxpass, ngiven, ldxgiven, nextra; + double border, maxchisq, mindeviation; + double *xgiven; + peakfinder_t peakfinder; + divonne() : + key1(47), key2(1), key3(1), maxpass(5), border(0), maxchisq(10), + mindeviation(.25), ngiven(0), ldxgiven(0), xgiven(NULL), + peakfinder(NULL), nextra(0) + {}; + }; + + struct vegas { + int nstart, nincrease, nbatch, gridno; + vegas() : + nstart(1000), nincrease(500), nbatch(1000), gridno(0) + {}; + }; + + struct cuhre { + int key; + cuhre() : + key(0) + {}; + }; + + suave Suave; + divonne Divonne; + vegas Vegas; + cuhre Cuhre; + }; + + extern const options DefaultOptions; + + int integrate (integrand_t integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, double integral[], + double error[], double prob[], + const unsigned int& verbosity=0, + void *userdata=NULL, const options& Options=DefaultOptions); + + int integrate (integrand_tc integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, std::complex integral[], + std::complex error[], std::complex prob[], + const unsigned int& verbosity=0, + void *userdata=NULL, const options& Options=DefaultOptions); +}; + +#endif diff --git a/cubaint_internal.hpp b/cubaint_internal.hpp new file mode 100644 index 0000000..e5ca94e --- /dev/null +++ b/cubaint_internal.hpp @@ -0,0 +1,59 @@ +#ifndef CUBAINT_INTERNAL_HPP +#define CUBAINT_INTERNAL_HPP + +namespace cubaint { + typedef void (*kernelmap_intvar)(const double& x01, double& xab, + const std::array& range); + + extern const kernelmap_intvar kernelmap_intvars[]; + + void kernelmap_intvar_numnum(const double& x01, double& xab, + const std::array& range); + void kernelmap_intvar_numinf(const double& x01, double& xab, + const std::array& range); + void kernelmap_intvar_infnum(const double& x01, double& xab, + const std::array& range); + void kernelmap_intvar_infinf(const double& x01, double& xab, + const std::array& range); + + typedef double (*kernelmap_intmulti)(const double& x01, + const std::array& range); + + double kernelmap_intmulti_numnum(const double& x01, + const std::array& range); + double kernelmap_intmulti_numinf(const double& x01, + const std::array& range); + double kernelmap_intmulti_infnum(const double& x01, + const std::array& range); + double kernelmap_intmulti_infinf(const double& x01, + const std::array& range); + + extern const kernelmap_intmulti kernelmap_intmultis[]; + + template + int masterintegrand(const int *ndim, const double x01[], + const int *ncomp, double f[], void *MetaUserData); + + template + struct metauserdata { + void *userdata; + const limits& Limits; + integrandtype integrand; + + metauserdata (void *_userdata, const limits& _Limits, + integrandtype _integrand) : + userdata(_userdata), Limits(_Limits), integrand(_integrand) + {}; + }; + + template + int _integrate (integrandtype integrand, const intmethod& IntMethod, + const limits& Limits, const int& ncomp, + int& nregions, + int& neval, int& fail, compreal integral[], + compreal error[], compreal prob[], + const unsigned int& verbosity=0, void *userdata=NULL, + const options& Options=DefaultOptions); +}; + +#endif diff --git a/demo.cpp b/demo.cpp new file mode 100644 index 0000000..93914bf --- /dev/null +++ b/demo.cpp @@ -0,0 +1,106 @@ +#include +#include + +#include "cubaint.hpp" + +int test_integrand (const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata) +{ + for (int icomp=0; icomp<*ncomp; icomp++) { + f[icomp] = 1; + for (int idim=0; idim<*ndim; idim++) + f[icomp] *= exp(-pow(x[idim],2)); + } + + return 0; +} + +int test_integrand_c (const int *ndim, const double x[], + const int *ncomp, std::complex f[], void *userdata) +{ + for (int icomp=0; icomp<*ncomp; icomp++) { + f[icomp] = 1.; + for (int idim=0; idim<*ndim; idim++) + f[icomp] *= exp(-pow(x[idim],2)); + } + + return 0; +} + +int main () +{ + int nregions, neval, fail; + double error, prob, integral; + + std::cout << "int_0^1 dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{0,1}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::cout << "int_0^INF dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{0,cubaint::INF}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::cout << "int_-INF^0 dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{cubaint::INF,0}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::cout << "int_-INF^INF dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{cubaint::INF,cubaint::INF}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::cout << "int_1^INF dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{1,cubaint::INF}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::cout << "int_-INF^1 dx e-x^2: " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand, cubaint::intmethod(itype), + {{cubaint::INF,1}}, 1, + nregions, + neval, fail, &integral, + &error, &prob, 0); + std::cout << integral << " +- " << error << " (" << itype << ")" << std::endl; + } + + std::complex error_c, prob_c, integral_c; + std::cout << "int_0^1 dx (e-x^2+i0): " << std::endl; + for (int itype=0; itype<4; itype++) { + cubaint::integrate (&test_integrand_c, cubaint::intmethod(itype), + {{0,1}}, 2, + nregions, + neval, fail, &integral_c, + &error_c, &prob_c, 0); + std::cout << integral_c << " +- " << error_c << " (" << itype << ")" << std::endl; + } + +} -- 2.39.5