/* C++ source code for EvolGenius v6.1 for Windows and Mac/Unix
Richard M. Kliman (2014) */

/* Floating point issue in summing of migrant genotype proportions led to erroneous warning that proportions
exceeded 1.0.  This only happened in Windows.  The code has been modified to work around the rounding error.
- Feb 16, 2015 */

#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <string>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <fstream>

using namespace std;

#define FMax 100
#define SMax 500
#define NMax 1000000
#define IMax 50000
#define GMax 9999999

/* FMax  = maximum length of a filename */
/* SMax = maximum length of general character string */
/* NMax  = maximum population size */
/* IMax  = maximum number of iterations */
/* GMax  = maximum number of generations allowed per iteration */

/* output file */
ofstream datafl;    /* .csv file with run data */
ofstream sumfl;     /* .txt file with run setting */
char onm[FMax];     /* output file prefix */
char datanm[FMax];  /* name of .csv file */
char sumnm[FMax];   /* name of .txt file */

/* general use variables */
char strng[10];	/* general use string variable */
char fchar;		/* general use character variable */
double rVal;	/* random real value b/n 0 and 1 */
int chk,chk2;	/* general use boolean */
long calcVal;   /* general use longint value */
double calcVal2;    /* general use double value */

/* program settings */
int rsd;		/* random number generator seed */
long nIter;		/* number of iterations */
long genMax;	/* maximum number of generations per iteration */
long resMax;	/* maximum number of generations for screen output resolution*/
long N[3][3];	/* number of individuals w/ each genotype */
long NTot;		/* total number of individuals in population */
double w[3][3];	/* fitnesses of genotypes */
int doMig;		/* boolean: 1=allow migration */
double migN;	/* average number of immigrants/generation */
double migP;	/* probability of migration (migN / NTot); */
double migF[3][3];	/* genotype proportions of immigrants */
double migCumF[3][3];	/* cumulative genotype proportions of immigrants */
double mutA[2];	/* 0: A>a mutation rate		1: a>A mutation rate */
double mutB[2];	/* 0: B>b mutation rate		1: b>B mutation rate */
double rec_cM;	/* recombination rate in cM (0-50) b/n A and B */
double rec;		/* recombination rate in freq (0-0.5) b/n A and B */
int mateRule;	/* mating rules:
					0=self-fert not allowed; 1=self-fert allowed; 2=enforced monogamy */
int trackCopies;/* boolean: 1=track gene copy histories */
int setAllFix[2];	/* check for allele-specific allele fixation
							element 0 = Alpha	0=any, 1=A, 2=a
							element 1 = Beta	0=any, 1=B, 2=b */
int stopAtFix;	/* boolean: 1=stop iteration upon allele and gene copy fixation*/

/* other program variables */
double trust;	/* used to determine if sum stats should be calculated */
long rep;		/* replicate being performed */
int contIter;	/* boolean: continue iteration */
int doProgAgain;/* boolean: 1=run the program again */
int doSeed;		/* boolean: 1=get random number seed */
long gen;		/* generation of iteration */
int firstRun;	/* boolean: 1=first run of program */
int fitCheck;	/* genotype differ in fitness */
int mutCheck;	/* some mutation is allowed */
int mateCheck;	/* nonrandom mating: 0 = random mating; 1 = postmating isolation; 2 = premating isolation */
int whichMating[3][3]={0,1,2,1,3,4,2,4,5};
		/* identifies which of the 6 possible matings has been chosen
					AA	Aa	aa
			  AA	0	1	2
			  Aa	1	3	4
			  aa	2	4	5	*/

int varyAInit;	/* set to 0 if population starts entirely AA or aa */
int varyBInit;	/* set to 0 if population starts entirely BB or bb */
long startPop[3][3];	/* number of individuals in starting population with
						the following genotypes:
						00=AABB		01=AABb		02=AAbb
						10=AaBB		11=AaBb		12=Aabb
						20=aaBB		21=aaBb		22=aabb */
double startF[3][3]; /* starting genotype proportions */
double matePref[6];	/* relative mating prefs (range from 0 to 1)
						matePref[0]: AA with AA
						matePref[1]: AA with Aa
						matePref[2]: AA with aa
						matePref[3]: Aa with Aa
						matePref[4]: Aa with aa
						matePref[5]: aa with aa */

/* data & results */
struct info { int AAll[2]; long AID[2]; int AHap[2];
				int BAll[2]; long BID[2]; int BHap[2];
				struct info *next; };
				/* AAll: 0=A, 1=a; BAll: 0=B, 1=b */
				/* each allele is part of an initial haplotype:
					AHap[], BHap[]: 0=AB 	1=Ab	2=aB	3=ab
					This must be recorded, because it can not be inferred
					for dihybrids */
struct info *pIND;
struct info *pNIND;
struct info *pPERMUTEIND;
double *pRand;		/* pointer to array of random numbers for permutation of starting pop */
struct info mIND[10]; /* used for immigration -- see RUNINIT() */
long *pMated;		/* used when monogamy is enforced;
						-1 if unmated; array pos of mate if mated */
double f[3][3];		/* proportions of genotypes */
int	fixedAAll;		/* 0: A fixed; 1: a fixed; 2: neither fixed */
int fixedBAll;		/* 0: B fixed; 1: b fixed; 2: neither fixed */
long totAFixed[3];	/* number of times A [0] or a [1] allele fixed; neither[2] */
long totBFixed[3];	/* number of times B [0] or b [1] allele fixed; neither[2] */
long fixedACopy;	/* ID of MRCA of fixed A copies; -1 if never fixed */
long fixedBCopy;	/* ID of MRCA of fixed B copies; -1 if never fixed */
long noFixedACopy; 	/* number of times no A copy fixed */
long noFixedBCopy;	/* number of times no B copy fixed */
long genFixedAAll;	/* generation at which A or a allele is fixed; genMax+1 if never */
long genFixedBAll;	/* generation at which B or b allele is fixed; genMax+1 if never */
long genFixedACopy;	/* generation at which A gene copy is fixed; genMax+1 if never */
long genFixedBCopy;	/* generation at which B gene copy is fixed; genMax+1 if never */
int ancAHap;		/* haplotype of MRCA for A gene at time of fixation
						0=AB	1=Ab	2=aB	3=ab*/
int ancBHap;		/* haplotype of MRCA for B gene - see above */
int iterExtinct;	/* extinction in iteration */
long nIterExtinct;	/* number of iterations leading to extinction */

/* statistics */
long sumGenFixedAAll;		/* cumulative sums of fixation generations*/
long sumGenFixedBAll;
long sumGenFixedACopy;
long sumGenFixedBCopy;
long sumGenFixedAAll2;	/* cumulative sums of (fixation generations)^2 */
long sumGenFixedBAll2;
long sumGenFixedACopy2;
long sumGenFixedBCopy2;

/* genotype arrays (for convenience) */
char aGType[3][3]={"AA","Aa","aa"};
char bGType[3][3]={"BB","Bb","bb"};
char haplo[4][3]={"AB","Ab","aB","ab"};

/* prototypes */
void GIVEUP(void);		/* halts program gracefully */
void SEEDRND(void);		/* assigns random number generator seed */
int RNDI(void);			/* used in random number generation	*/
double RANDREAL(void);	/* generates higher resolution random number between 0 and 1 */
int GF_INTCHK(void);	/* called to see that keyboard input is a valid integer */
int GF_REALCHK(void);	/* called to see that keyboard input is a valid real # */
void OUTFILE(void);		/* assigns output file names */
void STARTDEF(void);    /* sets default values for certain variables */
void START(void);		/* main menu for program settings */
void MEMALLOC(void);	/* allocates memory for population arrays */
void RUNINIT(void);		/* initializes variables for entire program run */
void ITERINIT(void);	/* initializes variables for a given iteration */
void SHUFFLE(void);     /* shuffles the starting population to minimize randnum biases */
void POPBUILD(void);	/* builds the population for generation 0 */
void MUTATE(void);		/* mutation algorithm */
void FREQCALC(void);	/* calculates proportions */
void EXCHANGE(void);	/* replaces parent generation with offspring */
void ADDNEW(void);		/* adds results to variables for summary stats */
void OUTWRITE(void);	/* writes iteration results to output file */
void CONTCHK(void);		/* checks to see if iteration can continue */
void COALESCE(void);	/* checks for gene copy fixation if trackCopies=1*/
void ALCOMP(void);		/* checks for allele fixation */
void GENBUILD(void);	/* builds new generation */
void AGAIN(void);		/* prompts user to run program again */
void SUMSTATS(void);	/* performs summary statistics calculations */


void GIVEUP(void)
{
	cout <<"\nEvolGenius 6.1 halted."<<endl;
	exit(0);
}

void SEEDRND(void)
{
	srand(rsd);
}

int RNDI(void) /* see RAN0() in Numerical Recipes in C; modified to allow for variable RAND_MAX */
{
	double static y,v[98], *pv;
	double dum;
	int j;
	double fCalc;

	if(doSeed) /* the first call to rndi (when doseed =1) fills the array */
	{
		pv=v;
		for(j=1;j<98;j++) dum=rand();
		for(j=1;j<98;j++) v[j]=rand();
		y=rand();
		doSeed=0;
	}
	j=(int)1+(97*y/RAND_MAX);
	y=*(pv+j);	/* this assigns a "random" random number from the array */
	*(pv+j)=rand(); /* this refills the array */
	/* need to ensure that random integer does not exceed 32767; RAND_MAX may be higher */
	fCalc=(double)y/RAND_MAX;
	y=(int)32768*fCalc;
	/* If RAND_MAX is >>32768 (e.g., as it is using the g++ compiler in Mac/Unix), y can only be 32768 when random
        number is RAND_MAX, which should be very rare.  Therefore, this result can be pooled with "real" 32767's.
        If RAND_MAX is 32767 (e.g., as it is in Visual C++ for Win), y can be 32766 or 32768 with equal probability,
        but not 32767.  Therefore, 32768 can be converted to 32767.
        NOTE: there may be values of RANDMAX for which this approach is unsuitable!!! */
	if(y==32768) y=32767;
	return((int)y);
}

