From c566b616fdc98544815081402c97fbde5ef3f9c8 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Fri, 30 Oct 2026 20:46:52 +0000 Subject: [PATCH 1/2] NNIDist(): score all pairs in one C++ call Renumber and postorder each tree once, compute each tree's splits once, and score every pair in cpp_nni_distance_all_pairs(). Output identical to CompareAll(trees, NNIDist). Co-Authored-By: Claude Opus 5.5 --- .AGENTS/memory/Optimizations.md | 31 +++ NEWS.md | 5 + R/RcppExports.R | 4 + R/tree_distance_nni.R | 38 ++- src/RcppExports.cpp | 13 + src/nni_distance.cpp | 320 ++++++++++++++++-------- tests/testthat/test-tree_distance_nni.R | 25 ++ 7 files changed, 327 insertions(+), 109 deletions(-) diff --git a/.AGENTS/memory/Optimizations.md b/.AGENTS/memory/Optimizations.md index 2a553898..62351c90 100644 --- a/.AGENTS/memory/Optimizations.md +++ b/.AGENTS/memory/Optimizations.md @@ -1020,3 +1020,34 @@ pay; what pays is what the parent overlaps let a candidate *skip*. double matrix once the fill is no longer dominant. - Next: depth-2 pair candidates (two leaves) by the same route; the unbounded MCI could use the Tier 1 form if tests tolerate ulp changes. + +## NNIDist() all pairs in one C++ call (2026-10-30) + +- `NNIDist(trees)` (one list, all pairs) went through `CompareAll()` -> + `.NNIDistSingle()`, which per pair ran `RenumberTips()`, two `Postorder()` + calls and `cpp_nni_distance()`; on GHA ~85% of the time was this R overhead. +- Now `.NNIDistAllPairs()` renumbers and postorders each tree once and calls + `cpp_nni_distance_all_pairs()`, which computes each tree's splits once + (under its own rootedness) and per pair runs only the matching and scoring. + Shared core: `nni_tree()`, `nni_edge_to_splits()`, `nni_rf_matching()`, + `nni_score()`, also used by `cpp_nni_distance()`. +- Identity: a pair is scored under the *first* tree's rootedness (as before), + so when two trees with equal edge counts differ in rootedness (polytomies), + the second tree's splits are recomputed for that pair. Output `identical()` + to `CompareAll(trees, NNIDist)` on all 945 7-leaf trees, unrooted 7-leaf, + random 60-200 leaves with shuffled labels, 3 leaves, polytomies, mixed + rootedness. +- Release installs, min of N, single thread, R-devel (2026-06-29) on Windows, + same machine/session (main via `git archive origin/main`): + + | Workload | main | PR | + |---|---|---| + | 945 × 7 leaves (`as.phylo(0:944, 7)`) | 136.7 s (n = 1) | 0.23 s (n = 3) | + | 100 × 50 leaves | 1.80 s (n = 3) | 0.040 s (n = 10) | + | 40 × 200 leaves | 0.86 s (n = 3) | 0.110 s (n = 10) | + + Main's 945-tree time is inflated by `CompareAll()` indexing a multiPhylo + (PR #227); GHA measured 46 s on R 4.6, against ~7 s for the bare C++ calls. +- Next: at 200 leaves the O(s^2 · bins) split matching now dominates. + Canonicalizing splits (tip-1 bit clear) and sorting each tree's splits once + would allow a linear merge per pair. diff --git a/NEWS.md b/NEWS.md index df06e482..e968e1f7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -35,6 +35,11 @@ ## Performance +- `NNIDist()` on a single list of trees prepares each tree once and scores + all pairs in one C++ call, rather than renumbering and reordering both + trees for every pair: c. 8–45× faster on 40–100 trees, and 945 seven-leaf + trees take 0.2 s. + - Trees with different leaf sets are pruned to their shared leaves in C++, without constructing pruned trees in R: `TreeDistance()` on 150 gene trees is c. 15× faster, with lower memory use (#123). diff --git a/R/RcppExports.R b/R/RcppExports.R index 2de3878f..2a76846f 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -98,6 +98,10 @@ cpp_nni_distance <- function(edge1, edge2, nTip) { .Call(`_TreeDist_cpp_nni_distance`, edge1, edge2, nTip) } +cpp_nni_distance_all_pairs <- function(edges, nTip) { + .Call(`_TreeDist_cpp_nni_distance_all_pairs`, edges, nTip) +} + min_subtree_entropy <- function(n) { .Call(`_TreeDist_min_subtree_entropy`, n) } diff --git a/R/tree_distance_nni.R b/R/tree_distance_nni.R index 7924a248..e2098e48 100644 --- a/R/tree_distance_nni.R +++ b/R/tree_distance_nni.R @@ -78,7 +78,40 @@ #' @family tree distances #' @export NNIDist <- function(tree1, tree2 = tree1) { - .TreeDistance(.NNIDistSingle, tree1, tree2) + if (!inherits(tree1, "phylo") && length(tree1) > 1 && + identical(tree1, tree2)) { + .NNIDistAllPairs(tree1) + } else { + .TreeDistance(.NNIDistSingle, tree1, tree2) + } +} + +.StopEdgeMismatch <- function() { + stop("Both trees must have the same number of edges. ", + "Is one rooted and the other unrooted?") +} + +#' @importFrom TreeTools Postorder RenumberTips TipLabels +#' @importFrom stats setNames +.NNIDistAllPairs <- function(trees) { + labelList <- lapply(trees, TipLabels) + .CheckLabelsSame(labelList) + tipLabels <- as.character(labelList[[1]]) + nTip <- length(tipLabels) + .CheckMaxTips(nTip) + edges <- lapply(trees, function(tree) { + Postorder(RenumberTips(tree, tipLabels)[["edge"]]) + }) + if (length(unique(vapply(edges, nrow, integer(1)))) != 1) { + .StopEdgeMismatch() + } + bounds <- cpp_nni_distance_all_pairs(edges, nTip) + + # Return: + lapply(setNames(nm = rownames(bounds)), function(bound) { + structure(bounds[bound, ], Size = length(trees), Labels = names(trees), + Diag = FALSE, Upper = TRUE, class = "dist") + }) } #' @importFrom TreeTools Postorder RenumberTips @@ -86,8 +119,7 @@ NNIDist <- function(tree1, tree2 = tree1) { .NNIDistSingle <- function(tree1, tree2, nTip, ...) { .CheckMaxTips(nTip) if (nrow(tree1[["edge"]]) != nrow(tree2[["edge"]])) { - stop("Both trees must have the same number of edges. ", - "Is one rooted and the other unrooted?") + .StopEdgeMismatch() } tree2 <- RenumberTips(tree2, as.character(tree1[["tip.label"]])) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b127ca30..e5acdc03 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -269,6 +269,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cpp_nni_distance_all_pairs +IntegerMatrix cpp_nni_distance_all_pairs(const List& edges, const IntegerVector& nTip); +RcppExport SEXP _TreeDist_cpp_nni_distance_all_pairs(SEXP edgesSEXP, SEXP nTipSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const List& >::type edges(edgesSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type nTip(nTipSEXP); + rcpp_result_gen = Rcpp::wrap(cpp_nni_distance_all_pairs(edges, nTip)); + return rcpp_result_gen; +END_RCPP +} // min_subtree_entropy Rcpp::List min_subtree_entropy(const int n); RcppExport SEXP _TreeDist_min_subtree_entropy(SEXP nSEXP) { @@ -896,6 +908,7 @@ static const R_CallMethodDef CallEntries[] = { {"_TreeDist_lapjv", (DL_FUNC) &_TreeDist_lapjv, 2}, {"_TreeDist_cpp_mast", (DL_FUNC) &_TreeDist_cpp_mast, 3}, {"_TreeDist_cpp_nni_distance", (DL_FUNC) &_TreeDist_cpp_nni_distance, 3}, + {"_TreeDist_cpp_nni_distance_all_pairs", (DL_FUNC) &_TreeDist_cpp_nni_distance_all_pairs, 2}, {"_TreeDist_min_subtree_entropy", (DL_FUNC) &_TreeDist_min_subtree_entropy, 1}, {"_TreeDist_nni_max_step", (DL_FUNC) &_TreeDist_nni_max_step, 2}, {"_TreeDist_cpp_mutual_clustering_all_pairs", (DL_FUNC) &_TreeDist_cpp_mutual_clustering_all_pairs, 3}, diff --git a/src/nni_distance.cpp b/src/nni_distance.cpp index 3a14ffd7..9d367ec0 100644 --- a/src/nni_distance.cpp +++ b/src/nni_distance.cpp @@ -2,9 +2,11 @@ #include #include +#include // for max #include #include // for unique_ptr, uninitialized_fill_n #include // for once_flag +#include #include "tree_distances.h" #include "li_diameters.h" @@ -177,20 +179,54 @@ void update_score(const int32_t subtree_edges, int32_t *lower_bound, } -// Edges must be listed in postorder -inline void nni_edge_to_splits(const IntegerMatrix& edge, +constexpr int32_t NOT_TRIVIAL = std::numeric_limits::max(); +constexpr int32_t N_NNI_BOUNDS = 6; + +// Edges in postorder, with 1-based node numbers +struct NNITree { + std::vector parent; + std::vector child; + int32_t root; + bool rooted; +}; + +inline NNITree nni_tree(const IntegerMatrix& edge) { + const int32_t n_edge = int32_t(edge.nrow()); + NNITree tree; + tree.parent.resize(n_edge); + tree.child.resize(n_edge); + for (int32_t i = 0; i < n_edge; ++i) { + tree.parent[i] = static_cast(edge(i, 0)); + tree.child[i] = static_cast(edge(i, 1)); + } + tree.root = tree.parent[n_edge - 1]; + tree.rooted = tree.parent[n_edge - 3] != tree.root; + return tree; +} + +// Node (0-based) whose split duplicates the root split in a rooted tree +inline int32_t nni_trivial_two(const NNITree& tree, const bool rooted) { + return rooted ? tree.child.back() - 1 : NOT_TRIVIAL; +} + +inline int32_t nni_n_splits(const int32_t n_edge, const int32_t n_tip, + const bool rooted) { + return n_edge - (rooted ? 1 : 0) - n_tip; +} + +inline void nni_edge_to_splits(const NNITree& tree, const int32_t n_tip, - const int32_t n_edge, - const int32_t n_node, const int32_t n_bin, - const int32_t trivial_origin, const int32_t trivial_two, - std::unique_ptr& splits, - std::unique_ptr& names) { - + uint64_t* splits, + int32_t* names) { + ASSERT((n_bin == n_tip + SL_BIN_SIZE - 1) / SL_BIN_SIZE); + const int32_t n_edge = int32_t(tree.parent.size()); + const int32_t n_node = n_edge + 1; + const int32_t trivial_origin = tree.root - 1; std::vector tmp_splits(n_node * n_bin, 0); - + for (int32_t i = 0; i < n_tip; ++i) { const size_t bin_idx = static_cast(i / SL_BIN_SIZE); const size_t idx = i * n_bin + bin_idx; @@ -198,15 +234,15 @@ inline void nni_edge_to_splits(const IntegerMatrix& edge, ASSERT((i % SL_BIN_SIZE) < 64); // 1 << 64 = 0 tmp_splits[idx] = static_cast(1) << (i % SL_BIN_SIZE); } - + for (int32_t i = 0; i < n_edge - 1; ++i) { /* final edge is second root edge */ - const ptrdiff_t parent_row_i = (edge(i, 0) - 1) * n_bin; - const ptrdiff_t child_row_i = (edge(i, 1) - 1) * n_bin; + const ptrdiff_t parent_row_i = (tree.parent[i] - 1) * n_bin; + const ptrdiff_t child_row_i = (tree.child[i] - 1) * n_bin; for (int32_t j = 0; j < n_bin; ++j) { tmp_splits[parent_row_i + j] |= tmp_splits[child_row_i + j]; } } - + int32_t n_trivial = 0; for (int32_t i = n_tip; i < n_node; ++i) { if (i == trivial_origin || i == trivial_two) { @@ -222,28 +258,28 @@ inline void nni_edge_to_splits(const IntegerMatrix& edge, } } } - -grf_match nni_rf_matching ( - const std::unique_ptr& a, - const std::unique_ptr& b, + +void nni_rf_matching ( + const uint64_t* a, + const uint64_t* b, const int32_t n_splits, const int32_t n_bins, - const int32_t n_tips) { - - ASSERT(n_splits > 0); + const int32_t n_tips, + grf_match& matching) { + ASSERT(n_tips > 3); ASSERT(n_tips <= NNI_MAX_TIPS); // Validated by exported wrapper - + const int32_t last_bin = n_bins - 1; - const int32_t unset_tips = (n_tips % SL_BIN_SIZE) ? + const int32_t unset_tips = (n_tips % SL_BIN_SIZE) ? SL_BIN_SIZE - n_tips % SL_BIN_SIZE : 0; - + const uint64_t unset_mask = ALL_ONES >> unset_tips; - - grf_match matching(n_splits, NA_INT32); + + matching.assign(n_splits, NA_INT32); HybridBuffer b_complement(n_splits * n_bins); - + for (int32_t i = 0; i < n_splits; ++i) { const int32_t row_i = i * n_bins; for (int32_t bin = 0; bin < last_bin; ++bin) { @@ -253,15 +289,15 @@ grf_match nni_rf_matching ( const int32_t last_cell = row_i + last_bin; b_complement[last_cell] = b[last_cell] ^ unset_mask; } - + for (int32_t ai = 0; ai < n_splits; ++ai) { const int32_t a_row = ai * n_bins; for (int32_t bi = 0; bi < n_splits; ++bi) { - + bool all_match = true; bool all_complement = true; const int32_t b_row = bi * n_bins; - + for (int32_t bin = 0; bin < n_bins; ++bin) { const int32_t a_cell = int32_t(a_row) + int32_t(bin); const int32_t b_cell = int32_t(b_row) + int32_t(bin); @@ -286,81 +322,32 @@ grf_match nni_rf_matching ( } } } - - return matching; } -// [[Rcpp::export]] -IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, - const IntegerMatrix& edge2, - const IntegerVector& nTip) { - - ASSERT(nTip[0] <= NNI_MAX_TIPS && "Cannot calculate NNI distance for trees with so many tips."); - const int32_t n_tip = static_cast(nTip[0]); +// Writes lower, best_upper, tight_upper, loose_upper, fack_upper and li_upper +// to `bounds`. `names_1` gives the node that defines each split of tree1; +// `match` maps these splits to those of the other tree. +void nni_score(const NNITree& tree1, const int32_t* names_1, + const grf_match& match, const int32_t n_tip, + int* bounds) { const int32_t node_0 = n_tip; const int32_t node_0_r = n_tip + 1; - const int32_t n_edge = int32_t(edge1.nrow()); - + const int32_t n_edge = int32_t(tree1.parent.size()); + const bool rooted = tree1.rooted; + int32_t lower_bound = 0; int32_t tight_score_bound = 0; - + int32_t loose_score_bound = 0; int32_t best_upper_bound = 0; int32_t fack_score_bound = 0; int32_t li_score_bound = 0; - - ASSERT(n_edge == int32_t(edge2.nrow()) - && "Both trees must have the same number of edges."); - if (n_tip < 4) { - return(IntegerVector::create(Named("lower") = 0, - _["best_upper"] = 0, - _["tight_upper"] = 0, - _["loose_upper"] = 0, - _["fack_upper"] = 0, - _["li_upper"] = 0)); - } - - const int32_t root_1 = static_cast(edge1(n_edge - 1, 0)); - const int32_t root_2 = static_cast(edge2(n_edge - 1, 0)); - bool rooted = static_cast(edge1(n_edge - 3, 0)) != root_1; - constexpr int32_t NOT_TRIVIAL = std::numeric_limits::max(); - const int32_t n_node = n_edge + 1; - const int32_t n_bin = static_cast(((n_tip - 1) / SL_BIN_SIZE) + 1); - const int32_t trivial_origin_1 = root_1 - 1; - const int32_t trivial_origin_2 = root_2 - 1; - const int32_t trivial_two_1 = rooted ? - static_cast(edge1(n_edge - 1, 1)) - 1 : - NOT_TRIVIAL; - const int32_t trivial_two_2 = rooted ? - static_cast(edge2(n_edge - 1, 1)) - 1 : - NOT_TRIVIAL; const int32_t n_distinct_edge = int32_t(n_edge - (rooted ? 1 : 0)); - const int32_t n_splits = n_distinct_edge - n_tip; - - if (n_splits < 1) { - ASSERT(false && "NNI distance is undefined for trees with no splits"); // #nocov - } - - std::unique_ptr splits1(new uint64_t[n_splits * n_bin]); - std::unique_ptr splits2(new uint64_t[n_splits * n_bin]); - // std::vector names_1; - // names_1.reserve(n_splits); - std::unique_ptr names_1(new int32_t[n_splits]); - ensure_cache(); - - if (n_edge != n_tip && n_tip > 3) { - nni_edge_to_splits(edge2, n_tip, n_edge, n_node, n_bin, - trivial_origin_2, trivial_two_2, splits2, names_1); - nni_edge_to_splits(edge1, n_tip, n_edge, n_node, n_bin, - trivial_origin_1, trivial_two_1, splits1, names_1); - } // else no internal nodes resolved - - grf_match match = nni_rf_matching(splits1, splits2, n_splits, n_bin, n_tip); - + HybridBuffer matched_1(n_tip); HybridBuffer unmatched_below(n_tip); - + const int32_t match_size = static_cast(match.size()); for (int32_t i = 0; i < match_size; ++i) { ASSERT(n_edge != n_tip && n_tip > 3); // else names_1 uninitialized @@ -372,11 +359,11 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, matched_1[node_i] = true; } } - + const int32_t edges_to_check = n_distinct_edge - (rooted ? 1 : 0); for (int32_t i = 0; i < edges_to_check; ++i) { - const int32_t parent_i = static_cast(edge1(i, 0)) - 1; - const int32_t child_i = static_cast(edge1(i, 1)) - 1; + const int32_t parent_i = tree1.parent[i] - 1; + const int32_t child_i = tree1.child[i] - 1; // If edge is unmatched, add one to subtree size. if (child_i >= n_tip) { if (!matched_1[child_i - node_0]) { @@ -388,13 +375,13 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, } } } - + // Root edges: - const int32_t root_node = root_1 - node_0_r; - + const int32_t root_node = tree1.root - node_0_r; + if (rooted) { - const int32_t root_child_1 = static_cast(edge1(n_edge - 1, 1)) - 1; - const int32_t root_child_2 = static_cast(edge1(n_edge - 2, 1)) - 1; + const int32_t root_child_1 = tree1.child[n_edge - 1] - 1; + const int32_t root_child_2 = tree1.child[n_edge - 2] - 1; const int32_t unmatched_1 = root_child_1 < n_tip ? 0 : unmatched_below[root_child_1 - node_0]; if (root_child_2 >= n_tip) { @@ -425,12 +412,133 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, &li_score_bound, &fack_score_bound); } - return IntegerVector::create( - Named("lower") = lower_bound, - _["best_upper"] = best_upper_bound, - _["tight_upper"] = tight_score_bound == NA_INT32 - ? NA_INTEGER : tight_score_bound, - _["loose_upper"] = loose_score_bound, - _["fack_upper"] = fack_score_bound, - _["li_upper"] = li_score_bound); + bounds[0] = lower_bound; + bounds[1] = best_upper_bound; + bounds[2] = tight_score_bound == NA_INT32 ? NA_INTEGER : tight_score_bound; + bounds[3] = loose_score_bound; + bounds[4] = fack_score_bound; + bounds[5] = li_score_bound; +} + +inline CharacterVector nni_bound_names() { + return CharacterVector::create("lower", "best_upper", "tight_upper", + "loose_upper", "fack_upper", "li_upper"); +} + +// [[Rcpp::export]] +IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, + const IntegerMatrix& edge2, + const IntegerVector& nTip) { + + ASSERT(nTip[0] <= NNI_MAX_TIPS && "Cannot calculate NNI distance for trees with so many tips."); + const int32_t n_tip = static_cast(nTip[0]); + const int32_t n_edge = int32_t(edge1.nrow()); + + ASSERT(n_edge == int32_t(edge2.nrow()) + && "Both trees must have the same number of edges."); + + IntegerVector ret(N_NNI_BOUNDS); + ret.names() = nni_bound_names(); + if (n_tip < 4) { + return ret; + } + + const NNITree tree1 = nni_tree(edge1); + const NNITree tree2 = nni_tree(edge2); + const bool rooted = tree1.rooted; + const int32_t n_bin = static_cast(((n_tip - 1) / SL_BIN_SIZE) + 1); + const int32_t n_splits = nni_n_splits(n_edge, n_tip, rooted); + + if (n_splits < 1) { + ASSERT(false && "NNI distance is undefined for trees with no splits"); // #nocov + } + + std::unique_ptr splits1(new uint64_t[n_splits * n_bin]); + std::unique_ptr splits2(new uint64_t[n_splits * n_bin]); + std::unique_ptr names_1(new int32_t[n_splits]); + ensure_cache(); + + if (n_edge != n_tip) { + nni_edge_to_splits(tree2, n_tip, n_bin, nni_trivial_two(tree2, rooted), + splits2.get(), names_1.get()); + nni_edge_to_splits(tree1, n_tip, n_bin, nni_trivial_two(tree1, rooted), + splits1.get(), names_1.get()); + } // else no internal nodes resolved + + grf_match match; + nni_rf_matching(splits1.get(), splits2.get(), n_splits, n_bin, n_tip, match); + nni_score(tree1, names_1.get(), match, n_tip, INTEGER(ret)); + + return ret; +} + +// One column per pair of trees, in the order of a `dist` object. +// Edges must be listed in postorder, with leaves numbered consistently. +// [[Rcpp::export]] +IntegerMatrix cpp_nni_distance_all_pairs(const List& edges, + const IntegerVector& nTip) { + const int32_t n_tip = static_cast(nTip[0]); + ASSERT(n_tip <= NNI_MAX_TIPS); + const R_xlen_t n_tree = edges.size(); + const R_xlen_t n_pair = n_tree * (n_tree - 1) / 2; + + IntegerMatrix ret(N_NNI_BOUNDS, n_pair); + rownames(ret) = nni_bound_names(); + if (n_tip < 4 || n_pair == 0) { + return ret; + } + ensure_cache(); + + std::vector trees; + trees.reserve(n_tree); + for (R_xlen_t i = 0; i < n_tree; ++i) { + trees.push_back(nni_tree(as(edges[i]))); + } + + const int32_t n_edge = int32_t(trees[0].parent.size()); + const int32_t n_bin = static_cast(((n_tip - 1) / SL_BIN_SIZE) + 1); + const size_t max_splits = size_t(std::max(n_edge - n_tip, 0)); + const size_t stride = max_splits * n_bin; + + // Splits of each tree, computed once under its own rootedness + std::vector splits(stride * n_tree); + std::vector names(max_splits * n_tree); + if (n_edge != n_tip) { + for (R_xlen_t i = 0; i < n_tree; ++i) { + nni_edge_to_splits(trees[i], n_tip, n_bin, + nni_trivial_two(trees[i], trees[i].rooted), + &splits[stride * i], &names[max_splits * i]); + } + } + + // A pair is scored under the first tree's rootedness, so a second tree + // that differs needs its splits recomputing + std::vector recast(stride); + std::vector recast_names(max_splits); + + grf_match match; + int* out = INTEGER(ret); + R_xlen_t pair = 0; + for (R_xlen_t i = 0; i < n_tree - 1; ++i) { + const NNITree& tree_i = trees[i]; + const int32_t n_splits = nni_n_splits(n_edge, n_tip, tree_i.rooted); + const uint64_t* splits_i = &splits[stride * i]; + for (R_xlen_t j = i + 1; j < n_tree; ++j) { + const uint64_t* splits_j = &splits[stride * j]; + if (trees[j].rooted != tree_i.rooted) { + nni_edge_to_splits(trees[j], n_tip, n_bin, + nni_trivial_two(trees[j], tree_i.rooted), + recast.data(), recast_names.data()); + splits_j = recast.data(); + } + nni_rf_matching(splits_i, splits_j, n_splits, n_bin, n_tip, match); + nni_score(tree_i, &names[max_splits * i], match, n_tip, + out + pair * N_NNI_BOUNDS); + if (++pair % 4096 == 0) { + Rcpp::checkUserInterrupt(); + } + } + } + + return ret; } diff --git a/tests/testthat/test-tree_distance_nni.R b/tests/testthat/test-tree_distance_nni.R index 09ec2f7d..2b18b271 100644 --- a/tests/testthat/test-tree_distance_nni.R +++ b/tests/testthat/test-tree_distance_nni.R @@ -230,6 +230,31 @@ test_that("NNI with lists of trees", { ) }) +test_that("NNIDist() all pairs matches pairwise calculation", { + .ExpectPairwise <- function(trees) { + expect_identical(NNIDist(trees), CompareAll(trees, NNIDist)) + } + set.seed(0) + .ExpectPairwise(as.phylo(0:20, 5)) + .ExpectPairwise(as.phylo(0:3, 3)) + .ExpectPairwise(c(as.phylo(0:3, 8), as.phylo(0:2, 8))) + .ExpectPairwise(lapply(as.phylo(0:20, 6), UnrootTree)) + shuffled <- lapply(1:5, function(i) { + RandomTree(sample(paste0("t", 1:60)), root = i %% 2 == 1) + }) + .ExpectPairwise(lapply(shuffled, UnrootTree)) + .ExpectPairwise(structure(lapply(shuffled, RootTree, "t1"), + names = letters[1:5])) + + # Same number of edges, but only some trees rooted + mixed <- ape::read.tree(text = c("((a, b, c), (d, e));", "(a, b, (c, (d, e)));", + "((a, d, c), (b, e));", "(a, e, (c, (d, b)));")) + .ExpectPairwise(mixed) + + expect_error(NNIDist(list(BalancedTree(5), RootOnNode(BalancedTree(5), 1))), + "same number of edges") +}) + test_that("NNIDiameter() is sane", { exacts <- NNIDiameter(3:12) From 4c12013d04e4200acda86a509f1096844f5ba1d8 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Mon, 2 Nov 2026 09:28:40 +0000 Subject: [PATCH 2/2] nni_distance.cpp: widen index products before multiplying Addresses CodeQL "multiplication result converted to larger type". Co-Authored-By: Claude Opus 5.5 --- src/nni_distance.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/nni_distance.cpp b/src/nni_distance.cpp index 9d367ec0..32f075d7 100644 --- a/src/nni_distance.cpp +++ b/src/nni_distance.cpp @@ -225,19 +225,19 @@ inline void nni_edge_to_splits(const NNITree& tree, const int32_t n_edge = int32_t(tree.parent.size()); const int32_t n_node = n_edge + 1; const int32_t trivial_origin = tree.root - 1; - std::vector tmp_splits(n_node * n_bin, 0); + std::vector tmp_splits(size_t(n_node) * n_bin, 0); for (int32_t i = 0; i < n_tip; ++i) { const size_t bin_idx = static_cast(i / SL_BIN_SIZE); - const size_t idx = i * n_bin + bin_idx; + const size_t idx = size_t(i) * n_bin + bin_idx; ASSERT(idx < tmp_splits.size()); ASSERT((i % SL_BIN_SIZE) < 64); // 1 << 64 = 0 tmp_splits[idx] = static_cast(1) << (i % SL_BIN_SIZE); } for (int32_t i = 0; i < n_edge - 1; ++i) { /* final edge is second root edge */ - const ptrdiff_t parent_row_i = (tree.parent[i] - 1) * n_bin; - const ptrdiff_t child_row_i = (tree.child[i] - 1) * n_bin; + const ptrdiff_t parent_row_i = ptrdiff_t(tree.parent[i] - 1) * n_bin; + const ptrdiff_t child_row_i = ptrdiff_t(tree.child[i] - 1) * n_bin; for (int32_t j = 0; j < n_bin; ++j) { tmp_splits[parent_row_i + j] |= tmp_splits[child_row_i + j]; } @@ -248,7 +248,7 @@ inline void nni_edge_to_splits(const NNITree& tree, if (i == trivial_origin || i == trivial_two) { ++n_trivial; } else { - const ptrdiff_t i_row = i * n_bin; + const ptrdiff_t i_row = ptrdiff_t(i) * n_bin; for (int32_t j = 0; j < n_bin; ++j) { const size_t idx = i - n_tip - n_trivial; ASSERT(idx >= 0 && idx < n_splits); @@ -278,7 +278,7 @@ void nni_rf_matching ( matching.assign(n_splits, NA_INT32); HybridBuffer - b_complement(n_splits * n_bins); + b_complement(size_t(n_splits) * n_bins); for (int32_t i = 0; i < n_splits; ++i) { const int32_t row_i = i * n_bins; @@ -453,8 +453,8 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, ASSERT(false && "NNI distance is undefined for trees with no splits"); // #nocov } - std::unique_ptr splits1(new uint64_t[n_splits * n_bin]); - std::unique_ptr splits2(new uint64_t[n_splits * n_bin]); + std::unique_ptr splits1(new uint64_t[size_t(n_splits) * n_bin]); + std::unique_ptr splits2(new uint64_t[size_t(n_splits) * n_bin]); std::unique_ptr names_1(new int32_t[n_splits]); ensure_cache();