From 310876108d76ce8a3217bc44be9b7aefbeca6bf9 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Mon, 19 Oct 2026 16:24:19 +0100 Subject: [PATCH 1/2] Score Rogue single-leaf candidates from parent split overlaps Most single-leaf candidates (96% on random pairs) need no further reduction and share no split once the leaf is removed. Their trees are then the parent trees less one split each, so the MCI cost matrix follows from the parent overlaps minus the leaf's membership (SingleDrop, mutual_clustering_bounded_counts()), skipping KeepAndReduce(), GetSplits() and FirstMatchingSplit(). Eligibility replays reduce_trees()' cherry, root and chain tests on ledgers with the leaf removed. The bounded fill also writes each cell as sum f(n_q) + row + column terms with f(k) = k lg k tabled, instead of four lookups and multiplications. Popcount was measured not to be the fill's cost. Min-of-9 vs main: 30 tips -56%, 80 -61%, 150 -60%; 15/30-move pairs -53%/-55%; depth 2 -3%. Argmin sets match the rebuild path on 13,833 score vectors; vs main only exact ties broken by rounding noise differ. Co-Authored-By: Claude Opus 5.5 --- .AGENTS/memory/Optimizations.md | 33 +++ NEWS.md | 7 + R/RcppExports.R | 4 +- dev/rogue_mci_bound_findings.md | 4 + src/RcppExports.cpp | 10 +- src/split_counts.h | 27 ++ src/spr.cpp | 327 +++++++++++++++++++++++- src/tree_distances.cpp | 139 +++++++++- tests/testthat/test-tree_distance_spr.R | 73 +++++- 9 files changed, 605 insertions(+), 19 deletions(-) create mode 100644 src/split_counts.h diff --git a/.AGENTS/memory/Optimizations.md b/.AGENTS/memory/Optimizations.md index aeaa5664..2a553898 100644 --- a/.AGENTS/memory/Optimizations.md +++ b/.AGENTS/memory/Optimizations.md @@ -987,3 +987,36 @@ row-filling skips ≤1.5% of fill. Next lever: build child cost matrices incrementally from the parent (97.6% of candidates drop one leaf); the fill is now 37–55% of runtime. Parent prices are computed lazily (first candidate past the row bound). + +## Rogue SPR: single-leaf candidates from parent overlaps (2026-10-19) +Measured first (rdtsc-instrumented release build, workloads in +`dev/rogue_bound/workloads.rds`): fill 35% of runtime at 80 tips, 54% at 150; +KeepAndReduce + GetSplits + FirstMatchingSplit 46% / 31%. **Popcount is not +the fill's cost**: in a standalone kernel (1-3 64-bit bins), replacing the +popcount with a lookup of parent counts saves 1 of ~14.5 cycles/cell. The +cost is the per-cell log arithmetic. So "O(1) quadrant counts" alone does not +pay; what pays is what the parent overlaps let a candidate *skip*. +- Tier 1 (all bounded MCI calls): cell information as + `sum_q f(n_q) + row_term + col_term`, `f(k) = k lg2 k` tabled + (`KLg2K()`), row/column terms hoisted. Kernel 14.5 -> 8.5 cycles/cell; + costs agree within 1 unit. Alone: -6..-16% (80/150 tips). +- Tier 2 (`SingleDrop` in `spr.cpp`, `mutual_clustering_bounded_counts()`): + for leaf l != 1 whose removal neither creates a reducible structure + (cherry/root/chain tests of `reduce_trees()` replayed on ledgers with l + removed) nor makes a pair of splits identical (precomputed: parent pairs + differing in l alone), the child matrix comes from parent overlaps minus l's + membership; the vanishing split per tree is par(l)'s row, or sib(l)'s if + par(l) is the root's child. 96% of candidates qualify (84% on 15-move + pairs). Predicate agreed with KeepAndReduce on all 589,686 candidates. +- Min of 9 interleaved runs vs main: rand30 -56%, rand80 -61%, rand150 -60%, + mid80 -53%, mid150 -55%; depth 2 only -3% (pair candidates still rebuild). +- Correctness: 13,833 score vectors (2.9 M candidates): argmin sets identical + to the non-incremental path in the same build; vs main, 5 differ, all exact + ties (within 1e-12) that main broke by rounding noise. 1/410 SPRDist + results changed (an 8-leaf tie). Max score difference 4e-13 bits. +- Dropped from the counts path: the initial entropy bound (never fired in + 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. diff --git a/NEWS.md b/NEWS.md index 4db05c58..c641aeff 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,6 +18,13 @@ problems and warm-starting the rest: c. 50% faster at 80 leaves and 65% at 150 leaves on unrelated trees. +- `SPRDist(method = "Rogue")` scores most single-leaf candidates from the + overlaps between the parent trees' splits, without rebuilding, reducing or + splitting the smaller trees, and fills cost matrices with fewer operations: + a further 55–60% faster at 30–150 leaves. Scores can differ from before in + the last few bits, which can break exact ties between candidates + differently. + - Improve performance when comparing many similar trees with `ClusteringInfoDistance()`, `MatchingSplitDistance()`, `InfoRobinsonFoulds()` and `JaccardRobinsonFoulds()`. diff --git a/R/RcppExports.R b/R/RcppExports.R index 5e1a70fa..0ebe6ad4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -208,8 +208,8 @@ keep_and_reduce <- function(tree1, tree2, keep) { .Call(`_TreeDist_keep_and_reduce`, tree1, tree2, keep) } -spr_rogue_scores <- function(tree1, tree2, depth) { - .Call(`_TreeDist_spr_rogue_scores`, tree1, tree2, depth) +spr_rogue_scores <- function(tree1, tree2, depth, incremental = TRUE) { + .Call(`_TreeDist_spr_rogue_scores`, tree1, tree2, depth, incremental) } spr_bounded_cid <- function(tree1, tree2, cutoff) { diff --git a/dev/rogue_mci_bound_findings.md b/dev/rogue_mci_bound_findings.md index ef5450be..9ef843ed 100644 --- a/dev/rogue_mci_bound_findings.md +++ b/dev/rogue_mci_bound_findings.md @@ -72,3 +72,7 @@ The LAPs that remain run 2.0×, 3.4× and 4.1× faster at 30, 80 and 150 tips. T - 97.6% of candidates drop a single leaf. A child cell's quadrant counts are the parent cell's counts minus that leaf's quadrant, so each cell costs O(1). - This would also skip KeepAndReduce and GetSplits (about 30% of runtime) for those candidates. - It is the largest remaining lever. + +## Follow-up: single-leaf candidates (2026-10-19) + +Done; see `.AGENTS/memory/Optimizations.md`, "Rogue SPR: single-leaf candidates from parent overlaps". Popcount turned out not to be the fill's cost (parent counts save ~1 of ~14.5 cycles per cell); the gains come from a cheaper cell formula and from skipping KeepAndReduce, GetSplits and FirstMatchingSplit for the 96% of candidates whose trees need no further reduction: 55–61% faster at 30–150 leaves. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 81edc3dd..973f1701 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -546,15 +546,16 @@ BEGIN_RCPP END_RCPP } // spr_rogue_scores -List spr_rogue_scores(const List& tree1, const List& tree2, const int depth); -RcppExport SEXP _TreeDist_spr_rogue_scores(SEXP tree1SEXP, SEXP tree2SEXP, SEXP depthSEXP) { +List spr_rogue_scores(const List& tree1, const List& tree2, const int depth, const bool incremental); +RcppExport SEXP _TreeDist_spr_rogue_scores(SEXP tree1SEXP, SEXP tree2SEXP, SEXP depthSEXP, SEXP incrementalSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const List& >::type tree1(tree1SEXP); Rcpp::traits::input_parameter< const List& >::type tree2(tree2SEXP); Rcpp::traits::input_parameter< const int >::type depth(depthSEXP); - rcpp_result_gen = Rcpp::wrap(spr_rogue_scores(tree1, tree2, depth)); + Rcpp::traits::input_parameter< const bool >::type incremental(incrementalSEXP); + rcpp_result_gen = Rcpp::wrap(spr_rogue_scores(tree1, tree2, depth, incremental)); return rcpp_result_gen; END_RCPP } @@ -828,7 +829,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } - // cpp_popcount_force_sw bool cpp_popcount_force_sw(const bool force); RcppExport SEXP _TreeDist_cpp_popcount_force_sw(SEXP forceSEXP) { @@ -885,7 +885,7 @@ static const R_CallMethodDef CallEntries[] = { {"_TreeDist_confusion", (DL_FUNC) &_TreeDist_confusion, 2}, {"_TreeDist_keep_and_reroot", (DL_FUNC) &_TreeDist_keep_and_reroot, 3}, {"_TreeDist_keep_and_reduce", (DL_FUNC) &_TreeDist_keep_and_reduce, 3}, - {"_TreeDist_spr_rogue_scores", (DL_FUNC) &_TreeDist_spr_rogue_scores, 3}, + {"_TreeDist_spr_rogue_scores", (DL_FUNC) &_TreeDist_spr_rogue_scores, 4}, {"_TreeDist_spr_bounded_cid", (DL_FUNC) &_TreeDist_spr_bounded_cid, 3}, {"_TreeDist_spr_first_shared_split", (DL_FUNC) &_TreeDist_spr_first_shared_split, 2}, {"_TreeDist_spr_keep_and_reduce", (DL_FUNC) &_TreeDist_spr_keep_and_reduce, 3}, diff --git a/src/split_counts.h b/src/split_counts.h new file mode 100644 index 00000000..ecb72f0a --- /dev/null +++ b/src/split_counts.h @@ -0,0 +1,27 @@ +#ifndef TREEDIST_SPLIT_COUNTS_H_ +#define TREEDIST_SPLIT_COUNTS_H_ + +#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). +struct SplitCounts { + int n_tip; // Leaves in each child tree + int n_split; + const int16_t *overlap; + int stride; + const int *row; + const int *row_in; + const int16_t *row_has; + const int *col; + const int *col_in; + const int16_t *col_has; +}; + +#endif diff --git a/src/spr.cpp b/src/spr.cpp index 387164f7..3f932823 100644 --- a/src/spr.cpp +++ b/src/spr.cpp @@ -7,6 +7,7 @@ #include #include "reduce_tree.h" #include "col_prices.h" +#include "split_counts.h" #include // for log2 #include // for sort, lower_bound, equal #include @@ -266,6 +267,9 @@ extern double mutual_clustering_bounded_score(const RawMatrix &x, const int32 n_tip, const double min_score, const ColPrices *prices); +extern double mutual_clustering_bounded_counts(const SplitCounts &sc, + const double min_score, + const ColPrices *prices); extern std::vector mutual_clustering_col_prices( const RawMatrix &x, const RawMatrix &y, const int32 n_tips); @@ -329,10 +333,12 @@ bool MapLeaves(const ParentPrices &pp, const CharacterVector &labels, // its non-tip child. Edges must be in postorder. // If `leaf_index` is given, also writes each split to `indexed` as // `n_word` 64-bit words in which leaf t sets bit leaf_index[t]. +// If `row_node` is given, also writes the node (1-based) of each split. TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip, const std::vector *leaf_index = nullptr, const int n_word = 0, - std::vector *indexed = nullptr) { + std::vector *indexed = nullptr, + std::vector *row_node = nullptr) { const int n_edge = edge.nrow(); const int n_vertex = n_edge + 1; const int n_bin = (n_tip - 1) / SPLIT_BIN + 1; @@ -392,6 +398,12 @@ TreeSplits GetSplits(const IntegerMatrix &edge, const int n_tip, } } const int n_split = rows.size(); + if (row_node) { + row_node->resize(n_split); + for (int r = 0; r != n_split; ++r) { + (*row_node)[r] = rows[r] + 1; + } + } RawMatrix ret(n_split, n_bin); if (leaf_index) { indexed->resize(size_t(n_split) * n_word); @@ -679,6 +691,295 @@ 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. +// +// 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 { +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) + : n_tip_(n_tip), n_word_(n_word), n_split_(int(node1.size())), + bits1_(bits1), bits2_(bits2) { + ready_ = Init(edge1, edge2, node1, node2); + } + + // 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)) { + 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; + } + 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; + } + 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; + // As GetSplits() + const double n_lg2_n = n * std::log2(double(n)); + double entropy1 = 0; + double entropy2 = 0; + for (int k = 0; k != n_child; ++k) { + const int a = row_in_[k]; + const int b = col_in_[k]; + entropy1 += n_lg2_n - a * std::log2(double(a)) - + (n - a) * std::log2(double(n - a)); + entropy2 += n_lg2_n - b * std::log2(double(b)) - + (n - b) * std::log2(double(n - b)); + } + // As BoundedCID() + const double info = entropy1 / n + entropy2 / n; + 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()}; + pp_ = pp; + const ColPrices prices{&SingleDrop::Prices, this}; + const double mci = mutual_clustering_bounded_counts( + sc, (info - cutoff) / 2 - tol, pp ? &prices : nullptr); + const double ret = info - mci - mci; + return ret < tol ? 0 : ret; + } + +private: + // False if no candidate can be scored here. + bool Init(const IntegerMatrix &edge1, const IntegerMatrix &edge2, + const std::vector &node1, const std::vector &node2) { + const int n_tip = n_tip_; + if (n_tip < 8 || n_split_ != n_tip - 3 || int(node2.size()) != n_split_) { + return false; + } + const int n_vertex = 2 * n_tip; + x_.Build(edge1, n_vertex); + y_.Build(edge2, n_vertex); + row_of1_.assign(n_vertex, -1); + row_of2_.assign(n_vertex, -1); + in1_.resize(n_split_); + in2_.resize(n_split_); + for (int r = 0; r != n_split_; ++r) { + row_of1_[node1[r]] = r; + row_of2_[node2[r]] = r; + 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 + overlap_.resize(size_t(n_split_) * n_split_); + shared_.assign(n_tip, 0); + 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) { + const uint64_t *a = &bits1_[size_t(i) * n_word_]; + for (int j = 0; j != n_split_; ++j) { + const uint64_t *b = &bits2_[size_t(j) * n_word_]; + int a_and_b = 0; + for (int w = 0; w != n_word_; ++w) { + a_and_b += TreeDist::count_bits(a[w] & b[w]); + } + overlap_[size_t(i) * n_split_ + j] = int16_t(a_and_b); + const int a_xor_b = in1_[i] + in2_[j] - a_and_b - a_and_b; + 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; + } + } + } + } + } + return true; + } + + // Tree structure, as in reduce_trees() + struct Ledger { + std::vector child, sibling, parent; + void Build(const IntegerMatrix &edge, const int n_vertex) { + child.assign(n_vertex, 0); + sibling.assign(n_vertex, 0); + parent.assign(n_vertex, 0); + for (int e = edge.nrow(); e--; ) { + const int p = edge(e, 0); + const int c = edge(e, 1); + parent[c] = p; + const int existing = child[p]; + if (existing) { + sibling[existing] = c; + sibling[c] = existing; + } else { + child[p] = c; + } + } + } + // As REMOVE_TIP() in reduce_trees(); `tip`'s parent is not the root. + void RemoveTip(const int tip) { + const int p = parent[tip]; + const int s = sibling[tip]; + parent[s] = parent[p]; + sibling[sibling[p]] = s; + sibling[s] = sibling[p]; + child[parent[s]] = sibling[s]; + sibling[tip] = 0; + } + }; + + int Count(const uint64_t *x) const { + int ret = 0; + for (int w = 0; w != n_word_; ++w) { + ret += TreeDist::count_bits(x[w]); + } + 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) { + const int parent = ledger.parent[tip]; + const int sibling = ledger.sibling[tip]; + if (row_of[parent] >= 0) { + if (merge) *merge = row_of[sibling]; + return row_of[parent]; + } + // `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]; + } + + // 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); + 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, + &xk = xc_.child, &ys = yc_.sibling, &yp = yc_.parent, &yk = yc_.child; + // Cherries + for (int i = 1; i <= n_tip; ++i) { + const int s = xs[i]; + if (s && IsTip(s) && s == ys[i]) return true; + } + // Root + int x_1 = xk[xs[1]]; + if (!IsTip(x_1)) x_1 = xs[x_1]; + if (IsTip(x_1)) { + const int y_1 = yk[ys[1]]; + if (x_1 == y_1 || x_1 == ys[y_1]) return true; + } + // Chains + for (int i = n_tip; i >= 2; --i) { + if (!xs[i]) continue; + const int aunt = xs[xp[i]]; + if (!aunt || !IsTip(aunt)) continue; + const int g_aunt = xs[xp[aunt]]; + if (!g_aunt || !IsTip(g_aunt)) continue; + const int gg_aunt = xs[xp[g_aunt]]; + if (aunt == ys[yp[i]] && g_aunt == ys[yp[aunt]]) { + if (gg_aunt && IsTip(gg_aunt) && gg_aunt == ys[yp[g_aunt]]) { + return true; + } + } else if (i == ys[yp[aunt]] && aunt == ys[yp[g_aunt]]) { + const int niece1 = yk[ys[g_aunt]]; + if (gg_aunt && IsTip(gg_aunt) && + (gg_aunt == niece1 || gg_aunt == ys[niece1] || + gg_aunt == ys[g_aunt])) { + return true; + } + } + } + return false; + } + + // Parent prices for the prepared candidate's splits of tree 2. A split + // that two parent splits restrict to takes the lesser price. + static const price_t *Prices(void *context) { + SingleDrop &sd = *static_cast(context); + const std::vector &parent = sd.pp_->Price(); + if (parent.empty()) return nullptr; // #nocov + const int n_child = sd.n_split_ - 1; + sd.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]; + } + return sd.price_.data(); + } + + const int n_tip_; + const int n_word_; + const int n_split_; + const std::vector &bits1_; + const std::vector &bits2_; + bool ready_ = false; + Ledger x_, y_, xc_, yc_; + std::vector row_of1_, row_of2_, in1_, in2_; + std::vector overlap_; + std::vector shared_; + // The prepared candidate + 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_; +}; + } // Scores for dropping each leaf from `tree1` and `tree2`, as .SPRRogue(); @@ -686,8 +987,11 @@ double ScoreWithout(const IntegerMatrix &pre1, const IntegerMatrix &pre2, // combn(), if `depth` > 1 and every single-leaf score is finite. // 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. // [[Rcpp::export]] -List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { +List spr_rogue_scores(const List &tree1, const List &tree2, const int depth, + const bool incremental = true) { const CharacterVector labels = tree1["tip.label"]; const int n_tip = labels.size(); const IntegerMatrix pre1 = reverse(tree1["edge"]); @@ -705,24 +1009,33 @@ List spr_rogue_scores(const List &tree1, const List &tree2, const int depth) { for (int t = 0; t != n_tip; ++t) { identity[t] = t; } + std::vector bits1; + std::vector node1, node2; { - const TreeSplits full1 = GetSplits(tree1["edge"], n_tip); + const TreeSplits full1 = GetSplits(tree1["edge"], n_tip, &identity, + pp.n_word, &bits1, &node1); const TreeSplits full2 = GetSplits(tree2["edge"], n_tip, &identity, - pp.n_word, &pp.bits); + pp.n_word, &pp.bits, &node2); pp.splits1 = full1.bits; pp.splits2 = full2.bits; 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); 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(); - keep[i] = false; - scores[i] = ScoreWithout(pre1, pre2, labels, keep, best, pp_ptr); - keep[i] = true; + if (single.Prepare(i)) { + scores[i] = single.Score(best, pp_ptr); + } else { + keep[i] = false; + scores[i] = ScoreWithout(pre1, pre2, labels, keep, best, pp_ptr); + keep[i] = true; + } if (!R_FINITE(scores[i])) { all_finite = false; break; diff --git a/src/tree_distances.cpp b/src/tree_distances.cpp index 7189bf9d..e937e5e4 100644 --- a/src/tree_distances.cpp +++ b/src/tree_distances.cpp @@ -7,6 +7,7 @@ #include #include "tree_distances.h" #include "col_prices.h" +#include "split_counts.h" using Rcpp::_; using Rcpp::IntegerVector; @@ -463,6 +464,19 @@ static cost_matrix &BoundedPool(const int i, const std::size_t dim) { return pool[i]; } +// k log2 k for 0 <= k <= n_tips, for the bounded mutual clustering fill. +static const double *KLg2K(const int32 n_tips) { + static thread_local std::vector table; + if (table.size() <= static_cast(n_tips)) { + const int32 old_size = static_cast(table.size()); + table.resize(n_tips + 1); + for (int32 k = old_size; k <= n_tips; ++k) { + table[k] = k * lg2_lookup(k); + } + } + return table.data(); +} + template List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, const int32 n_tips, const double min_score, @@ -519,6 +533,22 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, std::vector a_match(a.n_splits); std::unique_ptr b_match = std::make_unique(b.n_splits); + // Bounded only: write each cell's information as + // sum_q f(n_q) - f(na) - f(nA) - f(nb) - f(nB) + f(n), + // where f(k) = k log2 k and q ranges over the four quadrants, hoisting the + // row and column terms. This takes four table lookups and additions per + // cell, against four lookups and multiplications; the result agrees with + // the unbounded form to within a unit of the integer cost. + [[maybe_unused]] static thread_local std::vector col_term; + [[maybe_unused]] const double *f = nullptr; + if constexpr (Bounded) { + f = KLg2K(n_tips); + col_term.resize(b.n_splits); + for (split_int bi = 0; bi < b.n_splits; ++bi) { + col_term[bi] = -f[b.in_split[bi]] - f[n_tips - b.in_split[bi]]; + } + } + for (split_int ai = 0; ai < a.n_splits; ++ai) { if (!Bounded && (ai & 1023) == 0) Rcpp::checkUserInterrupt(); if (a_match[ai]) continue; @@ -526,9 +556,13 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, const split_int na = a.in_split[ai]; const split_int nA = n_tips - na; const auto *a_row = a.state[ai]; - const double offset_a = lg2_n - lg2_lookup(na); - const double offset_A = lg2_n - lg2_lookup(nA); + [[maybe_unused]] const double offset_a = lg2_n - lg2_lookup(na); + [[maybe_unused]] const double offset_A = lg2_n - lg2_lookup(nA); [[maybe_unused]] cost row_min = max_score; + [[maybe_unused]] double row_term = 0; + if constexpr (Bounded) { + row_term = f[n_tips] - f[na] - f[nA]; + } for (split_int bi = 0; bi < b.n_splits; ++bi) { @@ -556,6 +590,10 @@ List mutual_clustering_impl(const RawMatrix &x, const RawMatrix &y, && a_and_b == a_and_B && a_and_b == A_and_B) { score(ai, bi) = max_score; // Avoid rounding errors + } else if constexpr (Bounded) { + const double ic_sum = f[a_and_b] + f[a_and_B] + f[A_and_b] + + f[A_and_B] + row_term + col_term[bi]; + score(ai, bi) = max_score - static_cast(ic_sum * max_over_tips); } else { const double lg2_nb = lg2_lookup(nb); const double lg2_nB = lg2_lookup(nB); @@ -902,6 +940,103 @@ double mutual_clustering_bounded_score(const RawMatrix &x, const RawMatrix &y, mutual_clustering(x, y, n_tip, min_score, prices)["score"]); } +// 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) { + const split_int n_tips = sc.n_tip; + const split_int dim = sc.n_split; + constexpr cost max_score = BIG; + constexpr double over_max_score = 1.0 / static_cast(max_score); + const double n_tips_reciprocal = 1.0 / n_tips; + const double max_over_tips = static_cast(max_score) * + n_tips_reciprocal; + + // Row bound, as mutual_clustering(). The entropy bound that precedes + // it there is omitted: it rejects a candidate only if the trees' entropies + // differ by more than the best score so far, which no candidate met in + // 600,000 tested. + double upper_bound = 0; + for (split_int i = 0; i < dim; ++i) { + upper_bound += TreeDist::ic_matching(sc.row_in[i], n_tips - sc.row_in[i], + n_tips); + } + upper_bound *= n_tips_reciprocal; + + const double *f = KLg2K(n_tips); + static thread_local std::vector col_term; + col_term.resize(dim); + for (split_int j = 0; j < dim; ++j) { + col_term[j] = -f[sc.col_in[j]] - f[n_tips - sc.col_in[j]]; + } + cost_matrix &score = BoundedPool(0, dim); + for (split_int i = 0; i < dim; ++i) { + const split_int na = sc.row_in[i]; + const split_int nA = n_tips - na; + const int16_t *overlap = sc.overlap + std::size_t(sc.row[i]) * sc.stride; + const int16_t row_has = sc.row_has[i]; + const double row_term = f[n_tips] - f[na] - f[nA]; + 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 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; + const split_int A_and_B = nA - A_and_b; + ASSERT(a_and_B || A_and_b); + ASSERT(a_and_b || A_and_B); + const cost c = (a_and_b == A_and_b && a_and_b == a_and_B && + a_and_b == A_and_B) ? max_score : + max_score - static_cast( + (f[a_and_b] + f[a_and_B] + f[A_and_b] + f[A_and_B] + row_term + + col_term[j]) * max_over_tips); + out[j] = c; + if (c < row_min) row_min = c; + } + upper_bound += static_cast(max_score - row_min) * over_max_score - + TreeDist::ic_matching(na, nA, n_tips) * n_tips_reciprocal; + if (upper_bound < min_score) { + return upper_bound; + } + } + + const cost *col_price = (prices && dim <= SL_MAX_SPLITS / 2) ? + prices->get(prices->context) : nullptr; + if (col_price) { + cost lb = 0; + for (split_int j = 0; j < dim; ++j) { + lb += col_price[j]; + } + for (split_int i = 0; i < dim; ++i) { + const cost *row = score.row(i); + cost row_min = row[0] - col_price[0]; + for (split_int j = 1; j < dim; ++j) { + const cost reduced = row[j] - col_price[j]; + if (reduced < row_min) row_min = reduced; + } + lb += row_min; + } + const double bound = static_cast(max_score * dim - lb) * + over_max_score; + if (bound < min_score) { + return bound; + } + } + + static thread_local std::vector rowsol; + static thread_local std::vector colsol; + static thread_local LapScratch warm_scratch; + resize_uninitialized(rowsol, dim); + resize_uninitialized(colsol, dim); + const cost lap_cost = col_price ? + lap_warm(dim, score, rowsol, colsol, warm_scratch, col_price) : + lap(dim, score, rowsol, colsol, false); // #nocov: > 16k leaves + return static_cast(max_score * dim - lap_cost) * over_max_score; +} + // 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 cd9953ab..f39db1c7 100644 --- a/tests/testthat/test-tree_distance_spr.R +++ b/tests/testthat/test-tree_distance_spr.R @@ -259,10 +259,15 @@ test_that(".BoundedCID() matches ClusteringInfoDist()", { tree1 <- RandomTree(nTip, root = TRUE) tree2 <- RandomTree(nTip, root = TRUE) exact <- ClusteringInfoDist(tree1, tree2) - expect_equal(.BoundedCID(tree1, tree2), exact, tolerance = 0) - expect_equal(.BoundedCID(tree1, tree2, exact), exact, tolerance = 0) + # The bounded fill sums its terms in a different order, so agrees with + # ClusteringInfoDist() only to within rounding error + expect_equal(.BoundedCID(tree1, tree2), exact, tolerance = 1e-12) + expect_equal(.BoundedCID(tree1, tree2, exact), exact, tolerance = 1e-12) # Within numerical noise of the cutoff: must be exact, to preserve ties - expect_equal(.BoundedCID(tree1, tree2, exact - 1e-12), exact, tolerance = 0) + expect_equal(.BoundedCID(tree1, tree2, exact - 1e-12), exact, + tolerance = 1e-12) + expect_identical(.BoundedCID(tree1, tree2, exact - 1e-12), + .BoundedCID(tree1, tree2)) if (exact > 0) { bounded <- .BoundedCID(tree1, tree2, exact / 2) expect_gt(bounded, exact / 2) @@ -321,6 +326,68 @@ test_that("C++ Rogue scorer matches R implementation", { } }) +test_that("Single-leaf 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 + # win may receive different lower bounds + Agree <- function(x, y) { + expect_equal(is.finite(x), is.finite(y)) + fin <- is.finite(x) + if (!any(fin)) return() + best <- min(x[fin]) + expect_equal(min(y[fin]), best, tolerance = 1e-12) + tol <- 1e-9 * max(1, best) + expect_equal(which(fin & x <= best + tol), which(fin & y <= best + tol)) + exact <- fin & x <= best + 1e-6 & y <= best + 1e-6 + expect_equal(x[exact], y[exact], tolerance = 1e-12) + } + Compare <- function(tree1, tree2) { + reduced <- ReduceTrees(tree1, tree2) + if (is.null(reduced) || NTip(reduced[[1]]) < 8) return(FALSE) + for (depth in 1:2) { + fast <- spr_rogue_scores(reduced[[1]], reduced[[2]], depth) + slow <- spr_rogue_scores(reduced[[1]], reduced[[2]], depth, FALSE) + Agree(fast[["scores"]], slow[["scores"]]) + Agree(fast[["pair_scores"]], slow[["pair_scores"]]) + } + TRUE + } + + # Removing a leaf leaves a chain of leaves common to both trees, which + # reduce_trees() would shorten + chain <- ape::read.tree(text = c( + "(t1,((((((t2,t11),t10),t15),t13),((((((t3,t8),t7),t14),t4),t6),t5)),(t9,t12)));", + "(t1,(((t8,((t15,((t2,t11),t10)),t13)),(((((t3,(t7,t9)),t14),t4),t6),t5)),t12));" + )) + expect_true(Compare(chain[[1]], chain[[2]])) + + Moved <- function(tree, nMoves) { + for (i in seq_len(nMoves)) { + leaf <- sample(tree[["tip.label"]], 1) + pruned <- DropTip(tree, leaf) + tree <- AddTip(pruned, sample(pruned[["edge"]][, 2], 1), leaf) + } + tree + } + set.seed(4174) + nScored <- 0 + for (i in 1:60) { + nTip <- sample(8:40, 1) + tree1 <- RandomTree(nTip, root = TRUE) + tree2 <- if (i %% 2) RandomTree(nTip, root = TRUE) else + Moved(tree1, sample(2:5, 1)) + nScored <- nScored + Compare(tree1, tree2) + } + expect_gt(nScored, 40) + + # Trees that share a split use the full rebuild for every candidate + tree1 <- Postorder(RootTree(BalancedTree(12), 1)) + tree2 <- Postorder(RootTree(RenumberTips(PectinateTree(12), tree1), 1)) + expect_equal(spr_rogue_scores(tree1, tree2, 1L), + spr_rogue_scores(tree1, tree2, 1L, FALSE)) +}) + test_that("cpp_mutual_clustering_bounded() exits early", { library("TreeTools", quietly = TRUE) MCI <- function(x, y, minScore) { From 7339e468f3d46c361e463866c23ee5a2a122c047 Mon Sep 17 00:00:00 2001 From: ms609-agent <313734811+ms609-agent@users.noreply.github.com> Date: Mon, 19 Oct 2026 16:29:10 +0100 Subject: [PATCH 2/2] Use parent overlaps only for trees rooted on leaf 1 The reduction predicate replays reduce_trees(), which expects trees rooted on leaf 1; other rootings fall back to rebuilding each candidate. Co-Authored-By: Claude Opus 5.5 --- src/spr.cpp | 11 +++++++++++ tests/testthat/test-tree_distance_spr.R | 6 ++++++ 2 files changed, 17 insertions(+) diff --git a/src/spr.cpp b/src/spr.cpp index 3f932823..544d6b8e 100644 --- a/src/spr.cpp +++ b/src/spr.cpp @@ -794,6 +794,9 @@ class SingleDrop { const int n_vertex = 2 * n_tip; x_.Build(edge1, n_vertex); y_.Build(edge2, n_vertex); + if (!x_.RootedOnLeaf1(n_tip) || !y_.RootedOnLeaf1(n_tip)) { + return false; + } row_of1_.assign(n_vertex, -1); row_of2_.assign(n_vertex, -1); in1_.resize(n_split_); @@ -860,6 +863,14 @@ class SingleDrop { } } } + // Whether the root's children are leaf 1 and an internal node, as + // reduce_trees() expects. + bool RootedOnLeaf1(const int n_tip) const { + const int root = parent[1]; + if (parent[root]) return false; + const int other = child[root] == 1 ? sibling[1] : child[root]; + return sibling[1] == other && other > n_tip; + } // As REMOVE_TIP() in reduce_trees(); `tip`'s parent is not the root. void RemoveTip(const int tip) { const int p = parent[tip]; diff --git a/tests/testthat/test-tree_distance_spr.R b/tests/testthat/test-tree_distance_spr.R index f39db1c7..ab4c083f 100644 --- a/tests/testthat/test-tree_distance_spr.R +++ b/tests/testthat/test-tree_distance_spr.R @@ -362,6 +362,12 @@ test_that("Single-leaf candidates scored from parent splits", { )) expect_true(Compare(chain[[1]], chain[[2]])) + # Trees not rooted on leaf 1 use the full rebuild for every candidate + tree1 <- Postorder(RootTree(chain[[1]], "t5")) + tree2 <- Postorder(RenumberTips(chain[[2]], tree1)) + expect_equal(spr_rogue_scores(tree1, tree2, 1L), + spr_rogue_scores(tree1, tree2, 1L, FALSE)) + Moved <- function(tree, nMoves) { for (i in seq_len(nMoves)) { leaf <- sample(tree[["tip.label"]], 1)