/*
 ============================================================================
 Name        : importance-sampling-vanilla.c
 Author      : Stippinger Marcell
 Version     : 1.0
 Copyright   : BME-TTK FI SzamSzim gyak segedanyag
 Description : Demonstration of importance sampling in C, Ansi-style
 ============================================================================
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Lattice dimension */
#define Npoint 1E3
#define SIGMA 1

void brute_force(gsl_rng * r)
{
	int i;
	double a, sum=0, sum2=0, prob, weight=0;

	for (i=0; i<=Npoint; i++) {
		a = 6.0*SIGMA*gsl_rng_uniform(r)-3.0*SIGMA; /* UNI(-3SIGMA,3SIGMA) */
		prob = gsl_ran_gaussian_pdf(a,SIGMA);

		weight += prob;
		sum += prob*a;
		sum2 += prob*a*a;
	}

	printf("Brute force  mean = %lf, variance = %lf\n",sum/weight,(sum2/weight-pow(sum/weight,2)));
}

void importance(gsl_rng * r)
{
	int i;
	double a, sum=0, sum2=0, prob, weight=0;

	for (i=0; i<=Npoint; i++) {
		a = gsl_ran_gaussian(r,SIGMA); /* N(0,SIGMA^2) */
		prob = 1;

		weight += prob;
		sum += prob*a;
		sum2 += prob*a*a;
	}

	printf("Gaussian  mean = %lf, variance = %lf\n",sum/weight,(sum2/weight-pow(sum/weight,2)));
}

int main(int argc, char* argv[], char* envp[])
{
	/* GSL RNG */
	const gsl_rng_type * T;
	gsl_rng * r;
	gsl_rng_env_setup();
	T = gsl_rng_mt19937;
	r = gsl_rng_alloc (T);

	brute_force(r);
	importance(r);

	/* Free the memory */
	gsl_rng_free (r);

	return 0;
}