#include "drift_load.h"


const string Int2Str(const int x)
{
	ostringstream o;
	if (!(o << x))return"ERROR";
	return o.str();
}

// main function 

#if CLUSTER
int main(int argc,char* argv[])
{
	//Get the curretn directroy,
	char*buffer = getcwd(NULL, 0);
	dir = buffer;
	free(buffer);
	dir = dir + "/"; // current directory path
	dirout = dir + "Outputs/";

	para.simNr = std::atoi(argv[1]); // the column for which parameters to change when running in the cluster
	para.reps = std::atoi(argv[2]);
	para.gens = std::atoi(argv[3]);
	
	// track time of the line

	string name = dirout + "Sim" + Int2Str(para.simNr) + "_Para.txt";
	para.outPara(name);

	runmodel();
	std::cout << "Simulation completed" << endl;
	//cin.get();// stop befor finishing and keep console open
	return 0;

}

#else
int main(void)
{
	//Get the curretn directroy
	char *buffer = _getcwd(NULL, 0);
	dir = buffer;
	free(buffer);
	dir = dir + "\\"; // current directory path
	//dirout = dir + "Outputs\\";

	dirout = "C:\\Users\\r01mt18\\Documents\\Models\\purging_driftload\\Outputs\\";
	cout << dirout << endl;
	
	string name = dirout + "Sim" + Int2Str(para.simNr) + "_Para.txt";
	para.outPara(name);
	
	runmodel();

	
	std::cout << "Simulation completed" << endl;
	cin.get();// stop befor finishing and keep console open
	return 0;
}
#endif

// start function : initialises population with N equal to K individuals. 
// needs class population and individuals 



void runmodel(void)
{
	std::cout << "Simulation nr. " << para.simNr << endl;
	extime = clock();
	
	

	if (para.pop_slim) {
		out_pop_header_slim();
	}
	else {
		out_pop_header();
	}

	out_extinction_header();
	out_indsdeme_header();
	if (para.mutout_on) out_dload_header();
	if (para.mutout_on) outPopMut_header();
	if (para.out_Ne) outNe_header();


	for (int rep = 0; rep < para.reps; rep++)
	{
		std::cout << " ============================ REP = " << rep << "  ===============================  " << endl;
		cur_rep = rep;
		critW = 0;
		last_gen = false;
		max_dadID = 0;
		max_mumID =0;
		start();

		for (int g = 0; g < para.gens; g++)
		{
			cur_g = g;

			if (g % 1 == 0) {
				std::cout << " ======================== GEN = " << g << "  ==========================  " << endl;
				extime = clock() - extime;
				//std::cout << "time = " << (float)extime / CLOCKS_PER_SEC << " sec" << endl;
				extime = clock();
			}

			
			survival();
		
			create_deme();
			
			reproduction();
			
			
			if (!demes.empty()) demes.clear();

			if (last_gen)g = para.gens; // finish current replicate when popualtion is extinct


		}
		
		
		// clear population vectors and counters that were filled by reproduction phase
		if (!pop.females.empty()) pop.females.clear();
		if (!pop.males.empty()) pop.males.clear();

		pop.N = 0;
		pop.Nf = 0;
		pop.Nm = 0;

	}

	// close the ouput files
	extinction.close();
	pops.close();
	indsdeme.close();
	if (popmut.is_open()) popmut.close();
	if (dload.is_open()) dload.close();
	if (ne.is_open()) ne.close();

}


// ---> RUNMODEL COMPONENTS <---

void start() {
	

	//pop.N = para.K; // set N equal to K

	for (int i = 0; i <para.K / 2; i++) // initialise females in popualtion
	{
		pop.females.push_back(Individuals(true)); // boolean true (1) is female
		pop.females[i].init_genom(k,para); // give the individuals their properties, genome.. setgenom function
		pop.females[i].init_cond(para); // set condition
		pop.Nf++;
		pop.N++;
	}
	
	
	for (int i = 0; i< para.K / 2; i++)
	{
	pop.males.push_back(Individuals(false)); // do the same for the males.
	pop.males[i].init_genom(k,para);// setgenom 
	pop.males[i].init_cond(para); // set condition
	pop.Nm++;
	pop.N++;
	}

}


// individuals are in males and females vector

//--> LIFE CYCLE FUNCTIONS

