/* allocation simulation
*/

#include	<stdlib.h>
#include 	<stdio.h>
#include 	<time.h>

#define	MAXGENS	 	200
#define	NPLANTS		500
#define	NUMMATES	5
#define	MATEDIST	5
#define	SEASONLEN	100
#define	HERBHITS	11 /* 0 for no herbivory; 99 for fully hit system */
#define	HERBDIST	15
#define	HERBCUE		1.0 /* 0 for none, 1.0 for max */
#define	HERBRATE	0.05
#define	HERBSTART	30
#define	HERBEND		40
#define	HERBENDVAR	0.0 /* 0 for no variability, 1.0 for max */
#define RRR		0.02
#define	MUTATERATE	0.01
#define LZERO		1.0

typedef struct plant_def
{
	double	trait1[SEASONLEN];
	double	trait2[SEASONLEN];

	double	sizet[SEASONLEN];
	double	size;
	double	fitness;
} Plant;
Plant plants1[NPLANTS], plants2[NPLANTS];
int	herbhits[NPLANTS];

void InitSeed(void);
void InitPlants(Plant *plnt);
void ZeroPlants(Plant *plnt);
double HerbHits(void);
int wrapit(int ind);
void GrowPlants(Plant *plnt, double *ave, double *high);
void ReproducePlants(double avefit, double hifit, Plant *parents, Plant *babies);
void Mate(Plant *ma, Plant *pa, Plant *babe);
void Xover(double *mapnt, double *papnt, double *bapnt, double inher);
void PrintPlants(Plant *plnt);

FILE *fopen(const char *, const char *);

int main(void)
{
	int	ttt;
	double	avefit, highfit, avehits;
	Plant	*babes, *adults, *tmpy;
	FILE *fileid;

	adults = plants2;
	babes = plants1;
	InitSeed();
	InitPlants(babes);
	fileid = fopen("time.dat","w");

	for(ttt=0;ttt<MAXGENS;ttt++) /* measurement period */
	{
		tmpy = adults;
		adults = babes;
		babes = tmpy;

		ZeroPlants(adults);
		avehits = HerbHits();
		GrowPlants(adults,&avefit,&highfit);
		ReproducePlants(avefit,highfit,adults,babes);
		if(ttt%5==0)
		{
			printf("%d %f %f %f\n",ttt,avehits,avefit,highfit);
			fprintf(fileid,"%d %f %f %f\n",ttt,avehits,avefit,highfit);
		}
	}
	PrintPlants(adults);
	return(0);
}


void InitSeed(void)
{
	int	seed;
	time_t	nowtime;
	struct	tm *preztime;

	time(&nowtime);
	preztime = localtime(&nowtime);
	seed = (int)((preztime->tm_sec+1)*(preztime->tm_min+1)*
		(preztime->tm_hour+1)*(preztime->tm_year)*(preztime->tm_year));
	if(seed%2==0) seed++;

	srand48(seed);
}

void InitPlants(Plant *plnt)
{
	int i, j;

	for(i=0;i<NPLANTS;i++,plnt++)
	{
		for(j=0;j<SEASONLEN;j++)
		{
			plnt->trait1[j] = drand48();
			plnt->trait2[j] = drand48();
		}
	}
}

void ZeroPlants(Plant *plnt)
{
	int i;

	for(i=0;i<NPLANTS;i++,plnt++)
	{
		plnt->size = LZERO;
		plnt->fitness = 0.0;
		herbhits[i] = 0;
	}
}

double HerbHits(void)
{
	int i,j,k;
	double hits=0.0;

	if(HERBHITS==99)
	{
		for(i=0;i<NPLANTS;i++) herbhits[i] = 1;
		hits = 1.0;
	}
	else
	{
		for(i=0;i<HERBHITS;i++)
		{
			j = NPLANTS*drand48();
			for(k=j-HERBDIST;k<=j+HERBDIST;k++)
				herbhits[wrapit(k)] = 1;
		}
		for(i=0;i<NPLANTS;i++) hits += herbhits[i];
		hits /= (double)NPLANTS;
	}
	return(hits);
}

int wrapit(int ind)
{
	if(ind<0) return(ind+NPLANTS);
	else if(ind>NPLANTS-1) return(ind-NPLANTS);
	else return(ind);
}

