libsequence  1.9.5
Trajectories.tcc
1 // -*- C++ -*-
2 #ifndef __SEQUENCE_COALESCENT_BITS_TRAJECTORIES_TCC__
3 #define __SEQUENCE_COALESCENT_BITS_TRAJECTORIES_TCC__
4 #include <cassert>
5 #include <cmath>
6 #include <algorithm>
7 
8 namespace Sequence
9 {
10  namespace coalsim {
11  template< typename uni01_generator >
12  void ConditionalTraj_details(uni01_generator & uni01,
13  std::vector<double> * traj,
14  const unsigned & N,
15  const double & s,
16  const double & dt,
17  const double & initial_frequency,
18  const double & final_frequency)
19  /*!
20  Implementation details
21  */
22  {
23  assert( initial_frequency > 0. );
24  assert( final_frequency > 0. );
25  assert( final_frequency <= 1. );
26  assert( s>=0. );
27  double x = initial_frequency;
28  traj->erase(traj->begin(),traj->end());
29  traj->push_back(x);
30  while( x <= final_frequency )
31  {
32  double ux = (double(2*N)*s*x*(1.-x))/tanh(double(2*N)*s*x);
33  if( uni01() < 0.5 )
34  {
35  x += ( ux*dt + std::sqrt( x*(1.-x)*dt ) );
36  }
37  else
38  {
39  x += ( ux*dt - std::sqrt( x*(1.-x)*dt ) );
40  }
41  assert( std::isfinite(x) );
42  traj->push_back( ( (x<final_frequency)
43  ? x : final_frequency ) );
44  }
45  std::reverse(traj->begin(),traj->end());
46  }
47 
48  template<typename uni01_generator>
49  void ConditionalTrajNeutral_details(uni01_generator & uni01,
50  std::vector<double> * traj,
51  const double & dt,
52  const double & initial_freq,
53  const double & final_freq)
54  /*!
55  Implementation details
56  */
57  {
58  double x = initial_freq;
59  traj->erase(traj->begin(),traj->end());
60  traj->push_back(x);
61  while( x >= final_freq )
62  {
63  double ux = -1.*x;
64  if( uni01() < 0.5 )
65  {
66  x += ( ux*dt + std::sqrt( x*(1.-x)*dt ) );
67  }
68  else
69  {
70  x += ( ux*dt - std::sqrt( x*(1.-x)*dt ) );
71  }
72  assert( std::isfinite(x) );
73  traj->push_back( ( (x>final_freq) ? x : final_freq ) );
74  }
75  }
76 
77  template< typename uni01_generator >
78  void ConditionalTraj(uni01_generator & uni01,
79  std::vector<double> * traj,
80  const unsigned & N,
81  const double & s,
82  const double & dt,
83  const double & initial_frequency,
84  const double & final_frequency)
85  /*!
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
91  population of size 2N
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
97  */
98  {
99  ConditionalTraj_details(uni01,traj,N,s,dt,initial_frequency,final_frequency);
100  }
101 
102  template< typename uni01_generator >
103  void ConditionalTraj(const uni01_generator & uni01,
104  std::vector<double> * traj,
105  const unsigned & N,
106  const double & s,
107  const double & dt,
108  const double & initial_frequency,
109  const double & final_frequency)
110  /*!
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).
122  */
123  {
124  ConditionalTraj_details(uni01,traj,N,s,dt,initial_frequency,final_frequency);
125  }
126 
127  template<typename uni01_generator>
128  void ConditionalTrajNeutral(uni01_generator & uni01,
129  std::vector<double> * traj,
130  const double & dt,
131  const double & initial_freq,
132  const double & final_freq)
133 
134  /*!
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.
145  */
146  {
147  ConditionalTrajNeutral_details(uni01,traj,dt,initial_freq,final_freq);
148  }
149 
150  template<typename uni01_generator>
151  void ConditionalTrajNeutral(const uni01_generator & uni01,
152  std::vector<double> * traj,
153  const double & dt,
154  const double & initial_freq,
155  const double & final_freq)
156  /*!
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.
167  */
168  {
169  ConditionalTrajNeutral_details(uni01,traj,dt,initial_freq,final_freq);
170  }
171  }
172 }
173 #endif