/* metacoal.c - written by John Wakeley for Wakeley J and N Aliacar */
/* 2001 Gene genealogies in a metapopulation.  Genetics 159:893-905 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "Obaku:John:Cprograms:Headers:random.h" 
/* the line above includes a file with the random number generator I like; */
/* put your own or a system supplied one in unirandom()/sunirandom() below */

#define ANC_BASE 'A'
#define MUT_BASE 'G'
#define MAXN 11  /* maximum per-deme sample size */

struct node {  /* the nodes in the tree */ 
	int id;
	double time;
	struct node *desc1;
	struct node *desc2;
};

struct tipset { 
	int nprime;  /* size of lineage entering collecting phase */
	int tipids[MAXN];  /* labels of the lineages descendents in the sample */
};

void MyError(error_text)
char  error_text[];
{
	fprintf(stderr,"\nJW's run-time error...\n");
	fprintf(stderr,"\n\t%s\n",error_text);
	fprintf(stderr,"\n...now exiting to system...\n\n");
	exit(1);
}

/* TableSterlings() makes matrices of Stirling numbers of the second kind (s2)  */
/* and the coefficients in the paper (s3) REMEMBER: s3[i][j][k] = S^{(k)}_{j,i} */

double s2[MAXN][MAXN],s3[MAXN][MAXN][MAXN];

void TableSterlings()
{
	int x,y,z; 

	for(x=0;x<MAXN;x++) 
		for(y=0;y<MAXN;y++) {
			s2[x][y] = 0.0;
			for(z=0;z<MAXN;z++)
				s3[x][y][z] = 0.0;
		}
	s2[1][1] = 1.0;
	for(x=2;x<MAXN;x++) {
		for(y=1;y<=x;y++) {
			s2[x][y] = ((double)y)*s2[x-1][y] + s2[x-1][y-1]; 
		}
	}
	for(z=1;z<MAXN;z++) {
		s3[z][z][1] = 1.0;
		for(x=z+1;x<MAXN;x++) 
			for(y=1;y<=(x-z+1);y++) 
				s3[z][x][y] = (double)(x-1)*s3[z][x-1][y] + s3[z][x-1][y-1]; 
	}
}

double pdown(a,r)
double a;
int r;
{
	int x;
	double product;

	product = 1.0;
	for(x=0;x<r;x++)
		product *= (a-((double)x));
	return product;
}

double newprod(ni,j,twoM,twoE)
int ni,j;
double twoM,twoE;
{
	int x;
	double product;

	product = 1.0;
	for(x=j+1;x<=ni;x++)
		product *= ( twoM + ((double)(x-1)) + twoE/((double)x) );
	return product;
}

/* PME() is equation (2) in the paper */

double PME(x,j,ni,twoM,twoE) 
int x,j,ni;
double twoM,twoE;
{
	double denom,efactor;
	
	efactor = (twoE/((double)j))/( twoM + ((double)(j-1)) + twoE/((double)j) );
	denom = newprod(ni,j,twoM,twoE);
	if(j==1)
		return (s3[j][ni][x+1]*pow(twoM,(double)x)/denom);
	else
		return (efactor*s3[j][ni][x+1]*pow(twoM,(double)x)/denom);
}

/* BallBoxProb() is equation (5) in the paper */

double BallBoxProb(n,np,k) 
int n,np,k;
{
	return (s2[n][np]*pdown((double)k,np)/pow((double)k,(double)n));
}

/* ScatterProbME() is equation (6) in the paper */

double ScatterProbME(n,np,k,twoM,twoE) 
int n,np,k;
double twoM,twoE;
{
	int j,x;
	double sum;
	
	if( np > n || np < 1 || n < 1 || n >= MAXN ) MyError("Bad ScatterProbME()");
	
	sum = 0.0;
	for(j=1;j<=n;j++) 
		for(x=0;x<=n-j;x++)
			sum += PME(x,j,n,twoM,twoE)*BallBoxProb(j,np-x,k);
	return sum;
}

/* ES() is Watterson's (1975) expected number of segregating sites */

double ES(n,theta)
int n;
double theta;
{
	int x; 
	double a1;
	
	a1 = 0.0;
	for(x=1;x<n;x++)
		a1 += 1.0/((double)x);
	return theta*a1;
}

/* ESonedeme() is equation (25) in the paper, but for one deme only */

