#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"

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

#define NEED_sv_2pv_flags
#include "ppport.h"

struct params {
  SV* eqn;
};

double integrand(double *x, size_t size, void *params) {

  dSP;

  SV* eqn;
  int i, count, dim;
  double val;

  eqn = ((struct params *)params)->eqn;

  ENTER;
  SAVETMPS;

  PUSHMARK(SP);

  dim = (int)size;
  for(i=0; i<dim; i++) {
    XPUSHs(sv_2mortal(newSVnv(x[i])));
  }

  PUTBACK;

  count = call_sv(eqn, G_SCALAR);
  if (count != 1) 
    croak("Integrand closure did not return a value");

  SPAGAIN;

  val = POPn;

  PUTBACK;
  FREETMPS;
  LEAVE;

  return val;
}

AV* c_int_multi(SV* eqn, AV* lower, AV* upper, int calls) {

  int i, dim;
  double * xl;
  double * xu;
  double res, err;
  size_t dim_size, calls_size;
  const gsl_rng_type *T;
  gsl_rng *r;
  struct params myparams;
  gsl_monte_function F;
  AV* ret = newAV();

  /* deal with lower and upper limits */

  dim = av_len(lower) + 1;
  if((av_len(upper)+1) != dim)
    croak("dimension mismatch");

  Newx(xl, dim, double);
  if(xl == NULL) 
    croak ("Failed to allocate memory to 'xl' in 'c_int_multi'");
  Newx(xu, dim, double);
  if(xu == NULL) 
    croak ("Failed to allocate memory to 'xu' in 'c_int_multi'");

  for (i=0; i<dim; i++) {
    xl[i] = SvNV(*av_fetch(lower, i, 0));
    xu[i] = SvNV(*av_fetch(upper, i, 0));
  }

  /* store the equation to pass to mock function */
  myparams.eqn = eqn;

  /* setup the random number generator */
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  /* setup the integrator */
  dim_size = (size_t)dim;
  calls_size = (size_t)calls;

  F.f = &integrand;
  F.dim = dim_size;
  F.params = &myparams;
  gsl_monte_miser_state *s = gsl_monte_miser_alloc(dim_size);

  gsl_monte_miser_integrate(&F, xl, xu, dim_size, calls_size, r, s, &res, &err);

  av_push(ret, newSVnv(res));
  av_push(ret, newSVnv(err));

  /* cleanup */
  gsl_monte_miser_free(s);
  gsl_rng_free(r);
  Safefree(xl);
  Safefree(xu);

  return ret;
}

MODULE = PerlGSL::Integration::MultiDim	PACKAGE = PerlGSL::Integration::MultiDim	

PROTOTYPES: DISABLE

AV *
c_int_multi (eqn, lower, upper, calls)
	SV*	eqn
	AV*	lower
	AV*	upper
	int	calls