double RANDREAL(void)
{
	int a;
	long checkRand[2];
	long calcVal;
	double calcVal2;

    for(a=0;a<2;a++) checkRand[a]=RNDI();
    calcVal=32768*checkRand[0]+checkRand[1];
    calcVal2=(double)calcVal/1073741824;
	return((double)calcVal2);
}

int GF_INTCHK(char instrng[SMax])
{
	int a;

	for(a=0;a<strlen(instrng);a++) if(!isdigit(instrng[a])) return(0);
	return(1);
}

int GF_REALCHK (char instrng[SMax])
{
	int period,a;

	period=0;
	for(a=0;a<std::strlen(instrng);a++)
	{
		if(!isdigit(instrng[a]))	/* found something other than a digit */
		{
			if(instrng[a]=='.')		/* found a period */
			{
				if(period) return(0);	/* found a second period */
				else period=1;
			}
			else return(0);		/* character was neither a digit nor a period */
		}
	}
	return(1);
}

void OUTFILE(void)
{
	chk=0;
	while(!chk)
	{
		cout<<"Prefix of output files (.csv and .txt suffixes added): "<<flush;
		cin>>strng;
		strcpy(onm,strng);  /* assign filename prefix */
		strcpy(datanm,onm); /* copy prefix to .csv filename */
		strcpy(sumnm,onm);  /* copy prefix to .txt filename */
		strcat(datanm,".csv\0");
		strcat(sumnm,".txt\0");
		datafl.clear();
		datafl.open(datanm,ios::out);
		if(datafl.fail())
		{
			cout<<"\aTHE FILE "<<datanm<<" MAY ALREADY EXIST."<<endl
				<<"TYPE ""Q"" <enter> TO QUIT "<<endl
				<<"OR ""N"" <enter> FOR NEW FILENAME.> "<<flush;
			fchar='.';
			while(!(fchar=='Q' || fchar=='N'))
			{
				cin>>strng;
				fchar=toupper(strng[0]);
			}
			if(fchar=='Q') GIVEUP();
			datafl.clear();
		}
		else chk=1;
        sumfl.clear();
		sumfl.open(sumnm,ios::out);
		if(sumfl.fail())
		{
			cout<<"\aTHE FILE "<<sumnm<<" MAY ALREADY EXIST."<<endl
				<<"TYPE ""Q"" <enter> TO QUIT "<<endl
				<<"OR ""N"" <enter> FOR NEW FILENAME.> "<<flush;
			fchar='.';
			while(!(fchar=='Q' || fchar=='N'))
			{
				cin>>strng;
				fchar=toupper(strng[0]);
			}
			if(fchar=='Q') GIVEUP();
			sumfl.clear();
		}
		else chk=1;
	}
	datafl.close();
	sumfl.close();
}

void STARTDEF(void)
{
	int a,b;

	genMax=30000;		/* default: 30000 max generations */
	startPop[0][0]=10;	/* default: 1:2:1:2:4:2:1:2:1 genotype ratio */
	startPop[0][1]=20;
	startPop[0][2]=10;
	startPop[1][0]=20;
	startPop[1][1]=40;
	startPop[1][2]=20;
	startPop[2][0]=10;
	startPop[2][1]=20;
	startPop[2][2]=10;
	NTot=0;
	for(a=0;a<3;a++) for(b=0;b<3;b++)
	{
		NTot+=startPop[a][b];	/* calculate pop size */
		w[a][b]=1.0;			/* set all relative fitnesses to 1.0 */
	}

	/* calculate starting genotype proportions */
	for(a=0;a<3;a++) for(b=0;b<3;b++) startF[a][b]=(double)startPop[a][b]/NTot;

	for(a=0;a<2;a++)
	{
		mutA[a]=0.0;	/* default: no mutation */
		mutB[a]=0.0;
	}
	nIter=1;		/* default: 1 iteration */
	mateRule=1;		/* default: self-fertilization is allowed */
	for(a=0;a<6;a++) matePref[a]=1.0;	/* default: random mating */
	trackCopies=1;	/* default: track gene copy history */
	setAllFix[0]=0;	/* default: can fix either Alpha allele */
	setAllFix[1]=0;	/* default: can fix either Beta allele */
	stopAtFix=1;	/* default: stop iteration upon fixation */
}

