#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <fstream>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctime>
#include <random>  




time_t seconds;

//Created July 2, 2018
//Last updated November 2, 2018



using namespace std;


char FileRoot[50], TraitData[50] = "Traits_", BinData[50] = "Bins_", MetaData[50] = "Meta_", PosteriorData[50] = "Posterior_", ParameterData[50] = "Parameters_", TrueParameterData[50] = "TrueParameters_", AccuracyData[50] = "Accuracy_", MPDAccuracyData[50] = "MPDAccuracy_",VisualizeData[50]="Visualize_";
char fileinA[50], fileinP[50];
double hX,hY,Gx, Gy,Vx,Vy, mx, my,Nx,Ny, Mx[500][500], My[500][500], gx, gy,ETx,ETy,VTx,VTy, Tx[500], Ty[500],EAx,EAy,VAx,VAy,Ax[500],Ay[500];
double MetaMeanX, MetaMeanY, MetaVarX, MetaVarY, MetaCov, MetaSDX, MetaSDY, MetaRo;
double ThreshEx,ThreshEy,ThreshSDx,ThreshSDy,ThreshRo, DataMeanX, DataMeanY, DataSDX, DataSDY, DataRo;
int rep,MaxRuns,Hits,RunType,PosteriorSize,Runs, Succeeds;
double EAXStore[20000], EAYStore[20000], VAXStore[20000], VAYStore[20000];




// read X[P] as the population mean phenotype for trait X in population P. 
double x,y,X[500], Y[500], Xp[500], Yp[500];
int i,j,Npops,MigMod,G,test,Equilibrated, PopsSampled;
double Wx, Wy,EWx,EWy, Fx, Fy, WXB, WYB;
double XMin,XMax,XStep,YMin,YMax,YStep,SD,Acceptance;
double DataEAx, DataEAy,DataVAx,DataVAy, DataETx, DataETy,DatahX,DatahY,Datamx,Datamy,DataNx,DataNy,Datagx,Datagy,DataVTx,DataVTy;
double StoreMetaMeanX[2000], StoreMetaMeanY[2000], StoreMetaVarX[2000], StoreMetaVarY[2000], StoreMetaCov[2000], StoreMetaSDX[2000], StoreMetaSDY[2000], StoreMetaRo[2000];



void DrawPriorPars();
void MakeMigMat();
void MakeIndividuals();
void Move();
void Selection();
void Drift();
void TraitOutput();
void MetaOutput();
void Mate();
void MetaStats();
void TestSummary();
void ProcessPosterior();





ofstream out_traits;
ofstream out_meta;
ofstream out_posterior;
ofstream out_pars;
ofstream out_Bins;

#define _CRT_SECURE_NO_WARNINGS


double Big_Rand(double MIN, double MAX)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	uniform_real_distribution<> dist(MIN, MAX); // distribute results between 1 and 6 inclusive.  
	return dist(gen);
}

int Fish_Rand(double p)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	poisson_distribution<> distr(p); // distribute results between 1 and 6 inclusive.  
	return distr(gen);
}

double Norm_Rand(double Mu,double SD)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	normal_distribution<> distr(Mu,SD); // 
	return distr(gen);
}

double LogNorm_Rand(double p, double q)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	lognormal_distribution<> distr(p, q); // p is Alpha (shape) and q is Beta (rate)
	return distr(gen);
}

double Gam_Rand(double p, double q)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	gamma_distribution<> distr(p, q); // p is Alpha (shape) and q is Beta (rate)
	return distr(gen);
}

int main()

