2 #ifndef __SEQUENCE_COALESCENT_BITS_TRAJECTORIES_TCC__
3 #define __SEQUENCE_COALESCENT_BITS_TRAJECTORIES_TCC__
11 template< typename uni01_generator >
12 void ConditionalTraj_details(uni01_generator & uni01,
13 std::vector<double> * traj,
17 const double & initial_frequency,
18 const double & final_frequency)
20 Implementation details
23 assert( initial_frequency > 0. );
24 assert( final_frequency > 0. );
25 assert( final_frequency <= 1. );
27 double x = initial_frequency;
28 traj->erase(traj->begin(),traj->end());
30 while( x <= final_frequency )
32 double ux = (double(2*N)*s*x*(1.-x))/tanh(double(2*N)*s*x);
35 x += ( ux*dt + std::sqrt( x*(1.-x)*dt ) );
39 x += ( ux*dt - std::sqrt( x*(1.-x)*dt ) );
41 assert( std::isfinite(x) );
42 traj->push_back( ( (x<final_frequency)
43 ? x : final_frequency ) );
45 std::reverse(traj->begin(),traj->end());
48 template<typename uni01_generator>
49 void ConditionalTrajNeutral_details(uni01_generator & uni01,
50 std::vector<double> * traj,
52 const double & initial_freq,
53 const double & final_freq)
55 Implementation details
58 double x = initial_freq;
59 traj->erase(traj->begin(),traj->end());
61 while( x >= final_freq )
66 x += ( ux*dt + std::sqrt( x*(1.-x)*dt ) );
70 x += ( ux*dt - std::sqrt( x*(1.-x)*dt ) );
72 assert( std::isfinite(x) );
73 traj->push_back( ( (x>final_freq) ? x : final_freq ) );
77 template< typename uni01_generator >
78 void ConditionalTraj(uni01_generator & uni01,
79 std::vector<double> * traj,
83 const double & initial_frequency,
84 const double & final_frequency)
86 Stochastic trajectory of beneficial mutations, following
87 Coop and Griffiths (2004).
88 \param uni01 a random number generator returning a U(0,1].
89 \param traj For a diploid population of size N, this function will return trajectories
90 equivalent to what one would get from a Wright-Fisher simulation of a haploid
92 \param dt amount by which to change increment time during simulation.
93 \param initial_frequency Initial frequency of beneficial allele (i.e. 1/2N)
94 \param final_frequency Final frequency of beneficial allele (1 means fixation).
95 \return A vector of length L, such that,for dt = 1/(k*2N),
96 L/(k*2N) is the length of the sweep, in units of 2N generations
99 ConditionalTraj_details(uni01,traj,N,s,dt,initial_frequency,final_frequency);
102 template< typename uni01_generator >
103 void ConditionalTraj(const uni01_generator & uni01,
104 std::vector<double> * traj,
108 const double & initial_frequency,
109 const double & final_frequency)
111 Stochastic trajectory of beneficial mutations, following
112 Coop and Griffiths (2004).
113 \param uni01 a random number generator returning a U(0,1].
114 \param traj A vector of length L, such that,for dt = 1/(k*2N),
115 L/(k*2N) is the length of the sweep, in units of 2N generations
116 \note For a diploid population of size N, this function will return trajectories
117 equivalent to what one would get from a Wright-Fisher simulation of a haploid
118 population of size 2N
119 \param dt amount by which to change increment time during simulation.
120 \param initial_frequency Initial frequency of beneficial allele (i.e. 1/2N)
121 \param final_frequency Final frequency of beneficial allele (1 means fixation).
124 ConditionalTraj_details(uni01,traj,N,s,dt,initial_frequency,final_frequency);
127 template<typename uni01_generator>
128 void ConditionalTrajNeutral(uni01_generator & uni01,
129 std::vector<double> * traj,
131 const double & initial_freq,
132 const double & final_freq)
135 Stochastic trajectory of a neutral allele, following Coop & Griffiths (2004) TPB,
136 and Przeworski et al. (2005) Evolution. The simulation is backwards in time.
137 \param uni01 a random number generator returning a U(0,1].
138 \param traj A vector of length L, and for dt=1/(k*2N), the length of time,
139 in units of 2N generations, is given by L/(k*2N). The vector describes
140 the change in allele frequency from \a initial_freq to \a final_freq,
141 in jumps in time of \a dt.
142 \param dt amount by which to change increment time during simulation.
143 \param initial_freq Initial frequency of the neutral allele.
144 \param final_freq final frequency of the neutral allele.
147 ConditionalTrajNeutral_details(uni01,traj,dt,initial_freq,final_freq);
150 template<typename uni01_generator>
151 void ConditionalTrajNeutral(const uni01_generator & uni01,
152 std::vector<double> * traj,
154 const double & initial_freq,
155 const double & final_freq)
157 Stochastic trajectory of a neutral allele, following Coop & Griffiths (2004) TPB,
158 and Przeworski et al. (2005) Evolution. The simulation is backwards in time.
159 \param uni01 a random number generator returning a U(0,1].
160 \param traj A vector of length L, and for dt=1/(k*2N), the length of time,
161 in units of 2N generations, is given by L/(k*2N). The vector describes
162 the change in allele frequency from \a initial_freq to \a final_freq,
163 in jumps in time of \a dt.
164 \param dt amount by which to change increment time during simulation.
165 \param initial_freq Initial frequency of the neutral allele.
166 \param final_freq final frequency of the neutral allele.
169 ConditionalTrajNeutral_details(uni01,traj,dt,initial_freq,final_freq);