population.h

#ifndef DEF_POPULATION_H
#define DEF_POPULATION_H
#include <iterator>
#include <algorithm>
#include "humandna.h"
#include "demographic_function.h"
#include "boost/utility/enable_if.hpp"
 #include "boost/type_traits/is_base_of.hpp"

/* Population is a class that contains  2 array of human: female and male
// humans can be of type human, humanmt, humandna -> we use template on those type
this population evolve through non overlapping generation with cycles of Mutation - Reproduction*/
template <class T>
class Population // <> indicates a template specification defintion in the .cpp file
{
 private:
  std::ostream& print(std::ostream&)const;
  void random_shuffle_pop(); // handwritten shuffling of indiv of current generation: fast handy and efficient
  void initiate_pop(int N);
  void initiate_pop(int N,int mt,int xL, int xS,int y,int aL,int aS);
  void reinitiate_ID_table(); 
  std::vector<int> femaleIDtable; 
  std::vector<int> maleIDtable;  
  int nbHuman; 
  int nbFemale; 
  int generation; 
  int nbCrumbMt; 
  int nbCrumbA; 
  int nbCrumbY; 
  int nbCrumbX; 
  int nbCrumbNuc; 
  int nbCellPerA; 
  int nbCellPerX; 
  int nbChromosomesA; 
  int nbChromosomesX; 
  int sizeMt; 
  int sizeY; 
  int sizePerX; 
  int sizePerA; 
  double popsizeDemo; 
  T tmpmale; // handy and fast
  T tmpfemale; // handy and fast
  double muMt1;
  double muMt2;
  double muX1;
  double muX2;
  double muY1;
  double muY2;
  double muA1; // Transition A<->G C<->T i.e. change on 1 bit
  double muA2; // Transversion A<->C A<->T G<->C G<->T i.e. change on 2 bits
  std::vector<double> growthFactor; // exponential growth factor
  int matingSystem; 
  // 1 : completely random, monogamy 
  // 2 : polygamy, still poisson law with mean nb of child ~2 for the nb of child per women
  // 3 : polygyny, males can mate several female (without restriction), but female mate only one male
  // 4 : polyandry -> high variance in the nb of children / women
  // 5 : true random matig -> no more poisson in the nb of children
  int inbreeding;   // 1 : no inbreeding between brothers and sisters