void START(void)
{
	int a,b,c;
	int StepP=0;	/* to change starting population */
	int StepM=0;	/* to change maximum number of generations */
	int StepW=0;	/* to change fitness values */
	int StepL=0;	/* to change linkage map distance */
	int StepD=0;	/* to change migration settings */
	int Step1=0;	/* to change mutation rate of A to a */
	int Step2=0;	/* to change mutation rate of a to A */
	int Step3=0;	/* to change mutation rate of B to b */
	int Step4=0;	/* to change mutation rate of b to B */
	int StepI=0;	/* to change number of iterations */
	int StepF=0;	/* to toggle allow self-fertilization */
	int StepN=0;	/* to change relative mating preferences */
	int StepG=0;	/* to toggle tracking of gene copies */
	int StepA=0;	/* to toggle allele fixation setting */
	int StepV=0;	/* to toggle stop at fixation */
	int nfail;
	int startAGtype[3];	/* number of individuals w/ AA, Aa and aa gtypes */
	char chng='0';	/* menu option chosen */
	long tVal;		/* typed value at prompt (later converted to number) */
	double wMax,mateMax;	/* used to calculate relative values */
	double migCum;	/* used to calculate migCumF */

	/* check to be sure simulation can run given initial settings*/
	int fertChk;	/* 1 if enough fertile individuals to continue */
	int mateChk;	/* 1 if mating is possible */
	int migChk;		/* 1 if more immigrants than total population size */

	OUTFILE();
	chk=1;
	while(chk)
	{
		cout<<"Random number generator seed [integer, 1 - 32767]: "<<flush;
		cin>>strng;
		if(GF_INTCHK(strng))
		{
			tVal=atol(strng);
			if(tVal>0 && tVal<32767)
			{
				chk=0;
				rsd=(int)tVal;
			}
		}
		if(chk)	cout<<"\a\nBad seed.  Try again."<<endl;
	}

	do
	{
		for(a=0;a<25;a++) cout<<endl;
		while(StepP)
		{
			cout<<"Maximum population size is "<<NMax<<".\n"<<endl;
			cout<<"Enter Values for the Starting Population:"<<endl;
			for(a=0;a<3;a++) for(b=0;b<3;b++)
			{
				chk=0;
				do
				{
					cout<<"Number of "<<aGType[a]<<bGType[b]<<" individuals: "<<flush;
					cin>>strng;

					if(GF_INTCHK(strng))
					{
						tVal=atol(strng);
						if(tVal>=0)
						{
							chk=1;
							startPop[a][b]=(long)tVal;
						}
					}
					if(!chk) cout<<"\a\nBad number.  Try again."<<endl;
				} while(!chk);
			}
			NTot=0;
			for(a=0;a<3;a++) for(b=0;b<3;b++) NTot+=startPop[a][b];
			if(NTot>0 && NTot<=NMax) StepP=0;
			else
			{
				if(NTot==0) cout<<"\a\nPopulation size is zero."<<endl;
				if(NTot>NMax) cout<<"\a\nPopulation size exceeds "<<NMax<<"."<<endl;
			}
			if(!StepP)
				for(a=0;a<3;a++) for(b=0;b<3;b++) startF[a][b]=(double)startPop[a][b]/NTot;
		}

		while(StepM)
		{
			cout<<"Maximum number of generations/iteration [1 - "<<GMax<<"]: "<<flush;
			cin>>strng;
			if(GF_INTCHK(strng))
			{
				tVal=atol(strng);
				if(tVal>0 && tVal<=GMax)
				{
					StepM=0;
					genMax=(long)tVal;
				}
			}
			if(StepM) cout<<"\a\nBad number.  Try again."<<endl;
		}

		if(StepW)
		{
			for(a=0;a<3;a++) for(b=0;b<3;b++)
			{
				chk=0;
				do
				{
					cout<<"Fitness of "<<aGType[a]<<bGType[b]<<flush;
					cout<<" genotype [0.0 - 1.0]: "<<flush;
					cin>>strng;
					if(GF_REALCHK(strng))
					{
						w[a][b]=atof(strng);
						if(w[a][b]>=0.0 && w[a][b]<=1.0) chk=1;
					}
					if(!chk) cout<<"\a\nBad value.  Try again."<<endl;
				} while (!chk);
			}
			wMax=0.0;
			for(a=0;a<3;a++) for(b=0;b<3;b++) if(w[a][b]>wMax) wMax=w[a][b];
			for(a=0;a<3;a++) for(b=0;b<3;b++) w[a][b]=w[a][b]/wMax;
			StepW=0;
		}

		while(StepL)
		{
			cout<<"You must specify the map distance between"<<endl;
			cout<<"genes A and B."<<endl;
			cout<<"  A map distance of 0.0 cM allows no"<<endl;
			cout<<"  crossing-over between the two genes.\n"<<endl;
			cout<<"  A map distance of 50 cM allows for"<<endl;
			cout<<"  independent assortment of the two genes.\n"<<endl;
			cout<<"Map distance [0.0 - 50.0 cM]: "<<flush;
			cin>>strng;
			if(GF_REALCHK(strng))
			{
				rec_cM=atof(strng);
				if(rec_cM>=0.0 && rec_cM<=50.0)
				{
					rec=rec_cM/100;
					StepL=0;
				}
			}
			if(StepL) cout<<"\a\nBad value.  Try again."<<endl;
		}

		if(StepF)
		{
			mateRule++;
			if(mateRule==3) mateRule=0;
			StepF=0;
		}

		while(StepD)
		{
			if(!doMig)
			{
				chk2=0; /* used to make sure freqs add up to 1.0 */
				while(!chk2)
				{
					for(a=0;a<3;a++) for(b=0;b<3;b++)
					{
						chk=0;
						do
						{
							if(a==2 && b==2) /* calculate p(aabb) */
							{
								migF[a][b]=(double)1.0-(migF[0][0]+migF[0][1]+
									migF[0][2]+migF[1][0]+migF[1][1]+
									migF[1][2]+migF[2][0]+migF[2][1]);
                                if(migF[2][2]>-0.000000000000001 && migF[2][2]<0.0000000000001) migF[2][2]=0;
								chk=1;
							}
							else /* enter other proportions */
							{
								cout<<"Proportion of "<<aGType[a]<<bGType[b]<<flush;
								cout<<" genotype in immigrants [0.0 - 1.0]: "<<flush;
								cin>>strng;
								if(GF_REALCHK(strng))
                                {
									migF[a][b]=(double)(atof(strng));
									if(migF[a][b]>=0.0 && migF[a][b]<=1.0)
										chk=1;
								}
								if(!chk) cout<<"\a\nBad value.  Try again."<<endl;
							}
						} while (!chk);
					}
					if(migF[2][2]<0.0)
					{
						chk2=0;
						cout<<"\a\nGenotype proportions add up to a number"<<endl;
						cout<<"greater than 1.0.  Try again.\n"<<endl;
					}
					else chk2=1;
				}

				migCum=0.0;
				for(a=0;a<3;a++) for(b=0;b<3;b++)
				{
					migCum=migCum+migF[a][b];
					migCumF[a][b]=migCum;
				}

				cout<<"\nImmigrant genotype proportions:"<<endl;
				for(a=0;a<3;a++) for(b=0;b<3;b++)
				{
					cout<<"p("<<aGType[a]<<bGType[b]<<")="<<flush;
					cout<<setprecision(4)<<setiosflags(ios::fixed)<<migF[a][b]<<endl;
				}

				chk=0;
				while(!chk)
				{
					cout<<"\n\nAverage number of immigrants per generation."<<endl;
					cout<<"This number does not have to be an integer."<<endl;
					cout<<"It must be greater than zero and less than the"<<endl;
					cout<<"population size ("<<NTot<<")."<<endl;
					cout<<"Immigrants/generation [>0.0 - <"<<NTot<<"]: "<<flush;
					cin>>strng;
					if(GF_REALCHK(strng))
                    {
						migN=atof(strng);
						if(migN>0.0 && migN<NTot) chk=1;
					}
					if(!chk) cout<<"\a\nBad number.  Try again."<<endl;
				}
				doMig=1;
				StepD=0;
				if(trackCopies)
				{
					trackCopies=0;
					cout<<"\nNote: You can not track gene copies when"<<endl;
					cout<<"migration is allowed.  The program settings"<<endl;
					cout<<"have been changed to reflect this.\n"<<endl;
					cout<<"Type ""C""<enter> to continue "<<flush;
					fchar='.';
					while(fchar!='C')
					{
						cin>>strng;
						fchar=toupper(strng[0]);
					}
				}
			}
			else
			{
				doMig=0;
				migN=0.0;
				migP=0.0;
				StepD=0;
			}
		}

		while(Step1)
		{
			cout<<"Mutation rate of A to a [0.0 - 1.0]: "<<flush;
			cin>>strng;
			if(GF_REALCHK(strng))
			{
				mutA[0]=atof(strng);
				if(mutA[0]>=0.0 && mutA[0]<=1.0) Step1=0;
			}
			if(Step1) cout<<"\a\nBad value.  Try again."<<endl;
		}

		while(Step2)
		{
			cout<<"Mutation rate of a to A [0.0 - 1.0]: "<<flush;
			cin>>strng;
			if(GF_REALCHK(strng))
			{
				mutA[1]=atof(strng);
				if(mutA[1]>=0.0 && mutA[1]<=1.0) Step2=0;
			}
			if(Step2) cout<<"\a\nBad value.  Try again."<<endl;
		}

		while(Step3)
		{
			cout<<"Mutation rate of B to b [0.0 - 1.0]: "<<flush;
			cin>>strng;
			if(GF_REALCHK(strng))
			{
				mutB[0]=atof(strng);
				if(mutB[0]>=0.0 && mutB[0]<=1.0) Step3=0;
			}
			if(Step3) cout<<"\a\nBad value.  Try again."<<endl;
		}

		while(Step4)
		{
			cout<<"Mutation rate of b to B [0.0 - 1.0]: "<<flush;
			cin>>strng;
			if(GF_REALCHK(strng))
			{
				mutB[1]=atof(strng);
				if(mutB[1]>=0.0 && mutB[1]<=1.0) Step4=0;
			}
			if(Step4) cout<<"\a\nBad value.  Try again."<<endl;
		}

		while(StepI)
		{
			cout<<"Number of iterations using the same settings "<<flush;
			cout<<"[1 - "<<IMax<<"]: "<<flush;
			cin>>strng;
			if(GF_INTCHK(strng))
			{
				tVal=atol(strng);
				if(tVal>0 && tVal<IMax+1)
				{
					StepI=0;
					nIter=(long)tVal;
				}
			}
			if(StepI) cout<<"\a\nBad number.  Try again."<<endl;
		}

		if(StepN)
		{
			cout<<"Set relative mating preferences for the"<<endl;
			cout<<"six possible matings involving gene Alpha."<<endl;
			cout<<"Relative mating preferences will be scaled"<<endl;
			cout<<"to a maximum of 1.0.  Negative values are"<<endl;
			cout<<"not allowed.\n"<<endl;
			cout<<"To prevent a given mating, set the relative"<<endl;
			cout<<"mating preference to 0.0.\n"<<endl;
			for(a=0;a<6;a++)
			{
				chk=0;
				do
				{
					if(a==0) cout<<"AA with AA [0.0 - 1.0]: "<<flush;
					else if(a==1) cout<<"AA with Aa [0.0 - 1.0]: "<<flush;
					else if(a==2) cout<<"AA with aa [0.0 - 1.0]: "<<flush;
					else if(a==3) cout<<"Aa with Aa [0.0 - 1.0]: "<<flush;
					else if(a==4) cout<<"Aa with aa [0.0 - 1.0]: "<<flush;
					else cout<<"aa with aa [0.0 - 1.0]: "<<flush;
					cin>>strng;
					if(GF_REALCHK(strng))
					{
						matePref[a]=atof(strng);
						if(matePref[a]>=0.0 && matePref[a]<=1.0) chk=1;
					}
					if(!chk) cout<<"\a\nBad value.  Try again."<<endl;
				} while (!chk);
			}
			mateMax=0.0;
			for(a=0;a<6;a++) if(matePref[a]>mateMax) mateMax=matePref[a];
			for(a=0;a<6;a++) matePref[a]=matePref[a]/mateMax;

            mateCheck=0; /* default: there are no mating preferences */
			for(a=0;a<6 && !mateCheck;a++) if(matePref[a]<1)
            {
                chk=0;
                do
                {
                    cout<<"\nYou have indicated that there are mating preferences."<<endl;
                    cout<<"Type '1' for postmating isolaton; type '2' for premating isolation."<<endl;
                    cin>>strng;
                    if(GF_INTCHK(strng))
                    {
                        chk=1;
                        tVal=atol(strng);
                        if(tVal==1) mateCheck=1;
                        else if(tVal==2) mateCheck=2;
                        else chk=0;
                    }
                    if(!chk) cout<<"Invalid choice.  Try again."<<endl;
                } while (!chk);
            }
			StepN=0;
		}

		if(StepG)
		{
			if(trackCopies) trackCopies=0;
			else
			{
				if(doMig)
				{
					cout<<"Note: You can not track gene copies if"<<endl;
					cout<<"migration is allowed.\n"<<endl;
					cout<<"Type ""C""<enter> to continue. "<<flush;
					fchar='.';
					while(fchar!='C')
					{
						cin>>strng;
						fchar=toupper(strng[0]);
					}
				}
				else trackCopies=1;
			}
			StepG=0;
		}

		if(StepA)
		{
			cout<<"If you do not care which Alpha allele fixes, type 0."<<endl;
			cout<<"If you want to know when the A allele fixes, type 1."<<endl;
			cout<<"If you want to know when the a allele fixes, type 2."<<endl;
			chk=1;
			while(chk)
			{
				cout<<"Enter 0, 1 or 2: "<<flush;
				cin>>strng;
				if(GF_INTCHK(strng))
				{
					tVal=atol(strng);
					if(tVal>=0 && tVal<=2)
					{
						chk=0;
						setAllFix[0]=(int)tVal;
					}
				}
				if(chk)	cout<<"\a\nBad value.  Try again."<<endl;
			}
			cout<<"\nIf you do not care which Beta allele fixes, type 0."<<endl;
			cout<<"If you want to know when the B allele fixes, type 1."<<endl;
			cout<<"If you want to know when the b allele fixes, type 2."<<endl;
			chk=1;
			while(chk)
			{
				cout<<"Enter 0, 1 or 2: "<<flush;
				cin>>strng;
				if(GF_INTCHK(strng))
				{
					tVal=atol(strng);
					if(tVal>=0 && tVal<=2)
					{
						chk=0;
						setAllFix[1]=(int)tVal;
					}
				}
				if(chk)	cout<<"\a\nBad value.  Try again."<<endl;
			}
			if(setAllFix[0]>0 || setAllFix[1]>0)
			{
				cout<<"\nWARNING:"<<endl;
				cout<<"You have chosen to track fixation of a specific allele."<<endl;
				cout<<"In the absence of appropriate mutation and/or migration,"<<endl;
				cout<<"fixation of that allele may not occur before the specified"<<endl;
				cout<<"maximum number of generations.\n"<<endl;
				cout<<"Type ""C""<enter> to continue. "<<endl;
				fchar='.';
				while(fchar!='C')
				{
					cin>>strng;
					fchar=toupper(strng[0]);
				}
			}
			StepA=0;
		}

		if(StepV)
		{
			if(stopAtFix) stopAtFix=0;
			else stopAtFix=1;
			StepV=0;
		}

		for(a=0;a<25;a++) cout<<endl;
		cout<<"Population Size: "<<setw(7)<<NTot<<flush;
		cout<<"          MUTATION RATES"<<endl;

		cout<<"Type P to change POPULATION "<<flush;
		cout<<"      1 A>a: "<<flush;
		cout<<setprecision(8)<<setiosflags(ios::fixed)<<mutA[0]<<flush;
		cout<<"    2 a>A: "<<flush;
		cout<<setprecision(8)<<setiosflags(ios::fixed)<<mutA[1]<<endl;

		cout<<"Type W to change FITNESS    "<<flush;
		cout<<"      3 B>b: "<<flush;
		cout<<setprecision(8)<<setiosflags(ios::fixed)<<mutB[0]<<flush;
		cout<<"    4 b>B: "<<flush;
		cout<<setprecision(8)<<setiosflags(ios::fixed)<<mutB[1]<<endl;

		cout<<"           P            W       "<<endl;

		cout<<"AABB "<<setw(7)<<startPop[0][0]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[0][0]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[0][0]<<flush;
        if(mateCheck==0) cout<<"       N MATING PREFERENCES (none)"<<endl;
        else if(mateCheck==1)  cout<<"       N MATING PREFERENCES (postmating)"<<endl;
        else  cout<<"       N MATING PREFERENCES (premating)"<<endl;
		cout<<"AABb "<<setw(7)<<startPop[0][1]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[0][1]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[0][1]<<flush;
		cout<<"       AA:AA "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[0]<<"  "<<flush;
		cout<<"    AA:Aa "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[1]<<endl;

		cout<<"AAbb "<<setw(7)<<startPop[0][2]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[0][2]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[0][2]<<flush;
		cout<<"       AA:aa "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[2]<<"  "<<flush;
		cout<<"    Aa:Aa "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[3]<<endl;

		cout<<"AaBB "<<setw(7)<<startPop[1][0]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[1][0]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[1][0]<<flush;
		cout<<"       Aa:aa "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[4]<<"  "<<flush;
		cout<<"    aa:aa "<<setprecision(7)<<setiosflags(ios::fixed)<<matePref[5]<<endl;

		cout<<"AaBb "<<setw(7)<<startPop[1][1]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[1][1]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[1][1]<<endl;

		cout<<"Aabb "<<setw(7)<<startPop[1][2]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[1][2]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[1][2]<<flush;
		cout<<"       L LINKAGE DISTANCE:  "<<flush;
		cout<<setprecision(5)<<setiosflags(ios::fixed)<<rec_cM<<" cM"<<endl;

		cout<<"aaBB "<<setw(7)<<startPop[2][0]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[2][0]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[2][0]<<flush;
		cout<<"       F FERTILIZATION:     "<<flush;
		if(!mateRule) cout<<"Self-fert not allowed"<<endl;
		else if(mateRule==1) cout<<"Self-fert allowed"<<endl;
		else cout<<"Monogamy enforced"<<endl;

		cout<<"aaBb "<<setw(7)<<startPop[2][1]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[2][1]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[2][1]<<flush;
		cout<<"       I ITERATIONS:        "<<nIter<<endl;


		cout<<"aabb "<<setw(7)<<startPop[2][2]<<"  "<<flush;
		cout<<setprecision(3)<<setiosflags(ios::fixed)<<startF[2][2]<<"  "<<flush;
		cout<<setprecision(4)<<setiosflags(ios::fixed)<<w[2][2]<<flush;
		cout<<"       M MAX GENERATIONS:   "<<genMax<<endl<<endl;

		cout<<"D MIGRATION: "<<flush;
		if(!doMig) cout<<"none"<<endl;
		else
		{
			cout<<"Average number of immigrants/generation: "<<flush;
			cout<<setprecision(5)<<migN<<endl;
			cout<<"    p(AABB)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][0]<<" "<<flush;
			cout<<"  p(AABb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][1]<<" "<<flush;
			cout<<"  p(AAbb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][2]<<endl;
			cout<<"    p(AaBB)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][0]<<" "<<flush;
			cout<<"  p(AaBb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][1]<<" "<<flush;
			cout<<"  p(Aabb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][2]<<endl;
			cout<<"    p(aaBB)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][0]<<" "<<flush;
			cout<<"  p(aaBb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][1]<<" "<<flush;
			cout<<"  p(aabb)="<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][2]<<endl;
		}

		cout<<"G TRACK GENE COPY HISTORY? "<<flush;
		if(trackCopies) cout<<"Yes"<<endl;
		else cout<<"No"<<endl;
		cout<<"V STOP ITERATION AT FIXATION? "<<flush;
		if(stopAtFix) cout<<"Yes"<<endl;
		else cout<<"No"<<endl;
		cout<<"A ALLELE FIXATION: Alpha="<<flush;
		if(setAllFix[0]==0) cout<<"any | "<<flush;
		else if(setAllFix[0]==1) cout<<"A  | "<<flush;
		else cout<<"a   | "<<flush;
		cout<<"Beta="<<flush;
		if(setAllFix[1]==0) cout<<"any"<<endl;
		else if(setAllFix[1]==1) cout<<"B"<<endl;
		else cout<<"b"<<endl;

		trust=log(nIter)+2*(log(NTot));
		cout<<"\nType the letter or number <enter> for the field you wish "<<flush;
		cout<<"to change."<<endl;
		cout<<"Type ""C""<enter> to run the simulation(s).  Type ""Q""<enter> to quit. "<<flush;
		nfail=0;
		chk=0;
		while(!chk)
		{
			cin>>strng;
			chng=toupper(strng[0]);
			if(chng=='Q') GIVEUP();
			switch(chng)
			{
				case 'P': case 'M': case 'W': case 'L': case 'F':
				case '1': case '2': case '3': case '4': case 'I':
				case 'N': case 'G': case 'D': case 'A': case 'V':
				case 'C':	{chk=1; break;}
				default:
					nfail++;
					if(nfail==1) cout<<"\aTRY AGAIN "<<flush;
					else cout<<"\a> ";
			}
		}
		if(chng=='P') StepP=1; if(chng=='M') StepM=1; if(chng=='W') StepW=1;
		if(chng=='L') StepL=1; if(chng=='1') Step1=1; if(chng=='2') Step2=1;
		if(chng=='3') Step3=1; if(chng=='4') Step4=1; if(chng=='I') StepI=1;
		if(chng=='F') StepF=1; if(chng=='N') StepN=1; if(chng=='G') StepG=1;
		if(chng=='D') StepD=1; if(chng=='A') StepA=1; if(chng=='V') StepV=1;

		fertChk=1;
		c=0;	/* number of fertile individuals */
		for(a=0;a<3;a++) for(b=0;b<3;b++) if(w[a][b]>0) c+=startPop[a][b];
		if(mateRule && !c) /* either self-fert allowed or enforced monogamy, with no fertile inds */
		{
			for(a=0;a<25;a++) cout<<endl;
			cout<<"\n\nThe starting population has no fertile individuals.\n"<<endl;
			fertChk=0;
		}
		if((!mateRule || mateRule==2) && c<2) /* either no self-fert or enforced monogamy, with fewer than 2 fertile inds */
		{
			for(a=0;a<25;a++) cout<<endl;
			cout<<"\n\nThe starting population has fewer than two fertile individuals.\n"<<endl;
			fertChk=0;
		}
		if(!fertChk)
		{
			cout<<"Type ""C""<enter> to continue. "<<endl;
			fchar='.';
			while(fchar!='C')
			{
				cin>>strng;
				fchar=toupper(strng[0]);
			}
		}

		mateChk=0;
		startAGtype[0]=startPop[0][0]+startPop[0][1]+startPop[0][2];
		startAGtype[1]=startPop[1][0]+startPop[1][1]+startPop[1][2];
		startAGtype[2]=startPop[2][0]+startPop[2][1]+startPop[2][2];
		if(startAGtype[0]+mateRule>1 && matePref[0]>0.0) mateChk=1;
		if(!mateChk) if(startAGtype[1]+mateRule>1 && matePref[3]>0.0) mateChk=1;
		if(!mateChk) if(startAGtype[2]+mateRule>1 && matePref[5]>0.0) mateChk=1;
		if(!mateChk) if(startAGtype[0] && startAGtype[1] && matePref[1]>0.0) mateChk=1;
		if(!mateChk) if(startAGtype[0] && startAGtype[2] && matePref[2]>0.0) mateChk=1;
		if(!mateChk) if(startAGtype[1] && startAGtype[2] && matePref[4]>0.0) mateChk=1;
		if(!mateChk)
		{
			cout<<"\a\n\n\nThere are no allowable matings."<<endl;
			cout<<"Change either the starting population or"<<endl;
			cout<<"the mating preferences.\n"<<endl;
			cout<<"Type ""C""<enter> to continue. "<<endl;
			fchar='.';
			while(fchar!='C')
			{
				cin>>strng;
				fchar=toupper(strng[0]);
			}
		}

		migChk=1;
		if(migN>NTot)
		{
			migChk=0;
			cout<<"\a\n\nYou have changed your initial population size."<<endl;
			cout<<"You now have more immigrants/generation than allowed."<<endl;
			cout<<"Change either the initial population or the ";
			cout<<"migration settings.\n"<<endl;
			cout<<"Type ""C""<enter> to continue. "<<endl;
			fchar='.';
			while(fchar!='C')
			{
				cin>>strng;
				fchar=toupper(strng[0]);
			}
		}

	} while(chng!='C' || !fertChk || !mateChk || !migChk);

	SEEDRND();
	fitCheck=0;
	for(a=0;a<3;a++) for(b=0;b<3;b++) if(w[a][b]<1.0) fitCheck=1;
	mutCheck=0;
	if(mutA[0]>0.0 || mutA[1]>0.0 || mutB[0]>0.0 || mutB[1]>0.0) mutCheck=1;
	datafl.open(datanm,ios::out);
	sumfl.open(sumnm,ios::out);
	sumfl<<"EvolGenius v6.1 by R.M. Kliman (2014)"<<endl;
}

