Add fast RF cross-pairs via Day 1985 ClusterTable batch

RobinsonFoulds(trees1, trees2) previously fell through to per-pair
R dispatch via CalculateTreeDistance() -> .mapply(), which was ~27x
slower per pair than the all-pairs path.

Add robinson_foulds_cross_pairs() C++ function in day_1985.cpp using
the same Day 1985 algorithm and ClusterTable encoding as the existing
all-pairs function. Wire it into RobinsonFoulds() as a fast path when
both inputs are tree lists with matching tip labels.

Benchmark: 21x speedup for RobinsonFoulds(trees1, trees2).
Per-pair cost is now 1.84x all-pairs (down from 27x).

Tests: 5 new tests in test-batch_coverage.R covering distance,
similarity, normalization, single-tree fallback, and all-pairs
consistency.

Author: ms609

R/RcppExports.R

@@ -13,6 +13,10 @@ robinson_foulds_all_pairs <- function(tables) {
     .Call(`_TreeDist_robinson_foulds_all_pairs`, tables)
 }
 
+robinson_foulds_cross_pairs <- function(tables_a, tables_b) {
+    .Call(`_TreeDist_robinson_foulds_cross_pairs`, tables_a, tables_b)
+}
+
 consensus_info <- function(trees, phylo, p) {
     .Call(`_TreeDist_consensus_info`, trees, phylo, p)
 }

R/tree_distance_rf.R

