diff --git a/NEWS.md b/NEWS.md index a96ed54a9..d26ba3532 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,33 @@ # To integrate into 2.0.0 notes +- `inapplicable = "xform"` scores are now reported at a canonical rooting, so a + reported score is reproducible. The x-transformation's step matrix is + asymmetric -- a gain costs one more than the number of secondary characters it + brings into existence, against 1 for a loss -- which makes a tree's length + depend on where it is rooted, unlike parsimony under the symmetric criteria. + `MaximizeParsimony()` recorded its best score mid-search at whatever rooting + the replicate held, while returning trees re-rooted on the first taxon, so + `attr(result, "score")` did not match `TreeLength()` of the very tree returned + (measured: 178 reported against 183 returned on a 36-taxon matrix), and + re-rooting a returned tree changed its length again. Both boundaries now + canonicalise on the first taxon of `dataset`, so one topology has one length + and the two agree by construction. + + **X-transformation scores may therefore differ slightly from previous + versions**, and will not decrease: the reported value is the length of the tree + you are handed rather than of a rooting discarded during search. It is an + upper bound on the rooting-free minimum, exceeding it by at most the total + number of secondary characters across hierarchy blocks (attained exactly by + 87--98% of rootings in simulation). This changes reporting only -- what the + search optimises is untouched. + + `MaximizeParsimony()` now also warns when the trees it returns do not share a + length at that common rooting, which can happen because pool membership is + still decided on scores taken at differing rootings. Only the x-transformation + is affected; HSJ reporting is deliberately unchanged, since there + rooting-invariance is a property the method requires rather than a convention + to pick. + - `MaximizeParsimony(effort = )` replaces `strategy = `, which is removed (it was never released). `effort` is a **relative** offset, not an absolute level: `0` (the default) accepts the amount of search the dataset's size and diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 1ebc660b5..4116d84ef 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -673,7 +673,18 @@ #' \item{`"xform"`}{Step-matrix recoding approximating maximum homology #' via x-transformations #' \insertCite{Goloboff2021;textual}{TreeSearch}. Requires a -#' `hierarchy`.} +#' `hierarchy`. **Scores are rooting-sensitive**: the step matrix of +#' this recoding is asymmetric -- gaining the controlling character costs one more +#' than the number of secondaries it brings into existence, against 1 to +#' lose it -- so a tree's length depends on where its root sits, whereas +#' parsimony under the other methods does not. Lengths are therefore +#' reported at a canonical rooting, on the first taxon of `dataset`, which +#' is the rooting the returned trees carry; `TreeLength()` canonicalises +#' identically, so it reproduces the reported score and one topology has +#' one length. That value is an upper bound on the rooting-free minimum, +#' exceeding it by at most the total number of secondary characters across +#' hierarchy blocks, and attaining it for 87-98% of rootings in +#' simulation.} #' } #' @param hsj_alpha Numeric in \[0, 1\]: scaling parameter for secondary- #' character contributions under the HSJ method. 0 = secondaries ignored; @@ -867,6 +878,13 @@ #' the true optimum, so reporting them would look like erratic progress. Any #' score the heartbeat prints is therefore directly comparable with the final #' tree score. +#' @param .rung Internal. Pins a named entry of the effort ladder +#' (`"sprint"`, `"default"`, `"thorough"`, `"large"`), or `"none"` to apply no +#' preset at all; `NULL` (default) selects the rung from `effort` and the +#' dataset's size, which is what every ordinary call should do. Exists for +#' controlled experiments and the preset smoke tests, which need to name a rung +#' absolutely rather than relative to the automatic choice. Not part of the +#' stable interface: prefer `effort`. #' @param control A [`SearchControl`] object (or a named list) of low-level #' search parameters. Most users can rely on `effort` and #' ignore this argument; see [`SearchControl()`] for full documentation @@ -1601,7 +1619,9 @@ MaximizeParsimony <- function( # Contract zero-length (unsupported) branches into polytomies, à la TNT's # "collapse zero-length branches" -- done entirely in C++ (ts_collapse_pool) # to avoid a per-tree R surgery quagmire. The kernel re-roots each tree on - # tip 0 (so root-adjacent edges are trivial -> rooting-invariant collapse), + # tip 0 (so root-adjacent edges are trivial -> rooting-invariant *contraction*; + # note the LENGTH is not rooting-invariant under HSJ/XFORM, T-374, which is + # why the XFORM pool is rescored at this rooting below), # flags aggressive (min-length-0) internal edges in the *search's* scoring # mode, contracts them, and deduplicates on the collapsed topology. # @@ -1650,6 +1670,50 @@ MaximizeParsimony <- function( outTrees <- list(treeTpl) } + # --- XFORM: report the score of the tree we are actually returning --- + # `result$best_score` is recorded mid-search at whatever rooting the replicate + # held. XFORM's step matrix is asymmetric, so the score is rooting-dependent, + # and `ts_collapse_pool()` above hands back every tree re-rooted on tip 0. + # Reporting `best_score` therefore gives the user a number that `TreeLength()` + # of the returned tree does not reproduce -- measured at 178 reported against + # 183 returned (T-385; repro in dev/red-team/heavy-tests/). Rescore the + # returned pool at the canonical rooting instead: |pool| evaluations, negligible + # against a search, and `TreeLength()` canonicalises identically, so the two + # agree by construction. + # + # This deliberately does NOT change what the search optimises. The reported + # value stays a rooting-dependent upper bound on the min-over-rootings + # objective, exceeding it by at most the sum of `nSec` over hierarchy blocks + # (measured: attained by 87-98% of rootings, mean overstatement 0.02-0.17 + # steps). + # + # Min-over-rootings reporting -- the variant the plan calls better -- is NOT used. + # It reports a quantity the search never compared, but the deciding objection is + # cost: evaluated naively it is (2 * nTip - 3) x on the Sankoff term, so a + # 4000-tip pool of 100 trees would need ~800k evaluations at the boundary. + # Doing it affordably needs an all-rootings up-down DP, which is its own piece + # of work (Option 4), not a line in a reporting fix. + # See dev/plans/2026-07-29-t374b-xform-rooting-policy.md (Option 3). + bestScore <- result$best_score + if (useXform && length(outTrees) > 0L) { + canonicalScores <- TreeLength( + structure(outTrees, class = "multiPhylo"), + dataset, inapplicable = "xform", hierarchy = hierarchy + ) + bestScore <- min(canonicalScores) + if (diff(range(canonicalScores)) > sqrt(.Machine$double.eps)) { + # Pool membership is chosen on search-time scores taken at differing + # rootings (`result$scores` above), so trees held to be equally + # parsimonious can differ once scored at one rooting. Not silently + # averaged away: this is the open residue of T-374, and staying quiet about + # it is what let the reporting gap survive this long. + warning("Returned trees do not share a length at a common rooting (", + paste(signif(range(canonicalScores), 8), collapse = " to "), + "); reporting the smallest. The x-transformation's score is ", + "rooting-dependent -- see ?MaximizeParsimony.") + } + } + # --- Output --- if (verbosity > 0L) { total_s <- round(sum(unlist(result$timings), na.rm = TRUE) / 1000, 1) @@ -1658,7 +1722,7 @@ MaximizeParsimony <- function( else if (isTRUE(result$perturb_stop)) "perturbation limit" else "replicate limit" cli_alert_success(paste0( - "Search complete: score {.strong {signif(result$best_score, 7)}}, ", + "Search complete: score {.strong {signif(bestScore, 7)}}, ", "{result$replicates} replicate{?s} ", "(last improved: #{result$last_improved_rep}), ", "{result$hits_to_best} hit{?s} to best, ", @@ -1669,7 +1733,7 @@ MaximizeParsimony <- function( structure( outTrees, - score = result$best_score, + score = bestScore, replicates = result$replicates, hits_to_best = result$hits_to_best, n_topologies = nTopologies, diff --git a/R/RcppExports.R b/R/RcppExports.R index afce48989..dc9d764e1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -40,7 +40,6 @@ MaddisonSlatkin_clear_cache <- function() { #' @param nj Integer vector giving the block sizes of the second partition #' (also summing to `N`). #' @return The expected mutual information, in bits. -#' @references \insertAllCited{} #' @seealso [`SiteConcordance`] #' @examples #' # Expected MI between a 3|4 split and a 2|5 split of 7 items: diff --git a/R/recode_hierarchy.R b/R/recode_hierarchy.R index 34c1d3f92..2cc435a50 100644 --- a/R/recode_hierarchy.R +++ b/R/recode_hierarchy.R @@ -22,6 +22,24 @@ #' - **Present → present:** Hamming distance (number of secondaries with #' different states). #' +#' ## Rooting +#' +#' Because gain and loss cost differently, this cost matrix is **asymmetric** +#' whenever a block has at least one secondary character -- and the length of a +#' tree under an asymmetric step matrix depends on where the tree is rooted, +#' unlike ordinary parsimony. The asymmetry is the point of the recoding (the +#' first gain of the controlling character pays for the secondaries it brings +#' into existence), so this is intrinsic rather than a defect. +#' +#' `TreeSearch` treats topologies as unrooted, so [`TreeLength()`] and +#' [`MaximizeParsimony()`] both evaluate x-transformation lengths at a canonical +#' rooting -- on the first taxon of `dataset` -- giving one length per topology +#' and making a reported score reproducible. That length is an upper bound on +#' the rooting-free minimum, exceeding it by at most the total number of +#' secondary characters across blocks. If a rooting is biologically meaningful +#' to you, score the tree yourself with the block's `cost_matrix` rather than +#' relying on the canonical value. +#' #' @param dataset A [`phyDat`][phangorn::phyDat] object. #' @param hierarchy A [`CharacterHierarchy`] object. #' diff --git a/R/tree_length.R b/R/tree_length.R index 856a3c874..54191f270 100644 --- a/R/tree_length.R +++ b/R/tree_length.R @@ -188,7 +188,20 @@ TreeLength.phylo <- function(tree, dataset, concavity = Inf, .BuildTipLabels(dataset), .HSJAbsentState(dataset)) } else if (useXform) { - tree <- RenumberTips(Renumber(tree), names(dataset)) + # The x-transformation's step matrix is asymmetric -- a gain costs `nSec + 1` + # against 1 for a loss -- so the Sankoff term below is rooting-dependent, + # while the rest of the pipeline treats topologies as unrooted and moves the + # root freely (Wagner addition, fusing, sector search and `ts_collapse_pool` + # all reroot). Score at a canonical rooting: the dataset's first taxon, which + # is the tip-0 rooting `ts_collapse_pool()` already imposes on the trees + # `MaximizeParsimony()` returns. One unrooted topology then has one length + # whatever rooting the user's `phylo` happens to carry, and this agrees with + # the score `MaximizeParsimony()` reports (T-374 / T-385). Root by NAME and + # re-align afterwards, so tip i still indexes `tip_data` row i. + # Decision and the measured `nSec`-per-block bound: + # dev/plans/2026-07-29-t374b-xform-rooting-policy.md. + tree <- RenumberTips(Renumber(RootTree(tree, names(dataset)[[1]])), + names(dataset)) at <- attributes(dataset) contrast <- at$contrast tip_data <- matrix(unlist(dataset, use.names = FALSE), @@ -296,6 +309,18 @@ TreeLength.list <- function(tree, dataset, concavity = Inf, needRoot <- !vapply(tree, TreeIsRooted, logical(1L)) if (any(needRoot)) warning("Unrooted tree rooted on tip 1.") tree[] <- lapply(tree, function(tr) if (TreeIsRooted(tr)) tr else RootTree(tr, 1)) + if (useXform) { + # XFORM's score is rooting-dependent (asymmetric step matrix), so an + # *already*-rooted tree must be canonicalised too, not just an unrooted one: + # otherwise the same topology gets different lengths from different rootings + # and no length agrees with what `MaximizeParsimony()` reports. See the + # single-tree method above for the full rationale. No warning here -- unlike + # the unrooted case, nothing is being assumed about the user's intent; the + # rooting simply is not part of this criterion's input. + rootTaxon <- names(dataset)[[1]] + tree[] <- lapply(tree, function(tr) RootTree(tr, rootTaxon)) + tree[] <- RenumberTips(tree, dataset) + } nEdge <- unique(vapply(tree, function(tr) dim(tr[["edge"]])[1], integer(1))) if (length(nEdge) > 1L) { diff --git a/dev/red-team/heavy-tests/t385-diagnose-rooting.R b/dev/red-team/heavy-tests/t385-diagnose-rooting.R new file mode 100644 index 000000000..8d0e4c2f4 --- /dev/null +++ b/dev/red-team/heavy-tests/t385-diagnose-rooting.R @@ -0,0 +1,102 @@ +# T-385 diagnostic: what rooting are the returned trees actually at, and do the +# pool members agree under a genuinely common one? +# +# Motivated by an apparent contradiction in t385-xform-report-agreement.R: +# all 32 returned trees scored 183 AS RETURNED, but scored 183/182 after +# RootTree(t, t$tip.label[1]). Both were described as "a common rooting", so at +# most one of those descriptions is right. Settle it before designing the fix. + +lib <- if (dir.exists(".agent-t385")) ".agent-t385" else NULL +suppressPackageStartupMessages({ + library("TreeSearch", lib.loc = lib) + library("TreeTools") +}) + +if (!file.exists("dev/red-team/heavy-tests/t385-make-xform-data.R")) { + stop("Run this from the package root: the source() path below is relative.") +} +source("dev/red-team/heavy-tests/t385-make-xform-data.R") + +dat <- MakeXformData() +ds <- dat$dataset +h <- dat$hierarchy + +set.seed(11) +res <- suppressWarnings( + MaximizeParsimony(ds, inapplicable = "xform", hierarchy = h, + maxReplicates = 4L, verbosity = 0L) +) +cat("reported:", unique(attr(res, "score")), "| n trees:", length(res), "\n\n") + +Score <- function(tr) TreeLength(tr, ds, inapplicable = "xform", hierarchy = h) + +# Is tip 1 (the dataset's first taxon == the kernel's tip 0) a child of the root? +RootedAtTip1 <- function(tr) { + root <- length(tr$tip.label) + 1L + tip1 <- match(names(ds)[1], tr$tip.label) + tip1 %in% tr$edge[tr$edge[, 1] == root, 2] +} + +# Does the root have exactly 2 children (rooted) or 3+ (unrooted-as-stored)? +RootDegree <- function(tr) { + root <- length(tr$tip.label) + 1L + sum(tr$edge[, 1] == root) +} + +info <- data.frame( + i = seq_along(res), + asReturned = vapply(res, Score, numeric(1)), + rootDeg = vapply(res, RootDegree, integer(1)), + atTip1 = vapply(res, RootedAtTip1, logical(1)), + viaRootTree = vapply(res, function(tr) Score(RootTree(tr, names(ds)[1])), + numeric(1)) +) +info$rtRootDeg <- vapply(res, function(tr) RootDegree(RootTree(tr, names(ds)[1])), + integer(1)) + +cat("--- per-tree ---\n") +print(head(info, 12)) +cat("...\n\n") + +cat("as-returned scores :", paste(sort(unique(info$asReturned)), + collapse = " "), "\n") +cat("via RootTree(., taxon 1) :", paste(sort(unique(info$viaRootTree)), + collapse = " "), "\n") +cat("root degree as returned :", paste(sort(unique(info$rootDeg)), + collapse = " "), "\n") +cat("root degree via RootTree :", paste(sort(unique(info$rtRootDeg)), + collapse = " "), "\n") +cat("tip1 is a root child :", sum(info$atTip1), "/", nrow(info), "\n\n") + +# The crux: are the trees whose score MOVES under RootTree the same ones that +# were not already rooted at taxon 1? +moved <- info$asReturned != info$viaRootTree +cat("--- trees whose score moves under RootTree ---\n") +cat("n moved:", sum(moved), "\n") +if (any(moved)) { + cat("of those, already rooted at tip1:", sum(info$atTip1[moved]), "\n") + cat("root degree of movers (as returned):", + paste(sort(unique(info$rootDeg[moved])), collapse = " "), "\n") + cat("deltas:", paste(unique(info$viaRootTree[moved] - info$asReturned[moved]), + collapse = " "), "\n") +} + +# Are the MOVED trees topologically distinct from the others, or the same +# topology stored differently? If RootTree changes the score, and the score is +# supposed to depend only on the unrooted topology, then either the topologies +# differ or the scorer is rooting-sensitive (which is the finding). +cat("\n--- topology check ---\n") +# Canonical key: root every tree at the same taxon, sort, serialise. (Avoids +# as.Splits comparisons -- see the Splits/%in% dispatch trap.) +TopoKey <- function(tr) { + ape::write.tree(SortTree(RootTree(tr, names(ds)[1]))) +} +keys <- vapply(res, TopoKey, character(1)) +cat("distinct topologies among returned trees:", length(unique(keys)), + "of", length(res), "\n") +# Do trees sharing one topology share a score? +byTopo <- tapply(info$viaRootTree, keys, function(x) length(unique(x))) +cat("topologies whose duplicate copies disagree on score:", + sum(byTopo > 1L), "\n") +cat("score range within the pool at a common rooting:", + paste(range(info$viaRootTree), collapse = " - "), "\n") diff --git a/dev/red-team/heavy-tests/t385-find-small-case.R b/dev/red-team/heavy-tests/t385-find-small-case.R new file mode 100644 index 000000000..9c0a4dc66 --- /dev/null +++ b/dev/red-team/heavy-tests/t385-find-small-case.R @@ -0,0 +1,63 @@ +# Find a SMALL xform dataset + tree whose Sankoff term genuinely varies with the +# rooting, for use as a self-guarding regression test. Measured with +# ts_sankoff_test() directly -- the raw kernel, which the T-385 fix does not +# touch -- so the sensitivity is a property of the data, not of TreeLength(). + +lib <- if (dir.exists(".agent-t385")) ".agent-t385" else NULL +suppressPackageStartupMessages({ + library("TreeSearch", lib.loc = lib) + library("TreeTools") +}) + +KernelSpread <- function(ds, h, tree) { + recoded <- RecodeHierarchy(ds, h) + xf <- TreeSearch:::.PrepareXformArgs(recoded, length(ds)) + vals <- vapply(names(ds), function(taxon) { + tr <- RenumberTips(Renumber(RootTree(tree, taxon)), names(ds)) + TreeSearch:::ts_sankoff_test(tr[["edge"]], xf$n_states, xf$cost_matrices, + xf$tip_states, xf$forced_root, xf$combo_grids, + xf$tip_sec_known)$score + }, numeric(1)) + vals +} + +best <- NULL +for (seed in 1:400) { + set.seed(seed) + nTip <- 8L + tipNames <- paste0("t", seq_len(nTip)) + primary <- sample(c("0", "1"), nTip, replace = TRUE) + if (length(unique(primary)) < 2L) next + sec1 <- ifelse(primary == "0", "-", sample(c("0", "1"), nTip, replace = TRUE)) + sec2 <- ifelse(primary == "0", "-", sample(c("0", "1"), nTip, replace = TRUE)) + mat <- cbind(primary, sec1, sec2) + dimnames(mat) <- list(tipNames, NULL) + ds <- try(phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"), + ambiguity = "?"), silent = TRUE) + if (inherits(ds, "try-error")) next + h <- CharacterHierarchy("1" = 2:3) + tree <- RandomTree(tipNames, root = TRUE) + vals <- try(KernelSpread(ds, h, tree), silent = TRUE) + if (inherits(vals, "try-error")) next + sp <- diff(range(vals)) + if (sp > 0) { + cat(sprintf("seed %3d: spread %g values %s\n", seed, sp, + paste(vals, collapse = " "))) + if (is.null(best) || sp > best$spread) { + best <- list(seed = seed, spread = sp, mat = mat, vals = vals, + newick = ape::write.tree(tree)) + } + if (sp >= 2) break + } +} + +if (is.null(best)) { + cat("no rooting-sensitive 8-tip case found\n") +} else { + cat("\n=== BEST ===\nseed:", best$seed, " spread:", best$spread, "\n") + cat("newick:", best$newick, "\n") + cat("matrix:\n") + print(best$mat) + cat("\nvalues by root taxon:\n") + print(best$vals) +} diff --git a/dev/red-team/heavy-tests/t385-make-xform-data.R b/dev/red-team/heavy-tests/t385-make-xform-data.R new file mode 100644 index 000000000..f5a860127 --- /dev/null +++ b/dev/red-team/heavy-tests/t385-make-xform-data.R @@ -0,0 +1,41 @@ +# Shared synthetic XFORM dataset generator for the T-385 scripts. +# +# nBlock hierarchy blocks, each a controlling primary + nSec secondaries, plus +# some free (non-hierarchy) characters. The gain/loss asymmetry that makes XFORM +# rooting-sensitive needs nSec >= 1; the free characters carry no gaps, so +# `has_na` is driven only by the hierarchy blocks (the has_na = FALSE branch of +# the T-374b analysis, i.e. the Sankoff term stays out of the convergence sweep). + +MakeXformData <- function(nTip = 36L, nBlock = 6L, nSec = 2L, nFree = 10L, + seed = 1L) { + set.seed(seed) + tipNames <- paste0("t", seq_len(nTip)) + cols <- list() + hierArgs <- list() + + for (b in seq_len(nBlock)) { + primary <- sample(c("0", "1"), nTip, replace = TRUE) + priIdx <- length(cols) + 1L + cols[[length(cols) + 1L]] <- primary + secIdx <- integer(nSec) + for (s in seq_len(nSec)) { + # Secondary is inapplicable exactly where the primary is absent. + sec <- ifelse(primary == "0", "-", + sample(c("0", "1"), nTip, replace = TRUE)) + cols[[length(cols) + 1L]] <- sec + secIdx[s] <- length(cols) + } + hierArgs[[as.character(priIdx)]] <- secIdx + } + + for (f in seq_len(nFree)) { + cols[[length(cols) + 1L]] <- sample(c("0", "1"), nTip, replace = TRUE) + } + + mat <- do.call(cbind, cols) + dimnames(mat) <- list(tipNames, NULL) + ds <- phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"), + ambiguity = "?") + list(dataset = ds, hierarchy = do.call(CharacterHierarchy, hierArgs), + nSec = nSec, nBlock = nBlock) +} diff --git a/dev/red-team/heavy-tests/t385-xform-report-agreement.R b/dev/red-team/heavy-tests/t385-xform-report-agreement.R new file mode 100644 index 000000000..8f9e4d431 --- /dev/null +++ b/dev/red-team/heavy-tests/t385-xform-report-agreement.R @@ -0,0 +1,122 @@ +# T-385 (P1) — does MaximizeParsimony's reported XFORM score agree with +# TreeLength() on the tree it just returned? +# +# The finding (2026-07-31, found while building the T-378 benchmark panel): +# attr(res, "score") disagreed with TreeLength(res[[1]], ...) by 10 steps on a +# 36-tip synthetic XFORM dataset, and rerooting the *same* returned topology +# gave a third value again (386 reported / 396 as-returned / 394 rerooted). +# +# This script reproduces that against whatever tip it is run on, BEFORE any fix +# is written, because the tip has moved since the finding was recorded. It +# reports rather than asserts: it prints a verdict line and exits 0 when the +# three quantities agree, 1 when they do not. +# +# Run: Rscript dev/red-team/heavy-tests/t385-xform-report-agreement.R +# (expects the package in .agent-t385, i.e. the AGENTS.md tarball build) + +lib <- if (dir.exists(".agent-t385")) ".agent-t385" else NULL +suppressPackageStartupMessages({ + library("TreeSearch", lib.loc = lib) + library("TreeTools") +}) + +cat("TreeSearch:", as.character(packageVersion("TreeSearch")), "\n") + +# ---- A synthetic XFORM dataset ----------------------------------------------- +# Generator shared with t385-diagnose-rooting.R so the two cannot drift apart. +if (!file.exists("dev/red-team/heavy-tests/t385-make-xform-data.R")) { + stop("Run this from the package root: the source() path below is relative.") +} +source("dev/red-team/heavy-tests/t385-make-xform-data.R") + +dat <- MakeXformData() +ds <- dat$dataset +h <- dat$hierarchy +cat(sprintf("dataset: %d tips, %d chars, %d blocks x %d secondaries\n", + length(ds), attr(ds, "nr"), dat$nBlock, dat$nSec)) +cat(sprintf("bound on rooting overstatement (sum nSec per block): %d\n", + dat$nBlock * dat$nSec)) + +# ---- Search ------------------------------------------------------------------ +set.seed(11) +res <- suppressWarnings( + MaximizeParsimony(ds, inapplicable = "xform", hierarchy = h, + maxReplicates = 4L, verbosity = 0L) +) + +reported <- attr(res, "score") +cat("\n--- reported ---\n") +cat("attr(res, 'score'):", paste(unique(reported), collapse = " "), "\n") +cat("n trees returned:", length(res), "\n") + +# ---- Symptom 1: reported score vs TreeLength() of the returned tree --------- +asReturned <- vapply(seq_along(res), function(i) { + TreeLength(res[[i]], ds, inapplicable = "xform", hierarchy = h) +}, numeric(1)) + +cat("\n--- symptom 1: reported vs TreeLength(returned tree) ---\n") +cat("TreeLength as-returned:", paste(round(asReturned, 4), collapse = " "), "\n") +gap1 <- max(abs(asReturned - min(reported))) +cat(sprintf("max |TreeLength - reported| = %g\n", gap1)) + +# ---- Symptom 2: rerooting the SAME topology changes the score --------------- +# The engine's rooting is an artefact (RandomTree/ts_collapse_pool tip-0), so a +# score that moves under rerooting is the un-pinned reporting gap. +tree1 <- res[[1]] +rootings <- vapply(tree1$tip.label[seq_len(min(8L, length(tree1$tip.label)))], + function(tip) { + TreeLength(RootTree(tree1, tip), ds, + inapplicable = "xform", hierarchy = h) + }, numeric(1)) + +cat("\n--- symptom 2: same topology, different rootings ---\n") +print(round(rootings, 4)) +spread <- diff(range(rootings)) +cat(sprintf("spread across %d rootings = %g (documented bound = %d)\n", + length(rootings), spread, dat$nBlock * dat$nSec)) + +# ---- Symptom 3: is the returned MPT set internally consistent? ------------- +cat("\n--- symptom 3: MPT-set self-consistency ---\n") +if (length(res) > 1L) { + # Root every tree at the SAME NAMED taxon. Do NOT use res[[i]]$tip.label[1]: + # after Renumber() that is a different taxon for different trees, so it scores + # each tree at a different rooting and manufactures a spurious disagreement. + # (That error produced a retracted "2 distinct scores" claim on 2026-07-31.) + rootTaxon <- names(ds)[1] + canon <- vapply(seq_along(res), function(i) { + TreeLength(RootTree(res[[i]], rootTaxon), ds, + inapplicable = "xform", hierarchy = h) + }, numeric(1)) + cat("scores at a common (tip-1) rooting:", + paste(round(canon, 4), collapse = " "), "\n") + nDistinct <- length(unique(round(canon, 8))) + cat(sprintf("distinct scores among %d 'equally parsimonious' trees: %d\n", + length(canon), nDistinct)) +} else { + canon <- asReturned + nDistinct <- 1L + cat("single tree returned; not informative\n") +} + +# ---- Verdict ----------------------------------------------------------------- +tol <- 1e-8 +agree1 <- gap1 < tol +agree2 <- spread < tol +agree3 <- nDistinct == 1L + +cat("\n===== VERDICT =====\n") +cat(sprintf("reported == TreeLength(returned) : %s\n", + if (agree1) "AGREE" else sprintf("DISAGREE by %g", gap1))) +cat(sprintf("score invariant under rerooting : %s\n", + if (agree2) "INVARIANT" else sprintf("VARIES, spread %g", spread))) +cat(sprintf("MPT set shares one score : %s\n", + if (agree3) "CONSISTENT" + else sprintf("INCONSISTENT, %d distinct", nDistinct))) + +if (agree1 && agree2 && agree3) { + cat("\nT-385 does NOT reproduce on this tip.\n") + quit(status = 0) +} else { + cat("\nT-385 REPRODUCES on this tip.\n") + quit(status = 1) +} diff --git a/inst/WORDLIST b/inst/WORDLIST index e2d9eb350..940b38b7e 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -191,6 +191,7 @@ XPIWE XSS YU ZANOL +Zanol Zoologica abcd ac @@ -259,6 +260,7 @@ preorder preset's propto pscore +rootings rRNA rearranger reconverged diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 7fb1e5a75..4f458e822 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -114,7 +114,18 @@ hierarchy required.} \item{\code{"xform"}}{Step-matrix recoding approximating maximum homology via x-transformations \insertCite{Goloboff2021;textual}{TreeSearch}. Requires a -\code{hierarchy}.} +\code{hierarchy}. \strong{Scores are rooting-sensitive}: the step matrix of +this recoding is asymmetric -- gaining the controlling character costs one more +than the number of secondaries it brings into existence, against 1 to +lose it -- so a tree's length depends on where its root sits, whereas +parsimony under the other methods does not. Lengths are therefore +reported at a canonical rooting, on the first taxon of \code{dataset}, which +is the rooting the returned trees carry; \code{TreeLength()} canonicalises +identically, so it reproduces the reported score and one topology has +one length. That value is an upper bound on the rooting-free minimum, +exceeding it by at most the total number of secondary characters across +hierarchy blocks, and attaining it for 87-98\% of rootings in +simulation.} }} \item{hsj_alpha}{Numeric in [0, 1]: scaling parameter for secondary- @@ -317,6 +328,14 @@ collapse, so an enforced-but-unsupported clade (a zero-length branch) stays visible (a constraint encodes external evidence the matrix does not capture); unsupported non-constraint branches still collapse.} +\item{.rung}{Internal. Pins a named entry of the effort ladder +(\code{"sprint"}, \code{"default"}, \code{"thorough"}, \code{"large"}), or \code{"none"} to apply no +preset at all; \code{NULL} (default) selects the rung from \code{effort} and the +dataset's size, which is what every ordinary call should do. Exists for +controlled experiments and the preset smoke tests, which need to name a rung +absolutely rather than relative to the automatic choice. Not part of the +stable interface: prefer \code{effort}.} + \item{...}{Backward compatibility.} } \value{ diff --git a/man/RecodeHierarchy.Rd b/man/RecodeHierarchy.Rd index e8b8e8e7d..91e24135a 100644 --- a/man/RecodeHierarchy.Rd +++ b/man/RecodeHierarchy.Rd @@ -65,6 +65,25 @@ number of secondary characters. different states). } } + +\subsection{Rooting}{ + +Because gain and loss cost differently, this cost matrix is \strong{asymmetric} +whenever a block has at least one secondary character -- and the length of a +tree under an asymmetric step matrix depends on where the tree is rooted, +unlike ordinary parsimony. The asymmetry is the point of the recoding (the +first gain of the controlling character pays for the secondaries it brings +into existence), so this is intrinsic rather than a defect. + +\code{TreeSearch} treats topologies as unrooted, so \code{\link[=TreeLength]{TreeLength()}} and +\code{\link[=MaximizeParsimony]{MaximizeParsimony()}} both evaluate x-transformation lengths at a canonical +rooting -- on the first taxon of \code{dataset} -- giving one length per topology +and making a reported score reproducible. That length is an upper bound on +the rooting-free minimum, exceeding it by at most the total number of +secondary characters across blocks. If a rooting is biologically meaningful +to you, score the tree yourself with the block's \code{cost_matrix} rather than +relying on the canonical value. +} } \references{ \insertAllCited{} diff --git a/man/Resample.Rd b/man/Resample.Rd index becf48e81..e99630af0 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -124,7 +124,18 @@ hierarchy required.} \item{\code{"xform"}}{Step-matrix recoding approximating maximum homology via x-transformations \insertCite{Goloboff2021;textual}{TreeSearch}. Requires a -\code{hierarchy}.} +\code{hierarchy}. \strong{Scores are rooting-sensitive}: the step matrix of +this recoding is asymmetric -- gaining the controlling character costs one more +than the number of secondaries it brings into existence, against 1 to +lose it -- so a tree's length depends on where its root sits, whereas +parsimony under the other methods does not. Lengths are therefore +reported at a canonical rooting, on the first taxon of \code{dataset}, which +is the rooting the returned trees carry; \code{TreeLength()} canonicalises +identically, so it reproduces the reported score and one topology has +one length. That value is an upper bound on the rooting-free minimum, +exceeding it by at most the total number of secondary characters across +hierarchy blocks, and attaining it for 87-98\% of rootings in +simulation.} }} \item{hsj_alpha}{Numeric in [0, 1]: scaling parameter for secondary- diff --git a/man/TreeLength.Rd b/man/TreeLength.Rd index 4f76a73fb..767c11c9c 100644 --- a/man/TreeLength.Rd +++ b/man/TreeLength.Rd @@ -144,7 +144,18 @@ hierarchy required.} \item{\code{"xform"}}{Step-matrix recoding approximating maximum homology via x-transformations \insertCite{Goloboff2021;textual}{TreeSearch}. Requires a -\code{hierarchy}.} +\code{hierarchy}. \strong{Scores are rooting-sensitive}: the step matrix of +this recoding is asymmetric -- gaining the controlling character costs one more +than the number of secondaries it brings into existence, against 1 to +lose it -- so a tree's length depends on where its root sits, whereas +parsimony under the other methods does not. Lengths are therefore +reported at a canonical rooting, on the first taxon of \code{dataset}, which +is the rooting the returned trees carry; \code{TreeLength()} canonicalises +identically, so it reproduces the reported score and one topology has +one length. That value is an upper bound on the rooting-free minimum, +exceeding it by at most the total number of secondary characters across +hierarchy blocks, and attaining it for 87-98\% of rootings in +simulation.} }} \item{hsj_alpha}{Numeric in [0, 1]: scaling parameter for secondary- diff --git a/man/expected_mi.Rd b/man/expected_mi.Rd index 4adbde79f..9b455127e 100644 --- a/man/expected_mi.Rd +++ b/man/expected_mi.Rd @@ -35,9 +35,6 @@ the hypergeometric distribution of cell overlaps, and is returned in bits # Expected MI between a 3|4 split and a 2|5 split of 7 items: expected_mi(c(3L, 4L), c(2L, 5L)) } -\references{ -\insertAllCited{} -} \seealso{ \code{\link{SiteConcordance}} } diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index 05f25c6c3..3e02b6e00 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -47,7 +47,6 @@ inline double l2factorial(int n) { //' @param nj Integer vector giving the block sizes of the second partition //' (also summing to `N`). //' @return The expected mutual information, in bits. -//' @references \insertAllCited{} //' @seealso [`SiteConcordance`] //' @examples //' # Expected MI between a 3|4 split and a 2|5 split of 7 items: diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 0f53b8f94..914314b1d 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2188,8 +2188,13 @@ List ts_collapse_pool( ts::TreeState tree; tree.init_from_edge(&edge(0, 0), &edge(0, 1), edge.nrow(), ds); - // Root on tip 0 so root-adjacent edges are trivial (rooting-invariant - // collapse), then refresh state arrays for the flag computation. + // Root on tip 0 so root-adjacent edges are trivial, then refresh state + // arrays for the flag computation. This makes the CONTRACTION + // rooting-invariant, and it is also what fixes the rooting the returned + // trees are handed back at — which matters because HSJ/XFORM lengths are + // NOT rooting-invariant (T-374). MaximizeParsimony() rescores the XFORM + // pool at this same tip-0 rooting before reporting, so that the reported + // score is the score of the tree returned (T-385). ts::reroot_at_tip(tree, 0); tree.reset_states(ds); ts::score_tree(tree, ds); diff --git a/src/ts_tbr.cpp b/src/ts_tbr.cpp index a1eab3e9e..e03edfe62 100644 --- a/src/ts_tbr.cpp +++ b/src/ts_tbr.cpp @@ -120,9 +120,15 @@ static double full_rescore(TreeState& tree, const DataSet& ds) { } // Re-root the tree so tip `t` is a direct child of the root pseudo-node n_tip. -// Parsimony length is root-invariant, so this only changes the representation -// (which edges are clippable and where the root edge sits) — it lets the search -// reach moves the current rooting hides. Rebuilds postorder; does NOT refresh +// Parsimony length is root-invariant for the SYMMETRIC criteria (EW / IW / NA / +// profile), so for those this only changes the representation (which edges are +// clippable and where the root edge sits) — it lets the search reach moves the +// current rooting hides. It is NOT root-invariant under HSJ or XFORM: the +// x-transformation's step matrix is asymmetric (gain = nSec + 1 against loss = 1) +// and HSJ's alpha.d/m term is evaluated on a directional pick, so a reroot can +// change their scores by up to sum(nSec) over blocks (T-374; measured in +// dev/plans/2026-07-29-t374b-xform-rooting-policy.md). Callers under those modes +// must treat the score as rooting-relative. Rebuilds postorder; does NOT refresh // Fitch state arrays, so the caller must full_rescore() afterwards. // Generalises reroot_at_tip0() in ts_fuse.cpp to an arbitrary tip. // Declared in ts_tbr.h (used by the output-collapse kernel in ts_rcpp.cpp). @@ -3149,7 +3155,10 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, reroot_prev = best_score; reroot_at_tip(tree, reroot_tip); reroot_tip = (reroot_tip + 1) % tree.n_tip; - best_score = full_rescore(tree, ds); // root-invariant; refreshes states + // Refreshes states. Root-invariant for EW / IW / NA / profile only — under + // HSJ / XFORM the score is rooting-relative (T-374), so this value is the + // length at the rooting the tree currently carries, not an absolute. + best_score = full_rescore(tree, ds); score_fresh = true; if (!collapsed.empty()) { if (collapse_aggr) compute_collapsed_flags_aggressive(tree, ds, collapsed); diff --git a/tests/testthat/test-stop-patience.R b/tests/testthat/test-stop-patience.R index 95465510c..4024995c9 100644 --- a/tests/testthat/test-stop-patience.R +++ b/tests/testthat/test-stop-patience.R @@ -124,11 +124,21 @@ test_that("stopPatience also stops the parallel search", { expect_equal(attr(full, "replicates"), cap) # nothing else ends the search expect_lt(attr(short, "replicates"), cap) expect_true(attr(short, "perturb_stop")) - # NOT expect_equal: the 200 ms monitor poll makes the firing replicate (and thus which - # dry spell trips the rule) vary run to run, so demanding the SAME optimum from a - # 5-9-replicate stop as from the full 120-replicate run is a search-quality claim - # conditional on timing, not a property of the stopping rule. What the rule actually + # Deliberately NOT `expect_equal(short score, full score)`. That assertion held on + # this matrix but is not a property the rule guarantees, and it broke under covr: + # instrumentation slows every replicate, the parallel path evaluates the dry spell on a + # 200 ms monitor poll over replicates completed into the shared pool, so far fewer + # replicates land per poll and patience 3 fires genuinely earlier in the search + # (13.063 against 12.988). Stopping sooner is *allowed* to cost score -- that is the + # trade-off the parameter exists to offer -- so equality was asserting a coincidence of + # this machine's timing. What the rule does guarantee is asserted above: the search + # stops before the cap, and it stops for the no-improvement reason. The score is only + # checked for being a valid, finite improvement over a random start, which holds at any + # speed. See [[parallel-stop-rules-poll-granularity]]. + # What the rule actually # guarantees is that stopping early can't do BETTER than letting the search run on. + expect_true(is.finite(min(attr(short, "score")))) + expect_lt(min(attr(short, "score")), min(attr(full, "score")) * 1.5) expect_gte(min(attr(short, "score")), min(attr(full, "score"))) }) diff --git a/tests/testthat/test-ts-xform.R b/tests/testthat/test-ts-xform.R index da9f97d50..f5606918d 100644 --- a/tests/testthat/test-ts-xform.R +++ b/tests/testthat/test-ts-xform.R @@ -397,8 +397,19 @@ test_that("Xform search handles all-hierarchy data (zero Fitch words)", { expect_length(setdiff(seq_len(5L), HierarchyChars(h)), 0L) set.seed(42) - res <- MaximizeParsimony(ds, hierarchy = h, inapplicable = "xform", - maxReplicates = 4L, targetHits = 3L, verbosity = 0L) + # This all-hierarchy matrix provably triggers T-374's open residue: pool + # membership is decided on search-time scores taken at differing rootings, so + # the returned trees do not share a length at the common rooting they are + # reported at (measured: 7 to 9), and MaximizeParsimony() warns and reports the + # smallest. Pinned as an expectation rather than left as ambient noise in a + # green suite -- if this warning ever STOPS firing, the residue has been fixed + # (or masked) and that deserves to be noticed here. + res <- NULL + expect_warning( + res <- MaximizeParsimony(ds, hierarchy = h, inapplicable = "xform", + maxReplicates = 4L, targetHits = 3L, + verbosity = 0L), + "do not share a length") expect_s3_class(res[[1]], "phylo") expect_equal(length(res[[1]]$tip.label), 6L) for (tr in res) { @@ -406,3 +417,115 @@ test_that("Xform search handles all-hierarchy data (zero Fitch words)", { expect_true(TreeIsRooted(tr)) } }) + + +# ===== T-385: a reported length must be reproducible ========================== +# The x-transformation's step matrix is asymmetric (gain = nSec + 1, loss = 1), +# so a tree's length depends on where its root sits -- unlike parsimony under the +# symmetric criteria. The engine recorded `best_score` mid-search at whatever +# rooting the replicate held, while `ts_collapse_pool()` returns every tree +# re-rooted on tip 0, so `attr(res, "score")` did not reproduce under +# `TreeLength()` of the very tree returned (measured: 178 reported against 183 +# returned, 36 tips / 6 blocks). Both boundaries now score at a canonical +# rooting. See dev/plans/2026-07-29-t374b-xform-rooting-policy.md. + +test_that("TreeLength xform is rooting-invariant (T-385)", { + # 8-tip case whose Sankoff term genuinely varies with the rooting. That + # precondition is ASSERTED below rather than assumed, so this test cannot pass + # vacuously on data that happens to be rooting-insensitive. + mat <- matrix(c( + "0", "-", "-", + "1", "1", "1", + "0", "-", "-", + "0", "-", "-", + "1", "0", "0", + "0", "-", "-", + "0", "-", "-", + "0", "-", "-" + ), nrow = 8, byrow = TRUE, + dimnames = list(paste0("t", 1:8), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("1" = 2:3) + tree <- ape::read.tree(text = "(((t1,t3),((t2,t5),t7)),(t4,(t6,t8)));") + taxa <- names(ds) + + # Precondition: the raw Sankoff kernel IS rooting-sensitive here (3 to 5). + # ts_sankoff_test() is untouched by the fix, so this measures the data. + recoded <- RecodeHierarchy(ds, h) + xf <- TreeSearch:::.PrepareXformArgs(recoded, length(ds)) + kernel <- vapply(taxa, function(taxon) { + tr <- RenumberTips(Renumber(RootTree(tree, taxon)), taxa) + TreeSearch:::ts_sankoff_test(tr[["edge"]], xf$n_states, xf$cost_matrices, + xf$tip_states, xf$forced_root, xf$combo_grids, + xf$tip_sec_known)$score + }, numeric(1)) + expect_gt(diff(range(kernel)), 0) + + # Given that, TreeLength() must still return ONE length for ONE topology. + lengths <- vapply(taxa, function(taxon) { + TreeLength(RootTree(tree, taxon), ds, hierarchy = h, + inapplicable = "xform") + }, numeric(1)) + expect_equal(diff(range(lengths)), 0) + + # The multiPhylo method must canonicalise identically to the single-tree one: + # it previously rooted only trees that arrived unrooted. + multi <- TreeLength( + structure(lapply(taxa, function(taxon) RootTree(tree, taxon)), + class = "multiPhylo"), + ds, hierarchy = h, inapplicable = "xform" + ) + expect_equal(unname(multi), unname(lengths)) +}) + +test_that("MaximizeParsimony xform reports the length of the tree it returns (T-385)", { + # End-to-end guard on the report path. Larger than the test above so the + # search has room to end at a rooting other than the returned one. + set.seed(7) + nTip <- 16L + cols <- list() + hierArgs <- list() + for (b in 1:3) { + primary <- sample(c("0", "1"), nTip, replace = TRUE) + priIdx <- length(cols) + 1L + cols[[length(cols) + 1L]] <- primary + secIdx <- integer(2) + for (s in 1:2) { + cols[[length(cols) + 1L]] <- ifelse( + primary == "0", "-", sample(c("0", "1"), nTip, replace = TRUE)) + secIdx[s] <- length(cols) + } + hierArgs[[as.character(priIdx)]] <- secIdx + } + for (f in 1:6) { + cols[[length(cols) + 1L]] <- sample(c("0", "1"), nTip, replace = TRUE) + } + mat <- do.call(cbind, cols) + dimnames(mat) <- list(paste0("t", seq_len(nTip)), NULL) + ds <- make_dat(mat) + h <- do.call(CharacterHierarchy, hierArgs) + + set.seed(11) + res <- suppressWarnings( + MaximizeParsimony(ds, inapplicable = "xform", hierarchy = h, + maxReplicates = 3L, verbosity = 0L)) + + reported <- unique(unname(attr(res, "score"))) + expect_length(reported, 1L) + returned <- unname(vapply(res, function(tr) { + TreeLength(tr, ds, hierarchy = h, inapplicable = "xform") + }, numeric(1))) + + # The contract is that the reported score is the canonical length of the pool, + # reproducible from the trees returned -- NOT that every returned tree shares + # it. Pool membership is still decided on search-time scores taken at differing + # rootings, so on some data the pool genuinely spans a range and + # MaximizeParsimony() warns (all-hierarchy matrices do: see the zero-Fitch-words + # test above, which spans 7 to 9). That is the open residue of T-374 and is + # deliberately out of scope for the reporting fix. Assert the contract: + expect_equal(reported, min(returned)) + + # On THIS dataset the pool does happen to be self-consistent. Asserted as a + # property of the data, not as something the code promises. + expect_equal(diff(range(returned)), 0) +}) diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index 6e9af41eb..d77b371f5 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -269,6 +269,34 @@ program or Sankoff contribution) that an incremental Fitch delta cannot capture; for these modes an accepted move is therefore re-scored in full before the best score is updated. +### Reporting an x-transformation score + +The x-transformation's step matrix is asymmetric: gaining the controlling +character costs one more than the number of secondary characters it brings +into existence, whereas losing it costs one. +A tree's length under an asymmetric step matrix depends on where the tree is +rooted, whereas parsimony under the symmetric criteria does not -- and the +search moves the root freely, since Wagner addition, tree fusing, sector search +and the output collapse all re-root. + +The score the search compares is therefore the length at whichever rooting a +replicate happens to hold. +That is a sound quantity to search on: it is an upper bound on the rooting-free +minimum, and in simulation 87--98% of rootings attain that minimum exactly, with +the excess never exceeding the total number of secondary characters across +hierarchy blocks. +It is not, however, reproducible from the returned tree, whose rooting the +output stage has since canonicalised. +`MaximizeParsimony()` therefore re-scores the returned trees at that same +canonical rooting -- the first taxon of the dataset -- before reporting, and +`TreeLength()` canonicalises identically, so a reported x-transformation score +is the length of the tree in hand, and one topology has one length. +This affects reporting only: the quantity the search optimises is unchanged. + +Because pool membership is still decided on scores taken at differing rootings, +the returned trees need not all share that canonical length; +`MaximizeParsimony()` warns when they do not. + ### Zero-length edge skipping Following @Goloboff1996, TBR, SPR, and drift skip provably unproductive