cubaint.cpp

   1 #include "cubaint.hpp"
   2 #include "cubaint_internal.hpp"
   3
   4 #include <cmath>
   5 #include <iostream>
   6 #include <cassert>
   7
   8 int cubaint::integrate (integrand_t integrand, const intmethod& IntMethod,
   9                         const limits& Limits, const int& ncomp,
  10                         int& nregions,
  11                         int& neval, int& fail, double integral[],
  12                         double error[], double prob[],
  13                         const unsigned int& verbosity,
  14                         void *userdata, const options& Options)
  15 {
  16   assert (verbosity < 4);
  17
  18   return _integrate( integrand, IntMethod, Limits, ncomp,
  19                      nregions, neval, fail, integral,
  20                      error, prob, verbosity, userdata, Options );
  21 }
  22
  23 int cubaint::integrate (integrand_tc integrand, const intmethod& IntMethod,
  24                         const limits& Limits, const int& ncomp,
  25                         int& nregions,
  26                         int& neval, int& fail, std::complex<double> integral[],
  27                         std::complex<double> error[], std::complex<double> prob[],
  28                         const unsigned int& verbosity,
  29                         void *userdata, const options& Options)
  30 {
  31   assert (verbosity < 4);
  32
  33   return _integrate( integrand, IntMethod, Limits, 2*ncomp,
  34                      nregions, neval, fail, integral,
  35                      error, prob, verbosity, userdata, Options );
  36 }
  37
  38 template <typename integrandtype, typename compreal>
  39 int cubaint::_integrate (integrandtype integrand, const intmethod& IntMethod,
  40                          const limits& Limits, const int& ncomp,
  41                          int& nregions,
  42                          int& neval, int& fail, compreal integral[],
  43                          compreal error[], compreal prob[],
  44                          const unsigned int& verbosity,
  45                          void *userdata, const options& Options)
  46 {
  47   metauserdata<integrandtype> MetaUserData(userdata, Limits, integrand);
  48
  49   switch (IntMethod) {
  50   case SUAVE:
  51     Suave( Limits.size(), ncomp,
  52            &masterintegrand<integrandtype,compreal>, &MetaUserData, Options.nvec,
  53            Options.epsrel, Options.epsabs,
  54            Options.flags | verbosity, Options.seed,
  55            Options.mineval, Options.maxeval,
  56            Options.Suave.nnew, Options.Suave.nmin,
  57            Options.Suave.flatness,
  58            Options.statefile, Options.spin,
  59            &nregions, &neval, &fail,
  60            reinterpret_cast<double*>(integral),
  61            reinterpret_cast<double*>(error),
  62            reinterpret_cast<double*>(prob) );
  63     break;
  64   case DIVONNE:
  65     Divonne( Limits.size(), ncomp,
  66              &masterintegrand<integrandtype,compreal>,
  67              &MetaUserData, Options.nvec,
  68              Options.epsrel, Options.epsabs,
  69              Options.flags | verbosity, Options.seed,
  70              Options.mineval, Options.maxeval,
  71              Options.Divonne.key1, Options.Divonne.key2,
  72              Options.Divonne.key3, Options.Divonne.maxpass,
  73              Options.Divonne.border, Options.Divonne.maxchisq,
  74              Options.Divonne.mindeviation,
  75              Options.Divonne.ngiven, Options.Divonne.ldxgiven,
  76              Options.Divonne.xgiven,
  77              Options.Divonne.nextra, Options.Divonne.peakfinder,
  78              Options.statefile, Options.spin,
  79              &nregions, &neval, &fail,
  80              reinterpret_cast<double*>(integral),
  81              reinterpret_cast<double*>(error),
  82              reinterpret_cast<double*>(prob) );
  83     break;
  84   case VEGAS:
  85     Vegas( Limits.size(), ncomp,
  86            &masterintegrand<integrandtype,compreal>,
  87            &MetaUserData, Options.nvec,
  88            Options.epsrel, Options.epsabs,
  89            Options.flags | verbosity, Options.seed,
  90            Options.mineval, Options.maxeval,
  91            Options.Vegas.nstart, Options.Vegas.nincrease, Options.Vegas.nbatch,
  92            Options.Vegas.gridno, Options.statefile, Options.spin,
  93            &neval, &fail,
  94            reinterpret_cast<double*>(integral),
  95            reinterpret_cast<double*>(error),
  96            reinterpret_cast<double*>(prob) );
  97     break;
  98   case CUHRE:
  99     Cuhre( Limits.size(), ncomp,
 100            &masterintegrand<integrandtype,compreal>,
 101            &MetaUserData, Options.nvec,
 102            Options.epsrel, Options.epsabs,
 103            Options.flags, Options.mineval, Options.maxeval,
 104            Options.Cuhre.key,
 105            Options.statefile, Options.spin,
 106            &nregions, &neval, &fail,
 107            reinterpret_cast<double*>(integral),
 108            reinterpret_cast<double*>(error),
 109            reinterpret_cast<double*>(prob) );
 110     break;
 111   }
 112
 113   return 0;
 114 }
 115
 116 const cubaint::options cubaint::DefaultOptions;
 117
 118 void cubaint::kernelmap_intvar_numnum(const double& x01, double& xab,
 119                                       const std::array<limit,2>& range)
 120 {
 121   xab = range[0].number + x01*( range[1].number - range[0].number );
 122 }
 123
 124 void cubaint::kernelmap_intvar_numinf(const double& x01, double& xab,
 125                                       const std::array<limit,2>& range)
 126 {
 127   xab = range[0].number + (1-x01)/x01;
 128 }
 129
 130 void cubaint::kernelmap_intvar_infnum(const double& x01, double& xab,
 131                                       const std::array<limit,2>& range)
 132 {
 133   xab = range[1].number - (1-x01)/x01;
 134 }
 135
 136 void cubaint::kernelmap_intvar_infinf(const double& x01, double& xab,
 137                                       const std::array<limit,2>& range)
 138 {
 139   xab = (2*x01-1) / ( 1 - pow( 2*x01-1, 2 ) );
 140 }
 141
 142 double cubaint::kernelmap_intmulti_numnum(const double& x01,
 143                                           const std::array<limit,2>& range)
 144 {
 145   return range[1].number - range[0].number;
 146 }
 147
 148 double cubaint::kernelmap_intmulti_numinf(const double& x01,
 149                                           const std::array<limit,2>& range)
 150 {
 151   return pow(x01,-2);
 152 }
 153
 154 double cubaint::kernelmap_intmulti_infnum(const double& x01,
 155                                           const std::array<limit,2>& range)
 156 {
 157   return pow(x01,-2);
 158 }
 159
 160 double cubaint::kernelmap_intmulti_infinf(const double& x01,
 161                                           const std::array<limit,2>& range)
 162 {
 163   const double dings = pow( 2*x01-1, 2 );
 164   return 2 * (1+dings) * pow( 1-dings, -2 );
 165 }
 166
 167 const cubaint::kernelmap_intvar cubaint::kernelmap_intvars[] =
 168   {&cubaint::kernelmap_intvar_numnum, &cubaint::kernelmap_intvar_numinf,
 169    &cubaint::kernelmap_intvar_infnum, &cubaint::kernelmap_intvar_infinf};
 170
 171 const cubaint::kernelmap_intmulti cubaint::kernelmap_intmultis[] =
 172   {&cubaint::kernelmap_intmulti_numnum, &cubaint::kernelmap_intmulti_numinf,
 173    &cubaint::kernelmap_intmulti_infnum, &cubaint::kernelmap_intmulti_infinf};
 174
 175 cubaint::limit::limit(const double& _number) :
 176   LimType(limtype::NUM), number(_number)
 177 {
 178 }
 179
 180 cubaint::limit::limit(const limtype& _LimType) :
 181   LimType(_LimType), number(0)
 182 {
 183 }
 184
 185 template <typename integrandtype, typename compreal>
 186 int cubaint::masterintegrand(const int *ndim, const double x01[],
 187                              const int *ncomp, double f[], void *_MetaUserData)
 188 {
 189   metauserdata<integrandtype> *MetaUserData = (metauserdata<integrandtype>*)_MetaUserData;
 190
 191   double xab[*ndim];
 192   for (int idim=0; idim<*ndim; idim++)
 193     kernelmap_intvars[ MetaUserData->Limits[idim][0].LimType*2 +
 194                        MetaUserData->Limits[idim][1].LimType ]
 195       (x01[idim], xab[idim], MetaUserData->Limits[idim]);
 196
 197   MetaUserData->integrand( ndim, xab,
 198                            ncomp,
 199                            reinterpret_cast<compreal*>(f),
 200                            MetaUserData->userdata );
 201
 202   for (int idim=0; idim<*ndim; idim++)
 203     for (int icomp=0; icomp<*ncomp; icomp++)
 204       f[icomp] *=
 205         kernelmap_intmultis[ MetaUserData->Limits[idim][0].LimType*2 +
 206                              MetaUserData->Limits[idim][1].LimType ]
 207         (x01[idim], MetaUserData->Limits[idim]);
 208
 209   return 0;
 210 }