style: rename L to neighborhood_size (#22)

* Mainly this is to replace the two Ls in Matcher.hpp and ImputationMatcher.hpp, however I've applied to all other instances too.

Author: alexallmont

src/HMM.cpp

@@ -75,9 +75,9 @@ void HMM::compute_mutation_scores(std::vector<double> bp_sizes, double mutation_
 std::vector<int> HMM::breakpoints(std::vector<bool> observations, int start) {
   // Viterbi
   // Initialize
-  int L = observations.size();
-  std::vector<unsigned short> z(L);
-  int end = start + L;
+  int neighborhood_size = observations.size();
+  std::vector<unsigned short> z(neighborhood_size);
+  int end = start + neighborhood_size;
   for (int i = 0; i < num_states; i++) {
     double score = observations[0] ? het_score[start][i] : hom_score[start][i];
     trellis[start][i] = score;
@@ -86,7 +86,7 @@ std::vector<int> HMM::breakpoints(std::vector<bool> observations, int start) {
   // Main routine
   double score;
   unsigned short running_argmax;
-  for (int j = 1; j < L; j++) {
+  for (int j = 1; j < neighborhood_size; j++) {
     for (int i = 0; i < num_states; i++) {
       double running_max = 0;
       for (int k = 0; k < num_states; k++) {
@@ -107,26 +107,26 @@ std::vector<int> HMM::breakpoints(std::vector<bool> observations, int start) {
   }
 
   // Get best path
-  double running_max = trellis[start + L - 1][0];
+  double running_max = trellis[start + neighborhood_size - 1][0];
   unsigned short argmax = 0;
   for (int k = 1; k < num_states; k++) {
-    double s = trellis[start + L - 1][k];
+    double s = trellis[start + neighborhood_size - 1][k];
     if (s > running_max) {
       running_max = s;
       argmax = k;
     }
   }
 
   // Traceback
-  z[L - 1] = argmax;
-  for (int j = L - 1; j >= 1; j--) {
+  z[neighborhood_size - 1] = argmax;
+  for (int j = neighborhood_size - 1; j >= 1; j--) {
     z[j - 1] = pointers[j + start][z[j]];
   }
 
   // Break it up
   std::vector<int> breakpoints;
   breakpoints.push_back(0 + start);
-  for (int j = 1; j < L; j++) {
+  for (int j = 1; j < neighborhood_size; j++) {
     if (z[j - 1] != z[j]) {
       breakpoints.push_back(j + start);
     }

src/ImputationMatcher.cpp

@@ -11,9 +11,9 @@
 
 ImputationMatcher::ImputationMatcher(int _n_ref, int _n_target,
                                      const std::vector<double>& _genetic_positions,
-                                     double _query_interval_size, int _L)
+                                     double _query_interval_size, int _neighborhood_size)
     : num_reference(_n_ref), num_target(_n_target), genetic_positions(_genetic_positions),
-      query_interval_size(_query_interval_size), L(_L) {
+      query_interval_size(_query_interval_size), neighborhood_size(_neighborhood_size) {
   if (genetic_positions.size() <= 2) {
     throw std::runtime_error("Need at least 3 sites, found " +
                              std::to_string(genetic_positions.size()));
@@ -132,7 +132,7 @@ void ImputationMatcher::process_site(const std::vector<int>& genotype) {
     }
 
     for (auto& sorting_twople : target_sort) {
-      // get L-sized neighbourhood per target sample around target_sort[target_id] in ref_sorting
+      // get neighborhood-sized per target sample around target_sort[target_id] in ref_sorting
       int target_id = sorting_twople.first;
       int sorting_idx = sorting_twople.second;
       int insert_start;

src/ImputationMatcher.hpp

@@ -24,18 +24,17 @@ class ImputationMatcher {
   // adapted to be used for imputation.
   // TODO: include a second pass through data here to get divergence values and not do that using
   // Threads-fastLS
-  int L = 0;
+  int neighborhood_size = 0;
   std::vector<int> query_sites;
   int num_samples = 0;
   int num_reference = 0;
   int num_target = 0;
   int num_sites = 0;
   double query_interval_size = 0.0;
   std::unordered_map<int, std::unordered_set<int>> match_sets;
-  int neighborhood_size = 0;
   std::vector<double> genetic_positions;
   ImputationMatcher(int _n_ref, int _n_target, const std::vector<double>& _genetic_positions,
-                    double _query_interval_size, int _L);
+                    double _query_interval_size, int _neighborhood_size);
 
   // Do all the work
   void process_site(const std::vector<int>& genotype);

src/Matcher.cpp

@@ -75,7 +75,7 @@ void MatchGroup::filter_matches(int min_matches) {
     }
   }
 
-  // Then determine top 4 candidates for neighbouring groups
+  // Then determine top 4 candidates for neighboring groups
   top_four_maps.reserve(num_samples);
   for (int i = 0; i < num_samples; i++) {
     top_four_maps.emplace_back(std::min(4, (int) match_candidates.at(i).size()));
@@ -104,10 +104,10 @@ void MatchGroup::clear() {
 }
 
 Matcher::Matcher(int _n, const std::vector<double>& _genetic_positions, double _query_interval_size,
-                 double _match_group_interval_size, int _L, int _min_matches)
+                 double _match_group_interval_size, int _neighborhood_size, int _min_matches)
     : num_samples(_n), genetic_positions(_genetic_positions),
       query_interval_size(_query_interval_size),
-      match_group_interval_size(_match_group_interval_size), L(_L), min_matches(_min_matches) {
+      match_group_interval_size(_match_group_interval_size), neighborhood_size(_neighborhood_size), min_matches(_min_matches) {
   if (genetic_positions.size() <= 2) {
     throw std::runtime_error("Need at least 3 sites, found " +
                              std::to_string(genetic_positions.size()));
@@ -216,7 +216,7 @@ void Matcher::process_site(const std::vector<int>& genotype) {
   }
   sorting = next_sorting;
 
-  // Threading-neighbour queries
+  // Threading-neighbor queries
   if (match_group_idx < match_group_sites.size() - 1 &&
       sites_processed >= match_group_sites.at(match_group_idx + 1)) {
     match_group_idx++;
@@ -239,17 +239,17 @@ void Matcher::process_site(const std::vector<int>& genotype) {
     for (int i = 1; i < num_samples; i++) {
       std::vector<int> matches;
       int allele = genotype.at(i);
-      matches.reserve(L);
+      matches.reserve(neighborhood_size);
       auto iter = threaded.insert(permutation.at(i));
       auto iter_up = iter.first;
       auto iter_down = iter.first;
       // Check if genotypes are identical, just to be sure
-      while (matches.size() < L && (iter_down != threaded.begin() || iter_up != threaded.end())) {
+      while (matches.size() < neighborhood_size && (iter_down != threaded.begin() || iter_up != threaded.end())) {
         if (iter_down != threaded.begin()) {
           iter_down--;
           matches.push_back(sorting.at(*iter_down));
         }
-        if (matches.size() < L && iter_up != threaded.end()) {
+        if (matches.size() < neighborhood_size && iter_up != threaded.end()) {
           iter_up++;
           if (iter_up != threaded.end()) {
             matches.push_back(sorting.at(*iter_up));

src/Matcher.hpp

@@ -35,18 +35,17 @@ class Matcher {
   std::vector<int> permutation;
 
 public:
-  int L = 0;
+  int neighborhood_size = 0;
   std::vector<int> query_sites;
   std::vector<int> match_group_sites;
   int num_samples = 0;
   int num_sites = 0;
   double query_interval_size = 0.0;
   // matches in these groups are considered together in the hmm
   double match_group_interval_size = 0.0;
-  int neighborhood_size = 0;
   std::vector<double> genetic_positions;
   Matcher(int _n, const std::vector<double>& _genetic_positions, double _query_interval_size,
-          double _match_group_interval_size, int _L, int _min_matches);
+          double _match_group_interval_size, int _neighborhood_size, int _min_matches);
 
   // Do all the work
   void process_site(const std::vector<int>& genotype);

src/ThreadsFastLS.cpp

@@ -1017,12 +1017,12 @@ std::vector<std::tuple<int, std::vector<int>>> ThreadsFastLS::traceback(Tracebac
 }
 
 /**
- * Similar to normal traceback, but picks the up-to L best matches and stores their overlap with the
+ * Similar to normal traceback, but picks the up-to neighborhood_size best matches and stores their overlap with the
  * input sequence returns a list of a tuple of lists :-P I.e., a list of segments, and each segment
- * is a tuple (sample_IDs, starts, ends) all of equal length <= L
+ * is a tuple (sample_IDs, starts, ends) all of equal length <= neighborhood_size
  */
 std::vector<std::tuple<std::vector<int>, std::vector<int>, std::vector<int>>>
-ThreadsFastLS::traceback_impute(std::vector<bool>& genotypes, TracebackState* tb, Node* match, int L) {
+ThreadsFastLS::traceback_impute(std::vector<bool>& genotypes, TracebackState* tb, Node* match, int neighborhood_size) {
   std::vector<std::tuple<std::vector<int>, std::vector<int>, std::vector<int>>> imputation_path;
   int prev_end = num_sites;
   while (tb != nullptr) {
@@ -1084,7 +1084,7 @@ ThreadsFastLS::traceback_impute(std::vector<bool>& genotypes, TracebackState* tb
     std::vector<int> segment_starts;
     std::vector<int> sample_ids;
     std::vector<int> segment_ends;
-    for (int j = idx.size() - 1; j >= std::max(0, (int) (idx.size() - L)); j--) {
+    for (int j = idx.size() - 1; j >= std::max(0, (int) (idx.size() - neighborhood_size)); j--) {
       segment_starts.push_back(segment_start);
       segment_ends.push_back(segment_end);
       sample_ids.push_back(div_states[idx[j]]);
@@ -1344,7 +1344,6 @@ ThreadsFastLS::threads_ls(const std::vector<bool>& genotype) {
   return best_path;
 }
 
-// todo: L is no longer required here
 std::tuple<std::vector<int>, std::vector<std::vector<int>>, std::vector<double>, std::vector<int>>
 ThreadsFastLS::thread(const int new_sample_ID, const std::vector<bool>& genotype) {
   // Compute LS path
@@ -1415,13 +1414,13 @@ ThreadsFastLS::thread(const int new_sample_ID, const std::vector<bool>& genotype
   return remove_burn_in(bp_starts, target_IDs, segment_ages, het_sites);
 }
 
-std::vector<ImputationSegment> ThreadsFastLS::impute(std::vector<bool>& genotype, int L) {
+std::vector<ImputationSegment> ThreadsFastLS::impute(std::vector<bool>& genotype, int neighborhood_size) {
   // vector of sample_ids, seg_starts, seg_ends (buffered)
   std::vector<std::tuple<std::vector<int>, std::vector<int>, std::vector<int>>> best_path;
   Node* match;
   TracebackState* traceback;
   std::tie(traceback, match) = fastLS(genotype);
-  best_path = traceback_impute(genotype, traceback, match, L);
+  best_path = traceback_impute(genotype, traceback, match, neighborhood_size);
   traceback_states.clear();
 
   std::vector<double> seg_ages;

src/ThreadsFastLS.hpp

@@ -111,13 +111,13 @@ class ThreadsFastLS {
                  std::vector<double>& segment_ages, std::vector<int>& het_sites);
   std::vector<std::tuple<int, std::vector<int>>> threads_ls(const std::vector<bool>& genotype);
 
-  std::vector<ImputationSegment> impute(std::vector<bool>& genotype, int L);
+  std::vector<ImputationSegment> impute(std::vector<bool>& genotype, int neighborhood_size);
   std::pair<TracebackState*, Node*> fastLS(const std::vector<bool>& genotype,
                                            bool imputation = false);
   std::vector<std::tuple<int, std::vector<int>>> traceback(TracebackState* tb, Node* match,
                                                            bool return_all = false);
   std::vector<std::tuple<std::vector<int>, std::vector<int>, std::vector<int>>>
-  traceback_impute(std::vector<bool>& genotypes, TracebackState* tb, Node* match, int L);
+  traceback_impute(std::vector<bool>& genotypes, TracebackState* tb, Node* match, int neighborhood_size);
   std::array<std::pair<TracebackState*, Node*>, 2> fastLS_diploid(const std::vector<int>& genotype);
   std::array<std::vector<std::tuple<int, std::vector<int>>>, 2>
   diploid_ls(std::vector<int> unphased_genotypes);

src/threads_arg/__main__.py

@@ -394,13 +394,13 @@ def threads_infer(pgen, map_gz, recombination_rate, demography, mutation_rate, q
     # There are four params here to mess with:
     #   - query interval
     #   - match group size
-    #   - neighbourhood size (L)
+    #   - neighborhood size
     #   - min_matches
     # Keep min_matches low for small sample sizes, can be 2 up to ~1,000 but then > 3
     # Smaller numbers run faster and consume less memory
     MIN_MATCHES = 4
-    L = 4
-    matcher = Matcher(2 * num_samples, genetic_positions[ac_mask], query_interval, match_group_interval, L, MIN_MATCHES)
+    neighborhood_size = 4
+    matcher = Matcher(2 * num_samples, genetic_positions[ac_mask], query_interval, match_group_interval, neighborhood_size, MIN_MATCHES)
 
     # Matching step
     for b in range(n_batches):

src/threads_arg_pybind.cpp

@@ -63,7 +63,7 @@ PYBIND11_MODULE(threads_arg_python_bindings, m) {
   py::class_<Matcher>(m, "Matcher")
       .def(py::init<int, std::vector<double>, double, double, int, int>(), "Initialize",
            py::arg("num_samples"), py::arg("genetic_positions"), py::arg("query_interval_size"),
-           py::arg("match_group_interval_size"), py::arg("L") = 4, py::arg("min_matches") = 4)
+           py::arg("match_group_interval_size"), py::arg("neighborhood_size") = 4, py::arg("min_matches") = 4)
       .def_readonly("query_sites", &Matcher::query_sites)
       .def_readonly("genetic_positions", &Matcher::genetic_positions)
       .def_readonly("query_interval_size", &Matcher::query_interval_size)
@@ -81,7 +81,7 @@ PYBIND11_MODULE(threads_arg_python_bindings, m) {
   py::class_<ImputationMatcher>(m, "ImputationMatcher")
       .def(py::init<int, int, const std::vector<double>&, double, int>(), "Initialize",
            py::arg("num_reference"), py::arg("num_target"), py::arg("genetic_positions"),
-           py::arg("query_interval_size"), py::arg("L") = 8)
+           py::arg("query_interval_size"), py::arg("neighborhood_size") = 8)
       .def("process_site", &ImputationMatcher::process_site)
       .def("get_matches", &ImputationMatcher::get_matches);