/*****
 * Carlo C. Maley   July 19, 1994
 * 
 * These are the methods for the Organism, Population and Context Objects in a
 * model inspired by R.C. Lewontin.  (See Lew.cpp for description)
 *
 *****/

#include <stdio.h>
#include <stdlib.h>
#include <oops.h>
#include <math.h>
#include "LewObjects.h"
#include "windowUtils.h"
#include "random.h"
#include "poisson.h"

extern	WindowPtr	DrawingWindow;	// The drawing window is
									// defined in windowUtils.c

/*--------------------------------------------------------------------------*/

#define FIFTYFIFTY uniform(2)  /*FIFTYFIFTY gives a 50/50 chance of returning a 1.*/
#define GETPERCENT uniform(PERCENT) /*Gets a random number 0-99.*/

#define MUTRATE 5 /* The mutation gene mutates at this rate*/
double gene_probs[RANDPERCENT];

/* Useful procedures ------------------------------------------*/

Boolean check_gene(unsigned char gene)
{		
	return (Boolean) (knuth_random() < gene_probs[gene]);
}

void setup_probabilities()
{
	int i;
	
/*	for(i=1; i < PERCENT; i++)
		gene_probs[i] = 1 / pow((double) i, MUTSCALER);
	gene_probs[0] = 1.0;
	gene_probs[PERCENT] = 0.0;
*/

	/*This creates a set of poisson lambdas for the tag mutation from 0 to 16.*/
	for(i=0; i <= PERCENT; i++)
		gene_probs[i] = pow((double) i / POISSONSCALE, POISSONPOWER);
				
}

/* Methods for cOrganism --------------------------------------*/

cOrganism::cOrganism() /*initializes an organism at time of creation.*/
{
	longevity = (unsigned char) uniform(RANDPERCENT);
	mutation = (unsigned char) uniform(RANDPERCENT);
	pursue = (unsigned char) uniform(MAXCHAR);
	escape = (unsigned char) uniform(MAXCHAR);
	
	alive = TRUE;
	fitness = 0;
	age = 0;
	num_offspring = 0;
}

Boolean cOrganism::aging()  /*check for death by senescence.*/
{
	age++;
	
	if (uniform(PERCENT) >= longevity)/*(knuth_random() < gene_probs[longevity]) */
	{
		die();
		return TRUE;
	}
	else return FALSE;
}

/*
void cOrganism::die()		 
{
	alive = FALSE;
}
*/

void cOrganism::give_birth(cOrganism *progeny)	 /*Produce an offpring.*/
{
	progeny->be_born(repro_gene(longevity, mutation),
					 repro_gene(mutation, MUTRATE),  /***This used to be mutation.*/
					 repro_tag(pursue, mutation),
					 repro_tag(escape, mutation));
	num_offspring++;

}

  /*Initialize new born from parental genes.*/
void cOrganism::be_born(unsigned char longe, unsigned char mutat, unsigned char purs, unsigned char esc)  
{
	longevity = longe;
	mutation = mutat;
	pursue = purs;
	escape = esc;
	/* alive = TRUE; This happens in ruach() to prevent newborns from reproducing.*/
	age = 0;
	num_offspring = 0;
	fitness = 0;
}


  /*Score a double bonus for a match of my pursuit tag against my opponents escape
  tag.  But then subtract my opponent's pursuit tag match against my escape tag. 
  This gives a prisoner's dilemma when we consider a mutation of the escape tag to
  be a defection.  Other "pursreward" and "escreward" give different payoffs.*/

void cOrganism::interact(cOrganism *other, double pursreward, double escreward) 
{
	int mypursuit, myescape;
	unsigned char compare;
	int i;
	
	/*First score my pursuit vs. my opponent's escape...*/
	/*XOR pursue and other->escape*/
	compare = pursue ^ other->escape;
	
	/*Count the number of 0's.*/
	for(i=0, mypursuit = 0; i < TAGLENGTH; i++)
		if (((compare >> i) % 2) == 0) mypursuit++;
	
	/*Subtract half the tag length and that is the score.*/
	mypursuit -= TAGLENGTH / 2;
	
	fitness += pursreward * mypursuit;  /*A successful pursuit counts for double.*/


	/*Second score my opponent's pursuit vs. my escape...*/
	/*XOR other->pursue and escape*/
	compare = other->pursue ^ escape;
	
	/*Count the number of 0's.*/
	for(i=0, myescape = 0; i < TAGLENGTH; i++)
		if (((compare >> i) % 2) == 0) myescape++;
	
	/*Subtract half the tag length and that is the score.*/
	myescape -= TAGLENGTH / 2;
	
	fitness -= escreward * myescape;
	
	
}

/*returns the "alive" instance variable.
Boolean cOrganism::living() 
{
	return alive;
}
*/

void cOrganism::print_state(int id)
{
	printf("%d :: longe: %d mutat: %d purs: %x esc: %x fit: %d age: %d kids: %d\n", 
			id, longevity, mutation, pursue, escape, fitness, age, num_offspring);
}