@@ -148,10 +148,43 @@ RobinsonFoulds <- function(tree1, tree2 = NULL, similarity = FALSE,
                                 class = "dist")
     }
   } else {
-    unnormalized <- CalculateTreeDistance(RobinsonFouldsSplits, tree1, tree2,
-                                          reportMatching)
-    if (similarity) {
-      unnormalized <- .MaxValue(tree1, tree2, NSplits) - unnormalized
+    # Fast cross-pairs path: batch C++ via ClusterTable (Day 1985).
+    # Only applicable when both inputs are lists of trees with matching
+    # tip labels and reportMatching is not requested.
+    fast <- NULL
+    if (!reportMatching &&
+        !inherits(tree1, c("phylo", "Splits")) &&
+        !inherits(tree2, c("phylo", "Splits")) &&
+        is.null(getOption("TreeDist-cluster"))) {
+      lab1 <- TipLabels(tree1)
+      lab2 <- TipLabels(tree2)
+      if (is.list(lab1)) lab1 <- lab1[[1]]
+      if (is.list(lab2)) lab2 <- lab2[[1]]
+      if (setequal(lab1, lab2)) {
+        ct1 <- as.ClusterTable(tree1, tipLabels = lab1)
+        ct2 <- as.ClusterTable(tree2, tipLabels = lab1)
+        shared <- robinson_foulds_cross_pairs(
+          if (is.list(ct1)) ct1 else list(ct1),
+          if (is.list(ct2)) ct2 else list(ct2)
+        )
+        splits1 <- NSplits(tree1)
+        splits2 <- NSplits(tree2)
+        if (similarity) {
+          fast <- shared + shared
+        } else {
+          fast <- outer(splits1, splits2, "+") - shared - shared
+        }
+        dimnames(fast) <- list(names(tree1), names(tree2))
+      }
+    }
+    if (is.null(fast)) {
+      unnormalized <- CalculateTreeDistance(RobinsonFouldsSplits, tree1, tree2,
+                                            reportMatching)
+      if (similarity) {
+        unnormalized <- .MaxValue(tree1, tree2, NSplits) - unnormalized
+      }
+    } else {
+      unnormalized <- fast
     }
   }
 

src/RcppExports.cpp

@@ -44,6 +44,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// robinson_foulds_cross_pairs
+IntegerMatrix robinson_foulds_cross_pairs(const List& tables_a, const List& tables_b);
+RcppExport SEXP _TreeDist_robinson_foulds_cross_pairs(SEXP tables_aSEXP, SEXP tables_bSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const List& >::type tables_a(tables_aSEXP);
+    Rcpp::traits::input_parameter< const List& >::type tables_b(tables_bSEXP);
+    rcpp_result_gen = Rcpp::wrap(robinson_foulds_cross_pairs(tables_a, tables_b));
+    return rcpp_result_gen;
+END_RCPP
+}
 // consensus_info
 double consensus_info(const List trees, const LogicalVector phylo, const NumericVector p);
 RcppExport SEXP _TreeDist_consensus_info(SEXP treesSEXP, SEXP phyloSEXP, SEXP pSEXP) {
@@ -620,6 +632,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_TreeDist_binary_entropy_counts", (DL_FUNC) &_TreeDist_binary_entropy_counts, 2},
     {"_TreeDist_COMCLUST", (DL_FUNC) &_TreeDist_COMCLUST, 1},
     {"_TreeDist_robinson_foulds_all_pairs", (DL_FUNC) &_TreeDist_robinson_foulds_all_pairs, 1},
+    {"_TreeDist_robinson_foulds_cross_pairs", (DL_FUNC) &_TreeDist_robinson_foulds_cross_pairs, 2},
     {"_TreeDist_consensus_info", (DL_FUNC) &_TreeDist_consensus_info, 3},
     {"_TreeDist_HMI_xptr", (DL_FUNC) &_TreeDist_HMI_xptr, 2},
     {"_TreeDist_HH_xptr", (DL_FUNC) &_TreeDist_HH_xptr, 1},

src/day_1985.cpp

@@ -318,6 +318,110 @@ IntegerVector robinson_foulds_all_pairs(const List& tables) {
   return shared;
 }
 
+// Cross-pairs variant: compute RF shared splits between two collections.
+// tables_a and tables_b are lists of ClusterTable XPtr objects
+// (all trees must share the same tip labels).
+// Returns an nA × nB IntegerMatrix of shared-split counts.
+// [[Rcpp::export]]
+IntegerMatrix robinson_foulds_cross_pairs(const List& tables_a,
+                                          const List& tables_b) {
+  const int nA = static_cast<int>(tables_a.size());
+  const int nB = static_cast<int>(tables_b.size());
+  if (nA == 0 || nB == 0) return IntegerMatrix(nA, nB);
+
+  std::vector<ClusterTable*> tbl_a, tbl_b;
+  tbl_a.reserve(nA);
+  tbl_b.reserve(nB);
+  for (int i = 0; i < nA; ++i) {
+    Rcpp::XPtr<ClusterTable> xp = tables_a[i];
+    tbl_a.push_back(xp.get());
+  }
+  for (int j = 0; j < nB; ++j) {
+    Rcpp::XPtr<ClusterTable> xp = tables_b[j];
+    tbl_b.push_back(xp.get());
+  }
+
+  IntegerMatrix result(nA, nB);
+
+  const int32 n_tip = tbl_a[0]->N();
+  StackEntry* S_start;
+  std::array<StackEntry, ct_stack_threshold> S_stack;
+  std::vector<StackEntry> S_heap;
+  if (n_tip <= ct_stack_threshold) {
+    S_start = S_stack.data();
+  } else {
+    S_heap.resize(n_tip);
+    S_start = S_heap.data();
+  }
+
+  for (int i = 0; i < nA; ++i) {
+
+    ClusterTable* Xi = tbl_a[i];
+
+    for (int j = 0; j < nB; ++j) {
+
+      int32 v;
+      int32 w;
+      int32 n_shared = 0;
+
+      ClusterTable* Tj = tbl_b[j];
+
+      StackEntry* S_top = S_start;
+
+      Tj->TRESET();
+      Tj->NVERTEX_short(&v, &w);
+
+      while (v) {
+        if (Tj->is_leaf(v)) {
+          const auto enc_v = Xi->ENCODE(v);
+          *S_top++ = {enc_v, enc_v, 1, 1};
+        } else {
+          const StackEntry& entry = *--S_top;
+          int32 L = entry.L;
+          int32 R = entry.R;
+          int32 N = entry.N;
+          const int32 W_i = entry.W;
+          int32 W = 1 + W_i;
+
+          w -= W_i;
+
+          if (w) {
+            const StackEntry& entry = *--S_top;
+            const int32 W_i = entry.W;
+
+            L = std::min(L, entry.L);
+            R = std::max(R, entry.R);
+            N += entry.N;
+            W += W_i;
+            w -= W_i;
+
+            while (w) {
+              const StackEntry& entry = *--S_top;
+              const int32 W_i = entry.W;
+
+              L = std::min(L, entry.L);
+              R = std::max(R, entry.R);
+              N += entry.N;
+              W += W_i;
+              w -= W_i;
+            }
+          }
+
+          *S_top++ = {L, R, N, W};
+
+          if (N == R - L + 1) {
+            if (Xi->ISCLUST(L, R)) ++n_shared;
+          }
+        }
+        Tj->NVERTEX_short(&v, &w);
+      }
+      result(i, j) = n_shared;
+    }
+  }
+
+  return result;
+}
+
 // [[Rcpp::export]]
 double consensus_info(const List trees, const LogicalVector phylo,
                       const NumericVector p) {

tests/testthat/test-batch_coverage.R

@@ -363,6 +363,56 @@ test_that("Issue #162: CID by-hand matches function for 33-taxon trees", {
   expect_equal(d_fn, d_hand, tolerance = 1e-10)
 })
 
+# RobinsonFoulds cross-pairs (Day 1985 ClusterTable batch) ----
+
+test_that("RF cross-pairs fast path agrees with per-pair", {
+  cross <- RobinsonFoulds(tA, tB)
+  expect_equal(dim(cross), c(5L, 5L))
+  expect_equal(cross[1, 1],
+               RobinsonFoulds(tA[[1]], tB[[1]]),
+               tolerance = 0)
+  expect_equal(cross[3, 2],
+               RobinsonFoulds(tA[[3]], tB[[2]]),
+               tolerance = 0)
+  expect_equal(cross[5, 5],
+               RobinsonFoulds(tA[[5]], tB[[5]]),
+               tolerance = 0)
+})
+
+test_that("RF cross-pairs similarity mode agrees with per-pair", {
+  cross_sim <- RobinsonFoulds(tA, tB, similarity = TRUE)
+  expect_equal(cross_sim[2, 3],
+               RobinsonFoulds(tA[[2]], tB[[3]], similarity = TRUE),
+               tolerance = 0)
+})
+
+test_that("RF cross-pairs normalization agrees with per-pair", {
+  cross_norm <- RobinsonFoulds(tA, tB, normalize = TRUE)
+  expect_equal(cross_norm[1, 4],
+               RobinsonFoulds(tA[[1]], tB[[4]], normalize = TRUE),
+               tolerance = 1e-10)
+})
+
+test_that("RF cross-pairs handles single-tree inputs via fallback", {
+  # Single tree1 → falls back to CalculateTreeDistance
+  d <- RobinsonFoulds(tA[[1]], tB)
+  expect_equal(length(d), 5L)
+  expect_equal(d[1], RobinsonFoulds(tA[[1]], tB[[1]]), tolerance = 0)
+})
+
+test_that("RF cross-pairs all-pairs matches cross-pairs diagonal", {
+  all_dist <- as.matrix(RobinsonFoulds(trees20))
+  # Cross-pairs of trees20 with itself should match all-pairs
+  cross <- RobinsonFoulds(trees20, trees20)
+  expect_equal(dim(cross), c(10L, 10L))
+  for (i in 1:10) {
+    for (j in 1:10) {
+      expect_equal(cross[i, j], all_dist[i, j], tolerance = 0,
+                   info = paste0("i=", i, " j=", j))
+    }
+  }
+})
+
 test_that("Issue #162: batch path agrees with per-pair for 33-taxon trees", {
   tree1 <- ape::read.tree(text = "(B,A,((AG,AF),((((C,(E,D)),((F,G),H)),((K,I),J)),((Q,R),((((AE,(AC,AD)),(AB,(N,(P,O)))),(((Y,(Z,AA)),(W,X)),(V,(T,U)))),((M,S),L))))));")
   tree2 <- ape::read.tree(text = "(B,A,((AG,AF),((((C,(E,D)),((F,G),H)),(((Q,R),((((AE,(AC,AD)),AB),(N,(P,O))),((Y,(Z,AA)),((V,(T,U)),(W,X))))),((M,L),S))),((K,I),J))));")