#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>  


//ADD INTO THIS DECLINING SPECIES ABUNDNANCE AND SHIFTS IN EVENNESS (MORE ABUNDANT GETS MORE ABUNDNAT
//TEST TO SEE IF RATES OF SPECIES LOSS AND INTERACTION LOSS DEPEND ON COEVOLUTION VS. NO COEVOLUTION


time_t seconds;

//Created January 13, 2017


using namespace std;


char FileRoot[50], FigureData[50] = "Figs_", TraitData[50] = "Traits_", MetaData[50] = "Meta_", AbundanceData[50] = "Abundance_", SampledTraitData[50] = "SampledTraits_", AbundByHeritData[50] = "AbundByHerit_", BInteractionData[50] = "BInteractions_", NetworkSummary[50] = "NetworkSummary_", Parameters[50] = "Parameters_", TraitDegree[50] = "TraitDegree_", AnalyticalPredictions[50] = "AnalyticalPredictions_", SummaryStats[50] = "SummaryStats_", Phylogeny[50] = "Phylogeny_", K[50] = "K_", T[50] = "T_",EvoRate[50]="EvoRate_";
char fileinA[50], fileinP[50];
double EThetaA,EThetaP,VThetaA,VThetaP,EEpsA,EEpsP,VEpsA,VEpsP,EChiA,EChiP,VChiA,VChiP;
// read the below as the phenotype of [species][individual]
double ZA[100][30000],ZP[100][30000],ZAP[100][30000],ZPP[100][30000];
double WAbioPlant[100][30000],WAbioAnimal[100][30000], WBioPlant[100][30000], WBioAnimal[100][30000], WTotPlant[100][30000], WTotAnimal[100][30000], TrueMeanEpsA,TrueMeanEpsP;
double ThetaA[100],ThetaP[100],GammaA[100],GammaP[100],AlphaM,EpsA[100],EpsP[100],AlphaC,ChiA[100],ChiP[100],PBaseW,ABaseW,AnimalSampled[100],PlantSampled[100];
double SegVarAnimal,SegVarPlant;
double PMean[100],AMean[100],PVar[100],AVar[100];
double C,NODFA,NODFP;
double ComplementaritySampleA[10000000],ComplementaritySampleP[10000000];
int KA[100],KP[100], KAroot[100], KProot[100], NP[100],NA[100],NPP[100],NAP[100],NumIntObs,NPSamp,NASamp;
int IMAT[100][100],BIMAT[100][100], UniqueInteractions;
int PSpecies,ASpecies,PS,AS, Babies[100][10000];
int NInteract,Sampled,ModelType,TotalAnimals,TotalPlants,MAX;
int DegA[100],DegP[100];
int Runs,hRun,pRun,G,Gmax,PureNull,Large_Random;
double MeanEpsA,MeanEpsP,VarEpsA,VarEpsP,MeanKA,MeanKP,VarKA,VarKP,MeanGammaA,MeanGammaP,VarGammaA,VarGammaP,MeanThetaA,MeanThetaP,VarThetaA,VarThetaP,MeanChiA,MeanChiP,VarChiA,VarChiP,MeanSA,MeanSP;
double ENullC,ENullComplementarity,ENullNODFA,ENullNODFP,ENullVA,ENullVP,CorrectNODFA,CorrectNODFP,APercentGreater,APercentLesser,PPercentGreater,PPercentLesser,VAPercentile,VPPercentile,CPercentile,CorrectVA,CorrectVP,CorrectComp;
double PMetaMean,AMetaMean,PMetaVar,AMetaVar,PThetaMetaMean,AThetaMetaMean,PThetaMetaVar,AThetaMetaVar,PKMetaMean,PKMetaVar,AKMetaMean,AKMetaVar,PExpVar,AExpVar,AnSolPMetaMean,AnSolAMetaMean,AnSolAMetaVar,AnSolPMetaVar,SP,SA,GamSetA,GamSetP;
double ThetaVarA,ThetaVarP,AnSolCon;
double AnSolAVar,AnSolPVar,AnSolAMean,AnSolPMean,AnSolCor,IPComplementarity,IPVA,IPVP,IPMeanA,IPMeanP,IPCovAP,Efficiency,LNMean,LNVar;
double NullBIMAT[100][100],NullIMAT[100][100],NullComplementaritySampleA[10000000],NullComplementaritySampleP[10000000],NullC,NullComplementarity,NullNODFA,NullNODFP,NullDegA[100],NullDegP[100],NEA,NEP;
int z,zz,Test,AnimalReproTrials,PlantReproTrials,Encounters, SuccessfulEncounters;
double EvenA, EvenP, HA, HP,HInteract, HAmax, HPmax,PerPlantSuccesses,PerPlantEncounters;
int General[100][100];
double FIMAT[100][100],MeanGeneral,MeanADeg,MeanPDeg;
int TdivA[100][100], TdivP[100][100], TotalHabitatA, TotalHabitatP;

int BabySpecies[10000000], BabyMother[10000000], PlantSpecies[10000000], PlantIndividual[10000000], ExtantPlants, ExtantAnimals;

int test,ResCount,it,iii,AnimalWhack[200],PlantWhack[200];
int AExtinct[100], PExtinct[100],Import,EH;
int AE, PE,ObservableAnimals,ObservablePlants,HModel,DurModel;
double ATM, PTM,hARoot,hPRoot,hA[100],hP[100],PropLostA,PropLostP, ALossRate, PLossRate,AEChangeRate,PEChangeRate,AnimalChange[100],PlantChange[100];
int DestroyTime,AssemblyTime,coRun,Surrogate, DisturbType;
double PhyloThetaA[100], PhyloThetaP[100];
double AEvo[100], PEvo[100],ANullEvo[100],PNullEvo[100], AMetaMeanEvo, PMetaMeanEvo, AMetaMeanNullEvo, PMetaMeanNullEvo;


void DrawStaticPars();
void DrawDynamicPars();
void MakeIndividuals();
void Constraint();
void InterGuild();
void Mate();
void PopStats();
void Output();
void InitializeOutput();
//void OutputInteract();
void AnalyzeInteractions();
void NullNetwork();
void AnPred();
void SetK();
void SetTheta();
void BirthTheta();
void Immigrate();



ofstream out_traits;
ofstream out_meta;
ofstream out_abundance;
ofstream out_sampledtraits;
ofstream out_AbundByHerit;
ofstream out_Binteract;
ofstream out_network_summary;
ofstream out_pars;
ofstream out_trait_degree;
ofstream out_pred;
ofstream out_SummaryStats;
ofstream out_Phylogeny;
ofstream out_K; 
ofstream out_T;
ofstream out_Figs;
ofstream out_EvoRate;



int Big_Rand(int MAX)
{
	random_device rd;   // non-deterministic generator  
	mt19937 gen(rd());  // to seed mersenne twister.  
	uniform_int_distribution<> dist(0, 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);
}



int main()