void MEMALLOC(void)
{
	if((pIND=(struct info *)calloc(NTot,sizeof(struct info)))==NULL)
	{
		cout<<"\a\nNOT ENOUGH MEMORY."<<endl;
		GIVEUP();
	}
	if((pNIND=(struct info *)calloc(NTot,sizeof(struct info)))==NULL)
	{
		cout<<"\a\nNOT ENOUGH MEMORY."<<endl;
		GIVEUP();
	}
	if((pPERMUTEIND=(struct info *)calloc(NTot,sizeof(struct info)))==NULL)
	{
		cout<<"\a\nNOT ENOUGH MEMORY."<<endl;
		GIVEUP();
	}
	if((pRand=(double *)calloc(NTot,sizeof(double)))==NULL)
	{
		cout<<"\a\nNOT ENOUGH MEMORY."<<endl;
		GIVEUP();
	}

	if((pMated=(long *)calloc(NTot,sizeof(long)))==NULL)
	{
		cout<<"\a\nNOT ENOUGH MEMORY."<<endl;
		GIVEUP();
	}
}

void RUNINIT(void)
{
	int a,b;

	migP=(double)migN/NTot;
	sumGenFixedAAll=sumGenFixedBAll=0;
	sumGenFixedACopy=sumGenFixedBCopy=0;
	sumGenFixedAAll2=sumGenFixedBAll2=0;
	sumGenFixedACopy2=sumGenFixedBCopy2=0;
	totAFixed[0]=totAFixed[1]=totAFixed[2]=0;
	totBFixed[0]=totBFixed[1]=totBFixed[2]=0;
	noFixedACopy=noFixedBCopy=0;
	nIterExtinct=0;
	/* AABB */
	mIND[0].AAll[0]=0; mIND[0].AAll[1]=0; mIND[0].BAll[0]=0; mIND[0].BAll[1]=0;
	/* AABb */
	mIND[1].AAll[0]=0; mIND[1].AAll[1]=0; mIND[1].BAll[0]=0; mIND[1].BAll[1]=1;
	/* AAbb */
	mIND[2].AAll[0]=0; mIND[2].AAll[1]=0; mIND[2].BAll[0]=1; mIND[2].BAll[1]=1;
	/* AaBB */
	mIND[3].AAll[0]=0; mIND[3].AAll[1]=1; mIND[3].BAll[0]=0; mIND[3].BAll[1]=0;
	/* AB/ab */
	mIND[4].AAll[0]=0; mIND[4].AAll[1]=1; mIND[4].BAll[0]=0; mIND[4].BAll[1]=1;
	/* Ab/AB */
	mIND[5].AAll[0]=0; mIND[5].AAll[1]=1; mIND[5].BAll[0]=1; mIND[5].BAll[1]=0;
	/* Aabb */
	mIND[6].AAll[0]=0; mIND[6].AAll[1]=1; mIND[6].BAll[0]=1; mIND[6].BAll[1]=1;
	/* aaBB */
	mIND[7].AAll[0]=1; mIND[7].AAll[1]=1; mIND[7].BAll[0]=0; mIND[7].BAll[1]=0;
	/* aaBb */
	mIND[8].AAll[0]=1; mIND[8].AAll[1]=1; mIND[8].BAll[0]=0; mIND[8].BAll[1]=1;
	/* aabb */
	mIND[9].AAll[0]=1; mIND[9].AAll[1]=1; mIND[9].BAll[0]=1; mIND[9].BAll[1]=1;
	for(a=0;a<10;a++) for(b=0;b<2;b++)
	{
		mIND[a].AID[b]=-1;
		mIND[a].AHap[b]=-1;
		mIND[a].BID[b]=-1;
		mIND[a].BHap[b]=-1;
	}
}