void survival() {

	vector<Individuals>::iterator iter;
	std::map<double, mutation>::iterator iter_mut;

	double sum_vf = 0.0; // sum of all the females viability in a population
	double sum_vm = 0.0; // sum of all the males viability in a population
	double sp = 1.0; // probability of survival of an individual

	if (pop.N > 0)
	{

		if (para.dens_ind_sel) // if true,  selection is density independent and followed by density-dependent regulation
		{
			// first selection occurs for both females and males

			for (iter = pop.females.begin(); iter != pop.females.end(); iter++)
			{
				
				sp = fmin(std::pow(iter->w,para.strength_sel),1.0); // selection acts directly on genomic fitness MT 050421
				
			
				bernoulli_distribution selection_dist(sp);

				if (selection_dist(rgen))
				{
					pop.tsuv_females.push_back(*iter);
					pop.Ntf++;
					pop.Nt++;
				}
				else {
					iter->deleteInd(); // delete mutations that ind carries , 					
				}
				

			}


			for (iter = pop.males.begin(); iter != pop.males.end(); iter++)
			{
				
				sp = fmin(std::pow(iter->w, para.strength_sel),1.0); //selection acts directly on genomic fitness MT 050421

				bernoulli_distribution selection_dist(sp);

				if (selection_dist(rgen))
				{
					pop.tsuv_males.push_back(*iter);
					pop.Ntm++;
					pop.Nt++;
				}
				else {
					iter->deleteInd();
					
				}
				
			}


			// caclulate variance in female number of offspring after suvival-selection part

			vector<int> mothers_nr_surv_off; // number of surviving offspring for each mother

			if (cur_g > 0) {

				// fill vectors with 0s

				for (int i = 0; i < max_mumID; i++) {
					mothers_nr_surv_off.push_back(0);
				}

				// females in the previous breeding season (mothers)
				for (int i = 0; i < max_mumID; i++) {

					// male offspring

					for (iter = pop.tsuv_males.begin(); iter != pop.tsuv_males.end(); iter++) {

					
						if (iter->mum_ID == i) {
							
							mothers_nr_surv_off[i]++;
							//break;
						}
						else {
						}

					}

					// female offspring

					for (iter = pop.tsuv_females.begin(); iter != pop.tsuv_females.end(); iter++) {

						
						if (iter->mum_ID == i) {
							
							mothers_nr_surv_off[i]++;
							//break;
						}
						else {
						}
					}
				}

				

				pop.mothers_mean_OFF_postsel = ((double)pop.Ntf + (double)pop.Ntm) / (double)max_mumID; //


				double sumdev_mothers_OFF_postsel = 0.0;


				for (int z = 0; z < max_mumID; z++) {

					sumdev_mothers_OFF_postsel += (((double)mothers_nr_surv_off[z]) - ((double)pop.mothers_mean_OFF_postsel))*(((double)mothers_nr_surv_off[z]) - ((double)pop.mothers_mean_OFF_postsel));
				

				}


				pop.var_mothers_postsel = sumdev_mothers_OFF_postsel / (double)max_mumID;

				


			}


			// now density-dependent regulation occurs

			sp = fmin((double)para.K / (double)pop.Nt, 1.0); // global for females and males
			bernoulli_distribution surv_distr(sp);

			for (iter = pop.tsuv_females.begin(); iter != pop.tsuv_females.end(); iter++)
			{

				if (surv_distr(rgen))
				{
					
					iter->ind_ID = pop.Ntf2;
					pop.tsuv2_females.push_back(*iter);
					pop.sum_c += iter->c; // sum condition
					pop.sum_cf += iter->c; // sum condition for females
					pop.calc_popW(iter->w, iter->sex, iter->chromo.nMut); // enter ind info 
					
					pop.Ntf2++;
					pop.Nt2++;
					pop.n++;
					pop.nf++;

					

					
					
				
					// extract mutations from surviving individual 
					for (iter_mut = iter->chromo.mutations.begin(); iter_mut != iter->chromo.mutations.end(); iter_mut++) {

					
						//add mutation to population mutation map 
						if (para.mutout_on) {

							if (cur_g > para.dload_out_start - 1 && (cur_g % para.dload_interval == 0)) // driftload
								pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h,iter_mut->second.homol);

							else {

								if (critW == 1) { // extract when critW has been reached 

									pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);
								}
							}

						}



					}
					

				}
				else {
					iter->deleteInd();
					
				}

			}

			for (iter = pop.tsuv_males.begin(); iter != pop.tsuv_males.end(); iter++)
			{
				

				if (surv_distr(rgen))
				{
					
					iter->ind_ID = pop.Ntm2;
					pop.tsuv2_males.push_back(*iter);
					pop.sum_c += iter->c; // sum condition 
					pop.sum_cm += iter->c; // sum condition for males
					pop.calc_popW(iter->w, iter->sex, iter->chromo.nMut); // enter ind info 
					pop.Ntm2++;
					pop.Nt2++;
					pop.n++;
					pop.nm++;

					


			
					// extract mutations from surviving individual 
					for (iter_mut = iter->chromo.mutations.begin(); iter_mut != iter->chromo.mutations.end(); iter_mut++) {

						//add mutation to population mutation map 
						if (para.mutout_on) {

							if (cur_g > para.dload_out_start - 1 && (cur_g % para.dload_interval == 0)) // driftload
								pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);

							else {

								if (critW == 1) { // extract when critW has been reached 

									pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);
								}
							}

						}



					}
				}
				else {
					iter->deleteInd();
					
				}

			}

			
			// calculate Ne after offspring mortality has occured

			// get max ID of mum and dad ( last generation Nm-1,Nf-1)

			vector<int> males_nr_surv_moff,males_nr_surv_foff,females_nr_surv_moff,females_nr_surv_foff; // vectors to calculate variance and covariance

			if (cur_g > 0) {

				// fill vectors with 0s

			

				for (int i = 0; i < max_dadID; i++) {
					males_nr_surv_moff.push_back(0);
					males_nr_surv_foff.push_back(0);
				}
				for (int i = 0; i < max_mumID; i++) {
					females_nr_surv_moff.push_back(0);
					females_nr_surv_foff.push_back(0);
				}

				

				// males in the previous breeding season (dads)
				for (int i = 0; i < max_dadID; i++) {
;
					// male offspring

					for (iter = pop.tsuv2_males.begin(); iter != pop.tsuv2_males.end(); iter++) {

					
						if (iter->dad_ID == i) {
							
							males_nr_surv_moff[i]++;
							//break;
						}
						else{
						}

					}


					// female offspring

					for (iter = pop.tsuv2_females.begin(); iter != pop.tsuv2_females.end(); iter++) {
						
						
						if (iter->dad_ID == i) {
							
							males_nr_surv_foff[i] ++;
							//break;
						}
						else {}


					}
				}

				// females in the previous breeding season (mums)

				

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

					
					// male offspring

					for (iter = pop.tsuv2_males.begin(); iter != pop.tsuv2_males.end(); iter++) {
						
						
						if (iter->mum_ID == i) {
						
							females_nr_surv_moff[i] ++;
							//break;
						}
						else {}

					}


					// female offspring

					for (iter = pop.tsuv2_females.begin(); iter != pop.tsuv2_females.end(); iter++) {
						
						
						if (iter->mum_ID == i) {
							
							females_nr_surv_foff[i] ++;
							//break;
						}
						else {}


					}
				}
			
				// vectors with offspring numbers after suvrival phase done
				
				
				// calculate NE demographic including offspring mortality

				if (cur_g> 0) {

					pop.NE2_m_mean_moff = (double)pop.Ntm2 / (double)max_dadID; // mean male offspring 
					pop.NE2_m_mean_foff = (double)pop.Ntf2 / (double)max_dadID; // mean female offspring
					pop.NE2_f_mean_moff = (double)pop.Ntm2 / (double)max_mumID; // mean male offspring 
					pop.NE2_f_mean_foff = (double)pop.Ntf2 / (double)max_mumID; // mean female offspring

					pop.mothers_mean_OFF = ((double)pop.Ntf2 + (double)pop.Ntm2) / (double)max_mumID; // 
					pop.fathers_mean_OFF = ((double)pop.Ntf2 + (double)pop.Ntm2) / (double)max_dadID; //



					double NE2sumdev_f_foff = 0.0;
					double NE2sumdev_f_moff = 0.0;
					double NE2sumdev_f_fmoff = 0.0;
					
					double NE2sumdev_m_foff = 0.0;
					double NE2sumdev_m_moff = 0.0;
					double NE2sumdev_m_fmoff = 0.0;

					double sumdev_mothers_OFF = 0.0;
					double sumdev_fathers_OFF = 0.0;

					for (int z = 0; z < max_mumID; z++) {

						NE2sumdev_f_foff += ((double)females_nr_surv_foff[z] - pop.NE2_f_mean_foff)*((double)females_nr_surv_foff[z] - pop.NE2_f_mean_foff);//

						NE2sumdev_f_moff += ((double)females_nr_surv_moff[z] - pop.NE2_f_mean_moff)*((double)females_nr_surv_moff[z] - pop.NE2_f_mean_moff);

						NE2sumdev_f_fmoff += ((double)females_nr_surv_foff[z] - pop.NE2_f_mean_foff)*((double)females_nr_surv_moff[z] - pop.NE2_f_mean_moff);

						sumdev_mothers_OFF += (((double)females_nr_surv_foff[z] + (double)females_nr_surv_moff[z])-((double)pop.mothers_mean_OFF))*(((double)females_nr_surv_foff[z] + (double)females_nr_surv_moff[z]) - ((double)pop.mothers_mean_OFF));
					}

					for (int z = 0; z < max_dadID; z++) {

						NE2sumdev_m_foff += ((double)males_nr_surv_foff[z] - pop.NE2_m_mean_foff)*((double)males_nr_surv_foff[z] - pop.NE2_m_mean_foff);//

						NE2sumdev_m_moff += ((double)males_nr_surv_moff[z] - pop.NE2_m_mean_moff)*((double)males_nr_surv_moff[z] - pop.NE2_m_mean_moff);

						NE2sumdev_m_fmoff += ((double)males_nr_surv_foff[z] - pop.NE2_m_mean_foff)*((double)males_nr_surv_moff[z] - pop.NE2_m_mean_moff);

						sumdev_fathers_OFF += (((double)males_nr_surv_foff[z] + (double)males_nr_surv_moff[z]) - ((double)pop.fathers_mean_OFF))*(((double)males_nr_surv_foff[z] + (double)males_nr_surv_moff[z]) - ((double)pop.fathers_mean_OFF));
					}

					

					pop.NE2_varF_foff = NE2sumdev_f_foff / (double)max_mumID;
					pop.NE2_varF_moff = NE2sumdev_f_moff / (double)max_mumID;
					pop.NE2_covF_mfoff = NE2sumdev_f_fmoff / (double)max_mumID;

					pop.NE2_varM_foff = NE2sumdev_m_foff / (double)max_dadID;
					pop.NE2_varM_moff = NE2sumdev_m_moff / (double)max_dadID;
					pop.NE2_covM_mfoff = NE2sumdev_m_fmoff / (double)max_dadID;

					pop.var_mothers = sumdev_mothers_OFF / (double)max_mumID;
					pop.var_fathers = sumdev_fathers_OFF / (double)max_dadID;
				
					
				

					pop.NE2_Sk = (((double)max_mumID / ((double)max_dadID + (double)max_mumID))*(pop.NE2_varM_moff + 2.0 * ((double)max_dadID / (double)max_mumID)*pop.NE2_covM_mfoff + (((double)max_dadID / (double)max_mumID)*((double)max_dadID / (double)max_mumID))*pop.NE2_varM_foff)) + (((double)max_dadID / ((double)max_dadID + (double)max_mumID))*(pop.NE2_varF_foff + (2.0*((double)max_mumID / (double)max_dadID)*pop.NE2_covF_mfoff) + (((double)max_mumID / (double)max_dadID)*((double)max_mumID / (double)max_dadID))*pop.NE2_varF_moff));

					

					pop.NE2_d = ((16.0*(double)max_mumID*(double)max_dadID) / (((double)max_dadID + (double)max_mumID))) / (2 + pop.NE2_Sk);
					
					

					pop.NE2_varF_foff =0.0;
					pop.NE2_varF_moff = 0.0;
					pop.NE2_covF_mfoff = 0.0;

					pop.NE2_varM_foff = 0.0;
					pop.NE2_varM_moff = 0.0;
					pop.NE2_covM_mfoff = 0.0;

					
				}
				






			}
			// set counters maximum breeding individuals for next generation



			max_dadID = pop.Ntm2;
			max_mumID = pop.Ntf2; // check if this is the same as in other 
			
			// clear vectors for calculating NE
			if (!males_nr_surv_foff.empty()) males_nr_surv_foff.clear();
			if (!males_nr_surv_moff.empty()) males_nr_surv_moff.clear();
			if (!females_nr_surv_foff.empty()) females_nr_surv_foff.clear();
			if (!females_nr_surv_moff.empty()) females_nr_surv_moff.clear();

			// vector clearing for calculating variation after survival selection 
			if (!mothers_nr_surv_off.empty()) mothers_nr_surv_off.clear();
			

			//surviving individuals are in female and male vectors tsuv_2. Now clear male, female, both tsuv vectors and set counters to zero.

			if (!pop.females.empty()) pop.females.clear();
			if (!pop.tsuv_females.empty()) pop.tsuv_females.clear();
			if (!pop.males.empty()) pop.males.clear();
			if (!pop.tsuv_males.empty()) pop.tsuv_males.clear();

			pop.N = 0;
			pop.Nf = 0;
			pop.Nm = 0;
			pop.Nt = 0;
			pop.Ntf = 0;
			pop.Ntm = 0;





		}

		// if selection is density-dependent
		else
		{
			// first calculate sum of viability values
			double sum_v = 0.0;

			for (iter = pop.females.begin(); iter != pop.females.end(); iter++)
			{
				sum_v += std::pow(iter->w, para.strength_sel); // selection acts directly on genomic fitness MT 050421
				
			}

			for (iter = pop.males.begin(); iter != pop.males.end(); iter++)
			{
				sum_v += std::pow(iter->w, para.strength_sel);  // selection acts directly on genomic fitness MT 050421
				
			}

			
			
			for (iter = pop.females.begin(); iter != pop.females.end(); iter++)
			{
				
				sp = fmin((double)para.K*(std::pow(iter->w, para.strength_sel)/sum_v), 1.0); // selection acts directly on genomic fitness MT 050421

				
				bernoulli_distribution selection_dist(sp);

				if (selection_dist(rgen))
				{
					pop.tsuv2_females.push_back(*iter);
					pop.sum_c += iter->c; // sum condition
					pop.sum_cf += iter->c; // sum condition for females
					pop.calc_popW(iter->w, iter->sex, iter->chromo.nMut); // enter ind info 
					pop.Ntf2++;
					pop.Nt2++;

					// add mutation to mutmap for adult females
					for (iter_mut = iter->chromo.mutations.begin(); iter_mut != iter->chromo.mutations.end(); iter_mut++) {

						
						//add mutation to population mutation map 
						if (para.mutout_on) {

							if (cur_g > para.dload_out_start - 1 && (cur_g % para.dload_interval == 0)) // driftload
								pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);

							else {

								if (critW == 1) { // extract when critW has been reached 

									pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);
								}
							}

						}


					}



				}
				else {
					iter->deleteInd();
					
				}
				

			}


			for (iter = pop.males.begin(); iter != pop.males.end(); iter++)
			{
				sp = fmin(para.K*(std::pow(iter->w, para.strength_sel) / sum_v), 1.0); // selection acts directly on genomic fitness MT 050421

				
				bernoulli_distribution selection_dist(sp);

				if (selection_dist(rgen))
				{
					pop.tsuv2_males.push_back(*iter);
					pop.sum_c += iter->c; // sum condition 
					pop.sum_cm += iter->c; // sum condition for males
					pop.calc_popW(iter->w, iter->sex, iter->chromo.nMut); // enter ind info 
					pop.Ntm2++;
					pop.Nt2++;

					// add mutation to mutmap for adult males
					for (iter_mut = iter->chromo.mutations.begin(); iter_mut != iter->chromo.mutations.end(); iter_mut++) {

						
						//add mutation to population mutation map 
						if (para.mutout_on) {

							if (cur_g > para.dload_out_start - 1 && (cur_g % para.dload_interval == 0)) // driftload
								pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);

							else {

								if (critW == 1) { // extract when critW has been reached 

									pop.addMutation(iter_mut->first, iter_mut->second.s, iter_mut->second.h, iter_mut->second.homol);
								}
							}

						}



					}

				}
				else {
					iter->deleteInd();
					
				}
				
				
			}

			// output into pops
			//if (cur_g %para.out_pop_interval == 0 || (cur_g > para.out2_limit && cur_g%para.out_pop_interval2 == 0)) pop.outPop(cur_rep, cur_g, &pops);// 


			//surviving individuals are in female and male vectors tsuv_2. Now clear male, female, set counters to zero.

			if (!pop.females.empty()) pop.females.clear();
			if (!pop.males.empty()) pop.males.clear();

			pop.N = 0;
			pop.Nf = 0;
			pop.Nm = 0;
			pop.Nt = 0;
			pop.Ntf = 0;
			pop.Ntm = 0;
			
			

		}

	}
	else {
		
		pop.outExtinction(cur_rep, cur_g, &extinction);
		last_gen = true;

	}
}


