From f39e9f2fbbef11f6a5d2c250de84a35d06ed28ef Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Fri, 23 Oct 2026 09:50:35 +0100 Subject: [PATCH 1/4] Rogue SPR: score leaf-pair candidates from parent split overlaps Extends SingleDrop (now LeafDrop) to depth-2 candidates: removing leaves i and j replays two removals on the parent ledgers, recording the split that vanishes from each tree at each step; child overlaps are parent overlaps less the popcount of shared membership of {i, j} (2-bit masks; the single-leaf kernel is unchanged via a template). Pairs fall back to ScoreWithout() if either leaf is leaf 1, removing either leaf alone makes two splits identical, the trees have splits differing in exactly {i, j}, or reduce_trees() would reduce the smaller trees. Co-Authored-By: Claude Opus 5.5 --- src/split_counts.h | 16 +- src/spr.cpp | 276 +++++++++++++++--------- src/tree_distances.cpp | 21 +- tests/testthat/test-tree_distance_spr.R | 2 +- 4 files changed, 205 insertions(+), 110 deletions(-) diff --git a/src/split_counts.h b/src/split_counts.h index ecb72f0a..ec80c3ca 100644 --- a/src/split_counts.h +++ b/src/split_counts.h @@ -4,13 +4,14 @@ #include // Splits of a pair of trees, each obtained from a parent tree by removing -// one leaf, described through the parent's splits. The trees have `n_split` -// splits each. Child split i of the first tree restricts parent split -// row[i]; it has row_in[i] leaves, and row_has[i] is 1 if the parent split -// contained the removed leaf, else 0. col, col_in and col_has describe the -// second tree likewise. overlap[r * stride + c] counts the leaves shared by -// parent splits r and c, so the corresponding child splits share -// overlap[r * stride + c] - (row_has & col_has). +// one leaf (or two, if `pair`), described through the parent's splits. The +// trees have `n_split` splits each. Child split i of the first tree +// restricts parent split row[i]; it has row_in[i] leaves, and bit k of +// row_has[i] is set if the parent split contained removed leaf k. col, +// col_in and col_has describe the second tree likewise. +// overlap[r * stride + c] counts the leaves shared by parent splits r and c, +// so the corresponding child splits share +// overlap[r * stride + c] - popcount(row_has & col_has). struct SplitCounts { int n_tip; // Leaves in each child tree int n_split; @@ -22,6 +23,7 @@ struct SplitCounts { const int *col; const int *col_in; const int16_t *col_has; + bool pair; // Two leaves removed; else row_has and col_has are 0 or 1 }; #endif diff --git a/src/spr.cpp b/src/spr.cpp index 544d6b8e..26d94170 100644 --- a/src/spr.cpp +++ b/src/spr.cpp @@ -691,26 +691,27 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, return score1 + BoundedCID(sub2[0], sub2[1], cutoff - score1, pp); } -// Scores candidates that drop a single leaf without rebuilding the trees. +// Scores candidates that drop one or two leaves without rebuilding the trees. // // Removing leaf l from a binary tree suppresses its parent node, so exactly // one split vanishes from each tree, and every other split restricts to a -// distinct split of the smaller tree. If, once l is removed, the trees -// still cannot be reduced and still share no split, then the smaller trees -// are those that KeepAndReduce() would return (up to leaf numbering and -// rooting), so their mutual clustering cost matrix follows from the overlaps -// between the parent splits, without building, reducing or splitting trees. -class SingleDrop { +// distinct split of the smaller tree; removing a second leaf repeats this. +// If, once the leaves are removed, the trees still cannot be reduced and +// still share no split, then the smaller trees are those that +// KeepAndReduce() would return (up to leaf numbering and rooting), so their +// mutual clustering cost matrix follows from the overlaps between the parent +// splits, without building, reducing or splitting trees. +class LeafDrop { public: // `edge1` and `edge2`: postorder edges of the parent trees, binary and // rooted on leaf 1. `bits1` and `bits2`: their splits, as from // GetSplits() with identity leaf indexing; `node1` and `node2`: the node // of each split. Disabled if `n_tip` is zero. - SingleDrop(const IntegerMatrix &edge1, const IntegerMatrix &edge2, - const int n_tip, const int n_word, - const std::vector &bits1, - const std::vector &bits2, - const std::vector &node1, const std::vector &node2) + LeafDrop(const IntegerMatrix &edge1, const IntegerMatrix &edge2, + const int n_tip, const int n_word, + const std::vector &bits1, + const std::vector &bits2, + const std::vector &node1, const std::vector &node2) : n_tip_(n_tip), n_word_(n_word), n_split_(int(node1.size())), bits1_(bits1), bits2_(bits2) { ready_ = Init(edge1, edge2, node1, node2); @@ -719,44 +720,56 @@ class SingleDrop { // Prepares to score the removal of `leaf` (0-based); false if it must be // scored by ScoreWithout() instead. bool Prepare(const int leaf) { - const int tip = leaf + 1; - if (!ready_ || tip == 1 || shared_[leaf] || Reduces(tip)) { + if (!ready_ || leaf == 0 || shared_[leaf]) { return false; } - const int drop1 = Vanishing(x_, row_of1_, tip, nullptr); - drop2_ = Vanishing(y_, row_of2_, tip, &merge2_); - const int n_child = n_split_ - 1; - row_.resize(n_child); - row_in_.resize(n_child); - row_has_.resize(n_child); - col_.resize(n_child); - col_in_.resize(n_child); - col_has_.resize(n_child); - const int word = leaf >> 6; - const uint64_t bit = uint64_t(1) << (leaf & 63); - for (int r = 0, k = 0; r != n_split_; ++r) { - if (r == drop1) continue; - const int16_t has = (bits1_[size_t(r) * n_word_ + word] & bit) != 0; - row_[k] = r; - row_has_[k] = has; - row_in_[k] = in1_[r] - has; - ++k; + xc_ = x_; + yc_ = y_; + drop1_[0] = Remove(xc_, row_of1_, leaf + 1, nullptr); + drop2_[0] = Remove(yc_, row_of2_, leaf + 1, &merge2_[0]); + drop1_[1] = drop2_[1] = merge2_[1] = -1; + if (Reduces()) { + return false; } - for (int c = 0, k = 0; c != n_split_; ++c) { - if (c == drop2_) continue; - const int16_t has = (bits2_[size_t(c) * n_word_ + word] & bit) != 0; - col_[k] = c; - col_has_[k] = has; - col_in_[k] = in2_[c] - has; - ++k; + Fill(leaf, -1); + return true; + } + + // Prepares to score the removal of leaves `i` and `j` (0-based, i < j); + // false if they must be scored by ScoreWithout() instead. + bool Prepare(const int i, const int j) { + if (!ready_ || i == 0 || shared_[i] || shared_[j] || + std::binary_search(shared_pair_.begin(), shared_pair_.end(), + PairKey(i, j))) { + return false; } + if (base_leaf_ != i) { + // Trees less leaf i, reused for each j + xi_ = x_; + yi_ = y_; + base_drop1_ = Remove(xi_, row_of1_, i + 1, nullptr); + base_drop2_ = Remove(yi_, row_of2_, i + 1, &base_merge2_); + base_leaf_ = i; + } + xc_ = xi_; + yc_ = yi_; + drop1_[0] = base_drop1_; + drop2_[0] = base_drop2_; + merge2_[0] = base_merge2_; + drop1_[1] = Remove(xc_, row_of1_, j + 1, nullptr); + drop2_[1] = Remove(yc_, row_of2_, j + 1, &merge2_[1]); + if (Reduces()) { + return false; + } + Fill(i, j); return true; } // As ScoreWithout(), for the prepared candidate. double Score(const double cutoff, const ParentPrices *pp) { - const int n_child = n_split_ - 1; - const int n = n_tip_ - 1; + const int n_drop = drop1_[1] < 0 ? 1 : 2; + const int n_child = n_split_ - n_drop; + const int n = n_tip_ - n_drop; // As GetSplits() const double n_lg2_n = n * std::log2(double(n)); double entropy1 = 0; @@ -774,9 +787,10 @@ class SingleDrop { const double tol = std::max(1.0, info) * SQRT_EPS; const SplitCounts sc{n, n_child, overlap_.data(), n_split_, row_.data(), row_in_.data(), row_has_.data(), - col_.data(), col_in_.data(), col_has_.data()}; + col_.data(), col_in_.data(), col_has_.data(), + n_drop == 2}; pp_ = pp; - const ColPrices prices{&SingleDrop::Prices, this}; + const ColPrices prices{&LeafDrop::Prices, this}; const double mci = mutual_clustering_bounded_counts( sc, (info - cutoff) / 2 - tol, pp ? &prices : nullptr); const double ret = info - mci - mci; @@ -807,11 +821,12 @@ class SingleDrop { in1_[r] = Count(&bits1_[size_t(r) * n_word_]); in2_[r] = Count(&bits2_[size_t(r) * n_word_]); } - // Overlaps between parent splits; and, for each leaf, whether removing - // it makes two splits identical, which happens if they differ in that - // leaf alone + // Overlaps between parent splits; and the leaves whose removal makes two + // splits identical: those in which they differ, if they differ in one + // leaf (shared_) or two (shared_pair_) overlap_.resize(size_t(n_split_) * n_split_); shared_.assign(n_tip, 0); + shared_pair_.clear(); const int spare = n_word_ * 64 - n_tip; const uint64_t last_mask = spare ? ~uint64_t(0) >> spare : ~uint64_t(0); for (int i = 0; i != n_split_; ++i) { @@ -827,19 +842,29 @@ class SingleDrop { if (a_xor_b == 0 || a_xor_b == n_tip) { return false; // The trees share a split } - if (a_xor_b == 1 || a_xor_b == n_tip - 1) { - const uint64_t flip = a_xor_b == 1 ? 0 : ~uint64_t(0); - for (int w = 0; w != n_word_; ++w) { - const uint64_t diff = (a[w] ^ b[w] ^ flip) & - (w == n_word_ - 1 ? last_mask : ~uint64_t(0)); - if (diff) { - shared_[w * 64 + __builtin_ctzll(diff)] = 1; - break; - } + const int n_diff = std::min(a_xor_b, n_tip - a_xor_b); + if (n_diff > 2) { + continue; + } + // The leaves in which the splits differ + const uint64_t flip = a_xor_b == n_diff ? 0 : ~uint64_t(0); + int diff_leaf[2]; + int n_found = 0; + for (int w = 0; w != n_word_ && n_found != n_diff; ++w) { + uint64_t diff = (a[w] ^ b[w] ^ flip) & + (w == n_word_ - 1 ? last_mask : ~uint64_t(0)); + for (; diff && n_found != n_diff; diff &= diff - 1) { + diff_leaf[n_found++] = w * 64 + __builtin_ctzll(diff); } } + if (n_diff == 1) { + shared_[diff_leaf[0]] = 1; + } else { + shared_pair_.push_back(PairKey(diff_leaf[0], diff_leaf[1])); + } } } + std::sort(shared_pair_.begin(), shared_pair_.end()); return true; } @@ -883,6 +908,11 @@ class SingleDrop { } }; + static uint64_t PairKey(const int i, const int j) { + return i < j ? (uint64_t(i) << 32) | uint64_t(j) : + (uint64_t(j) << 32) | uint64_t(i); + } + int Count(const uint64_t *x) const { int ret = 0; for (int w = 0; w != n_word_; ++w) { @@ -891,31 +921,34 @@ class SingleDrop { return ret; } - // The split that vanishes from `ledger`'s tree when `tip` is removed; - // `merge`, if given, receives the split that it then duplicates, or -1. - // `tip`'s parent is not the root. - static int Vanishing(const Ledger &ledger, const std::vector &row_of, - const int tip, int *merge) { + // Removes `tip` from `ledger`'s tree (a parent tree, perhaps already less + // a leaf), returning the parent split that vanishes from it; `merge`, if + // given, receives the parent split that the vanishing split then + // duplicates, or -1. The root's children must be leaf 1 and a node with + // more than three leaves. + static int Remove(Ledger &ledger, const std::vector &row_of, + const int tip, int *merge) { const int parent = ledger.parent[tip]; const int sibling = ledger.sibling[tip]; - if (row_of[parent] >= 0) { + int ret; + if (ledger.parent[parent] != ledger.parent[1]) { if (merge) *merge = row_of[sibling]; - return row_of[parent]; + ret = row_of[parent]; + } else { + // `parent` is the root's child, whose split omits leaf 1 alone; the + // sibling's split will omit only leaf 1 too + if (merge) *merge = -1; + ret = row_of[sibling]; } - // `parent` is the root's child, whose split omits leaf 1 alone; the - // sibling's split will omit only leaf 1 too - if (merge) *merge = -1; - return row_of[sibling]; + ASSERT(ret >= 0); + ledger.RemoveTip(tip); + return ret; } - // Whether reduce_trees() would remove further leaves once `tip` is gone: - // the tests of its cherry, root and chain rules, applied to the parent - // trees with `tip` removed. - bool Reduces(const int tip) { - xc_ = x_; - yc_ = y_; - xc_.RemoveTip(tip); - yc_.RemoveTip(tip); + // Whether reduce_trees() would remove further leaves from the prepared + // trees: the tests of its cherry, root and chain rules, applied to the + // parent trees with the candidate's leaves removed. + bool Reduces() const { const int n_tip = n_tip_; const auto IsTip = [n_tip](const int i) { return i <= n_tip; }; const std::vector &xs = xc_.sibling, &xp = xc_.parent, @@ -956,20 +989,62 @@ class SingleDrop { return false; } + // Child splits of the prepared candidate, which removes leaf `i` and, if + // `j` >= 0, leaf `j`. + void Fill(const int i, const int j) { + const int n_child = n_split_ - (j < 0 ? 1 : 2); + row_.resize(n_child); + row_in_.resize(n_child); + row_has_.resize(n_child); + col_.resize(n_child); + col_in_.resize(n_child); + col_has_.resize(n_child); + const int word_i = i >> 6; + const uint64_t bit_i = uint64_t(1) << (i & 63); + const int word_j = j < 0 ? 0 : j >> 6; + const uint64_t bit_j = j < 0 ? 0 : uint64_t(1) << (j & 63); + const auto Has = [&](const uint64_t *split) { + return int16_t(((split[word_i] & bit_i) != 0) | + (((split[word_j] & bit_j) != 0) << 1)); + }; + for (int r = 0, k = 0; r != n_split_; ++r) { + if (r == drop1_[0] || r == drop1_[1]) continue; + const int16_t has = Has(&bits1_[size_t(r) * n_word_]); + row_[k] = r; + row_has_[k] = has; + row_in_[k] = in1_[r] - (has & 1) - (has >> 1); + ++k; + } + for (int c = 0, k = 0; c != n_split_; ++c) { + if (c == drop2_[0] || c == drop2_[1]) continue; + const int16_t has = Has(&bits2_[size_t(c) * n_word_]); + col_[k] = c; + col_has_[k] = has; + col_in_[k] = in2_[c] - (has & 1) - (has >> 1); + ++k; + } + } + // Parent prices for the prepared candidate's splits of tree 2. A split - // that two parent splits restrict to takes the lesser price. + // that several parent splits restrict to takes the least of their prices. static const price_t *Prices(void *context) { - SingleDrop &sd = *static_cast(context); - const std::vector &parent = sd.pp_->Price(); + LeafDrop &ld = *static_cast(context); + const std::vector &parent = ld.pp_->Price(); if (parent.empty()) return nullptr; // #nocov - const int n_child = sd.n_split_ - 1; - sd.price_.resize(n_child); + // The split into which the first vanishing split merges may itself + // vanish with the second leaf, merging in turn into the second's target + const int merge0 = ld.merge2_[0] >= 0 && ld.merge2_[0] == ld.drop2_[1] ? + ld.merge2_[1] : ld.merge2_[0]; + const int n_child = int(ld.col_.size()); + ld.price_.resize(n_child); for (int k = 0; k != n_child; ++k) { - const int c = sd.col_[k]; - sd.price_[k] = c == sd.merge2_ ? - std::min(parent[c], parent[sd.drop2_]) : parent[c]; + const int c = ld.col_[k]; + price_t p = parent[c]; + if (c == merge0) p = std::min(p, parent[ld.drop2_[0]]); + if (c == ld.merge2_[1]) p = std::min(p, parent[ld.drop2_[1]]); + ld.price_[k] = p; } - return sd.price_.data(); + return ld.price_.data(); } const int n_tip_; @@ -978,15 +1053,20 @@ class SingleDrop { const std::vector &bits1_; const std::vector &bits2_; bool ready_ = false; - Ledger x_, y_, xc_, yc_; + Ledger x_, y_, xi_, yi_, xc_, yc_; std::vector row_of1_, row_of2_, in1_, in2_; std::vector overlap_; std::vector shared_; - // The prepared candidate + std::vector shared_pair_; // Sorted PairKey()s + // Trees less leaf base_leaf_ + int base_leaf_ = -1; + int base_drop1_ = -1, base_drop2_ = -1, base_merge2_ = -1; + // The prepared candidate: the parent split that vanishes from each tree + // as each leaf is removed (-1 if only one leaf is removed), and the split + // of tree 2 that it then duplicates + int drop1_[2] = {-1, -1}, drop2_[2] = {-1, -1}, merge2_[2] = {-1, -1}; std::vector row_, row_in_, col_, col_in_; std::vector row_has_, col_has_; - int drop2_ = -1; - int merge2_ = -1; const ParentPrices *pp_ = nullptr; std::vector price_; }; @@ -999,7 +1079,7 @@ class SingleDrop { // Candidates that cannot beat the best so far may receive a lower bound // rather than their exact score. // `incremental = FALSE` scores every candidate by rebuilding its trees -// (SingleDrop is not used), for testing. +// (LeafDrop is not used), for testing. // [[Rcpp::export]] List spr_rogue_scores(const List &tree1, const List &tree2, const int depth, const bool incremental = true) { @@ -1032,16 +1112,16 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth, pp.n_split = full2.bits.nrow(); } const ParentPrices *pp_ptr = pp.n_split ? &pp : nullptr; - SingleDrop single(tree1["edge"], tree2["edge"], incremental ? n_tip : 0, - pp.n_word, bits1, pp.bits, node1, node2); + LeafDrop drop(tree1["edge"], tree2["edge"], incremental ? n_tip : 0, + pp.n_word, bits1, pp.bits, node1, node2); NumericVector scores(n_tip); double best = R_PosInf; bool all_finite = true; LogicalVector keep(n_tip, true); for (int i = 0; i != n_tip; ++i) { Rcpp::checkUserInterrupt(); - if (single.Prepare(i)) { - scores[i] = single.Score(best, pp_ptr); + if (drop.Prepare(i)) { + scores[i] = drop.Score(best, pp_ptr); } else { keep[i] = false; scores[i] = ScoreWithout(pre1, pre2, labels, keep, best, pp_ptr); @@ -1065,10 +1145,14 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth, Rcpp::checkUserInterrupt(); keep[i] = false; for (int j = i + 1; j != n_tip; ++j) { - keep[j] = false; - pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best, - pp_ptr); - keep[j] = true; + if (drop.Prepare(i, j)) { + pair_scores[pair] = drop.Score(best, pp_ptr); + } else { + keep[j] = false; + pair_scores[pair] = ScoreWithout(pre1, pre2, labels, keep, best, + pp_ptr); + keep[j] = true; + } if (!R_FINITE(pair_scores[pair])) { return List::create(Named("scores") = scores, _["pair_scores"] = pair_scores); diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp index e937e5e4..a3460d29 100644 --- a/src/tree_distances.cpp +++ b/src/tree_distances.cpp @@ -941,11 +941,11 @@ double mutual_clustering_bounded_score(const RawMatrix &x, const RawMatrix &y, } // As mutual_clustering_bounded_score(), for trees obtained by removing one -// leaf from a parent pair whose splits overlap by `sc.overlap`, when the -// resulting trees share no split. See split_counts.h. -double mutual_clustering_bounded_counts(const SplitCounts &sc, - const double min_score, - const ColPrices *prices) { +// or two leaves from a parent pair whose splits overlap by `sc.overlap`, when +// the resulting trees share no split. See split_counts.h. +template +double bounded_counts_impl(const SplitCounts &sc, const double min_score, + const ColPrices *prices) { const split_int n_tips = sc.n_tip; const split_int dim = sc.n_split; constexpr cost max_score = BIG; @@ -981,7 +981,9 @@ double mutual_clustering_bounded_counts(const SplitCounts &sc, cost *out = score.row(i); cost row_min = max_score; for (split_int j = 0; j < dim; ++j) { - const split_int a_and_b = overlap[sc.col[j]] - (row_has & sc.col_has[j]); + const int16_t both = row_has & sc.col_has[j]; + const split_int a_and_b = overlap[sc.col[j]] - + (Pair ? (both & 1) + (both >> 1) : both); const split_int nb = sc.col_in[j]; const split_int a_and_B = na - a_and_b; const split_int A_and_b = nb - a_and_b; @@ -1037,6 +1039,13 @@ double mutual_clustering_bounded_counts(const SplitCounts &sc, return static_cast(max_score * dim - lap_cost) * over_max_score; } +double mutual_clustering_bounded_counts(const SplitCounts &sc, + const double min_score, + const ColPrices *prices) { + return sc.pair ? bounded_counts_impl(sc, min_score, prices) : + bounded_counts_impl(sc, min_score, prices); +} + // Optimal column duals of the mutual clustering LAP between the splits `x` // and `y` (identical splits included), one per split in `y`, in cost units, // shifted so that the least is zero and capped at BIG: prices for diff --git a/tests/testthat/test-tree_distance_spr.R b/tests/testthat/test-tree_distance_spr.R index ab4c083f..c2112bec 100644 --- a/tests/testthat/test-tree_distance_spr.R +++ b/tests/testthat/test-tree_distance_spr.R @@ -326,7 +326,7 @@ test_that("C++ Rogue scorer matches R implementation", { } }) -test_that("Single-leaf candidates scored from parent splits", { +test_that("Leaf and leaf-pair candidates scored from parent splits", { library("TreeTools", quietly = TRUE) # Scores must agree with those from rebuilding the reduced trees, to within # numerical error, with the same best candidates; candidates that cannot From 5f5498ca4bce4d469632dc5ad3ef14bb5676fea7 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Fri, 23 Oct 2026 09:59:10 +0100 Subject: [PATCH 2/4] Record leaf-pair Rogue results in NEWS and Optimizations.md Co-Authored-By: Claude Opus 5.5 --- .AGENTS/memory/Optimizations.md | 36 +++++++++++++++++++++++++++++++-- NEWS.md | 4 ++++ 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/.AGENTS/memory/Optimizations.md b/.AGENTS/memory/Optimizations.md index 2a553898..e8a02d0a 100644 --- a/.AGENTS/memory/Optimizations.md +++ b/.AGENTS/memory/Optimizations.md @@ -1018,5 +1018,37 @@ pay; what pays is what the parent overlaps let a candidate *skip*. 600k candidates). Rejected: storing parent cell information and updating it by an O(1) delta (5.3 cycles/cell vs 8.5) - not worth a d^2 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. +- Next: depth-2 pair candidates (two leaves) by the same route (done, see + below); the unbounded MCI could use the Tier 1 form if tests tolerate ulp + changes. + +## Rogue SPR: leaf-pair candidates from parent overlaps (2026-10-23) +`SingleDrop` became `LeafDrop`, with `Prepare(i, j)` for depth-2 candidates. +- Two removals are replayed on copies of the parent ledgers (leaf i once per + outer loop, then j per candidate), recording the vanishing split per tree + at each step. Covers the cherry, sib(i) = par(j), mirror and root-child + cases without enumerating them; the one trap is that the root-child test + must be structural (`parent[parent[tip]] == parent[1]`), not + `row_of[parent] < 0`: after removing i, sib(i) can become the root's + child while still owning a row. +- Overlap correction is popcount(row_has & col_has) on 2-bit masks; + `mutual_clustering_bounded_counts()` dispatches on `SplitCounts::pair` to + a template so the depth-1 kernel is unchanged. +- Shared-split test: parent pairs whose polarized XOR is two leaves are + recorded as sorted pair keys; a pair also falls back if either leaf alone + makes a shared split. Price merge chains: the first vanishing split's + merge target can itself vanish second, then it merges into the second's + target (prices only affect speed). +- Fallback rate on 295 pairs (9-40 leaves, half near): 75% of incremental + pair candidates scored from overlaps; fallbacks are leaf 1 or a disabled + tree (14%), a single-leaf shared split (7%), a pair-shared split (1.6%), + reduction (1%). +- Min of 3 interleaved runs vs main, depth 2: rand30 3.86 -> 1.65 s (-57%), + mid80 5.47 -> 3.10 s (-43%), near80 0.39 -> 0.21 s (-46%). Depth 1 + unchanged within noise (rand80 1.03-1.09 vs 1.07-1.11 s over 4 rounds). + SPRDist sums identical on all workloads. +- Correctness: 295 tree pairs, 70,142 pair candidates, with ASSERTs on: + argmin sets of single and pair score vectors identical to + `incremental = FALSE`; max difference among near-best scores 7e-15. +- Remaining at depth 2: candidates involving leaf 1 (~2/n) could be scored + by rerooting; ScoreWithout() for shared-split pairs. diff --git a/NEWS.md b/NEWS.md index df06e482..622063d6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -51,6 +51,10 @@ the last few bits, which can break exact ties between candidates differently. +- With `options(sprDepth = 2)`, `SPRDist(method = "Rogue")` likewise scores + most pairs of leaves from the parent trees' splits: 45–55% faster at 30–80 + leaves. + - Improve performance when comparing many similar trees with `ClusteringInfoDistance()`, `MatchingSplitDistance()`, `InfoRobinsonFoulds()` and `JaccardRobinsonFoulds()`. From 67ee0fb0dc9c881070ea3277f188ab46971b2855 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Thu, 29 Oct 2026 06:55:07 +0000 Subject: [PATCH 3/4] CI: unblock R 4.1 and benchmark legs The 7-leaf NNIDist() bounds test (#222) called NNIDist() on all 945 trees: 446,040 pairs through CompareAll()'s R loop, c. 80 s locally and over an hour on the R 4.1 leg, which then hit its 60-minute limit. The test now calls cpp_nni_distance() on pre-ordered edges (c. 7 s); the 6-leaf case still checks that NNIDist() agrees. benchmark and rev-dep-check ran on ubuntu-latest, now 26.04, where setup-r apt-installs ~175 packages and timed out at 15 minutes; their RSPM serves noble binaries anyway, so pin ubuntu-24.04. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/R-CMD-check.yml | 6 ++++-- tests/testthat/test-tree_distance_nni.R | 23 ++++++++++++++++++----- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml index bbdc2526..03024173 100644 --- a/.github/workflows/R-CMD-check.yml +++ b/.github/workflows/R-CMD-check.yml @@ -338,7 +338,9 @@ jobs: !cancelled() && needs.sense-check.result == 'success' && contains(fromJSON('["success", "skipped"]'), needs.codemeta.result) && (github.event_name == 'pull_request' || github.event_name == 'workflow_dispatch') - runs-on: ubuntu-latest + # Pinned: ubuntu-latest is now 26.04, on which setup-r apt-installs ~175 + # packages and has timed out; RSPM below serves noble (24.04) binaries. + runs-on: ubuntu-24.04 # A hung setup or a wedged benchmark otherwise runs to the 6-hour default, # holding a runner the whole time. timeout-minutes: 45 @@ -672,7 +674,7 @@ jobs: # and codemeta is skipped outside pull requests, so without this the # backstop never ran on push, schedule or workflow_dispatch. if: ${{ !cancelled() && needs.core.result == 'success' }} - runs-on: ubuntu-latest + runs-on: ubuntu-24.04 # As benchmark timeout-minutes: 120 # revdepchecks are genuinely long name: revdepcheck (${{ matrix.config.deps }}) diff --git a/tests/testthat/test-tree_distance_nni.R b/tests/testthat/test-tree_distance_nni.R index 09ec2f7d..b33c69d2 100644 --- a/tests/testthat/test-tree_distance_nni.R +++ b/tests/testthat/test-tree_distance_nni.R @@ -186,13 +186,26 @@ test_that("NNIDist() lower bounds are valid", { testthat::expect_equal(max(exact), NNIDiameter(nTip)[, "exact"], ignore_attr = TRUE) - bounds <- lapply(NNIDist(trees), as.matrix) + # Bounds for each pair, in the order of exact[lower.tri(exact)]. + # Calling the C++ kernel directly avoids NNIDist()'s per-pair R overhead, + # which makes 7 leaves (446,040 pairs) take minutes, and an hour on R 4.1; + # tip labels are already in a common order. + pairs <- which(lower.tri(exact), arr.ind = TRUE) + edges <- lapply(trees, function(tr) Postorder(tr[["edge"]])) + bounds <- vapply(seq_len(nrow(pairs)), function(k) { + cpp_nni_distance(edges[[pairs[k, 1]]], edges[[pairs[k, 2]]], nTip) + }, integer(6)) + if (nTip == 6) { + viaR <- vapply(NNIDist(trees), function(d) as.matrix(d)[pairs], + double(nrow(pairs))) + testthat::expect_equal(t(bounds), viaR, ignore_attr = TRUE) + testthat::expect_equal(colnames(viaR), rownames(bounds)) + } exact <- exact[lower.tri(exact)] - .Bound <- function(field) bounds[[field]][lower.tri(bounds[[field]])] - testthat::expect_true(all(.Bound("lower") <= exact)) + testthat::expect_true(all(bounds["lower", ] <= exact)) for (upper in c("tight_upper", "best_upper", "loose_upper", "fack_upper", "li_upper")) { - testthat::expect_true(all(.Bound(upper) >= exact)) + testthat::expect_true(all(bounds[upper, ] >= exact)) } } @@ -201,7 +214,7 @@ test_that("NNIDist() bounds bracket exact distances: 6 leaves", { }) test_that("NNIDist() bounds bracket exact distances: 7 leaves", { - skip_on_cran() # c. 70 s + skip_on_cran() # c. 7 s .ExpectBoundsBracketExact(7) }) From bb046d216b7ae261747c6b46d86727035a2e4b91 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Thu, 29 Oct 2026 07:16:04 +0000 Subject: [PATCH 4/4] CompareAll(): index a plain list, not a pair-length multiPhylo On R 4.1, ape's [[.multiPhylo (class(x) <- NULL) copies the whole list on each access. CompareAll() built two multiPhylo objects with one element per pair and indexed each once per pair, so its time grew with the square of the number of pairs: 446,040 pairs of 7-leaf trees took over 50 minutes on the R 4.1 CI leg, against c. 80 s on R 4.6, where the access does not copy. It now indexes a plain list of the trees by pair indices. Co-Authored-By: Claude Opus 5.5 --- NEWS.md | 4 +++ R/tree_distance_utilities.R | 21 ++++++++------ tests/testthat/_problems/test-parallel-19.R | 29 ++++++++++++++++++++ tests/testthat/testthat-problems.rds | Bin 0 -> 4400 bytes 4 files changed, 45 insertions(+), 9 deletions(-) create mode 100644 tests/testthat/_problems/test-parallel-19.R create mode 100644 tests/testthat/testthat-problems.rds diff --git a/NEWS.md b/NEWS.md index df06e482..718c7c07 100644 --- a/NEWS.md +++ b/NEWS.md @@ -35,6 +35,10 @@ ## Performance +- `CompareAll()` no longer takes time quadratic in the number of pairs on + older versions of R (e.g. 4.1), where each element access to a + `multiPhylo` object copied the whole list. + - 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/tree_distance_utilities.R b/R/tree_distance_utilities.R index 518b7bba..4cb4976d 100644 --- a/R/tree_distance_utilities.R +++ b/R/tree_distance_utilities.R @@ -535,28 +535,31 @@ Ntropy <- function(...) { CompareAll <- function(x, Func, FUN.VALUE = Func(x[[1]], x[[1]], ...), ...) { nTree <- length(x) + labels <- names(x) countUp <- seq_len(nTree - 1) - i <- x[rep(countUp, rev(countUp))] - j <- x[unlist(sapply(countUp, function(n) n + seq_len(nTree - n)))] - + i <- rep(countUp, rev(countUp)) + j <- unlist(lapply(countUp, function(n) n + seq_len(nTree - n))) + # A plain list: on R 4.1, each `[[` on a multiPhylo copies the whole list + x <- lapply(seq_len(nTree), function(k) x[[k]]) + cluster <- getOption("TreeDist-cluster") ret <- if (is.null(cluster)) { cli_progress_bar("Comparing", total = length(i)) vapply(seq_along(i), function(k) { cli_progress_update(1, .envir = parent.frame(2)) - Func(i[[k]], j[[k]], ...) + Func(x[[i[k]]], x[[j[k]]], ...) }, FUN.VALUE) } else { - do.call("cbind", - parLapply(cluster, seq_along(i), - function(k, Func) Func(i[[k]], j[[k]]), + do.call("cbind", + parLapply(cluster, seq_along(i), + function(k, Func) Func(x[[i[k]]], x[[j[k]]]), Func = Func)) } - + .WrapReturn <- function(dists) { structure(dists, Size = nTree, - Labels = names(x), + Labels = labels, Diag = FALSE, Upper = TRUE, class = "dist") diff --git a/tests/testthat/_problems/test-parallel-19.R b/tests/testthat/_problems/test-parallel-19.R new file mode 100644 index 00000000..37e11cd1 --- /dev/null +++ b/tests/testthat/_problems/test-parallel-19.R @@ -0,0 +1,29 @@ +# Extracted from test-parallel.R:19 + +# setup ------------------------------------------------------------------------ +library(testthat) +test_env <- simulate_test_env(package = "TreeDist", path = "..") +attach(test_env, warn.conflicts = FALSE) + +# test ------------------------------------------------------------------------- +library("TreeTools") +trees <- as.phylo(0:20, 20) +suppressMessages({ + expect_equal(options("TreeDist-cluster"), StartParallel(2)) + cl <- getOption("TreeDist-cluster") + + parallel <- ClusteringInfoDistance(trees) + SetParallel(NULL) + expect_null(GetParallel()) + expect_equal(ClusteringInfoDistance(trees), parallel) + + expect_equal(options("TreeDist-cluster"), SetParallel(cl)) + expect_equal(cl, getOption("TreeDist-cluster")) + expect_equal(cl, GetParallel()) + parallel <- CompareAll(trees, NNIDist) + SetParallel(NULL) + expect_equal(CompareAll(trees, NNIDist), parallel) + expect_false(StopParallel()) + SetParallel(cl) + expect_true(StopParallel()) + }) diff --git a/tests/testthat/testthat-problems.rds b/tests/testthat/testthat-problems.rds new file mode 100644 index 0000000000000000000000000000000000000000..47bbda638a15b888ebfcce0bb55321c0e999cdb9 GIT binary patch literal 4400 zcmW-hcRUpSAIHO2LS)OxRm#jxwo{TlN>;{|%#h7x+)<&-?Cp#=BhDTPb8ccFEx56~>85=a6dMmUYlVDD z_a0~z-v0i+R?40w)f6HuT4)yjYYMJ$vg0zKQ{T_TPph`kP0F`yZ{U^1Ls9Wc>mu8YGcxU~4CRr7eI`vTLh(iLrm!WyHj8 zW8y?I(A;C@-j>6R5KAt>_*+G-L+{;|yzJeAbG4-g>n>Qyrkk0$O^q@Yd!4rm^Tjfl z3lr~kFiU3x@(XhtGI8J?-dM$w^jQm8n5NbWAM#eAZpV)0wWbBE>{MKU)>Fi{WV!UF z9Ibh^jA^!&5Z#XZw;;cGL0KauGa?dg`oE18H(WP_2M>w@xIEN5roWQ!`WqFS-=f^B zDEwW=XGhnP$aew9V@&WcBvfAHAFl5#UUtnKeaXQ*XUT=skHyO}qwHhxl2tq5E?{57 zxo}r7$dI5dfNHyO>;xpcVC&KrN>z{M#>SpmeleF!EEgZH6A4sLZ@C}6h1=(V1eqHy z>X)}_2*}jpO1l(M&HV|71IJzEjv)bi!P&neW;S+`oOcv-5>X#qr*?HJ^GAgqEl^;F<*=7SHsa{-oFE3-mq;5G`_!_o^A0n% z2s(b1aj#8>8CArt!;HkgYrj`_UlKvkuVJavC}eo7YjULa%CvGczgpwgXkXHndXRtF z=bdq>mnK+8`@W=g+I^`*jCfr?un_;ZKj6T^VcVcimtP~89U*U*h2uvFHMT!!>|oz50hkMs>>uwY%%p2gMl7`>-{YPa^Y z{?g08Lo%g=6ND=pzuI#M{_OD9Of|9dv-MO(A4SzU@3Mgtm)}?qnV#BY_m4h>uDvhF z91b2SPz(HAgV^&^gt~;32Fj*uc_9=uQhGoQ{o9G+(@k2AIPUs_g1lr`Ovq}D(kGYi1GaUu*`GyX~f`@jv0mG z-=4OkX=M~Kc9-Xe_tevJId_+Z5!V@TyOZ9smwq6B-3h!6_^~RvZ_u3S;Gv47pDx*$St=OQc)d%{13@18
sWEniAFrX_)^aep>M7u*;v?55sjIN<^_Vgjc{ z*G#=?Y_AkRx9@4RV)FO7pO zwDPs{M~2FkZ33+nL*tb8-xlaw-x(->K^*XG?|wd^v`}Uy!G3ne1VoZF^3NnMu^)20 zzxLSTLTGX#`z$vm-I$cYCiiE|r|L3qpRQDu*@1JjC%k{X-Pou4GO@$H$iS&yEUs2f zaIr3z%MaGX-bYk$pG0VPv#?hCgRj0`ySsetz5cd|Ks1~Oyes5iUx(s5acc&7sZZhH$SY&fh6?GsWqHO{q6O{mffqKS1eY~h`S6wec$<4im`}{Rfv=!T77~X z=tx_j^8BBnc{So#%nxlO$OyMnD_iu{vKRYuV5L#@4D8i!r;PuBSSQQ7(<#9w=plmJ z`3&@dj?lxo@2V^2CUok&;6`-nhK+*kw(hi0ew0Z}5#e4N?@EmT{~3u5j`)CLJ`r>S z?YG$Tgf&80dsLLnD1iyrA<6TD`;+WjxzI7D5Q+V_^Ir%;m*HnkmZ|e>7w0{BTIH~z z&!~2<)xX#XI+1qhHg#VyEKQ=i{EitJi?)vxB1M;iB6nhXFajf&Iqo}sBjXwVq!7dU z2dL*pBUv<0<}v%vmTPh;_R;Lc$oKw@RE6Rm7R@PO_rBMMHo>_AD2E5fipFL5VW;UE zanc(Wsz%P4e9gRBe-%{Z*=fjfxE&>ho~;ZC^RK3va??Bl+)g7CEvbC2JCE1dDn!ZW zmNH-|(&S6>WM*qLC6W7#`@u8nZw%q+8MWpPvH>wp=*O7TCzT;t)RGX+WMKv)NrDTs zV8iPJ;^|3Tcs-~7iBZ2Km0QQ12iyp7KRv$;rZ$mawk|g}*o zLyQvT&a#FrIZ%iHAu8Q}5YNa~v<-TEJU_5+!*j59%A|n4*3;Q$d3K*7LLYW{B{_3RH;XMw%NWGEZtF`nNYV;7O z4FwDC=?608;U7^@=l3rTZVA9oBd@IM(FP3cp^eWzm*RU)D}m^fm`?Gf*Fw+;bkBuo z1Rzn{TX=gNObIfH=C)0lmRP}8sL%uuvxyy1oduzP>=%ANsy{|C?b8p23&4+2h%f0> zex4J&>nVRjPTnH7#Ox zE`%*W8+s}*DlIYpMRf(gCa5_D)mwz6pb?FyIG2JyNBk!9$;j-V)NCC@2{MHS3?gdI zVqt*lWrF5c8$oZM?RgbrAEN%f`2OqSXIdAbU0*g;x>cvBwGV#-kI7r6+t7nbxb}mS z6NXWHiM9JDC!(Vy8nDoWjh@)Hq6#vZdEF-iK?{e+{7WEzLnw9P)TH;wra06B7408S z^^M*8Zr`6y3C)3s63u$H=Ales{41T8;bf{V4RO2_8Si|0WJ~#cS)4?Dw2V$zod0&n zR;KMu>&3hdpnir6P9`(4Z_#L^(oJyh9(HZ!G_FVszz++t3u7pW)lH#P-~GS*J{-j@ zKgIeLAPuq1lM)0Mbo_>CzXH0Hoe`YA=w%y`FCZRakP18`8?SF4+C_(bnIfy7>tMNE z8X*vH2hHB?AYE@tgC)tt9Luf+5nE3sD6;GGi!=9lmZ+75|T^NAGeFf&DSqmM3XV$4n0J@y5|4uEO zFseVzH9>+_w0FQnTx!!vxQ-_1PaLr&F+H-~=6o`(Tqs}yp_@#>>;V1bF3!XUr2B8e zNHheY(i3f0Q{+i)bXV@Zl3tqBenU$d_^}d@4`43LM_);_B!ZB?7PEVwwx(otgiukGEe(<#M8 zMbqg;gyJR`e$Rc)YKXeVMkUdWbBhEWl9&oC+Y89R6!%m(u%brJ>k!A9Qf3?@&; ziyN9B!wYn&o$CY3w?i1NIkD{XpR5FG2mUfjh)a29pJ1>40IwA?{e7%K!|dWE5rCk8 zh>#|DUqWD3WGXo>^H)snt**o?S`h5-h3|xQ|2BvGY6m}3IZkELxTabdSQjv}V22x< z#;aNy)FS>VOkigGHiZwF*KD_a*Jh&jA+j3yINA^=%1K{UkY2uB!O0jH5XzsB21fmp+O6=8r zkI~kDK`guYxXprs$@Nr|AwUl>iB9lGs}$M_E{q7MG^rtWt`;wL*K2w{fHpFFUc>R!AIS^L z&&*q~*_u4vNNQK{k-YUy?4(bBP`Iurdj>r5`JVe{Uf4YSh|)pup6R%7OREHkC^Z=u znRP>&ZEtS(iKdeZKGU-Zo9+8?iR%eO%Qt zh!7HF^B-euRVmi!xZ)MsZcp&NT)R|Z`RjZ5{6Mdz7vob$!&{EQte7jDamqtC_blHm z@Ti}q+J-5*x_UnKgOy7aA*{_cnk)wXjQ4u0?>IJS8(e*<8s3!o{GMXT2U*h!_(Iyy z^ZTi<-I}R)Ux(ivB=~8E`%6U9^x(5G86D&8G$x*HjuGRM+{&@PQWTnXTN~q4K&mcI z3!e(ue{p%vC)N#4_u!d*EKZFxkiUiumL!3ckr8@OG|}MHo%fIm)LZ9L1f+B%97_2c zR(NtcMQZ07n|=Hqnv_)Zz60|+Yx(J-o1ylL#@r5e;&Xh=>1}c@3`W)2DOxKL(wOB4 zoqHG>G+mrlm+vi=rOlR|w5E|^dQ&_i?J4ASHD})27I*#YTV&CsnxP(<^Ro<0z$DW} zFPzn2HmjK}j>kGPBKjmeK3s_CT+`Bdn8VGvI$&ge`d?