double ESonedeme(n,k,twoM,twoE,theta) 
int n,k;
double twoM,twoE,theta;
{
	int np,x;
	double sum;
	
	if( n < 1 || n >= MAXN ) MyError("Bad n in ESonedeme()");
	
	sum = 0.0;
	for(np=2;np<=n;np++) 
			sum += ScatterProbME(n,np,k,twoM,twoE)*ES(np,theta);
	return sum;
}

/* EStwodemes() is equation (25) in the paper, but for two demes only */

double EStwodemes(n1,n2,k1,k2,twoM1,twoM2,twoE1,twoE2,theta) 
int n1,n2,k1,k2;
double twoM1,twoM2,twoE1,twoE2,theta;
{
	int np1,np2;
	double sum;
	
	if( n1 < 1 || n1 >= MAXN ) MyError("Bad n1 in EStwodemes()");
	if( n2 < 1 || n2 >= MAXN ) MyError("Bad n2 in EStwodemes()");
	
	sum = 0.0;
	for(np1=1;np1<=n1;np1++) 
		for(np2=1;np2<=n2;np2++) 
			sum += ScatterProbME(n1,np1,k1,twoM1,twoE1)*ScatterProbME(n2,np2,k2,twoM2,twoE2)*ES(np1+np2,theta);
	return sum;
}

/*	SUNIRANDOM - Seed your random number generator	*/

void sunirandom(seed)
unsigned int seed;
{
	/* put your own srandom here */
	srandom(seed);
} 

/*	UNIRANDOM - Uniform [0,1) random deviates.	*/

#define IM 2147483648.0

double unirandom()
{
	/* put your own random number generator here */
	return (double)random()/IM; 
}

#undef IM

/*	EXPONENTIALDEV - mean 1, deviates. For mean m, multiple by m.	*/

double exponentialdev() 
{ 
	double dum;
 
	do 
		dum = unirandom();
	while( dum == 0.0 ); 
	return -log(dum); 
} 

/*	LOCALGAMMLN - Natural log of the gamma function.	*/

double localgammln(xx)
double xx;
{
	double x,y,tmp,ser;
	static double cof[6]={76.18009172947146,-86.50532032941677,
		24.01409824083091,-1.231739572450155,
		0.1208650973866179e-2,-0.5395239384953e-5};
	int j;

	y=x=xx;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}

/*	POISSONDEV - Poisson deviates with mean xm.	*/

#define PI 3.141592653589793 
 
double poissondev(xm) 
double xm;
{ 
	double localgammln(); 
	double unirandom(); 
	static double sq,alxm,g,oldm=(-1.0); 
	double em,t,y; 
 
	if (xm < 12.0) { 
		if (xm != oldm) { 
			oldm=xm; 
			g=exp(-xm); 
		} 
		em = -1; 
		t=1.0; 
		do { 
			++em; 
			t *= unirandom(); 
		} while (t > g); 
	} else { 
		if (xm != oldm) { 
			oldm=xm; 
			sq=sqrt(2.0*xm); 
			alxm=log(xm); 
			g=xm*alxm-localgammln(xm+1.0); 
		} 
		do { 
			do { 
				y=tan(PI*unirandom()); 
				em=sq*y+xm; 
			} while (em < 0.0); 
			em=floor(em); 
			t=0.9*(1.0+y*y)*exp(em*alxm-localgammln(em+1.0)-g); 
		} while (unirandom() > t); 
	} 
	return em; 
} 

#undef PI 

/* MakeTipsOneDeme() simulates the scattering phase, returns nprime */
/* for the deme and stores the sizes of these lineages in sizes[]   */