{
	srand(time(NULL));

	cout << "Enter root file name\n";
	cin >> FileRoot;

	out_traits.open(strcat(TraitData, FileRoot));
	out_abundance.open(strcat(AbundanceData, FileRoot));
	out_sampledtraits.open(strcat(SampledTraitData, FileRoot));
	out_AbundByHerit.open(strcat(AbundByHeritData, FileRoot));
	out_Binteract.open(strcat(BInteractionData, FileRoot));
	out_network_summary.open(strcat(NetworkSummary, FileRoot));
	out_pars.open(strcat(Parameters, FileRoot));
	out_trait_degree.open(strcat(TraitDegree, FileRoot));
	out_pred.open(strcat(AnalyticalPredictions, FileRoot));
	out_meta.open(strcat(MetaData, FileRoot));
	out_SummaryStats.open(strcat(SummaryStats, FileRoot));
	out_Phylogeny.open(strcat(Phylogeny, FileRoot));
	out_K.open(strcat(K, FileRoot));
	out_T.open(strcat(T, FileRoot));
	out_Figs.open(strcat(FigureData, FileRoot));
	out_EvoRate.open(strcat(EvoRate, FileRoot));




	//cout << "Matching (0) Escalation (1) or ?\n";
	//cin >> ModelType;
	ModelType = 0;
	//cout << "Enter the number of Animal species\n";
	//cin >> ASpecies;
	ASpecies = 80;
	//cout << "Enter the number of Plant species\n";
	//cin >> PSpecies;
	PSpecies = 80;
	cout << "Set Alpha\n";
	cin >> AlphaM;
	cout << "Habitat destruction (1) or climate change (2)?\n";
	cin >> DisturbType;
	cout<<"Heritability Model (0), (1), or (2)\n";
	cin >> HModel;
	cout << "Coevolutionary Duration (0), (1), or (2)\n";
	cin >> DurModel;




	out_SummaryStats << "DisturbanceType,EvoHist,PRun,HRun,AssemblyTime,hA,hP,alpha,gA,gP,AssemblyTime,DestroyTime,PropLostA,PropLostP,ExtantPlants(BD),ExtantAnimals(BD),UniqueInteractions(BD),HA(BD),HP(BD),HI(BD),Efficiency(BD),ADegree(BD),PDegree(BD),,ExtantPlants(AD),ExtantAnimals(AD),UniqueInteractions(AD),HA(AD),HP(AD),HI(AD),Efficiency(AD),ADegree(AD),PDegree(AD)\n";
	out_EvoRate << "DisturbanceType,EvoHist,PRun,HRun,AssemblyTime,hA,hP,alpha,gA,gP,AssemblyTime,DestroyTime,PropLostA,PropLostP,G,AnimalEvo,PlantEvo,AnimalNullEvo,PlantNullEvo\n";
	InitializeOutput();
	//for (EH = 0; EH < 5; EH++)//run five evolutionary histories for each parameter combination
	for (EH = 0; EH < 10; EH++)//run 10 evolutionary histories for each parameter combination
	{
		BirthTheta(); //make a birth tree and use it to simulate theta


		/**
		cout<<"Enter EpsA\n";
		cin>>TrueMeanEpsA;
		cout<<"Enter EpsP\n";
		cin>>TrueMeanEpsP;*/
		TrueMeanEpsA = 5.0;//expected number of babies per successful interaction
		TrueMeanEpsP = 5.0;//expected number of babies per successful interaction
		//cout << "Enter Alpha\n";
		//cin >> AlphaM;

		//another way to disturb the world would be to reduce the expected number of plants an animal encounters in proportion to plant density...

		//AlphaM=5.0;


		InitializeOutput();


		for (AS = 0; AS < ASpecies; AS++)
		{
			AExtinct[AS] = -1;
		}
		for (PS = 0; PS < PSpecies; PS++)
		{
			PExtinct[PS] = -1;
		}

		DrawStaticPars();
		DrawDynamicPars();
		ResCount = 0;

		SetK();



		
		//for (pRun = 0; pRun <= 3; pRun++)//run over different levels of habitat destruction or environmental change
		for (pRun = 0; pRun <= 10; pRun++)//run over different levels of habitat destruction or environmental change Feb 16, 2018
		//for (pRun = 0; pRun < 1; pRun++)//run over different levels of habitat destruction or environmental change August 1, 2017
		{
			//for (hRun = 0; hRun <= 3; hRun++)//run over different heritabilties
			//for (hRun = 0; hRun <= 1; hRun++)//run over different heritabilties
			for (hRun = 0; hRun < 1; hRun++)//heritability is set manually on entry
			{
				//for (coRun = 0; coRun <= 2; coRun++)//run over different durations of coevolution
				for (coRun = 0; coRun <1; coRun++)//Set manually on entry Feb 16, 2018
				//for (coRun = 0; coRun <= 10; coRun++)//run over different durations of coevolution August 1, 2017
				{
					//reset thetas to their ancestral values
					for (PS = 0; PS < PSpecies; PS++)
					{
						ThetaP[PS] = PhyloThetaP[PS];
					}
					for (AS = 0; AS < ASpecies; AS++)
					{
						ThetaA[AS] = PhyloThetaA[AS];
					}


					SegVarAnimal = 0.01;
					SegVarPlant = 0.01;

					
					if (DurModel == 0)
					{
						AssemblyTime = 10; //Feb 16, 2018
					}
					if (DurModel == 1)
					{
						AssemblyTime = 100; //Feb 16, 2018
					}
					if (DurModel == 2)
					{
						AssemblyTime = 300; //Feb 16, 2018
					}
					//AssemblyTime = 0 + coRun * 20;//August 1, 2017
					DestroyTime = 50;
					Gmax = AssemblyTime + DestroyTime;

					//Set habitat loss target as a % of total habitat... 
					if (DisturbType == 1)//Habitat destruction
					{
						//PropLostA = 0.2 + 0.2*pRun;
						//PropLostP = 0.2 + 0.2*pRun;
						PropLostA = 0.4 + 0.05*pRun;//Feb 16, 2018
						PropLostP = 0.4 + 0.05*pRun;//Feb 16, 2018
						//PropLostA = 0.6;//August 1, 2017
						//PropLostP = 0.6;//August 1, 2017
						ALossRate = PropLostA / (1.0*DestroyTime);
						PLossRate = PropLostP / (1.0*DestroyTime);

						for (PS = 0; PS < PSpecies; PS++)
						{
							PlantWhack[PS] = PLossRate*(0.5 + 1.0*rand() / (1.0*RAND_MAX))*KProot[PS];
						}
						for (AS = 0; AS < ASpecies; AS++)
						{
							AnimalWhack[AS] = ALossRate*(0.5 + 1.0*rand() / (1.0*RAND_MAX))*KAroot[AS];
						}
					}

					if (DisturbType == 2)//Environmental change
					{
						//PropLostA = 1.0*sqrt(SegVarAnimal) + 2.0*sqrt(SegVarAnimal)*pRun;//I call this prop lost for economy of variables but it actually meausures total environmental change
						//PropLostP = 1.0*sqrt(SegVarPlant) + 2.0*sqrt(SegVarPlant)*pRun;//I call this prop lost for economy of variables but it actually meausures total environmental change
						PropLostA = 1.0*sqrt(SegVarAnimal) + 1.0*sqrt(SegVarAnimal)*pRun;//I call this prop lost for economy of variables but it actually meausures total environmental change
						PropLostP = 1.0*sqrt(SegVarPlant) + 1.0*sqrt(SegVarPlant)*pRun;//I call this prop lost for economy of variables but it actually meausures total environmental change
						AEChangeRate = PropLostA / (1.0*DestroyTime); //This is script C in the paper
						PEChangeRate = PropLostP / (1.0*DestroyTime);

						for (PS = 0; PS < PSpecies; PS++)
						{
							PlantChange[PS] = PEChangeRate*(0.5 + 1.0*rand() / (1.0*RAND_MAX)); //this is sc in the paper

						}
						for (AS = 0; AS < ASpecies; AS++)
						{
							AnimalChange[AS] = AEChangeRate*(0.5 + 1.0*rand() / (1.0*RAND_MAX));
						}
					}


					//Set carrying capacity at the start of each run based on predefined root values
					for (PS = 0; PS < PSpecies; PS++)
					{
						KP[PS] = KProot[PS];
					}
					for (AS = 0; AS < ASpecies; AS++)
					{
						KA[AS] = KAroot[AS];
					}




					MakeIndividuals();

					//Set heritability
					//hARoot = hRun / (4.0);
					//hPRoot = hRun / (4.0);
					//hARoot = hRun / (2.0);
					//hPRoot = hRun / (2.0);
				//	hARoot = hRun / (4.0/3.0);
					//hPRoot = hRun / (4.0/3.0);
					//hARoot =0.75;//August 1, 2017
				//	hPRoot = 0.75;//August 1, 2017
					if (HModel == 0)//No heritability at all
					{
						hARoot = 0.0;//Feb 16, 2018
						hPRoot = 0.0;//Feb 16, 2018
					}
					if (HModel == 1)//Unilateral evolution
					{
						hARoot = 0.0;//Feb 16, 2018
						hPRoot = 0.75;//Feb 16, 2018
					}
					if (HModel == 2)//Coevolution
					{
						hARoot = 0.75;//Feb 16, 2018
						hPRoot = 0.75;//Feb 16, 2018
					}


					//Done setting root values of heritbility
					//Set heritabilities for individual species
					for (PS = 0; PS < PSpecies; PS++)
					{
						hP[PS] = hPRoot + 0.2*rand() / (1.0*RAND_MAX) - 0.1;
						if (hP[PS] < 0)
						{
							hP[PS] = 0;
						}
						if (hP[PS] > 1)
						{
							hP[PS] = 1;
						}
						if (hPRoot == 0)
						{
							hP[PS] = 0;
						}
						if (hPRoot == 1)
						{
							hP[PS] = 1;
						}
					}
					for (AS = 0; AS < ASpecies; AS++)
					{
						hA[AS] = hARoot + 0.2*rand() / (1.0*RAND_MAX) - 0.1;
						if (hA[AS] < 0)
						{
							hA[AS] = 0;
						}
						if (hA[AS] > 1)
						{
							hA[AS] = 1;
						}
						if (hARoot == 0)
						{
							hA[AS] = 0;
						}
						if (hARoot == 1)
						{
							hA[AS] = 1;
						}
					}


					for (G = 0; G <= Gmax; G++)
					{
						//DRAW SOME SPECIES AT RANDOM TO COLONIZE and add 100 individuals to their numbers (NEED TO SET ALL ABUNDANCES to ZERO TO START AND ALL ASPECIES AND PSPECIES CALLS TO RUN THROUGH ALL POSSIBLE SPECIES)


						if (G == 0)//must run interactions prior to gen zero output 
						{
							//Immigrate();
							InterGuild();
							PopStats();
						}

						if (ResCount == 0 || ResCount == 1)
						{
							ResCount = 0;
							AnalyzeInteractions();
							Output();
						}

						if (G == AssemblyTime)
						{
							out_SummaryStats << DisturbType << "," <<EH<<","<< pRun << "," << hRun << "," << coRun << "," << hARoot << "," << hPRoot << "," << AlphaM << "," << GamSetA << "," << GamSetP << "," << AssemblyTime << "," << DestroyTime << "," << PropLostP << "," << PropLostA << "," << ExtantPlants << "," << ExtantAnimals << "," << UniqueInteractions << "," << HA << "," << HP << "," << HInteract << "," << Efficiency << "," << MeanADeg << "," << MeanPDeg << ",,";
						}

						if (G > AssemblyTime)
						{
							if (DisturbType == 1)//Habitat destruction
							{
								//Impose the disturbance
								for (PS = 0; PS < PSpecies; PS++)
								{
									if (PlantWhack[PS] > 0)
									{
										KP[PS] = KP[PS] - rand() % (2 * PlantWhack[PS]);
									}
								}

								//The error is in the for loop ABOVE
								for (AS = 0; AS < ASpecies; AS++)
								{
									if (AnimalWhack[AS] > 0)
									{
										KA[AS] = KA[AS] - rand() % (2 * AnimalWhack[AS]);
									}
								}
							}

							if (DisturbType == 2)//Environmental change
							{
								//Impose the disturbance
								for (PS = 0; PS < PSpecies; PS++)
								{
									ThetaP[PS] = ThetaP[PS] - PlantChange[PS] * (0.9 + 0.2*rand() / (1.0*RAND_MAX));//set to decrease theta to mimic climate change

								}

								
								for (AS = 0; AS < ASpecies; AS++)
								{
									ThetaA[AS] = ThetaA[AS] - AnimalChange[AS] * (0.9 + 0.2*rand() / (1.0*RAND_MAX));//set to decrease theta to mimic climate change
								}
							}

						}


						//Immigrate();
						Constraint();
						InterGuild();
						Mate();
						PopStats();

						out_EvoRate << DisturbType << "," << EH << "," << pRun << "," << hRun << "," << coRun << "," << hARoot << "," << hPRoot << "," << AlphaM << "," << GamSetA << "," << GamSetP << "," << AssemblyTime << "," << DestroyTime << "," << PropLostP << "," << PropLostA << "," << G<<","<<AMetaMeanEvo<<","<<PMetaMeanEvo<< "," << AMetaMeanNullEvo << "," << PMetaMeanNullEvo << "\n";

						//	PopStats();

						cout <<EH<<","<<pRun << "," << hRun << "," << coRun << "," << G << "\n";
						ResCount++;

						if (G == Gmax)
						{
							out_SummaryStats << ExtantPlants << "," << ExtantAnimals << "," << UniqueInteractions << "," << HA << "," << HP << "," << HInteract << "," << Efficiency << "," << MeanADeg << "," << MeanPDeg << "\n";

							out_AbundByHerit << "Plant Species,";
							for (PS = 0; PS < PSpecies; PS++)
							{
								out_AbundByHerit << PS << ",";
							}
							out_AbundByHerit << "\n";
							out_AbundByHerit << "Heritibility,";
							for (PS = 0; PS < PSpecies; PS++)
							{
								out_AbundByHerit << hP[PS] << ",";
							}
							out_AbundByHerit << "\n";
							out_AbundByHerit << "Abundance,";
							for (PS = 0; PS < PSpecies; PS++)
							{
								out_AbundByHerit << NP[PS] << ",";
							}
							out_AbundByHerit << "\n";

							out_AbundByHerit << "Animal Species,";
							for (AS = 0; AS < ASpecies; AS++)
							{
								out_AbundByHerit << AS << ",";
							}
							out_AbundByHerit << "\n";
							out_AbundByHerit << "Heritibility,";
							for (AS = 0; AS < ASpecies; AS++)
							{
								out_AbundByHerit << hA[AS] << ",";
							}
							out_AbundByHerit << "\n";
							out_AbundByHerit << "Abundance,";
							for (AS = 0; AS < ASpecies; AS++)
							{
								out_AbundByHerit << NA[AS] << ",";
							}
							out_AbundByHerit << "\n\n";

						}

						out_SummaryStats.flush();
						out_AbundByHerit.flush();

					}//End of main generational Loop

				}//end of coevolve time loop
			}//end of heritability runs
		}//end of proportion disturbed runs
	}//end of evolutionary history runs



	
	return 0;
}


void DrawStaticPars()
{
	double SumGammaP,SumGammaPP;
	double SumGammaA,SumGammaAA;	

	LNMean=6.5;
	LNVar=0.3;

	ThetaVarA=2.0;
	ThetaVarP=2.0;

	//SegVarAnimal = 0.05;// 0.05;// 1.0;
	//SegVarPlant = 0.05;// 0.05;// 1.0;

	GamSetA = 1.0 + 2.0*rand() / (1.0*RAND_MAX);
	GamSetP = 1.0 + 2.0*rand() / (1.0*RAND_MAX);

	//Sampled=200+400*(rand()%3);
	Sampled=10.0*ASpecies*PSpecies;
	
	AlphaC=0.0;

	ABaseW=1.0;
	PBaseW=1.0;

	PlantReproTrials = 1;
	AnimalReproTrials = 1;

	SumGammaP=0;
	SumGammaPP=0;

	for(PS=0;PS<PSpecies;PS++)
	{
		//Gamma is assumed fixed across species and is drawn at random at the guild level for every run
		GammaP[PS]=GamSetP;
		SumGammaP=SumGammaP+GammaP[PS];
		SumGammaPP=SumGammaPP+GammaP[PS]*GammaP[PS];
	}

	MeanGammaP=SumGammaP/(1.0*PSpecies);
	VarGammaP=SumGammaPP/(1.0*PSpecies)-MeanGammaP*MeanGammaP;
	
	SumGammaA=0;
	SumGammaAA=0;

	for(AS=0;AS<ASpecies;AS++)
	{
		//Gamma is assumed fixed across species and is drawn at random at the guild level for every run
		GammaA[AS]=GamSetA;
		SumGammaA=SumGammaA+GammaA[AS];
		SumGammaAA=SumGammaAA+GammaA[AS]*GammaA[AS];
	}

	MeanGammaA=SumGammaA/(1.0*ASpecies);
	VarGammaA=SumGammaAA/(1.0*ASpecies)-MeanGammaA*MeanGammaA;

}