void SHUFFLE(void)
{
    /* This is a routine to shuffle the members of NIND as protection against possible biases
        in the random number generator.  While all subsequent generations are filled from random
        matings, the starting generation is filled in order of genotype.  Thus, minor biases could
        influence the fate of rare genotypes.  By shuffling the starting generation, the biases
        would be uncorrelated with genotype. */

    long a,b;
    long lowPos;
   	double lowRand;		/* low remaining value in pRAND array during permutation */

   	for(a=0;a<NTot;a++) pRand[a]=RANDREAL();	/* fill array of random numbers for permutation */
	lowPos=0;
	for(a=0;a<NTot;a++) /* fill permuted array of individuals */
	{
		lowRand=2;  /* a value that must exceed all original values in random number array */
		for(b=0;b<NTot;b++) /* find lowest remaining value in array of random numbers */
		{
			if(pRand[b]<lowRand)
			{
				lowPos=b;
				lowRand=pRand[b];
			}
		}
        *(pPERMUTEIND+a)=*(pNIND+b); /* assigns individual b in NIND to individual a in PERMUTEIND */
		pRand[lowPos]=3;    /* convert lowest remaining value in array of random numbers to a value >2 */
	}
    for(a=0;a<NTot;a++) *(pNIND+a)=*(pPERMUTEIND+a); /* re-fills NIND to feed into simulation */
}

void POPBUILD(void)
{
	long a,b;

	for(a=0;a<3;a++) for(b=0;b<3;b++)
	{
		N[a][b]=startPop[a][b];
		f[a][b]=startF[a][b];
	}
	a=0;
	for(b=0;b<N[0][0];b++) /* AABB */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=0;
		pIND[a].BAll[0]=0; pIND[a].BAll[1]=0;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=0; pIND[a].BHap[0]=0; /* haplotype = AB */
		pIND[a].AHap[1]=0; pIND[a].BHap[1]=0; /* haplotype = AB */
		a++;
	}
	for(b=0;b<N[0][1];b++) /* AABb */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=0;
		pIND[a].BAll[0]=0; pIND[a].BAll[1]=1;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=0; /* haplotype = AB */
		pIND[a].AHap[1]=pIND[a].BHap[1]=1; /* haplotype = Ab */
		a++;
	}
	for(b=0;b<N[0][2];b++) /* AAbb */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=0;
		pIND[a].BAll[0]=1; pIND[a].BAll[1]=1;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=1; /* haplotype = Ab */
		pIND[a].AHap[1]=pIND[a].BHap[1]=1; /* haplotype = Ab */
		a++;
	}
	for(b=0;b<N[1][0];b++) /* AaBB */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=1;
		pIND[a].BAll[0]=0; pIND[a].BAll[1]=0;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=0; /* haplotype = AB */
		pIND[a].AHap[1]=pIND[a].BHap[1]=2; /* haplotype = aB */
		a++;
	}
	for(b=0;b<N[1][1];b++) /* AaBb */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=1;
		if(RANDREAL()<0.5)
		{
			pIND[a].BAll[0]=0; pIND[a].BAll[1]=1;
			pIND[a].AHap[0]=pIND[a].BHap[0]=0; /* haplotype = AB */
			pIND[a].AHap[1]=pIND[a].BHap[1]=3; /* haplotype = ab */
		}
		else
		{
			pIND[a].BAll[0]=1; pIND[a].BAll[1]=0;
			pIND[a].AHap[0]=pIND[a].BHap[0]=1; /* haplotype = Ab */
			pIND[a].AHap[1]=pIND[a].BHap[1]=2; /* haplotype = aB */
		}
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		a++;
	}
	for(b=0;b<N[1][2];b++) /* Aabb */
	{
		pIND[a].AAll[0]=0; pIND[a].AAll[1]=1;
		pIND[a].BAll[0]=1; pIND[a].BAll[1]=1;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=1; /* haplotype = Ab */
		pIND[a].AHap[1]=pIND[a].BHap[1]=3; /* haplotype = ab */
		a++;
	}
	for(b=0;b<N[2][0];b++) /* aaBB */
	{
		pIND[a].AAll[0]=1; pIND[a].AAll[1]=1;
		pIND[a].BAll[0]=0; pIND[a].BAll[1]=0;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=2; /* haplotype = aB */
		pIND[a].AHap[1]=pIND[a].BHap[1]=2; /* haplotype = aB */
		a++;
	}
	for(b=0;b<N[2][1];b++) /* aaBb */
	{
		pIND[a].AAll[0]=1; pIND[a].AAll[1]=1;
		pIND[a].BAll[0]=0; pIND[a].BAll[1]=1;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=2; /* haplotype = aB */
		pIND[a].AHap[1]=pIND[a].BHap[1]=3; /* haplotype = ab */
		a++;
	}
	for(b=0;b<N[2][2];b++) /* aabb */
	{
		pIND[a].AAll[0]=1; pIND[a].AAll[1]=1;
		pIND[a].BAll[0]=1; pIND[a].BAll[1]=1;
		pIND[a].AID[0]=pIND[a].BID[0]=2*a;
		pIND[a].AID[1]=pIND[a].BID[1]=2*a+1;
		pIND[a].AHap[0]=pIND[a].BHap[0]=3; /* haplotype = ab */
		pIND[a].AHap[1]=pIND[a].BHap[1]=3; /* haplotype = ab */
		a++;
	}
}