void create_deme(void) {

	int n_demes = 0;

	
	vector<int> fem_index, male_index;

	for (int i = 0; i < pop.Ntf2; i++) fem_index.push_back(i);
	for (int i = 0; i < pop.Ntm2; i++) male_index.push_back(i);

	random_shuffle(fem_index.begin(), fem_index.end());
	random_shuffle(male_index.begin(), male_index.end());

	int index;
	///

	//std::random_shuffle(pop.tsuv2_females.begin(), pop.tsuv2_females.end());
	//std::random_shuffle(pop.tsuv2_males.begin(), pop.tsuv2_males.end());

	// adjust number of demes to number of males
	// females will then be divided in proportion to males in demes
	// this will result in all but one deme having specified deme size, deme with less males will get proportionally less females
	// if sex ratios are not too unequal, this should result in roughly equal sex ratios within demes
	
	
	n_demes = (int)(pop.Ntm2 / para.deme_size); // deme_size number of males in a deme now
	
	
	
	// test if one deme has to be added because demes and individuals are integers

	double t_ntm,t_ntf, t_ndems, t_demsize; // number of males (as double), number of females (as double), nr of demes, demsize (double, females per male, females per deme

	t_ntm = (double)pop.Ntm2;
	t_demsize = (double)para.deme_size;
	t_ndems = t_ntm / t_demsize;

	

	if (n_demes < t_ndems)
	{
		n_demes = n_demes + 1;
	}

	
	
	// ceate demes according to n_demes

	int count=0;
	int update_deme_size=0;

	for (int i = 0; i < n_demes; i++) {
		demes.push_back(Deme());
		
		// if there are enough individuals in the the deme vector, continue to fill demes accroding to deme_size
		

		if (count < para.deme_size*n_demes - para.deme_size)
		{
				for (int j = 0; j < para.deme_size; j++) {
						
						index = male_index[count];
						demes[i].males.push_back(pop.tsuv2_males[index]);

						
						demes[i].Ndm++;
						demes[i].Nd++;
						count++;
						

						demes[i].sumNmales += pow(exp(demes[i].males[demes[i].Ndm-1].c), para.alpha); // sum the male trait values
						
				}

		}
        
		// if there are less individuals than defined by deme_size, update deme_size size
		else
		{
			
			update_deme_size = pop.Ntm2 - count;
		


			 for (int j = 0; j < update_deme_size; j++) {
				 
				 index = male_index[count];
				 demes[i].males.push_back(pop.tsuv2_males[index]);

			
				 demes[i].Ndm++;
				 demes[i].Nd++;
				 count++;

				 
				 demes[i].sumNmales += pow(exp(demes[i].males[demes[i].Ndm-1].c), para.alpha);
			 }
		}
	}

	count = 0;
	update_deme_size = 0;
	
	// assign females to vectors according to the number of males within a deme, shoudl result in the one deme having less males receiving proportionally less females
		double fpm=0.0, fdm=0.0,db_count2=0.0;
		int count2=0;
		t_ntf = (double)pop.Ntf2; // type double, females in population
		fpm = t_ntf / t_ntm; // type double, females per male
		
		

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

			t_ntf = (double)pop.Ntf2 - (double)count; // update the number of females
			fpm = t_ntf / (t_ntm - db_count2); // type double, updated females per male
			count2 += demes[i].Ndm;
			db_count2 = (double)count2;
		
			
			
			double Ndm2 = (double) demes[i].Ndm; // type double, males per deme[i]
			fdm = fpm * Ndm2; // females tp be assigned to deme[i]
			
			int fdm_integer = (int)fdm;
			double fdm_recast_double = (double)fdm_integer;


			
			if ((fdm_recast_double - fdm) < 0) {
				fdm_integer++;
			}
			
			

			if (count < pop.Ntf2 - fdm_integer)
			{
				
				for (int j = 0; j < fdm_integer; j++) {

					
					index = fem_index[count];
					demes[i].females.push_back(pop.tsuv2_females[index]);

					//demes[i].females.push_back(pop.tsuv2_females[count]);
					demes[i].Ndf++;
					demes[i].Nd++;
					count++;
					
				

				}

			}
				// if there are less individuals than defined by deme_size, update deme_size size
			else
			{
				
				update_deme_size = pop.Ntf2 - count;
				

				for (int j = 0; j < update_deme_size; j++) {
					
					index = fem_index[count];
					demes[i].females.push_back(pop.tsuv2_females[index]);

					//demes[i].females.push_back(pop.tsuv2_females[count]);
					demes[i].Ndf++;
					demes[i].Nd++;
					count++;
					
				}
			}
			
		}

		

	
		fem_index.clear();
		male_index.clear();

	// clear post survival vectors and counters
		if (!pop.tsuv2_females.empty()) pop.tsuv2_females.clear();
		if (!pop.tsuv2_males.empty()) pop.tsuv2_males.clear();

		pop.Nt2 = 0;
		pop.Ntf2 = 0;
		pop.Ntm2 = 0;
}


