#include <random>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <csmethodsserial.hpp>

using namespace std;

template<class AnyClass, class SystemInfo> /*Algorithm can be used for a class that defines any system*/
void cloning(double tmax, int nc, vector<int> dim, double dt, double s, char obs, char proc, vector<double> rates, int iseed, double repeats, int clonemethod)
{
        void (*selectclones)(int*, double*, int, double, SystemInfo**, int);
        if(clonemethod==1)
        {
          selectclones = csm1serial<SystemInfo*>;
        }
	if(clonemethod==2)
	{
	  selectclones = csm2serial<SystemInfo*>;
	}

        int sLL = int(round(s*double(dim[0])*double(dim[0])));
	int logt = int(round(log10(tmax)));
	int lognc = int(round(log10(nc)));
	int dtprint = int(round(dt));
	
	FILE * pFile;
	char FILENAME[50];
	sprintf(FILENAME,"%dSerial%d,%d,%d,%d,%d,%d.txt",iseed,clonemethod,logt,dim[0],lognc,sLL,dtprint);
	pFile = fopen(FILENAME,"w");
	timespec start, end;
        double repcount;
	for(repcount=0;repcount<repeats;repcount++)
	{
	clock_gettime(CLOCK_REALTIME, &start);  

	int iter;

	/*Initialise Variables*/
	int i, k; //For loop indices
	double j;
	double r; //Random number for picking lattices

	/*Generate seeds for each system's internal twister*/
	int *seeds; //Seeds for lattices
	seeds = new int[nc]();
	mt19937 twister(iseed); //Twister for generating seeds
	std::uniform_int_distribution<int> dist(0, pow(10, 7)); //Distribution for twister seeds
	for (i = 0; i < nc; i++) //Initialise seeds to be in order
	{
		seeds[i] = dist(twister);
	}
	///***Variables for cloning process***///
	double *upsilon;
	upsilon = new double[nc](); //Holds cloning/pruning rates on lattices local to this processor
	int *c;
	c = new int[nc](); //Indices of systems to clone

	//Note: Two sets of clones prevent over-writing during cloning
	int arrswitch = 0; //Switch between two sets of lattices so one one not overwritten during cloning
	AnyClass *systems[2 * nc]; //Vector that holds two sets of lattices
	SystemInfo *infos[2 * nc]; //Vector that holds information about observables
	mt19937 *twisters[nc];
	double totups = 0; //Holds total of cloning/pruning factors
	double X = 1; //Holds product of cloning factors
	double lnX = 0; //Log value easier to store

	/*'Key' is ued in the binary search*/
	double *key; //Holds cumulative upsilon values to be used as a key in a binary search
	key = new double[nc + 1]();
	key[0] = 0; //First element of key is always 0

	/*Generate Systems, Info Classes and Random Number Generators*/
	for (i = 0; i < nc; i++)
	{
       		twisters[i] = new mt19937(seeds[i]);
       		systems[i] = new AnyClass(dim); //Fill vector a with newly constructed lattices
       		systems[i + nc] = new AnyClass(dim); //Fill vector b with newly constructed lattices
       		infos[i] = new SystemInfo(dim, s, obs, proc, rates, twisters[i]); //Fill vector a with newly constucted classes that hold information about processes
       		infos[i + nc] = new SystemInfo(dim, s, obs, proc, rates, twisters[i]); //Fil vector b with info classes
    	}
        
       	///***Initialise each lattice***///
	for (i = 0; i < nc; i++)
	{
	        infos[i]->init(systems[i]);
		//Note: Only lattices top row are initialised. Lattices in set B carry on process from final particle positions of particles in lattices in set A
	}

        for (j = 0; j < tmax / dt; j++) //At each time step
	{
	        totups = 0; //Reset the sum of upsilons to 0
		for (i = 0; i < nc; i++) //For each lattice
		{
		               (infos[(arrswitch*nc)+i]->*(infos[(arrswitch*nc)+i]->dynamics))(dt, systems[(arrswitch*nc)+i]); //Run Dynamics for one cloning interval
			       (infos[(arrswitch*nc)+i]->*(infos[(arrswitch*nc)+i]->calcups))(&upsilon[i]); //Calculate Cloning Factor Upsilon Y
	    	}

		        for (i = 0; i < nc; i++)
			{
				key[i + 1] = key[i] + upsilon[i]; //Generate key for selecting clones
			}
			totups = key[nc]; //Extract total of upsilons
			X = double(totups) / double(nc); //Calculate cloning factor and store as log to avoid infinitely large values
			lnX = lnX + log(X);

		selectclones(c, key, nc, totups, infos, arrswitch); 

		for (i=0; i<nc; i++)
		{
		  *systems[(1 - arrswitch)*nc + i] = *systems[(arrswitch*nc) + c[i]]; //Clone System
		}

		arrswitch = 1 - arrswitch; //Set the other set of systems to be used in next time loop		 
	}

	/*Calculate psi for this value of s*/
	double psi;
	psi = double(lnX) / double(j*dt);
	fprintf(pFile,"\n%f", psi);
	clock_gettime(CLOCK_REALTIME, &end);
	fprintf(pFile,"\n%f", float((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/float(pow(10.0,9.0))));

	/*Delete Arrays*/
	delete upsilon;
	delete key;
	delete seeds;

	for(i=0;i<nc;i++)
	{
	  delete infos[i];
	  delete infos[nc+i];
	  delete systems[i];
	  delete systems[nc+i];
	  delete twisters[i];
	}

	iseed+=1;
        }
}

int main(int argc, char** argv)
{
	int T = atoi(argv[1]);
	int nc = atoi(argv[2]);
	std::vector<int> dim;
	int L = atoi(argv[3]);
	dim.push_back(L);
	int N = 0.5*L;
	dim.push_back(N);
	double dt = atof(argv[4]);
        double s = double(atof(argv[5]))/double(L*L);
	std::vector<double> rates;
	rates.push_back(1);
	rates.push_back(1);
	double seed = double(atof(argv[6]));
	double repeats = atof(argv[7]);
	int clonemethod = atoi(argv[8]);
	printf("Serial Code: %d units of time, %d clones %d sites %d particles %f repeats, s = %f and Cloning Interval of %f units of time",T,nc,L,N,repeats,s,dt);
        cloning<lattice, latticeinfo>(T, nc, dim, dt, s, 'h', 'A', rates, seed, repeats, clonemethod);
}
