Optimize multi-tree comparisons

R/RcppExports.R

@@ -57,6 +57,10 @@ cpp_msi_distance <- function(x, y, nTip) {
     .Call(`_TreeDist_cpp_msi_distance`, x, y, nTip)
 }
 
+cpp_all_pairs_mci <- function(x, nTip) {
+    .Call(`_TreeDist_cpp_all_pairs_mci`, x, nTip)
+}
+
 cpp_mutual_clustering <- function(x, y, nTip) {
     .Call(`_TreeDist_cpp_mutual_clustering`, x, y, nTip)
 }

R/tree_distance_info.R

@@ -297,8 +297,15 @@ ExpectedVariation <- function (tree1, tree2, samples = 1e+4) {
 #' @export
 MutualClusteringInfo <- function (tree1, tree2 = NULL, normalize = FALSE,
                                   reportMatching = FALSE, diag = TRUE) {
-  unnormalized <- CalculateTreeDistance(MutualClusteringInfoSplits, tree1, tree2,
-                                        reportMatching)
+  unnormalized <- if (is.null(tree2) && is.list(tree1) &&
+                      !inherits(tree1, 'phylo')) {
+    structure(cpp_all_pairs_mci(as.Splits(tree1, tree1[[1]]), NTip(tree1[[1]])),
+              class = "dist",
+              Size = length(tree1), diag = FALSE, upper = FALSE)
+  } else {
+    CalculateTreeDistance(MutualClusteringInfoSplits, tree1, tree2,
+                          reportMatching)
+  }
   if (diag && is.null(tree2)) {
     unnormalized <- as.matrix(unnormalized)
     diag(unnormalized) <- ClusteringEntropy(tree1)

src/RcppExports.cpp

@@ -177,6 +177,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// cpp_all_pairs_mci
+NumericVector cpp_all_pairs_mci(const List x, const IntegerVector nTip);
+RcppExport SEXP _TreeDist_cpp_all_pairs_mci(SEXP xSEXP, SEXP nTipSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const List >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const IntegerVector >::type nTip(nTipSEXP);
+    rcpp_result_gen = Rcpp::wrap(cpp_all_pairs_mci(x, nTip));
+    return rcpp_result_gen;
+END_RCPP
+}
 // cpp_mutual_clustering
 List cpp_mutual_clustering(const RawMatrix x, const RawMatrix y, const IntegerVector nTip);
 RcppExport SEXP _TreeDist_cpp_mutual_clustering(SEXP xSEXP, SEXP ySEXP, SEXP nTipSEXP) {
@@ -219,6 +231,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_TreeDist_cpp_matching_split_distance", (DL_FUNC) &_TreeDist_cpp_matching_split_distance, 3},
     {"_TreeDist_cpp_jaccard_similarity", (DL_FUNC) &_TreeDist_cpp_jaccard_similarity, 5},
     {"_TreeDist_cpp_msi_distance", (DL_FUNC) &_TreeDist_cpp_msi_distance, 3},
+    {"_TreeDist_cpp_all_pairs_mci", (DL_FUNC) &_TreeDist_cpp_all_pairs_mci, 2},
     {"_TreeDist_cpp_mutual_clustering", (DL_FUNC) &_TreeDist_cpp_mutual_clustering, 3},
     {"_TreeDist_cpp_shared_phylo", (DL_FUNC) &_TreeDist_cpp_shared_phylo, 3},
     {NULL, NULL, 0}

src/SplitList.cpp

@@ -1,7 +1,7 @@
 #include <Rcpp.h>
 using namespace Rcpp;
-#include "tree_distances.h" // Before SplitList.h, for int16
 #include <stdint.h>
+#include "tree_distances.h" // Before SplitList.h, for int16
 #include "SplitList.h"
 
 #define INLASTBIN(n, size) ((size) - ((size) - ((n) % (size))) % (size))
@@ -52,3 +52,131 @@ SplitList::SplitList(RawMatrix x) {
     }
   }
 }
+
+SplitList::SplitList(RawMatrix x, int16 n_tips) {
+  n_splits = x.rows();
+  const int16 
+    n_input_bins = x.cols(),
+    input_bins_per_bin = BIN_SIZE / R_BIN_SIZE
+  ;
+
+  n_bins = (n_input_bins + R_BIN_SIZE - 1) / input_bins_per_bin;
+
+  if (n_bins < 1) throw std::invalid_argument("No leaves present.");
+  if (n_splits < 0) throw std::invalid_argument("Negative number of splits!?");
+  if (n_bins > MAX_BINS) {
+    throw std::length_error("This many leaves cannot be supported. "            // # nocov
+                            "Please contact the TreeDist maintainer if you "    // # nocov
+                            "need to use more!");                               // # nocov
+  }
+  // Rcout << "\n=== SplitList::SplitList(x, n_tips) ===\n";
+
+  for (int16 split = 0; split != n_splits; split++) {
+    int16 last_bin = n_bins - 1;
+    const int16 raggedy_bins = INLASTBIN(n_input_bins, R_BIN_SIZE);
+
+    // Rcout << n_input_bins << " bins in; " << raggedy_bins << " raggedy bins\n";
+    assert(raggedy_bins <= BIN_SIZE / R_BIN_SIZE);
+    assert(raggedy_bins > 0);
+    state[split][last_bin] = INSUBBIN(last_bin, 0);
+    // Rcout << " State[" << split << "][" << last_bin << "] = " << state[split][last_bin] << ".\n";
+    for (int16 input_bin = 1; input_bin != raggedy_bins; input_bin++) {
+      // Rcout << " state [" << split << "][" << last_bin
+      //       << "] increased from " << state[split][last_bin];
+      state[split][last_bin] += INBIN(input_bin, last_bin);
+      // Rcout << " by " << INBIN(input_bin, last_bin) << "\n";
+    }
+
+    const bool invert = state[split][last_bin] & splitbit(1);
+    if (invert) {
+      const int16 last_bin_tips = INLASTBIN(n_tips, BIN_SIZE);
+      // Rcout << "Last bin tips: " << last_bin_tips << ".\n";
+      const splitbit last_mask = last_bin_tips == sizeof(splitbit) * 8 ?
+        splitbit(0) - 1 :
+        (splitbit(1) << INLASTBIN(n_tips, BIN_SIZE)) - 1;
+      // Rcout << "Last mask: " << last_mask << ".\n";
+
+      state[split][last_bin] = ~state[split][last_bin] & last_mask;
+    }
+
+    in_split[split] = count_bits(state[split][last_bin]);
+
+    for (int16 bin = 0; bin != n_bins - 1; bin++) {
+
+      // Rcout << "Split " << split << ", bin << " << bin << ".\n";
+      state[split][bin] = INSUBBIN(bin, 0);
+
+      for (int16 input_bin = 1; input_bin != input_bins_per_bin; input_bin++) {
+        // Rcout << "   Adding " << INBIN(input_bin, bin) << " = " 
+        //       << (splitbit (x(split, (bin * input_bins_per_bin) + input_bin))) << " << "
+        //       << (R_BIN_SIZE * input_bin) << " to state [" << split << "]["
+        //       << bin << "], was " << state[split][bin] << "\n";
+
+        state[split][bin] += INBIN(input_bin, bin);
+      }
+
+      if (invert) {
+        state[split][bin] = ~state[split][bin];
+      }
+
+      in_split[split] += count_bits(state[split][bin]);
+    }
+  }
+}
+
+void SplitList::swap(int16 a, int16 b) {
+  for (int16 bin = n_bins; bin--; ) {
+    const splitbit tmp = state[a][bin];
+    state[a][bin] = state[b][bin];
+    state[b][bin] = tmp;
+  }
+  const int16 tmp = in_split[a];
+  in_split[a] = in_split[b];
+  in_split[b] = tmp;
+}
+
+bool SplitList::less_than(int16 a, int16 b) {
+  for (int16 bin = 0; bin != n_bins; ++bin) {
+    if (state[a][bin] < state[b][bin]) {
+      return true;
+    }
+  }
+  return false;
+}
+
+bool SplitList::greater_than(int16 a, int16 b) {
+  for (int16 bin = 0; bin != n_bins; ++bin) {
+    if (state[a][bin] > state[b][bin]) {
+      return true;
+    }
+  }
+  return false;
+}
+
+int16 SplitList::partition(int16 lo, int16 hi) {
+  const int16 pivot = (hi + lo) / 2;
+  for (int16 i = lo - 1, j = hi + 1; ; ) {
+    do {
+      ++i;
+    } while (less_than(i, pivot));
+    do {
+      --j;
+    } while (greater_than(j, pivot));
+    if (i >= j) {
+      return j;
+    }
+    swap(i, j);
+  }
+}
+
+void SplitList::quicksort(int16 lo, int16 hi) {
+  if (lo < hi) {
+    const int16 p = partition(lo, hi);
+    quicksort(lo, p);
+    quicksort(p + 1, hi);
+  }
+}
+
+void SplitList::quicksort() {
+  quicksort(0, n_splits - 1);
+}

src/SplitList.h

@@ -3,16 +3,26 @@
 
 #include <stdint.h>
 #include <Rcpp.h>
+#include "ints.h"
+#include "tree_distances.h"
+
 using namespace Rcpp;
 
 const int16 R_BIN_SIZE = 8;
 
 class SplitList {
+  int16 partition(int16 lo, int16 hi);
+  void swap(int16 a, int16 b);
+  bool less_than(int16 a, int16 b);
+  bool greater_than(int16 a, int16 b);
 public:
   int16 n_splits, n_bins;
   int16 in_split[MAX_SPLITS];
   splitbit state[MAX_SPLITS][MAX_BINS];
-  SplitList(RawMatrix);
+  SplitList(RawMatrix x);
+  SplitList(RawMatrix x, int16 n_tips);
+  void quicksort();
+  void quicksort(int16 lo, int16 hi);
 };
 
 #endif