void cOrganism::save_state(FILE *endfp)
{
	fprintf(endfp, "%3d\t %3d\t %.2x\t %.2x\t %3d\t %3d\t %3d\n", longevity, mutation,
			pursue, escape, fitness, age, num_offspring);
}

/*
void cOrganism::new_turn()
{
	fitness = 0;
}
*/

/*
int cOrganism::get_fitness()
{
	return fitness;
}
*/

	/* Protected methods for Organisms -------------------------*/
	
  /* At reproduction gene_probs[mutat] number of bits are flipped on average
     (with a poisson distribution around that average).  This basically means
     0 to 8 bits.*/
	
unsigned char cOrganism::repro_tag(unsigned char tag, unsigned char mutat)
{
	int num_mut;
	int loc;	
	
	if (gene_probs[mutat] > 0) num_mut = poisson(gene_probs[mutat]);
	else num_mut = 0;
	
	/* printf("Old tag: %x ", tag);*/

	for(; num_mut > 0; num_mut--)
	{
		loc = uniform(TAGLENGTH);
		tag = (1 << loc) ^ tag;
	}
	
	
/*	
	int i;
	unsigned char heritability;
	
	heritability = PERCENT - mutat;
	
	
	for(i=0; i < TAGLENGTH; i++)
	{
		if (uniform(PERCENT) < mutat) 
			if (FIFTYFIFTY) tag = (1 << i) ^ tag;
	}
	
*/
	
	return tag;
}


  /* The new gene probably ought to vary from the parent's gene in a poisson 
  distribution with a width determined by the mutation gene.  The boundaries
  of 0 and 100 cut off this distribution.*/

unsigned char cOrganism::repro_gene(unsigned char gene, unsigned char mutat)
{
	int newgene;
	
	if (mutat <= 0) newgene = gene;
	else newgene = poisson((double) mutat) - mutat + gene;
	
	/*newgene = ((randgene * (PERCENT - mutat)) + (gene * mutat)) / PERCENT;*/
	
	if (newgene > 100) newgene = 100;
	if (newgene < 0) newgene = 0;
	
	/*printf("Old gene: %d  New gene: %d\n", gene, newgene);*/
			
	return newgene;

}


/* Methods for cPopulation -------------------------------------------------*/

void cPopulation::update_pop(void)
{
	/*printf("updating the population.\n");*/
	
	/*death_by_senescence();*/
	death_by_competition();
	
	/*printf("The remaining population: \n");
	print_pop();*/
	
	reproduction();
}

void cPopulation::init_pop(cContext *context) /*Let everyone live for one time step.*/
{
	int i;
	
	next_corpse = 0;
	numorgs = context->pop_size;
	
	corpses = (unsigned int *) calloc(numorgs, sizeof(int));
	for(i=0; i < numorgs; i++)
		corpses[i] = 0;
		
	orgs = (cOrganism **) calloc(numorgs, sizeof(cOrganism*));
	for(i = 0; i < numorgs; i++)
		orgs[i] = new cOrganism;
		
	ave_fitness = 0.0;
	radius = context->radius;
	
	next_parent = 0;
	next_competitor = 0;
	max_corpses = (int) (numorgs * GENGAP);
}


cPopulation::~cPopulation()
{
	int i;
	
	free(corpses);
	
	for(i=0; i < numorgs; i++)
		delete orgs[i];
	free(orgs);
}

void cPopulation::print_pop()
{
	int i;
	
	for(i=0; i < numorgs; i++)
		if (orgs[i]->living()) orgs[i]->print_state(i);
}


 /*Creates a new random population (once init_pop has been done).*/

void cPopulation::reset_pop()
{
	int i;
	
	for(i=0; i < numorgs; i++)
		orgs[i]->be_born((unsigned char) uniform(RANDPERCENT), (unsigned char) uniform(RANDPERCENT),
						 (unsigned char) uniform(MAXCHAR), (unsigned char) uniform(MAXCHAR));

}

void cPopulation::save_pop(FILE *endfp)
{
	int i;
	
	fprintf(endfp, "Average Fitness: %f\n", ave_fitness);
	fprintf(endfp, "Lon\t Mut\t Pr\t Es\t Fit\t Age\t Kids\n");
	
	for(i=0; i < numorgs; i++)
		orgs[i]->save_state(endfp);

}


	/* Protected Population methods ---------------------------------------*/
		
void cPopulation::competition(cContext *context)
{
	int i, j;
	int rad;
	
	ave_fitness = 0.0;
	
	for(i=0; i < numorgs; i++)
	{
		orgs[i]->new_turn();
		orgs[i]->aging(); /*because it used to happen in senescence.*/
		for(rad = radius; rad > 0; rad--)
		{
			/*printf("%d interacts with %d\n", i, (i + rad) % numorgs);*/
			orgs[i]->interact(orgs[(i + rad) % numorgs], 
							  (double) context->pursuit_factor,
							  (double) context->escape_factor);
		}
		for(rad = -1 * radius; rad < 0; rad++)
		{
			/*printf("%d interacts with %d\n", i, (i + rad + numorgs) % numorgs);*/
			orgs[i]->interact(orgs[(i + rad + numorgs) % numorgs], 
							  (double) context->pursuit_factor,
							  (double) context->escape_factor);
		}
		
		ave_fitness += orgs[i]->get_fitness();
	}
	
	ave_fitness = ave_fitness / numorgs;
	/*printf("Average fitness this time step = %f\n", ave_fitness);*/
}