void reproduction(void) {

	vector<Individuals>::iterator iter_m; // iterator of males in deme i
	vector<Individuals>::iterator iter_f; // iterator of female in deme i

	int counter_m = 0;
	int counter_f = 0;

	// male competition
	double sumNmales; // sum of condition of all males in deme
	double cumprob; // upper boundary for ID male 
	double choice; // choice of ID male
	int IDmale; // ID for female mate in male deme vector

	//Offspring production 
	int noffspring; // number of offspring produced by female
	bool sex; // sex of offspring
	Individuals *ind; //iterator for offspring in inheritance function

	double sumdev_m_foff = 0.0;
	double sumdev_m_moff = 0.0;
	double sumdev_m_fmoff = 0.0;

	double sumdev_f_foff = 0.0;
	double sumdev_f_moff = 0.0;
	double sumdev_f_fmoff = 0.0;

	// rel reproductive succes

	double rel_rs_m = 0.0;
	double sum_m_rs = 0.0;
	double sqr_sum_m_rs = 0.0;
	double rel_rs_f = 0.0;
	double sum_f_rs = 0.0;
	double sqr_sum_f_rs = 0.0;

	// mating success
	double var_ms_m=0.0; //variance in mating success
	double sum_ms_squared=0.0; // sums of squares in mating success
	double sum_ms = 0.0; // square of sums in mating success

	// Mutation
	int ndel = 0; // number of deleterous mutations for each offspring
	int nback = 0;// number of back mutations for each offspring
	int nben = 0; // number of beneficial mutations

	// vectors and distributions
	vector<double> cumdistr; // is the cumulative distribution of males prob of being chosen cumdistr, used to identify ID male
	vector<double> m_foff; // vector of male number of female offspring
	vector<double> m_moff; // vector of male number of male offspring
	vector<double> f_foff; // vector of female number of female offspring
	vector<double> f_moff; // vector of female number of male offspring
	vector<double> m_noff;
	vector<double> f_noff;
	vector<double> m_c, f_c; // adult vectors for condition
	 

	uniform_real_distribution<> mate(0.0, 1.0); // determine the ID of a female mate
	poisson_distribution<> offdistr(para.r); // determine number of offspring when no fecundity selection applies
	bernoulli_distribution sexdist(0.5); // determine sex of offspring

	std::poisson_distribution<> n_dmut(para.Ud); // distribution from which to draw number of deleterious mutations from per newly produced offspring
	std::poisson_distribution<> n_benmut(para.Ub);// distribution to draw beneficial mutations
	bernoulli_distribution neut_mutation(para.neut_mu); // determine if neutral mutation occurs according to mut probability (mu)

	pop.N = 0; // counter of past survival phase is set to 0 again, new offpsring will fill that counter again
	



	// caculate mean condition of population in reproduction phase

	pop.mean_c_f = pop.sum_cf / (double)pop.nf; // mean among females
	pop.mean_c_m = pop.sum_cm / (double)pop.nm; // mean among males
	pop.mean_c = pop.sum_c / (double)pop.n; // whole population

	// calculate mean popualtion fitness and  mean number of mutations per individual in population during reproduction phase


	pop.mean_W = pop.sum_W / (double)pop.n; // mean fitness overall
	pop.mean_fW = pop.sum_fW / (double)pop.nf; // mean fitness females
	pop.mean_mW = pop.sum_mW / (double)pop.nm; // mean fitness males
	pop.mean_nmut = pop.sum_nmut / (double)pop.n; // mean number of mut per ind ??

	int n_demes = (int)demes.size();

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


		if (demes[i].Ndf > 0) {


			cumprob = 0.0; // cumulative probability of being the mate


			for (iter_m = demes[i].males.begin(); iter_m != demes[i].males.end(); iter_m++)
			{
				iter_m->mating_prob = (pow(exp(iter_m->c), para.alpha) / demes[i].sumNmales);
				cumprob += iter_m->mating_prob;
				cumdistr.push_back(cumprob);
				
			}

			demes[i].sumNmales = 0.0;
			int counter = 1;
			// now loop through females and determine mate
			for (iter_f = demes[i].females.begin(); iter_f != demes[i].females.end(); iter_f++)
			{
				
			

				choice = mate(rgen);

			

				IDmale = 0;

				for (int z = 0; z < (int)cumdistr.size(); z++) // iterating over vector that contains the added probabilities
				{

					if (choice <= cumdistr[z]) 
					{
						IDmale = z;
						break;
					}

				}




				// male mate identified in IDmale

				// male mate is identified for each female in deme (at position IDmale), now offspring production


				if (para.fecundity_selection)
				{
					// fecundity-selection function
					double fec = para.r*fmin(pow(iter_f->c, para.fec_par), 1.0); // assuming that para.r is physiological optimum, beneficial mutations can just compensate for del


					poisson_distribution<> offdistr1(fec);
					noffspring = offdistr1(rgen);

				

				}
				else
				{
					if (para.noff_det) {
						noffspring = (int)para.r; // offspring deterministic constant
					}
					else {
						noffspring = offdistr(rgen); // offspring random varialbe drawn from poisson

					}

				}

				// assign (mating and) reproductive sucess to mum and dad in terms of number of offspring 
				iter_f->rep_success += noffspring;
				iter_f->mat_success = 1;
				demes[i].males[IDmale].rep_success += noffspring; // update male reproductive success
				demes[i].males[IDmale].mat_success++; // update male mating success 


				if (noffspring > 0)
				{
					for (int z = 0; z < noffspring; z++)
					{

						sex = sexdist(rgen);


						if (sex)
						{

							// create new individual
							ind = new Individuals(sex);
							inheritance2(ind, *iter_f, demes[i].males[IDmale]); // where fitness is calulated
							iter_f->sum_expWoff += ind->w;// assign mother expected offsrping fitness
							demes[i].males[IDmale].sum_expWoff += ind->w;// assign father expected offsrpring fitness

							iter_f->foff_success++; // update female offpsring count produced for mother


							demes[i].males[IDmale].foff_success++; // update female offspring count produced for father

							
							ind->dad_ID = demes[i].males[IDmale].ind_ID; // give offspring dad id
							ind->mum_ID = iter_f->ind_ID; // give offspring mum id

						

#if BACKMUT
							double mean_backmut= para.Ur * ind->chromo.nMut;
							std::poisson_distribution<> n_backmut(mean_backmut); // distribution from which to draw number of reverse-mutations
							nback = n_backmut(rgen);
							

							if (nback > 0 && ind->chromo.nMut >= nback) ind->back_mutation(nback);

#endif

#if BENMUT
							// beneficial mutations
							nben = n_benmut(rgen);
							if (nben > 0) ind->benef_mutation(nben, para);

#endif

#if DELMUT
							// deleterious mutations
							ndel = n_dmut(rgen);
							
							if (ndel > 0) ind->deleterious_mutation(ndel, k, para);

#endif

							// add mutation on neutral loci // represents mutation rate per individual 
							// is 2*per allele mutation rate

							if (neut_mutation(rgen))
							{
								ind->add_neut_mut(para);

							}

							// add neutral loci value of kid to popualtion in order to calculate Ne from Vneut
							pop.sum_neutral += ind->chromo.neutM[0];
							pop.sum_neutral += ind->chromo.neutM[1];

							

							double squared_neutral = (ind->chromo.neutM[0] * ind->chromo.neutM[0]) + (ind->chromo.neutM[1] * ind->chromo.neutM[1]);

							pop.sum_squared_neutral += squared_neutral;


							// offspring has recombined genom and new mutations occured, now set condtion 

							ind->set_cond(para);



							// Push offspring into population female vector
							pop.females.push_back(*ind);
							pop.Nf++;
							pop.N++;

							ind->deleteInd();
							delete ind;

						}

						else {


							// create new individual
							ind = new Individuals(sex);
							inheritance2(ind, *iter_f, demes[i].males[IDmale]);
							iter_f->sum_expWoff += ind->w; // assign mother expected offspring fitness
							demes[i].males[IDmale].sum_expWoff += ind->w; // assign father expected offspring fitness before mutations occur

							iter_f->moff_success++; // update female offpsring count produced for mother

							demes[i].males[IDmale].moff_success++; // update female offspring count produced for father


							ind->dad_ID = demes[i].males[IDmale].ind_ID; // give offspring dad id
							ind->mum_ID = iter_f->ind_ID; // give offspring mum id

						// add mutations	
						//reverse (back) mutations

#if BACKMUT
							double mean_backmut = para.Ur * ind->chromo.nMut;
							std::poisson_distribution<> n_backmut(mean_backmut); // distribution from which to draw number of reverse-mutations
							nback = n_backmut(rgen);
							
							
						
				

							if (nback > 0 && ind->chromo.nMut >= nback) ind->back_mutation(nback);

#endif

#if BENMUT
							// beneficial mutations
							nben = n_benmut(rgen);
							if (nben > 0) ind->benef_mutation(nben, para);

#endif

#if DELMUT
							// deleterious mutations
							ndel = n_dmut(rgen);
							
							if (ndel > 0) ind->deleterious_mutation(ndel, k, para);

#endif




							// add mutation on neutral loci // represents mutation rate per individual 
							// is 2*per allele mutation rate
							if (neut_mutation(rgen))
							{
								ind->add_neut_mut(para);

							}


							// add neutral loci value of kid to popualtion in order to calculate Ne from Vneut
							pop.sum_neutral += ind->chromo.neutM[0];
							pop.sum_neutral += ind->chromo.neutM[1];

							

							double squared_neutral = (ind->chromo.neutM[0] * ind->chromo.neutM[0]) + (ind->chromo.neutM[1] * ind->chromo.neutM[1]);

							pop.sum_squared_neutral += squared_neutral;

							

							// offspring has recombined genom and new mutations occured, now set condtion 

							ind->set_cond(para);



							// Push offspring into population male vector
							pop.males.push_back(*ind);
							pop.Nm++;
							pop.N++;

							ind->deleteInd();
							delete ind;
						}


					}


					
				}

				else {
					
					iter_f->sum_expWoff = 0.0;
					demes[i].males[IDmale].sum_expWoff += 0.0;
				}

				// 1a.INDIVIDUAL OUTPUT :  output female individual information (adult females in demes, went through natural selection)
				int dnr = i;


				if (cur_g < para.out2_limit && cur_g%para.out_ind_interval == 0) {
					
					iter_f->outindsdeme(cur_rep, cur_g, dnr, &indsdeme);
				}
				else
				{
					if (cur_g > para.out2_limit && cur_g%para.out_ind_interval2 == 0) iter_f->outindsdeme(cur_rep, cur_g, dnr, &indsdeme);

				}

				//pop.calc_popW(iter_f->w, iter_f->sex, iter_f->chromo.nMut); // enter ind info for calculating pop level


			}

			cumdistr.clear();// reset cumulative male values for next deme

		}
		else { 

		}

		// 
		for (iter_m = pop.males.begin(); iter_m != pop.males.end(); iter_m++)
		{
			


		}


		// 1b.INDIVIDUAL OUTPUT :output male individual information,  (adult males in demes, went through natural selection and sexual selection)

		if (cur_g < para.out2_limit && cur_g%para.out_ind_interval == 0) {


			for (iter_m = demes[i].males.begin(); iter_m != demes[i].males.end(); iter_m++)
			{

				
				int dnr = i;
				iter_m->outindsdeme(cur_rep, cur_g, dnr, &indsdeme);


			}
		}
		else
		{
			if (cur_g > para.out2_limit && cur_g%para.out_ind_interval2 == 0)
			{
				int dnr = i;
				iter_m->outindsdeme(cur_rep, cur_g, dnr, &indsdeme);

			}

		}



		// clear deme (adult) vectors and iterators
		for (iter_m = demes[i].males.begin(); iter_m != demes[i].males.end(); iter_m++)
		{
			counter_m++;
			m_foff.push_back(iter_m->foff_success);
			m_moff.push_back(iter_m->moff_success);
			m_noff.push_back(iter_m->rep_success);
			
			sum_ms += (double)iter_m->mat_success;
			sum_ms_squared += (double)iter_m->mat_success*(double)iter_m->mat_success;
			
			// caluclate variance in fitness and condtion //


			pop.sum_sqr_dev_cm += (((double)iter_m->c - pop.mean_c_m)*((double)iter_m->c - pop.mean_c_m));
			pop.sum_sqr_dev_c += (((double)iter_m->c - pop.mean_c)*((double)iter_m->c - pop.mean_c));

			pop.sum_sqr_dev_mW += (((double)iter_m->w - pop.mean_mW) * ((double)iter_m->w - pop.mean_mW));
			pop.sum_sqr_dev_W += (((double)iter_m->w - pop.mean_W) * ((double)iter_m->w - pop.mean_W));

			
			iter_m->deleteInd();

		}

		for (iter_f = demes[i].females.begin(); iter_f != demes[i].females.end(); iter_f++)
		{
			counter_f++;
			f_foff.push_back(iter_f->foff_success);
			f_moff.push_back(iter_f->moff_success);
			f_noff.push_back(iter_f->rep_success);
			//f_c.push_back(iter_f->c); 

			// calculate variance in genomic fitness and condtion //

			pop.sum_sqr_dev_cf += (((double)iter_f->c - pop.mean_c_f)*((double)iter_f->c - pop.mean_c_f));
			pop.sum_sqr_dev_c += (((double)iter_f->c - pop.mean_c)*((double)iter_f->c - pop.mean_c));

			pop.sum_sqr_dev_fW += (((double)iter_f->w - pop.mean_fW) * ((double)iter_f->w - pop.mean_fW));
			pop.sum_sqr_dev_W += (((double)iter_f->w - pop.mean_W) * ((double)iter_f->w - pop.mean_W));

			
			iter_f->deleteInd();
		}

		demes[i].females.clear();
		demes[i].males.clear();

		demes[i].Nd = 0;
		demes[i].Ndf = 0;
		demes[i].Ndm = 0;

	}


	//2. CALCULATE POPULATION LEVEL OUTPUT

	// calculate mean number of offspring per sex 

	pop.mean_noff_f = (double)pop.N / (double)pop.nf; // number of newly pushed offspring divided by number of parents
	pop.mean_noff_m = (double)pop.N / (double)pop.nm;

	var_ms_m= (sum_ms_squared - ((sum_ms*sum_ms)/(double)pop.nm)) / (double)pop.nm;
	pop.Is = var_ms_m / (((double)pop.nf / (double)pop.nm)*((double)pop.nf / (double)pop.nm));



	// MUTATION MAP IN Parent generation 
	// used to be population level mutations of offspring (generation +1) produced, before selection acts on offspring
	

	if (para.mutout_on) {
		if (critW == 1)
		{

			pop.outMutations(pop.n, cur_rep, cur_g, &popmut);
			critW++;
		}

		else {
			if (cur_g > para.full_mutmap_start - 1 && (cur_g % para.full_mutmap_interval == 0))
			{

				pop.outMutations(pop.n, cur_rep, cur_g, &popmut);
			}

		}
	}



		
	// Calculate drift load 

	if (para.mutout_on) {
		if (cur_g > para.dload_out_start - 1 && (cur_g % para.dload_interval == 0)) {

			pop.calc_driftload(pop.n);
			pop.out_fixload(cur_rep, cur_g, pop.dload,pop.mean_W, &dload);

		}
		else {
			if (critW == 2) {
			
				pop.calc_driftload(pop.n);
				pop.out_fixload(cur_rep, cur_g, pop.dload,pop.mean_W, &dload);
				critW++;
			}

		}


	}

	if (pop.mean_W < para.W_critical) critW++;
	
	//3a. NE DEMOGRAPHIC: calcualte demopgraphic ne measure before offspring selection

	pop.f_mean_foff = (double)pop.Nf / (double)pop.nf; // number of female offpsring divided by number of mothers
	pop.f_mean_moff = (double)pop.Nm / (double)pop.nf; // number of male offpsring divided by number of mothers
	pop.m_mean_foff = (double)pop.Nf / (double)pop.nm; // number of female offpsring divided by number of fathers
	pop.m_mean_moff = (double)pop.Nm / (double)pop.nm; // number of male offpsring divided by number of fathers




	if (n_demes > 0) {

		// calcualte variance and covariance in offspring for mothers
		// and relative rep success

		if (counter_f > 0) {
			for (int z = 0; z < counter_f; z++) {

				sumdev_f_foff += ((double)f_foff[z] - pop.f_mean_foff)*((double)f_foff[z] - pop.f_mean_foff);//

				sumdev_f_moff += ((double)f_moff[z] - pop.f_mean_moff)*((double)f_moff[z] - pop.f_mean_moff);

				sumdev_f_fmoff += ((double)f_moff[z] - pop.f_mean_moff)*((double)f_foff[z] - pop.f_mean_foff);

				// calc relative reproductive success
				rel_rs_f = ((double)(f_noff[z] / pop.mean_noff_f)); //  relative reproductive success in f
				sum_f_rs += rel_rs_f;
				sqr_sum_f_rs += rel_rs_f * rel_rs_f;
			}

		}



		pop.varF_foff = sumdev_f_foff / (double)pop.nf;
		pop.varF_moff = sumdev_f_moff / (double)pop.nf;
		pop.covF_mfoff = sumdev_f_fmoff / (double)pop.nf;



		// calculate variance in males offspring
		if (counter_m > 0) {
			for (int z = 0; z < counter_m; z++) {

				sumdev_m_foff += (m_foff[z] - pop.m_mean_foff)*(m_foff[z] - pop.m_mean_foff);
				sumdev_m_moff += (m_moff[z] - pop.m_mean_moff)*(m_moff[z] - pop.m_mean_moff);

				sumdev_m_fmoff += (m_moff[z] - pop.m_mean_moff)*(m_foff[z] - pop.m_mean_foff);

				rel_rs_m = ((double)(m_noff[z] / pop.mean_noff_m)); // 
				sum_m_rs += rel_rs_m;
				sqr_sum_m_rs += rel_rs_m * rel_rs_m;


			}


		}



	}
	pop.var_rel_rs_f = (sqr_sum_f_rs - ((sum_f_rs*sum_f_rs) / pop.nf)) / pop.nf;
	pop.var_rel_rs_m = (sqr_sum_m_rs - ((sum_m_rs*sum_m_rs) / pop.nm)) / pop.nm;



	pop.varM_foff = sumdev_m_foff / (double)pop.nm;
	pop.varM_moff = sumdev_m_moff / (double)pop.nm;
	pop.covM_mfoff = sumdev_m_fmoff / (double)pop.nm;


	// calculate variance in condition and genomic fitness
	pop.var_c = pop.sum_sqr_dev_c / (double)pop.n;
	pop.var_c_f = pop.sum_sqr_dev_cf / (double)pop.nf;
	pop.var_c_m = pop.sum_sqr_dev_cm / (double)pop.nm;

	pop.var_W = pop.sum_sqr_dev_W / (double)pop.n;
	pop.var_fW = pop.sum_sqr_dev_fW / (double)pop.nf;
	pop.var_mW = pop.sum_sqr_dev_mW / (double)pop.nm;

	//  calculation of NE according to Wang et al. 2016 (equ. 2 and 3), see also Nomura (2002) and Hill (1979)

	pop.Sk = (((double)pop.nf / ((double)pop.nm + (double)pop.nf))*(pop.varM_moff + (2.0 * ((double)pop.nm / (double)pop.nf)*pop.covM_mfoff) + (((double)pop.nm / (double)pop.nf)*((double)pop.nm / (double)pop.nf))*pop.varM_foff)) + (((double)pop.nm / ((double)pop.nm + (double)pop.nf))*(pop.varF_foff + (2.0*((double)pop.nf / (double)pop.nm)*pop.covF_mfoff) + (((double)pop.nf / (double)pop.nm)*((double)pop.nf / (double)pop.nm))*pop.varF_moff));


	pop.Ne_d = ((16.0*(double)pop.nf*(double)pop.nm) / (((double)pop.nm + (double)pop.nf))) / (2.0 + pop.Sk);


	//POPULATION OUTPUT: output population level information during reproduction phase

	if (cur_g < para.out2_limit && cur_g%para.out_pop_interval == 0) {

		if (para.pop_slim) {
			
			pop.outPop_slim(cur_rep, cur_g, &pops);//

		}
		else {
			pop.outPop(cur_rep, cur_g, &pops);//
		}
		

	}
	else
	{
		if (cur_g > para.out2_limit && cur_g%para.out_pop_interval2 == 0) pop.outPop(cur_rep, cur_g, &pops);

	}

	// NE GENETIC: output genetic effective population size (measured in offspring generation)

	if (para.out_Ne) {
		double mean_neutral_effect = 0.0;
		double sum_of_square = 0.0;

		mean_neutral_effect = pop.sum_neutral / (2 * pop.N); // calcualte mean 

		
		pop.Vneut = ((pop.sum_squared_neutral - (pop.sum_neutral*pop.sum_neutral) / (2 * pop.N)) / (2 * pop.N));


		// output Ne 

		if (cur_g >= para.out_Ne_start && cur_g%para.out_Ne_interval == 0) {

			pop.outNe(cur_rep, cur_g, pop.Vneut, pop.N, &ne);


		}


	}


	// clear population record and  mean genomic fitness and number of mutations, and caclualtions of Ne

	// Population record:

	counter_m = 0;
	counter_f = 0;
	pop.n = 0;
	pop.nm = 0;
	pop.nf = 0;
	pop.sum_nmut = 0;
	
	// genomic fitness
	pop.sum_W = 0.0;
	pop.sum_mW = 0.0;
	pop.sum_fW=0.0;
	pop.mean_W = 0.0;
	pop.mean_fW = 0.0;
	pop.mean_mW = 0.0;
	
	pop.sum_sqr_dev_W = 0.0;
	pop.sum_sqr_dev_mW = 0.0;
	pop.sum_sqr_dev_fW = 0.0;

	pop.var_W = 0.0;
	pop.var_mW = 0.0;
	pop.var_fW = 0.0;
	
	pop.mean_nmut = 0.0;

	pop.dload = 0.0;

	// NE genetic:
	pop.sum_neutral = 0.0;
	pop.sum_squared_neutral = 0.0;
	pop.Vneut = 0.0;

	// NE demographic:
	pop.varF_foff = 0.0;
	pop.varM_foff = 0.0;
	pop.varF_moff = 0.0;
	pop.varM_moff = 0.0;

	pop.covF_mfoff = 0.0;
	pop.covM_mfoff = 0.0;

	pop.f_mean_foff = 0.0;
	pop.f_mean_moff = 0.0;
	pop.m_mean_foff = 0.0;
	pop.m_mean_moff = 0.0;

	 sumdev_m_foff = 0.0;
	 sumdev_m_moff = 0.0;
	 sumdev_m_fmoff = 0.0;

	 sumdev_f_foff = 0.0;
	 sumdev_f_moff = 0.0;
	 sumdev_f_fmoff = 0.0;

	 pop.mean_noff_f = 0.0;
	 pop.mean_noff_m = 0.0;
	 
	 // NE demogrpahic post selection (gen-1) and var in reprductive success of parents (gen-1)
	 pop.NE2_d = 0.0;
	 pop.var_mothers = 0.0;
	 pop.var_fathers = 0.0;

	 // variance in condition during reproduction phase:
	 pop.sum_c = 0.0;
	 pop.sum_cf = 0.0;
	 pop.sum_cm = 0.0;
	 pop.sum_sqr_dev_cf = 0.0;
	 pop.sum_sqr_dev_cm = 0.0;
	 pop.sum_sqr_dev_c = 0.0;

	pop.mean_c= 0.0;
	pop.mean_c_f = 0.0; 
	pop.mean_c_m = 0.0;

	pop.var_c = 0.0;
	pop.var_c_f = 0.0;
	pop.var_c_m = 0.0;

	// mating success
	var_ms_m = 0.0; 
	sum_ms_squared = 0.0; 
	sum_ms = 0.0; 
	pop.Is = 0.0;

	// relative rep success
	rel_rs_m = 0.0;
	sum_m_rs = 0.0;
	sqr_sum_m_rs = 0.0;
	rel_rs_f = 0.0;
	sum_f_rs = 0.0;
	sqr_sum_f_rs = 0.0;

	pop.var_rel_rs_f = 0.0;
	pop.var_rel_rs_m = 0.0;




	 // empty vectors
	if (!f_foff.empty()) f_foff.clear();
	if (!f_moff.empty()) f_moff.clear();
	if (!m_foff.empty()) m_foff.clear();
	if (!m_moff.empty()) m_moff.clear();
	if (!f_noff.empty()) f_noff.clear();
	if (!m_noff.empty()) m_noff.clear();
	if (!m_c.empty()) m_c.clear();
	if (!f_c.empty()) f_c.clear();


	

	//clear population mutation map of next generation 
	if (!pop.popMuts.empty()) pop.popMuts.clear();

	map<double, pop_muts>::iterator iter_map;

	for (iter_map = pop.popMuts.begin(); iter_map != pop.popMuts.end(); iter_map++) {

		cout << " gen =" << cur_g << " s = " << iter_map->second.s << "h = " << iter_map->second.h << endl;
	}
	
}