void DrawDynamicPars()
{
	double SumThetaP,SumThetaPP,SumGammaP,SumGammaPP,SumChiP,SumChiPP,SumKP,SumKPP,SumEpsP,SumEpsPP,SumSP,SumSA;
	double SumThetaA,SumThetaAA,SumGammaA,SumGammaAA,SumChiA,SumChiAA,SumKA,SumKAA,SumEpsA,SumEpsAA;
	double u,v,r2,scale,Noise;
	double LNMean,LNVar;
	int SeedA,SeedP;

			
	SumEpsP=0;
	SumEpsPP=0;
	
	for(PS=0;PS<PSpecies;PS++)
	{
		//Draw eps from a normal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 *rand()/(1.0*RAND_MAX);
			v = -1 + 2 *rand()/(1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		}
		while (r2 > 1.0 || r2 == 0);
		scale = sqrt (-2.0 * log (r2) / r2);
		Noise = sqrt(0.5) * u * scale;
		EpsP[PS]=TrueMeanEpsP+Noise;
		if(TrueMeanEpsP==0)
		{
			EpsP[PS]=0;
		}
		//Eliminate negative values as they are biologically non-sensical
		if(EpsP[PS]<0)
		{
			EpsP[PS]=0;
		}

		SumEpsP=SumEpsP+EpsP[PS];
		SumEpsPP=SumEpsPP+EpsP[PS]*EpsP[PS];		
	}

	
	MeanEpsP=SumEpsP/(1.0*PSpecies);
	VarEpsP=SumEpsPP/(1.0*PSpecies)-MeanEpsP*MeanEpsP;

	//CALCULATE ANALYTICAL PARAMETER VALUES
	if(ModelType==0)// matching
	{
		SP=(2.0*AlphaM*MeanEpsP)/(1.0+MeanEpsP);// This is only for one kind of model dipshit!
	}
	if(ModelType==1)// escalation
	{
		SP=(AlphaM*MeanEpsP)/(2.0*(2.0+MeanEpsP));// This is only for one kind of model dipshit!
	}

	SumEpsA=0;
	SumEpsAA=0;
	
	for(AS=0;AS<ASpecies;AS++)
	{
		//Draw Eps from a normal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 *rand()/(1.0*RAND_MAX);
			v = -1 + 2 *rand()/(1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		}
		while (r2 > 1.0 || r2 == 0);
		scale = sqrt (-2.0 * log (r2) / r2);
		Noise = sqrt(0.5) * u * scale;//I changed this from .05 1/28/2011 upon realizing that it was too much variance (overwhelming the mean)
		EpsA[AS]=TrueMeanEpsA+Noise;
		//Truncate any values below zero as that is biologically non-sensical
		if(TrueMeanEpsA==0)
		{
			EpsA[AS]=0;
		}
		if(EpsA[AS]<0)
		{
			EpsA[AS]=0;
		}

		SumEpsA=SumEpsA+EpsA[AS];
		SumEpsAA=SumEpsAA+EpsA[AS]*EpsA[AS];
	
		
	}

	
	MeanEpsA=SumEpsA/(1.0*ASpecies);
	VarEpsA=SumEpsAA/(1.0*ASpecies)-MeanEpsA*MeanEpsA;

	//CALCULATE ANALYTICAL PARAMETER VALUES
	if(ModelType==0)// matching
	{
		SA=(2.0*AlphaM*MeanEpsA)/(1.0+MeanEpsA);// This is only for one kind of model dipshit!
	}
	if(ModelType==1)// escalation
	{
		SA=(AlphaM*MeanEpsA)/(2.0*(2.0+MeanEpsA));
	}





	//Now output the detailed parameters for each run to a special parameter file
	out_pars<<"Run #"<<pRun<<"."<<hRun<<",KA,";
	for(AS=0;AS<ASpecies;AS++)
	{
		out_pars<<KAroot[AS]<<",";
	}
	out_pars<<"\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",ThetaA,";
	for(AS=0;AS<ASpecies;AS++)
	{
		out_pars<<ThetaA[AS]<<",";
	}
	out_pars<<"\n";
	out_pars << "Run #" << pRun << "." << hRun << ",hA,";
	for (AS = 0; AS<ASpecies; AS++)
	{
		out_pars << hA[AS] << ",";
	}
	out_pars << "\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",EpsA,";
	for(AS=0;AS<ASpecies;AS++)
	{
		out_pars<<EpsA[AS]<<",";
	}
	out_pars<<"\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",GamA,";
	for(AS=0;AS<ASpecies;AS++)
	{
		out_pars<<GammaA[AS]<<",";
	}
	out_pars<<"\n";


	out_pars<<"Run #"<< pRun << "." << hRun <<",KP,";
	for(PS=0;PS<PSpecies;PS++)
	{
		out_pars<<KProot[PS]<<",";
	}
	out_pars<<"\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",ThetaP,";
	for(PS=0;PS<PSpecies;PS++)
	{
		out_pars<<ThetaP[PS]<<",";
	}
	out_pars<<"\n";
	out_pars << "Run #" << pRun << "." << hRun << ",hP,";
	for (PS = 0; PS<PSpecies; PS++)
	{
		out_pars << hP[PS] << ",";
	}
	out_pars << "\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",EpsP,";
	for(PS=0;PS<PSpecies;PS++)
	{
		out_pars<<EpsP[PS]<<",";
	}
	out_pars<<"\n";
	out_pars<<"Run #"<< pRun << "." << hRun <<",GamP,";
	for(PS=0;PS<PSpecies;PS++)
	{
		out_pars<<GammaP[PS]<<",";
	}
	out_pars<<"\n";

	out_pars<<"\n";	

	out_pars.flush();
}

void MakeIndividuals()
{
	int i;
	double u, v, r2, scale, Noise;
	
	for (PS = 0; PS < PSpecies; PS++)
	{
		NP[PS] = KP[PS];

		for (i = 0; i < KP[PS]; i++)
		{
			//Draw starting variance from a normal
			do
			{
				// choose x,y in uniform square (-1,-1) to (+1,+1) 
				u = -1 + 2 * rand() / (1.0*RAND_MAX);
				v = -1 + 2 * rand() / (1.0*RAND_MAX);
				// see if it is in the unit circle 
				r2 = u * u + v * v;
			} while (r2 > 1.0 || r2 == 0);
			scale = sqrt(-2.0 * log(r2) / r2);
			Noise = sqrt(SegVarPlant) * u * scale;
			ZP[PS][i] = ThetaP[PS] + Noise;
		}
	}
	
	for (AS = 0; AS < ASpecies; AS++)
	{
		NA[AS] = KA[AS];
	
		for (i = 0; i < KA[AS]; i++)
		{
			//Draw starting varuance from a normal
			do
			{
				// choose x,y in uniform square (-1,-1) to (+1,+1) 
				u = -1 + 2 * rand() / (1.0*RAND_MAX);
				v = -1 + 2 * rand() / (1.0*RAND_MAX);
				// see if it is in the unit circle 
				r2 = u * u + v * v;
			} while (r2 > 1.0 || r2 == 0);
			scale = sqrt(-2.0 * log(r2) / r2);
			Noise = sqrt(SegVarAnimal) * u * scale;
			ZA[AS][i] = ThetaA[AS] + Noise;
		}
	}

}

void Constraint()
{
	//first, lets implement stabilizing selection
	int k;
		
	//First do plants
	for(PS=0;PS<PSpecies;PS++)
	{
		for(k=0;k<NP[PS];k++)
		{	
			WAbioPlant[PS][k]=exp(-GammaP[PS]*pow(ZP[PS][k]-ThetaP[PS],2.0));


		}
	}
	//Animal
	for(AS=0;AS<ASpecies;AS++)
	{
		for(k=0;k<NA[AS];k++)
		{
			WAbioAnimal[AS][k]=exp(-GammaA[AS]*pow(ZA[AS][k]-ThetaA[AS],2.0));
		
		}
	}
	

}

void InterGuild() 
{
	//TRY THIS INSTEAD: MAKE AN UNFOLDED LIST OF ALL PLANTS. THEN CHOOSE EACH POLLINATOR AND HAVE IT VISIT n RANDOMLY SELECTED PLANTS. 
	int i,j,BreakPointA[200],BreakPointP[200],RanP,RanA,Samp,k,RandomEncounters,WalkPlants,Winner,Successes;
	int InteractingAnimalSpecies,InteractingPlantSpecies,InteractingAnimalIndividual,InteractingPlantIndividual,TotalHabitat;
	double PI,ExpectedVisits,FractionOccupied;
	// the crux challenge here is how to draw individuals at random in proportion to their frequency...

	// First clear out the biotic fitness vectors
	for (PS = 0; PS < PSpecies; PS++)
	{
		for (k = 0; k < NP[PS]; k++)
		{
			WBioPlant[PS][k] = 0;
		}
	}
	for (AS = 0; AS < ASpecies; AS++)
	{
		for (k = 0; k < NA[AS]; k++)
		{
			WBioAnimal[AS][k] = 0;
		}
	}
	//initialize interaction matrix
	for(PS=0;PS<PSpecies;PS++)
	{
		for(AS=0;AS<ASpecies;AS++)
		{
			IMAT[PS][AS]=0;
		}
	}

	//NEW CODE GOES IN HERE
	WalkPlants = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		for (k = 0; k < NP[PS]; k++)
		{
			PlantSpecies[WalkPlants] = PS;
			PlantIndividual[WalkPlants] = k;
			WalkPlants++;
		}
	}
	FractionOccupied = WalkPlants / (1.0*TotalHabitatP);
	//cout << WalkPlants << "\n";
	//Check to make sure we don't overload array
	if (WalkPlants > 1000000)
	{
		cout << "You Busted your ARRAY\n";
		cin >> Test;
	}

	//PUT A CAP ON THIS... A LIFETIME MAXIMUM IF YOU WILL. HOW FAR CAN THE BEE GO?
	
	ExpectedVisits = 2.00*FractionOccupied;//Each animal can visit this NUMBER OF PLANTS WHICH SATURATES


	//done checking for array overload
	Encounters = 0;
	SuccessfulEncounters = 0;
	Samp = 0;
	Successes = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		for (k = 0; k < NA[AS]; k++)
		{
			
			if (ExpectedVisits > 0)
			{
				RandomEncounters = Fish_Rand(ExpectedVisits);
			}
			else
			{
				RandomEncounters = 0;
			}
			while(RandomEncounters>0)
			{
				Encounters++;
				//Choose a random Plant for this animal to interact with
				if (WalkPlants > 1)
				{
					Winner = Big_Rand(WalkPlants - 1);
				}
				else
				{
					Winner = 0;
				}
				

				//Now determine whether an interaction occurs between these individuals
				if (ModelType == 0)//do matching
				{
					PI = exp(-AlphaM*pow((ZP[PlantSpecies[Winner]][PlantIndividual[Winner]] - ZA[AS][k]), 2.0));
				}
				if (ModelType == 1)//do escalation
				{
					//This function assumes plants want to be smaller than animules. Changed on 12/25/2010
					PI = 1.0 / (1.0 + exp(-AlphaM*(ZA[AS][k] - ZP[PlantSpecies[Winner]][PlantIndividual[Winner]])));
				}
				//then modify fitness if the interaction occured
				if ((rand() / (1.0*RAND_MAX)) < PI)
				{
					//then assign a change in fitness if the interaction occurred. Since multiple encouonters are possible, this should really be changed to be additive with respect to interactions (but multiplicative across abiotic and biotic components)

					//NEED TO SATURATE FITNESS SO THAT INFINITE FITNESS DOES NOT RESULT FROM LOTS OF VISITS??

					WBioPlant[PlantSpecies[Winner]][PlantIndividual[Winner]] = WBioPlant[PlantSpecies[Winner]][PlantIndividual[Winner]] + EpsP[PlantSpecies[Winner]];
					WBioAnimal[AS][k] = WBioAnimal[AS][k] + EpsA[AS];

					Successes++;
					SuccessfulEncounters++;
					//and record the fact that this interaction did occur
				
					IMAT[PlantSpecies[Winner]][AS]++;
					ComplementaritySampleA[Samp] = ZA[AS][k];//Note that sampling is no longer random with respect to animals so this needs to be repaired if it is to be used
					ComplementaritySampleP[Samp] = ZP[PlantSpecies[Winner]][PlantIndividual[Winner]];//Note that sampling is no longer random with respect to animals so this needs to be repaired if it is to be used
					Samp++;				
				}
				RandomEncounters--;//count down the number of encounters for each animal
			}//closes the number of random encounters per animal loop
		}//closes the animal individual loop
	}//closes the animal species loop
	Efficiency = (1.0*Successes)/ (1.0*Encounters);
	PerPlantSuccesses = SuccessfulEncounters / (1.0*WalkPlants);
	PerPlantEncounters = Encounters / (1.0*WalkPlants);
}//end interactions loop