int MakeTipsOneDeme(ni,propagulesize,sizes,twoM,twoE) 
int ni,propagulesize,sizes[MAXN];
double twoM,twoE;
{
	int x,y,n,nprime,numoccupied,pick1,pick2,counter,picks[MAXN],occupancies[MAXN],coalsizes[MAXN];
	double prob,pcoal,pmig,denom;
	
	n = ni;
	if(n<1 || n>=MAXN) 
		MyError("\n... bad n in MakeTipsOneDeme() ...\n\n");
	if(n==1) {
		nprime = 1;
		sizes[0] = 1;
	} else {
		for(x=0;x<n;x++) {
			sizes[x] = 0;
			coalsizes[x] = 1;
		}
		nprime = 0;
		while(n > 0) {
			denom = ((double)(n-1)) + twoM + twoE/((double)n);
			pcoal = ((double)(n-1))/denom;
			pmig = twoM/denom;
			prob = unirandom();
			if(prob < pcoal) {
				pick1 = n*unirandom();
				n -= 1;
				pick2 = n*unirandom();
				if(pick2 >= pick1)
					pick2 += 1;
				coalsizes[pick1] += coalsizes[pick2];
				for(x=pick2;x<n;x++)
					coalsizes[x] = coalsizes[x+1];
			} else 
				if(prob < (pcoal+pmig)) {
					pick1 =  n*unirandom();
					sizes[nprime] = coalsizes[pick1];
					n -= 1;
					for(x=pick1;x<n;x++)
						coalsizes[x] = coalsizes[x+1];
					nprime += 1;
				} else { 
					for(x=0;x<n;x++) {
						picks[x] = propagulesize*unirandom();
						occupancies[x] = 0;
					}
					numoccupied = 0;
					for(x=0;x<n;x++) {
						pick1 = picks[x];
						if(pick1>=0) {
							for(y=0;y<n;y++)
								if(picks[y]==pick1) {
									occupancies[numoccupied] += 1;
									picks[y] = -1;
								}
							numoccupied += 1;
						}
					}
					counter = 0;
					for(x=0;x<numoccupied;x++) 
						for(y=0;y<occupancies[x];y++) {
							sizes[nprime+x] += coalsizes[counter];
							counter += 1;
						}
					nprime += numoccupied; 
					n = 0;
				} 
		}
	}
	return nprime;
}

/* MakeTips() does the scattering phase for all sampled demes */

int MakeTips(numdemes,demesams,propagulesizes,twoMs,twoEs,scattertips)
int numdemes,*demesams,*propagulesizes;
double *twoMs,*twoEs;
struct tipset *scattertips;
{
	int i,x,y,nprime,niprime,tipcounter1,tipcounter2,sizes[MAXN];
	
	nprime = 0;
	tipcounter1 = 0;
	tipcounter2 = 0;
	for(i=0;i<numdemes;i++) {
		niprime = MakeTipsOneDeme(demesams[i],propagulesizes[i],sizes,twoMs[i],twoEs[i]);
		for(x=0;x<niprime;x++) {
			scattertips[tipcounter1].nprime = sizes[x];
			for(y=0;y<sizes[x];y++) {
				scattertips[tipcounter1].tipids[y] = tipcounter2;
				tipcounter2 += 1;
			}
			tipcounter1 += 1;
		}
		nprime += niprime;
	}
	return nprime; 
}

void InitTree(numseqs,tree)
int numseqs;
struct node *tree;
{
	int x;

	for(x=0;x<(2*numseqs);x++) {
		(tree + x)->id = x;
		(tree + x)->time = 0.0;
		(tree + x)->desc1 = NULL;
		(tree + x)->desc2 = NULL;
	}
}

/* MakeTree() standard coalescent tree building routine */

void MakeTree(numseqs,theta,tree)
int numseqs;
double theta;
struct node *tree;
{
	int x,pick,cnode;
	double t;
	struct node **list;
	
	list = (struct node **)malloc( (size_t)numseqs*sizeof(struct node *) );
	if(list == NULL) MyError("problem with list memory allocation");
	for(x=0;x<numseqs;x++)
		list[x] = tree + x;
	t = 0.0;
	for(x=numseqs;x>1;x--) {
		cnode = 2*numseqs - x;
		pick = x*unirandom();
		tree[cnode].desc1 = list[pick];
		list[pick] = list[x-1];
		pick = (x-1)*unirandom();
		tree[cnode].desc2 = list[pick];
		list[pick] = tree + cnode; 
		t += exponentialdev()*theta/(double)(x*(x-1));
		tree[cnode].time = t;
	}
	free( (void *)list );
}

/* SizeChange() lineages before "time" become "factor" times longer */

void SizeChange(samplesize,time,factor,tree)
int samplesize;
double time,factor;
struct node *tree;
{
	int x;
	
	for(x=samplesize;x<(2*samplesize-1);x++) {
		if( tree[x].time > time )
			tree[x].time = factor*(tree[x].time - time) + time;
	}
}

double TreeLength(numseqs,tree)
int numseqs;
struct node *tree;
{
	int x;
	double Ttot;
	
