libsequence  1.9.5
NeutralSample.hpp
2 #include <Sequence/Coalescent/Coalesce.hpp>
4 #include <Sequence/Coalescent/Coalesce.hpp>
5 #include <Sequence/Coalescent/Recombination.hpp>
6 #include <Sequence/Coalescent/Mutation.hpp>
7 #include <Sequence/SimData.hpp>
8 #include <utility>
9 
10 namespace Sequence
11 {
12  namespace coalsim {
13  template<typename uniform_generator,
14  typename uniform01_generator,
15  typename exponential_generator,
16  typename poisson_generator>
17  Sequence::SimData neutral_sample( uniform_generator & uni,
18  uniform01_generator & uni01,
19  exponential_generator & expo,
20  poisson_generator & poiss,
21  const double & theta,
22  const double & rho,
23  const int & nsites,
24  const int & nsam,
25  std::vector<chromosome> * sample,
26  arg * sample_history,
27  unsigned * max_chromosomes = NULL,
28  const unsigned & max_chromosomes_inc = 0)
62  {
63  int NSAM = nsam;
64 
65  //this is rho = 4Nr/site
66  double littler = rho/(double(nsites-1));
67 
68  //a chromosome with nsites sites has nsites-1 positions ("links")
69  //at which crossovers can occur, so the total number of
70  //links at the start of the simulation is:
71  int nlinks = nsam*(nsites-1);
72  double t = 0.;
73  while(NSAM>1)
74  {
75  double rcoal = double(NSAM*(NSAM-1));
76  double rrec = (rho>0.) ? littler*double(nlinks) : 0.;
77 
78  //note--the function calls below scale time in units of 4Ne
79  double tcoal = expo(1./rcoal);
80  double trec = expo(1./rrec);
81  if ( trec < tcoal ) //crossover event
82  {
83  t+=trec;
84  std::pair<int,int> pos_rec = pick_uniform_spot(uni01(),
85  nlinks,
86  sample->begin(),NSAM);
87  assert( pos_rec.second >= 0 );
88  assert( pos_rec.second >= (sample->begin()+pos_rec.first)->first() );
89  assert( pos_rec.second <= ((sample->begin()+pos_rec.first)->last() ) );
90 
91  assert( (sample->begin()+pos_rec.first)->links()>0 );
92  nlinks -= crossover(NSAM,pos_rec.first,pos_rec.second,
93  sample,sample_history);
94  NSAM++;
95  }
96  else //common ancestor event
97  {
98  t+=tcoal;
99  std::pair<int,int> two = pick2(uni,NSAM);
100  NSAM -= coalesce(t,nsam,NSAM,two.first,two.second,nsites,
101  &nlinks,sample,sample_history);
102  }
103  if (unsigned(NSAM) < sample->size()/5)
104  {
105  sample->erase(sample->begin()+NSAM+1,sample->end());
106  }
107  }
108  if (max_chromosomes != NULL && sample->size() > *max_chromosomes)
109  *max_chromosomes += max_chromosomes_inc;
110 
111  //As we have scaled time in units of 4Nr generations, we pass theta
112  //to the mutation function. If we had used the following code to
113  //generate times:
114  // double tcoal = expo(2./rcoal);
115  // double trec = expo(2./rrec);
116  //Then time would be scaled in units of 2Ne generations,
117  //and if our simulation considers theta=4Neu, we'd pass theta/2
118  //to the infinite_sites routine
119 
120  SimData gametes_obj = infinite_sites_sim_data(poiss,uni,
121  nsites,*sample_history,theta);
122  return gametes_obj;
123  }
124  }
125 }
std::pair< int, int > pick_uniform_spot(const double &random_01, const int &nlinks, std::vector< chromosome >::const_iterator sample_begin, const unsigned &current_nsam)
Pick a crossover point for the model where recombination rates are constant across a recion...
std::list< marginal > arg
Ancestral Recombination Graph.
Definition: SimTypes.hpp:217
Sequence::SimData neutral_sample(uniform_generator &uni, uniform01_generator &uni01, exponential_generator &expo, poisson_generator &poiss, const double &theta, const double &rho, const int &nsites, const int &nsam, std::vector< chromosome > *sample, arg *sample_history, unsigned *max_chromosomes=NULL, const unsigned &max_chromosomes_inc=0)
A simple function to generate samples under a neutral equilibrium model.
The namespace in which this library resides.
declaration of types for coalescent simulation
Declaration of Sequence::SimData, a class representing polymorphism data from coalescent simulations ...
Data from coalescent simulations.
int coalesce(const double &time, const int &ttl_nsam, const int &current_nsam, const int &c1, const int &c2, const int &nsites, int *nlinks, std::vector< chromosome > *sample, arg *sample_history)
Common ancestor routine for coalescent simulation. Merges chromosome segments and updates marginal tr...
int crossover(const int &current_nsam, const int &chromo, const int &pos, std::vector< chromosome > *sample, arg *sample_history)
Recombination function.