Integrator_VEGAS.h

#ifndef QNANO_NEW_INTEGRATOR_VEGAS_DEFINED_H
#define QNANO_NEW_INTEGRATOR_VEGAS_DEFINED_H

#include "tools/Integrator.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>


class Integrator_VEGAS: public Integrator{
public:

  virtual double integrate(const Integrable_Function& ifct, const Integrator_Discretization& discr) const{
    if(discr.size()!=ifct.dim){
      std::cerr<<"Integrator_VEGAS::integrate: discr.size()!=ifct.dim"<<std::endl;
      exit(1);
    }
    gsl_monte_vegas_state *vst=gsl_monte_vegas_alloc(ifct.dim);
    double res,err;

    gsl_rng_env_setup();
    const gsl_rng_type *T=gsl_rng_default;
    gsl_rng *r=gsl_rng_alloc(T);

    gsl_monte_function G = { ifct.fct, ifct.dim, (void *)ifct.params };

    size_t Ntot=discr[0].N;

    double *lmin=new double[discr.size()];
    for(size_t i=0; i<discr.size(); i++)lmin[i]=discr[i].lim_a;
    double *lmax=new double[discr.size()];
    for(size_t i=0; i<discr.size(); i++)lmax[i]=discr[i].lim_b;

    gsl_monte_vegas_integrate(&G, lmin, lmax, ifct.dim, Ntot, r, vst, &res, &err);
    
    delete [] lmax;
    delete [] lmin;

    gsl_monte_vegas_free(vst);
    return res;
  }

};

#endif