void Mate() 
{

	int i, j, k, Mom, Dad, TotalBabies, RemainingBabies, Winner, Seeker, Niche, SpacesRemaining[100],OffProd;
	double u, v, r2, scale, Noise, Luck, PSurvive,MeanOffA,MeanOffP,MidP;

	//First, calculate the total fitness of each individual by multiplying abiotic and biotic components of fitness
	for (PS = 0; PS < PSpecies; PS++)
	{
		for (k = 0; k < NP[PS]; k++)
		{
			WTotPlant[PS][k]= WAbioPlant[PS][k]*WBioPlant[PS][k];
		}
	}
	for (AS = 0; AS < ASpecies; AS++)
	{
		for (k = 0; k < NA[AS]; k++)
		{
			WTotAnimal[AS][k] = WAbioAnimal[AS][k]* WBioAnimal[AS][k];
			//cout << WTotAnimal[AS][k] << "\n";
			//cin >> Test;
		}
	}
	
	
	//For accuracy, calculate mean phenotypes of each species (used in heritability calcs)
	for (PS = 0; PS < PSpecies; PS++)
	{
		PMean[PS] = 0; 
		for (k = 0; k < NP[PS]; k++)
		{
			PMean[PS] = PMean[PS] + ZP[PS][k];			
		}
		PMean[PS] = PMean[PS] / NP[PS];
	}
	for (AS = 0; AS < ASpecies; AS++)
	{
		AMean[AS] = 0;
		for (k = 0; k < NA[AS]; k++)
		{
			AMean[AS] = AMean[AS] + ZA[AS][k];
		}
		AMean[AS] = AMean[AS] / NA[AS];
	}
	
	
	for (i = 0; i < ASpecies; i++)
	{
		SpacesRemaining[i] = KA[i];
	}

	TotalBabies = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		for (i = 0; i < NA[AS]; i++)//Step through all moms
		{
			//USE FITNESS OF MOM TO DETERMINE HOW MANY BABIES SHE MAKES... CAP TOTAL # OFFSPRING PER FEMALE at 10
			if (WTotAnimal[AS][i] < 10)
			{
				MeanOffA = WTotAnimal[AS][i];
			}
			else
			{
				MeanOffA = 10;
			}

			if (MeanOffA > 0)
			{
				OffProd = Fish_Rand(MeanOffA);
				for (j = 0; j < OffProd; j++)
				{
					BabySpecies[TotalBabies] = AS;
					BabyMother[TotalBabies] = i;
					TotalBabies++;
				}
			}
		}
	}
	//OK. I now know how many babies will be produced by each mother within each species. I will use this information to populate the next generation.


	//cout << TotalBabies << "\n";
	//cin >> Test;

	//Initialize the next generation number counter
	for (AS = 0; AS < ASpecies; AS++)
	{
		NAP[AS] = 0;
	}
	//Now, draw mom at random but in proportion to her offspring production. Then mate her with a random father, produce an offspring and drop it into a niche to see if it survives.
	RemainingBabies = TotalBabies;//ONE MIGHT WORRY THIS IS SAMPLING WITH REPLACEMENT, BUT THE BABIES DON"T ACTUALLY EXIST, WE JUST USE THEIR NUMBERS TO SET BABIES PROPORTIONAL TO FITNESS
	while (RemainingBabies > 0)
	{
		if (TotalBabies > 1)
		{
			Winner = Big_Rand(TotalBabies - 1);//Pick a random mother  // could speed this up but use more memory if I just stored the baby info a 3d array...
		}
		else
		{
			Winner = 0;
		}
	//	cout << Winner << "\n";
													
		
		if (SpacesRemaining[BabySpecies[Winner]] > 0)//spaces remain in the chosen species' niche
		{
			Niche = BabySpecies[Winner];
			PSurvive = 1.0;
		}
		else
		{
			Niche = (ASpecies-1) * rand() / (RAND_MAX);//choose a random niche
			if (SpacesRemaining[Niche] > 0)//check to be sure it is not full
			{
				PSurvive = 0.20;
			}
			else
			{
				PSurvive = 0.0;
			}
		}
		Luck = 1.0* rand() / (1.0*RAND_MAX);
		if (Luck < PSurvive)
		{
			//draw a random father (selfing allowed)
		//	MAX = NA[AS] - 1;
			if (NA[BabySpecies[Winner]] > 0)
			{
				Dad = Big_Rand(NA[BabySpecies[Winner]]);//why was there a -1 in this?
			}
			else
			{
				Dad = 0;
			}
			//produce a baby
			//draw some random gaussian noise with mean zero and variance SegVar
			do
			{
				// choose x,y in uniform square (-1,-1) to (+1,+1) 
				u = -1 + 2 * rand() / (1.0*RAND_MAX);
				v = -1 + 2 * rand() / (1.0*RAND_MAX);
				// see if it is in the unit circle 
				r2 = u * u + v * v;
			} while (r2 > 1.0 || r2 == 0);
			scale = sqrt(-2.0 * log(r2) / r2);
			Noise = sqrt(SegVarAnimal) * u * scale;
			MidP = (ZA[BabySpecies[Winner]][BabyMother[Winner]] + ZA[BabySpecies[Winner]][Dad]) / 2.0;

			ZAP[BabySpecies[Winner]][NAP[BabySpecies[Winner]]] = AMean[BabySpecies[Winner]]+hA[BabySpecies[Winner]]*(MidP- AMean[BabySpecies[Winner]]) + Noise;
			NAP[BabySpecies[Winner]]++;
			SpacesRemaining[Niche]--;

	
			
		}//END IF LUCK
		RemainingBabies--;


		//cout << RemainingBabies << "\n";


	}//Closes the while loop which runs until a number of babies equal to the total offspring production of the community has been produced
	//cout << "TO HERE\n";
	//cin >> Test;
	//We must finish by updating the phenotype and abundance arrays
	for(AS=0; AS < ASpecies;AS++)
	{
		for (i = 0; i < NAP[AS]; i++)
		{
			ZA[AS][i] = ZAP[AS][i];
			//cout << ZA[AS][i] << "\n";
			//cin >> Test;
		}
		NA[AS] = NAP[AS];	
		if (NA[AS] == 0 && AExtinct[AS] == -1)
		{
			AExtinct[AS] = G;
		}
	}
	//In principle, all is now complete FOR ANIMALS
	//cout << TotalBabies << "\n";
	//cout << TempCountForDebug1 << "\n";
	//cout << TempCountForDebug2 << "\n";
	//cin >> Test;



	// NOW DO PLANTS
	for (i = 0; i < PSpecies; i++)
	{
		SpacesRemaining[i] = KP[i];
	}

	TotalBabies = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		for (i = 0; i < NP[PS]; i++)//Step through all moms
		{
			//USE FITNESS OF MOM TO DETERMINE HOW MANY BABIES SHE MAKES... CAP TOTAL # OFFSPRING PER FEMALE at 10
			if (WTotPlant[PS][i] < 10)
			{
				MeanOffP = WTotPlant[PS][i];
			}
			else
			{
				MeanOffP = 10;
			}

			if (MeanOffP > 0)
			{
				OffProd = Fish_Rand(MeanOffP);
				for (j = 0; j < OffProd; j++)
				{
					BabySpecies[TotalBabies] = PS;
					BabyMother[TotalBabies] = i;
					TotalBabies++;
				}
			}
		}
	}
	//OK. I now know how many babies will be produced by each mother within each species. I will use this information to populate the next generation.


	//cout << TotalBabies << "\n";
	//cin >> Test;

	//Initialize the next generation number counter
	for (PS = 0; PS < PSpecies; PS++)
	{
		NPP[PS] = 0;
	}
	//Now, draw mom at random but in proportion to her offspring production. Then mate her with a random father, produce an offspring and drop it into a niche to see if it survives.
	RemainingBabies = TotalBabies;//ONE MIGHT WORRY THIS IS SAMPLING WITH REPLACEMENT, BUT THE BABIES DON"T ACTUALLY EXIST, WE JUST USE THEIR NUMBERS TO SET BABIES PROPORTIONAL TO FITNESS
	while (RemainingBabies > 0)
	{
		if (TotalBabies > 1)
		{
			Winner = Big_Rand(TotalBabies - 1);//Pick a random mother  // could speed this up but use more memory if I just stored the baby info a 3d array...
		}
		else
		{
			Winner = 0;
		}
		if (Winner < 0)
		{
			cout << TotalBabies << "," << (TotalBabies - 1) << "," << Winner << "\n";
			cout<<"Fuck Me, it is\n";
			cin >> Test;
		}
		if (SpacesRemaining[BabySpecies[Winner]] > 0)//spaces remain in the chosen species' niche
		{
			Niche = BabySpecies[Winner];
			PSurvive = 1.0;
		}
		else
		{
			Niche = PSpecies * rand() / (RAND_MAX);//choose a random niche
			if (SpacesRemaining[Niche] > 0)//check to be sure it is not full
			{
				PSurvive = 0.20;
			}
			else
			{
				PSurvive = 0.0;
			}
		}
		Luck = 1.0* rand() / (1.0*RAND_MAX);
		if (Luck < PSurvive)
		{
			if (NP[BabySpecies[Winner]] > 1)
			{
				Dad = Big_Rand(NP[BabySpecies[Winner]]);//why was a 1 subtracted here?
			}
			else
			{
				Dad = 0;
			}
			//produce a baby
			//draw some random gaussian noise with mean zero and variance SegVar
			do
			{
				// choose x,y in uniform square (-1,-1) to (+1,+1) 
				u = -1 + 2 * rand() / (1.0*RAND_MAX);
				v = -1 + 2 * rand() / (1.0*RAND_MAX);
				// see if it is in the unit circle 
				r2 = u * u + v * v;
			} while (r2 > 1.0 || r2 == 0);
			scale = sqrt(-2.0 * log(r2) / r2);
			Noise = sqrt(SegVarPlant) * u * scale;
			MidP = (ZP[BabySpecies[Winner]][BabyMother[Winner]] + ZP[BabySpecies[Winner]][Dad]) / 2.0;

			ZPP[BabySpecies[Winner]][NPP[BabySpecies[Winner]]] = PMean[BabySpecies[Winner]]+ hP[BabySpecies[Winner]]*(MidP- PMean[BabySpecies[Winner]]) + Noise;
			NPP[BabySpecies[Winner]]++;
			SpacesRemaining[Niche]--;

		}
		RemainingBabies--;


	}//Closes the while loop which runs until a number of babies equal to the total offspring production of the community has been produced
	 //We must finish by updating the phenotype and abundance arrays
	for (PS=0; PS < PSpecies; PS++)
	{
		for (i = 0; i < NPP[PS]; i++)
		{
			ZP[PS][i] = ZPP[PS][i];
		}
		NP[PS] = NPP[PS];
		if (NP[PS] == 0 && PExtinct[PS] == -1)
		{
			PExtinct[PS] = G;
		}
	}
	//In principle, all is now complete for PLANTS
			
			

}