	Ttot = 0.0;
	for(x=(2*numseqs-2);x>=numseqs;x--) {
		Ttot += (tree[x].time - (tree+x)->desc1->time); 
		Ttot += (tree[x].time - (tree+x)->desc2->time);   
	}
	return Ttot;
}

int PickBranch(numseqs,Ttot,tree)
int numseqs;
double Ttot;
struct node *tree;
{
	int x;
	double Tsum,Tpick;

	Tsum = 0.0;
	Tpick = Ttot*unirandom();
	for(x=(2*numseqs-2);x>=numseqs;x--) {
		Tsum += (tree[x].time - (tree+x)->desc1->time); 
		if( Tpick < Tsum ) 
			return ( (tree+x)->desc1->id ); 
		Tsum += (tree[x].time - (tree+x)->desc2->time);   
		if( Tpick < Tsum ) 
			return ( (tree+x)->desc2->id ); 
	}
}

void MutateTips(numscat,site,datamatrix,thisnode,scattertips)
int numscat,site;
char **datamatrix;
struct node *thisnode;
struct tipset *scattertips;
{
	int x;
	
	if( thisnode->id >= numscat ) {
		MutateTips(numscat,site,datamatrix,thisnode->desc1,scattertips);
		MutateTips(numscat,site,datamatrix,thisnode->desc2,scattertips);
	} else 
		for(x=0;x<scattertips[thisnode->id].nprime;x++)
			datamatrix[scattertips[thisnode->id].tipids[x]][site] = MUT_BASE;
}

void MakeDataTips(numscat,numseqs,S,datamatrix,Ttot,tree,scattertips)
int numscat,numseqs,S;
char **datamatrix;
double Ttot;
struct node *tree;
struct tipset *scattertips;
{
	int x,y,nodenum;

	for(x=0;x<numseqs;x++) {
		for(y=0;y<S;y++) 
			datamatrix[x][y] = ANC_BASE;
		datamatrix[x][S] = '\0';
	}
	for(x=0;x<S;x++) { 
		nodenum = PickBranch(numscat,Ttot,tree); 
		MutateTips(numscat,x,datamatrix,tree+nodenum,scattertips);
	}
}

/* CountFreqs() site-frequency distribution in the entire sample */

void CountFreqs(numseqs,numpats,S,fullfreqs,foldfreqs,datamatrix)
int numseqs,numpats,S;
double *fullfreqs,*foldfreqs;
char **datamatrix;
{
	int x,y,sitefreq;
	double count;

	for(x=1;x<numseqs;x++)
		fullfreqs[x] = 0.0;
	for(x=1;x<=numpats;x++)
		foldfreqs[x] = 0.0;
	for(x=0;x<S;x++) {
		sitefreq = 0;
		for(y=0;y<numseqs;y++)
			if(datamatrix[y][x] == MUT_BASE)
				sitefreq += 1;
		fullfreqs[sitefreq] += 1.0;
		for(y=1;y<=numpats;y++) 
			if( sitefreq == y || sitefreq == (numseqs-y) ) 
				foldfreqs[y] += 1.0;
	}
}

void PrintData(numseqs,S,datamatrix,fptr)
int numseqs,S;
char **datamatrix;
FILE *fptr;
{
	int x;

	fprintf(fptr," %d %d \n",numseqs,S);
	for(x=0;x<numseqs;x++) {
		fputs(datamatrix[x],fptr);
		fputc('\n',fptr);
	}
	fputc('\n',fptr);
	fflush(fptr);
}

void WriteTime(start,finish,fptr)
long start,finish;
FILE *fptr;
{
	double seconds,hours,minutes;
	
	seconds = (double)(finish-start);
	modf(seconds/3600.0,&hours);
	seconds -= 3600.0*hours;
	modf(seconds/60.0,&minutes);
	seconds -= 60.0*minutes;
	fprintf(fptr,"\n\nTotal Running Time = %1.0f:%1.0f:%1.0f\n\n",hours,minutes,seconds);
	fflush(fptr);
}