{
	srand(time(NULL));

	cout << "Enter root file name\n";
	cin >> FileRoot;
//	The code is sufficiently general to run any possible model of migration, however, we have currently studied performance only under the island model. Thus, the option below is commented out.
//	cout << "Enter migration model: (1) for Island Model or (2) for Stepping Stone\n";
//	cin >> MigMod;
	MigMod = 1;

	out_traits.open(strcat(TraitData, FileRoot));
	out_meta.open(strcat(MetaData, FileRoot));
	out_posterior.open(strcat(PosteriorData, FileRoot));
	out_pars.open(strcat(ParameterData, FileRoot));
	out_Bins.open(strcat(BinData, FileRoot));


	out_posterior << "Run,EAX,EAY,VAX,VAY,gx,gy,hx,hy,mx,my,Vx,Vy,Nx,Ny,ETx,ETy,VTx,VTy\n";
	out_pars << "Run,PopsSampled,EAX,EAY,VAX,VAY,gx,gy,hx,hy,mx,my,Vx,Vy,Nx,Ny,ETx,ETy,VTx,VTy\n";
	out_meta << "Run, G , MetaMeanX ,MetaMeanY,MetaVarX, MetaVarY , MetaCov\n";


	//Here I have manually entered the summary statistics describing the Toju and Sota data
	DataMeanX = 10.7;
	DataMeanY = 8.37;
	DataSDX = 1.38;
	DataSDY = 2.51;
	DataRo = 0.87;
	//Done SETTING THE TARGET DATA

	


	Npops = 13;//This is set to the number of populations within the Honshu Clade from Toju and Sota
	PopsSampled = Npops;

	//Set responsive thresholds as a function of the data set
	ThreshEx = 0.05 + 0.1*DataMeanX;
	ThreshEy = 0.05 + 0.1*DataMeanY;
	ThreshSDx = 0.05 + 0.2*DataSDX;
	ThreshSDy = 0.05 + 0.2*DataSDY;
	ThreshRo = 0.05 + 0.2*DataRo;
	//Done setting responsive thresholds

	Hits = 0;
	PosteriorSize = 2000;
	MaxRuns = 604000;
	Runs = 0;
	while (Hits < PosteriorSize&&Runs<=MaxRuns)//run until the posterior contains PosteriorSize points. Stop if more than MaxRuns have been executed or if no hits occur within the first X tries. 
	{
		cout <<Runs<<","<< Hits << "\n";
			


		DrawPriorPars();
		out_pars << "Prior(" << Runs << ")," << PopsSampled << "," << EAx << "," << EAy << "," << VAx << "," << VAy << "," << gx << "," << gy << "," << hX << "," << hY << "," << mx << "," << my << "," << Vx << "," << Vy << "," << Nx << "," << Ny << "," << ETx << "," << ETy << "," << VTx << "," << VTy << "\n";
		MakeIndividuals();
		MakeMigMat();

		Equilibrated = 0;
		for (G = 0; G <=500; G++)//Main generational recursion
		{
			Selection();
			Drift();
			Move();
			Mate();//NOTE: MATING MUST be the final event of the life cycle!
			MetaStats();

			if (G == 500)
			{
				cout << "Equilibriation failure\n";
			}

		}//End generational loop

		
		MetaOutput();
		TestSummary();
		Runs++;
		
	}//End of posterior generating loop
	out_meta.flush();
	out_pars.flush();
	out_posterior.flush();

	
	if (Runs < MaxRuns)
	{
		Succeeds++;
	}
	ProcessPosterior();





}


void DrawPriorPars()
{
	double k, M, Mx, My;//k defines the shape and M the mode of the Gamma Distributions

	//Draw phenotypic variances from a Gamma distribution with Shape k and mode M
	k = 50.0;
	M = 0.71;
	Vx = Gam_Rand(k, M / (k - 1.0));
	k = 50;
	M = 3.72;
	Vy = Gam_Rand(k, M / (k - 1.0));


	//Draw heritabilities from a uniform distribution since this is a {0,1} quantity
	hX = Big_Rand(0.312, 0.512);
	hY = Big_Rand(0.6375, 0.8375);


	//Draw the strengths of stabilizing selection from a gamma 
	k = 3.0;// 10.0;
	M = 0.05;
	gx = Gam_Rand(k, M / (k - 1.0));

	k = 3.0;// 10.0;
	M = 0.05;
	gy = Gam_Rand(k, M / (k - 1.0));

	//Draw local population sizes from a Gamma distribution with Shape k and mode M
	k = 200.0;// 10.0;
	M = 30870;
	Nx = Gam_Rand(k, M / (k - 1.0));

	k = 100.0;// 10.0;
	M = 1818;
	Ny = Gam_Rand(k, M / (k - 1.0));

	//Draw rates of movement from a Gamma distribution with Shape k and mode M
	k = 2.0;
	M = 0.00000307;
	mx = Gam_Rand(k, M / (k - 1.0));

	k = 2.0;
	M = 0.000001;
	my = Gam_Rand(k, M / (k - 1.0));

	//Draw the spatial average strengths of coevolutionary selection
	EAx = Big_Rand(0.0, 3.0);
	EAy = Big_Rand(0.0, 3.0);

	//Draw the spatial average phenotypic optima favored by stabilizing selection from a Gamma distribution with Shape k and mode M
	k = 30.0;
	M = 5.911818;
	ETx = Gam_Rand(k, M / (k - 1.0));

	k = 30.0;
	M = 6.248235;
	ETy = Gam_Rand(k, M / (k - 1.0));

	VAx = 0;
	VAy = 0;
	VTx = 0;
	VTy = 0;

	//Although the simulation code can accomodate spatial variation in coevoltuionary sensitivities, these are set to zero by default.
	if (EAx>0)
	{
		VAx = 0.0;// Big_Rand(0.0, 0.02);
	}
	if (EAy>0)
	{
		VAy = 0.0;// Big_Rand(0.0, 0.02);
	}

	//Draw the spatial variance in phenotypic optima favored by stabilizing selection from a Gamma distribution with Shape k and mode M
	if (ETx > 0)
	{
		k = 1.5;
		Mx = 0.15;
		VTx = Gam_Rand(k, Mx / (k - 1.0));
	}
	else
	{
		VTx = 0.0;
	}
	if (ETy > 0)
	{

		k = 2;
		My = 1.27;
		VTy = Gam_Rand(k, My / (k - 1.0));
	}
	else
	{
		VTy = 0.0;
	}


	for (i = 0; i < Npops; i++)
	{
		if (VTx > 0)
		{
			Tx[i] = -1;
			while (Tx[i] < 0)//Make sure the Abiotic optimum is positive
			{
				Tx[i] = Norm_Rand(ETx, sqrt(VTx));
			}
		}
		else
		{
			Tx[i] = ETx;
		}

		if (VTy > 0)
		{
			Ty[i] = -1;
			while (Ty[i] < 0)//Make sure the abiotic optimum is positive
			{
				Ty[i] = Norm_Rand(ETy, sqrt(VTy));
			}
		}
		else
		{
			Ty[i] = ETy;
		}

		if (VAx > 0)
		{
			Ax[i] = -1;
			while (Ax[i] < 0)//Make sure the strength of coevolution is positive
			{
				Ax[i] = Norm_Rand(EAx, sqrt(VAx));
			}
		}
		else
		{
			Ax[i] = EAx;
		}

		if (VAy > 0)
		{
			Ay[i] = -1;
			while (Ay[i] < 0)//Make sure the strength of coevolution is positive
			{
				Ay[i] = Norm_Rand(EAy, sqrt(VAy));
			}
		}
		else
		{
			Ay[i] = EAy;
		}

	}
}