void GrowPlants(Plant *plnt, double *ave, double *high)
{
	int i, j, herbbegin, herbend;
	double growth, fitness, *sked;

	*high = fitness = 0.0;
	for(i=0;i<NPLANTS;i++,plnt++)
	{
		herbbegin = herbend = SEASONLEN;
		sked = plnt->trait1;
		if(herbhits[i]==1)
		{
			herbbegin = HERBSTART;
			herbend = herbbegin + (HERBEND-herbbegin+1)
						*(1.0+HERBENDVAR*(2.0*drand48()-1.0));
			if(drand48()<HERBCUE)
					sked = plnt->trait2;
		}

		for(j=0;j<herbbegin;j++)
		{
			growth = RRR*plnt->size;
			plnt->size += (1.0-sked[j])*growth;
			plnt->sizet[j] = plnt->size;
			plnt->fitness += sked[j]*growth;
		}

		sked = plnt->trait2;
		for(j=herbbegin;j<herbend;j++)
		{
			growth = RRR*plnt->size;
			plnt->size += (1.0-sked[j])*growth;
			plnt->size *= (1.0-HERBRATE);
			plnt->sizet[j] = plnt->size;
			plnt->fitness += sked[j]*growth;
		}
		for(j=herbend;j<SEASONLEN;j++)
		{
			growth = RRR*plnt->size;
			plnt->size += (1.0-sked[j])*growth;
			plnt->sizet[j] = plnt->size;
			plnt->fitness += sked[j]*growth;
		}

		fitness += plnt->fitness;
		if(plnt->fitness>*high) *high = plnt->fitness;
	}
	*ave = fitness/(double)NPLANTS;
}

void ReproducePlants(double avefit, double hifit, Plant *parents, Plant *babies)
{
	int	i, j, k, momid[NPLANTS], totmoms, momind;
	Plant	*mom, *dad, *newdad;

/* Reproduction */
	for(i=0;i<NPLANTS;i++) momid[i] = i; /* all plants produce one offspring */
	totmoms = NPLANTS;			/* but they land in random places. */
	for(i=0;i<NPLANTS;i++,babies++)
	{
		k = totmoms*drand48();
		momind = momid[k];
		momid[k] = momid[--totmoms];
		mom = &parents[momind];

		k = momind - MATEDIST + (2*MATEDIST+1)*drand48();
		dad = &parents[wrapit(k)];
		for(j=0;j<NUMMATES-1;j++)
		{
			k = momind - MATEDIST + (2*MATEDIST+1)*drand48();
			newdad = &parents[wrapit(k)];
			if(newdad->fitness>dad->fitness) dad = newdad;
		}
		Mate(mom, dad, babies);
	}
}

void Mate(Plant *ma, Plant *pa, Plant *babe)
{
	double *mapnt, *papnt, *bapnt;

	mapnt = ma->trait1;
	papnt = pa->trait1;
	bapnt = babe->trait1;
	Xover(mapnt,papnt,bapnt,0.5);

	mapnt = ma->trait2;
	papnt = pa->trait2;
	bapnt = babe->trait2;
	Xover(mapnt,papnt,bapnt,0.5);
}

void Xover(double *mapnt, double *papnt, double *bapnt, double inher)
{
	int j;
	double xxx;

	for(j=0;j<SEASONLEN;j++,mapnt++,papnt++,bapnt++)
	{
		xxx = drand48();
		if(xxx<MUTATERATE)
			*bapnt = xxx/MUTATERATE;
		else if((xxx-MUTATERATE)/(1.0-MUTATERATE)<inher)
			*bapnt = *mapnt;
		else	*bapnt = *papnt;
	}
}

void PrintPlants(Plant *plnt)
{
	int i, j;
	double avetrait1[SEASONLEN], avetrait2[SEASONLEN], avesize[SEASONLEN];
	FILE *fileid;

	fileid = fopen("skedz.dat","w");
	for(j=0;j<SEASONLEN;j++) avetrait1[j] = avetrait2[j] = avesize[j] = 0.0;
	for(i=0;i<NPLANTS;i++)
	{
		for(j=0;j<SEASONLEN;j++)
		{
			avetrait1[j] += plnt[i].trait1[j];
			avetrait2[j] += plnt[i].trait2[j];
			avesize[j] += plnt[i].sizet[j];
		}
	}
	for(j=0;j<SEASONLEN;j++)
	{
		printf("%d %f %f %f\n",j,avesize[j]/(double)NPLANTS,
			avetrait1[j]/(double)NPLANTS, avetrait2[j]/(double)NPLANTS);
		fprintf(fileid,"%d %f %f %f\n",j,avesize[j]/(double)NPLANTS,
			avetrait1[j]/(double)NPLANTS, avetrait2[j]/(double)NPLANTS);
	}
}