/*What about organisms that fall right on the average?  Give them a 50/50 chance.*/
/*  The old algorithm.  If there are a few bad competitors, only they will die:
	for(i=0; i < numorgs; i++)
	{
		if (orgs[i]->living())
		{
			if (orgs[i]->get_fitness() < ave_fitness ||
				(orgs[i]->get_fitness() == ave_fitness && FIFTYFIFTY))
			{
					orgs[i]->die();
					bury(i);
			}
		}
	}
*/

/* The carrying capacity of the enivoronment is popsize * (1 - GENGAP).  
   Every time step, popsize * GENGAP must die due to resource limitations.*/
   
void cPopulation::death_by_competition()
{
	int i;
	double cutoff;
	
	for(cutoff = ave_fitness; next_corpse < max_corpses; cutoff++)
	{
		for(i=0; 
			i < numorgs && next_corpse < max_corpses; 
			i++, next_competitor = (next_competitor + 1) % numorgs)
		{
			if (orgs[next_competitor]->living())
			{
				if (orgs[next_competitor]->get_fitness() < cutoff)
				{
					orgs[next_competitor]->die();
					bury(next_competitor);
				}
			}
		}
	}
	
}

void cPopulation::death_by_senescence()
{
	int i;
	
	for(i=0; i < numorgs; i++)
		if (orgs[i]->living())
			if (orgs[i]->aging()) bury(i);
}

  /*If there is no one left alive, generate a new random population.  Otherwise,
    keep cycling around*/
void cPopulation::reproduction()
{
	int i, offspring, num_dead;
	
	if (next_corpse >= numorgs) reset_pop();
	else
	{
		num_dead = next_corpse;
		for(i=next_parent, offspring = get_next_corpse(); 
			offspring >= 0; 
			offspring = get_next_corpse(), i = (i+1) % numorgs)
		{
			/*find the next living organism.*/
			while(!orgs[i]->living())
				i = (i+1) % numorgs;
				
			/*let it reproduce into the next corpse position in the stack.*/
			orgs[i]->give_birth(orgs[offspring]);	
		}
		next_parent = (i + 1) % numorgs;
		for(i=0; i < num_dead; i++)
			orgs[corpses[i]]->ruach();  /*The breath of life - sets alive = TRUE.*/
	}

}


	/*Returns -1 if there are no more corpses, otherwise it pops the next corpse off the 
	  top of the corpses stack.*/

int cPopulation::get_next_corpse()
{
	if (next_corpse <= 0) return -1;
	else
	{
		next_corpse--;
		return corpses[next_corpse];
	}
}

	/*Pushes a new corpse onto the corpses stack.*/								

void cPopulation::bury(int corpse_id) 
{
	if (next_corpse >= numorgs) 
		printf("Error in bury: there are more dead than in the total population.\n");
	else
	{
		corpses[next_corpse] = corpse_id;
		next_corpse++;
	}
}

/* Methods for cContext --------------------------------------------------*/

/*The constructor: asks the user for pop_size and max_iterations.*/
cContext::cContext()	
{
	printf("\nHow many trials would you like to run? ");
	scanf("%d", &num_trials);
	fflush(stdin);

	printf("\nWhat is the starting trial number? ");
	scanf("%d", &current_trial);
	fflush(stdin);

	printf("\nWhat population size would you like to work with? ");
	scanf("%d", &pop_size);
	fflush(stdin);
	
	printf("What will be the radius of interaction? ");
	scanf("%d", &radius);
	fflush(stdin);
	
	printf("Reward multiplier for pursuit? ");
	scanf("%f", &pursuit_factor);
	fflush(stdin);
	
	printf("Reward multiplier for escape? ");
	scanf("%f", &escape_factor);
	fflush(stdin);
	
	printf("How many time steps do you want to run? ");
	scanf("%d", &max_iterations);
	fflush(stdin);
	
	printf("Random seed (negative = clock)? ");
	scanf("%d", &randseed);
	fflush(stdin);
	randseed = seed_random(randseed);
	/*printf("The random number seed = %d\n", randseed);*/
	
	printf("Prefix for the files to save the final states? ");
	scanf("%s", prefix);
	fflush(stdin);
	

	
	time_step = 0;  /*Set the clock to zero.*/
	endfp = NULL;
	
	setup_probabilities();
}


/*Returns TRUE if time_step < max_iterations, while updating time_step.*/

Boolean cContext::running()  
{
	if (time_step < max_iterations) return TRUE;
	else return FALSE;
}