2 #ifndef __SEQUENCE_COALESCENT_BITS_RECOMBINATION_TCC__
3 #define __SEQUENCE_COALESCENT_BITS_RECOMBINATION_TCC__
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)
24 while( (sample_begin+recombinant) <
25 (sample_begin+current_nsam) && ! flag)
27 sum += reclens[std::vector<double>::size_type(recombinant)];
28 if ( ran <= sum/total_reclen )
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)
39 if ( ran2 <= sum2/cum_prob )
50 return std::make_pair(recombinant,pos);
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)
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)
81 return pick_spot_details(uni01,total_reclen,reclens,sample_begin,current_nsam,rec_map);
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)
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)
111 return pick_spot_details(uni01,total_reclen,reclens,sample_begin,current_nsam,rec_map);
114 }//namespace Sequence