void MakeIndividuals()
{
	for (i = 0; i < Npops; i++)
	{
		X[i] = Tx[i] +1.0*rand() / (1.0*RAND_MAX);
		Y[i] = Ty[i] +1.0*rand() / (1.0*RAND_MAX);
	}
}


void MakeMigMat()
{
	if (MigMod == 1)//Run Island Model
	{
		for (i = 0; i < Npops; i++)
		{
			for (j = 0; j < Npops; j++)
			{
				if (i != j)
				{
					Mx[i][j] = mx/(1.0*Npops-1);
					My[i][j] = my/(1.0*Npops-1);
				}
				if (i == j)
				{
					Mx[i][j] = 1.0-mx ;
					My[i][j] = 1.0-my ;
				}
			}
		}
	}

	if (MigMod == 2)//Run Stepping Stone Model
	{
		for (i = 0; i < Npops; i++)
		{
			for (j = 0; j < Npops; j++)
			{
				if (i == 0)
				{
					if (j == 0)
					{
						Mx[i][j] = 1.0-mx/2.0;
						My[i][j] = 1.0-my/2.0;
					}
					if (j == 1)
					{
						Mx[i][j] = mx/2.0;
						My[i][j] = my/2.0;

					}
					if (j >1)
					{
						Mx[i][j] = 0.0;
						My[i][j] = 0.0;

					}
				}
				if (i == Npops-1)
				{
					if (j == Npops - 1)

					{
						Mx[i][j] = 1.0 - mx / 2.0;
						My[i][j] = 1.0 - my / 2.0;
					}
					if (j == Npops - 2)
					{
						Mx[i][j] = mx /2.0;
						My[i][j] = my /2.0;
					}
					if (j < Npops - 2)
					{
						Mx[i][j] = 0.0;
						My[i][j] = 0.0;
					}
				}
				if (i >0&&i<Npops-1)
				{
					if (i == j)
					{
						Mx[i][j] = 1.0 - mx;
						My[i][j] = 1.0 - my;
					}
					if (j==i - 1)
					{
						Mx[i][j] = mx / 2.0;
						My[i][j] = my / 2.0;
					}
					if (j==i + 1)
					{
						Mx[i][j] = mx / 2.0;
						My[i][j] = my / 2.0;
					}
					if (j > i + 1)
					{
						Mx[i][j] = 0.0;
						My[i][j] = 0.0;
					}
					if (j < i - 1)
					{
						Mx[i][j] = 0.0;
						My[i][j] = 0.0;
					}
				}
			}
		}
	}
}


void Move()
{
	double XPM[500], YPM[500];
	for (i = 0; i < Npops; i++)
	{
		XPM[i] = 0;
		YPM[i] = 0;
		for (j = 0; j<Npops ; j++)
		{
			XPM[i] = XPM[i] + Mx[i][j] * Xp[j];//Move individuals according to the migration matrix Mx[Target][Source];
			YPM[i] = YPM[i] + My[i][j] * Yp[j];//Move individuals according to the migration matrix My[Target][Source];
		}
	}

	for (i = 0; i < Npops; i++)
	{
		Xp[i] = XPM[i];
		Yp[i] = YPM[i];
	}

}

