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   _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   _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.flatness,
  57            Options.statefile,
  58            &nregions, &neval, &fail,
  59            reinterpret_cast<double*>(integral),
  60            reinterpret_cast<double*>(error),
  61            reinterpret_cast<double*>(prob) );
  62     break;
  63   case DIVONNE:
  64     Divonne( Limits.size(), ncomp,
  65              &masterintegrand<integrandtype,compreal>,
  66              &MetaUserData, Options.nvec,
  67              Options.epsrel, Options.epsabs,
  68              Options.flags | verbosity, Options.seed,
  69              Options.mineval, Options.maxeval,
  70              Options.Divonne.key1, Options.Divonne.key2,
  71              Options.Divonne.key3, Options.Divonne.maxpass,
  72              Options.Divonne.border, Options.Divonne.maxchisq,
  73              Options.Divonne.mindeviation,
  74              Options.Divonne.ngiven, Options.Divonne.ldxgiven,
  75              Options.Divonne.xgiven,
  76              Options.Divonne.nextra, Options.Divonne.peakfinder,
  77              Options.statefile,
  78              &nregions, &neval, &fail,
  79              reinterpret_cast<double*>(integral),
  80              reinterpret_cast<double*>(error),
  81              reinterpret_cast<double*>(prob) );
  82     break;
  83   case VEGAS:
  84     Vegas( Limits.size(), ncomp,
  85            &masterintegrand<integrandtype,compreal>,
  86            &MetaUserData, Options.nvec,
  87            Options.epsrel, Options.epsabs,
  88            Options.flags | verbosity, Options.seed,
  89            Options.mineval, Options.maxeval,
  90            Options.Vegas.nstart, Options.Vegas.nincrease, Options.Vegas.nbatch,
  91            Options.Vegas.gridno, Options.statefile,
  92            &neval, &fail,
  93            reinterpret_cast<double*>(integral),
  94            reinterpret_cast<double*>(error),
  95            reinterpret_cast<double*>(prob) );
  96     break;
  97   case CUHRE:
  98     Cuhre( Limits.size(), ncomp,
  99            &masterintegrand<integrandtype,compreal>,
 100            &MetaUserData, Options.nvec,
 101            Options.epsrel, Options.epsabs,
 102            Options.flags, Options.mineval, Options.maxeval,
 103            Options.Cuhre.key,
 104            Options.statefile,
 105            &nregions, &neval, &fail,
 106            reinterpret_cast<double*>(integral),
 107            reinterpret_cast<double*>(error),
 108            reinterpret_cast<double*>(prob) );
 109     break;
 110   }
 111
 112   return 0;
 113 }
 114
 115 const cubaint::options cubaint::DefaultOptions;
 116
 117 void cubaint::kernelmap_intvar_numnum(const double& x01, double& xab,
 118                                       const std::array<limit,2>& range)
 119 {
 120   xab = range[0].number + x01*( range[1].number - range[0].number );
 121 }
 122
 123 void cubaint::kernelmap_intvar_numinf(const double& x01, double& xab,
 124                                       const std::array<limit,2>& range)
 125 {
 126   xab = range[0].number + (1-x01)/x01;
 127 }
 128
 129 void cubaint::kernelmap_intvar_infnum(const double& x01, double& xab,
 130                                       const std::array<limit,2>& range)
 131 {
 132   xab = range[1].number - (1-x01)/x01;
 133 }
 134
 135 void cubaint::kernelmap_intvar_infinf(const double& x01, double& xab,
 136                                       const std::array<limit,2>& range)
 137 {
 138   xab = (2*x01-1) / ( 1 - pow( 2*x01-1, 2 ) );
 139 }
 140
 141 double cubaint::kernelmap_intmulti_numnum(const double& x01,
 142                                           const std::array<limit,2>& range)
 143 {
 144   return range[1].number - range[0].number;
 145 }
 146
 147 double cubaint::kernelmap_intmulti_numinf(const double& x01,
 148                                           const std::array<limit,2>& range)
 149 {
 150   return pow(x01,-2);
 151 }
 152
 153 double cubaint::kernelmap_intmulti_infnum(const double& x01,
 154                                           const std::array<limit,2>& range)
 155 {
 156   return pow(x01,-2);
 157 }
 158
 159 double cubaint::kernelmap_intmulti_infinf(const double& x01,
 160                                           const std::array<limit,2>& range)
 161 {
 162   const double dings = pow( 2*x01-1, 2 );
 163   return 2 * (1+dings) * pow( 1-dings, -2 );
 164 }
 165
 166 const cubaint::kernelmap_intvar cubaint::kernelmap_intvars[] =
 167   {&cubaint::kernelmap_intvar_numnum, &cubaint::kernelmap_intvar_numinf,
 168    &cubaint::kernelmap_intvar_infnum, &cubaint::kernelmap_intvar_infinf};
 169
 170 const cubaint::kernelmap_intmulti cubaint::kernelmap_intmultis[] =
 171   {&cubaint::kernelmap_intmulti_numnum, &cubaint::kernelmap_intmulti_numinf,
 172    &cubaint::kernelmap_intmulti_infnum, &cubaint::kernelmap_intmulti_infinf};
 173
 174 cubaint::limit::limit(const double& _number) :
 175   LimType(limtype::NUM), number(_number)
 176 {
 177 }
 178
 179 cubaint::limit::limit(const limtype& _LimType) :
 180   LimType(_LimType), number(0)
 181 {
 182 }
 183
 184 template <typename integrandtype, typename compreal>
 185 int cubaint::masterintegrand(const int *ndim, const double x01[],
 186                              const int *ncomp, double f[], void *_MetaUserData)
 187 {
 188   metauserdata<integrandtype> *MetaUserData = (metauserdata<integrandtype>*)_MetaUserData;
 189
 190   double xab[*ndim];
 191   for (int idim=0; idim<*ndim; idim++)
 192     kernelmap_intvars[ MetaUserData->Limits[idim][0].LimType*2 +
 193                        MetaUserData->Limits[idim][1].LimType ]
 194       (x01[idim], xab[idim], MetaUserData->Limits[idim]);
 195
 196   MetaUserData->integrand( ndim, xab,
 197                            ncomp,
 198                            reinterpret_cast<compreal*>(f),
 199                            MetaUserData->userdata );
 200
 201   for (int idim=0; idim<*ndim; idim++)
 202     for (int icomp=0; icomp<*ncomp; icomp++)
 203       f[icomp] *=
 204         kernelmap_intmultis[ MetaUserData->Limits[idim][0].LimType*2 +
 205                              MetaUserData->Limits[idim][1].LimType ]
 206         (x01[idim], MetaUserData->Limits[idim]);
 207
 208   return 0;
 209 }