void inheritance2(Individuals *kid, Individuals mum, Individuals dad) {
	

	int n_crossovers;
	double cross;
	int hom;
	int hom_neut;

	std::set<double> recomSites;
	std::set<double>::iterator itercross;
	std::map<double, mutation>::iterator iter, iter2;


	std::poisson_distribution<> crossn(para.R);// sample number of crossovers
	std::uniform_real_distribution<> position(0.0, para.R);//Deleterious crossovers positions)
	std::bernoulli_distribution Bern(0.5); // bernoulli trial with prob 0.5
	
	//Recombination (see Roze & Rousset 2009, JEB - Appendix 2)
	//kids homologue 1 comes from the mother----------------------------------------------

	if (mum.chromo.nMut > 0) {

		//sample starting homologue
		hom = Bern(rgen);
		hom_neut = hom;

		//sample no. of crossovers
		n_crossovers = crossn(rgen);

		//sample crossover positions
		for (int i = 0; i < n_crossovers; i++) {
			cross = position(rgen);
			recomSites.insert(cross);
		}

		// NEUTRAL LOCUS
		for (itercross = recomSites.begin(); itercross != recomSites.end(); ++itercross) {
		
			

			if (*itercross < (para.R/(double)2.0)) {
				if (hom_neut == 0) { 
					hom_neut++;
					
				}
				else hom_neut--;
			}
			
		}
		
		kid->chromo.neutM[1] = mum.chromo.neutM[hom_neut]; // kid neut mut inherited from mother is on kid chrom 1

		// MUTATIONS WITH FITNESS EFFECT 
		itercross = recomSites.begin(); //iterator through crossover positions
		iter = mum.chromo.mutations.begin();


		//no mutations before cross-overs positions //
		while (n_crossovers > 0 && *itercross < iter->first) {
			itercross++;
			n_crossovers--;
		}

		for (iter = mum.chromo.mutations.begin(); iter != mum.chromo.mutations.end(); iter++) {
			
																		//cross-overs
			while (n_crossovers > 0 && *itercross < iter->first) {
				if (hom == 0) hom++;
				else hom--;
				itercross++;
				n_crossovers--;
			}
			//if mutation is on the right homologue inherit it, otherwise ignore it
			if (iter->second.homol == hom || iter->second.homol == 2) {
				kid->chromo.mutations[iter->first] = iter->second;
				kid->chromo.mutations[iter->first].homol = 0; //inherit first homologue from mom
				kid->chromo.nMut++;


				//calculate fitness considering the mutation as heterozygote //

				kid->w *= (1.0 - iter->second.h * iter->second.s);


			

			}
		}


		if (!recomSites.empty()) recomSites.clear();

	}
	else {
		// only neutral locus

		//sample starting homologue
		hom = Bern(rgen);
		

		kid->chromo.neutM[1] = mum.chromo.neutM[hom]; // kid neut mut inherited from mother is on kid chrom 1

	}


	if (dad.chromo.nMut > 0) {

		//sample starting homologue
		hom = Bern(rgen);
		hom_neut = hom;

		//sample no. of crossovers
		n_crossovers = crossn(rgen);
		//sample crossover positions
		for (int i = 0; i < n_crossovers; i++) {
			cross = position(rgen);
			recomSites.insert(cross);
		}
		
		// NEUTRAL LOCUS
		for (itercross = recomSites.begin(); itercross != recomSites.end(); ++itercross) {

			

			if (*itercross < (para.R / (double)2.0)) {
				if (hom_neut == 0)
				{
					hom_neut++;
					
				}
				else { 
					hom_neut--;
					
				}
			}

		}

		kid->chromo.neutM[0] = dad.chromo.neutM[hom_neut]; // kid neut mut inherited from dad is on kid chrom 0

		// MUTATIONS WITH FITNESS EFFECT 
		
		
		itercross = recomSites.begin(); //iterator through crossover positions
		iter = dad.chromo.mutations.begin();
		
		
		//no mutations before cross-overs positions
		while (n_crossovers > 0 && *itercross < iter->first) {
			itercross++;
			n_crossovers--;
		}

		for (iter = dad.chromo.mutations.begin(); iter != dad.chromo.mutations.end(); iter++) {
			//crossovers
			while (n_crossovers > 0 && *itercross < iter->first) {
				if (hom == 0) hom++;
				else hom--;
				itercross++;
				n_crossovers--;
			}

			//if mutation is on the right homologue inherit it, otherwise ignore it
			if (iter->second.homol == hom || iter->second.homol == 2) {

				iter2 = kid->chromo.mutations.find(iter->first);
				//if mutation is already present --> it is homozygous
				if (iter2 != kid->chromo.mutations.end()) {
					iter2->second.homol = 2; //mutation is homozygote
					kid->chromo.Nho++;
					//change fitness effect

				

					kid->w /= (1.0 - iter2->second.h * iter2->second.s);

					
					kid->w *= (1.0 - iter2->second.s);


				}
				else { //mutation is heterozygote
					kid->chromo.mutations[iter->first] = iter->second;
					kid->chromo.mutations[iter->first].homol = 1;
					kid->chromo.nMut++;

					//fitness effect as a heterzygot 
					kid->w *= (1.0 - iter->second.h * iter->second.s);

				}
				
				

				
			}

		}


		if (!recomSites.empty()) recomSites.clear();
	}

	else {
		// only neutral locus

		//sample starting homologue
		hom = Bern(rgen);
		
		kid->chromo.neutM[0] = dad.chromo.neutM[hom]; // kid neut mut inherited from dad is on kid chrom 0
	}



}