void PopStats()
{
	double Sum,SumSq,OtherSum,FP[100],FA[100],pA,pP,PSD,ASD;
	int i;

	ExtantPlants = 0;
	ObservablePlants = 0;
	for(PS=0;PS<PSpecies;PS++)
	{
		if (NP[PS] > 50)
		{
			ObservablePlants++;
		}
		if (NP[PS] > 0)
		{
			ExtantPlants++;
			Sum = 0;
			SumSq = 0;
			for (i = 0; i < NP[PS]; i++)
			{
				Sum = Sum + ZP[PS][i];
				SumSq = SumSq + ZP[PS][i] * ZP[PS][i];
			}

			PEvo[PS]= fabs(log(PMean[PS])-log(Sum / (1.0*NP[PS])));//calculates the change in mean phenotype that occured over the generation...
			PNullEvo[PS] = sqrt((hP[PS]*PVar[PS]+SegVarPlant) / NP[PS])*sqrt(2.0 / 3.141592);
			PMean[PS] = Sum / (1.0*NP[PS]);
			PVar[PS] = SumSq / (1.0*NP[PS]) - PMean[PS] * PMean[PS];
		}
		else
		{
			PMean[PS] = 0;
			PVar[PS] = 0;
			PEvo[PS] = 0;

		}

	}

	ExtantAnimals = 0;
	ObservableAnimals = 0;
	for(AS=0;AS<ASpecies;AS++)
	{
		if (NA[AS] > 50)
		{
			ObservableAnimals++;
		}
		if (NA[AS] > 0)
		{
			ExtantAnimals++;
			Sum = 0;
			SumSq = 0;
			for (i = 0; i < NA[AS]; i++)
			{
				Sum = Sum + ZA[AS][i];
				SumSq = SumSq + ZA[AS][i] * ZA[AS][i];
			}
			AEvo[AS] = fabs(log(AMean[AS]) - log(Sum / (1.0*NA[AS])));
			ANullEvo[AS] = sqrt((hA[AS]*AVar[AS]+SegVarAnimal) / NA[AS])*sqrt(2.0 / 3.141592);
			AMean[AS] = Sum / (1.0*NA[AS]);
			AVar[AS] = SumSq / (1.0*NA[AS]) - AMean[AS] * AMean[AS];
		}
		else
		{
			AMean[AS] = 0;
			AVar[AS] = 0;
			AEvo[AS] = 0;
		}
	}

	//now calculate meta statistics for EXTANT plants
	TotalPlants = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		TotalPlants = TotalPlants + NP[PS];
	}
	for (PS = 0; PS < PSpecies; PS++)
	{
		FP[PS] =1.0*NP[PS]/(1.0*TotalPlants);
	}

	PMetaMean = 0;
	PThetaMetaMean = 0;
	PKMetaMean = 0;
	PMetaMeanEvo = 0;
	PMetaMeanNullEvo = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		PMetaMeanEvo = PMetaMeanEvo + FP[PS] * PEvo[PS] ;
		PMetaMeanNullEvo= PMetaMeanNullEvo +  FP[PS] * PNullEvo[PS];
		PMetaMean = PMetaMean + FP[PS] * PMean[PS];
		PThetaMetaMean = PThetaMetaMean +FP[PS]* ThetaP[PS];
		PKMetaMean = PKMetaMean + FP[PS] * KP[PS];
	}

	PMetaVar = 0;
	PThetaMetaVar = 0;
	PKMetaVar = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		PMetaVar = PMetaVar + FP[PS] * (PMean[PS] - PMetaMean) * (PMean[PS] - PMetaMean);
		PThetaMetaVar = PThetaMetaVar + FP[PS] * (ThetaP[PS] - PThetaMetaMean) * (ThetaP[PS] - PThetaMetaMean);
		PKMetaVar = PKMetaVar + FP[PS] * (KP[PS] - PKMetaMean)*(KP[PS] - PKMetaMean);
	}



	//Now do animals
	TotalAnimals = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		TotalAnimals = TotalAnimals + NA[AS];
	}
	for (AS = 0; AS < ASpecies; AS++)
	{
		FA[AS] = 1.0*NA[AS] / (1.0*TotalAnimals);
	}

	AMetaMean = 0;
	AThetaMetaMean = 0;
	AKMetaMean = 0;
	AMetaMeanEvo = 0;
	AMetaMeanNullEvo = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		AMetaMeanEvo = AMetaMeanEvo + FA[AS] * AEvo[AS];
		AMetaMeanNullEvo = AMetaMeanNullEvo +  FA[AS] * ANullEvo[AS];
		AMetaMean = AMetaMean + FA[AS] * AMean[AS];
		AThetaMetaMean = AThetaMetaMean + FA[AS] * ThetaA[AS];
		AKMetaMean = AKMetaMean + FA[AS] * KA[AS];
	}

	AMetaVar = 0;
	AThetaMetaVar = 0;
	AKMetaVar = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		AMetaVar = AMetaVar + FA[AS] * (AMean[AS] - AMetaMean) * (AMean[AS] - AMetaMean);
		AThetaMetaVar = AThetaMetaVar + FA[AS] * (ThetaA[AS] - AThetaMetaMean) * (ThetaA[AS] - AThetaMetaMean);
		AKMetaVar = AKMetaVar + FA[AS] * (KA[AS] - AKMetaMean)*(KA[AS] - AKMetaMean);
	}


	//Diversity Indices
	HA = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		if (FA[AS] > 0)
		{
			HA = HA - FA[AS] * log(FA[AS]);
		}
	}
	HP = 0;

	for (PS = 0; PS < PSpecies; PS++)
	{
		if (FP[PS] > 0)
		{
			HP = HP - FP[PS] * log(FP[PS]);
		}
	}
	//Theoretical maxima for perfectly even communities of current number of species
	HAmax = 0;
	pA = ((1.0*TotalAnimals) / (1.0*ASpecies))/(1.0*TotalAnimals);
	for (AS = 0; AS < ASpecies; AS++)
	{
		HAmax = HAmax - pA* log(pA);
	}
	HPmax = 0;

	pP = ((1.0*TotalPlants) / (1.0*PSpecies)) / (1.0*TotalPlants);
	for (PS = 0; PS < PSpecies; PS++)
	{
		HPmax = HPmax - (pP) * log(pP);
	}
	//Evenness
	EvenA = HA / HAmax;
	EvenP = HP / HPmax;

	//Total Number of Unique Interactions


}

void InitializeOutput()
{
	int Samp,i;

	out_meta << "DisturbType,EvoHist,PropDestroyedP,PropDestroyedA,HARoot,HPRoot,AssemblyTime,G,ZbarP,ZvarP,ZbarA,ZvarA,,TbarP,TvarP,TbarA,TvarA,,KbarP,KvarP,KbarA,KvarA,,Encounters/Plant,Successful Encounters/Plant,Connectance,Efficiency,UniqueInteractions,Total # Successful Interactions,MeanADeg,MeanPDeg,MeanGeneralization,InteractionDiversity,,ExtantAnimals,ExtantPlants,Total Animals,TotalPlants,HA,HP,EvenA,EvenP\n";


	out_traits<<"PRun,HRun,G"<<",";
	for(PS=0;PS<PSpecies;PS++)
	{
		out_traits<<"PS#"<<PS<<",";
	}

	for(AS=0;AS<ASpecies;AS++)
	{
		out_traits<<"AS#"<<AS<<",";
	}
	out_traits<<"\n";

	out_abundance << "PRun,HRun,G" << ",";
	for (PS = 0; PS<PSpecies; PS++)
	{
		out_abundance << "PS#" << PS << ",";
	}

	for (AS = 0; AS<ASpecies; AS++)
	{
		out_abundance << "AS#" << AS << ",";
	}
	out_abundance << "\n";

	out_K << "PRun,HRun,G" << ",";
	for (PS = 0; PS<PSpecies; PS++)
	{
		out_K << "PS#" << PS << ",";
	}

	for (AS = 0; AS<ASpecies; AS++)
	{
		out_K << "AS#" << AS << ",";
	}
	out_K << "\n";

	out_T << "PRun,HRun,G" << ",";
	for (PS = 0; PS<PSpecies; PS++)
	{
		out_T << "PS#" << PS << ",";
	}

	for (AS = 0; AS<ASpecies; AS++)
	{
		out_T << "AS#" << AS << ",";
	}
	out_T << "\n";

	/*
	out_sampledtraits<<"Run,G"<<",";
	for(i=0;i<PSpecies;i++)
	{
		if(PlantSampled[i]==1)
		{
			out_sampledtraits<<"PS#"<<i<<",";
		}
	}

	for(i=0;i<ASpecies;i++)
	{
		if(AnimalSampled[i]==1)
		{
			out_sampledtraits<<"AS#"<<Samp<<",";
		}
	}
	out_sampledtraits<<"\n";
	*/


	/*
	out_Binteract<<"Run,G,";
	for(AS=0;AS<ASpecies;AS++)
	{
		out_Binteract<<"AS#"<<AS<<",";
	}
	out_Binteract<<"\n";
	*/
	
	/*
	out_trait_degree<<"Run,G,Animal Species"<<","<<"Population Size"<<","<<"Trait Mean"<<","<<"Degree"<<","<<"Plant Species"<<","<<"Population Size"<<","<<"Trait Mean"<<","<<"Degree"<<"\n";
	out_network_summary<<"Model,NullModel,Run,G,Sampled,AlphaM,AlphaC,GammaA,GammaP,TrueMeanEpsA,TrueMeanEpsP,MeanEpsA,MeanEpsP,VarEpsA,VarEpsP,Heat,MeanKA,MeanKP,VarKA,VarKP,MeanThetaA,MeanThetaP,VarThetaA,VarThetaP,MeanSA,MeanSP,NullConnectance,NullConvergeA,NullConvergeP,NullComplementarity,NullNODFA,NullNODFP,ObservedConnectance,Efficiency,ObservedConvergeA,CorrectedConvergeA,ConvergeA_Percentile,ObservedConvergeP,CorrectedConvergeP,ConvergeP_Percentile,ObservedComplementarity,CorrectedComplementarity,Complementarity_Percentile,ObservedNODFA,ObservedNODFP,CorrectedNODFA,CorrectedNODFP,NODFA_Percentile,NODFP_Percentile\n";
	out_pred<<"Run,G,GammaA,GammaP,SA,SP,MeanThetaA,MeanThetaP,VarThetaA,VarThetaP,AverageVarWithinA,AverageVarWithinP,AnalyticalMeanA,ObservedMeanA,AnalyticalMeanP,ObservedMeanP,AnalyticalVarA,ObservedVarA,AnalyticalVarP,ObservedVarP,AnalyticalConnectance,Efficiency,ObservedConnectance,AnalyticalIPMeanA,ObservedIPMeanA,AnalyticalIPMeanP,ObservedIPMeanP,AnalyticalIPVA,ObservedIPVA,AnalyticalIPVP,ObservedIPVP,AnalyticalComplementarity,ObservedComplementarity\n";
	*/
}
void Output()
{
	
	out_meta << DisturbType<<","<<EH<<","<<PropLostP << "," << PropLostA <<","<<hARoot<<","<< hPRoot << "," << AssemblyTime << "," << G << "," << PMetaMean<< "," <<PMetaVar<< "," << AMetaMean << "," << AMetaVar << ",," << PThetaMetaMean << "," << PThetaMetaVar << "," << AThetaMetaMean << "," << AThetaMetaVar << ",," << PKMetaMean << "," << PKMetaVar << "," << AKMetaMean << "," << AKMetaVar<<",,"<<PerPlantEncounters<<","<<PerPlantSuccesses<<","<<C<<","<<Efficiency<<","<< UniqueInteractions<<","<<SuccessfulEncounters<<","<<MeanADeg<<","<<MeanPDeg<<","<<MeanGeneral<<","<<HInteract<<",,"<<ObservableAnimals<<","<<ObservablePlants<<","<< TotalAnimals<<","<< TotalPlants<<","<< HA<<","<< HP<<","<< EvenA<<","<< EvenP<<"\n";

	out_meta.flush();
	out_traits<<pRun<<","<<hRun<<","<<G<<",";
	
	for(PS=0;PS<PSpecies;PS++)
	{
		if (NP[PS] > 0)
		{
			out_traits << PMean[PS] << ",";
		}
		else
		{
			out_traits <<  ",";
		}
	}

	for(AS=0;AS<ASpecies;AS++)
	{
		if (NA[AS] > 0)
		{
			out_traits << AMean[AS] << ",";
		}
		else
		{
			out_traits << ",";
		}
	}
	out_traits<<"\n";

	out_traits.flush();

	out_abundance << pRun << "," << hRun  << "," << G << ",";

	for (PS = 0; PS<PSpecies; PS++)
	{
		out_abundance << NP[PS] << ",";
	}

	for (AS = 0; AS<ASpecies; AS++)
	{
		out_abundance << NA[AS] << ",";
	}
	out_abundance << "\n";

	out_abundance.flush();


	out_K << pRun << "," << hRun << "," << G << ",";

	for (PS = 0; PS<PSpecies; PS++)
	{
		
		out_K << KP[PS] << ",";
	
	}

	for (AS = 0; AS<ASpecies; AS++)
	{

		out_K << KA[AS] << ",";
	
	}
	out_K << "\n";

	out_K.flush();

	out_T << pRun << "," << hRun << "," << G << ",";

	for (PS = 0; PS<PSpecies; PS++)
	{

		out_T << ThetaP[PS] << ",";

	}

	for (AS = 0; AS<ASpecies; AS++)
	{

		out_T << ThetaA[AS] << ",";

	}
	out_T << "\n";

	out_T.flush();

	
	/*
	out_sampledtraits<<Run<<","<<G<<",";
	for(i=0;i<PSpecies;i++)
	{
		if(PlantSampled[i]==1)
		{
			out_sampledtraits<<PMean[i]<<",";
		}
	}
	for(i=0;i<ASpecies;i++)
	{
		if(AnimalSampled[i]==1)
		{
			out_sampledtraits<<AMean[i]<<",";
		}
	}
	out_sampledtraits<<"\n";

	out_sampledtraits.flush();
	*/
}

