#include "individuals.h"

Parameters para2;

#if DET 
std::mt19937 gen(para2.seed); // how to get para.seed recognized without connecting drift_load.h to individuals.h ? should work
//std::random_device rd2;
//std::mt19937 gen(rd2());
#else
std::random_device rd2;
std::mt19937 gen(rd2());
#endif

std::bernoulli_distribution findhom(0.5); //for sampling the homologue

Individuals::Individuals(bool sx) {
	sex = sx;
	w = 1.0;
	c = 1.0;
	mt = 1.0;
	rep_success = 0;
	moff_success = 0;
	foff_success = 0;
	sum_expWoff = 0.0;
	mat_success = 0;
	dad_ID = 0;
	mum_ID = 0;
	ind_ID = 1;
}


//--------------------------------------------------------------------------
Individuals::~Individuals() {

}


// function to initialise genom with deleterious mutations -> init_genom


void Individuals::init_genom(double k, Parameters para)
{
	int hom;
	double pos;
	double ss;
	double hh;

	std::gamma_distribution<> finds(para.shape_mut,(para.mean_mut/para.shape_mut)); //
	std::uniform_real_distribution<> findpos (0,para.R); // 

	if (para.initial_nMut > 0) {
		

		for (int i = 0; i < para.initial_nMut; i++)  // for each deleterious mutation 
		{
			//sample homologue
			hom = findhom(gen);
			
			//sample selection coefficient
			ss = fmin(finds(gen),1.0);

			
			
			// sample dominance coefficient
			#if DOM // if dominance is modelled,
			
			std::uniform_real_distribution<> findh(0.0, std::exp(-k * ss)); // sample dominance according to emperical results 
																			 
			hh = findh(gen);

			#else 
			hh = 0.5; // codominance
			#endif
			
			//sample position of mutation on chromosome
			pos = findpos(gen);

			//add mutation
			w *= chromo.addDelMutation(hom, pos, ss, hh);
			

		}
		
		
	}
	else
	{
		w = 1;
	}


}


void Individuals::init_cond(Parameters para)
{
	double error;

	if (para.add_error)
	{
		std::normal_distribution<> finder(0.0,para.sd_error);
		error = finder(gen);
	}
	
	else
	{
		error = 0.0;
	}
	
	c = fmax(0.0,w + error); //20/09/20 

	//mt = fmax(0.0, c + error);
	//if (mt > 1.0)mt = 1.0; // adding mating trait  12/03/21
										
}

void Individuals::set_cond(Parameters para)
{
	double error;

	if (para.add_error)
	{
		std::normal_distribution<> finder(0.0, para.sd_error);
		error = finder(gen);
	}

	else
	{
		error = 0.0;
	}

	c = fmax(0.0, w + error);

	//mt = fmax(0.0, c+error);
	//if (mt > 1.0)mt = 1.0;

}


void Individuals::add_neut_mut(Parameters para) {
	
	int hom;
	double neut_effect=0.0;
	std::normal_distribution<> neut_muteffect(0.0, para.neut_Veffect);
	
	neut_effect = neut_muteffect(gen);
	hom = findhom(gen);

	if (hom) chromo.neutM[1] += neut_effect;
	else chromo.neutM[0] += neut_effect;

	

}


void Individuals::deleterious_mutation(int ndel,double k,Parameters para)
{
	int hom;
	double pos;
	double ss;
	double hh;
	

	std::gamma_distribution<> finds(para.shape_mut,(para.mean_mut/para.shape_mut)); //  gamma distribution of fitness effects of del muts
	std::uniform_real_distribution<> findpos(0.0, para.R); // 


	for (int i = 0; i < ndel; i++)  // for each deleterious mutation 
	{
		//sample homologue
		hom = findhom(gen);

		//sample selection coefficient

		if (para.s_const) {

			ss = para.mean_mut; 
		}
		else {
			ss = fmin(finds(gen), 1.0);
		}

		


		// sample dominance coefficient
#if DOM // if dominance is modelled,

		std::uniform_real_distribution<> findh(0.0, std::exp(-k * ss)); // sample dominance according to emperical results 
																		 
																		 
																		 
		hh = findh(gen);

#else 
		hh = 0.5; //  co-dominance
#endif

//sample position of mutation on chromosome
		pos = findpos(gen);

		

		//add mutation
		w *= chromo.addDelMutation(hom, pos, ss, hh);
		
	}



}


void Individuals::back_mutation(int nmut)
{
	int rdn;
	double w_effect=0.0;
	std::map<double, mutation, std::less<double>>::iterator it;

	std::bernoulli_distribution hom(0.5); // from which homolog to delete if homozygote

	for (int i = 0; i < nmut; i++) {

		std::uniform_int_distribution<> samp(0, (int)chromo.mutations.size() - 1);
		rdn = samp(gen);
		it = chromo.mutations.begin();
		std::advance(it, rdn);

		

		if (it->second.homol == 2) { //mutation is homozygote
			
			it->second.homol = hom(gen);
			
			w_effect = (1.0 - it->second.h * it->second.s); // calculate fitness effect of single mutation 
			chromo.Nho--; // delete mutation on one chromosome
			
			//update fitness
			w = w/w_effect;

		}

		else { //mutation is heterozygote
			
			w_effect = (1.0 - it->second.h * it->second.s);
			chromo.mutations.erase(it);
			chromo.nMut--;
			
			//update fitness
		
			w = w/w_effect;
		}

	}

}

void Individuals::benef_mutation(int nmut, Parameters para) {
	int hom;
	double pos;
	double ss;
	std::bernoulli_distribution findhom(0.5); //for sampling the homologue
	std::uniform_real_distribution<> findpos(0, para.R); //
	std::gamma_distribution<> finds(para.shape_mut, (para.mean_mut/para.shape_mut));

	for (int i = 0; i < nmut; i++) {
		//sample homologue
		hom = findhom(gen);
		//sample selection coefficient
		ss = (-1)*fmin(finds(gen),1);
		//sample position
		pos = findpos(gen);
		//add mutation
		w *= chromo.addDelMutation(hom, pos, ss, 1.0); //assume complete dominance of beneficial mutations
	}
}

void Individuals::deleteInd(void) {
	

	
#if DELMUT
	chromo.deleteChromo();
#endif


}

void Individuals::outindsdeme(int repl, int gen, int dnr, std::ofstream *out)
{

	*out << repl << "\t" << gen << "\t" << dnr << "\t" << sex << "\t" << c << "\t" << w << "\t" << mat_success << "\t" << rep_success <<"\t"<< foff_success << "\t"<<moff_success << "\t" << sum_expWoff<< endl;

}