void Selection() //employ a five sigma rule for integration range... 
{
	int k,NSteps;

	
	
	//Currently set to run a phenotypic differences or arms race model. This is approprite for the weevils and camellias. Running other forms of interaction model is straightforward and just involves replacing the escalation functions below with any other sort of function (e.g., matching)
	NSteps = 50;

	for (i = 0; i < Npops; i++)
	{
		XMin = X[i] - 5 * (sqrt(Vx));
		XMax = X[i] + 5 * (sqrt(Vx));
		XStep = (XMax - XMin) /(1.0*NSteps);
		YMin = Y[i] - 5 * (sqrt(Vy));
		YMax = Y[i] + 5 * (sqrt(Vy));
		YStep = (YMax - YMin) / (1.0*NSteps);
	
		//Calculate mean fitness for species X
		WXB = 0;
		y = YMin;
		while(y<YMax)
		{
			x = XMin;
			while (x<XMax)
			{
				Wx = (exp(-gx*(pow(x - Tx[i], 2.0)))) / (1.0 + exp(-Ax[i]*(x - y)));
				Wy = (exp(-gy*(pow(y - Ty[i], 2.0)))) / (1.0 + exp(-Ay[i] *(y - x)));
				Fx = (exp(-pow(x - X[i], 2.0) / (2.0*Vx))) / (sqrt(2.0*3.1415*Vx));
				Fy = (exp(-pow(y - Y[i], 2.0) / (2.0*Vy))) / (sqrt(2.0*3.1415*Vy));

				//cout << Wx << "," << Wy << "," << Fx << "," << Fy << "\n";
				//cin >> test;

				WXB = WXB + XStep*YStep*Fx*Fy*Wx;
				x = x + XStep;
			}
			y = y + YStep;
		}
		//Calculate post-selection mean phenotype for species X

		Xp[i] = 0;
		x = XMin;
		while(x<XMax)
		{
			y = YMin;
			EWx = 0;
			while(y<YMax)
			{
				Fy = (exp(-pow(y - Y[i], 2.0) / (2.0*Vy))) / (sqrt(2.0*3.1415*Vy));
				Wx = (exp(-gx*(pow(x - Tx[i], 2.0)))) / (1.0 + exp(-Ax[i] *(x - y)));
				EWx = EWx + YStep*Fy*Wx;
				y = y + YStep;
			}

			Fx = (exp(-pow(x - X[i], 2.0) / (2.0*Vx))) / (sqrt(2.0*3.1415*Vx));
			Xp[i] = Xp[i] + XStep*x*Fx*EWx / WXB;
			x = x + XStep;
		
		}


		//Calculate mean fitness for species Y
		WYB = 0;

		y = YMin;
		while(y<YMax)
		{
			x = XMin;
			while(x<XMax)
			{
				Wx = (exp(-gx*(pow(x - Tx[i], 2.0)))) / (1.0 + exp(-Ax[i] *(x - y)));
				Wy = (exp(-gy*(pow(y - Ty[i], 2.0)))) / (1.0 + exp(-Ay[i] *(y - x)));
				Fx = (exp(-pow(x - X[i], 2.0) / (2.0*Vx))) / (sqrt(2.0*3.1415*Vx));
				Fy = (exp(-pow(y - Y[i], 2.0) / (2.0*Vy))) / (sqrt(2.0*3.1415*Vy));



				WYB = WYB + XStep*YStep*Fx*Fy*Wy;
				x = x + XStep;
			}
			y = y + YStep;
		}
		//Calculate post-selection mean phenotype for species Y

		Yp[i] = 0;
		y = YMin;
		while(y<YMax)
		{
			x = XMin;
			EWy = 0;
			while(x<XMax)
			{
				Fx = (exp(-pow(x - X[i], 2.0) / (2.0*Vx))) / (sqrt(2.0*3.1415*Vx));
				Wy = (exp(-gy*(pow(y - Ty[i], 2.0)))) / (1.0 + exp(-Ay[i] *(y - x)));
				EWy = EWy + XStep*Fx*Wy;
				x = x + XStep;
			}
			Fy = (exp(-pow(y - Y[i], 2.0) / (2.0*Vy))) / (sqrt(2.0*3.1415*Vy));
			Yp[i] = Yp[i] + YStep*y*Fy*EWy / WYB;
			y = y + YStep;
		}



	}
}

void Drift()//drift results from sampling some number N of the parents who survive selection. The mean of these parents after drift is the same as before drift, and the variance in their mean depends on the number sampled, N, and the phenotypic variance, Vx
{

	for (i = 0; i < Npops; i++)
	{
		Xp[i] = Xp[i] + Norm_Rand(0.0,sqrt(Vx/(1.0*Nx)));
		Yp[i] = Yp[i] + Norm_Rand(0.0,sqrt(Vy/(1.0*Ny)));
	}

}

void Mate()
{

	for (i = 0; i < Npops; i++)
	{
		X[i] = hX*Xp[i] +(1-hX)*X[i];
		Y[i] = hY*Yp[i] +(1-hY)*Y[i];
	}

}

