libsequence  1.9.5
Recombination.tcc
1 // -*- C++ -*-
2 #ifndef __SEQUENCE_COALESCENT_BITS_RECOMBINATION_TCC__
3 #define __SEQUENCE_COALESCENT_BITS_RECOMBINATION_TCC__
4 #include <utility>
5 #include <numeric>
6 
7 namespace Sequence
8 {
9  namespace coalsim {
10 #ifndef DOXYGEN_SKIP
11  template<typename uniform01_generator>
12  std::pair<int,int> pick_spot_details(uniform01_generator & uni01,
13  const double & total_reclen,
14  const std::vector<double> & reclens,
15  std::vector<chromosome>::const_iterator sample_begin,
16  const unsigned & current_nsam,
17  const double * rec_map)
18  {
19  double sum=0.;
20  int recombinant=0;
21  int pos=0;
22  double ran = uni01();
23  bool flag = false;
24  while( (sample_begin+recombinant) <
25  (sample_begin+current_nsam) && ! flag)
26  {
27  sum += reclens[std::vector<double>::size_type(recombinant)];
28  if ( ran <= sum/total_reclen )
29  {
30  int beg = (sample_begin+recombinant)->segs->beg;
31  int end = ((sample_begin+recombinant)->segs+(sample_begin+recombinant)->nsegs-1)->end;
32  double cum_prob = std::accumulate(rec_map+beg,rec_map+end,0.);
33  double ran2 = uni01(),sum2=0.;
34  //we only want to simulate crossover spots
35  //from the density that covers the interval beg to end-1
36  for(int i=beg ; i < end ; ++i)
37  {
38  sum2+= *(rec_map+i);
39  if ( ran2 <= sum2/cum_prob )
40  {
41  pos = i;
42  flag = true;
43  break;
44  }
45  }
46  }
47  if(flag) break;
48  ++recombinant;
49  }
50  return std::make_pair(recombinant,pos);
51  }
52 #endif
53 
54  template<typename uniform01_generator>
55  std::pair<int,int> pick_spot( uniform01_generator & uni01,
56  const double & total_reclen,
57  const std::vector<double> & reclens,
58  std::vector<chromosome>::const_iterator sample_begin,
59  const unsigned & current_nsam,
60  const double * rec_map)
61  /*!
62  Picks a positions amongst all chromosomes at which a recombination event
63  will occur, based on an arbitrary genetic map
64  \param uni01 a function/object which takes no arguments and can return a U[0,1)
65  \param total_reclen the total recombination length of all chromosomes in the sample
66  \param reclens a vector of the proportion of \a total_reclen contributed by each chromosome.
67  This needs to be ordered in the same order as \a sample_begin to (\a sample_begin + \a current_nsam - 1)
68  \param sample_begin an iterator pointing to the beginning of the sample
69  \param current_nsam the current sample size in the simulation
70  \param rec_map an array of probabilities describing the recombination map. The map is completely up to
71  the programmer, and it is not checked for sanity at all in this function. For a region of k sites,
72  indexes 0 to k-2 of this array should be filled. The i-th element should contain the probability
73  that a crossover occurs between position i and i+1. The sum of all elements should be 1, such
74  that the array describes the recombination map in terms of a probability distribution function. An
75  example of how to do this is in the file examples/msbeta.cc that comes with the source for this library.
76  \return a pair of integers containing the index of the recombinant chromosome (.first),
77  and the position at which the crossover will occur (.second)
78  \ingroup coalescent
79  */
80  {
81  return pick_spot_details(uni01,total_reclen,reclens,sample_begin,current_nsam,rec_map);
82  }
83 
84  template<typename uniform01_generator>
85  std::pair<int,int> pick_spot( const uniform01_generator & uni01,
86  const double & total_reclen,
87  const std::vector<double> & reclens,
88  std::vector<chromosome>::const_iterator sample_begin,
89  const unsigned & current_nsam,
90  const double * rec_map)
91  /*!
92  Picks a positions amongst all chromosomes at which a recombination event
93  will occur, based on an arbitrary genetic map
94  \param uni01 a function/object which takes no arguments and can return a U[0,1)
95  \param total_reclen the total recombination length of all chromosomes in the sample
96  \param reclens a vector of the proportion of \a total_reclen contributed by each chromosome.
97  This needs to be ordered in the same order as \a sample_begin to (\a sample_begin + \a current_nsam - 1)
98  \param sample_begin an iterator pointing to the beginning of the sample
99  \param current_nsam the current sample size in the simulation
100  \param rec_map an array of probabilities describing the recombination map. The map is completely up to
101  the programmer, and it is not checked for sanity at all in this function. For a region of k sites,
102  indexes 0 to k-2 of this array should be filled. The i-th element should contain the probability
103  that a crossover occurs between position i and i+1. The sum of all elements should be 1, such
104  that the array describes the recombination map in terms of a probability distribution function. An
105  example of how to do this is in the file examples/msbeta.cc that comes with the source for this library.
106  \return a pair of integers containing the index of the recombinant chromosome (.first),
107  and the position at which the crossover will occur (.second)
108  \ingroup coalescent
109  */
110  {
111  return pick_spot_details(uni01,total_reclen,reclens,sample_begin,current_nsam,rec_map);
112  }
113  } //ns coalsim
114 }//namespace Sequence
115 #endif