int main() 
{
	int analytic,simulation,seed;
	int x,y,z,S,numseqs,numpats,numsites,numscat,numdemes,numreps,*demesams,*propagulesizes;
	double totave,totexp,theta,Q,T,Ttot,*twoMs,*twoEs,*fullfreqs,*foldfreqs,*fullaves,*foldaves,*nprimes;
	char **datamatrix;
	struct node *tree;
	struct tipset *scattertips;
	time_t start,finish;
	int i,j,imax=7,jmax=7;  /* use maxs to vary 2M and 2E; must be <= size of the following two arrays */ 
	double M2s[7]={0.1,0.31623,1.0,3.1623,10.0,31.623,100.0}; /* index i below */
	double E2s[7]={0.1,0.31623,1.0,3.1623,10.0,31.623,100.0}; /* index j below */
	FILE *fptr;
			
	TableSterlings();
	
	/* fptr = fopen("meta_data","w"); */
	fptr = stdout;
	numdemes = 2;
	theta = 10.0;  
	Q = 1.0;  /* instantaneous change in Ne, Q=NeA/Ne (wasn't in paper) */
	T = 1.0;  /* when the instantaneous change in Ne happened, in units of Ne generations */
	numreps = 10000;
	analytic=1;
	simulation=1;
	
	/* This loop prints out a matrix of expected numbers of segregating sites in a single */
	/* deme, over the ranges of 2M and 2E specified above; format is for importing into   */
	/* Mathematica. Right now it can print out the expected number in one or two demes.  */
	if(analytic==1) {
		fputc('{',fptr);
		for(i=0;i<imax;i++) {
			fputc('{',fptr);
			for(j=0;j<jmax;j++) {
				demesams = (int *)malloc( (size_t)numdemes*sizeof(int) );
				if(demesams == NULL) MyError("problem getting memory for demesams in main()");
				propagulesizes = (int *)malloc( (size_t)numdemes*sizeof(int) );
				if(propagulesizes == NULL) MyError("problem getting memory for propagulesizes in main()");
				twoMs = (double *)malloc( (size_t)numdemes*sizeof(double) );
				if(twoMs == NULL) MyError("problem getting memory for twoMs in main()");
				twoEs = (double *)malloc( (size_t)numdemes*sizeof(double) );
				if(twoEs == NULL) MyError("problem getting memory for twoEs in main()");
				numseqs = 0;
				for(x=0;x<numdemes;x++) {
					twoMs[x] = M2s[i];
					twoEs[x] = E2s[j];
					propagulesizes[x] = 1;
					demesams[x] = 5;
					numseqs += demesams[x];
				}
				if(numdemes==1)
					fprintf(fptr,"%lf",ESonedeme(demesams[0],propagulesizes[0],twoMs[0],twoEs[0],theta)); 
				if(numdemes==2)
					fprintf(fptr,"%lf",EStwodemes(demesams[0],demesams[1],propagulesizes[0],propagulesizes[1],twoMs[0],twoMs[1],twoEs[0],twoEs[1],theta)); 
				if(numdemes!=1 && numdemes!=2)
					MyError("Analytical loop can do more than two demes: set analytical=1 in main");
				free( (void *)demesams );
				free( (void *)propagulesizes );
				free( (void *)twoMs );
				free( (void *)twoEs );
				if(j<(jmax-1)) 
					fputc(',',fptr);
				fflush(fptr);	
			}
			fputc('}',fptr);
			if(i<(imax-1)) {
				fputc(',',fptr);
				fputc('\n',fptr);
			}
		}
		fputc('}',fptr);
		fputc('\n',fptr);
		fputc('\n',fptr);
	} 

	/* This loop simulates data.  Right now it is set up to compute the average site-frequency */
	/* distribution for a sample from one deme, over the ranges of 2M and 2E specified above,  */
	/* but the only thing it prints out now is the average total number of segregating sites.  */
	/* If numdemes = 1 or 2, and same propagule sizes, this and the analytic loop should agree */
	if(simulation==1) {
		for(i=0;i<imax;i++) {
		fputc('{',fptr);
			for(j=0;j<jmax;j++) {
				demesams = (int *)malloc( (size_t)numdemes*sizeof(int) );
				if(demesams == NULL) MyError("problem getting memory for demesams in main()");
				propagulesizes = (int *)malloc( (size_t)numdemes*sizeof(int) );
				if(propagulesizes == NULL) MyError("problem getting memory for propagulesizes in main()");
				twoMs = (double *)malloc( (size_t)numdemes*sizeof(double) );
				if(twoMs == NULL) MyError("problem getting memory for twoMs in main()");
				twoEs = (double *)malloc( (size_t)numdemes*sizeof(double) );
				if(twoEs == NULL) MyError("problem getting memory for twoEs in main()");
				numseqs = 0;
				for(x=0;x<numdemes;x++) {
					twoMs[x] = M2s[i];
					twoEs[x] = E2s[j];
					propagulesizes[x] = 1;
					demesams[x] = 5;
					numseqs += demesams[x];
				}
				numpats = (int)floor(((double)numseqs)/2.0);
				nprimes = (double *)malloc( (size_t)(numseqs+1)*sizeof(double) );
				if(nprimes == NULL) MyError("problem getting memory for nprimes in main()"); 
				fullfreqs = (double *)malloc( (size_t)numseqs*sizeof(double) );
				if(fullfreqs == NULL) MyError("problem getting memory for fullfreqs in main()");
				fullaves = (double *)malloc( (size_t)numseqs*sizeof(double) );
				if(fullfreqs == NULL) MyError("problem getting memory for fullaves in main()");
				foldfreqs = (double *)malloc( (size_t)(numpats+1)*sizeof(double) );
				if(foldfreqs == NULL) MyError("problem getting memory for foldfreqs in main()");
				foldaves = (double *)malloc( (size_t)(numpats+1)*sizeof(double) );
				if(foldaves == NULL) MyError("problem getting memory for foldaves in main()");
				for(y=1;y<=numseqs;y++) 
					nprimes[y] = 0.0; 
				for(y=1;y<numseqs;y++) 
					fullaves[y] = 0.0;
				for(y=1;y<=numpats;y++) 
					foldaves[y] = 0.0;
				time(&start); 
				sunirandom( (unsigned int)start ); 
				for(x=0;x<numreps;x++) {
					scattertips = (struct tipset *)malloc( (size_t)numseqs*sizeof(struct tipset) );
					if(scattertips == NULL) MyError("problem getting memory for scattertips in main()");
					numscat = MakeTips(numdemes,demesams,propagulesizes,twoMs,twoEs,scattertips);
					nprimes[numscat] += 1.0; 
					tree = (struct node *)malloc( (size_t)2*numscat*sizeof(struct node) );
					if(tree == NULL) MyError("problem getting memory for tree in main()");
					InitTree(numscat,tree);
					MakeTree(numscat,theta,tree);
					SizeChange(numscat,T*theta,Q,tree);
					Ttot = TreeLength(numscat,tree);
					S = (int)poissondev(Ttot);
					if(S>0) {  
						datamatrix = (char **)malloc( (size_t)numseqs*sizeof(char *) );
						if( datamatrix == NULL ) MyError("prob w/datamatrix in main()");
						for(y=0;y<numseqs;y++) {
							datamatrix[y] = (char *)malloc( (size_t)(S+1)*sizeof(char) );
							if( datamatrix[y] == NULL ) MyError("prob w/datamatrix[y] in main()");
						}
						MakeDataTips(numscat,numseqs,S,datamatrix,Ttot,tree,scattertips);
						/* PrintData(numseqs,S,datamatrix,fptr); */
						CountFreqs(numseqs,numpats,S,fullfreqs,foldfreqs,datamatrix); 
						for(y=1;y<numseqs;y++) 
							fullaves[y] += fullfreqs[y];
						for(y=1;y<=numpats;y++) 
							foldaves[y] += foldfreqs[y];
						for(y=0;y<numseqs;y++) 
							free( (void *)datamatrix[y] );
						free( (void *)datamatrix ); 
					} 
					free( (void *)tree ); 
					free( (void *)scattertips );
				} 
				/* time(&finish);
				WriteTime(start,finish,stdout); */
				totave = 0.0;
				totexp = 0.0;
				for(y=1;y<numseqs;y++) {
					totave += fullaves[y];
					totexp += theta/((double)y);
				}
				fprintf(fptr,"%lf",totave/((double)numreps));
				if(j<(jmax-1)) 
					fputc(',',fptr);
				fflush(fptr);
				free( (void *)nprimes ); 
				free( (void *)fullfreqs );
				free( (void *)fullaves );
				free( (void *)foldfreqs );
				free( (void *)foldaves );
				free( (void *)demesams );
				free( (void *)propagulesizes );
				free( (void *)twoMs );
				free( (void *)twoEs );
			}
			fputc('}',fptr);
			if(i<(imax-1)) {
				fputc(',',fptr);
				fputc('\n',fptr);
			}
		}
		fputc('}',fptr);
	  	fputc('\n',fptr);
	  	fputc('\n',fptr);
	} 
	
	fclose(fptr);
	return 0;
}