#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);

  return _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);

  return _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.nmin,
	   Options.Suave.flatness,
	   Options.statefile, Options.spin,
	   &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, Options.spin,
	     &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, Options.spin,
	   &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, Options.spin,
	   &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;
}