 public:
  // *** Public attributes ***
  std::vector<T> female[2]; 
  std::vector<T> male[2]; 
  // *** Constructor / destructor ***
  Population() // dull
    { initiate_pop(0);}
  explicit Population(int N,int mt,int xL, int xS,int y,int aL,int aS) 
  { 
    initiate_pop(N,mt,xL,xS,y,aL,aS);
  } // assume mu=0
  Population(Population const& a) 
    { 
      initiate_pop(a.nbHuman,a.sizeMt,a.nbChromosomesX,a.sizePerX,a.sizeY,a.nbChromosomesA/2,a.sizePerA);
      generation = a.generation;
      nbFemale = a.nbFemale;
      popsizeDemo = a.popsizeDemo;
      matingSystem = a.matingSystem;
      inbreeding = a.inbreeding;
      muA1 = a.muA1;
      muA2 = a.muA2;
      muX1 = a.muX1;
      muX2 = a.muX2;
      muY1 = a.muY1;
      muY2 = a.muY2;
      muMt1 = a.muMt1;
      muMt2 = a.muMt2;
      female[0] = a.female[0];
      male[0] = a.male[0];
      female[1] = a.female[1];
      male[1] = a.male[1];
      femaleIDtable = a.femaleIDtable;
      maleIDtable = a.maleIDtable;
      growthFactor = a.growthFactor;
    }
  explicit Population(std::string const& filename);  // <>
  ~Population(){
  }; 
  // *** std get/set ***
  std::vector<T> get_pop(int current) const; 
  Population& operator = (Population const&);
  friend std::ostream& operator << (std::ostream& os, Population const&a)
  {
    a.print(os);
    return os;
  }
  void reset(int N);
  void set_mating_system(int); 
  int get_matingSystem() const; 
  void set_inbreeding(int); 
  int get_inbreeding() const; 
  void set_mu(double,double,double,double); 
  // mt x y a
  void set_mu(double,double,double,double,double,double,double,double); 
  void set_growthFactor(double,double,double); 
  void set_growthFactor(std::vector<double>); 
  void resize(int); 
  double get_muMt1() const; 
  double get_muMt2() const; 
  double get_muA1() const; 
  double get_muA2() const; 
  double get_muX1() const; 
  double get_muX2() const; 
  double get_muY1() const; 
  double get_muY2() const; 
  int get_generation() const; 
  int get_nbHuman() const; 
  int get_nbFemale() const; 
  int get_nbMale() const; 
  double get_popsizeDemo() const; 
  int get_sizeMt() const;
  int get_sizeY() const;
  int get_sizePerX() const;
  int get_sizePerA() const;
  int get_crumbA() const;
  int get_crumbY() const;
  int get_crumbX() const;
  int get_crumbMt() const;
  int get_nbCellPerA() const;
  int get_nbCellPerX() const;
  int get_nbChromosomesA() const;
  int get_nbChromosomesX() const;
  std::vector<double> get_growthFactor() const; 
  std::vector<int> get_nb_mates(); 
  std::vector<int> get_nb_children(); 
  // *** Sub sample ***
  Population subsample(int sampleSize, int nbWomen); 
  Population subsample_female();  // fast and easy, if want to subsample a # of women use subsample(#,#)
  Population subsample_male();  // idem subsample(#,0)
  void analyse();  // TODO
  // *** Evolution functions ***
  void evolve(int gen); 
  void evolve_gene(int gen,std::string filename); 
  // *** Outputs ***
  void write_txt_pop(std::string const&) const; 
  void output_pop_size() const; 
  void write_nexus(std::string const& file_name) const;  // <>
  void write_fasta(std::string const& file_name) const;  // <>
  void output_genealogy(std::ofstream & fs);  // for the current generation output id of mother and father in a file, format : generation \t  mother \t father \n one line per indiv
  // private:
  void write_smartFormat(std::string const& filename) const;  // <>
  void output_diversity();
  void output_arlequin(std::string const& filename); 
  void output_nb_child(); 
  void output_nb_mates(); 
  void test(int); // dummy but <> can test on humanmt ror humanddna some stuff
  // *** Evolution functions ***
  int burnin(int gen,double facMu);  // <> // output the number of gen run under the burnin
  int burnin(double diversityTreshold,int gen,double facMu); // <> // output the number of gen run under the burnin
  void evolve_one(); 
  typeDNA polySites(int nuclear) const;  // <>
  typeDNA polySites_mtdna() const;  // <>
  double pairwise_distance_mtdna() const;  // <>
  double diversity_A() const;  // <>
  double diversity_Y() const;  // <>
  double diversity_X() const;  // <>
  double diversity_mtdna() const;  // <>

 private:
  void mutation(); 
  void reproduction(); 
  void birth_boy(int indexfemale,int indexmale,int nbBabyMale);  // <> 
  void birth_girl(int indexfemale,int indexmale,int nbBabyFemale);  // <>
  void selection(); 
  // *** Diversity estimators ***
};