void ITERINIT(void)
{
	int a;
	long b;

	contIter=1; gen=0;
	cout<<"\nReplicate #"<<rep<<" "<<flush;
	/* check for initial variation in A and B loci */
	varyAInit=1; varyBInit=1;
	b=0;
	for(a=0;a<3;a++) b+=N[0][a];
	if(b==NTot && setAllFix[0]!=2)
	{
		varyAInit=0;
		genFixedAAll=0;
		fixedAAll=0;
	}
	b=0;
	for(a=0;a<3;a++) b+=N[2][a];
	if(b==NTot && setAllFix[0]!=1)
	{
		varyAInit=0;
		genFixedAAll=0;
		fixedAAll=1;
	}
	b=0;
	for(a=0;a<3;a++) b+=N[a][0];
	if(b==NTot && setAllFix[1]!=2)
	{
		varyBInit=0;
		genFixedBAll=0;
		fixedBAll=0;
	}
	b=0;
	for(a=0;a<3;a++) b+=N[a][2];
	if(b==NTot && setAllFix[1]!=1)
	{
		varyBInit=0;
		genFixedBAll=0;
		fixedBAll=1;
	}
	if(varyAInit)
	{
		fixedAAll=2;
		genFixedAAll=genMax+1;
	}
	if(varyBInit)
	{
		fixedBAll=2;
		genFixedBAll=genMax+1;
	}
	fixedACopy=-1;
	fixedBCopy=-1;
	genFixedACopy=genMax+1;
	genFixedBCopy=genMax+1;
	iterExtinct=0;
}

void MUTATE(void)
{
	long a;
	int b;

	for(a=0;a<NTot;a++) for(b=0;b<2;b++)
	{
		/* A mutation */
		if(pNIND[a].AAll[b]==0)
		{
			if(RANDREAL()<mutA[0]) pNIND[a].AAll[b]=1;
		}
		else
		{
			if(RANDREAL()<mutA[1]) pNIND[a].AAll[b]=0;
		}
		/* B mutation */
		if(pNIND[a].BAll[b]==0)
		{
			if(RANDREAL()<mutB[0]) pNIND[a].BAll[b]=1;
		}
		else
		{
			if(RANDREAL()<mutB[1]) pNIND[a].BAll[b]=0;
		}
	}
}

void FREQCALC(void)
{
	long a;
	int b;
	int gtA,gtB;

	for(a=0;a<3;a++) for(b=0;b<3;b++) N[a][b]=0;
	for(a=0;a<NTot;a++)
	{
		gtA=pNIND[a].AAll[0]+pNIND[a].AAll[1];
		gtB=pNIND[a].BAll[0]+pNIND[a].BAll[1];
		N[gtA][gtB]++;
	}
	for(a=0;a<3;a++) for(b=0;b<3;b++) f[a][b]=(double)N[a][b]/NTot;
}

void EXCHANGE(void)
{
	long a;

	for(a=0;a<NTot;a++) *(pIND+a)=*(pNIND+a);
}

void GENBUILD(void)
{
	long a,imm;
	int b,c;
	int parent[2];
	int gTypeParA[2],gTypeParB[2];/* 0=GG, 1=Gg, 2=gg,
										where G = A or B
										and g = a or b */
	int letMate;

	gen++;
	imm=0;
	if(doMig)
	{
		for(a=0;a<NTot;a++)
		{
			rVal=RANDREAL();
			if(RANDREAL()<migP)
			{
				if(rVal<migCumF[0][0]) b=0;			/* AABB */
				else if(rVal<migCumF[0][1]) b=1;	/* AABb */
				else if(rVal<migCumF[0][2]) b=2;	/* AAbb */
				else if(rVal<migCumF[1][0]) b=3;	/* AaBB */
				else if(rVal<migCumF[1][1])
				{
					if(RANDREAL()<0.5) b=4;			/* AB/ab */
					else b=5;						/* Ab/aB */
				}
				else if(rVal<migCumF[1][2]) b=6;	/* Aabb */
				else if(rVal<migCumF[2][0]) b=7;	/* aaBB */
				else if(rVal<migCumF[2][1]) b=8;	/* aaBb */
				else b=9;							/* aabb */
				pNIND[imm].AAll[0]=mIND[b].AAll[0];
				pNIND[imm].AAll[1]=mIND[b].AAll[1];
				pNIND[imm].AID[0]=mIND[b].AID[0];
				pNIND[imm].AID[1]=mIND[b].AID[1];
				pNIND[imm].AHap[0]=mIND[b].AHap[0];
				pNIND[imm].AHap[1]=mIND[b].AHap[1];
				pNIND[imm].BAll[0]=mIND[b].BAll[0];
				pNIND[imm].BAll[1]=mIND[b].BAll[1];
				pNIND[imm].BID[0]=mIND[b].BID[0];
				pNIND[imm].BID[1]=mIND[b].BID[1];
				pNIND[imm].BHap[0]=mIND[b].BHap[0];
				pNIND[imm].BHap[1]=mIND[b].BHap[1];
				imm++;
			}
		}
	}

	if(mateRule==2) for(a=0;a<NTot;a++) pMated[a]=-1; /* enforced monogamy; sets all individuals to unmated */
	for(a=imm;a<NTot;a++)
	{
		letMate=0;
		while(!letMate)
		{
			for(b=0;b<2;b++) /* choosing parents */
			{
				chk=0;
				while(!chk)
				{
					chk=1;
					parent[b]=NTot;
					while(parent[b]==NTot) parent[b]=(int)(RANDREAL()*NTot);
					if(mateRule!=1 && b==1)
					{
						/* If monogamy and parent[0] has mated before, assign prior mate as parent[1] */
						if(mateRule==2)	if(pMated[parent[0]]>=0) parent[1]=pMated[parent[0]];
						/* Self-mating is not allowed; also if monogamy and parent[0] needs a mate */
						if(parent[1]==parent[0] || (mateRule==2 && pMated[parent[0]]==-1)) parent[1]=NTot;
    					{
							while(parent[1]==NTot)
							{
								parent[1]=(int)(RANDREAL()*NTot);
                                if(parent[1]==parent[0]) parent[1]=NTot;
								/* if monogamy and parent[0] has not mated before, be sure mate has also not mated before */
								if(mateRule==2) if(pMated[parent[1]]>=0) parent[1]=NTot;
							}
						}
					}
					/* get parent[b] A genotype */
					gTypeParA[b]=pIND[parent[b]].AAll[0]+pIND[parent[b]].AAll[1];
					/* get parent[b] B genotype */
					gTypeParB[b]=pIND[parent[b]].BAll[0]+pIND[parent[b]].BAll[1];

                    /* if premating isolation and second parent is being selected */
                    if(b==1 && mateCheck==2)
           				if(RANDREAL()>matePref[whichMating[gTypeParA[0]][gTypeParA[1]]])
                            chk=0;
					if(chk && fitCheck) /* selection */
						if(RANDREAL()>w[gTypeParA[b]][gTypeParB[b]]) chk=0;
				}
			}

			/* if postmating isolation and both parents have been selected, make sure mating is allowed */
			letMate=1;
			if(mateCheck==1)
				if(RANDREAL()>matePref[whichMating[gTypeParA[0]][gTypeParA[1]]])
					letMate=0;
			if(letMate && mateRule==2)
			{
				if(pMated[parent[0]]==-1) pMated[parent[0]]=parent[1];
				if(pMated[parent[1]]==-1) pMated[parent[1]]=parent[0];
			}
		}

		/* building offspring genotype */
		for(b=0;b<2;b++)
		{
			if(RANDREAL()<0.5) c=0; else c=1;
			pNIND[a].AAll[b]=pIND[parent[b]].AAll[c];
			pNIND[a].AID[b] =pIND[parent[b]].AID[c];
			pNIND[a].AHap[b]=pIND[parent[b]].AHap[c];
			/* opportunity for recombination if recombination is allowed */
			if(rec>0)
			{
				rVal=RANDREAL();
				if(rVal<rec)
				{
					if(c==0) c=1;
					else c=0;
				}
			}
			pNIND[a].BAll[b]=pIND[parent[b]].BAll[c];
			pNIND[a].BID[b]=pIND[parent[b]].BID[c];
			pNIND[a].BHap[b]=pIND[parent[b]].BHap[c];
		}
	}

	/* opportunity for mutation if mutation is allowed */
	if(mutCheck) MUTATE();
	FREQCALC();
	EXCHANGE();
}

void ALCOMP(void)
{
	if(fixedAAll==2)
	{
		if(N[0][0]+N[0][1]+N[0][2]==NTot && setAllFix[0]!=2)
		{
			fixedAAll=0;
			genFixedAAll=gen;
		}
		else if(N[2][0]+N[2][1]+N[2][2]==NTot && setAllFix[0]!=1)
		{
			fixedAAll=1;
			genFixedAAll=gen;
		}
	}

	if(fixedBAll==2)
	{
		if(N[0][0]+N[1][0]+N[2][0]==NTot && setAllFix[1]!=2)
		{
			fixedBAll=0;
			genFixedBAll=gen;
		}
		else if(N[0][2]+N[1][2]+N[2][2]==NTot && setAllFix[1]!=1)
		{
			fixedBAll=1;
			genFixedBAll=gen;
		}
	}
}

void COALESCE(void)	/* only used if trackCopies=1 */
{
	long a;

	if(fixedACopy==-1)
	{
		fixedACopy=pIND[0].AID[0];
		for(a=0;a<NTot;a++)
		{
			if(pIND[a].AID[0]!=pIND[0].AID[0])
			{
				fixedACopy=-1;
				break;
			}
			if(fixedACopy!=-1) if(pIND[a].AID[1]!=pIND[0].AID[0])
			{
				fixedACopy=-1;
				break;
			}
		}
		if(fixedACopy!=-1)
		{
			genFixedACopy=gen;
			ancAHap=pIND[0].AHap[0];
		}
	}
	if(fixedBCopy==-1)
	{
		fixedBCopy=pIND[0].BID[0];
		for(a=0;a<NTot;a++)
		{
			if(pIND[a].BID[0]!=pIND[0].BID[0])
			{
				fixedBCopy=-1;
				break;
			}
			if(fixedBCopy!=-1) if(pIND[a].BID[1]!=pIND[0].BID[0])
			{
				fixedBCopy=-1;
				break;
			}
		}
		if(fixedBCopy!=-1)
		{
			genFixedBCopy=gen;
			ancBHap=pIND[0].BHap[0];
		}
	}
}

