2 #include <Sequence/Coalescent/Coalesce.hpp> 4 #include <Sequence/Coalescent/Coalesce.hpp> 5 #include <Sequence/Coalescent/Recombination.hpp> 6 #include <Sequence/Coalescent/Mutation.hpp> 13 template<
typename uniform_generator,
14 typename uniform01_generator,
15 typename exponential_generator,
16 typename poisson_generator>
18 uniform01_generator & uni01,
19 exponential_generator & expo,
20 poisson_generator & poiss,
25 std::vector<chromosome> * sample,
27 unsigned * max_chromosomes = NULL,
28 const unsigned & max_chromosomes_inc = 0)
66 double littler = rho/(double(nsites-1));
71 int nlinks = nsam*(nsites-1);
75 double rcoal = double(NSAM*(NSAM-1));
76 double rrec = (rho>0.) ? littler*double(nlinks) : 0.;
79 double tcoal = expo(1./rcoal);
80 double trec = expo(1./rrec);
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() ) );
91 assert( (sample->begin()+pos_rec.first)->links()>0 );
92 nlinks -=
crossover(NSAM,pos_rec.first,pos_rec.second,
93 sample,sample_history);
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);
103 if (
unsigned(NSAM) < sample->size()/5)
105 sample->erase(sample->begin()+NSAM+1,sample->end());
108 if (max_chromosomes != NULL && sample->size() > *max_chromosomes)
109 *max_chromosomes += max_chromosomes_inc;
120 SimData gametes_obj = infinite_sites_sim_data(poiss,uni,
121 nsites,*sample_history,theta);
std::pair< int, int > pick_uniform_spot(const double &random_01, const int &nlinks, std::vector< chromosome >::const_iterator sample_begin, const unsigned ¤t_nsam)
Pick a crossover point for the model where recombination rates are constant across a recion...
std::list< marginal > arg
Ancestral Recombination Graph.
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 ¤t_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 ¤t_nsam, const int &chromo, const int &pos, std::vector< chromosome > *sample, arg *sample_history)
Recombination function.