void TraitOutput()
{
	out_traits << G << ",";
	for (i = 0; i < Npops; i++)
	{
		out_traits << X[i] << "," << Y[i] << ",";
	}
	out_traits << "\n";

}
void MetaOutput()
{
	out_meta <<Runs<<","<<G<<","<< MetaMeanX << "," << MetaMeanY << "," << MetaVarX << "," << MetaVarY << "," << MetaCov <<"\n";
}
void MetaStats()
{
	double XB, YB, XB2, YB2, XBYB,Sum1,Sum2,Sum3,Sum4,Sum5;


	//First do means
	XB = 0;
	YB = 0;
	for (i = 0; i < Npops; i++)
	{
		XB = XB + X[i];
		YB = YB + Y[i];
	}

	MetaMeanX = XB / (1.0*Npops);
	MetaMeanY = YB / (1.0*Npops);

	//Then higher moments
	XB2 = 0;
	YB2 = 0;
	XBYB = 0;

	for (i = 0; i < Npops; i++)
	{
		XB2 = XB2 + (X[i]- MetaMeanX)*(X[i] - MetaMeanX);
		YB2 = YB2 + (Y[i] - MetaMeanY)*(Y[i] - MetaMeanY);
		XBYB = XBYB + (X[i] - MetaMeanX)*(Y[i] - MetaMeanY);
	}

	MetaVarX = XB2 / (1.0*Npops-1.0);
	MetaVarY = YB2 / (1.0*Npops-1.0);
	MetaCov = XBYB / (1.0*Npops-1.0);
	MetaSDX = sqrt(MetaVarX);
	MetaSDY = sqrt(MetaVarY);
	MetaRo= MetaCov / (sqrt(MetaVarX)*sqrt(MetaVarY));


	StoreMetaMeanX[G] = MetaMeanX;
	StoreMetaMeanY[G] = MetaMeanY;
	StoreMetaVarX[G] = MetaVarX;
	StoreMetaVarY[G] = MetaVarY;
	StoreMetaCov[G] = MetaCov;


	StoreMetaSDX[G] = sqrt(MetaVarX);
	StoreMetaSDY[G] = sqrt(MetaVarY);
	StoreMetaRo[G] = MetaCov/(sqrt(MetaVarX)*sqrt(MetaVarY));



	if (G>0&&G%100 == 0)//Check to see if an approximate equilibrium as been reached every 100 generations
	{

		Sum1 = 0;
		Sum2 = 0;
		Sum3 = 0;
		Sum4 = 0;
		Sum5 = 0;
		for (i = 0; i < 10; i++)
		{
			Sum1 = Sum1 + (fabs(StoreMetaMeanX[G - i] - StoreMetaMeanX[G - i - 1])/ StoreMetaMeanX[G - i - 1]) /10.0;
			Sum2 = Sum2 + (fabs(StoreMetaMeanY[G - i] - StoreMetaMeanY[G - i - 1])/ StoreMetaMeanY[G - i - 1]) / 10.0;
			Sum3 = Sum3 + (fabs(StoreMetaSDX[G - i] - StoreMetaSDX[G - i - 1])/ StoreMetaSDX[G - i - 1]) / 10.0;
			Sum4 = Sum4 + (fabs(StoreMetaSDY[G - i] - StoreMetaSDY[G - i - 1])/ StoreMetaSDY[G - i - 1]) / 10.0;
			Sum5 = Sum5 + fabs((fabs(StoreMetaRo[G - i] - StoreMetaRo[G - i - 1])/ StoreMetaRo[G - i - 1])) / 10.0;
		}

	

		if (Sum1 < 0.01&&Sum2<0.01&&Sum3<0.05&&Sum4<0.05&&Sum5<0.05)
		{
			G = 10000;
			Equilibrated = 1;
		}
	}
}

void TestSummary()
{
	int Match;
	double DifMetaMeanX, DifMetaMeanY, DifMetaSDX, DifMetaSDY, DifMetaRo;


	Match = 0;
	DifMetaMeanX = fabs(MetaMeanX - DataMeanX);
	DifMetaMeanY = fabs(MetaMeanY - DataMeanY);
	DifMetaSDX = fabs(MetaSDX - DataSDX);
	DifMetaSDY = fabs(MetaSDY - DataSDY);
	DifMetaRo = fabs(MetaRo - DataRo);

	//output diagnostics
	cout << MetaMeanX - DataMeanX << "," << MetaMeanY - DataMeanY << "," << MetaSDX - DataSDX << "," << MetaSDY - DataSDY << "," << MetaRo - DataRo << "\n";

	if (DifMetaMeanX > ThreshEx)
	{
		Match++;
	}
	if (DifMetaMeanY > ThreshEy)
	{
		Match++;
	}
	if (DifMetaSDX > ThreshSDx)
	{
		Match++;
	}
	if (DifMetaSDY > ThreshSDy)
	{
		Match++;
	}
	if (DifMetaRo > ThreshRo)
	{
		Match++;
	}




	if (Match == 0)
	{
		
		if (Equilibrated = 1)//Only report to posterior if an equilibrium was reached
		{
			Hits++;
			EAXStore[Hits] = EAx;
			EAYStore[Hits] = EAy;
			VAXStore[Hits] = VAx;
			VAYStore[Hits] = VAy;
			out_posterior << Runs << "," << EAx << "," << EAy << "," << VAx << "," << VAy << "," << gx << "," << gy << "," << hX << "," << hY << "," << mx << "," << my << "," << Vx << "," << Vy << "," << Nx << "," << Ny << "," << ETx << "," << ETy << "," << VTx << "," << VTy <<"\n";
		}
	}

}

