Merge pull request #8 from smithlabcode/RemoveStateSeq

Remove state seq
smithlabcode · Feb 25, 2020 · ba639df · ba639df
2 parents 3068aa0 + 00214ed
commit ba639df
Show file tree

Hide file tree

Showing 28 changed files with 471 additions and 752 deletions.
diff --git a/src/common/EndCondSampling.cpp b/src/common/EndCondSampling.cpp
diff --git a/src/common/EpiEvoModel.cpp b/src/common/EpiEvoModel.cpp
@@ -20,7 +20,6 @@
  */
 
 #include "EpiEvoModel.hpp"
-#include "StateSeq.hpp"
 #include "epievo_utils.hpp"
 
 #include <iostream>

diff --git a/src/common/EpiEvoModel.hpp b/src/common/EpiEvoModel.hpp
@@ -25,6 +25,7 @@
 #include <string>
 #include <vector>
 #include <random>
+
 #include "epievo_utils.hpp"
 
 struct EpiEvoModel {
@@ -42,7 +43,7 @@ struct EpiEvoModel {
   two_by_two init_T;  // horizontal transition probs (initial)
   two_by_two Q;       // pair-wise potential densities (stationary)
   bool use_init_T;    // whether to use initial prior
-  
+
   std::vector<double> triplet_rates; // rates for triples
 
   /*********************************************************************

diff --git a/src/common/GlobalJump.cpp b/src/common/GlobalJump.cpp
@@ -21,8 +21,6 @@
 
 #include "GlobalJump.hpp"
 
-#include "StateSeq.hpp"
-
 #include <iostream>
 #include <fstream>
 #include <string>
@@ -72,15 +70,15 @@ operator<<(std::ostream &os, const GlobalJump &s) {
 
 void
 write_root_to_pathfile_global(const string &pathfile, const string &root_name,
-                              const StateSeq &root) {
+                              const state_seq &root) {
   std::ofstream of;
-  if (!pathfile.empty()) of.open(pathfile.c_str());
+  if (!pathfile.empty()) of.open(pathfile);
   std::ostream outpath(pathfile.empty() ? std::cout.rdbuf() : of.rdbuf());
   if (!outpath)
     throw std::runtime_error("bad output file: " + pathfile);
 
   outpath << ROOT_TAG << ':' << root_name << '\n';
-  copy(root.seq.begin(), root.seq.end(),
+  copy(begin(root), end(root),
        std::ostream_iterator<bool>(outpath));
   outpath << '\n';
 }
@@ -89,24 +87,24 @@ void
 append_to_pathfile_global(const string &pathfile, const string &node_name,
                           const vector<GlobalJump> &the_path) {
   std::ofstream of;
-  if (!pathfile.empty()) of.open(pathfile.c_str(), std::ofstream::app);
+  if (!pathfile.empty()) of.open(pathfile, std::ofstream::app);
   std::ostream outpath(pathfile.empty() ? std::cout.rdbuf() : of.rdbuf());
   if (!outpath)
     throw std::runtime_error("bad output file: " + pathfile);
-  
+
   outpath << NODE_TAG << ':' << node_name << '\n';
   for (size_t i = 0; i < the_path.size(); ++i)
     outpath << the_path[i] << '\n';
 }
 
 void
-read_pathfile_global(const string &pathfile, StateSeq &root,
+read_pathfile_global(const string &pathfile, state_seq &root,
                      vector<string> &node_names,
                      vector<vector<GlobalJump> > &the_paths) {
 
   the_paths.clear();
 
-  std::ifstream in(pathfile.c_str());
+  std::ifstream in(pathfile);
   if (!in)
     throw std::runtime_error("cannot read: " + pathfile);
 
@@ -120,7 +118,7 @@ read_pathfile_global(const string &pathfile, StateSeq &root,
 
   getline(in, buffer);
   for (size_t i = 0; i < buffer.length(); ++i)
-    root.seq.push_back(buffer[i] == '1');
+    root.push_back(buffer[i] == '1');
 
   the_paths.push_back(vector<GlobalJump>()); // empty on purpose
 

diff --git a/src/common/GlobalJump.hpp b/src/common/GlobalJump.hpp
@@ -23,7 +23,7 @@
 #include <string>
 #include <vector>
 
-struct StateSeq;
+#include "epievo_utils.hpp"
 
 struct GlobalJump {
   /* GlobalJump: represents a change in the boolean state within a
@@ -32,7 +32,8 @@ struct GlobalJump {
      deduced from some ground state (i.e. the root sequence).
    */
   GlobalJump() {}
-  GlobalJump(const double tp, const size_t pos) : timepoint(tp), position(pos) {}
+  GlobalJump(const double tp, const size_t pos) :
+    timepoint(tp), position(pos) {}
   double timepoint; // time of change
   size_t position; // position of change
   bool operator<(const GlobalJump &other) const {
@@ -49,14 +50,14 @@ operator>>(std::istream &is, GlobalJump &s);
 void
 write_root_to_pathfile_global(const std::string &pathfile,
                               const std::string &root_name,
-                              const StateSeq &root);
+                              const std::vector<bool> &root);
 
 void
 append_to_pathfile_global(const std::string &pathfile,
                           const std::string &node_name,
                           const std::vector<GlobalJump> &the_path);
 
 void
-read_pathfile_global(const std::string &pathfile, StateSeq &root,
+read_pathfile_global(const std::string &pathfile, std::vector<bool> &root,
                      std::vector<std::string> &node_names,
                      std::vector<std::vector<GlobalJump> > &the_paths);
diff --git a/src/common/IntervalSampler.cpp b/src/common/IntervalSampler.cpp
@@ -19,6 +19,8 @@
  * 02110-1301 USA
  */
 
+#include "IntervalSampler.hpp"
+
 #include <string>
 #include <vector>
 #include <cassert>
@@ -29,10 +31,10 @@
 #include <functional>
 #include <iomanip>
 
-#include "IntervalSampler.hpp"
+#include "EndCondSampling.hpp"
+#include "Segment.hpp"
 #include "PhyloTreePreorder.hpp"
 #include "Path.hpp"
-#include "StateSeq.hpp"
 
 using std::vector;
 using std::endl;

diff --git a/src/common/IntervalSampler.hpp b/src/common/IntervalSampler.hpp
@@ -24,8 +24,6 @@
 #include "Path.hpp"
 #include "EpiEvoModel.hpp"
 #include "TreeHelper.hpp"
-#include "EndCondSampling.hpp"
-#include "Segment.hpp"
 
 bool
 Metropolis_Hastings_interval(const EpiEvoModel &the_model, const TreeHelper &th,
@@ -36,7 +34,7 @@ Metropolis_Hastings_interval(const EpiEvoModel &the_model, const TreeHelper &th,
 
 /* CASE 1: Path ending in non-bifurcating nodes
   - end-conditioned uniformization sampling
- 
+
    CASE 2: Path involving bifurcating nodes
   - calculate posterior distribution at bifurcating node
   - end-conditioned uniformization sampling along involving intervals

diff --git a/src/common/ParamEstimation.cpp b/src/common/ParamEstimation.cpp
@@ -33,11 +33,12 @@
 
 #include "PhyloTreePreorder.hpp"
 #include "Path.hpp"
-#include "StateSeq.hpp"
 #include "EpiEvoModel.hpp"
 #include "TreeHelper.hpp"
 #include "TripletSampler.hpp" /* forward simulation */
 
+#include "epievo_utils.hpp"
+
 using std::vector;
 using std::endl;
 using std::cerr;
@@ -104,17 +105,14 @@ get_sufficient_statistics(const vector<vector<Path> > &all_paths,
 }
 
 void
-get_root_frequences(const vector<vector<Path> > &all_paths,
-                    vector<vector<double> > &counts) {
-
+get_root_frequencies(const vector<vector<Path>> &all_paths, two_by_two &counts) {
   assert(all_paths.size() >= 2 && !all_paths[1].empty());
 
-  counts = vector<vector<double> >(2, vector<double>(2, 0.0));
-
+  counts.reset();
   size_t prev = all_paths[1].front().init_state;
   for (size_t i = 1; i < all_paths[1].size(); ++i) {
     const size_t curr = all_paths[1][i].init_state;
-    counts[prev][curr]++;
+    counts(prev, curr)++;
     prev = curr;
   }
 }
@@ -398,7 +396,7 @@ compute_estimates_rates_and_branches(const bool VERBOSE,
   candidate_branches(J, D, updated_rates, branch_scale);
   std::transform(branch_scale.begin(), branch_scale.end(),
                  th.branches.begin(), updated_branches.begin(),
-                 std::multiplies<double>() );
+                 std::multiplies<double>());
 
   set_one_change_per_site_per_unit_time(updated_rates, updated_branches);
   the_model.rebuild_from_triplet_rates(updated_rates);
@@ -414,7 +412,7 @@ compute_estimates_rates_and_branches(const bool VERBOSE,
   for (size_t b = 1; b < th.n_nodes; b++)
     for (size_t i = 0; i < D[b].size(); i++) {
       J_collapsed[i] += J[b][i];
-      D_collapsed[i] += branch_scale[b] * D[b][i];
+      D_collapsed[i] += branch_scale[b]*D[b][i];
     }
 
   return log_likelihood(J_collapsed, D_collapsed, updated_rates);
@@ -448,14 +446,11 @@ compute_estimates_rates_and_branches(const bool VERBOSE,
 
 
 void
-estimate_root_distribution(const vector<vector<double> > &counts,
-                           EpiEvoModel &the_model) {
-
-  the_model.init_T[0][0] = counts[0][0] / (counts[0][0] + counts[0][1]);
-  the_model.init_T[0][1] = 1 - the_model.init_T[0][0];
-
-  the_model.init_T[1][1] = counts[1][1] / (counts[1][1] + counts[1][0]);
-  the_model.init_T[1][0] = 1 - the_model.init_T[1][1];
+estimate_root_distribution(const two_by_two &counts, EpiEvoModel &the_model) {
+  the_model.init_T[0][0] = counts(0, 0)/(counts(0, 0) + counts(0, 1));
+  the_model.init_T[0][1] = 1.0 - the_model.init_T[0][0];
+  the_model.init_T[1][1] = counts(1, 1)/(counts(1, 1) + counts(1, 0));
+  the_model.init_T[1][0] = 1.0 - the_model.init_T[1][1];
 }
 
 void
@@ -464,7 +459,7 @@ estimate_root_distribution(const vector<vector<Path> > &all_paths,
 
   assert(all_paths.size() >= 2 && !all_paths[1].empty());
 
-  vector<vector<double> > counts;
-  get_root_frequences(all_paths, counts);
+  two_by_two counts;
+  get_root_frequencies(all_paths, counts);
   estimate_root_distribution(counts, the_model);
 }
diff --git a/src/common/ParamEstimation.hpp b/src/common/ParamEstimation.hpp
@@ -27,6 +27,8 @@
 #include "TreeHelper.hpp"
 #include "EpiEvoModel.hpp"
 
+#include "epievo_utils.hpp"
+
 #include <vector>
 
 double
@@ -59,8 +61,7 @@ compute_estimates_rates_and_branches(const bool VERBOSE,
                                      EpiEvoModel &the_model);
 
 void
-estimate_root_distribution(const std::vector<std::vector<double> > &counts,
-                           EpiEvoModel &the_model);
+estimate_root_distribution(const two_by_two &counts, EpiEvoModel &the_model);
 
 void
 estimate_root_distribution(const std::vector<std::vector<Path> > &all_paths,
@@ -78,10 +79,10 @@ get_sufficient_statistics(const std::vector<std::vector<Path> > &all_paths,
                           std::vector<std::vector<double> > &J,
                           std::vector<std::vector<double> > &D);
 
-// count duplet frequences from root sequence
+// count duplet frequencies from root sequence
 void
-get_root_frequences(const std::vector<std::vector<Path> > &all_paths,
-                    std::vector<std::vector<double> > &counts);
+get_root_frequencies(const std::vector<std::vector<Path> > &all_paths,
+                     two_by_two &counts);
 
 // scale jump times according to updated branch lengths
 void

diff --git a/src/common/Path.cpp b/src/common/Path.cpp
@@ -32,7 +32,7 @@
 
 #include "smithlab_utils.hpp"
 
-#include "StateSeq.hpp"
+#include "epievo_utils.hpp"
 
 using std::vector;
 using std::string;
@@ -65,8 +65,8 @@ operator<<(std::ostream &os, const Path &p) {
 }
 
 void
-initialize_paths(const std::vector<bool> &seq, const double tot_time,
-                 std::vector<Path> &paths) {
+initialize_paths(const vector<bool> &seq, const double tot_time,
+                 vector<Path> &paths) {
   paths.resize(seq.size());
   for (size_t i = 0; i < seq.size(); ++i)
     paths[i] = Path(seq[i], tot_time);
@@ -177,12 +177,12 @@ void get_seq_init(const vector<Path> &paths, vector<bool> &seq) {
 ////////////////////////////////////////////////////////////////////////////////
 
 struct TriplePath {
-  std::vector<size_t> states; // triplet states, length k
-  std::vector<double> breaks; // start is first jump, end is total_time, length k
-  std::vector<size_t> jump_context_freq; // context freq. of jumps at middle site
+  vector<size_t> states; // triplet states, length k
+  vector<double> breaks; // start is first jump, end is total_time, length k
+  vector<size_t> jump_context_freq; // context freq. of jumps at middle site
 
   TriplePath(const Path &l, const Path &m, const Path &r);
-  void time_by_context(std::vector<double> &tbc) const;
+  void time_by_context(vector<double> &tbc) const;
 };
 
 TriplePath::TriplePath(const Path &l, const Path &m, const Path &r) {

diff --git a/src/common/Path.hpp b/src/common/Path.hpp
@@ -28,8 +28,6 @@
 #include <sstream>
 #include <cassert>
 
-#include "smithlab_utils.hpp"
-
 struct Path {
 
   Path() : init_state(false), tot_time(0.0) {}

diff --git a/src/common/Segment.cpp b/src/common/Segment.cpp
@@ -25,7 +25,7 @@
 #include <exception>
 
 #include "Segment.hpp"
-#include "StateSeq.hpp"
+#include "epievo_utils.hpp"
 
 using std::vector;
 

diff --git a/src/common/Segment.hpp b/src/common/Segment.hpp
@@ -37,10 +37,8 @@ struct SegmentInfo {
   double len;
 };
 
-
 void
 collect_segment_info(const std::vector<double> &rates,
                      const Path &l, const Path &r,
                      std::vector<SegmentInfo> &seg_info);
-
 #endif
diff --git a/src/common/SingleSiteSampler.cpp b/src/common/SingleSiteSampler.cpp
@@ -33,7 +33,9 @@
 #include "PhyloTreePreorder.hpp"
 #include "Path.hpp"
 #include "ParamEstimation.hpp"
-#include "StateSeq.hpp"
+#include "EndCondSampling.hpp"
+
+#include "epievo_utils.hpp"
 
 using std::vector;
 using std::endl;

diff --git a/src/common/SingleSiteSampler.hpp b/src/common/SingleSiteSampler.hpp
@@ -25,7 +25,6 @@
 #include "Segment.hpp"
 #include "EpiEvoModel.hpp"
 #include "TreeHelper.hpp"
-#include "EndCondSampling.hpp"
 
 struct FelsHelper {
   FelsHelper() {}
@@ -41,7 +40,7 @@ pruning(const TreeHelper &th, const size_t site_id,
         const std::vector<std::vector<double> > &emit,
         const std::vector<std::vector<SegmentInfo> > &seg_info,
         std::vector<FelsHelper> &fh);
- 
+
 void
 downward_sampling(const EpiEvoModel &the_model, const TreeHelper &th,
                   const size_t site_id,