#include <algorithm>
#include <iostream>
#include "rqmc.hpp"

double pos (double x){
    return x>0?x:0;
}

struct asian_option : public std::unary_function<std::vector<double>, double>
{
    asian_option(double r, double T, double S0, double V, int d, double K)
        : r(r), T(T), S0(S0), V(V), d(d), K(K) {};
    
    double operator()(std::vector<double> X) const {
        std::vector<double> S(d);
        S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[0]);
        for(int i=1;i<d;i++){
            S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*X[i]);
        }
        double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
        return exp(-r*T)*pos(temp-K);
        };
        
    private:
        double r;
        double T;
        double S0;
        double V;
        int d;
        double K;
    };
        
        

struct asian_option_qmc : public var_alea<double>
{
	asian_option_qmc(double r, double T, double S0, double V, int d, double K)
        : r(r), T(T), S0(S0), V(V), d(d), K(K), G(d), U(0,1), seed(d) {
            for(int i=0; i<d; i++){
                seed[i]=U();
            }
        };
        
	double operator()() {
        std::vector<double> S(d);
        std::vector<double> sob(d); 
        sob = G();
        S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(seed[0]+sob[0]), 1));
        for(int i=1;i<d;i++){
            S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(seed[i]+sob[i]), 1));
        }
        double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
        return exp(-r*T)*pos(temp-K);
    };
    
	private:
        double r;
        double T;
        double S0;
        double V;
        int d;
        double K;
        sobol G;
        uniform U;
        std::vector<double> seed;
};


int main(){
    init_alea(1);
    asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
    //~ int M= 1000000;
    int N= 10000;
    //~ int I= 100;
    
    //~ std::vector<double> meanvar = monte_carlo(M, A);
    //~ std::cout<<"espérance "<<meanvar[0] <<" IC "<<meanvar[2]<<std::endl;
    //~ std::vector<double> temp;
    //~ double m = 0;
    //~ double s = 0;
    //~ for(int i=0;i<I;i++){
        //~ asian_option_qmc B(0.05, 1.0, 50.0, 0.1, 16, 45);
        //~ temp = monte_carlo(N,B);
        //~ m+=temp[0];
        //~ s+=temp[0]*temp[0];
        //~ }
    //~ m = m/I;
    //~ s = s/I - m*m;
    //~ std::cout<<"espérance "<<m <<" IC "<<sqrt(s)*1.96/10<<std::endl;
    int d = 16;
    std::vector<double> test(d,1);
    std::cout<<A(test)<<std::endl;
    gaussian_d G(16,0,1);
    
    std::vector<double> r = monte_carlo(N, A, G);
    for(int i=0; i<3; i++){
        std::cout<<r[i]<<std::endl;
    }

    
    return 0;
}