void ProcessPosterior() 
{
	int Wrong,SuccessE,SuccessV,Hi,Lo,Med,Bins,Included,Max,EAXBin[1000],EAYBin[1000], VAXBin[1000], VAYBin[1000], i;
	double Temp1, Temp2;
	double S1Lo, S1Hi, BinWidthEAX, BinWidthEAY, BinWidthVAX, BinWidthVAY;
	int LoEAX, HiEAX, LoEAY, HiEAY, TempLoEAX, TempHiEAX, TempLoEAY, TempHiEAY, ModeEAX, ModeEAY, LoVAX, HiVAX, LoVAY, HiVAY, TempLoVAX, TempHiVAX, TempLoVAY, TempHiVAY, ModeVAX, ModeVAY,ToInclude,IntStart, Shortest,Length;
	
	Hi = int(0.975*PosteriorSize);
	Lo = int(0.025*PosteriorSize);
	Med = int(0.500*PosteriorSize);

	//SORT AND PROCESS EXPECTED VALUES OF COEVOLUTION
	//then sort X
	Wrong = 1;
	while (Wrong > 0)
	{
		Wrong = 0;
		for (i = 0; i < Hits; i++)
		{
			if (EAXStore[i] > EAXStore[i + 1])
			{
				Temp1 = EAXStore[i];
				Temp2 = EAXStore[i + 1];
				EAXStore[i] = Temp2;
				EAXStore[i + 1] = Temp1;
				Wrong++;
			}
		}
	}

	//then sort Y
	Wrong = 1;
	while (Wrong > 0)
	{
		Wrong = 0;
		for (i = 0; i < Hits; i++)
		{
			if (EAYStore[i] > EAYStore[i + 1])
			{
				Temp1 = EAYStore[i];
				Temp2 = EAYStore[i + 1];
				EAYStore[i] = Temp2;
				EAYStore[i + 1] = Temp1;
				Wrong++;
			}
		}
	}

	SuccessE = 0;
	//then ask whether the true value lies between the low 2.5% and the high 2.5%
	if (DataEAx >= EAXStore[Lo] && DataEAx <= EAXStore[Hi]&& DataEAy >= EAYStore[Lo] && DataEAy <= EAYStore[Hi])
	{
		SuccessE=1;
	}
	//DONE PROCESSING POSTERIORS FOR EXPECTED STRENGTH OF COEVOLUTION

	//NOW PROCESS POSTERIORS FOR THE VARIANCE IN THE STRENGTH OF COEVOLUTION
	//then sort X
	Wrong = 1;
	while (Wrong > 0)
	{
		Wrong = 0;
		for (i = 0; i < Hits; i++)
		{
			if (VAXStore[i] > VAXStore[i + 1])
			{
				Temp1 = VAXStore[i];
				Temp2 = VAXStore[i + 1];
				VAXStore[i] = Temp2;
				VAXStore[i + 1] = Temp1;
				Wrong++;
			}
		}
	}

	//then sort Y
	Wrong = 1;
	while (Wrong > 0)
	{
		Wrong = 0;
		for (i = 0; i < Hits; i++)
		{
			if (VAYStore[i] > VAYStore[i + 1])
			{
				Temp1 = VAYStore[i];
				Temp2 = VAYStore[i + 1];
				VAYStore[i] = Temp2;
				VAYStore[i + 1] = Temp1;
				Wrong++;
			}
		}
	}

	SuccessV = 0;
	//then ask whether the true value lies between the low 2.5% and the high 2.5%
	if (DataVAx >= VAXStore[Lo] && DataVAx <= VAXStore[Hi] && DataVAy >= VAYStore[Lo] && DataVAy <= VAYStore[Hi])
	{
		SuccessV = 1;
	}
	//DONE PROCESSING POSTERIORS FOR VARIANCE IN THE STRENGTH OF COEVOLUTION
	


	//*****************************************************NOW APPROACH THE PROBLEM USING THE MPD AND MODE***************************************************************************
	// TO BE STRICTLY AN HPD, THE BELOW ASSUMES THE POSTERIOR IS UNIMODAL. I MAKE THIS ASSUMPTION BECAUSE OTHERWISE WE HAVE MULTI-INTERVAL 
	// HPD's DUE TO DISCRETE NATURE OF POSTERIOR WITH SMALL POSTERIOR SIZE EVEN WHEN TRULY UNIMODAL
	// Updated May 10, 2018
	//*******************************************************************************************************************************************************************************
	
	//SET THE NUMBER OF BINS
	Bins = 20;

	//start with SPECIES X EXPECTED VALUE OF COEVOLUTION

	//zero the bins
	for (j = 0; j < Bins; j++)
	{
		EAXBin[j] = 0;
	}
	//Done zeroing bins
	BinWidthEAX = (EAXStore[Hits-1]- EAXStore[0])/(1.0*Bins);
	for (i = 0; i<Hits; i++)
	{

	
		S1Lo = 0;
		S1Hi = S1Lo + BinWidthEAX;
		for (j = 0; j < Bins; j++)
		{

			if (EAXStore[i] > S1Lo&&EAXStore[i] <= S1Hi)
			{
				EAXBin[j]++;
			}
			S1Lo = S1Lo + BinWidthEAX;
			S1Hi = S1Hi + BinWidthEAX;
		}


	}


	//Now find mode
	Max = 0;
	for (j = 0; j<Bins; j++)
	{
		if (EAXBin[j]>Max)
		{
			ModeEAX = j;
			Max = EAXBin[j];
		}
	}


	//Now find HPD REVISED MAY 10, 2018
	ToInclude = int(Hits*0.95);
	Shortest = Bins;//Set the starting minimum at the full length
	for (IntStart = 0; IntStart < Bins; IntStart++)//TRY ALL POSSIBLE STARTING POINTS FOR THE MPD INTERVAL
	{
		TempLoEAX = IntStart;
		TempHiEAX = IntStart;
		Included = 0;
		i = 0;
		while (Included < ToInclude&&TempHiEAX<Bins)//THE NUMBER IN HERE DEPENDS ON THE CONFIDENCE INTERVAL YOU WISH TO ACHIEVE AND THE NUMBER OF HITS IN THE POSTERIOR
		{
			Included = Included + EAXBin[IntStart + i];
			TempHiEAX = IntStart + i;
			i++;
		}

		//CHECK TO SEE IF NARROWER. IF SO, REPLACE.
		Length = TempHiEAX - TempLoEAX;
		if (Length < Shortest&&Included>=ToInclude)
		{
			LoEAX = TempLoEAX;
			HiEAX = TempHiEAX;
			Shortest = Length;
		}
	}
	//Done finding MPD


	//THEN DO SPECIES X VARIANCE IN VALUE OF COEVOLUTION

	//zero the bins
	for (j = 0; j < Bins; j++)
	{
		VAXBin[j] = 0;
	}
	//Done zeroing bins
	BinWidthVAX = (VAXStore[Hits - 1] - VAXStore[0]) / (1.0*Bins);
	for (i = 0; i<Hits; i++)
	{
	
		S1Lo = 0;
		S1Hi = S1Lo + BinWidthVAX;
		for (j = 0; j<Bins; j++)
		{

			if (VAXStore[i]>S1Lo&&VAXStore[i] <= S1Hi)
			{
				VAXBin[j]++;
			}
			S1Lo = S1Lo + BinWidthVAX;
			S1Hi = S1Hi + BinWidthVAX;
		}
	}
	//Now find mode
	Max = 0;
	for (j = 0; j<Bins; j++)
	{
		if (VAXBin[j]>Max)
		{
			ModeVAX = j;
			Max = VAXBin[j];
		}
	}

	//Now find HPD REVISED MAY 10, 2018
	ToInclude = int(Hits*0.95);
	Shortest = Bins;//Set the starting minimum at the full length
	for (IntStart = 0; IntStart < Bins; IntStart++)//TRY ALL POSSIBLE STARTING POINTS FOR THE MPD INTERVAL
	{
		TempLoVAX = IntStart;
		TempHiVAX = IntStart;
		Included = VAXBin[IntStart];
		i = 1;
		while (Included < ToInclude&&TempHiVAX<Bins)//THE NUMBER IN HERE DEPENDS ON THE CONFIDENCE INTERVAL YOU WISH TO ACHIEVE AND THE NUMBER OF HITS IN THE POSTERIOR
		{
			Included = Included + VAXBin[IntStart + i];
			TempHiVAX = IntStart + i;
			i++;
		}

		//CHECK TO SEE IF NARROWER. IF SO, REPLACE.
		Length = TempHiVAX - TempLoVAX;
		if (Length < Shortest&&Included >= ToInclude)
		{
			LoVAX = TempLoVAX;
			HiVAX = TempHiVAX;
			Shortest = Length;
		}
	}
	//Done finding MPD


	//THEN DO SPECIES Y EXPECTED VALUES OF COEVOLUTION
	//zero the bins
	for (j = 0; j < Bins; j++)
	{
		EAYBin[j] = 0;
	}
	//Done zeroing bins
	BinWidthEAY = (EAYStore[Hits - 1] - EAYStore[0]) / (1.0*Bins);
	for (i = 0; i<Hits; i++)
	{
	
		S1Lo = 0;
		S1Hi = S1Lo + BinWidthEAY;
		for (j = 0; j<Bins; j++)
		{

			if (EAYStore[i]>S1Lo&&EAYStore[i] <= S1Hi)
			{
				EAYBin[j]++;
			}
			S1Lo = S1Lo + BinWidthEAY;
			S1Hi = S1Hi + BinWidthEAY;
		}

	}
	//Now find mode
	Max = 0;
	for (j = 0; j<Bins; j++)
	{
		if (EAYBin[j]>Max)
		{
			ModeEAY = j;
			Max = EAYBin[j];
		}
	}
	//Now find HPD REVISED MAY 10, 2018
	ToInclude = int(Hits*0.95);
	Shortest = Bins;//Set the starting minimum at the full length
	for (IntStart = 0; IntStart < Bins; IntStart++)//TRY ALL POSSIBLE STARTING POINTS FOR THE MPD INTERVAL
	{
		TempLoEAY = IntStart;
		TempHiEAY = IntStart;
		Included = EAYBin[IntStart];
		i = 1;
		while (Included < ToInclude&&TempHiEAY<Bins)//THE NUMBER IN HERE DEPENDS ON THE CONFIDENCE INTERVAL YOU WISH TO ACHIEVE AND THE NUMBER OF HITS IN THE POSTERIOR
		{
			Included = Included + EAYBin[IntStart + i];
			TempHiEAY = IntStart + i;
			i++;
		}

		//CHECK TO SEE IF NARROWER. IF SO, REPLACE.
		Length = TempHiEAY - TempLoEAY;
		if (Length < Shortest&&Included >= ToInclude)
		{
			LoEAY = TempLoEAY;
			HiEAY = TempHiEAY;
			Shortest = Length;
		}
	}
	//Done finding MPD


	//AND FINALLY, SPECIES Y VARIANCE IN THE VALUE OF COEVOLUTION
	//zero the bins
	for (j = 0; j < Bins; j++)
	{
		VAYBin[j] = 0;
	}
	//Done zeroing bins
	BinWidthVAY = (VAYStore[Hits - 1] - VAYStore[0]) / (1.0*Bins);
	for (i = 0; i<Hits; i++)
	{
		
		S1Lo = 0;
		S1Hi = S1Lo + BinWidthVAY;
		for (j = 0; j<Bins; j++)
		{

			if (VAYStore[i]>S1Lo&&VAYStore[i] <= S1Hi)
			{
				VAYBin[j]++;
			}
			S1Lo = S1Lo + BinWidthVAY;
			S1Hi = S1Hi + BinWidthVAY;
		}
	}
	//Now find mode
	Max = 0;
	for (j = 0; j<Bins; j++)
	{
		if (VAYBin[j]>Max)
		{
			ModeVAY = j;
			Max = VAYBin[j];
		}
	}
	//Now find HPD REVISED MAY 10, 2018
	ToInclude = int(Hits*0.95);
	Shortest = Bins;//Set the starting minimum at the full length
	for (IntStart = 0; IntStart < Bins; IntStart++)//TRY ALL POSSIBLE STARTING POINTS FOR THE MPD INTERVAL
	{
		TempLoVAY = IntStart;
		TempHiVAY = IntStart;
		Included = VAYBin[IntStart];
		i = 1;
		while (Included < ToInclude&&TempHiVAY<Bins)//THE NUMBER IN HERE DEPENDS ON THE CONFIDENCE INTERVAL YOU WISH TO ACHIEVE AND THE NUMBER OF HITS IN THE POSTERIOR
		{
			Included = Included + VAYBin[IntStart + i];
			TempHiVAY = IntStart + i;
			i++;
		}

		//CHECK TO SEE IF NARROWER. IF SO, REPLACE.
		Length = TempHiVAY - TempLoVAY;
		if (Length < Shortest&&Included >= ToInclude)
		{
			LoVAY = TempLoVAY;
			HiVAY = TempHiVAY;
			Shortest = Length;
		}
	}
	//Done finding MPD

	SuccessE = 0;
	//then ask whether the true value lies between the low 2.5% and the high 2.5%
	if (DataEAx >= LoEAX*BinWidthEAX  && DataEAx <= HiEAX*BinWidthEAX  && DataEAy >= LoEAY*BinWidthEAY  && DataEAy <= HiEAY*BinWidthEAY)
	{
		SuccessE = 1;
	}

	SuccessV = 0;
	//then ask whether the true value lies between the low 2.5% and the high 2.5%
	if (DataVAx >= LoVAX*BinWidthVAX  && DataVAx <= HiVAX*BinWidthVAX  && DataVAy >= LoVAY*BinWidthVAY  && DataVAy <= HiVAY*BinWidthVAY)
	{
		SuccessV = 1;
	}
	

	//output summary stats
	if (Equilibrated = 1)//Only report to posterior if an equilibrium was reached
	{
	

		//output EAX Bin labels and bin counts
		out_Bins << "EAX (Bins):,";
		for (i = 0; i < Bins; i++)
		{
			out_Bins << i*BinWidthEAX << ",";
		}
		out_Bins << "\n";
		out_Bins << "EAX (Counts):,";
		for (i = 0; i < Bins; i++)
		{
			out_Bins << EAXBin[i] << ",";
		}
		out_Bins << "\n";

		//output EAY Bin labels and bin counts
		out_Bins << "EAY (Bins):,";
		for (i = 0; i < Bins; i++)
		{
			out_Bins << i*BinWidthEAY << ",";
		}
		out_Bins << "\n";
		out_Bins << "EAY (Counts):,";
		for (i = 0; i < 20; i++)
		{
			out_Bins << EAYBin[i] << ",";
		}
		out_Bins << "\n";

		//output VAX Bin labels and bin counts
		out_Bins << "VAX (Bins):,";
		for (i = 0; i <Bins; i++)
		{
			out_Bins << i*BinWidthVAX << ",";
		}
		out_Bins << "\n";
		out_Bins << "VAX (Counts):,";
		for (i = 0; i < Bins; i++)
		{
			out_Bins << VAXBin[i] << ",";
		}
		out_Bins << "\n";

		//output EAX Bin labels and bin counts
		out_Bins << "VAY (Bins):,";
		for (i = 0; i <Bins; i++)
		{
			out_Bins << i*BinWidthVAY << ",";
		}
		out_Bins << "\n";
		out_Bins << "VAY (Counts):,";
		for (i = 0; i < Bins; i++)
		{
			out_Bins << VAYBin[i] << ",";
		}
		out_Bins << "\n\n";
	}
}