simulation.h

/***  SMARTPOP, Simulating Mating Alliances as a Reproductive Tactic for POPuplations
    Copyright (C) 2013 E.G.Guillot

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
    
    Contact: e.guillot@massey.ac.nz
***/

#ifndef DEF_SIMULATION
#define DEF_SIMULATION
#include "genetics.h"


class Simulation
{
 public:
  Simulation(); 
  explicit Simulation(unsigned randomSeed); 
  Simulation(Simulation const &); 
  ~Simulation(); 
  Simulation& operator = (Simulation const&);
  // main operations
  int run(); 
  int print_help(); 
  void load_file(std::string); 
  void load_directory(std::string); 
  void export_file(std::string filename); 
  //  get /set
  void set_seed(int); 
  void set_popsize(int); 
  void set_sample(int); 
  void set_step(int); 
  void set_nstep(int); 
  void set_matsys(int); 
  void set_inbreeding(int); 
  void set_nSimu(int); 
  void set_sizeMt(int);
  void set_sizeY(int);
  void set_sizePerX(int);
  void set_sizePerA(int);
  void set_nbLociA(int);
  void set_nbLociX(int);
  void set_mu(double,double); 
  void set_mu(double,double,double,double); 
  void set_mu(double,double,double,double,double,double,double,double); 
  void set_growthFactor(double,double,double); 
  void set_burninT(double); 
  void set_dnaComplexity(int a); 
  void set_diversityX(int a); 
  void set_fileoutput(std::string); 
  void set_verbose_on(); 
  void set_fasta_on(); 
  void set_save_on(std::string dirName); 
  void set_arlequin_on(); 
  void set_genealogy_on(); 
  void set_header_on(); 
  int get_nSimu(); 
  bool get_verbose(); 
  void header_write(); 
  void check_parameters(); 
 private:
  unsigned seed; 
  int popsize; 
  int sample; 
  int step; 
  int nstep; 
  int matsys; 
  int inbreeding; 
  int nSimu; 
  int counter; 
  int sizeMt;
  int sizeY;
  int sizePerA;
  int sizePerX;
  int nbLociA;
  int nbLociX;
  double muMt1; 
  double muMt2; 
  double muX1; 
  double muX2; 
  double muY1; 
  double muY2; 
  double muA1; 
  double muA2; 
  std::vector<double> growthFactor; 
  double burninT; 
  double dnaComplexity;  // if 0 complete DNA,  if 1, mitochondrial, if 2 only human (-> genealogy)
  std::vector<int> diversityX;  // if 0 complete DNA, if 1 mtDNA, if 3 Y chr, if 2 X, if 4 autosome
  std::string fileoutput; 
  std::string directoryLoad; 
  std::string directorySave; 
  // booleans
  bool verbose;
  bool genealogy;
  bool fasta; 
  bool smartSave;
  bool arlequin;
  bool header;
  bool load;
  double maxpop; 
  // hidden functions
  int classic(); 
  int successive_classic(); 
  int mito(); 
  int successive_mito(); 
  int indiv_nodna(); 
  int successive_indiv_nodna(); 
  int output_main(int X, Population<Humandna> & epsilon); 
  int output_main(int X, Population<Humanmt> & epsilon); 
  int output_main(int X, Population<Human> & epsilon); 
  int verbose_start(); 
  int verbose_end(); 
  template <class T>
    void save(Population<T> epsilon); 
  template <class T>
    void evolve(Population<T> &epsilon); 
};

template <class T>
void Simulation::save(Population<T> epsilon)
{
  epsilon.write_smartFormat(directorySave + "_" + to_string(counter) + ".sim");
}

template <class T>
void Simulation::evolve(Population<T>& epsilon)
{
  std::string fileGenealogy = fileoutput + ".gene";
  epsilon.set_mating_system(matsys);
  epsilon.set_growthFactor(growthFactor);
  if(burninT>0)
    {
      epsilon.burnin(burninT,100,10);
    }
  for(int j(0);j<nstep;++j)
    {
      if(genealogy)
    epsilon.evolve_gene(step,fileGenealogy);
      else
    epsilon.evolve(step);
      for (unsigned k(0); k<diversityX.size();++k)
    if(diversityX[k]>-1)
      output_main(diversityX[k],epsilon);
    }
  if (smartSave)
    {
      save(epsilon);
    }
}

#endif