/*
void OutputInteract()
{
	out_Qinteract << Run<<","<<G<<",";
	for (AS = 0; AS<ASpecies; AS++)
	{
		for (PS = 0; PS < PSpecies; PS++)
		{
			out_Qinteract << FIMAT[PS][AS]<<",";
		}
	}
	out_Qinteract << ",,";

	for (AS = 0; AS<ASpecies; AS++)
	{
		for (PS = 0; PS < PSpecies; PS++)
		{
			out_Qinteract << General[PS][AS]<<",";
		}
	}
	out_Qinteract << "\n";

	out_Qinteract.flush();

	/*
	for(PS=0;PS<PSpecies;PS++)
	{
		out_Qinteract<<Run<<","<<G<<",";
		for(AS=0;AS<ASpecies;AS++)
		{
			out_Qinteract<<IMAT[PS][AS]<<",";
		}
		out_Qinteract<<"\n";
	}


	
	for(PS=0;PS<PSpecies;PS++)
	{
		out_Binteract<<Run<<","<<G<<",";
		for(AS=0;AS<ASpecies;AS++)
		{
			out_Binteract<<BIMAT[PS][AS]<<",";
		}
		out_Binteract<<"\n";
	}
	

	//The below outputs summary stats only

	out_network_summary<<ModelType<<","<<PureNull<<","<<Run<<","<<G<<","<<Sampled<<","<<AlphaM<<","<<AlphaC<<","<<MeanGammaA<<","<<MeanGammaP<<","<<TrueMeanEpsA<<","<<TrueMeanEpsP<<","<<MeanEpsA<<","<<MeanEpsP<<","<<VarEpsA<<","<<VarEpsP<<","<<MeanEpsA*MeanEpsP*AlphaM<<","<<AKMetaMean<<","<<PKMetaMean<<","<<AKMetaVar<<","<<PKMetaVar<<","<<AThetaMetaMean<<","<<PThetaMetaMean<<","<<AThetaMetaVar<<","<<PThetaMetaVar<<","<<MeanSA<<","<<MeanSP<<","<<ENullC<<","<<ENullVA<<","<<ENullVP<<","<<ENullComplementarity<<","<<ENullNODFA<<","<<ENullNODFP<<","<<C<<","<<Efficiency<<","<<IPVA<<","<<CorrectVA<<","<<VAPercentile<<","<<IPVP<<","<<CorrectVP<<","<<VPPercentile<<","<<IPComplementarity<<","<<CorrectComp<<","<<CPercentile<<","<<NODFA<<","<<NODFP<<","<<CorrectNODFA<<","<<CorrectNODFP<<","<<APercentLesser<<","<<PPercentLesser<<"\n";
	

	//the below outputs data on the relationship between species mean trait value and its degree
	if(ASpecies>=PSpecies)
	{
		for(AS=0;AS<ASpecies;AS++)
		{
			out_trait_degree<<Run<<","<<G<<","<<AS<<","<<KA[AS]<<","<<AMean[AS]<<","<<DegA[AS]<<",";
			if(AS<PSpecies)
			{
				out_trait_degree<<AS<<","<<KP[AS]<<","<<PMean[AS]<<","<<DegP[AS]<<"\n";
			}
			else
			{
				out_trait_degree<<""<<","<<""<<","<<""<<","<<""<<"\n";//output blanks
			}
		}
	}
	if(ASpecies<PSpecies)
	{
		for(PS=0;PS<PSpecies;PS++)
		{
			out_trait_degree<<Run<<","<<G<<","<<PS<<","<<KP[PS]<<","<<PMean[PS]<<","<<DegP[PS]<<",";
			if(PS<ASpecies)
			{
				out_trait_degree<<PS<<","<<KA[PS]<<","<<AMean[PS]<<","<<DegA[PS]<<"\n";
			}
			else
			{
				out_trait_degree<<""<<","<<""<<","<<""<<","<<""<<"\n";//output blanks
			}
		}
	}

	out_Binteract.flush();
	out_Qinteract.flush();
	out_network_summary.flush();
	out_trait_degree.flush();
	

}*/

void AnalyzeInteractions()
{
	int ASP,PSP,Samp,Extants;
	double NODINC,SumA,SumP,SumAP,SumAA,SumPP;
	//first, calculate frequency of interactions
	NumIntObs = 0;
	for (PS = 0; PS<PSpecies; PS++)
	{
		for (AS = 0; AS<ASpecies; AS++)
		{
			NumIntObs = NumIntObs + IMAT[PS][AS];
		}
	}

	for (PS = 0; PS<PSpecies; PS++)
	{
		for (AS = 0; AS<ASpecies; AS++)
		{
			FIMAT[PS][AS] = IMAT[PS][AS] /(1.0*NumIntObs);
		}
	}

	//CALCULATE AN INDEX OF INTERACTION DIVERSITY USING SHANNON WEINER ANALOGUE
	//Diversity Indices
	HInteract = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		for (PS = 0; PS < PSpecies; PS++)
		{
			if (FIMAT[PS][AS] > 0)
			{
				HInteract = HInteract - FIMAT[PS][AS] * log(FIMAT[PS][AS]);
			}
		}
	}


	


	//Next, make a binary matrix and calculate connectance
	C=0;
	UniqueInteractions = 0;
	NumIntObs=0;
	Extants = 0;
	for(PS=0;PS<PSpecies;PS++)
	{
		for(AS=0;AS<ASpecies;AS++)
		{
			if (NP[PS] > 0 && NA[AS] > 0)
			{
				Extants++;
			}
			BIMAT[PS][AS]=0;
			//AnimalSampled[AS]=0;
			//PlantSampled[PS]=0;
			if(IMAT[PS][AS]>0&& NP[PS] > 0 && NA[AS] > 0)
			{
				BIMAT[PS][AS]=1;
				//AnimalSampled[AS]=1;
				//PlantSampled[PS]=1;
				UniqueInteractions++;
				C++;
				NumIntObs++;
			}
		}
	}
	C=C/(1.0*Extants); //WARNING WARNING... AREN'T THESE CALCULATIONS SCREWED DUE TO EXTINCTION?  ALSO FOR GENERALIZATION?

	//next, calculate degree distribution for animals and plants
	MeanPDeg = 0;
	for(PS=0;PS<PSpecies;PS++)
	{
		DegP[PS]=0;
		for(AS=0;AS<ASpecies;AS++)
		{
			DegP[PS]=DegP[PS]+BIMAT[PS][AS];			
		}
		MeanPDeg = MeanPDeg + DegP[PS];
	}
	MeanPDeg = MeanPDeg / (1.0*ExtantPlants);

	MeanADeg = 0;
	for(AS=0;AS<ASpecies;AS++)
	{
		DegA[AS]=0;
		for(PS=0;PS<PSpecies;PS++)
		{
			DegA[AS]=DegA[AS]+BIMAT[PS][AS];			
		}
		MeanADeg = MeanADeg + DegA[AS];
	}
	MeanADeg = MeanADeg / (1.0*ExtantAnimals);

	//Now calculate generalization of interactions sensu Aizen
	MeanGeneral = 0;
	Extants= 0;
	for (PS = 0; PS<PSpecies; PS++)
	{
		for (AS = 0; AS<ASpecies; AS++)
		{
			if (NP[PS] > 0 && NA[AS] > 0)
			{
				Extants++;
				General[PS][AS] = (DegP[PS] + DegA[AS]) / 2.0;
				MeanGeneral = MeanGeneral + General[PS][AS];
			}
		}
	}
	MeanGeneral = MeanGeneral / (1.0*Extants);
	/*
	//next, calculate nestedness using NODF
	//First calculate animal nestedness
	NODFA=0;
	for(AS=0;AS<ASpecies;AS++)
	{
		for(ASP=0;ASP<ASpecies;ASP++)
		{
			if(DegA[AS]>DegA[ASP]&&DegA[ASP]!=0)
			{
				NODINC=0;
				for(PS=0;PS<PSpecies;PS++)
				{
					if(BIMAT[PS][ASP]==1&&BIMAT[PS][AS]==1)
					{
						NODINC++;
					}
				}
				NODFA=NODFA+100.0*NODINC/(1.0*DegA[ASP]);
			}
		}
	}
	NODFA=NODFA/(0.5*ASpecies*(ASpecies-1));
	//Next, calculate plant nestedness
	NODFP=0;
	for(PS=0;PS<PSpecies;PS++)
	{
		for(PSP=0;PSP<PSpecies;PSP++)
		{
			if(DegP[PS]>DegP[PSP]&&DegP[PSP]!=0)
			{
				NODINC=0;
				for(AS=0;AS<ASpecies;AS++)
				{
					if(BIMAT[PS][AS]==1&&BIMAT[PSP][AS]==1)
					{
						NODINC++;
					}
				}
				NODFP=NODFP+100.0*NODINC/(1.0*DegP[PSP]);
			}
		}
	}
	NODFP=NODFP/(0.5*PSpecies*(PSpecies-1));

	//Next, calculate complementarity
	SumA=0;
	SumP=0;
	SumAA=0;
	SumPP=0;
	SumAP=0;
	for(Samp=0;Samp<Sampled;Samp++)
	{
		SumA=SumA+ComplementaritySampleA[Samp];
		SumP=SumP+ComplementaritySampleP[Samp];
		SumAA=SumAA+ComplementaritySampleA[Samp]*ComplementaritySampleA[Samp];
		SumPP=SumPP+ComplementaritySampleP[Samp]*ComplementaritySampleP[Samp];
		SumAP=SumAP+ComplementaritySampleA[Samp]*ComplementaritySampleP[Samp];			
	}
	IPMeanA=SumA/(1.0*Sampled);
	IPMeanP=SumP/(1.0*Sampled);
	IPVA=SumAA/(1.0*Sampled)-(SumA/(1.0*Sampled))*(SumA/(1.0*Sampled));
	IPVP=SumPP/(1.0*Sampled)-(SumP/(1.0*Sampled))*(SumP/(1.0*Sampled));
	IPCovAP=SumAP/(1.0*Sampled)-(SumA/(1.0*Sampled))*(SumP/(1.0*Sampled));
	IPComplementarity=IPCovAP/(sqrt(IPVA)*sqrt(IPVP));
	*/

	

}

