/*	gyak2_rejection_sampling.c
 *
 *	created by: racze
 *	date: 2011-09-25
 *
 * 	This program is an example for rejection sampling.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_rng.h>
/*!
 * Generates a sample of size <b>length</b> with density function
 * \f[
 *	f(x; x_0, \alpha, \delta) = \frac{\alpha+\delta}{\delta}
 *	\frac{\alpha}{x_0}
 *	\left(
 *	\frac{x}{x_0}
 *	\right)^{-(\alpha+1)}
 *	\cdot
 *	\left[
 *	1-
 *	\left(
 *	\frac{x}{x_0}
 *	\right)^{-\delta}
 *	\right]
 * \f]
 *The complementary distribution function is
 *\f[
 * 	F(x; x_0, \alpha, \delta) \equiv \mathbf{P}(X \geq x) = \frac{\alpha+\delta}{\delta}
 *      \left(
 *      	\frac{x}{x_0}
 *      \right)^{-\alpha}
 *      \cdot
 *      \left[
 *		1-\frac{\alpha}{\alpha+\delta}
 *      	\left(
 *			\frac{x}{x_0}
 *      	\right)^{-\delta}
 *      \right]
 *\f]
 *
 *\return None.
 */
void Generate_Corrected_Pareto(double* pointer, size_t length, double alpha, double delta, double xmin, const gsl_rng *rng);

/*!A comparator for the C library quicksort function (stdlib.h/qsort).
 *Descending order for doubles.
 */
int Hi2Lo(const void* a, const void* b);

int main(void){
	size_t N = 10000;
	size_t i;
	double data[N];

	/*Initialize for the GSL random number generator.*/
	const gsl_rng_type * T;
	gsl_rng * r;
     
	gsl_rng_env_setup();
     	/*Choose rng type.*/
	T = gsl_rng_default;
	/*Allocate memory to the rng.*/
	r = gsl_rng_alloc (T);
	/*Set the seed of the random number generator.*/
	gsl_rng_set(r, 123456);

	Generate_Corrected_Pareto(data, N, 2.0, 1.0, 10.0, r);

/*	for(i = 0; i < N; ++i){
		fprintf(stdout, "%g\n" , data[i]);
	}*/
	
	/*Free the memory allocated to the rng.*/
	gsl_rng_free(r);

	/*
 *	HOMEWORK: Check whether rejection sampling actually works:
 *	1) analytically
 *	and/or
 *	2) by constructing the empirical distribution function of a larger sample (say N=1000)
 *	and comparing it to the analytic function
 */

	qsort(data, N, sizeof(double), Hi2Lo);

	for(i = 0; i < N; ++i){
		fprintf(stdout, "%g\t%g\n", data[i], ((double) i+1)/((double) N));
	}


	return 0;	
}


void Generate_Corrected_Pareto(double* data, size_t length, double alpha, double delta, double xmin, const gsl_rng *rng){
        size_t index = 0;
        short accepted;
        double x, u;
        while(index < length){
                accepted = 0;
		fprintf(stderr, "------------- index = %zu -------------", index);
                while(!accepted){
                        /*1) Generate a random number distributed according to (x/xmin)^(-alpha)*/
                        u = gsl_rng_uniform(rng);
                        x = xmin*pow(u, -1/alpha);
                        /*2) Generate another random number to decide whether 
 *                         to accept or reject the previously generated value*/
                        u = gsl_rng_uniform(rng);
                        accepted = (u <  1.0 - pow(x/xmin, -delta));
			// accepted = 1;
			fprintf(stderr, "accepted = %d\n", accepted);
                }
		data[index++] = x;
        }
}

int Hi2Lo(const void* a, const void* b){
	double vala = *(double*) a;
	double valb = *(double*) b;

	if(vala > valb){
		return -1;
	}else if(vala < valb){
		return 1;
	}else{
		return 0;
	}
}