/* VEGAS BUG */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

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

struct my_function_struct {
    double (*function)(double * x, size_t size, void *params);
    void  *  params;
};

typedef struct my_function_struct my_function;

struct i_params { double * x; size_t size; my_function * F;};
struct gaussian_params { double a; double b; };

/* Proto-types */
double i_transform(double * t, size_t size, void * params);
double gaussian(double * x, size_t size, void * params);

/* Functions */

/* Tranform integrand to the (-1,1) domain */
double i_transform(double * t , size_t size, void * params)
{
    struct i_params * transform_params = (struct i_params * ) params;
    my_function * F = transform_params->F;
    double  * x     = transform_params->x;
    size_t i;
    double theta, y, dx;

    dx = 1.0;
    for(i=0;i<size;i++) {
        theta = M_PI_2*t[i];
        x[i]  = tan(theta);
        dx   *= (M_PI_2)*(1.0/cos(theta)/cos(theta));
    }

    y = F->function(x,size,F->params);
    y *= dx;

    return y;
}

    
double gaussian(double * t, size_t size, void * params)
{
    size_t i;
    double y, sum;
    struct gaussian_params *g_params = ( struct gaussian_params * ) params; 

    sum = 0.0;
    for(i=0;i<size;i++) {
        sum+=(t[i]*t[i]);
    }

    y = g_params->b*exp(-1.0*g_params->a*sum);

    return y;

}
  
void display_results(const char * title, double res, double est_error, double  ans)
{
    printf("%s ====================\n", title);
    printf("result     = %1.9E\n", res);
    printf("est error  = %1.9E\n", est_error);
    printf("answer     = %1.9E\n", ans);
    printf("true error = %1.9E\n", fabs(res-ans));
    printf("sigma      = %.3f\n", fabs(res-ans)/est_error);
}
    
int main (int argv, char * argc[] ) 
{

    size_t  i,lcnt;

    size_t  num_dim = 1;
    double  ans     = 0.5000000000000000;
    double  norm    = pow(M_SQRTPI,num_dim);
    double  *xl, *xu;
    double  result, error;

    size_t calls = num_dim*500000;
    size_t MAX_LOOP_CNT = 10;

    struct i_params        trans_params;
    struct gaussian_params F_params;
    my_function  F;

    gsl_monte_function gsl_F;

    const gsl_rng_type *type;
    gsl_rng      *rng;

    gsl_monte_vegas_state * state;

    /* Set up lower and upper limits of integration */
    xl = (double *)malloc(sizeof(double)*num_dim);
    xu = (double *)malloc(sizeof(double)*num_dim);
    for ( i=0; i<num_dim; i++ ) {
        xl[i] = 0.0;
        xu[i] = 1.0;
    }

    /* Set up the RNG */
    gsl_rng_env_setup();
    type = gsl_rng_default;
    rng =  gsl_rng_alloc(type);
    
    /* Set up the function */
    F_params.a    = 1.0;
    F_params.b    = 1.0/norm;
    
    F.function    = &gaussian;
    F.params      = &F_params;

    trans_params.x    = (double *) malloc(sizeof(double)*num_dim);
    trans_params.size = num_dim;
    trans_params.F    = &F;

    gsl_F.f       = &i_transform;
    gsl_F.params  = &trans_params;
    gsl_F.dim     = num_dim;
    
    state = gsl_monte_vegas_alloc(num_dim);
    state->verbose = 1;
    /* Warm-up */
    gsl_monte_vegas_integrate(&gsl_F,xl,xu,num_dim,calls/20,rng,state,&result,&error);
    display_results("After warm up",result,error,ans);
    printf("chisq/dof = %.3f\n", state->chisq);
    lcnt = 0;
    do {
        gsl_monte_vegas_integrate(&gsl_F,xl,xu,num_dim,calls,rng,state,&result,&error);
        printf("result = %1.9E sigma = %1.9E chisq/dof = %.3f\n", result, error, state->chisq);
        lcnt++;
    }
    while ( fabs(1.0-state->chisq) > 0.5 && lcnt < MAX_LOOP_CNT );

    display_results("Final results", result, error,ans);

    /* Free up memory */
    free(trans_params.x);
    free(xl);
    free(xu);
    gsl_monte_vegas_free(state);

    return 0;
}