void expect_offnr(void) {





}

// --> OUTPUT FUNCTIONS 


void out_pop_header(void) {
	string name;

	name = dirout + "Sim" + Int2Str(para.simNr) + "_Pops.txt";
	pops.open(name.c_str());

	pops << "rep\tgen\tnf\tnm\tn\twf\twm\tw\tnmut\tne_d\tne2_d\tif\tim\tis\tvc\tvcf\tvcm\tvGW\tvmGW\tvfGW\tvmRS\tvfRS\tvfRS_s" << endl;


} // population metric during reproduction

void out_pop_header_slim(void) {
	string name;

	name = dirout + "Sim" + Int2Str(para.simNr) + "_Pops.txt";
	pops.open(name.c_str());

	pops << "rep\tgen\tnf\tnm\tw\tnmut\tne2_d\tis\tvcm\tvmGW\tvmRS\tvfRS\tvfRS_s" << endl;


} // population metric during reproduction


void out_indsdeme_header(void) {
	string name;

	name = dirout + "Sim" + Int2Str(para.simNr) + "_Indsdeme.txt";
	indsdeme.open(name.c_str());

	indsdeme << "rep\tgen\tdnr\tsex\tc\tw\tmat_success\trep_success\tfoff\tmoff\tsum_expWoff" << endl; // change this to include moff and foff


}// individual output during reproduction phase in demes

void outPopMut_header(void)
{
	string name;
	name = dirout + "Sim" + Int2Str(para.simNr) + "_PopMut.txt";
	popmut.open(name.c_str());

	popmut << "rep\tgen\ts\th\tfreq" << endl;
} // summarises distribution of mutational effects and h in offsrping generation (gen+1)

void outNe_header(void)
{
	string name;
	name = dirout + "Sim" + Int2Str(para.simNr) + "_Ne.txt";
	ne.open(name.c_str());

	ne << "rep\tgen\tVneut\tNc" << endl;
} // summarises distribution of mutational effects and h in offsrping generation (gen+1)

void out_dload_header(void) {
	string name;
	name = dirout + "Sim" + Int2Str(para.simNr) + "_dload.txt";
	dload.open(name.c_str());

	dload << "rep\tgen\tw\tdload" << endl;


}

void out_extinction_header(void) {
	string name;
	name = dirout + "Sim" + Int2Str(para.simNr) + "_extinction.txt";
	extinction.open(name.c_str());

	dload << "rep\tgen" << endl;

}