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

#define N 100  // population size
#define L 100 // max memory size
#define G 10000 // max number of different genes in the world

double s=0.1; //selection strength
double v=1.0; //prob of inventing new gene per generation
double p=1.0; //prob of dependency of memes
int nl=100; //number of opportunities to learn
int ngen=10000; //number of generations
int iseed=9095;
double wmax,wbar,nbar;
int gmax;

void one_generation(int meme[N][L],int nmeme[N],double w[N],double h[N],int dep[G],int copy[G]);
void fitness(int meme[N][L],int nmeme[N],double w[N],int dep[G],int copy[G]); 
void print_stuff(int meme[N][L],int nmeme[N],double w[N],double h[N],int dep[G],int copy[G]);

int main(){
	int i,j,igen;
	int meme[N][L],nmeme[N];
	int dep[G],copy[G];
	double w[N],h[N];
	
	srand48(iseed); //seed random number generator
	
	for(i=0;i<N;i++) {
	for(j=0;j<L;j++) {
		meme[i][j]=0; //this is the list of memes in each individual
	}}
	for(i=0;i<G;i++) {
		dep[i]=-1; //if dep[i]>0 meme i depends on meme dep[i]
			// if dep[i]=-1, meme i is independent
		copy[i]=0;
	}
	copy[1]=N; //N copies of meme 1
	
	for(i=0;i<N;i++) {
		meme[i][0]=1; //everyone has meme 1 to start
		nmeme[i]=1; //everyone has 1 meme
		w[i]=1.0+s;
		h[i]=1.0;
	}
	wmax=1.0+s;
	gmax=1;
	
	for(igen=1;igen<=ngen;igen++) {
	one_generation(meme,nmeme,w,h,dep,copy);
	printf("%d %f %f\n",igen,wbar,nbar);
	//print_stuff(meme,nmeme,w,h,dep,copy);
	}
	
	exit(0);
}
//-----------------------------------------------------------------------
void one_generation(int meme[N][L],int nmeme[N],double w[N],double h[N],int dep[G],int copy[G]) 
{
	int i,j,j1,j2,k,ipar,jpar,same,icul,new;
	int meme1[N][L],nmeme1[N];
	double r,w1[N],h1[N],vert;
	

	struct i_parent{
		int group;
		int value;
	};

	struct i_parent a1[1000];

	int flip_value = 0;


	for(i=0;i<N;i++) {
		//choose parent with prob proportional to w[ipar]
		do{
			ipar=drand48()*N;
			flip_value = drand48()*N;

			a1[N].value = ipar;

			if(flip_value > 0.5){
				a1[N].group = 1;
			}else if(flip_value < 0.5){
				a1[N].group = 0;
			}
			r=drand48();
		} while(r>w[ipar]/wmax); 
		do{
			flip_value = drand48()*N;
			
			if(flip_value > 0.5){
				flip_value = 1;
			}else flip_value = 0;
			
			if(a1[N].group == 0){ 
				jpar=drand48()*N;
			}
			if(a1[N].group == 1){
				
				if(flip_value == 1){
					jpar=drand48()*N;
				}
				else {
					while(flip_value!=1){
						flip_value = drand48()*N;
						if(flip_value == 1){
							jpar=drand48()*N;
						}
					}
				}

			}

			r=drand48();
		} while((r>w[jpar]/wmax)||(jpar==ipar)); //two parents must be different
		
		h1[i]=h[ipar];
		//h = prob of horizontal descent of meme (learn from parent)
		//initially this probability is copied without changing from ipar.
		
		//now have M chances of learning a meme
		//with prob vert/2 learn from ipar, prob vert/2 learn from jpar, prob h=1-vert learn from anyone
		j1=0;
		k=0;
		do{
			r=drand48();
			vert=(1.0-h[i])/2.0;
			if(r<vert) { //learn from ipar
				j2=drand48()*nmeme[ipar];
				new=meme[ipar][j2];
			}
			else if(r<vert*2.0) { //learn from jpar
				j2=drand48()*nmeme[jpar];
				new=meme[jpar][j2];
			}
			else { //learn from a random individual
				do{
					icul=drand48()*N;
				} while((drand48()>w[icul]/wmax)||(icul==ipar)||(icul==jpar));
				j2=drand48()*nmeme[icul];
				new=meme[icul][j2];
			}
			//now check whether the new meme is already known
			//if not, learn it
			same=0;
			for(j=0;j<j1;j++) if(meme1[i][j]==new) same=1;
			if(same==0) {
				meme1[i][j1]=new;
				j1++;
			}
			k++;
		} while((k<nl)&&(j1<L));
		nmeme1[i]=j1;
	}
	
	//v is possibility of inventing one new meme per generation
	if(drand48()<v) {
		gmax=1;
		while((copy[gmax]>0)&&(gmax<G)) gmax++;
		if(gmax==G) {
			printf("reached max meme array size\n");
			exit(0);
		}
		i=drand48()*N;
		if(drand48()<p) {
			j=drand48()*nmeme1[i];
			dep[gmax]=meme1[i][j];
		}
		else dep[gmax]=-1;
		
		if(nmeme1[i]==L) {
			meme1[i][L-1]=gmax;//prevents memory overflow
		}
		else {
			meme1[i][nmeme1[i]]=gmax;
			nmeme1[i]++;
		}
	}		
	
	//copy offspring into main array - offspring grow up
	for(i=0;i<N;i++) {
		nmeme[i]=nmeme1[i];
		for(j=0;j<nmeme[i];j++) meme[i][j]=meme1[i][j];
		h[i]=h1[i];
	}
	fitness(meme,nmeme,w,dep,copy);
}
//-------------------------------------------------------------------------
void fitness(int meme[N][L],int nmeme[N],double w[N],int dep[G],int copy[G]) 
{
	int i,j,j1,k,count,g;
	
	wmax=1.0;
	wbar=0.0;
	nbar=0.0;
	
	for(i=0;i<N;i++) {
		if(nmeme[i]==0) w[i]=0.0;
		else {
			w[i]=1.0;
			for(j=0;j<nmeme[i];j++) {
			g=meme[i][j];
			if(dep[g]<0) w[i]*=(1.0+s);
			else {
			for(j1=0;j1<nmeme[i];j1++) if(meme[i][j1]==dep[g]) w[i]*=(1.0+s);
			}
			}
		}
		if(w[i]>wmax) wmax=w[i];
		wbar+=w[i];
		nbar+=nmeme[i];
	}
	wbar/=N;
	nbar/=N;
	
	for(i=1;i<G;i++) copy[i]=0;
	for(i=0;i<N;i++) {
		for(j=0;j<nmeme[i];j++) copy[meme[i][j]]++;
	}
}
//--------------------------------------------------------------------
void print_stuff(int meme[N][L],int nmeme[N],double w[N],double h[N],int dep[G],int copy[G]) 
{
	//prints out lots of details for checking purposes
	int i,j;
	
	for(i=1;i<G;i++) {
	if(copy[i]>0) printf("meme %4d depends on %4d -- copies %d\n",i,dep[i],copy[i]);
	}
	
	for(i=0;i<N;i++) {
	printf("ind %4d fitness %f h %f has memes",i,w[i],h[i]);
	for(j=0;j<nmeme[i];j++) printf(" %4d",meme[i][j]);
	printf("\n");
	}
	
	printf("\n");
}
//--------------------------------------------------------------------