T-385: report the x-transformation score of the tree actually returned

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

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,

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:

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.
 #'

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) {

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")

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)
+}