void CONTCHK(void)
{
	int a,b,c;
	int AGtype[3];

	contIter=0;
	if(fixedAAll==2 || fixedBAll==2) contIter=1;
	if(trackCopies)
		if(fixedACopy==-1 || fixedBCopy==-1) contIter=1;
	if(!stopAtFix) contIter=1;
	if(gen==genMax) contIter=0;
	if(!contIter) cout<<" "<<gen<<" generations"<<endl;
	else if(gen%1000==0) cout<<"."<<flush;

	/* checking for fit parents to produce next generation*/
	c=0;
	for(a=0;a<3;a++) for(b=0;b<3;b++) if(w[a][b]>0) c+=N[a][b];
	if((mateRule==1 && !c) || (mateRule!=1 && c<2))
	{
		contIter=0;
		iterExtinct=1;
	}
	/* checking for allowable matings, assuming there are fit parents */
	if(mateCheck && !iterExtinct)
	{
		AGtype[0]=N[0][0]+N[0][1]+N[0][2];
		AGtype[1]=N[1][0]+N[1][1]+N[1][2];
		AGtype[2]=N[2][0]+N[2][1]+N[2][2];
		c=0;
		if(AGtype[0]>0 && matePref[0]>0.0) c=AGtype[0];
		if(mateRule!=1) c--;
		if(!c)
		{
			if(AGtype[1]>0 && matePref[3]>0.0) c=AGtype[1];
			if(mateRule==1) c--;
		}
		if(!c)
		{
			if(AGtype[2]>0 && matePref[5]>0.0) c=AGtype[2];
			if(mateRule==1) c--;
		}
		if(!c) if(AGtype[0] && AGtype[1] && matePref[1]>0.0) c=1;
		if(!c) if(AGtype[0] && AGtype[2] && matePref[2]>0.0) c=1;
		if(!c) if(AGtype[1] && AGtype[2] && matePref[4]>0.0) c=1;
		if(!c)
		{
			contIter=0;
			iterExtinct=1;
		}
	}
}

void ADDNEW(void)
{
	if(!iterExtinct)
	{
		if(fixedAAll<2)
		{
			totAFixed[fixedAAll]++;
			sumGenFixedAAll+=genFixedAAll;
			sumGenFixedAAll2+=(genFixedAAll*genFixedAAll);
		}
		else totAFixed[2]++;
		if(fixedBAll<2)
		{
			totBFixed[fixedBAll]++;
			sumGenFixedBAll+=genFixedBAll;
			sumGenFixedBAll2+=(genFixedBAll*genFixedBAll);
		}
		else totBFixed[2]++;
		if(trackCopies)
		{
			if(genFixedACopy<genMax+1)
			{
				sumGenFixedACopy+=genFixedACopy;
				sumGenFixedACopy2+=(genFixedACopy*genFixedACopy);
			}
			else noFixedACopy++;
			if(genFixedBCopy<genMax+1)
			{
				sumGenFixedBCopy+=genFixedBCopy;
				sumGenFixedBCopy2+=(genFixedBCopy*genFixedBCopy);
			}
			else noFixedBCopy++;
		}
	}
	else
	{
		nIterExtinct++;
		if(trackCopies)
		{
			noFixedACopy++;
			noFixedBCopy++;
		}
	}
}

void OUTWRITE(void)
{
	long a,b,c;
	long cumN;
	double cumF;
	int gtA,gtB;
	long nHap[4];	/* number of each haplotype: AB, Ab, aB and ab */
	double hVal;	/* used to calculate haplotype proportion */

	if(rep==1)
	{
	    /* Write settings to .txt file */
		sumfl<<"\n\nInitial settings\n"<<endl;
		sumfl<<"OUTPUT FILE: "<<sumnm<<endl;
		sumfl<<"RANDOM NUMBER SEED: "<<rsd<<"\n"<<endl;
		sumfl<<"                               Relative"<<endl;
		sumfl<<"Genotype   Number  Proportion  Fitness"<<endl;
		for(a=0;a<3;a++) for(b=0;b<3;b++)
		{
			sumfl<<aGType[a]<<bGType[b]<<"      "<<flush;
			sumfl<<setw(7)<<startPop[a][b]<<"  "<<flush;
            sumfl<<setprecision(8)<<setiosflags(ios::fixed)<<startF[a][b]<<"  "<<flush;
            sumfl<<setprecision(8)<<setiosflags(ios::fixed)<<w[a][b]<<endl;
		}
		sumfl<<"\nTotal Population Size: "<<NTot<<"\n"<<endl;
		sumfl<<"MUTATION RATES:"<<endl;
		sumfl<<"A to a: "<<setprecision(8)<<setiosflags(ios::fixed)<<mutA[0]<<endl;
		sumfl<<"a to A: "<<setprecision(8)<<setiosflags(ios::fixed)<<mutA[1]<<endl;
		sumfl<<"B to b: "<<setprecision(8)<<setiosflags(ios::fixed)<<mutB[0]<<endl;
		sumfl<<"b to B: "<<setprecision(8)<<setiosflags(ios::fixed)<<mutB[1]<<"\n"<<endl;
        if(mateCheck==0) sumfl<<"MATING PREFERENCES (none)"<<endl;
        else if(mateCheck==1) sumfl<<"MATING PREFERENCES (postmating)"<<endl;
        else  sumfl<<"MATING PREFERENCES (premating)"<<endl;
		sumfl<<"AA:AA "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[0]<<"  "<<flush;
		sumfl<<"AA:Aa "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[1]<<"  "<<flush;
		sumfl<<"AA:Aa "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[2]<<endl;
		sumfl<<"Aa:Aa "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[3]<<"  "<<flush;
		sumfl<<"Aa:aa "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[4]<<"  "<<flush;
		sumfl<<"aa:aa "<<setprecision(8)<<setiosflags(ios::fixed)<<matePref[5]<<"\n"<<endl;
        calcVal2=100*rec;
		sumfl<<"LINKAGE BETWEEN A & B GENES: "<<flush;
		sumfl<<setprecision(7)<<setiosflags(ios::fixed)<<calcVal2<<" cM"<<endl;
        sumfl<<"MIGRATION: "<<flush;
        if(!doMig) sumfl<<"none"<<endl;
		else
		{
            sumfl<<"Average number of immigrants/generation: "<<flush;
            sumfl<<setprecision(5)<<setiosflags(ios::fixed)<<migN<<endl;
            sumfl<<"    p(AABB)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][0]<<" "<<flush;
			sumfl<<"  p(AABb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][1]<<" "<<flush;
            sumfl<<"  p(AAbb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[0][2]<<endl;
            sumfl<<"    p(AaBB)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][0]<<" "<<flush;
			sumfl<<"  p(AaBb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][1]<<" "<<flush;
            sumfl<<"  p(Aabb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[1][2]<<endl;
            sumfl<<"    p(aaBB)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][0]<<" "<<flush;
			sumfl<<"  p(aaBb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][1]<<" "<<flush;
            sumfl<<"  p(aabb)="<<flush;
            sumfl<<setprecision(4)<<setiosflags(ios::fixed)<<migF[2][2]<<endl;
		}
        sumfl<<"FERTILIZATION: "<<flush;
		if(!mateRule) sumfl<<" Self-fertilization not allowed\n"<<endl;
		else if(mateRule==1) sumfl<<" Self-fertilization allowed\n"<<endl;
		else sumfl<<" Monogamy enforced\n"<<endl;
		sumfl<<"MAXIMUM NUMBER OF GENERATIONS PER ITERATION: "<<genMax<<endl;
		sumfl<<"ALPHA ALLELE FIXATION TRACKED: "<<flush;
		if(setAllFix[0]==0) sumfl<<"either A or a allele must fix"<<endl;
		else if(setAllFix[0]==1) sumfl<<"A allele must fix"<<endl;
        else sumfl<<"a allele must fix"<<endl;
		sumfl<<"BETA ALLELE FIXATION TRACKED: "<<flush;
		if(setAllFix[1]==0) sumfl<<"either B or b allele must fix"<<endl;
		else if(setAllFix[1]==1) sumfl<<"B allele must fix"<<endl;
        else sumfl<<"b allele must fix"<<endl;
        sumfl<<"ITERATIONS STOPPED AFTER FIXATION: "<<flush;
		if(stopAtFix) sumfl<<"Yes\n"<<endl;
		else sumfl<<"No\n"<<endl;
		/* finish writing settings to .txt file */
		/* write header to .csv file */
		datafl<<"iter,gens,"<<flush;
		datafl<<"A_all_fixed,gens_fix_A_all,finalp_A,gens_fix_A_copy,A_origin_ht,A_origin_gt,"<<flush;
		datafl<<"B_all_fixed,gens_fix_B_all,finalp_B,gens_fix_B_copy,B_origin_ht,B_origin_gt,"<<flush;
        datafl<<"p_AB_AB,p_AB_Ab,p_Ab_Ab,p_AB_aB,p_AB_ab,p_Ab_aB,p_Ab_ab,p_aB_aB,p_aB_ab,p_ab_ab,"<<flush;
        datafl<<"p_AB,p_Ab,p_aB,p_ab\n"<<flush;
        /* finish writing header to .csv file */
	}
	if(!iterExtinct)
	{
		datafl<<rep<<","<<flush;    /* write iteration */
        datafl<<gen<<","<<flush;    /* write number of gens completed */
        /* write allele of gene Alpha fixed, along with generations to fixation */
        if(fixedAAll==0) datafl<<"A,"<<genFixedAAll<<","<<flush;
        if(fixedAAll==1) datafl<<"a,"<<genFixedAAll<<","<<flush;
        if(fixedAAll==2) datafl<<"none,NA,"<<flush;
        /* write final proportion of A allele */
        cumF=0.0;
        for(b=0;b<3;b++) cumF=cumF+f[0][b]+0.5*f[1][b];
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<cumF<<","<<flush;
        /* write data for Alpha gene copy fixation */
        if(trackCopies)
		{
			if(fixedACopy!=-1)
			{
                datafl<<genFixedACopy<<","<<flush;
                datafl<<haplo[ancAHap]<<","<<flush;
				cumN=0;
				chk=1;
				for(a=0;(a<3 && chk);a++) for(b=0;b<3;b++)
				{
					cumN+=startPop[a][b];
					if(fixedACopy<=2*cumN)
					{
						gtA=a; gtB=b; chk=0; break;
					}
					else chk=1;
				}
                datafl<<aGType[gtA]<<bGType[gtB]<<","<<flush;
			}
			else datafl<<"none,NA,NA,"<<flush;
		}
		else datafl<<"NA,NA,NA,"<<flush;
		/* write allele of gene Beta fixe, along with generations to fixation */
        if(fixedBAll==0) datafl<<"B,"<<genFixedBAll<<","<<flush;
        if(fixedBAll==1) datafl<<"b,"<<genFixedBAll<<","<<flush;
        if(fixedBAll==2) datafl<<"none,NA,"<<flush;
        /* write final proportion of B allele */
		cumF=0.0;
		for(a=0;a<3;a++) cumF=cumF+f[a][0]+0.5*f[a][1];
		datafl<<setprecision(7)<<setiosflags(ios::fixed)<<cumF<<","<<flush;
		/* write data for Beta gene copy fixation */
		if(trackCopies)
		{
			if(fixedBCopy!=-1)
			{
                datafl<<genFixedBCopy<<","<<flush;
                datafl<<haplo[ancBHap]<<","<<flush;
				cumN=0;
				chk=1;
				for(a=0;(a<3 && chk);a++) for(b=0;b<3;b++)
				{
					cumN+=startPop[a][b];
					if(fixedBCopy<2*cumN)
					{
						gtA=a; gtB=b; chk=0; break;
					}
					else chk=1;
				}
                datafl<<aGType[gtA]<<bGType[gtB]<<","<<flush;
			}
			else datafl<<"none,NA,NA,"<<flush;
		}
		else datafl<<"NA,NA,NA,"<<flush;
		/* write final genotype proportions */
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[0][0]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[0][1]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[0][2]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[1][0]<<","<<flush;
        c=0;
        for(a=0;a<NTot;a++)
        {
            if(pIND[a].AAll[0]==0 && pIND[a].AAll[1]==1)
                if(pIND[a].BAll[0]==0 && pIND[a].BAll[1]==1) c++;
            if(pIND[a].AAll[0]==1 && pIND[a].AAll[1]==0)
                if(pIND[a].BAll[0]==1 && pIND[a].BAll[1]==0) c++;
        }
        calcVal2=(double)c/NTot;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<calcVal2<<","<<flush;
        calcVal2=(double)(N[1][1]-c)/NTot;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<calcVal2<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[1][2]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[2][0]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[2][1]<<","<<flush;
        datafl<<setprecision(7)<<setiosflags(ios::fixed)<<f[2][2]<<","<<flush;
        /* write final haplotype proportions */
		for(a=0;a<4;a++) nHap[a]=0;
		for(a=0;a<NTot;a++)
		{
			for(b=0;b<2;b++)
			{
				if(pIND[a].AAll[b]==0)
				{
					if(pIND[a].BAll[b]==0) nHap[0]++;
					else nHap[1]++;
				}
				else
				{
					if(pIND[a].BAll[b]==0) nHap[2]++;
					else nHap[3]++;
				}
			}
		}
		for(a=0;a<4;a++)
		{
			hVal=(double)nHap[a]/(2*NTot);
			datafl<<setprecision(7)<<setiosflags(ios::fixed)<<hVal<<flush;
			if(a<3) datafl<<","<<flush;
			else datafl<<endl;
		}
	}
}

