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

int main( int argc, char **argv){
	if(argc != 5){
		cout << "Usage :\n";
		cout << "\t Nx # of intervals ?\n";
		cout << "\t Correlation length?\n";
		cout << "\t Exponent?\n";
		cout << "\t KL Truncation level?\n";
		exit(0);
	}
	/* Read the problem paremater */
	unsigned N = atoi(argv[1]);		//# of spatial Cells
	double   L = atof(argv[2]);		//Correlation length
	double   P = atof(argv[3]);		//Exponent of the Covariance kernel
	unsigned T = atoi(argv[4]);		//Truncature of the KL expansion
	

	/* Create a pseudo random number generator and a normal distribution */
	std::default_random_engine generator;
	std::normal_distribution<double> normal(0.0,1.0);

	double dx = 1./(double)(N); 	// Spatial step
	MatrixXd K(N,N);				//Covariance matrix
	for(unsigned i=0; i<N; i++){
		for(unsigned j=0; j<N; j++){
			double d = dx*fabs(i-j)/L;			//Scaled distance
			K(i,j) = exp( - pow(d,P)/P )*dx;	//Two points correlation	
		}
	}
	SelfAdjointEigenSolver<MatrixXd> Eig(K);		//Decompose K
	MatrixXd Ev = Eig.eigenvectors()/sqrt(dx);		//Retrieve eigenvectors (and scale them) 
	VectorXd Lm = Eig.eigenvalues();				//Retrieve eigenvalues
	Lm.array() = Lm.array().abs();					//Move eigenvalues to their absolute value to ensure positivity even for close to zero value (round-off error)
	for(unsigned i=0; i<N/2; i++){		//Swap the eigenpairs to have the largest first
		VectorXd V = Ev.col(i);
		Ev.col(i) = Ev.col(N-1-i);
		Ev.col(N-1-i) = V;
		double l = fabs(Lm(i));
		Lm(i) = fabs(Lm(N-1-i));
		Lm(N-1-i) = l;
	}

	FILE *out = fopen("modes.gnu","w");				//Record the eigenmodes on disk
	for(unsigned i=0; i<N; i++){
		for(unsigned j=0; j<N; j++){
			fprintf(out,"%e %e %e\n",j*dx,Ev(j,i),Ev(j,i)*sqrt(Lm(i)));
			fprintf(out,"%e %e %e\n",(j+1)*dx,Ev(j,i),Ev(j,i)*sqrt(Lm(i)));
		}
		fprintf(out,"\n\n");
	}
	fclose(out);

	out = fopen("lambda.gnu","w");					//Record the eigenvalues on disk
	for(unsigned i=0; i<N; i++){
		fprintf(out,"%d %e\n",i,Lm(i));
	}
	fclose(out);

	out = fopen("realizations.gnu","w");			//Generate 100 realizations (truncated and non-truncated ones) and record on disk
	for(unsigned i=0; i<100; i++){	
		VectorXd U = VectorXd::Zero(N);
		for(unsigned j=0; j<T; j++) U += (Ev.col(j)*sqrt(fabs(Lm(j)))*normal(generator));  //truncated sum
		VectorXd V = U;
		for(unsigned j=T; j<N; j++) V += (Ev.col(j)*sqrt(fabs(Lm(j)))*normal(generator));	//continue for non-truncated sum
		for(unsigned j=0; j<N; j++){
			fprintf(out,"%e %e %e\n",j*dx,U(j),V(j));
			fprintf(out,"%e %e %e\n",(j+1)*dx,U(j),V(j));
		}
		fprintf(out,"\n\n");
	}
	fclose(out);

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