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..544d6b8e 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,306 @@ 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); + 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_); + 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; + } + } + } + // 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]; + 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 +998,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 +1020,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..ab4c083f 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,74 @@ 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]])) + + # 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) + 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) {