void SUMSTATS(void)
{
	double meanAAll,meanBAll,meanACopy,meanBCopy;

	if(trust)
	{
		/* calculate means if more than 1 value */
		if((totAFixed[0]+totAFixed[1])>1)
			meanAAll=(double)sumGenFixedAAll/(totAFixed[0]+totAFixed[1]);
		if((totBFixed[0]+totBFixed[1])>1)
			meanBAll=(double)sumGenFixedBAll/(totBFixed[0]+totBFixed[1]);
		if(trackCopies)
		{
			if(((rep-1)-noFixedACopy)>1)
				meanACopy=(double)sumGenFixedACopy/((rep-1)-noFixedACopy);
			if(((rep-1)-noFixedBCopy)>1)
				meanBCopy=(double)sumGenFixedBCopy/((rep-1)-noFixedBCopy);
		}
		sumfl<<endl;
		sumfl<<"SUMMARY STATISTICS\n"<<endl;
		sumfl<<"FIXATION TIMES"<<endl;
		sumfl<<"                    number   mean     "<<endl;
		sumfl<<"                    -------  ---------"<<endl;
		sumfl<<"Alpha Allele        "<<flush;
		calcVal=totAFixed[0]+totAFixed[1];
		sumfl<<setw(7)<<calcVal<<"  "<<flush;
		if((totAFixed[0]+totAFixed[1])>1)
            sumfl<<setw(9)<<setprecision(1)<<setiosflags(ios::fixed)<<meanAAll<<endl;
        else sumfl<<"nc"<<endl;
		if(trackCopies)
		{
			sumfl<<"Alpha Gene Copy     "<<flush;
            calcVal=(rep-1)-noFixedACopy;
            sumfl<<setw(7)<<calcVal<<"  "<<flush;
			if(((rep-1)-noFixedACopy)>1)
                sumfl<<setw(9)<<setprecision(1)<<setiosflags(ios::fixed)<<meanACopy<<endl;
			else sumfl<<"nc"<<endl;
		}
       	sumfl<<"Beta Allele         "<<flush;
		calcVal=totBFixed[0]+totBFixed[1];
		sumfl<<setw(7)<<calcVal<<"  "<<flush;
		if((totBFixed[0]+totBFixed[1])>1)
            sumfl<<setw(9)<<setprecision(1)<<setiosflags(ios::fixed)<<meanBAll<<endl;
        else sumfl<<"nc"<<endl;
		if(trackCopies)
		{
			sumfl<<"Beta Gene Copy      "<<flush;
            calcVal=(rep-1)-noFixedBCopy;
            sumfl<<setw(7)<<calcVal<<"  "<<flush;
			if(((rep-1)-noFixedBCopy)>1)
                sumfl<<setw(9)<<setprecision(1)<<setiosflags(ios::fixed)<<meanBCopy<<endl;
			else sumfl<<"nc"<<endl;
		}
	}
	sumfl<<"\n"<<endl;
	sumfl<<"ALLELE FIXATION -- OUTPUT REFLECTS SPECIFIED ALLELE FIXATION OPTIONS"<<endl;
	sumfl<<"ALPHA FIXATION"<<endl;
	if(setAllFix[0]!=2) sumfl<<"A allele fixed "<<totAFixed[0]<<" times."<<endl;
	if(setAllFix[0]!=1) sumfl<<"a allele fixed "<<totAFixed[1]<<" times."<<endl;
	sumfl<<"Fixation did not occur "<<totAFixed[2]<<" times.\n"<<endl;
	sumfl<<"BETA FIXATION"<<endl;
	if(setAllFix[1]!=2) sumfl<<"B allele fixed "<<totBFixed[0]<<" times."<<endl;
	if(setAllFix[1]!=1) sumfl<<"b allele fixed "<<totBFixed[1]<<" times."<<endl;
	sumfl<<"Fixation did not occur "<<totBFixed[2]<<" times.\n"<<endl;
	sumfl<<"Iteration halted due to extinction "<<nIterExtinct<<" times."<<endl;
}

void AGAIN(void)
{
	cout<<"\n"<<endl;
	cout<<"Do you want to run the program again? [Y/N] "<<flush;
	cin>>strng;
	fchar=toupper(strng[0]);
	cout<<endl; chk=0;
	if(fchar!='N' && fchar!='Y')
		while(!chk)
		{
			cout<<"\aTry again. [Y/N]? "<<flush;
			cin>>strng;
			fchar=toupper(strng[0]);
			if(fchar=='N' || fchar=='Y') chk=1;
		}
	if(fchar=='Y') doProgAgain=1; else doProgAgain=0;
	if(doProgAgain) firstRun=0;
	cout<<endl;
}

int main()
{
	/* initialize pointers */
	pIND=(struct info *)NULL;
	pNIND=(struct info *)NULL;
	pPERMUTEIND=(struct info *)NULL;
	pMated=(long *)NULL;
	pRand=(double *)NULL;

	doProgAgain=1;
	doSeed=1;
	STARTDEF();		/* sets defaults for menu setting upon startup */
	while(doProgAgain)
	{
		START();  	/* gets menu choices */
		MEMALLOC(); /* allocates memory for population arrays */
		RUNINIT();	/* initializes summary variables for program run */
		for(rep=1;rep<=nIter;rep++)
		{
			/* build and shuffle generation 0 population */
			POPBUILD();
			SHUFFLE();
			ITERINIT();	/* initializes iteration variables */
			while(contIter)
			{
				/* build new generation */
				GENBUILD();
				if(fixedAAll==2 || fixedBAll==2) ALCOMP();
				if(trackCopies) if(fixedACopy==-1 || fixedBCopy==-1) COALESCE();
				CONTCHK();
			}
			OUTWRITE();
			ADDNEW();
		}
		if((rep-1)>1) SUMSTATS();
		datafl.close();
		sumfl.close();
		AGAIN();
	}
	cout<<"\n\nEvolGenius 6.1 done."<<endl;
	return 0;
}