/* allocation simulation
cc simpall.c -o bsimp -lm
*/

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

#define	MAXGENS	 	200
#define	NPLANTS		500
#define	SEASONLEN	100
#define RRR		0.02
#define	MUTATERATE	0.01
#define LZERO		1.0

typedef struct plant_def
{
	double	trait[SEASONLEN];
	double	size;
	double	fitness;
} Plant;
Plant plants1[NPLANTS], plants2[NPLANTS];

void InitSeed(void);
void InitPlants(Plant *plnt);
void ZeroPlants(Plant *plnt);
void GrowPlants(Plant *plnt, double *ave, double *high);
void ReproducePlants(double avefit, Plant *parents, Plant *babies);
void MakeBrakes(double avefit, double *totr, double *breaks, Plant *parents);
Plant *GetParent(double totr, double *breaks, Plant *parents);
void Mate(Plant *ma, Plant *pa, Plant *babe);
void PrintPlants(Plant *plnt);

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

	adults = plants2;
	babes = plants1;
	InitSeed();
	InitPlants(babes);

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

		ZeroPlants(adults);
		GrowPlants(adults,&avefit,&highfit);
		ReproducePlants(avefit,adults,babes);
		if(ttt%5==0)
			printf("%d %f %f\n",ttt,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->trait[j] = drand48();
}

void ZeroPlants(Plant *plnt)
{
	int i;

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

void GrowPlants(Plant *plnt, double *ave, double *high)
{
	int i, j;
	double growth, fitness=0.0;

	*high=0.0;
	for(i=0;i<NPLANTS;i++,plnt++)
	{
		for(j=0;j<SEASONLEN;j++)
		{
			growth = RRR*plnt->size;
			plnt->size += (1.0-plnt->trait[j])*growth;
			plnt->fitness += plnt->trait[j]*growth;
		}
		fitness += plnt->fitness;
		if(plnt->fitness>*high) *high = plnt->fitness;
	}
	*ave = fitness/(double)NPLANTS;
}

void ReproducePlants(double avefit, Plant *parents, Plant *babies)
{
	int	i;
	double	totrep;
	double	brakes[NPLANTS];
	Plant	*mom, *dad;

/* measure the total reproduction in the habitat */
	MakeBrakes(avefit, &totrep, brakes, parents);

/* Reproduction */
	for(i=0;i<NPLANTS;i++,babies++)
	{
		mom = GetParent(totrep, brakes, parents);
		dad = GetParent(totrep, brakes, parents);
		Mate(mom, dad, babies);
	}
}

void MakeBrakes(double avefit, double *totr, double *breaks, Plant *parents)
{
	int i;

	*totr = 0.0;
	for(i=0;i<NPLANTS;i++,parents++)
	{
		if(parents->fitness>avefit)
			*totr += parents->fitness;
		breaks[i] = *totr;
	}
}

Plant *GetParent(double totr, double *breaks, Plant *parents)
{
	int j, rnum;

	rnum =  drand48() * totr;
	j = (int)((rnum*NPLANTS)/totr);

	if(rnum<breaks[0]) j = 0;
	else if(breaks[j]<rnum)
		while(rnum > breaks[j]) j++;
	else
	{
		while(rnum < breaks[j]) j--;
		j++;
	}
	return(&parents[j]);
}

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

	mapnt = ma->trait;
	papnt = pa->trait;
	bapnt = babe->trait;
	for(j=0;j<SEASONLEN;j++,mapnt++,papnt++,bapnt++)
	{
		xxx = drand48();
		if(xxx<MUTATERATE)
			*bapnt = xxx/MUTATERATE;
		else if(xxx<0.5+MUTATERATE/2.0)
			*bapnt = *mapnt;
		else	*bapnt = *papnt;
	}
}

void PrintPlants(Plant *plnt)
{
	int i, j, biggie=0;
	double hifit=0.0, avetrait[NPLANTS];

	for(j=0;j<SEASONLEN;j++) avetrait[j] = 0.0;
	for(i=0;i<NPLANTS;i++)
	{
		if(plnt[i].fitness>hifit)
		{
			biggie = i;
			hifit = plnt[i].fitness;
		}
		for(j=0;j<SEASONLEN;j++)
			avetrait[j] += plnt[i].trait[j];
	}
	for(j=0;j<SEASONLEN;j++)
		printf("%d %f %f\n",j,
			avetrait[j]/(double)NPLANTS,
			plnt[biggie].trait[j]);
}