#include "low_discrepancy.hpp"
#include <vector>
#include <gsl/gsl_cdf.h>
#include "var_alea.hpp"


double frac_part(double x);

//fonctions de compositions de monte_carlo.hpp

template <class _Result>
struct generator 
{
	typedef _Result result_type;
};

template <typename Fct, typename VA>
struct compose_t : public generator< typename Fct::result_type >
{
	compose_t(Fct f, VA X) : f(f), X(X) {};
	typename Fct::result_type operator()() { 
		return f(X()); 
	};
	private:
		Fct f; VA X;
};

template <typename Fct, typename VA>
inline compose_t<Fct, VA>
compose(Fct f, VA X) {
    return compose_t<Fct, VA>(f, X);
};

template <typename VA>
struct mean_t: public generator<double>
{
    mean_t(int n, VA X): n(n), X(X) {};
    double operator()(){
        double sum = 0;
        for(int i=0; i<n; i++){
            sum+=X();
        }
        return sum/n;
    }
    private :
        int n; VA X;
    };
    
template <typename VA>
inline mean_t<VA>
mean(int n, VA X){
    return mean_t<VA>(n, X);
};


//Les classes de monte-Carlo

template <typename L>
std::vector<double> quasi_mean(int n, L X, double p=0.05)
{
	std::vector<double> result(3,0);
	double x;
	for (int j = 0; j < n; j++) {
		x = X();
		result[0] += x;
		result[1] += x*x;
	}
	result[0] /= (double) n;
	result[1] = (result[1] - n*result[0]*result[0])/(double)(n-1);
	result[2] = gsl_cdf_gaussian_Pinv(1-p/2,1)*sqrt(result[1]/(double) n);
	return result;
};

template <typename L, typename Fct>
std::vector<double> quasi_mean(int n, Fct f, L X, double p=0.05)
{
	return quasi_mean(n, compose(f, X), p); 
};

//Les classes de quasi-mean

template <typename LDS>
struct quasi_gaussian : public var_alea<std::vector<double> >
{
    quasi_gaussian (int d, double mean = 0, double std =1)
        : d(d), mean(mean), std(std), s(d), U(0,1), seed(d) {
            for(int i=0; i<d; i++) {seed[i]=U();}
                };
    std::vector<double> operator ()(){
        std::vector<double> result = s();
        for(int i=0;i<d; i++){
            result[i] = mean + gsl_cdf_gaussian_Pinv(frac_part(result[i]+seed[i]), std);
        }
        return value = result;
    }
    private: 
    int d;
    double mean, std;
    LDS s;
    uniform U;
    std::vector<double> seed;
};



struct gaussian_d : public var_alea<std::vector<double> >
{
    gaussian_d(int d, double mean = 0, double std =1)
        : d(d), mean(mean), std(std), G(mean, std) {};
    std::vector<double> operator ()(){
        std::vector<double> result(d);
        for(int i=0;i<d; i++){
            result[i] = G();
        }
        return value = result;
    }
    private: 
    int d;
    double mean, std;
    gaussian G;
};