diff --git a/dev/profiling/drivers/tbr-accept-ab-cell.R b/dev/profiling/drivers/tbr-accept-ab-cell.R new file mode 100644 index 000000000..668b19580 --- /dev/null +++ b/dev/profiling/drivers/tbr-accept-ab-cell.R @@ -0,0 +1,41 @@ +# One A/B cell for the TBR-rerooting accept path — issue #38 (T-300) +# +# Prints ONE line so a shell loop can interleave the arms (this machine is +# shared, so alternating arms is what stops a drifting background load being +# read as an effect). The reported numbers come from the DLL's own +# steady_clock counters around the accept branch, not from R-level timing. +# +# Usage: TS_NA_TIMING=1 Rscript dev/profiling/drivers/tbr-accept-ab-cell.R \ +# [nCycles] + +args <- commandArgs(trailingOnly = TRUE) +libDir <- args[[1]] +arm <- args[[2]] +dsName <- args[[3]] +mode <- args[[4]] +seed <- as.integer(args[[5]]) +nCycles <- if (length(args) >= 6L) as.integer(args[[6]]) else 8L + +suppressMessages(library(TreeSearch, lib.loc = libDir)) +stopifnot(nzchar(Sys.getenv("TS_NA_TIMING"))) + +dataset <- TreeSearch::inapplicable.phyData[[dsName]] +at <- attributes(dataset) +tipData <- matrix(unlist(dataset, use.names = FALSE), + nrow = length(dataset), byrow = TRUE) +weight <- TreeSearch:::.ScaleWeight(at$weight) +concavity <- if (mode == "EW") -1 else 10 + +set.seed(seed) +tr <- ape::rtree(length(dataset), tip.label = names(dataset), rooted = FALSE) +startEdge <- ape::root(tr, 1L, resolve.root = TRUE)$edge + +set.seed(seed) +res <- TreeSearch:::ts_ratchet_search( + edge = startEdge, contrast = at$contrast, tip_data = tipData, + weight = weight, levels = at$levels, + nCycles = nCycles, perturbProb = 0.04, maxHits = 1L, concavity = concavity) + +cat(sprintf("%s,%s,%s,%d,%.4f,%.2f,%.2f,%.0f\n", + arm, dsName, mode, seed, res$score, + res$na_t_total_ms, res$na_t_accept_ms, res$na_n_accept)) diff --git a/dev/profiling/drivers/tbr-accept-reroot-audit.R b/dev/profiling/drivers/tbr-accept-reroot-audit.R new file mode 100644 index 000000000..0d46349b5 --- /dev/null +++ b/dev/profiling/drivers/tbr-accept-reroot-audit.R @@ -0,0 +1,96 @@ +# Score-identity battery for the TBR-rerooting accept path — issue #38 (T-300) +# +# Extending the dirty-set incremental rescore from SPR accepts to TBR-rerooting +# accepts claims to be SCORE-IDENTICAL: the incremental score must equal what +# full_rescore would have returned, so every accept/reject decision — and hence +# the whole search trajectory — must be unchanged. Run this against the +# pre-patch and post-patch libraries and diff the CSVs. +# +# Pair it with TS_TBR_ACCEPTCHK=1, which makes the DLL cross-check every +# incremental accept against full_rescore in-flight and abort on drift. That is +# the guard the reverted first attempt (b7303ee5, systematic delta = -3) lacked. +# +# Usage: +# TS_AUDIT_OUT=base.csv Rscript dev/profiling/drivers/tbr-accept-reroot-audit.R .agent-i38 +# TS_TBR_ACCEPTCHK=1 TS_AUDIT_OUT=patched.csv \ +# Rscript dev/profiling/drivers/tbr-accept-reroot-audit.R .agent-i38b + +args <- commandArgs(trailingOnly = TRUE) +libDir <- if (length(args) >= 1L) args[[1]] else ".agent-i38b" +library(TreeSearch, lib.loc = libDir) +library(TreeTools, quietly = TRUE) + +MakeData <- function(dataset) { + at <- attributes(dataset) + list( + contrast = at$contrast, + tipData = matrix(unlist(dataset, use.names = FALSE), + nrow = length(dataset), byrow = TRUE), + weight = at$weight, + levels = at$levels, + nTip = length(dataset) + ) +} + +# Weak-signal random matrices accept long chains of moves, which is what drives +# reroot accepts; the real matrices add NA blocks and realistic state counts. +cases <- list() +for (nTip in c(12L, 18L, 25L)) { + set.seed(1000 + nTip) + mat <- matrix(sample(0:3, nTip * 8L, replace = TRUE), nrow = nTip, + dimnames = list(paste0("t", seq_len(nTip)), NULL)) + cases[[paste0("rand", nTip)]] <- MatrixToPhyDat(mat) +} +data("inapplicable.phyData", package = "TreeSearch") +for (nm in c("Longrich2010", "Vinther2008", "Sansom2010", "DeAssis2011")) { + cases[[nm]] <- inapplicable.phyData[[nm]] +} + +rows <- list() +for (nm in names(cases)) { + dataset <- cases[[nm]] + d <- MakeData(dataset) + minSteps <- as.integer(MinimumLength(dataset, compress = TRUE)) + for (mode in c("EW", "IW")) { + searchConcavity <- if (mode == "EW") -1 else 10 + scoreConcavity <- if (mode == "EW") Inf else 10 + ms <- if (mode == "EW") integer(0) else minSteps + for (start in c(1, 17, 88, 256, 777)) { + tree <- as.phylo(start, d$nTip) + set.seed(start) + res <- TreeSearch:::ts_tbr_search( + tree$edge, d$contrast, d$tipData, d$weight, d$levels, + maxHits = 50L, min_steps = ms, concavity = searchConcavity) + independent <- TreeSearch:::ts_fitch_score( + res$edge, d$contrast, d$tipData, d$weight, d$levels, + min_steps = ms, concavity = scoreConcavity) + rows[[length(rows) + 1L]] <- data.frame( + case = nm, mode = mode, start = start, + score = res$score, independent = independent, + nAccepted = res$n_accepted, nEvaluated = res$n_evaluated, + stringsAsFactors = FALSE) + } + } +} + +tab <- do.call(rbind, rows) +tab$drift <- tab$score - tab$independent +bad <- tab[abs(tab$drift) > 1e-9, , drop = FALSE] + +cat(sprintf("cells: %d | total accepts: %d | in-flight audit: %s\n", + nrow(tab), sum(tab$nAccepted), + if (nzchar(Sys.getenv("TS_TBR_ACCEPTCHK"))) "ON" else "off")) +if (nrow(bad)) { + cat("SCORE DRIFT vs independent recomputation:\n") + print(bad, row.names = FALSE) +} else { + cat("all reported scores match an independent full recomputation\n") +} + +outFile <- Sys.getenv("TS_AUDIT_OUT", unset = "") +if (nzchar(outFile)) { + write.csv(tab[, c("case", "mode", "start", "score", "nAccepted", "nEvaluated")], + outFile, row.names = FALSE) + cat("wrote", outFile, "\n") +} +if (nrow(bad)) quit(status = 1L) diff --git a/dev/profiling/findings.md b/dev/profiling/findings.md index 8da0ecc59..486dd7b44 100644 --- a/dev/profiling/findings.md +++ b/dev/profiling/findings.md @@ -13,6 +13,7 @@ Tags: | ID-suggest | P? | Status | Depends | Headline | Detail (% time, mechanism, verified Δ, micro-bench path) | |------------|----|--------|---------|----------|---------------------------------------------------------| | T-300 | P1 | DONE | — | [Optimise] `full_rescore` after accepted TBR move (ts_tbr.cpp:1138): replace with incremental rescore | LANDED (commits f531bbcd EW + 014ccdea NA dirty-set). 19.2 % of NA-path DLL CPU; 15.2 % wall speedup on Zhu2013 NA (3.88→3.29 s). | +| #38 | P3 | DONE / AT-LIMIT | T-300 | [AT-LIMIT] the residue T-300 left: `full_rescore` on TBR-**rerooting** accepts | The rerooting arm is now incremental too (third dirty seed at `clip_node`), but the lever is spent. Post-T-300 the WHOLE accept branch is **0.18–0.51 %** of `tbr_search` wall (Vinther2008 / Agnarsson2004 / Zanol2014 × EW,IW; in-DLL `na_t_accept_ms`, 15 paired cells, `dev/profiling/drivers/tbr-accept-ab-cell.R`). The patch cuts that slice to ×0.901 median (12/15 cells) ⇒ **≈0.03 % of wall e2e — undetectable**. Merge case is codepath unification, not speed. Rerooting accepts are ~18 % of accepts, so this is not a coverage artefact. | ## Round 3 (2026-06-16) — standard-Fitch TNT-parity path (Zhu2013 `-`→`?`, auto→thorough) diff --git a/dev/profiling/focus-areas.md b/dev/profiling/focus-areas.md index f6f9aaaac..71d8cd457 100644 --- a/dev/profiling/focus-areas.md +++ b/dev/profiling/focus-areas.md @@ -21,7 +21,7 @@ further wins — skip unless code changes), `SKIPPED` (out of rotation). | 1 | NNI-perturb in driven pipeline | `src/ts_nni_perturb.cpp`, `src/ts_driven.cpp` (perturb call sites) | Disabled in thorough preset via T-274 (`nniPerturbCycles=0L` in `R/MaximizeParsimony.R`) — code on path only when caller sets `nni_perturb_per > 0` | T-274 filed; disabled at R level; re-evaluate only if default changes | 2026-05-18 | AT-LIMIT | | 2 | Ratchet inner loop | `src/ts_ratchet.cpp`, `src/ts_tbr.cpp` (called from ratchet) | 62 % of inner-loop search time (verbosity=2, Zhu2013 thorough, 2026-05-18); TBR dominates (perturbation overhead < 2 %) | 2.80 s/rep median (Zhu2013 thorough ×1 rep, nThreads=1); T-300 (`full_rescore`) is pending fix | 2026-05-18 | PROFILED | | 3 | RSS / sector search | `src/ts_sector.cpp`, `src/ts_prune_reinsert.cpp` | THROUGHPUT at-limit by inheritance (R6 2026-06-20): ~96 % is inner+global tbr_search (at-limit kernel); sector scaffolding ≤2 %. Banked T-S6c byte-identical ~2.8 %; T-S6d per-clip getenv ~22 % (TBR-wide). Efficiency axis (work-to-target) untouched. | inner tbr_search-dominated; see findings R6 | 2026-06-20 | AT-LIMIT | -| 4 | TBR full-rescore at acceptance | `src/ts_tbr.cpp:1138` (`full_rescore` after every accepted move) | T-300 RESOLVED — dirty-set incremental rescore landed for SPR accept (EW path `fitch_dirty_*`, NA path `fitch_na_dirty_*`); GHA-green; 15.2 % wall-time speedup on Zhu2013 (3.88 s → 3.29 s) confirmed via dev/profiling/t300_na_bench.R 2026-05-19 | resolved | 2026-05-19 | DONE | +| 4 | TBR full-rescore at acceptance | `src/ts_tbr.cpp` accept branch (`full_rescore` after an accepted move) | T-300 RESOLVED — dirty-set incremental rescore landed for SPR accept (EW path `fitch_dirty_*`, NA path `fitch_na_dirty_*`); GHA-green; 15.2 % wall-time speedup on Zhu2013 (3.88 s → 3.29 s) confirmed via dev/profiling/t300_na_bench.R 2026-05-19. Issue #38 then closed the rerooting residue (third dirty seed) and **measured the branch at-limit: 0.18–0.51 % of `tbr_search` wall, patch ×0.901 on that slice ⇒ ≈0.03 % e2e**. Do not re-profile. | at-limit | 2026-08-04 | AT-LIMIT | | 5 | quartet_concordance.cpp allocation | `src/quartet_concordance.cpp` | T-298 active PR #242 — matrix allocation hoist already benchmarked; re-profile after merge | hoist-fix in flight | 2026-05-12 | PROFILED | | 6 | CSS / XSS sector pipeline | `src/ts_sector.cpp`, `src/ts_simplify.cpp` (`ts_simplify_diag` entry) | Same verdict as #3 (R6): XSS uses search_sector (=RSS scaffolding ≤2 % + inner tbr_search); CSS uses sector-masked tbr_search directly — both inner-tbr-dominated ⇒ THROUGHPUT at-limit by inheritance. T-S6c levers + T-S6d getenv apply to all three modes. | inner tbr_search-dominated | 2026-06-20 | AT-LIMIT | | 7 | Hierarchical resampling parallelism | `src/ts_resample.cpp`, `R/Resample.R` | HSJ/XFORM hierarchical resampling 2-thread speedup 1.1× (vs Brazeau 2.5×) — serial R loop | known limitation (2026-03-19 Agent A) | — | NEW | diff --git a/src/ts_data.h b/src/ts_data.h index 5c078fda2..89e0d76c7 100644 --- a/src/ts_data.h +++ b/src/ts_data.h @@ -268,6 +268,20 @@ struct DataSet { mutable long long na_t_vroot_ns = 0; // vroot_cache build / compute_from_above mutable long long na_t_accept_ns = 0; // accept-path NA dirty rescores mutable long long na_n_accept = 0; + // Moves that rerooted the clipped fragment AND were rescored by the dirty-set + // accept path (i.e. incremental_ok held), counted ALWAYS (not only under + // TS_NA_TIMING) for the same reason as na_n_evs above: it is one increment per + // rescore, and it is the only evidence that the reroot arm of that path was + // reached at all. Without it the regression test for that arm cannot tell a + // correct rescore from an unexercised one — so it must count the arm, NOT the + // move class: a rerooting move under HSJ/XFORM takes full_rescore instead, and + // counting it here would let the test report coverage of code that never ran. + // + // Counted at the rescore, i.e. BEFORE the accept/reject decision: a move + // counted here may still be rejected by the constraint check, the tabu test or + // the score comparison. It is therefore an upper bound on accepted rerootings + // — read it against the arm it witnesses, not against n_accepted. + mutable long long n_reroot_accepts = 0; // Per-pattern step scratch for the weighted (IW/profile) full-rescore path // (fitch_score_ew). Lives on DataSet for the SAME reason as evs_false_cache diff --git a/src/ts_fitch.cpp b/src/ts_fitch.cpp index 10e491497..7e8a19485 100644 --- a/src/ts_fitch.cpp +++ b/src/ts_fitch.cpp @@ -287,9 +287,15 @@ void fitch_incremental_uppass(TreeState& tree, const DataSet& ds, // each affected node exactly once in postorder, reading current children's // prelims — which are guaranteed correct because postorder processes // children before parents. +// +// A TBR rerooting additionally rewrites the children of every node on +// clip_node..reroot_parent; passing clip_node as start_c covers them (see +// ts_fitch.h). Off-path nodes inside the moved fragment keep both their +// children and their whole subtree, so their prelim and local_cost are +// untouched and the returned delta stays exact. int fitch_dirty_downpass(TreeState& tree, const DataSet& ds, - int start_a, int start_b) { + int start_a, int start_b, int start_c) { std::vector dirty(tree.n_node, 0); // Mark the rootward path from `node` up to (and including) the root. @@ -304,6 +310,7 @@ int fitch_dirty_downpass(TreeState& tree, const DataSet& ds, }; mark_path(start_a); mark_path(start_b); + if (start_c >= 0) mark_path(start_c); int length_delta = 0; @@ -353,7 +360,7 @@ int fitch_dirty_downpass(TreeState& tree, const DataSet& ds, } void fitch_dirty_uppass(TreeState& tree, const DataSet& ds, - int start_a, int start_b) { + int start_a, int start_b, int start_c) { // Step 1: root final_ = prelim (root prelim may have changed in downpass). int root = tree.n_tip; size_t root_base = static_cast(root) * tree.total_words; @@ -375,6 +382,7 @@ void fitch_dirty_uppass(TreeState& tree, const DataSet& ds, }; mark_path(start_a); mark_path(start_b); + if (start_c >= 0) mark_path(start_c); // Step 3: reverse postorder — visit any node whose parent is dirty_up. // If that node's final_ changes, propagate the flag to it. diff --git a/src/ts_fitch.h b/src/ts_fitch.h index 9fcb394e9..c569ae7e4 100644 --- a/src/ts_fitch.h +++ b/src/ts_fitch.h @@ -56,13 +56,20 @@ int fitch_incremental_downpass(TreeState& tree, const DataSet& ds, void fitch_incremental_uppass(TreeState& tree, const DataSet& ds, int start_node); -// Dirty-set rescore after an SPR move (T-300). +// Dirty-set rescore after an SPR or TBR move (T-300). // -// Recomputes prelim and local_cost for every node on the union of paths -// start_a -> root and start_b -> root, visiting each node exactly once in -// postorder. start_a and start_b are the two clip endpoints whose children -// changed after apply_tbr_move (typically nz = clip grandparent and -// nx = regraft point). +// Recomputes prelim and local_cost for every node on the union of the rootward +// paths from start_a, start_b and (optionally) start_c, visiting each node +// exactly once in postorder. start_a and start_b are the two clip endpoints +// whose children changed after apply_tbr_move (typically nz = clip grandparent +// and nx = regraft point). +// +// start_c: optional third dirty seed (-1 = unused), needed only when the move +// rerooted the clipped fragment. apply_tbr_move reverses the parent/child +// links along clip_node..reroot_parent, so every node on that path gains new +// children; after the reversal that path IS clip_node's rootward chain, so +// seeding at clip_node marks exactly those nodes (and nothing else new). +// Mirrors fitch_na_dirty_downpass's start_c. // // Caller must call tree.build_postorder_prealloc() first so that // tree.postorder reflects the post-move topology. @@ -71,21 +78,22 @@ void fitch_incremental_uppass(TreeState& tree, const DataSet& ds, // For IW/profile, ignore the return value and use extract_char_steps + // compute_weighted_score after this call (local_cost is correct). int fitch_dirty_downpass(TreeState& tree, const DataSet& ds, - int start_a, int start_b); + int start_a, int start_b, int start_c = -1); // Companion uppass for fitch_dirty_downpass. Recomputes final_ for nodes // whose ancestor's final_ may have changed, seeded from the same start // points. Propagates downward. void fitch_dirty_uppass(TreeState& tree, const DataSet& ds, - int start_a, int start_b); + int start_a, int start_b, int start_c = -1); // --- NA-aware dirty-set incremental rescore (T-300 NA variant) --- // // Same dirty-set approach as fitch_dirty_downpass / fitch_dirty_uppass but // handles inapplicable-bearing blocks via the NA-aware Pass 1 / Pass 2 -// logic. Used for the SPR accept path under has_inapplicable to avoid -// full_rescore. The return value is the EW length delta for standard -// blocks only — NA block step counts require Pass 3, so call +// logic. Used for the SPR and TBR-rerooting accept paths under +// has_inapplicable to avoid full_rescore, and by exact_verify_sweep's +// incremental candidate rescore. The return value is the EW length delta for +// standard blocks only — NA block step counts require Pass 3, so call // fitch_na_pass3_score(tree, ds) on the updated state to obtain the // authoritative score. // start_c: optional third dirty seed (-1 = unused). The TBR-reroot dirty region diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 8ff590920..1a035fff1 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -924,7 +924,8 @@ List ts_tbr_search( Named("na_t_vroot_ms") = ds.na_t_vroot_ns / 1e6, Named("na_t_accept_ms") = ds.na_t_accept_ns / 1e6, Named("na_n_accept") = static_cast(ds.na_n_accept), - Named("n_candidates") = static_cast(ds.n_candidates_evaluated) + Named("n_candidates") = static_cast(ds.n_candidates_evaluated), + Named("n_reroot_accepts") = static_cast(ds.n_reroot_accepts) ); } @@ -985,7 +986,8 @@ List ts_ratchet_search( Named("na_t_vroot_ms") = ds.na_t_vroot_ns / 1e6, Named("na_t_accept_ms") = ds.na_t_accept_ns / 1e6, Named("na_n_accept") = static_cast(ds.na_n_accept), - Named("n_candidates") = static_cast(ds.n_candidates_evaluated) + Named("n_candidates") = static_cast(ds.n_candidates_evaluated), + Named("n_reroot_accepts") = static_cast(ds.n_reroot_accepts) ); } diff --git a/src/ts_tbr.cpp b/src/ts_tbr.cpp index 9dded3351..a912e1778 100644 --- a/src/ts_tbr.cpp +++ b/src/ts_tbr.cpp @@ -1663,6 +1663,8 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, refresh_collapsed_all_zero(); const bool revert_check = std::getenv("TS_REVERT_CHECK") != nullptr; const bool iw_scanchk = std::getenv("TS_IW_SCANCHK") != nullptr; + // Oracle for the dirty-set accept path (see TS_TBR_ACCEPTCHK below). + const bool acceptchk = std::getenv("TS_TBR_ACCEPTCHK") != nullptr; // TS_PHYS_REROOT selects the legacy physical-reroot reference path; it is read // once per outer reroot-loop iteration below (>=1/call), so hoist it too. const bool phys_reroot = std::getenv("TS_PHYS_REROOT") != nullptr; @@ -2790,21 +2792,32 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, tree.build_postorder_prealloc(work_stack); - // T-300: dirty-set incremental rescore for SPR moves. The two - // affected nodes after apply_tbr_move are nz (clip grandparent, - // children changed: nx -> ns) and nx (regraft point, children - // changed to {clip_node, below}). fitch_dirty_downpass updates - // every node on the union of paths nz->root and nx->root exactly - // once in postorder; sums correctly with no shared-ancestor - // ambiguity. TBR moves with non-trivial rerooting and NA - // datasets fall back to full_rescore. + // T-300: dirty-set incremental rescore. The affected nodes after + // apply_tbr_move are nz (clip grandparent, children changed: + // nx -> ns), nx (regraft point, children changed to + // {new_subtree_root, below}) and — for a TBR rerooting only — every + // node on clip_node..reroot_parent, whose parent/child links the move + // reverses. After that reversal the path IS clip_node's rootward + // chain, so a third seed at clip_node marks exactly those nodes; + // `third` stays -1 for SPR, leaving the two-seed set untouched. + // fitch_dirty_downpass updates every node on the union of the seeds' + // rootward paths exactly once in postorder; sums correctly with no + // shared-ancestor ambiguity. Only scoring modes whose total equals + // the Fitch/IW result (see incremental_ok) take this path; HSJ and + // XFORM still fall back to full_rescore. bool is_spr = (best_reroot_parent < 0 || best_reroot_parent == clip_node); + const int third = (!is_spr && clip_node >= tree.n_tip) ? clip_node : -1; + // Count the ARM, not the move class: under HSJ/XFORM a rerooting move + // still falls through to full_rescore below, and counting it here would + // let the regression test in test-ts-tbr-dirty-rescore.R report coverage + // of a dirty-set rescore that never ran (src/ts_data.h). + if (!is_spr && incremental_ok) ++ds.n_reroot_accepts; double actual; const auto _t_acc = na_timing ? std::chrono::steady_clock::now() : std::chrono::steady_clock::time_point{}; - if (is_spr && !has_na && incremental_ok) { - int delta = fitch_dirty_downpass(tree, ds, nz, nx); - fitch_dirty_uppass(tree, ds, nz, nx); + if (!has_na && incremental_ok) { + int delta = fitch_dirty_downpass(tree, ds, nz, nx, third); + fitch_dirty_uppass(tree, ds, nz, nx, third); if (use_iw) { std::fill(divided_steps.begin(), divided_steps.end(), 0); extract_char_steps(tree, ds, divided_steps); @@ -2812,14 +2825,15 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, } else { actual = best_score + static_cast(delta); } - } else if (is_spr && has_na && incremental_ok) { + } else if (has_na && incremental_ok) { // T-300 NA variant: dirty-set Pass 1 + Pass 2 instead of full // rescore. Pass 3 still runs over the full tree because it // populates internal down2 (read by extract_char_steps) and // counts NA-block steps directly. Savings come from skipping - // Pass 1 + Pass 2 on off-dirty nodes. - fitch_na_dirty_downpass(tree, ds, nz, nx); - fitch_na_dirty_uppass(tree, ds, nz, nx); + // Pass 1 + Pass 2 on off-dirty nodes. The same three-seed dirty + // region already backs the exact_verify_sweep incremental path. + fitch_na_dirty_downpass(tree, ds, nz, nx, third); + fitch_na_dirty_uppass(tree, ds, nz, nx, third); int ew_total = fitch_na_pass3_score(tree, ds); if (use_iw) { std::fill(divided_steps.begin(), divided_steps.end(), 0); @@ -2833,19 +2847,38 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, actual = static_cast(ew_total) + ds.ew_offset; } } else { - // Non-trivial TBR rerooting, or a scoring mode whose incremental - // delta is not exact (HSJ/XFORM, see incremental_ok): recompute the - // authoritative score via score_tree(). + // A scoring mode whose incremental delta is not exact (HSJ/XFORM, + // see incremental_ok): recompute the authoritative score via + // score_tree(). actual = full_rescore(tree, ds); } // Accept-path rescore: the price of ACCEPTING a move, as distinct from // scanning candidates. On NA this is the dirty down/uppass plus a - // full-tree Pass 3, or an outright full_rescore for a TBR rerooting. + // full-tree Pass 3; HSJ/XFORM still pay an outright full_rescore. if (na_timing) { ds.na_t_accept_ns += ns_since(_t_acc); ++ds.na_n_accept; } + // AUDIT (env TS_TBR_ACCEPTCHK): cross-check the incremental accept + // score against full_rescore and abort on any drift. This is the + // oracle for the dirty-set accept path — an earlier incremental + // attempt shipped a systematic delta of -3 (b7303ee5) precisely + // because no such check existed. full_rescore leaves prelim/final_ + // coherent for the whole tree, so running it here is state-neutral. + // No-op unless the env var is set. + if (acceptchk && incremental_ok) { + const double audit = full_rescore(tree, ds); + if (std::fabs(actual - audit) > 1e-6) { + Rcpp::stop("TS_TBR_ACCEPTCHK mismatch mode=%s reroot=%d clip=%d " + "incr=%.6f full=%.6f diff=%.6f", + has_na ? (use_iw ? "NA+IW" : "NA+EW") + : (use_iw ? "IW" : "EW"), + is_spr ? 0 : 1, clip_node, actual, audit, + actual - audit); + } + } + // DIAGNOSTIC (env TS_IW_SCANCHK): compare the scan's predicted // best_candidate against the authoritative post-apply score for EVERY // scorer. Pure EW should be 0 (its indirect length is exact); a diff --git a/tests/testthat/test-ts-tbr-dirty-rescore.R b/tests/testthat/test-ts-tbr-dirty-rescore.R index ccd0329f8..421b396a1 100644 --- a/tests/testthat/test-ts-tbr-dirty-rescore.R +++ b/tests/testthat/test-ts-tbr-dirty-rescore.R @@ -120,6 +120,61 @@ test_that("TBR dirty-set rescore matches full rescore (NA-IW dataset, many accep } }) +test_that("dirty-set rescore matches full rescore on TBR-REROOTING accepts", { + # Issue #38: the dirty-set accept path originally covered SPR-classified + # accepts only; accepts that rerooted the clipped fragment fell back to + # full_rescore. Extending it adds a third dirty seed at clip_node, because + # apply_tbr_move reverses the parent/child links along + # clip_node..reroot_parent and so gives every node on that path new children. + # + # The four tests above cannot guard this arm: they assert score identity but + # have no way to tell whether a rerooting accept ever occurred, so they would + # pass just as happily if the arm were never entered. `n_reroot_accepts` + # (src/ts_data.h) is what makes this one non-vacuous -- it is asserted + # positive, so losing coverage fails the test rather than silently voiding it. + data("inapplicable.phyData", package = "TreeSearch") + minSteps <- function(dataset) { + as.integer(MinimumLength(dataset, compress = TRUE)) + } + + set.seed(6273) + mat <- matrix(sample(0:3, 20 * 8, replace = TRUE), + nrow = 20, dimnames = list(paste0("t", 1:20), NULL)) + random20 <- MatrixToPhyDat(mat) + vinther <- inapplicable.phyData[["Vinther2008"]] + + cases <- list( + list(label = "EW", dataset = random20, concavity = -1, score_conc = Inf), + list(label = "IW", dataset = random20, concavity = 10, score_conc = 10), + list(label = "NA", dataset = vinther, concavity = -1, score_conc = Inf), + list(label = "NA-IW", dataset = vinther, concavity = 10, score_conc = 10) + ) + + for (case in cases) { + ds <- make_ts_data(case$dataset) + n_tip <- length(case$dataset) + ms <- if (is.finite(case$score_conc)) minSteps(case$dataset) else integer(0) + n_reroot <- 0 + + for (start in c(3, 29, 131, 512, 900)) { + tree <- as.phylo(start, n_tip) + set.seed(7000 + start) + result <- ts_tbr(tree, ds, maxHits = 50L, concavity = case$concavity, + min_steps = ms) + n_reroot <- n_reroot + result$n_reroot_accepts + + rt <- result_tree(result, tree) + independent <- ts_score(rt, ds, concavity = case$score_conc, + min_steps = ms) + expect_equal(result$score, independent, tolerance = 1e-10, + info = paste(case$label, "start =", start)) + validate_result(result, n_tip) + } + + expect_gt(n_reroot, 0) # coverage: the rerooting arm was actually entered + } +}) + test_that("XPIWE x4 + dirty-region opts are byte-identical to opts-off (port guard)", { # Regression guard for the IW->XPIWE opt port (src/ts_tbr.cpp `iw_family` # gate): the x4 reroot batch + extract_char_steps dirty-region must produce