#include "define.h"
#include "Neur_Class.h"
#include "Corr_Stat.h"
#include <time.h>

int
main(int ac, char **av)
{
	time_t          start, end1, end2;
	double          dG, dT, gmutualoldpre, gmutualnewpre, gmutualoldpost,
	                gmutualnewpost;
	int             maxsz, minsz;
	int              maxold=0, maxnew, minold=0, minnew;
	time(&start);
	int32           seed = time(NULL);
	TRandomMersenne rg(seed);
	string          s1, s2, s3, s4, s5;

	if (ac <= 14) {
		puts("NN");
		puts("endtime");
		puts("Iinput external current I_ext Pre neuron");
		puts("Baseline current which is kept fixed, to decide the regime Idc post neuron");
		puts("self inhibition strength");
		puts("mutual inhibition strength");
		puts("decay constant for inhibition");
		puts("do you want the time series data");
		puts("noise level");
		puts("1 if Learning present, O otherwise");
		puts("output file name");
		puts("Self inhibition strength for pre neuron: Heterogeneity through gs");
		puts("Tau Rise time");
		puts("Spike time output file name");
		puts("Iteration Number for Noise Study");
		exit(0);
	}
	int             NN = atoi(av[1]);
	double          endtime = atof(av[2]);
	double          I_ext = atof(av[3]);
	double          Idc = atof(av[4]);
//	int             Synapse_Type = atoi(av[5]);
	double          gself = atof(av[5]);
	double          gmutual = atof(av[6]);
	double          taus = atof(av[7]);
	int             scale = atoi(av[8]);
	double          noise = atof(av[9]);
	int             Learnpresent = atoi(av[10]);
	string          out = av[11];
	double          gselfnew = atof(av[12]);
	double          tauR = atof(av[13]);
    int NoiseNum=atoi(av[14]);
	double          tauD = taus;

	double          tim=0.0, timestep = epsilon, Threshold = 0;

	double gelec=0.0;
	vector < double >x1, x2, w, *T;
	vector < double >x1p, x2p, wp, *TP;
	T = new vector < double >[3000];
	TP = new vector < double >[3000];
	
	int            *bol, *bolp;
	bol = new int   [NN + 10];
	bolp = new int  [NN + 10];

	

	for (int i = 0; i < NN; i++) {
		x1.push_back(0.);
		x2.push_back(0.);
		w.push_back(0.);
		x1p.push_back(0.);
		x2p.push_back(0.);
		wp.push_back(0.);
	}

	ofstream        Outfile, Outfile1, Outfile2;
Outfile.open(out.c_str(), ios::app);
	s1 = "Timeseries";
	s2 = DoubleToStdStr(I_ext);
	s3 = DoubleToStdStr(noise);
	s4 = IntToString(Learnpresent);
	s1 = s1 + "_" + s2 + "_" + s3 + "_" + s4 + ".dat";


	s5 = "Rates";
	s5 = s5 + "_" + s2 + "_" + s3 + "_" + s4 + ".dat";
	s5 = "RatesNetwork/" + s5;
//Outfile2.open((char *) s5.c_str(), ios: :out);

	HHneuron      **neurpre, **neurpost;
	neurpre = new HHneuron *[NN];
	neurpost = new HHneuron *[NN];
	for (int i = 0; i < NN; i++) {
		neurpre[i] = new HHneuron(.1, 1);
		neurpre[i]->parameter[11] = -65;
		neurpre[i]->parameter[7] = 0.;
		neurpre[i]->x[0] = -65;
		neurpost[i] = new HHneuron(.1, 1);
		neurpost[i]->parameter[11] = -65;
		neurpost[i]->parameter[7] = 0;
		neurpost[i]->x[0] = -55.23;
	}



	Insynapse     **Inputpre, **Inputpost;
	Inputpre = new Insynapse *[NN];
	Inputpost = new Insynapse *[NN];

	for (int i = 0; i < NN; i++) {
		Inputpre[i] = new Insynapse(neurpre[i], 0.0);
		Inputpost[i] = new Insynapse(neurpost[i], 0.0);
	}

	Couplingcurrent ***coupAB, ***coupBA;
	coupAB = new Couplingcurrent **[NN];
	coupBA = new Couplingcurrent **[NN];
	for (int i = 0; i < NN; i++) {
		coupAB[i] = new Couplingcurrent *[NN];
		coupBA[i] = new Couplingcurrent *[NN];
		for (int j = 0; j < NN; j++) {
			coupAB[i][j] =
				new Couplingcurrent(neurpre[i], neurpost[j], gelec / NN);
			coupBA[i][j] =
				new Couplingcurrent(neurpost[i], neurpre[j], gelec / NN);
		}

	}

	TwoDsynapse  ***slowself, ***slowselfA, ***slowmutual, ***slowmutualA;
	slowself = new TwoDsynapse **[NN];
	slowselfA = new TwoDsynapse **[NN];
	slowmutual = new TwoDsynapse **[NN];
	slowmutualA = new TwoDsynapse **[NN];
	for (int i = 0; i < NN; i++) {
		slowself[i] = new TwoDsynapse *[NN];
		slowselfA[i] = new TwoDsynapse *[NN];
		slowmutual[i] = new TwoDsynapse *[NN];
		slowmutualA[i] = new TwoDsynapse *[NN];
		slowself[i][i] = new TwoDsynapse(neurpre[i], neurpre[i], gselfnew, 0., 0);
		slowselfA[i][i] = new TwoDsynapse(neurpost[i], neurpost[i], gself, 0., 0);
		slowself[i][i]->parameter[1] = -82;
		slowselfA[i][i]->parameter[1] = -82;
		slowself[i][i]->parameter[2] = tauR;
		slowself[i][i]->parameter[3] = tauD;
		slowselfA[i][i]->parameter[2] = tauR;
		slowselfA[i][i]->parameter[3] = tauD;


		for (int j = 0; j < NN; j++) {
			slowmutual[i][j] = new TwoDsynapse(neurpre[i], neurpost[j], 1.0*gmutual / NN, 0., 0);
			slowmutualA[i][j] = new TwoDsynapse(neurpost[i], neurpre[j],  gmutual / NN, 0., 0);
			slowmutual[i][j]->parameter[1] = -82;
			slowmutualA[i][j]->parameter[1] = -82;

			slowmutual[i][j]->parameter[2] = tauR;
			slowmutual[i][j]->parameter[3] = tauD;
			slowmutualA[i][j]->parameter[2] = tauR;
			slowmutualA[i][j]->parameter[3] = tauD;

		}
	}
	int             count = 0;

	double         *curr_pre, *curr_post, sumprecur=0, sumpostcur=0, Het;
	curr_pre = new double[NN + 10];
	curr_post = new double[NN + 10];

	for (int i = 0; i < NN; i++) {
		curr_pre[i] = I_ext + 0.0 * rg.Random();
		curr_post[i] = Idc + 0.0 * rg.Random();
		sumprecur += curr_pre[i];
		sumpostcur += curr_post[i];
	}
	Het = 100 * (sumprecur - sumpostcur) / (sumprecur + sumpostcur);

	/*******************parameters for the TwoDsynapse to work********************/

	double         *ref_time_pre, *ref_time_post;
	double         *vnew_pre, *vnew_post, *vold_pre, *vold_post;
	double         *diff_pre, *diff_post;
	int            *bol_pre, *bol_post, *bol_postpre, *bol_prepost;
	ref_time_pre = new double[NN];
	ref_time_post = new double[NN];
	vnew_pre = new double[NN];
	vnew_post = new double[NN];
	vold_pre = new double[NN];
	vold_post = new double[NN];
	diff_pre = new double[NN];
	diff_post = new double[NN];
	bol_pre = new int[NN];
	bol_post = new int[NN];
	bol_prepost = new int[NN];
	bol_postpre = new int[NN];

	for (int i = 0; i < NN; i++) {
		vnew_pre[i] = 0;
		vnew_post[i] = 0;
		vold_pre[i] = neurpre[i]->x[0];
		vold_post[i] = neurpost[i]->x[0];
		bol_pre[i] = bol_post[i] = 0;
		bol_prepost[i] = bol_postpre[i] = 0;
	}

	/*****************************************************************************/
	gmutualoldpre = gmutualoldpost = gmutual;


double Hetgself=100*(gselfnew-gself)/(gselfnew+gself);


	
	int szpreold=0,szpostold=0;
	double sumcurr=0,sumcurrpost=0;
	vector<double> PreCurr,PostCurr;
	double RuleNoise=0;
	
	while (tim < endtime) {

		for (int i = 0; i < NN; i++) {
			Inputpre[i]->set_I(0);
			Inputpost[i]->set_I(0);
		}


		if (tim > 199 && tim < endtime) {
			for (int i = 0; i < NN; i++) {
				Inputpre[i]->set_I(curr_pre[i]);
				Inputpost[i]->set_I(curr_post[i]);
			}
		}
		/*if (tim > 199.25 && tim < 200.) {
			for (int i = 0; i < NN; i++) {
				Inputpre[i]->set_I(50.0);
				Inputpost[i]->set_I(50.0);
			}
		}*/
		if (noise != 0.0) {

			for (int i = 0; i < NN; i++) {
				do {
					x1[i] = 2.0 * rg.Random() - 1.0;
					x2[i] = 2.0 * rg.Random() - 1.0;
					w[i] = x1[i] * x1[i] + x2[i] * x2[i];
				}
				while (w[i] >= 1.0);

				do {
					x1p[i] = 2.0 * rg.Random() - 1.0;
					x2p[i] = 2.0 * rg.Random() - 1.0;
					wp[i] = x1p[i] * x1p[i] + x2p[i] * x2p[i];
				}
				while (wp[i] >= 1.0);

				neurpre[i]->AddNoise(noise, 0., 1., x1[i], x2[i], w[i],
						     timestep);
				neurpost[i]->AddNoise(noise, 0., 1., x1p[i], x2p[i], wp[i],
						      timestep);
			}
		}
		
		for (int i = 0; i < NN; i++) {
			HH_Run(&tim, neurpre[i], timestep);
			HH_Run(&tim, neurpost[i], timestep);
			TwoD_Run(&tim,slowself[i][i],neurpre[i]->x[0], timestep);
			TwoD_Run(&tim, slowselfA[i][i],neurpost[i]->x[0], timestep);
			for (int j = 0; j < NN; j++) {
				TwoD_Run(&tim,slowmutual[i][j],neurpre[i]->x[0], timestep);
				TwoD_Run(&tim, slowmutualA[i][j],neurpost[i]->x[0], timestep);
			}
		}

		if (tim > 200) {
			for (int i = 0; i < NN; i++) {
				spiketimes(&tim, &neurpre[i]->x[0], Threshold, bol[i], T[i]);
				spiketimes(&tim, &neurpost[i]->x[0], Threshold, bolp[i],
					   TP[i]);
			}
		}
		
		/*Getting the total current per cycle entering the neuron*/
		
		if (szpreold!=T[0].size())
		{
			//cout<<"Entered loop "<<T[0].size()-1<<endl;
		    PreCurr.push_back(sumcurr);
		    sumcurr=0;
		}
		sumcurr+=I_ext+slowself[0][0]->Isyn()+slowmutualA[0][0]->Isyn();
		szpreold=T[0].size();
		
				if (szpostold!=TP[0].size())
		{
			//cout<<"Entered loop "<<TP[0].size()-1<<endl;
		    PostCurr.push_back(sumcurrpost);
		    sumcurrpost=0;
		}
		sumcurrpost+=Idc+slowselfA[0][0]->Isyn()+slowmutual[0][0]->Isyn();
		szpostold=TP[0].size();

				
		
		/* Learning Rule On Line */
		if (Learnpresent == 1 ) {
			dG = 0;
			
			
			for (int u = 0; u < NN; u++)
				for (int v = 0; v < NN; v++) {

			
			
					maxsz = max(T[u].size(), TP[v].size());
					minsz = min(T[u].size(), TP[v].size());

					maxnew = maxsz;
					minnew = minsz;

					if (maxnew != maxold || minnew != minold) {
                          // if(TP[v].size()>1)
            			
            			if (TP[v].size()>0 && T[u].size()>0)				
                         dT = TP[v][TP[v].size() - 1] - T[u][T[u].size() - 1];
                         else
                         dT=-1000;
						//RuleNoise=rg.Random();
				//	dG = Rule(dT, 10.0, 1.0, 10.0, 1.0, .01);
						dG=Rule(dT,5,.25,5,.25,.015);
						

					}
					gmutualnewpre = gmutualoldpre + dG;
					if (gmutualnewpre < 0)
						gmutualnewpre = 0;
					gmutualnewpost = gmutualoldpost - dG;
					if (gmutualnewpost < 0)
						gmutualnewpost = 0;
					slowmutual[u][v]->parameter[0] = 1.0*gmutualnewpre;
					slowmutualA[u][v]->parameter[0] = gmutualnewpost;
				
				}

			maxold = maxnew;
			minold = minnew;

		} else;
		gmutualoldpre = 1.0*gmutualnewpre;
		gmutualoldpost = gmutualnewpost;
		count += 1;
		
		if (scale==2)
		{
			if (tim>0)
	Outfile<<tim<<" "<<neurpre[0]->x[0]<<" "<<neurpost[0]->x[0]<<" "<<gmutualnewpost<<" "<<gmutualnewpre<<endl;
		}
		
		tim += epsilon;
	}
 

 
cout<<tim<<" "<<neurpre[0]->x[0]<<" "<<neurpost[0]->x[0]<<" "<<gmutualnewpre<<" "<<gmutualnewpost<<" "<<slowmutualA[0][0]->x[0]<<endl;

 
        vector < double >Periodpre, Periodpost;
	double          Avgperiodpre, Avgperiodpost;
	double          sumpree = 0, sumposte = 0;
	double        **kura, sumkura = 0, avgkura;
	kura = new double *[NN + 10];
	for (int j = 0; j < NN + 10; j++)
		kura[j] = new double[NN + 10];

	for (int i = 0; i < NN; i++) {
		Periodpre.push_back(ISI(T[i]));
		Periodpost.push_back(ISI(TP[i]));
	}

	for (int i = 0; i < NN; i++) {
		sumpree += Periodpre[i];
		sumposte += Periodpost[i];
	}
	Avgperiodpre = sumpree / NN;
	Avgperiodpost = sumposte / NN;
	
//cout<<I_ext<<" "<<Avgperiodpre<<" "<<Avgperiodpost<<endl;

/*
Outfile1.open((char *)SpikeTimeFile.c_str(), ios::out);
if (scale!=0)
{
	for (int i=0;i<PreCurr.size();i++)
	Outfile1<<PreCurr[i]<<endl;
	Outfile1<<-1000<<endl;
	for (int i=0;i<PostCurr.size();i++)
	Outfile1<<PostCurr[i]<<endl;
	
for (int i=0;i<T[0].size();i++)
Outfile1<<T[0][i]<<endl;
Outfile1<<-1000<<endl;
for (int i=0;i<TP[0].size();i++)
{Outfile1<<TP[0][i]<<endl;}
}
Outfile1.close();
*/	
    
	for (int i = 0; i < NN; i++) {
		int             findpre = isnan(Periodpre[i]);
		int             findpreinf = isinf(1. / Periodpre[i]);
		for (int j = 0; j < NN; j++) {
			int             findpost = isnan(Periodpost[j]);
			if (findpre != 1 && findpost != 1 && findpreinf != 1
			    && Periodpost[j] != 0.) {
				kura[i][j] = Kuramoto(T[i], TP[j], Avgperiodpost);
				sumkura += kura[i][j];
			}
		}
	}
	avgkura = sumkura / (NN * NN);

	int             MaxTime, MinTime;

	MaxTime = max(T[0].size(), TP[0].size());
	MinTime = min(T[0].size(), TP[0].size());
	//cout << MaxTime << " " << MinTime << endl;
	double          diffspike = MaxTime - 2 * MinTime;
	//Outfile<<I_ext<<" "<<Avgperiodpre<<endl;
	
	if (diffspike < 100 && MinTime > 100) {
		cout<< NoiseNum<<" "<<gmutual<<" "<<I_ext << " " << Het << " " <<gmutualnewpre << " "<<gmutualnewpost<<" " <<avgkura<<" "<< Avgperiodpre << " " <<Avgperiodpost << " "  << Avgperiodpost / Avgperiodpre << endl;
	} 
	
	

	double         *rateA, *rateB;
	rateA = new double[NN + 10];
	rateB = new double[NN + 10];

	Outfile.close();
	Outfile1.close();
	
time(&end1);
double diff = difftime (end1, start);
cout<<"Computation time for "<< I_ext<<" is "<<diff<<" sec"<<endl;

	delete         *neurpre; delete *neurpost;delete *Inputpre;delete *Inputpost;delete **coupAB;
	   delete           **coupBA; delete **slowself; delete **slowselfA; delete **slowmutual;
	      delete        **slowmutualA;
	delete[] neurpre;
	delete[] neurpost;
	delete[] Inputpre;
	delete[] Inputpost;
	delete[] coupAB;
	delete[] coupBA;
	delete[] slowself;
	delete[] slowselfA;
	delete[] slowmutual;
	delete[] slowmutualA;
	delete[] ref_time_pre;
	delete[] ref_time_post;
	delete[] vnew_pre;
	delete[] vnew_post;
	delete[] vold_pre;
	delete[] vold_post;
	delete[] diff_pre;
	delete[] diff_post;
	delete[] bol_pre;
	delete[] bol_post;
	delete[] bol_postpre;
	delete[] bol_prepost;
        delete[] T;
        delete []TP;
        delete []bol;
        delete []bolp;
	return 0;
}