void NullNetwork() // This is the "new" null which Pedro and Jordi wanted. Bottom line, it is the standard null.
{
	int i,j,BreakPointA[100],BreakPointP[100],TotalPlants,TotalAnimals,RanP,RanA,Samp;
	int InteractingAnimalSpecies,InteractingPlantSpecies,InteractingAnimalIndividual,InteractingPlantIndividual;
	double PI;
	int ASP,PSP,NullSample;
	double NODINC,SumA,SumP,SumAP,SumAA,SumPP,NullCovAP,NullVA,NullVP;
	int ABigger,ASmaller,PBigger,PSmaller,BiggerC,BiggerVA,BiggerVP;
	

	ENullC=0;
	ENullComplementarity=0;
	ENullNODFA=0;
	ENullNODFP=0;
	ENullVA=0;
	ENullVP=0;

	ABigger=0;
	ASmaller=0;
	PBigger=0;
	PSmaller=0;
	BiggerVA=0;
	BiggerVP=0;
	BiggerC=0;
	NullSample=100;
	for(i=0;i<NullSample;i++)
	{
		//First, generate a random interaction matrix
		//initialize interaction matrix
		for(PS=0;PS<PSpecies;PS++)
		{
			for(AS=0;AS<ASpecies;AS++)
			{
				NullIMAT[PS][AS]=0;
			}
		}
		//first, clear breakpoint vectors
		for(PS=0;PS<PSpecies;PS++)
		{
			BreakPointP[PS]=0;
		}
		for(AS=0;AS<ASpecies;AS++)
		{
			BreakPointA[AS]=0;
		}
		// First, establish break points?
		TotalPlants=KP[0];
		BreakPointP[0]=KP[0];
		for(PS=1;PS<PSpecies;PS++)
		{
			TotalPlants=TotalPlants+KP[PS];
			BreakPointP[PS]=BreakPointP[PS-1]+KP[PS];
		}
		TotalAnimals=KA[0];
		BreakPointA[0]=KA[0];
		for(AS=1;AS<ASpecies;AS++)
		{
			TotalAnimals=TotalAnimals+KA[AS];
			BreakPointA[AS]=BreakPointA[AS-1]+KA[AS];
		}


		//Now, draw an animal and plant at random and interact them. Do this for the observed number of interactions in the real network. This takes into account variation in density!
		
		Samp=0;
		j=0;
		while(j<NumIntObs) //this calculates the null interaction matrix
		{

			//first draw a random plant
			//MAX=TotalPlants-1;
			if (TotalPlants > 1)
			{
				RanP = Big_Rand(TotalPlants - 1);
			}
			else
			{
				RanP = 0;
			}
	
			//which species is this?
			for(PS=0;PS<1;PS++)
			{
				if(RanP<BreakPointP[PS]&&RanP>=0)
				{
					InteractingPlantSpecies=PS;
				}
			}

			for(PS=1;PS<PSpecies;PS++)
			{
				if(RanP<BreakPointP[PS]&&RanP>=BreakPointP[PS-1])
				{
					InteractingPlantSpecies=PS;
				}
			}
			//which individual is this?
			if(InteractingPlantSpecies==0)
			{
				InteractingPlantIndividual=RanP;
			}
			else
			{
				InteractingPlantIndividual=RanP-BreakPointP[InteractingPlantSpecies-1];
			}

			//then draw a random animal
			//MAX=TotalAnimals-1;
			RanA=Big_Rand(TotalAnimals - 1);
			//which species is this?
			for(AS=0;AS<1;AS++)
			{
				if(RanA<BreakPointA[AS]&&RanA>=0)
				{
					InteractingAnimalSpecies=AS;
				}
			}
			for(AS=1;AS<ASpecies;AS++)
			{
				if(RanA<BreakPointA[AS]&&RanA>=BreakPointA[AS-1])
				{
					InteractingAnimalSpecies=AS;
				}
			}
			//which individual is this?
			if(InteractingAnimalSpecies==0)
			{
				InteractingAnimalIndividual=RanA;
			}
			else
			{
				InteractingAnimalIndividual=RanA-BreakPointA[InteractingAnimalSpecies-1];
			}
		
			if(NullIMAT[InteractingPlantSpecies][InteractingAnimalSpecies]==0)
			{
				NullIMAT[InteractingPlantSpecies][InteractingAnimalSpecies]++;
				j++;
				
			}
		}
		
		//then calculate the null complementarity and convergence
		Samp=0;
	
		while(Samp<Sampled) 
		{

			//first draw a random plant
			//MAX=TotalPlants-1;
			RanP=Big_Rand(TotalPlants - 1);
	
			//which species is this?
			for(PS=0;PS<1;PS++)
			{
				if(RanP<BreakPointP[PS]&&RanP>=0)
				{
					InteractingPlantSpecies=PS;
				}
			}

			for(PS=1;PS<PSpecies;PS++)
			{
				if(RanP<BreakPointP[PS]&&RanP>=BreakPointP[PS-1])
				{
					InteractingPlantSpecies=PS;
				}
			}
			//which individual is this?
			if(InteractingPlantSpecies==0)
			{
				InteractingPlantIndividual=RanP;
			}
			else
			{
				InteractingPlantIndividual=RanP-BreakPointP[InteractingPlantSpecies-1];
			}

			//then draw a random animal
			//MAX=TotalAnimals-1;
			RanA=Big_Rand(TotalAnimals - 1);
			//which species is this?
			for(AS=0;AS<1;AS++)
			{
				if(RanA<BreakPointA[AS]&&RanA>=0)
				{
					InteractingAnimalSpecies=AS;
				}
			}
			for(AS=1;AS<ASpecies;AS++)
			{
				if(RanA<BreakPointA[AS]&&RanA>=BreakPointA[AS-1])
				{
					InteractingAnimalSpecies=AS;
				}
			}
			//which individual is this?
			if(InteractingAnimalSpecies==0)
			{
				InteractingAnimalIndividual=RanA;
			}
			else
			{
				InteractingAnimalIndividual=RanA-BreakPointA[InteractingAnimalSpecies-1];
			}
		
			
			NullComplementaritySampleA[Samp]=ZA[InteractingAnimalSpecies][InteractingAnimalIndividual];
			NullComplementaritySampleP[Samp]=ZP[InteractingPlantSpecies][InteractingPlantIndividual];
			Samp++;
		}
	

	


		//now calculate network metrics for this random replicate
		//first, make a binary matrix and calculate connectance
		NullC=0;
		for(PS=0;PS<PSpecies;PS++)
		{
			for(AS=0;AS<ASpecies;AS++)
			{
				NullBIMAT[PS][AS]=0;
				if(NullIMAT[PS][AS]>0)
				{
					NullBIMAT[PS][AS]=1;
					NullC++;
				}
			}
		}
		NullC=NullC/(1.0*PSpecies*ASpecies);
		

		//next, calculate degree distribution for animals and plants
		for(PS=0;PS<PSpecies;PS++)
		{
			NullDegP[PS]=0;
			for(AS=0;AS<ASpecies;AS++)
			{
				NullDegP[PS]=NullDegP[PS]+NullBIMAT[PS][AS];			
			}
		}
		for(AS=0;AS<ASpecies;AS++)
		{
			NullDegA[AS]=0;
			for(PS=0;PS<PSpecies;PS++)
			{
				NullDegA[AS]=NullDegA[AS]+NullBIMAT[PS][AS];			
			}
		}

		//next, calculate nestedness using NODF
		//First calculate animal nestedness
		NullNODFA=0;
		for(AS=0;AS<ASpecies;AS++)
		{
			for(ASP=0;ASP<ASpecies;ASP++)
			{
				if(NullDegA[AS]>NullDegA[ASP]&&NullDegA[ASP]!=0)
				{
					NODINC=0;
					for(PS=0;PS<PSpecies;PS++)
					{
						if(NullBIMAT[PS][ASP]==1&&NullBIMAT[PS][AS]==1)
						{
							NODINC++;
						}
					}
					NullNODFA=NullNODFA+100.0*NODINC/(1.0*NullDegA[ASP]);
				}
			}
		}
		NullNODFA=NullNODFA/(0.5*ASpecies*(ASpecies-1));
		//Next, calculate plant nestedness
		NullNODFP=0;
		for(PS=0;PS<PSpecies;PS++)
		{
			for(PSP=0;PSP<PSpecies;PSP++)
			{
				if(NullDegP[PS]>NullDegP[PSP]&&NullDegP[PSP]!=0)
				{
					NODINC=0;
					for(AS=0;AS<ASpecies;AS++)
					{
						if(NullBIMAT[PS][AS]==1&&NullBIMAT[PSP][AS]==1)
						{
							NODINC++;
						}
					}
					NullNODFP=NullNODFP+100.0*NODINC/(1.0*NullDegP[PSP]);
				}
			}
		}
		NullNODFP=NullNODFP/(0.5*PSpecies*(PSpecies-1));

		//Next, calculate complementarity
		SumA=0;
		SumP=0;
		SumAA=0;
		SumPP=0;
		SumAP=0;
		for(Samp=0;Samp<Sampled;Samp++)
		{
			SumA=SumA+NullComplementaritySampleA[Samp];
			SumP=SumP+NullComplementaritySampleP[Samp];
			SumAA=SumAA+NullComplementaritySampleA[Samp]*NullComplementaritySampleA[Samp];
			SumPP=SumPP+NullComplementaritySampleP[Samp]*NullComplementaritySampleP[Samp];
			SumAP=SumAP+NullComplementaritySampleA[Samp]*NullComplementaritySampleP[Samp];			
		}
		NullVA=SumAA/(1.0*Sampled)-(SumA/(1.0*Sampled))*(SumA/(1.0*Sampled));
		NullVP=SumPP/(1.0*Sampled)-(SumP/(1.0*Sampled))*(SumP/(1.0*Sampled));
		NullCovAP=SumAP/(1.0*Sampled)-(SumA/(1.0*Sampled))*(SumP/(1.0*Sampled));
		NullComplementarity=NullCovAP/(sqrt(NullVA)*sqrt(NullVP));
	
		ENullVA=ENullVA+NullVA;
		ENullVP=ENullVP+NullVP;
		ENullC=ENullC+NullC;
		ENullComplementarity=ENullComplementarity+NullComplementarity;
		ENullNODFA=ENullNODFA+NullNODFA;
		ENullNODFP=ENullNODFP+NullNODFP;

		if(NullNODFA>NODFA)
		{
			ABigger++;
		}
		if(NullNODFA<NODFA)
		{
			ASmaller++;
		}
		if(NullNODFP>NODFP)
		{
			PBigger++;
		}
		if(NullNODFP<NODFP)
		{
			PSmaller++;
		}
		if(IPVA>NullVA)
		{
			BiggerVA++;
		}
		if(IPVP>NullVP)
		{
			BiggerVP++;
		}
		if(IPComplementarity>NullComplementarity)
		{
			BiggerC++;
		}
	}
	ENullVA=ENullVA/(1.0*NullSample);
	ENullVP=ENullVP/(1.0*NullSample);

	ENullC=ENullC/(1.0*NullSample);

	ENullComplementarity=ENullComplementarity/(1.0*NullSample);
	ENullNODFA=ENullNODFA/(1.0*NullSample);
	ENullNODFP=ENullNODFP/(1.0*NullSample);
	APercentGreater=ABigger/(1.0*NullSample);
	APercentLesser=ASmaller/(1.0*NullSample);
	PPercentGreater=PBigger/(1.0*NullSample);
	PPercentLesser=PSmaller/(1.0*NullSample);

	CPercentile=BiggerC/(1.0*NullSample);
	VAPercentile=BiggerVA/(1.0*NullSample);
	VPPercentile=BiggerVP/(1.0*NullSample);



	//Finally, calculate corrected NODF scores by subtracting the mean of the randomized NodfScores
	CorrectNODFA=NODFA-ENullNODFA;
	CorrectNODFP=NODFP-ENullNODFP;
	CorrectComp=IPComplementarity-ENullComplementarity;
	CorrectVA=IPVA-ENullVA;
	CorrectVP=IPVP-ENullVP;


}


