genetics.h

#ifndef DEF_GENETICS
#define DEF_GENETICS
#include "boost/multi_array.hpp"
#include <cassert>
#include "population.h"

// file to make genetic stuff


typedef boost::multi_array<typeDNA, 2> array_type;
// site frequency spectrum SFS

class DNA
{
 public:
  // *** constructor destructor *
  DNA(); 
  template <class T> 
    explicit DNA(Population<T> *pop)
    {
      N = pop -> get_nbHuman();
      nbXcell = pop->get_crumbMt();
      nbS = pop->get_sizeMt();
      assert(nbXcell*32 == nbS);
      if(N==1)
    tab.resize(boost::extents[2][nbXcell]); // pb with multui_array
      else
    tab.resize(boost::extents[N][nbXcell]);
      int nbFemale = pop -> get_nbFemale();
      int currentG = pop -> get_generation();
      currentG &= 1;
      for(int i(0); i<N; ++i)
    {
      for(int j(0) ;j<nbXcell ;++j)
        {
          if(i < nbFemale)
        tab[i][j] = pop -> female[currentG][i].mtDna[j];
          else
        tab[i][j] = pop -> male[currentG][i-nbFemale].mtDna[j];
        }
    }  
    }
  DNA(Population<Humandna>* pop, int tdna); 
  DNA(DNA const & a); 
  DNA(DNA const & a, DNA const & b); 
  DNA(typeDNA* input,int size);
  DNA(array_type const & input,int sizeN, int sizeX);
  DNA(std::vector< std::vector<typeDNA> > input);
  //  DNA(typeDNA* input);
  ~DNA();
  friend std::ostream &operator<<(std::ostream&,DNA const&);
  friend DNA operator+(DNA const&,DNA const&); 
  DNA& operator=(DNA const&);
  int get_N() const; 
  int get_nbXcell() const; 
  int get_nbS() const;   
  long poly_sites() const; 
  double prop_poly_sites() const; 
  double mean_pw_diff() const; 
  double theta_waterson() const; //** Return Watterson's Theta for the population */
  double theta_hom() const; //** Return a theta computed from the homozygiosity */
  double theta_pi() const; 
  int nb_of_haplotypes() const; 
  DNA haplotypes() const; 
  std::vector<double> hap_frequency() const; 
  std::vector<double> site_frequency(int site) const; 
  double heterozigosity_hap() const; 
  double heterozigosity_sites() const; 
  double Neff(double theta, double mu) const; 
  double tajima_D() const; 
  std::vector<double> site_frequency_DNA(int site) const; 
  typeDNA dist_two(int a, int b) const; 
  void SFS(std::string filename) const; 
  void SFS_DNA(std::string filename) const; 
  void output_diversity(std::ofstream&) const; 
  void output_diversity_with_time(std::ofstream&,int Ntheo, int timeG) const; 
  void output_diversity(std::ofstream & fs, int Ntheo) const; 
  void output_arlequin(std::string filename) const; 
  void output_fasta(std::string filename) const; 
 private:
  std::ostream & print(std::ostream&) const;
  int N; 
  int nbXcell; 
  int nbS; 
  array_type tab; 
};
std::vector< std::vector<typeDNA> >  input_fasta(std::string filename);
void output_diversity_headline(std::ofstream &); 
void output_diversity_headline_with_time(std::ofstream &); 
#endif