template <class T>
void Population<T>::reinitiate_ID_table()
{
  int k = 0;
  for (std::vector<int>::iterator it = femaleIDtable.begin(); it != femaleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  k=nbFemale;
  for (std::vector<int>::iterator it = maleIDtable.begin(); it != maleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
}

template <class T>
void Population<T>::set_mating_system(int a)
{
  matingSystem = a;
}

template <class T>
int Population<T>::get_matingSystem() const
{ 
  return matingSystem;
}

template <class T>
void Population<T>::set_inbreeding(int a)
{
  inbreeding = a;
}

template <class T>
int Population<T>::get_inbreeding() const
{ 
  return inbreeding;
}

template <class T>
void Population<T>::set_mu(double a, double b, double c,double d)
{
  muMt1 = a;
  muMt2 = a;
  muX1 = b;
  muX2 = b;
  muY1 = c;
  muY2 = c;
  muA1 = d;
  muA2 = d;
}

template <class T>
void Population<T>::set_mu(double a, double b, double c,double d,double e, double f, double g,double h)
{
  muMt1 = a;
  muMt2 = b;
  muX1 = c;
  muX2 = d;
  muY1 = e;
  muY2 = f;
  muA1 = g;
  muA2 = h;
}

template <class T>
void Population<T>::set_growthFactor(double a,double b,double c)
{
  growthFactor[0] = a;
  growthFactor[1] = b;
  growthFactor[2] = c;
}

template <class T>
void Population<T>::set_growthFactor(std::vector<double> a)
{
  growthFactor = a;
}

template <class T>
double Population<T>::get_muA1() const
{ 
  return muA1;
}

template <class T>
double Population<T>::get_muA2() const
{ 
  return muA2;
}

template <class T>
double Population<T>::get_muX1() const
{ 
  return muX1;
}

template <class T>
double Population<T>::get_muX2() const
{ 
  return muX2;
}

template <class T>
double Population<T>::get_muY1() const
{ 
  return muY1;
}

template <class T>
double Population<T>::get_muY2() const
{ 
  return muY2;
}

template <class T>
double Population<T>::get_muMt1() const
{ 
  return muMt1;
}

template <class T>
double Population<T>::get_muMt2() const
{ 
  return muMt2;
}


template <class T>
int Population<T>::get_generation() const
{ 
  return generation;
}

template <class T>
int Population<T>::get_nbHuman() const
{
  return nbHuman;
}

template <class T>
int Population<T>::get_nbFemale() const
{
  return nbFemale;
}

template <class T>
int Population<T>::get_nbMale() const
{
  return nbHuman-nbFemale;
}

template <class T>
double Population<T>::get_popsizeDemo() const
{
  return popsizeDemo;
}

template <class T>
int Population<T>::get_crumbMt() const
{
  return nbCrumbMt;
}

template <class T>
int Population<T>::get_crumbA() const
{
  return nbCrumbA;
}

template <class T>
int Population<T>::get_crumbY() const
{
  return nbCrumbY;
}

template <class T>
int Population<T>::get_crumbX() const
{
  return nbCrumbX;
}

template <class T>
int Population<T>::get_nbCellPerX() const
{
  return nbCellPerX;
}

template <class T>
int Population<T>::get_nbCellPerA() const
{
  return nbCellPerA;
}

template <class T>
int Population<T>::get_nbChromosomesA() const
{
  return nbChromosomesA;
}

template <class T>
int Population<T>::get_nbChromosomesX() const
{
  return nbChromosomesX;
}

template <class T>
int Population<T>::get_sizeMt() const
{
  return sizeMt;
}

template <class T>
int Population<T>::get_sizePerA() const
{
  return sizePerA;
}

template <class T>
int Population<T>::get_sizePerX() const
{
  return sizePerX;
}

template <class T>
int Population<T>::get_sizeY() const
{
  return sizeY;
}

template <class T>
std::vector<int> Population<T>::get_nb_mates()
{
  int parentG = !(generation&1);
  std::vector<int> mates(nbHuman);
  for(int i=0; i<nbFemale; ++i)
    mates[i] = female[parentG][i].getNbMate(); // careful we have to look at the parent generation
  for(int i=0; i<(nbHuman-nbFemale); ++i)
    mates[i+nbFemale] = male[parentG][i].getNbMate();
  return mates;
}

template <class T>
std::vector<int> Population<T>::get_nb_children()
{
  int parentG = (generation&1)^1;
  std::vector<int> children(nbHuman);
  for(int i=0; i<nbFemale; ++i)
      children[i] = female[parentG][i].getNbChild(); // careful we have to look at the parent generation
  for(int i=0; i<(nbHuman-nbFemale); ++i)
    children[i+nbFemale] = male[parentG][i].getNbChild();
  return children;
}

template<class T> 
std::vector<double> Population<T>::get_growthFactor() const
{
  return growthFactor;
}

template <class T>
Population<T> Population<T>::subsample(int sampleSize, int nbWomen)
{
  // Subsample the Population for analysis
  // Create a new object of he class Population
  int currentG = generation&1;
  if(sampleSize>nbHuman)
    {
      std::cerr << "ERROR : attempt to subsample more individuals than the Population size\n";
      sampleSize = nbHuman;
    }
  if(nbWomen>nbFemale)
    {
      std::cerr << "ERROR : attempt to subsample more female than the Population size "<< nbWomen<< " versus" << nbFemale<< "\n";
      nbWomen = nbFemale;
    }  
  Population sample(sampleSize,sizeMt,nbChromosomesX,sizePerX,sizeY,nbChromosomesA,sizePerA);
  // NB : the Population sample will start at generation 0 -> must initiate at 0 and not currentG
  random_shuffle_pop();
  sample.nbFemale=nbWomen;
  sample.matingSystem=matingSystem;
  for(int i(0);i<nbWomen;++i)
    {
      sample.female[0][i].reborn(female[currentG][i]);
    }
  for(int i(0);i<sampleSize-nbWomen;++i)
    {      
      sample.male[0][i].reborn(male[currentG][i]);
    }
  return sample;
}
  
  template <class T>
Population<T> Population<T>::subsample_female()
{
  // Subsample the Population for analysis
  // Create a new object of he class Population
  int currentG = generation&1;
  Population sample(2*nbFemale,sizeMt,nbChromosomesX,sizePerX,sizeY,nbChromosomesA,sizePerA); // Default attribute for parameters - 2times female because we need female tab big enough
  sample.nbHuman = nbFemale; // then restore population size equals to female -> 'kill' the men
  for(int i(0);i<nbFemale;++i)
    {
      sample.female[0][i].reborn(female[currentG][i]);
      // To not resample twice the same individual we swap the sampled ones and check somewhere else
    }
  return sample;
}

template <class T>
Population<T> Population<T>::subsample_male()
{
  int currentG = generation&1;
  int nbMale = nbHuman-nbFemale;
  Population sample(2*nbMale,sizeMt,nbChromosomesX,sizePerX,sizeY,nbChromosomesA,sizePerA); // Default attribute for parameters -2times male because we need male tab big enough
  sample.nbHuman = nbMale; // then we 'kill' the female
  sample.nbFemale = 0;
  for(int i(0);i<(nbHuman-nbFemale);++i)
    {
      sample.male[0][i].reborn(male[currentG][i]);
      // To not resample twice the same individual we swap the sampled ones and check somewhere else
    }
  return sample;
}


template <class T>
void Population<T>::evolve(int gen)
{
  for(int i(0);i<gen;++i){
    if(nbHuman>0)
      evolve_one();
  }
  if(nbHuman==0)
    std::cout << "Population is extinct" << std::endl;
}

template <class T>
void Population<T>::evolve_gene(int gen,std::string filename)
{
  std::ofstream fs;
  fs.open(filename.c_str(),std::ios::trunc);
  for(int i(0);i<gen;++i){
    evolve_one();
    output_genealogy(fs);
  }
  fs.close();
}

template <class T>
void Population<T>::evolve_one()
{
  random_shuffle_pop();
  mutation();
  selection();
  reproduction();
  ++generation;
}


template <class T>
void Population<T>::reproduction()
{  
  int currentG = (generation&1);
  int nbMating(0),nbBabyToMake(floor(popsizeDemo));
  int randomNbChildren,nbBabyFemale(0),nbBabyMale(0);
  int indexmale(-1),indexfemale(-1);// force the segfault if not initialized
  int randomIndivSwap(0),testloop(0);
  bool sexChild;
  if((nbFemale>0)&&(nbHuman-nbFemale>0))
    { // if not no mating necessary population go to extinction
      switch(matingSystem)
    {
    case 1: // random monogamy
      {
        if(nbFemale>nbHuman/2){
          nbMating = nbHuman-nbFemale; }
        else{
          nbMating = nbFemale;} 
        for(int i(0);i<nbMating;++i) 
          {      
        indexmale = i;
        indexfemale = i;
        if(inbreeding==1)
          if(generation!=0){
            testloop = 0;
            while((female[currentG][indexfemale].shared_parents2(male[currentG][indexmale])) && (testloop<100))
              // WARNING : to prevent infinite loop -> test loop
              {
            randomIndivSwap = myrand(nbHuman-nbFemale-i);
            tmpmale.reborn(male[currentG][randomIndivSwap + i]);
            male[currentG][randomIndivSwap + i].reborn(male[currentG][indexmale]);
            male[currentG][indexmale].reborn(tmpmale);
            ++testloop;
              }}
        randomNbChildren = rBinom(nbBabyToMake,1/(double)(nbMating-i));
        if(TRACK_MATE_AND_CHILDREN)
          {
            female[currentG][indexfemale].newMate();
            male[currentG][indexmale].newMate();
            female[currentG][indexfemale].hadChild(randomNbChildren);
            male[currentG][indexmale].hadChild(randomNbChildren);
          }
        for(int j(0);j<randomNbChildren;++j)
          {
            sexChild = rand2;
            if(sexChild == true)
              {         
            birth_girl(indexfemale,indexmale,nbBabyFemale);
            ++nbBabyFemale;
              }
            else
              {
            birth_boy(indexfemale,indexmale,nbBabyMale);
            ++nbBabyMale;
              }     
            --nbBabyToMake;
          }
          }
        break;
      }
    case 2: // polygamy
      {
        nbMating = nbFemale; 
        for(int i(0);i<nbMating;++i) // first suppose we have as many male as female, not any more
          // if male and female have been well sorted (or randomly shuffled) gives us the proper way ...
          {      
        indexfemale = i;
        randomNbChildren = rBinom(nbBabyToMake,1/(double)(nbMating-i));
        if(TRACK_MATE_AND_CHILDREN)
          {
            female[currentG][indexfemale].hadChild(randomNbChildren);
          }
        for(int j(0);j<randomNbChildren;++j)
          {   
            {
              indexmale = myrand(nbHuman-nbFemale);
              if(inbreeding==1)
            if(generation!=0){
              testloop = 0;
              while((female[currentG][indexfemale].shared_parents2(male[currentG][indexmale])) && (testloop<100))
                {
                  randomIndivSwap = myrand(nbHuman-nbFemale);
                  tmpmale.reborn(male[currentG][randomIndivSwap + i]);
                  male[currentG][randomIndivSwap + i].reborn(male[currentG][indexmale]);
                  male[currentG][indexmale].reborn(tmpmale);
                  ++testloop;
                }}
              if(TRACK_MATE_AND_CHILDREN)
            {
              female[currentG][indexfemale].newMate();
              male[currentG][indexmale].newMate();
              male[currentG][indexmale].hadChild(1);
            }
              if(indexmale == nbHuman-nbFemale)
            {
              std::cerr << "RANDOM FAILURE" << std::endl;
              exit(2);
            }
              sexChild = rand2;
              if(sexChild == true)
            {
              birth_girl(indexfemale,indexmale,nbBabyFemale);
              ++nbBabyFemale;
            }
              else
            {
              birth_boy(indexfemale,indexmale,nbBabyMale);
              ++nbBabyMale;
            }       
            }
            --nbBabyToMake;
          }      
          }
        break;
      }
    case 3: // polygyny, male polygames but not women
      {
        nbMating = nbFemale; 
        for(int i(0);i<nbMating;++i) // first suppose we have as many male as female, not any more
          // if male and female have been well sorted (or randomly shuffled) gives us the proper way ...
          {
        indexfemale = i;          
        randomNbChildren = rBinom(nbBabyToMake,1/(double)(nbMating-i));
        if(TRACK_MATE_AND_CHILDREN)
          {
            female[currentG][indexfemale].newMate();
            female[currentG][indexfemale].hadChild(randomNbChildren);
          }
        for(int j(0);j<randomNbChildren;++j)
          {
            indexmale = myrand(nbHuman-nbFemale);
            if(inbreeding==1)
            if(generation!=0){
              testloop = 0;
              while((female[currentG][indexfemale].shared_parents2(male[currentG][indexmale])) && (testloop<100))
                {
                  randomIndivSwap = myrand(nbHuman-nbFemale-i);
                  tmpmale.reborn(male[currentG][randomIndivSwap + i]);
                  male[currentG][randomIndivSwap + i].reborn(male[currentG][indexmale]);
                  male[currentG][indexmale].reborn(tmpmale);
                  ++testloop;
                }}
            if(TRACK_MATE_AND_CHILDREN)
              {         
            male[currentG][indexmale].newMate();
            male[currentG][indexmale].hadChild(1);
              }
            sexChild = rand2;
            if(sexChild == true)
              {
            birth_girl(indexfemale,indexmale,nbBabyFemale);
            ++nbBabyFemale;
              }
            else
              {
            birth_boy(indexfemale,indexmale,nbBabyMale);
            ++nbBabyMale;
              }
            --nbBabyToMake;     
          }      
          }
        break;
      }
    case 4: // polyandry, female polygames but not men 
      { 
        nbMating = nbHuman - nbFemale;
        for(int i(0);i<nbMating;++i)
          {
        indexmale = i;
        randomNbChildren = rBinom(nbBabyToMake,1/(double)(nbHuman-nbFemale));
        if(TRACK_MATE_AND_CHILDREN)
          {
            male[currentG][indexmale].newMate();
            male[currentG][indexmale].hadChild(randomNbChildren);
          }
        for(int j(0);j<randomNbChildren;++j)
          {
            indexfemale = myrand(nbFemale);
            if(inbreeding==1)
              if(generation!=0){
            testloop = 0;
            while((female[currentG][indexfemale].shared_parents2(male[currentG][indexmale])) && (testloop<100))
              {
                randomIndivSwap = myrand(nbFemale);
                tmpfemale.reborn(female[currentG][randomIndivSwap + i]);
                female[currentG][randomIndivSwap + i].reborn(female[currentG][indexfemale]);
                female[currentG][indexfemale].reborn(tmpfemale);
                ++testloop;
              }}
            if(TRACK_MATE_AND_CHILDREN)
              {
            
            female[currentG][indexfemale].newMate();
            female[currentG][indexfemale].hadChild(1);
              }
            sexChild = rand2;
            if(sexChild == true)
              {
            birth_girl(indexfemale,indexmale,nbBabyFemale);
            ++nbBabyFemale;
              }
            else
              {
            birth_boy(indexfemale,indexmale,nbBabyMale);
            ++nbBabyMale;
              }
            --nbBabyToMake;
          } 
          }
        break;
      }
    case 5: // random mating
      {
        nbMating = std::floor(popsizeDemo); 
        for(int i(0);i<nbMating;++i)
          { // 1 baby at a time, as many as nneded to fill up to popsizeDemo
        indexfemale = myrand(nbFemale);
        indexmale = myrand(nbHuman-nbFemale);
        if(TRACK_MATE_AND_CHILDREN)
          {
            female[currentG][indexfemale].hadChild(1);
            male[currentG][indexmale].hadChild(1);
            female[currentG][indexfemale].newMate();
            male[currentG][indexmale].newMate();
          }
        if(indexmale == nbHuman-nbFemale)
          {
            std::cerr << "ERROR in the creation of the new generation" << std::endl;
            exit(2);
          }
        sexChild = rand2;
        if(sexChild == true)
          {
            birth_girl(indexfemale,indexmale,nbBabyFemale);
            ++nbBabyFemale;
          }
        else
          {
            birth_boy(indexfemale,indexmale,nbBabyMale);
            ++nbBabyMale;
          }     
        --nbBabyToMake;
          }
        break;
      }
    default:
      {
        std::cout << "Non existing mating system !" << std::endl;
        exit(2);
      }
    }    
      nbFemale = nbBabyFemale;
      nbHuman = nbBabyMale+nbBabyFemale;
    }
  else
    { // exctinction
      nbFemale = 0;
      nbHuman = 0;
    }
  if(nbHuman > popsizeDemo)
    {
      std::cerr << "ERROR : reproduction scheme is not working, too many baby were created" << std::endl;
      exit(2);
    }
  reinitiate_ID_table();
}

template <class T>
void Population<T>::selection() // in this step we just kill people to follow our demographic function - shoudl not append but safeguard
{
  if(nbHuman>popsizeDemo)
    {
      std::cerr << "TROUBLE Nb of children not working !" << std::endl;      
      int nbDead,nbFemaleDead;
      nbDead = nbHuman - std::floor(popsizeDemo); // nb people that we should kill to get that our pop is out growing
      nbFemaleDead = rBinom(nbDead,(double)nbFemale/(double)nbHuman); // random distribution of the dead between female & male, can have diff probability to die
      nbFemale -= nbFemaleDead;
      nbHuman -= nbDead;
    }
  // update popsizeDemo as a function of time -> demographic model
  popsizeDemo = demographic_function(popsizeDemo,growthFactor[0],growthFactor[1],growthFactor[2]);
}

template <class T>
Population<T>& Population<T>::operator = (Population<T> const& a)
{
  female[0] = a.female[0];
  male[0] = a.male[0];
  female[1] = a.female[1];
  male[1] = a.male[1];
  femaleIDtable = a.femaleIDtable;
  maleIDtable = a.maleIDtable;
  popsizeDemo = a.popsizeDemo;
  nbHuman = a.nbHuman;
  nbFemale = a.nbFemale;
  generation = a.generation;
  growthFactor = a.growthFactor;
  matingSystem = a.matingSystem;
  nbCrumbMt = a.nbCrumbMt;
  nbCrumbA = a.nbCrumbA;
  nbCrumbY = a.nbCrumbY;
  nbCrumbNuc = a.nbCrumbNuc;
  nbCellPerA = a.nbCellPerA;
  nbCellPerX = a.nbCellPerX;
  nbCrumbX = a.nbCrumbX;
  nbChromosomesA= a.nbChromosomesA;
  nbChromosomesX= a.nbChromosomesX;
  sizeMt= a.sizeMt;
  sizeY = a.sizeY;
  sizePerA = a.sizePerA;
  sizePerX = a.sizePerX;
  muA1 = a.muA1;
  muA2 = a.muA2;
  muX1 = a.muX1;
  muX2 = a.muX2;
  muY1 = a.muY1;
  muY2 = a.muY2;
  muMt1 = a.muMt1;
  muMt2 = a.muMt2;
  return *this;
}
  

template <class T>
std::ostream& Population<T>::print(std::ostream & os) const
{
  int currentG = (generation&1);
  os << "Pop with " << nbHuman << "individuals with a sex ratio of " << nbFemale/nbHuman << std::endl;
  for(int i(0);i<nbFemale;++i)
    {
      os << female[currentG][i];
    }
  for(int i(0);i<(nbHuman-nbFemale);++i)
    {
      os << male[currentG][i];
    }  
  return os;    
}


template <class T>
void Population<T>::write_txt_pop(std::string const& file_name) const
{
  int currentG = (generation&1);
  std::ofstream fs(file_name.std::string::c_str());
  fs << "generation " << generation << " -- nb humans " << nbHuman << " -- Sex ratio " << nbFemale/nbHuman << std::endl;
  if(nbHuman>0)
    {
      for(int i(0);i<nbFemale;++i)
    {
      fs << female[currentG][i];
    }
      for(int i(0);i<(nbHuman-nbFemale);++i)
    {
      fs << male[currentG][i];
    }  
      fs << "---" << std::endl;    
    }
}

template <class T>
void Population<T>::output_pop_size() const
{
  std::cout << "generation " << generation << " -- nb humans " << nbHuman << " -- Sex ratio " << (double)nbFemale/(double)nbHuman << std::endl;
}

template <class T>
void Population<T>::output_genealogy(std::ofstream & fs)
{
  int currentG = (generation&1);
  for(int i(0);i<nbFemale;++i)
    {
      fs << generation << "\t" << femaleIDtable[i] << "\t" << female[currentG][i].getMother() << "\t" << female[currentG][i].getFather() << "\n";
    }  
  for(int i(0);i<(nbHuman-nbFemale);++i)
    {
      fs << generation << "\t"  << maleIDtable[i] << "\t" << male[currentG][i].getMother() << "\t" << male[currentG][i].getFather() << "\n";
    }  
}

template <class T>
void Population<T>::initiate_pop(int N)
{
  nbHuman = N;
  generation = 0;
  nbFemale = N/2;
  popsizeDemo = N;
  matingSystem = 0;
  inbreeding = 0;
  muA1 = 0;
  muA2 = 0;
  muX1 = 0;
  muX2 = 0;
  muY1 = 0;
  muY2 = 0;
  muMt1 = 0;
  muMt2 = 0;
  nbCrumbMt = 0;
  nbCrumbA = 0;
  nbCrumbY = 0;
  nbCrumbX = 0;
  nbCrumbNuc = 0;
  nbChromosomesX = 0;
  nbChromosomesA = 0;
  nbCellPerA = 0;
  nbCellPerX = 0;
  sizeMt = 0;
  sizePerX = 0;
  sizeY = 0;
  sizePerA = 0;
  growthFactor.resize(3);
  growthFactor[0]=0;
  growthFactor[1]=1; // by default popsize (t+1) = popsize (t)
  growthFactor[2]=0;
  femaleIDtable.resize(std::ceil(popsizeDemo));
  maleIDtable.resize(std::ceil(popsizeDemo));
  int k = 0;
  for (std::vector<int>::iterator it = femaleIDtable.begin(); it != femaleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  k=nbFemale;
  for (std::vector<int>::iterator it = maleIDtable.begin(); it != maleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  T tmp(nbCrumbMt,nbCrumbNuc);
  tmpfemale = tmp;
  T tmp2(tmp);
  tmpmale = tmp2;
  for(int i(0);i<popsizeDemo;++i)
    {
      T people(tmpfemale);
      female[0].push_back(people);
      T people2(tmpfemale);
      female[1].push_back(people2);
    }
  for(int i(0);i<popsizeDemo;++i)
    {
      T people(tmpmale);
      male[0].push_back(people);
      T people2(tmpmale);
      male[1].push_back(people2);
    }
}

template <class T>
void Population<T>::reset(int N)
{
  nbHuman = N;
  generation = 0;
  nbFemale = N/2;
  popsizeDemo = N;
  int k = 0;
  for (std::vector<int>::iterator it = femaleIDtable.begin(); it != femaleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  k=nbFemale;
  for (std::vector<int>::iterator it = maleIDtable.begin(); it != maleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  for(int i(0);i<std::ceil(N/2.0);++i)
    {      
      male[0][i].reset();
    }
  for(int i(0);i<N/2;++i)
    {      
      female[0][i].reset();
    }
}

template <class T>
void Population<T>::initiate_pop(int N,int mt,int xL, int xS,int y,int aL,int aS)
{
  nbHuman = N;
  generation = 0;
  nbFemale = N/2;
  popsizeDemo = N;
  matingSystem = 0;
  inbreeding = 0;
  muA1 = 0;
  muA2 = 0;
  muX1 = 0;
  muX2 = 0;
  muY1 = 0;
  muY2 = 0;
  muMt1 = 0;
  muMt2 = 0;
  nbCrumbMt = std::ceil(mt*1.0/CRUMB);
  sizeMt = std::ceil(mt*1.0/CRUMB)*CRUMB;
  sizePerA = std::ceil(aS*1.0/CRUMB)*CRUMB; // forced to be long as a 32 bits chunks
  nbChromosomesA = aL * 2;
  nbCellPerA = sizePerA/CRUMB;
  nbCrumbA = nbChromosomesA * sizePerA /CRUMB;
  sizeY = std::ceil(y*1.0/CRUMB)*CRUMB;
  nbCrumbY = sizeY/CRUMB;
  sizePerX = std::ceil(xS*1.0/CRUMB)*CRUMB;
  nbCellPerX = sizePerX /CRUMB;
  nbChromosomesX = xL;
  nbCrumbX = nbChromosomesX * sizePerX /CRUMB;
  int nbCrumbSex = std::max(nbCrumbX,nbCrumbY);
  nbCrumbNuc = nbCrumbA + nbCrumbX  + nbCrumbSex;
  growthFactor.resize(3);
  growthFactor[0]=0;
  growthFactor[1]=1; // by default popsize = popsize
  growthFactor[2]=0;
  femaleIDtable.resize(std::ceil(popsizeDemo));
  maleIDtable.resize(std::ceil(popsizeDemo));
  int k = 0;
  for (std::vector<int>::iterator it = femaleIDtable.begin(); it != femaleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  k=nbFemale;
  for (std::vector<int>::iterator it = maleIDtable.begin(); it != maleIDtable.end(); ++it)
    {
      *it=k;
      k++;
    }
  T tmp(nbCrumbMt,nbCrumbNuc);
  tmpfemale = tmp;
  T tmp2(tmp);
  tmpmale = tmp2;
  for(int i(0);i<popsizeDemo;++i)
    {
      T people(tmpfemale);
      female[0].push_back(people);
      T people2(tmpfemale);
      female[1].push_back(people2);
    }
  for(int i(0);i<popsizeDemo;++i)
    {
      T people(tmpmale);
      male[0].push_back(people);
      T people2(tmpmale);
      male[1].push_back(people2);
    }
}

template <class T>
void Population<T>::resize(int a)
{
  int oldsize = male[0].size();
  for(int i(0); i<(a-oldsize); i++)
    {
      T newH0(tmpmale);
      male[0].push_back(newH0);
      T newH1(tmpmale);
      male[1].push_back(newH1);
      T newF0(tmpfemale);
      female[0].push_back(newF0);
      T newF1(tmpfemale);
      female[1].push_back(newF1);
    }
  femaleIDtable.resize(a);
  maleIDtable.resize(a);
}

template <class T>
void Population<T>::random_shuffle_pop() // knuth fisher yates algorithm stolen on http://www.codinghorror.com/blog/2007/12/the-danger-of-naivete.html
{
  int currentG = (generation%2);
  int n;
  int tmpIndex = -1;
  for (int i = nbFemale - 1; i > 0; --i)
    {
      n = myrand(i-1);
      tmpfemale.reborn(female[currentG][i]);
      female[currentG][i].reborn(female[currentG][n]);
      female[currentG][n].reborn(tmpfemale);
      tmpIndex = femaleIDtable[i];
      femaleIDtable[i] = femaleIDtable[n];
      femaleIDtable[n] = tmpIndex;
    }
  for (int i = nbHuman - nbFemale - 1; i > 0; --i)
    {
      n = myrand(i-1);
      tmpmale.reborn(male[currentG][i]);
      male[currentG][i].reborn(male[currentG][n]);
      male[currentG][n].reborn(tmpmale);
      tmpIndex = maleIDtable[i];
      maleIDtable[i] = maleIDtable[n];
      maleIDtable[n] = tmpIndex;
    }      
}

template <class T>
void Population<T>::output_diversity()
{  
  std::cout <<  generation << "\t" << generation << "\t" << diversity_mtdna() <<  "\t"<< pairwise_distance_mtdna() << " \t " << diversity_A() << "\t" << diversity_X() << "\t" << diversity_Y() << std::endl;
}

template <class T>
void Population<T>::output_nb_child()
{
  std::vector<int> nbChildren = get_nb_children();
  for(int i(nbFemale);i<nbHuman-nbFemale;++i)
    std::cout << "\n" << nbChildren[i];
  std::cout<<std::endl;
}

template <class T>
void Population<T>::output_nb_mates()
{
  std::vector<int> nbMate = get_nb_mates();
  for(int i(0);i<nbFemale;++i)
    std::cout<<"\n"<<nbMate[i];
  std::cout<<std::endl;
}


#endif