#include <iostream>
#include <iomanip>
#include <fstream>
#include <Eigen/Dense>
#include <random>
using namespace std;
using namespace Eigen;

int main( int argc, char **argv){
	if(argc != 3){
		cout << "Usage :\n";
		cout << "\t Max # of samples ?\n";
		cout << "\t order of H?\n";
		exit(0);
	}

	/* Read the problem paremater */
	unsigned N = atoi(argv[1]);				//Number of MC samples
	unsigned P = atoi(argv[2]); P++;		//Max order P of the Taylor expansion Hp of exp(x)
 
	std::default_random_engine generator;					//Create pseudo random number generator
	std::uniform_real_distribution<double> Unif(0,1.);		//and uniform distribution in [0,1]

	VectorXd EofHP_ex = VectorXd::Zero(P);						//Vector of exact means of Hp for p=0, P
	EofHP_ex(0) = 1;
	double xx = 1.;
	for(unsigned p=1; p<P; p++)	{
		xx /= (double)(p+1); 
		EofHP_ex(p) = EofHP_ex(p-1) + xx;
	}
	cout << "E(H_P) : " << EofHP_ex.transpose() << endl;
	
	
	double EofF_ex = exp(1)-1.;								//Exact expectation of F
	VectorXd MCSum = VectorXd::Zero(P);						//Accumulator of samples of (F-Hp) for p=0, P

	FILE *out = fopen("res.gnu","w");
	unsigned next =10;

	for(unsigned i=0; i<N; i++){
		double x = Unif(generator);						//Draw x from [0,1]
		double F = exp(x);								//evaluate F
		VectorXd Hp(P);									//Compute the Hp at x
		Hp(0) = 1;
		double xp = 1;
		for(unsigned p=1; p<P; p++){
			xp *= x/(double)(p);
			Hp(p) = Hp(p-1) + xp;
		}

		MCSum += (VectorXd::Ones(P)*F - Hp);			//Add sample value of (F-Hp) to the accumulator		
		if(i>next){											//Report convergence
			VectorXd Error = MCSum/(double)(i) + EofHP_ex;	//MC estimator E(F-Hp) + E(Hp)	
			fprintf(out,"%d ",i);
			for(unsigned p=0; p<P; p++) fprintf(out," %e ",fabs(Error(p)-EofF_ex)); 	//MC estimator error

			/* case of uniform discretization of [0,1] with i intervals */
			double dx = 1./(double)(i);		//interval size
			double sum_grd = 0;				//initialize integral
			for(unsigned j=0; j<i; j++){							//loop over integration grid
				sum_grd += exp(dx*(j+.5));								//compute integral
			}
			sum_grd *=dx;						//scale the result	

			/* case of the LHS (select one point at random over each of the i interval*/
			double sum_lhs = 0;										//initialize sum
			for(unsigned j=0; j<i; j++){							//loop over intervals
				sum_lhs += exp(dx*(double(j+Unif(generator))));		//compute sum
			}
			sum_lhs /=(double)(i);									//LHS estimate (mean value)

			fprintf(out," %e %e ",fabs(sum_grd-EofF_ex),fabs(sum_lhs-EofF_ex));		//Report results
			fprintf(out,"\n");
			fflush(out);				
			next *=1.3;										//set next iteration to report convergence
		}	
	}

	fclose(out);

  	cout << "Progam done " << endl;
};