void AnPred()
{
		if(ModelType==0)//do matching
		{
				
		}
		if(ModelType==1)//do escalation
		{
			

		}



		out_pred << pRun<<","<<hRun << "," << G << "," << "\n";


}


void SetK()
{
	double SumKP, SumKPP;
	double SumKA, SumKAA;
	double u, v, r2, scale, Noise;

	SumKP = 0;
	SumKPP = 0;

	TotalHabitatP = 0;
	for (PS = 0; PS<PSpecies; PS++)
	{
		//draw K from a lognormal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 * rand() / (1.0*RAND_MAX);
			v = -1 + 2 * rand() / (1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		} while (r2 > 1.0 || r2 == 0);
		scale = sqrt(-2.0 * log(r2) / r2);
		Noise = sqrt(LNVar) * u * scale;
		KProot[PS] = int(exp(LNMean + Noise));
		if (KProot[PS]>9998)
		{
			KProot[PS] = 9999;
		}
		TotalHabitatP = TotalHabitatP + KProot[PS];
	}


	TotalHabitatA = 0;
	for (AS = 0; AS<ASpecies; AS++)
	{
		//draw K from a lognormal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 * rand() / (1.0*RAND_MAX);
			v = -1 + 2 * rand() / (1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		} while (r2 > 1.0 || r2 == 0);
		scale = sqrt(-2.0 * log(r2) / r2);
		Noise = sqrt(LNVar) * u * scale;
		KAroot[AS] = int(exp(LNMean + Noise));
		if (KAroot[AS]>9998)
		{
			KAroot[AS] = 9999;
		}
		TotalHabitatA = TotalHabitatA + KAroot[AS];
	}
}

void SetTheta()
{
	double SumThetaP, SumThetaPP;
	double SumThetaA, SumThetaAA;
	double u, v, r2, scale, Noise;
	int TotKA, TotKP;
	double FracKA[100], FracKP[100], RoA, RoP;

	RoA = 0.0;
	RoP = 0.0;
	//calculate fractional carrying capacities
	TotKP = 0;
	for (PS = 0; PS < PSpecies; PS++)
	{
		TotKP = TotKP + KP[PS];
	}

	for (PS = 0; PS < PSpecies; PS++)
	{
		FracKP[PS]= (1.0*KP[PS])/ (1.0*TotKP);
	}

	//calculate fractional carrying capacities
	TotKA = 0;
	for (AS = 0; AS < ASpecies; AS++)
	{
		TotKA = TotKA + KA[AS];
	}

	for (AS = 0; AS < ASpecies; AS++)
	{
		FracKA[AS] = (1.0*KA[AS]) / (1.0*TotKA);
	}



	for (PS = 0; PS<PSpecies; PS++)
	{
		//Draw theta from a normal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 * rand() / (1.0*RAND_MAX);
			v = -1 + 2 * rand() / (1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		} 
		while (r2 > 1.0 || r2 == 0);
		scale = sqrt(-2.0 * log(r2) / r2);
		Noise = sqrt(ThetaVarP) * u * scale;
		ThetaP[PS] = 2.0 + 2.0*rand() / (1.0*RAND_MAX);;// 5 + RoP*(2.0*FracKP[PS]) + (1 - RoP)* Noise;
	}

	for (AS = 0; AS<ASpecies; AS++)
	{
		//Draw theta from a normal
		do
		{
			// choose x,y in uniform square (-1,-1) to (+1,+1) 
			u = -1 + 2 * rand() / (1.0*RAND_MAX);
			v = -1 + 2 * rand() / (1.0*RAND_MAX);
			// see if it is in the unit circle 
			r2 = u * u + v * v;
		} while (r2 > 1.0 || r2 == 0);
		scale = sqrt(-2.0 * log(r2) / r2);
		Noise = sqrt(ThetaVarA) * u * scale;
		ThetaA[AS] =  2.0 +2.0*rand() / (1.0*RAND_MAX);// 5 + RoA*(2.0*FracKA[AS]) + (1 - RoA)* Noise;
	}
}


void BirthTheta()
{
	int tcount,i,j,SpeciesExtant,Parent,counter;
	double Chance,bA,bP,Tnot,TotalTimeA,TotalTimeP,u,v,r2,scale,Noise, ADriftRate, PDriftRate;

	

	//First, generate a birth tree
	bA = 0.00001; //sets species birth rate
	bP = 0.00001;
	TotalTimeA = (ASpecies / bA);
	TotalTimeP = (PSpecies / bP);
	ADriftRate = 5.0/TotalTimeA;
	PDriftRate = 5.0 / TotalTimeP;

	SpeciesExtant = 1;
	tcount = 0;
	TdivA[0][0] = 0;
	PhyloThetaA[0] = 10;
	

	while (SpeciesExtant <= ASpecies)
	{
		for (i = 0; i < SpeciesExtant; i++)
		{
			//Draw evolutionary change from a normal
			do
			{
				// choose x,y in uniform square (-1,-1) to (+1,+1) 
				u = -1 + 2 * rand() / (1.0*RAND_MAX);
				v = -1 + 2 * rand() / (1.0*RAND_MAX);
				// see if it is in the unit circle 
				r2 = u * u + v * v;
			} while (r2 > 1.0 || r2 == 0);
			scale = sqrt(-2.0 * log(r2) / r2);
			Noise = sqrt(ADriftRate) * u * scale;
			PhyloThetaA[i] = PhyloThetaA[i]+Noise;
		}
		Chance = (1.0*rand() / (1.0*RAND_MAX));
		if (Chance < bA)//do we form a new species?
		{
			//if so, who is the parent?
			if (SpeciesExtant == 1)
			{
				TdivA[0][1] = tcount;
				TdivA[1][0] = tcount;
				PhyloThetaA[1] = PhyloThetaA[0];
			}
			else
			{
				Parent = Big_Rand(SpeciesExtant - 1);
				TdivA[Parent][SpeciesExtant] = tcount;
				TdivA[SpeciesExtant][Parent] = tcount;
				PhyloThetaA[SpeciesExtant] = PhyloThetaA[Parent];

				for (counter = 0; counter < SpeciesExtant; counter++)
				{
					if (counter != Parent)
					{
						TdivA[counter][SpeciesExtant] = TdivA[Parent][counter];
						TdivA[SpeciesExtant][counter] = TdivA[Parent][counter];
					}

				}
			}
			SpeciesExtant++;

		}
		tcount++;
	}

		//Now do plants
		//First, generate a birth tree
		SpeciesExtant = 1;
		tcount = 0;
		TdivP[0][0] = 0;
		PhyloThetaP[0] = 10;
		while (SpeciesExtant <= PSpecies)
		{
			for (i = 0; i < SpeciesExtant; i++)
			{
				//Draw evolutionary change from a normal
				do
				{
					// choose x,y in uniform square (-1,-1) to (+1,+1) 
					u = -1 + 2 * rand() / (1.0*RAND_MAX);
					v = -1 + 2 * rand() / (1.0*RAND_MAX);
					// see if it is in the unit circle 
					r2 = u * u + v * v;
				} while (r2 > 1.0 || r2 == 0);
				scale = sqrt(-2.0 * log(r2) / r2);
				Noise = sqrt(PDriftRate) * u * scale;
				PhyloThetaP[i] = PhyloThetaP[i] + Noise;
			}
			Chance = (1.0*rand() / (1.0*RAND_MAX));
			if (Chance < bP)//do we form a new species?
			{
				//if so, who is the parent?
				if (SpeciesExtant == 1)
				{
					TdivP[0][1] = tcount;
					TdivP[1][0] = tcount;
					PhyloThetaP[1] = PhyloThetaP[0];
				}
				else
				{
					Parent = Big_Rand(SpeciesExtant - 1);
					TdivP[Parent][SpeciesExtant] = tcount;
					TdivP[SpeciesExtant][Parent] = tcount;
					PhyloThetaP[SpeciesExtant] = PhyloThetaP[Parent];

					for (counter = 0; counter < SpeciesExtant; counter++)
					{
						if (counter != Parent)
						{
							TdivP[counter][SpeciesExtant] = TdivP[Parent][counter];
							TdivP[SpeciesExtant][counter] = TdivP[Parent][counter];
						}

					}
				}
				SpeciesExtant++;

			}
			tcount++;
		}


	for (i = 0; i< ASpecies; i++)
	{
		for (j = 0; j < ASpecies; j++)
		{

			out_Phylogeny << TdivA[i][j] << ",";
		}
		out_Phylogeny << "\n";
	}
	out_Phylogeny << "\n\n";

	for (i = 0; i< PSpecies; i++)
	{
		for (j = 0; j < PSpecies; j++)
		{

			out_Phylogeny << TdivP[i][j] << ",";
		}
		out_Phylogeny << "\n";
	}

	out_Phylogeny.flush();
	

}

void Immigrate()
{
	int i,j,Immigrants;
	double m;

	m = 0.3;//set immigration rate. This will be the mean of the Poisson

	for (i = 0; i < ASpecies; i++)
	{
		Immigrants= Fish_Rand(m);
		for (j = 0; j < Immigrants; j++)
		{
			ZA[i][j + NA[i]] = ThetaA[i];
		}
		NA[i] = NA[i] + Immigrants;


	}


	for (i = 0; i < PSpecies; i++)
	{
		Immigrants = Fish_Rand(m);
		for (j = 0; j < Immigrants; j++)
		{
			ZP[i][j + NP[i]] = ThetaP[i];
		}
		NP[i] = NP[i] + Immigrants;
	//	cout << NP[i] << "\n";
		//cin >> test;
	}

}