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/.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/NEWS.md b/NEWS.md index df06e482..cdac18dd 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). @@ -51,6 +55,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()`. 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/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/_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/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) }) 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 diff --git a/tests/testthat/testthat-problems.rds b/tests/testthat/testthat-problems.rds new file mode 100644 index 00000000..47bbda63 Binary files /dev/null and b/tests/testthat/testthat-problems.rds differ