Initial commit.

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
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 (file)
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 (file)
index 0000000..939da69
--- /dev/null
+++ b/cubaint.cpp
@@ -0,0 +1,209 @@
+#include "cubaint.hpp"
+#include "cubaint_internal.hpp"
+
+#include <cmath>
+#include <iostream>
+#include <cassert>
+
+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<double> integral[], 
+                       std::complex<double> error[], std::complex<double> 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 <typename integrandtype, typename compreal>
+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<integrandtype> MetaUserData(userdata, Limits, integrand);
+  
+  switch (IntMethod) {
+  case SUAVE:    
+    Suave( Limits.size(), ncomp,
+          &masterintegrand<integrandtype,compreal>, &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<double*>(integral), 
+          reinterpret_cast<double*>(error), 
+          reinterpret_cast<double*>(prob) );
+    break;
+  case DIVONNE:
+    Divonne( Limits.size(), ncomp,
+            &masterintegrand<integrandtype,compreal>, 
+            &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<double*>(integral), 
+            reinterpret_cast<double*>(error), 
+            reinterpret_cast<double*>(prob) );
+    break;
+  case VEGAS:
+    Vegas( Limits.size(), ncomp,
+          &masterintegrand<integrandtype,compreal>, 
+          &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<double*>(integral), 
+          reinterpret_cast<double*>(error), 
+          reinterpret_cast<double*>(prob) );
+    break;
+  case CUHRE:
+    Cuhre( Limits.size(), ncomp,
+          &masterintegrand<integrandtype,compreal>, 
+          &MetaUserData, Options.nvec,
+          Options.epsrel, Options.epsabs,
+          Options.flags, Options.mineval, Options.maxeval,
+          Options.Cuhre.key,
+          Options.statefile,
+          &nregions, &neval, &fail,
+          reinterpret_cast<double*>(integral), 
+          reinterpret_cast<double*>(error), 
+          reinterpret_cast<double*>(prob) );
+    break;
+  }
+  
+  return 0;
+}
+
+const cubaint::options cubaint::DefaultOptions;
+
+void cubaint::kernelmap_intvar_numnum(const double& x01, double& xab,
+                                     const std::array<limit,2>& 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<limit,2>& range)
+{
+  xab = range[0].number + (1-x01)/x01;
+}
+
+void cubaint::kernelmap_intvar_infnum(const double& x01, double& xab,
+                                     const std::array<limit,2>& range)
+{
+  xab = range[1].number - (1-x01)/x01;
+}
+
+void cubaint::kernelmap_intvar_infinf(const double& x01, double& xab,
+                                     const std::array<limit,2>& range)
+{
+  xab = (2*x01-1) / ( 1 - pow( 2*x01-1, 2 ) );
+}
+
+double cubaint::kernelmap_intmulti_numnum(const double& x01,
+                                         const std::array<limit,2>& range)
+{
+  return range[1].number - range[0].number;
+}
+
+double cubaint::kernelmap_intmulti_numinf(const double& x01,
+                                         const std::array<limit,2>& range)
+{
+  return pow(x01,-2);
+}
+
+double cubaint::kernelmap_intmulti_infnum(const double& x01,
+                                         const std::array<limit,2>& range)
+{
+  return pow(x01,-2);
+}
+
+double cubaint::kernelmap_intmulti_infinf(const double& x01,
+                                         const std::array<limit,2>& 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 <typename integrandtype, typename compreal>
+int cubaint::masterintegrand(const int *ndim, const double x01[],
+                            const int *ncomp, double f[], void *_MetaUserData)
+{
+  metauserdata<integrandtype> *MetaUserData = (metauserdata<integrandtype>*)_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<compreal*>(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 (file)
index 0000000..103ea86
--- /dev/null
+++ b/cubaint.hpp
@@ -0,0 +1,100 @@
+#ifndef CUBAINT_HPP
+#define CUBAINT_HPP
+
+#include <vector>
+#include <array>
+#include <complex>
+
+#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<double> f[], void *userdata);
+
+  typedef std::vector<std::array<limit,2>> 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<double> integral[], 
+                std::complex<double> error[], std::complex<double> 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 (file)
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<limit,2>& range);
+
+  extern const kernelmap_intvar kernelmap_intvars[];
+
+  void kernelmap_intvar_numnum(const double& x01, double& xab,
+                              const std::array<limit,2>& range);
+  void kernelmap_intvar_numinf(const double& x01, double& xab,
+                              const std::array<limit,2>& range);
+  void kernelmap_intvar_infnum(const double& x01, double& xab,
+                              const std::array<limit,2>& range);
+  void kernelmap_intvar_infinf(const double& x01, double& xab,
+                              const std::array<limit,2>& range);
+
+  typedef double (*kernelmap_intmulti)(const double& x01,
+                                      const std::array<limit,2>& range);
+  
+  double kernelmap_intmulti_numnum(const double& x01,
+                                  const std::array<limit,2>& range);
+  double kernelmap_intmulti_numinf(const double& x01,
+                                  const std::array<limit,2>& range);
+  double kernelmap_intmulti_infnum(const double& x01,
+                                  const std::array<limit,2>& range);
+  double kernelmap_intmulti_infinf(const double& x01,
+                                  const std::array<limit,2>& range);
+
+  extern const kernelmap_intmulti kernelmap_intmultis[];
+
+  template <typename integrandtype, typename compreal>
+  int masterintegrand(const int *ndim, const double x01[],
+                     const int *ncomp, double f[], void *MetaUserData);
+
+  template <typename integrandtype>
+  struct metauserdata {
+    void *userdata;
+    const limits& Limits;
+    integrandtype integrand;
+
+    metauserdata (void *_userdata, const limits& _Limits, 
+                 integrandtype _integrand) :
+      userdata(_userdata), Limits(_Limits), integrand(_integrand)
+    {};
+  };
+
+  template <typename integrandtype, typename compreal>
+  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 (file)
index 0000000..93914bf
--- /dev/null
+++ b/demo.cpp
@@ -0,0 +1,106 @@
+#include <cmath>
+#include <iostream>
+
+#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<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 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<double> 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;
+  }
+
+}