Merge pull request #1 from FaridKaveh/mutation_model

Mutation model

MoranModel/MoranModel.cpp

@@ -32,13 +32,12 @@ void MoranProcess::generateTree(){
 
 void MoranProcess::generateMuts(){
 
-    mutations.clear();
+    mutations.assign(population, std::vector<double> ());
 
     std::random_device rd;
     std::default_random_engine engine(rd());
-
-    double total_time = std::accumulate(event_times.begin(), event_times.end(), 0.0); 
-    double rate = population* total_time * THETA/2; 
+
+    double rate = events * THETA/(2*population); 
 
     std::poisson_distribution<int> mut_dist(rate);
     mut_number = mut_dist(engine); 
@@ -50,14 +49,15 @@ void MoranProcess::generateMuts(){
         weights.at(i) = i > 0 ? weights.at(i-1)+event_times.at(i) : event_times.at(i); 
     }  
 
-    std::uniform_real_distribution<> mut_drop(0, total_time); 
+    std::uniform_int_distribution<> mut_drop(0, events-1); 
     double draw;
     int box;
     std::vector<int> allocations (events, 0);
 
     for (int i = 0; i < mut_number; i++){ 
-        draw = mut_drop(engine); 
-        box = binarySearch(draw, weights);
+
+        box = mut_drop(engine); 
+
         try{ 
         ++allocations.at(box);
         }
@@ -67,17 +67,24 @@ void MoranProcess::generateMuts(){
         }
     }
 
-    mutations.insert(mutations.end(), population*events, 0); 
+
     std::uniform_int_distribution<> pick_line (0,population-1); 
+    std::uniform_real_distribution<> mutant_id (0,1); 
     int line; 
+    unsigned id = 0; 
 
-    for (int i = 0; i < events; ++i) { 
-        while (allocations.at(i) > 0){
-        line = pick_line(engine);
-        ++mutations.at(population*i+line);
+    for (unsigned i = 0; i < events; ++i) { 
+
+        line = event_history.at(2*i + 1);
 
+        while (allocations.at(i) > 0){
+
+        mutations.at(line).push_back(++id);
         --allocations.at(i);
+
         }
+
+        mutations.at(event_history.at(2*i)) = mutations.at(event_history.at(2*i+1));
     }
 }
 
@@ -121,16 +128,17 @@ int MoranProcess::calculateFamilyHistories(bool draw){
     return i;
 }
 
-int MoranProcess::calculateNumberOfMutationEvents(){ 
-    return std::accumulate(mutations.begin(), mutations.end(), 0); 
-}
 
 std::vector<int> MoranProcess::calcualteSegregatingSites(){ 
 
     std::vector<int> segregating_sites (population, 0);
     return segregating_sites;
 }
 
+std::vector<int> MoranProcess::calculateSiteFrequencySpectrum(){ 
+    return std::vector<int> ();
+}
+
 std::vector < std::vector<int> > MoranProcess::buildCoalescentTree(int level){ 
     if (level == -1){ 
         level = events; 

MoranModel/MoranModel.h

@@ -17,7 +17,7 @@ stores a vector of arrays that records the history. The statistics can all be in
 
     std::vector<double> event_times;
     std::vector <int> event_history; 
-    std::vector <int> mutations;
+    std::vector < std::vector<double> > mutations;
 
     void generatePath(){
         generateTree();
@@ -41,13 +41,16 @@ stores a vector of arrays that records the history. The statistics can all be in
 
     std::vector<double> getEventTimes(){return this -> event_times;}
     std::vector<int> getEventHistory(){return this -> event_history;}
-    std::vector<int> getMutations(){return this -> mutations;}
+    std::vector< std::vector<double> > getMutations(){return this -> mutations;}
 
     int calculateFamilyHistories(bool draw = false);
-    int calculateNumberOfMutationEvents();
     std::vector<int> calcualteSegregatingSites();
+    std::vector<int> calculateSiteFrequencySpectrum(); 
+
+
     std::vector< std::vector<int> > buildCoalescentTree(int level = -1);
 
+
 };
 
 #endif

main.cpp

@@ -45,7 +45,9 @@ int main(int argc, char* argv[]){
 
     int average_mutations = average_quotient + average_remainder / runs; 
     std::cout << "<S> = " << average_mutations << std::endl; 
-    // printVector(moran.getMutations());
+
+    printVectorOfVectors(moran.getMutations());
+    // printVectorOfVectors(moran.getMutations());
 
     return 0;
 }