From a0939a856bbb4e034ecd052a6d2309203427accb Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Tue, 4 Aug 2026 16:56:19 +0100
Subject: [PATCH 1/6] fix: derive XFORM secondary state spaces from contrasts,
 not token strings

A hierarchy block's secondary levels were the distinct token STRINGS observed
in that column, minus "-" and "?".  An ambiguity token such as "{01}" survived
that filter and became a level of its own: one Hamming step from both "0" and
"1" rather than matching either, so a polymorphic cell cost a step no
resolution of it needs, and each such cell multiplied the combination count
((k + 1)^m present states rather than k^m).  Read the levels from the phyDat
contrast matrix instead, as the HSJ path already does, and carry an
unresolved secondary as a bit mask of admissible levels so an ambiguity
narrower than the full state space still constrains the combination.

A secondary with no observed level -- reachable by dropping taxa from a
dataset that validated, since ValidateHierarchy runs before both subsetting
sites -- gave a zero-width state space, so no tip admitted any state and the
block's length was Inf.  Carry it as one unobserved level: it adds nothing to
any tree, and the gain cost still counts every secondary the primary controls.

MaximizeParsimony() rescores its own returned pool at the canonical rooting.
That call could not survive a polytomous pool (TreeLength scores the topology
it is given, and a contracted tree's length is that of its best resolution) or
a pool with no finite length (diff(range(.)) is NaN, which `if` aborts on).
Score the binary pool the returned trees were contracted from, and reduce the
lengths through a guard that reports rather than branches on non-finite ones.

Fixes #11
Fixes #12
Fixes #17

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 R/MaximizeParsimony.R                     |  89 +++++++--
 R/recode_hierarchy.R                      | 101 +++++++---
 man/RecodeHierarchy.Rd                    |  19 +-
 src/ts_rcpp.cpp                           |  26 +--
 tests/testthat/test-ts-xform-statespace.R | 222 ++++++++++++++++++++++
 5 files changed, 404 insertions(+), 53 deletions(-)
 create mode 100644 tests/testthat/test-ts-xform-statespace.R

diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R
index 3ccce1f44..1fe339d70 100644
--- a/R/MaximizeParsimony.R
+++ b/R/MaximizeParsimony.R
@@ -1709,22 +1709,18 @@ MaximizeParsimony <- function(
   # 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
+    # The second argument is evaluated only if `outTrees` turns out not to be
+    # binary, so the ordinary path builds nothing extra.
+    bestScore <- .XformPoolScore(
+      outTrees,
+      lapply(resultTrees[result$scores == result$best_score],
+             function(edgeMat) {
+               tr <- treeTpl
+               tr[["edge"]] <- edgeMat
+               Renumber(tr)
+             }),
+      dataset, hierarchy, result$best_score
     )
-    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 ---
@@ -1770,6 +1766,69 @@ MaximizeParsimony <- function(
   )
 }
 
+# Reduce a returned pool's canonical-rooting lengths to the single score
+# `MaximizeParsimony()` reports under `inapplicable = "xform"`, falling back to
+# `fallback` (the search's own best score) where the pool cannot supply one.
+#
+# `collapse = TRUE` contracts unsupported branches, and `TreeLength()` scores
+# the topology it is given -- the length of a polytomy is that of its best
+# resolution, which the Sankoff kernel does not compute, so a contracted pool
+# reports a number that is too small (T-401).  Nothing but the HSJ/XFORM no-op
+# in `compute_collapsed_flags_aggressive()` currently keeps `outTrees` binary,
+# and this must not depend on that staying in place: score instead the binary
+# pool the trees were contracted from, whose lengths are the same because only
+# zero-length branches are removed.  Kept separate from its caller because a
+# default search cannot reach that path, so this is the only place it can be
+# exercised.
+.XformPoolScore <- function(pool, binaryPool, dataset, hierarchy, fallback) {
+  nEdgeBinary <- 2L * length(dataset) - 2L
+  .Binary <- function(trees) {
+    vapply(trees, function(tr) dim(tr[["edge"]])[[1]], integer(1)) == nEdgeBinary
+  }
+  if (!all(.Binary(pool))) {
+    pool <- binaryPool[.Binary(binaryPool)]
+  }
+  if (length(pool) == 0L) {
+    return(fallback)
+  }
+  .ReportXformScore(
+    TreeLength(structure(pool, class = "multiPhylo"), dataset,
+               inapplicable = "xform", hierarchy = hierarchy),
+    fallback)
+}
+
+# Reduce a pool's canonical-rooting lengths to the one number reported.
+# `fallback` covers a pool with no finite length: `diff(range(.))` is then
+# `NaN`, which `if` cannot branch on -- the abort a degenerate hierarchy block
+# used to produce (T-394).
+.ReportXformScore <- function(canonicalScores, fallback) {
+  finite <- is.finite(canonicalScores)
+  if (!all(finite)) {
+    warning(sum(!finite), " of ", length(canonicalScores),
+            " returned trees have no ",
+            "finite x-transformation length; ",
+            if (any(finite)) "reporting the shortest of the rest."
+            else "reporting the search's own score.")
+  }
+  if (!any(finite)) {
+    return(fallback)
+  }
+  canonicalScores <- canonicalScores[finite]
+
+  if (diff(range(canonicalScores)) > sqrt(.Machine$double.eps)) {
+    # Pool membership is chosen on search-time scores taken at differing
+    # rootings (`result$scores` in the caller), 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.")
+  }
+  min(canonicalScores)
+}
+
 #' Launch tree search graphical user interface
 #'
 #' Opens a "shiny" app for interactive parsimony tree search and results
diff --git a/R/recode_hierarchy.R b/R/recode_hierarchy.R
index 2cc435a50..f6c43c470 100644
--- a/R/recode_hierarchy.R
+++ b/R/recode_hierarchy.R
@@ -14,6 +14,18 @@
 #' secondary character states (where \eqn{k_i} is the number of informative
 #' states of secondary character \eqn{i}).
 #'
+#' The informative levels of a secondary character are read from the dataset's
+#' `contrast` matrix, not from the token strings it carries.  An ambiguity token
+#' such as `"{01}"` therefore denotes the *set* of states it contrasts against,
+#' as it does under `inapplicable = "hsj"`, rather than becoming a level of its
+#' own; a token that admits every applicable state (`"?"`, or an ambiguity
+#' spanning them all) shows that no state is established and so contributes
+#' none.  A secondary whose levels are all unobserved -- reachable by dropping
+#' taxa from a dataset that validated -- is carried as a single unobserved
+#' level: it adds nothing to any tree's length, but still counts towards the
+#' block's gain cost, which is a property of the hierarchy rather than of the
+#' taxa sampled.
+#'
 #' ## Cost matrix
 #'
 #' - **Absent → present (gain):** cost = \eqn{n + 1}, where \eqn{n} is the
@@ -61,9 +73,10 @@
 #'         row \code{i} giving the 1-based level index of each secondary for
 #'         present-state \code{i + 1}.}
 #'       \item{`tip_sec_known`}{Integer matrix (\code{n_tip × n_secondary}).
-#'         For tips with \code{tip_states == -2}, column \code{s} holds the
-#'         1-based level index of secondary \code{s} if it was observed, or
-#'         0 if it was unknown; used to constrain the admissible states of a
+#'         For tips with \code{tip_states == -2}, column \code{s} holds a
+#'         bit mask of the levels secondary \code{s} may take at that tip
+#'         (bit \code{i - 1} set = level \code{i} admissible), or 0 where it
+#'         is unconstrained; used to restrict the admissible states of a
 #'         partially-known combination.}
 #'     }
 #'   }
@@ -82,6 +95,8 @@ RecodeHierarchy <- function(dataset, hierarchy) {
 
   idx <- attr(dataset, "index")
   allLevels <- attr(dataset, "allLevels")
+  levels <- attr(dataset, "levels")
+  contrast <- attr(dataset, "contrast")
   nChar <- length(idx)
   nTip <- length(dataset)
 
@@ -89,6 +104,21 @@ RecodeHierarchy <- function(dataset, hierarchy) {
   origMat <- do.call(rbind, lapply(dataset, function(x) {
     allLevels[x[idx]]
   }))
+  # ... and as `contrast` row indices, which is what a secondary's state space
+  # must be read from.  Deriving it from the token strings instead made an
+  # ambiguity token such as "{01}" a level of its own -- one Hamming step from
+  # both "0" and "1" rather than matching either, and one more factor in the
+  # combination count (T-393).  `tokenLevels` is the R-side counterpart of
+  # `DataSet::token_states`, which the HSJ path already reads.
+  tokenMat <- do.call(rbind, lapply(dataset, function(x) x[idx]))
+  applicable <- which(levels != "-")
+  tokenLevels <- lapply(seq_len(nrow(contrast)), function(tk) {
+    applicable[contrast[tk, applicable] > 0]
+  })
+  # A token admitting every applicable state establishes no state at all; this
+  # is the role "?" played under the old string test, and an ambiguity spanning
+  # the whole state space says exactly as much.
+  tokenGeneric <- lengths(tokenLevels) == length(applicable)
 
   .RecodeBlock <- function(node) {
     ctrl <- node$controlling
@@ -99,13 +129,26 @@ RecodeHierarchy <- function(dataset, hierarchy) {
            "Block controlled by character ", ctrl, " has sub-hierarchies.")
     }
 
-    # Informative levels for each secondary (exclude "-" and "?")
+    # Informative levels for each secondary, as state indices into `levels`
     secLevels <- lapply(deps, function(d) {
-      sort(setdiff(unique(origMat[, d]), c("-", "?")))
+      tokens <- unique(tokenMat[, d])
+      sort(unique(unlist(tokenLevels[tokens[!tokenGeneric[tokens]]])))
     })
     secNStates <- vapply(secLevels, length, integer(1))
+    if (any(secNStates > 31L)) {
+      stop("Secondary character ", deps[which.max(secNStates)],
+           " has more than 31 informative states; the x-transformation ",
+           "cannot recode it.")
+    }
+    # A secondary with no informative level -- every tip gap or fully ambiguous
+    # -- is carried as one unobserved level rather than dropped, so that a
+    # present primary still has a state to take (a zero-width state space made
+    # every tip cost infinite, T-394) and the block's gain cost still reflects
+    # how many secondaries the primary controls, not how many the sampled taxa
+    # happen to resolve.
+    secNLevels <- pmax(secNStates, 1L)
 
-    nPresent <- prod(secNStates)
+    nPresent <- prod(secNLevels)
     nStates <- nPresent + 1L
     nSec <- length(deps)
 
@@ -120,7 +163,7 @@ RecodeHierarchy <- function(dataset, hierarchy) {
     # All present-state combinations (expand.grid: first dim varies fastest)
     if (nSec > 0L) {
       comboGrid <- as.matrix(expand.grid(
-        lapply(secLevels, seq_along)
+        lapply(secNLevels, seq_len)
       ))
     } else {
       # No secondaries: 2 states (absent + one present)
@@ -145,13 +188,14 @@ RecodeHierarchy <- function(dataset, hierarchy) {
     }
 
     # --- Tip states ---
-    # `tipSecKnown[t, s]` records, per tip and per secondary, the 1-based
-    # level index of that secondary IF it was observed for this tip, or 0 if
-    # it was unknown ("-"/"?"/unrecognised token). Only consulted when
-    # `tipStates[t] == -2` (present, but not every secondary was resolvable):
-    # it lets the admissible-state set be restricted to combinations
-    # consistent with whichever secondaries WERE observed, rather than
-    # freeing every present state (T-379).
+    # `tipSecKnown[t, s]` records, per tip and per secondary, a bit mask of the
+    # levels that secondary may take at this tip (bit i - 1 = level i), or 0
+    # where it is unconstrained. Only consulted when `tipStates[t] == -2`
+    # (present, but not every secondary was resolvable): it lets the
+    # admissible-state set be restricted to combinations consistent with
+    # whatever the secondaries WERE observed to be, rather than freeing every
+    # present state (T-379). A mask rather than a single level index because a
+    # polymorphic token narrows a secondary without resolving it (T-393).
     tipStates <- integer(nTip)
     tipSecKnown <- matrix(0L, nrow = nTip, ncol = nSec)
     for (t in seq_len(nTip)) {
@@ -171,28 +215,37 @@ RecodeHierarchy <- function(dataset, hierarchy) {
         next
       }
 
-      secVals <- origMat[t, deps]
+      secVals <- tokenMat[t, deps]
       anyUnknown <- FALSE
       levelIndices <- integer(nSec)
+      secMasks <- integer(nSec)
       known <- logical(nSec)
 
       for (s in seq_len(nSec)) {
-        if (secVals[s] %in% c("-", "?")) {
+        if (tokenGeneric[[secVals[s]]]) {
           anyUnknown <- TRUE
           next
         }
-        mi <- match(secVals[s], secLevels[[s]])
-        if (is.na(mi)) {
-          anyUnknown <- TRUE
+        # Positions, within this secondary's levels, that its token admits
+        pos <- match(tokenLevels[[secVals[s]]], secLevels[[s]])
+        pos <- pos[!is.na(pos)]
+        if (length(pos) == 1L) {
+          levelIndices[s] <- pos
+          known[s] <- TRUE
           next
         }
-        levelIndices[s] <- mi
-        known[s] <- TRUE
+        anyUnknown <- TRUE
+        # Admitting every level (or none of them) constrains nothing, and 0
+        # says so more cheaply than the equivalent full mask.
+        if (length(pos) > 0L && length(pos) < secNStates[[s]]) {
+          secMasks[s] <- sum(bitwShiftL(1L, pos - 1L))
+        }
       }
 
       if (anyUnknown) {
-        tipStates[t] <- -2L  # present, one or more secondaries unknown
-        tipSecKnown[t, known] <- levelIndices[known]
+        tipStates[t] <- -2L  # present, one or more secondaries unresolved
+        secMasks[known] <- bitwShiftL(1L, levelIndices[known] - 1L)
+        tipSecKnown[t, ] <- secMasks
         next
       }
 
@@ -201,7 +254,7 @@ RecodeHierarchy <- function(dataset, hierarchy) {
       multiplier <- 1L
       for (s in seq_len(nSec)) {
         rowIdx <- rowIdx + (levelIndices[s] - 1L) * multiplier
-        multiplier <- multiplier * secNStates[s]
+        multiplier <- multiplier * secNLevels[s]
       }
       tipStates[t] <- rowIdx  # 1-based present state = Sankoff state index
     }
diff --git a/man/RecodeHierarchy.Rd b/man/RecodeHierarchy.Rd
index 91e24135a..05097a161 100644
--- a/man/RecodeHierarchy.Rd
+++ b/man/RecodeHierarchy.Rd
@@ -30,9 +30,10 @@ row-major: \code{cost_matrix[from, to]}.}
 row \code{i} giving the 1-based level index of each secondary for
 present-state \code{i + 1}.}
 \item{\code{tip_sec_known}}{Integer matrix (\code{n_tip × n_secondary}).
-For tips with \code{tip_states == -2}, column \code{s} holds the
-1-based level index of secondary \code{s} if it was observed, or
-0 if it was unknown; used to constrain the admissible states of a
+For tips with \code{tip_states == -2}, column \code{s} holds a
+bit mask of the levels secondary \code{s} may take at that tip
+(bit \code{i - 1} set = level \code{i} admissible), or 0 where it
+is unconstrained; used to restrict the admissible states of a
 partially-known combination.}
 }
 }
@@ -54,6 +55,18 @@ State 0 represents "primary absent".
 States \eqn{1 \ldots \prod k_i} represent all possible combinations of
 secondary character states (where \eqn{k_i} is the number of informative
 states of secondary character \eqn{i}).
+
+The informative levels of a secondary character are read from the dataset's
+\code{contrast} matrix, not from the token strings it carries.  An ambiguity token
+such as \code{"{01}"} therefore denotes the \emph{set} of states it contrasts against,
+as it does under \code{inapplicable = "hsj"}, rather than becoming a level of its
+own; a token that admits every applicable state (\code{"?"}, or an ambiguity
+spanning them all) shows that no state is established and so contributes
+none.  A secondary whose levels are all unobserved -- reachable by dropping
+taxa from a dataset that validated -- is carried as a single unobserved
+level: it adds nothing to any tree's length, but still counts towards the
+block's gain cost, which is a property of the hierarchy rather than of the
+taxa sampled.
 }
 
 \subsection{Cost matrix}{
diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp
index fcb6678c0..3d6d2e3e6 100644
--- a/src/ts_rcpp.cpp
+++ b/src/ts_rcpp.cpp
@@ -1921,8 +1921,9 @@ static void unpack_xform(Nullable<List> xformConfig,
       }
       int ns = ns_vec[ch];
       // Only needed to resolve state == -2 (present, secondaries partially
-      // unknown); combo_grid is n_present x n_sec, tip_sec_known is
-      // n_tip x n_sec (see RecodeHierarchy()).
+      // unknown); combo_grid is n_present x n_sec holding 1-based level
+      // indices, tip_sec_known is n_tip x n_sec holding a per-secondary bit
+      // mask of admissible levels, 0 = unconstrained (see RecodeHierarchy()).
       IntegerMatrix combo_grid = as<IntegerMatrix>(rc["combo_grid"]);
       IntegerMatrix tip_sec = as<IntegerMatrix>(rc["tip_sec_known"]);
       int n_sec = combo_grid.ncol();
@@ -1933,15 +1934,16 @@ static void unpack_xform(Nullable<List> xformConfig,
         if (state == -1) {
           for (int s = 0; s < ns; ++s) tip_ptr[s] = 0.0;
         } else if (state == -2) {
-          // Present, but one or more secondaries were unknown for this tip.
-          // Restrict admissible present-states to those consistent with the
-          // secondaries that WERE observed (T-379); previously this freed
+          // Present, but one or more secondaries were unresolved for this tip.
+          // Restrict admissible present-states to those consistent with what
+          // the secondaries WERE observed to be (T-379); previously this freed
           // every present state regardless of any known secondaries.
           for (int s = 1; s < ns; ++s) {
             bool admissible = true;
             for (int d = 0; d < n_sec; ++d) {
               int known = tip_sec(t, d);
-              if (known != 0 && combo_grid(s - 1, d) != known) {
+              if (known != 0 &&
+                  (known & (1 << (combo_grid(s - 1, d) - 1))) == 0) {
                 admissible = false;
                 break;
               }
@@ -3294,10 +3296,11 @@ List ts_sankoff_test(
                tip_states_r.nrow(), n_tip, n_tip);
   }
 
-  // combo_grids_r[ch] (n_present x n_sec) and tip_sec_known_r[ch]
-  // (n_tip x n_sec) resolve state == -2 to the states consistent with
-  // whichever secondaries WERE observed (T-379); absent (NULL), -2 falls
-  // back to freeing every present state, as before.
+  // combo_grids_r[ch] (n_present x n_sec, 1-based level indices) and
+  // tip_sec_known_r[ch] (n_tip x n_sec, per-secondary bit mask of admissible
+  // levels, 0 = unconstrained) resolve state == -2 to the states consistent
+  // with what the secondaries WERE observed to be (T-379); absent (NULL), -2
+  // falls back to freeing every present state, as before.
   bool have_combo = combo_grids_r.isNotNull() && tip_sec_known_r.isNotNull();
   List combo_grids, tip_sec_knowns;
   if (have_combo) {
@@ -3328,7 +3331,8 @@ List ts_sankoff_test(
           bool admissible = true;
           for (int d = 0; d < n_sec; ++d) {
             int known = tip_sec(t, d);
-            if (known != 0 && combo_grid(s - 1, d) != known) {
+            if (known != 0 &&
+                (known & (1 << (combo_grid(s - 1, d) - 1))) == 0) {
               admissible = false;
               break;
             }
diff --git a/tests/testthat/test-ts-xform-statespace.R b/tests/testthat/test-ts-xform-statespace.R
new file mode 100644
index 000000000..e352e4d81
--- /dev/null
+++ b/tests/testthat/test-ts-xform-statespace.R
@@ -0,0 +1,222 @@
+# Tier 2: skipped on CRAN; see tests/testing-strategy.md
+skip_on_cran()
+
+# Regressions for red-team T-393, T-394 and T-401: how a hierarchy block's
+# secondary state space is derived, what happens when a secondary has no
+# observed state, and what `MaximizeParsimony()` may score its own output with.
+
+library("TreeTools")
+
+# `MatrixToPhyDat()` rather than the `phangorn::phyDat()` helper the rest of the
+# xform suite uses: only it puts an ambiguity token such as "{01}" in the
+# contrast matrix as the state SET it denotes, which is the whole subject here.
+xss_dat <- function(mat) MatrixToPhyDat(mat)
+
+xss_tree <- function() ape::read.tree(text = "(((t1,t2),(t3,t4)),(t5,t6));")
+
+
+# ===== T-393: an ambiguity token is a state set, not a state =================
+# Deriving a secondary's levels from the observed token STRINGS admitted "{01}"
+# as a level of its own, one Hamming step from both "0" and "1" rather than
+# matching either -- so a polymorphic cell cost a step no resolution of it
+# needs, and every polymorphic cell multiplied the combination count.
+
+test_that("Polymorphic secondary costs no more than its best resolution", {
+  mat <- matrix(c(
+    "1", "0",
+    "1", "1",
+    "1", "{01}",
+    "1", "0",
+    "0", "-",
+    "1", "1"
+  ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
+  ds <- xss_dat(mat)
+  h <- CharacterHierarchy("1" = 2L)
+  tree <- xss_tree()
+
+  # "{01}" is a subset of {"0", "1"}, so the block's length is the smallest any
+  # concrete resolution attains -- never more.
+  resolved <- vapply(c("0", "1"), function(state) {
+    m <- mat
+    m[3, 2] <- state
+    TreeLength(tree, xss_dat(m), hierarchy = h, inapplicable = "xform")
+  }, numeric(1))
+  expect_equal(
+    TreeLength(tree, ds, hierarchy = h, inapplicable = "xform"),
+    min(resolved)
+  )
+
+  # The binary secondary has two levels, not three.
+  expect_equal(RecodeHierarchy(ds, h)$sankoff_chars[[1]]$n_states, 3)
+})
+
+
+test_that("Polymorphic cells do not inflate the state space", {
+  # Four binary secondaries = 2^4 + 1 = 17 states.  Reading a state space off
+  # the token strings made each polymorphic cell a third level, giving 82 and a
+  # spurious "> 32 states" warning (with the quadratic-per-node cost to match).
+  mat <- matrix(c(
+    "1", "0",    "0",    "0",    "0",
+    "1", "1",    "1",    "1",    "1",
+    "1", "{01}", "0",    "1",    "0",
+    "1", "0",    "{01}", "0",    "1",
+    "1", "1",    "0",    "{01}", "1",
+    "1", "0",    "1",    "0",    "{01}",
+    "0", "-",    "-",    "-",    "-"
+  ), nrow = 7, byrow = TRUE, dimnames = list(paste0("t", 1:7), NULL))
+  ds <- xss_dat(mat)
+  h <- CharacterHierarchy("1" = 2:5)
+
+  expect_silent(recoded <- RecodeHierarchy(ds, h))
+  expect_equal(recoded$sankoff_chars[[1]]$n_states, 17)
+})
+
+
+test_that("Polymorphism narrows a multistate secondary without freeing it", {
+  # With three levels available, "{01}" is neither resolved nor unconstrained:
+  # t2 must not be allowed to take state "2" to match its sister t1.
+  mat <- matrix(c(
+    "1", "2",
+    "1", "{01}",
+    "1", "2",
+    "1", "0",
+    "0", "-",
+    "1", "1"
+  ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
+  ds <- xss_dat(mat)
+  h <- CharacterHierarchy("1" = 2L)
+  tree <- xss_tree()
+
+  resolved <- vapply(c("0", "1", "2"), function(state) {
+    m <- mat
+    m[2, 2] <- state
+    TreeLength(tree, xss_dat(m), hierarchy = h, inapplicable = "xform")
+  }, numeric(1))
+  # Precondition: resolving to "2" is strictly cheaper here, so treating the
+  # token as wholly unknown would be a measurable under-count.
+  expect_lt(resolved[["2"]], min(resolved[c("0", "1")]))
+
+  expect_equal(
+    TreeLength(tree, ds, hierarchy = h, inapplicable = "xform"),
+    min(resolved[c("0", "1")])
+  )
+})
+
+
+# ===== T-394: a secondary with no observed state ============================
+# `ValidateHierarchy()` runs before both taxon-subsetting sites, so dropping a
+# taxon can leave a secondary all-gap/all-missing in a dataset that validated.
+# A zero-width state space then admitted no state at any tip: every tip cost
+# was infinite and the block's length `Inf`.
+
+xss_degenerate <- function() {
+  # Character 3 is resolved at s5 alone; scoring any tree over {s1..s4} drops
+  # s5, leaving char 3 with nothing observed.
+  matrix(c(
+    "1", "0", "?",
+    "1", "1", "?",
+    "1", "0", "?",
+    "0", "-", "-",
+    "1", "1", "0"
+  ), nrow = 5, byrow = TRUE, dimnames = list(paste0("s", 1:5), NULL))
+}
+
+test_that("Unobserved secondary leaves a finite length after taxon dropping", {
+  ds <- xss_dat(xss_degenerate())
+  h <- CharacterHierarchy("1" = 2:3)
+  tree <- ape::read.tree(text = "((s1,s2),(s3,s4));")
+
+  # Precondition: the dataset as supplied validates and recodes normally.
+  expect_equal(RecodeHierarchy(ds, h)$sankoff_chars[[1]]$n_states, 3)
+
+  expect_true(is.finite(
+    TreeLength(tree, ds, hierarchy = h, inapplicable = "xform")))
+})
+
+
+test_that("Search over a subset with an unobserved secondary completes", {
+  ds <- xss_dat(xss_degenerate())
+  h <- CharacterHierarchy("1" = 2:3)
+  tree <- ape::read.tree(text = "((s1,s2),(s3,s4));")
+
+  res <- suppressWarnings(
+    MaximizeParsimony(ds, tree = tree, hierarchy = h, inapplicable = "xform",
+                      maxReplicates = 2L, verbosity = 0L))
+  expect_true(is.finite(attr(res, "score")))
+})
+
+
+test_that("A pool with no finite length is reported, not branched on", {
+  # `diff(range(c(Inf, Inf)))` is `NaN`, which aborted the reporting `if` with
+  # "missing value where TRUE/FALSE needed" whatever produced the infinities.
+  expect_warning(reported <- TreeSearch:::.ReportXformScore(c(Inf, Inf), 7),
+                 "no finite x-transformation length")
+  expect_equal(reported, 7)
+
+  # A pool that is only partly infinite still reports the shortest finite
+  # length (and, these two differing, warns about that too).
+  expect_warning(
+    expect_warning(mixed <- TreeSearch:::.ReportXformScore(c(Inf, 3, 5), 7),
+                   "no finite x-transformation length"),
+    "do not share a length")
+  expect_equal(mixed, 3)
+})
+
+
+# ===== T-401: the report block must not score a contracted tree =============
+# `MaximizeParsimony()` rescores its own output at the canonical rooting.  With
+# `collapse = TRUE` that output may be polytomous, and `TreeLength()` scores the
+# topology it is given: the length of a polytomy is that of its best resolution,
+# which the Sankoff kernel does not compute.  Only the HSJ/XFORM no-op in
+# `compute_collapsed_flags_aggressive()` keeps the returned trees binary today.
+
+test_that("Contracted trees are not scored as though binary", {
+  # The T-330 fixture, whose equal-weights arm genuinely collapses (10 -> 8
+  # edges) while the hierarchy character supports the contracted clade.
+  mat <- matrix(c(
+    "0", "0",
+    "0", "0",
+    "1", "0",
+    "1", "1",
+    "1", "1",
+    "1", "-"
+  ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
+  ds <- phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"),
+                         ambiguity = "?")
+  h <- CharacterHierarchy("2" = integer(0))
+  binary <- Preorder(RenumberTips(
+    ape::read.tree(text = "(((t1,t2),(t3,(t4,t5))),t6);"), names(ds)))
+
+  at <- attributes(ds)
+  collapsed <- TreeSearch:::ts_collapse_pool(
+    list(binary[["edge"]]), at$contrast,
+    matrix(unlist(ds, use.names = FALSE), nrow = length(ds), byrow = TRUE),
+    as.integer(TreeSearch:::.NonHierarchyWeights(ds, h)), at$levels,
+    list(min_steps = integer(0), concavity = Inf, xpiwe = FALSE,
+         xpiwe_r = 0.5, xpiwe_max_f = 5.0, obs_count = integer(0),
+         infoAmounts = NULL),
+    NULL, NULL, NULL)
+  polytomous <- Renumber(structure(
+    list(edge = collapsed$trees[[1]],
+         Nnode = max(collapsed$trees[[1]]) - length(ds),
+         tip.label = names(ds)),
+    class = "phylo"))
+  expect_lt(dim(polytomous[["edge"]])[[1]], dim(binary[["edge"]])[[1]])
+
+  reference <- TreeLength(binary, ds, hierarchy = h, inapplicable = "xform")
+
+  # Precondition: handing the contracted tree straight to `TreeLength()` does
+  # not reproduce that length.  Tolerant of an error, because the kernel-level
+  # boundary check on non-binary edge matrices (T-400) turns the silent form of
+  # this into a loud one.
+  direct <- tryCatch(
+    TreeLength(structure(list(polytomous), class = "multiPhylo"), ds,
+               hierarchy = h, inapplicable = "xform"),
+    error = function(e) NA_real_)
+  expect_false(isTRUE(all.equal(direct, reference)))
+
+  # The reporting path scores the binary pool the tree was contracted from.
+  expect_equal(
+    TreeSearch:::.XformPoolScore(list(polytomous), list(binary), ds, h, -1),
+    reference)
+})

From d7049de5b821277175c748028d1d1430e50f863a Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:01:12 +0100
Subject: [PATCH 2/6] docs: record the x-transformation state-space changes

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .AGENTS/memory/feature-inapplicable.md | 21 +++++++++++++++++++
 NEWS.md                                | 28 ++++++++++++++++++++++++++
 2 files changed, 49 insertions(+)

diff --git a/.AGENTS/memory/feature-inapplicable.md b/.AGENTS/memory/feature-inapplicable.md
index 7c776b65c..9e2e3f26f 100644
--- a/.AGENTS/memory/feature-inapplicable.md
+++ b/.AGENTS/memory/feature-inapplicable.md
@@ -94,3 +94,24 @@ secondaries supported (state count = ∏k_i + 1). Nested hierarchies deferred.
 Integration complete: `ScoringMode::XFORM` in `score_tree()` dispatches
 Fitch(non-hierarchy) + Sankoff(recoded). `MaximizeParsimony()` accepts
 `inapplicable = "xform"`. End-to-end search verified.
+
+### State spaces come from `contrast`, never from token strings (T-393/T-394)
+
+A secondary's levels are the applicable states its column's tokens contrast
+against — the R-side counterpart of `DataSet::token_states`, which HSJ reads.
+Do not reach for `unique()` on the token strings: an ambiguity token such as
+`"{01}"` is then a level of its own (a Hamming step from both `"0"` and `"1"`,
+and one more factor in `prod(k_i)`), and level ordering follows locale
+collation. A token admitting every applicable state establishes none.
+
+Two consequences to preserve when touching `RecodeHierarchy()`:
+
+- `tip_sec_known` holds a **bit mask** of admissible levels per (tip,
+  secondary), 0 = unconstrained — not a single level index. Both consumers are
+  in `src/ts_rcpp.cpp` (`unpack_xform` and `ts_sankoff_test`); change them
+  together.
+- A secondary with no observed level is carried as one unobserved level, not
+  dropped. `ValidateHierarchy()` runs before both taxon-subsetting sites
+  (`R/MaximizeParsimony.R`, `R/tree_length.R`), so dropping a taxon can leave a
+  block degenerate in a dataset that validated. Keeping it in `nSec` keeps the
+  gain cost a property of the hierarchy, not of the taxa sampled.
diff --git a/NEWS.md b/NEWS.md
index 3492cf2d4..8c92f96dc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,33 @@
 # To integrate into 2.0.0 notes
 
+- `inapplicable = "xform"` no longer treats an ambiguity token as a state of its
+  own.  A hierarchy block took its secondary characters' states from the distinct
+  tokens observed in each column, so a polymorphic cell such as `{01}` was
+  admitted alongside `0` and `1` and sat one step from both rather than matching
+  either: a tree paid for a step that no resolution of the polymorphism requires.
+  Each such cell also multiplied the block's combination count -- four binary
+  secondary characters with one polymorphic cell apiece produced 82 states where
+  17 suffice, tripping the "large state space" warning and paying its quadratic
+  per-node cost -- and left the ordering of states at the mercy of the locale's
+  string collation.  States are now read from the dataset's contrast matrix, as
+  `inapplicable = "hsj"` already did.
+
+  **X-transformation lengths of data containing polymorphic or partly ambiguous
+  secondary characters will therefore change, and will not increase.**  A tree's
+  length is now the smallest that any resolution of its ambiguity attains, which
+  is what the criterion means by it.  Data coded only with unambiguous tokens,
+  `?` and `-` are unaffected.
+
+- `inapplicable = "xform"` no longer returns an infinite length when a secondary
+  character has no observed state.  Scoring a tree drops any taxon the tree does
+  not bear, and dropping the only taxon at which a secondary character was
+  resolved left that character with an empty state space: no state was
+  admissible at any tip, so `TreeLength()` returned `Inf` without comment and
+  `MaximizeParsimony()` stopped with "missing value where TRUE/FALSE needed".
+  Such a character is now carried as a single unobserved state, contributing
+  nothing to any tree's length while still counting towards its block's gain
+  cost -- which the hierarchy fixes, not the taxa that happen to be sampled.
+
 - `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

From ca2baaa5572fbfd05421587f1d78fd8d27cbc967 Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:25:14 +0100
Subject: [PATCH 3/6] fix: warn when a contracted pool is scored; tighten docs
 and tests

Addresses review of the x-transformation state-space fix: reporting the
length of the binary pool a contracted tree came from is itself the T-385
discrepancy, so say so; the collapse flags are decided over the Fitch term
alone, so the two pools' lengths are equal only if a future collapse fix is
hierarchy-aware, which the comment no longer assumes.

Also: state the 31-level ceiling and the max(k, 1) state count in the docs and
vignettes, pin the state space and the multi-bit mask that the length
assertions alone did not, and drive both mask consumers end to end.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 .AGENTS/memory/feature-inapplicable.md    |  10 +-
 NEWS.md                                   |   4 +-
 R/MaximizeParsimony.R                     |  27 ++++--
 R/recode_hierarchy.R                      |  17 ++--
 man/RecodeHierarchy.Rd                    |  11 ++-
 tests/testthat/test-ts-xform-statespace.R | 110 +++++++++++++++++-----
 vignettes/inapplicable.Rmd                |  10 ++
 vignettes/search-algorithm.Rmd            |   6 ++
 8 files changed, 147 insertions(+), 48 deletions(-)

diff --git a/.AGENTS/memory/feature-inapplicable.md b/.AGENTS/memory/feature-inapplicable.md
index 9e2e3f26f..f322ee1a1 100644
--- a/.AGENTS/memory/feature-inapplicable.md
+++ b/.AGENTS/memory/feature-inapplicable.md
@@ -111,7 +111,9 @@ Two consequences to preserve when touching `RecodeHierarchy()`:
   in `src/ts_rcpp.cpp` (`unpack_xform` and `ts_sankoff_test`); change them
   together.
 - A secondary with no observed level is carried as one unobserved level, not
-  dropped. `ValidateHierarchy()` runs before both taxon-subsetting sites
-  (`R/MaximizeParsimony.R`, `R/tree_length.R`), so dropping a taxon can leave a
-  block degenerate in a dataset that validated. Keeping it in `nSec` keeps the
-  gain cost a property of the hierarchy, not of the taxa sampled.
+  dropped. `ValidateHierarchy()` asks that a secondary be coded inapplicable
+  where its primary is absent, never that any state of it remain observed, so
+  the taxon subsetting in `MaximizeParsimony()` / `TreeLength()` can leave a
+  block degenerate in a dataset that validated — and re-validating after the
+  subset would not catch it either. Keeping it in `nSec` keeps the gain cost a
+  property of the hierarchy, not of the taxa sampled.
diff --git a/NEWS.md b/NEWS.md
index 8c92f96dc..c2c17421a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,7 +16,9 @@
   secondary characters will therefore change, and will not increase.**  A tree's
   length is now the smallest that any resolution of its ambiguity attains, which
   is what the criterion means by it.  Data coded only with unambiguous tokens,
-  `?` and `-` are unaffected.
+  `?` and `-` are unaffected by this change (though they may be affected by the
+  next).  A single secondary character may now take at most 31 states, and one
+  exceeding that is reported as an error rather than recoded incorrectly.
 
 - `inapplicable = "xform"` no longer returns an infinite length when a secondary
   character has no observed state.  Scoring a tree drops any taxon the tree does
diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R
index 1fe339d70..ba9f82d7b 100644
--- a/R/MaximizeParsimony.R
+++ b/R/MaximizeParsimony.R
@@ -1776,10 +1776,16 @@ MaximizeParsimony <- function(
 # reports a number that is too small (T-401).  Nothing but the HSJ/XFORM no-op
 # in `compute_collapsed_flags_aggressive()` currently keeps `outTrees` binary,
 # and this must not depend on that staying in place: score instead the binary
-# pool the trees were contracted from, whose lengths are the same because only
-# zero-length branches are removed.  Kept separate from its caller because a
+# pool the trees were contracted from.  Kept separate from its caller because a
 # default search cannot reach that path, so this is the only place it can be
 # exercised.
+#
+# That substitution reports the length of a tree the user is not handed, which
+# is the discrepancy T-385 was filed for, so it warns.  Whoever lifts the
+# HSJ/XFORM no-op can retire the warning only by making the collapse itself
+# hierarchy-aware: the flags are decided over `ds.blocks[]`, the Fitch term
+# alone, so an XFORM-blind unblocking would contract branches that are not
+# zero-length under the Sankoff term and the two pools would genuinely differ.
 .XformPoolScore <- function(pool, binaryPool, dataset, hierarchy, fallback) {
   nEdgeBinary <- 2L * length(dataset) - 2L
   .Binary <- function(trees) {
@@ -1787,6 +1793,10 @@ MaximizeParsimony <- function(
   }
   if (!all(.Binary(pool))) {
     pool <- binaryPool[.Binary(binaryPool)]
+    warning("Returned trees contain polytomies, whose x-transformation length ",
+            "is that of their best resolution; reporting the length of the ",
+            "binary trees they were contracted from, which `TreeLength()` of a ",
+            "returned tree need not reproduce.", call. = FALSE)
   }
   if (length(pool) == 0L) {
     return(fallback)
@@ -1800,7 +1810,9 @@ MaximizeParsimony <- function(
 # Reduce a pool's canonical-rooting lengths to the one number reported.
 # `fallback` covers a pool with no finite length: `diff(range(.))` is then
 # `NaN`, which `if` cannot branch on -- the abort a degenerate hierarchy block
-# used to produce (T-394).
+# used to produce (T-394).  A non-finite member is reported here and then
+# dropped, so it does not also raise the T-374 residue warning below; silence
+# there means "the finite lengths agree", not "the residue is fixed".
 .ReportXformScore <- function(canonicalScores, fallback) {
   finite <- is.finite(canonicalScores)
   if (!all(finite)) {
@@ -1808,7 +1820,7 @@ MaximizeParsimony <- function(
             " returned trees have no ",
             "finite x-transformation length; ",
             if (any(finite)) "reporting the shortest of the rest."
-            else "reporting the search's own score.")
+            else "reporting the search's own score.", call. = FALSE)
   }
   if (!any(finite)) {
     return(fallback)
@@ -1818,13 +1830,12 @@ MaximizeParsimony <- function(
   if (diff(range(canonicalScores)) > sqrt(.Machine$double.eps)) {
     # Pool membership is chosen on search-time scores taken at differing
     # rootings (`result$scores` in the caller), 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.
+    # parsimonious can differ once scored at one rooting.  Not silently averaged
+    # away: this is the open residue of T-374.
     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.")
+            "rooting-dependent -- see ?MaximizeParsimony.", call. = FALSE)
   }
   min(canonicalScores)
 }
diff --git a/R/recode_hierarchy.R b/R/recode_hierarchy.R
index f6c43c470..538fd9eb4 100644
--- a/R/recode_hierarchy.R
+++ b/R/recode_hierarchy.R
@@ -4,15 +4,16 @@
 #' \insertCite{Goloboff2021;textual}{TreeSearch}.
 #' Each hierarchy block (one controlling primary character plus \eqn{n}
 #' secondary characters) is combined into a single step-matrix character
-#' with \eqn{\prod k_i + 1} states and an asymmetric cost matrix.
+#' with \eqn{\prod \max(k_i, 1) + 1} states and an asymmetric cost matrix.
 #'
 #' @details
 #' ## State encoding
 #'
 #' State 0 represents "primary absent".
-#' States \eqn{1 \ldots \prod k_i} represent all possible combinations of
-#' secondary character states (where \eqn{k_i} is the number of informative
-#' states of secondary character \eqn{i}).
+#' States \eqn{1 \ldots \prod \max(k_i, 1)} represent all possible combinations
+#' of secondary character states (where \eqn{k_i} is the number of informative
+#' states of secondary character \eqn{i}; a secondary with none contributes a
+#' single unobserved state, as below).
 #'
 #' The informative levels of a secondary character are read from the dataset's
 #' `contrast` matrix, not from the token strings it carries.  An ambiguity token
@@ -24,7 +25,7 @@
 #' taxa from a dataset that validated -- is carried as a single unobserved
 #' level: it adds nothing to any tree's length, but still counts towards the
 #' block's gain cost, which is a property of the hierarchy rather than of the
-#' taxa sampled.
+#' taxa sampled.  A single secondary may take at most 31 levels.
 #'
 #' ## Cost matrix
 #'
@@ -226,9 +227,11 @@ RecodeHierarchy <- function(dataset, hierarchy) {
           anyUnknown <- TRUE
           next
         }
-        # Positions, within this secondary's levels, that its token admits
+        # Positions, within this secondary's levels, that its token admits.
+        # `secLevels[[s]]` is the union over the non-generic tokens of this very
+        # column, so a non-generic token's states are all levels of it and the
+        # match cannot fail.
         pos <- match(tokenLevels[[secVals[s]]], secLevels[[s]])
-        pos <- pos[!is.na(pos)]
         if (length(pos) == 1L) {
           levelIndices[s] <- pos
           known[s] <- TRUE
diff --git a/man/RecodeHierarchy.Rd b/man/RecodeHierarchy.Rd
index 05097a161..81d5d62df 100644
--- a/man/RecodeHierarchy.Rd
+++ b/man/RecodeHierarchy.Rd
@@ -46,15 +46,16 @@ Implements the x-transformation recoding of
 \insertCite{Goloboff2021;textual}{TreeSearch}.
 Each hierarchy block (one controlling primary character plus \eqn{n}
 secondary characters) is combined into a single step-matrix character
-with \eqn{\prod k_i + 1} states and an asymmetric cost matrix.
+with \eqn{\prod \max(k_i, 1) + 1} states and an asymmetric cost matrix.
 }
 \details{
 \subsection{State encoding}{
 
 State 0 represents "primary absent".
-States \eqn{1 \ldots \prod k_i} represent all possible combinations of
-secondary character states (where \eqn{k_i} is the number of informative
-states of secondary character \eqn{i}).
+States \eqn{1 \ldots \prod \max(k_i, 1)} represent all possible combinations
+of secondary character states (where \eqn{k_i} is the number of informative
+states of secondary character \eqn{i}; a secondary with none contributes a
+single unobserved state, as below).
 
 The informative levels of a secondary character are read from the dataset's
 \code{contrast} matrix, not from the token strings it carries.  An ambiguity token
@@ -66,7 +67,7 @@ none.  A secondary whose levels are all unobserved -- reachable by dropping
 taxa from a dataset that validated -- is carried as a single unobserved
 level: it adds nothing to any tree's length, but still counts towards the
 block's gain cost, which is a property of the hierarchy rather than of the
-taxa sampled.
+taxa sampled.  A single secondary may take at most 31 levels.
 }
 
 \subsection{Cost matrix}{
diff --git a/tests/testthat/test-ts-xform-statespace.R b/tests/testthat/test-ts-xform-statespace.R
index e352e4d81..07bdff242 100644
--- a/tests/testthat/test-ts-xform-statespace.R
+++ b/tests/testthat/test-ts-xform-statespace.R
@@ -10,9 +10,9 @@ library("TreeTools")
 # `MatrixToPhyDat()` rather than the `phangorn::phyDat()` helper the rest of the
 # xform suite uses: only it puts an ambiguity token such as "{01}" in the
 # contrast matrix as the state SET it denotes, which is the whole subject here.
-xss_dat <- function(mat) MatrixToPhyDat(mat)
+XssDat <- function(mat) MatrixToPhyDat(mat)
 
-xss_tree <- function() ape::read.tree(text = "(((t1,t2),(t3,t4)),(t5,t6));")
+XssTree <- function() ape::read.tree(text = "(((t1,t2),(t3,t4)),(t5,t6));")
 
 
 # ===== T-393: an ambiguity token is a state set, not a state =================
@@ -30,16 +30,16 @@ test_that("Polymorphic secondary costs no more than its best resolution", {
     "0", "-",
     "1", "1"
   ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
-  ds <- xss_dat(mat)
+  ds <- XssDat(mat)
   h <- CharacterHierarchy("1" = 2L)
-  tree <- xss_tree()
+  tree <- XssTree()
 
   # "{01}" is a subset of {"0", "1"}, so the block's length is the smallest any
   # concrete resolution attains -- never more.
   resolved <- vapply(c("0", "1"), function(state) {
     m <- mat
     m[3, 2] <- state
-    TreeLength(tree, xss_dat(m), hierarchy = h, inapplicable = "xform")
+    TreeLength(tree, XssDat(m), hierarchy = h, inapplicable = "xform")
   }, numeric(1))
   expect_equal(
     TreeLength(tree, ds, hierarchy = h, inapplicable = "xform"),
@@ -64,7 +64,7 @@ test_that("Polymorphic cells do not inflate the state space", {
     "1", "0",    "1",    "0",    "{01}",
     "0", "-",    "-",    "-",    "-"
   ), nrow = 7, byrow = TRUE, dimnames = list(paste0("t", 1:7), NULL))
-  ds <- xss_dat(mat)
+  ds <- XssDat(mat)
   h <- CharacterHierarchy("1" = 2:5)
 
   expect_silent(recoded <- RecodeHierarchy(ds, h))
@@ -72,6 +72,17 @@ test_that("Polymorphic cells do not inflate the state space", {
 })
 
 
+test_that("A secondary beyond the mask's width is reported, not mis-recoded", {
+  # One bit per level, so 31 is the most a secondary can carry.
+  tokens <- c(0:9, LETTERS)[1:32]
+  mat <- cbind(c("0", rep("1", 32)), c("-", tokens))
+  rownames(mat) <- paste0("t", seq_len(33))
+
+  expect_error(RecodeHierarchy(XssDat(mat), CharacterHierarchy("1" = 2L)),
+               "more than 31 informative states")
+})
+
+
 test_that("Polymorphism narrows a multistate secondary without freeing it", {
   # With three levels available, "{01}" is neither resolved nor unconstrained:
   # t2 must not be allowed to take state "2" to match its sister t1.
@@ -83,14 +94,14 @@ test_that("Polymorphism narrows a multistate secondary without freeing it", {
     "0", "-",
     "1", "1"
   ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
-  ds <- xss_dat(mat)
+  ds <- XssDat(mat)
   h <- CharacterHierarchy("1" = 2L)
-  tree <- xss_tree()
+  tree <- XssTree()
 
   resolved <- vapply(c("0", "1", "2"), function(state) {
     m <- mat
     m[2, 2] <- state
-    TreeLength(tree, xss_dat(m), hierarchy = h, inapplicable = "xform")
+    TreeLength(tree, XssDat(m), hierarchy = h, inapplicable = "xform")
   }, numeric(1))
   # Precondition: resolving to "2" is strictly cheaper here, so treating the
   # token as wholly unknown would be a measurable under-count.
@@ -100,16 +111,51 @@ test_that("Polymorphism narrows a multistate secondary without freeing it", {
     TreeLength(tree, ds, hierarchy = h, inapplicable = "xform"),
     min(resolved[c("0", "1")])
   )
+
+  recoded <- RecodeHierarchy(ds, h)$sankoff_chars[[1]]
+  # Three levels plus absent, not four plus absent: the length assertion above
+  # holds under the old encoding too, so this is what pins the state space.
+  expect_equal(recoded$n_states, 4)
+  # A mask of more than one bit, which is what distinguishes the mask encoding
+  # from the single level index it replaced.
+  expect_equal(recoded$tip_sec_known[2, 1], 3L)
+})
+
+
+test_that("Search and TreeLength agree on a multi-bit secondary mask", {
+  # The mask is read in two places -- `unpack_xform()` for the search and
+  # `ts_sankoff_test()` for `TreeLength()` -- and a divergence between them
+  # would mis-guide the search while the reported score stayed self-consistent.
+  # Only a mask of more than one bit tells the two encodings apart, so drive
+  # both over data that produces one.
+  mat <- matrix(c(
+    "1", "2",
+    "1", "{01}",
+    "1", "2",
+    "1", "0",
+    "0", "-",
+    "1", "1"
+  ), nrow = 6, byrow = TRUE, dimnames = list(paste0("t", 1:6), NULL))
+  ds <- XssDat(mat)
+  h <- CharacterHierarchy("1" = 2L)
+  expect_equal(RecodeHierarchy(ds, h)$sankoff_chars[[1]]$tip_sec_known[2, 1], 3L)
+
+  res <- MaximizeParsimony(ds, tree = XssTree(), hierarchy = h,
+                           inapplicable = "xform", maxReplicates = 3L,
+                           verbosity = 0L)
+  expect_equal(attr(res, "score"),
+               min(TreeLength(res, ds, hierarchy = h, inapplicable = "xform")))
 })
 
 
 # ===== T-394: a secondary with no observed state ============================
-# `ValidateHierarchy()` runs before both taxon-subsetting sites, so dropping a
-# taxon can leave a secondary all-gap/all-missing in a dataset that validated.
-# A zero-width state space then admitted no state at any tip: every tip cost
-# was infinite and the block's length `Inf`.
+# Validation asks that a secondary be coded inapplicable where its primary is
+# absent, never that any state of it remain observed, so dropping a taxon can
+# leave a secondary all-gap/all-missing in a dataset that validated.  A
+# zero-width state space then admitted no state at any tip: every tip cost was
+# infinite and the block's length `Inf`.
 
-xss_degenerate <- function() {
+XssDegenerate <- function() {
   # Character 3 is resolved at s5 alone; scoring any tree over {s1..s4} drops
   # s5, leaving char 3 with nothing observed.
   matrix(c(
@@ -122,27 +168,33 @@ xss_degenerate <- function() {
 }
 
 test_that("Unobserved secondary leaves a finite length after taxon dropping", {
-  ds <- xss_dat(xss_degenerate())
+  ds <- XssDat(XssDegenerate())
   h <- CharacterHierarchy("1" = 2:3)
   tree <- ape::read.tree(text = "((s1,s2),(s3,s4));")
 
   # Precondition: the dataset as supplied validates and recodes normally.
   expect_equal(RecodeHierarchy(ds, h)$sankoff_chars[[1]]$n_states, 3)
 
-  expect_true(is.finite(
-    TreeLength(tree, ds, hierarchy = h, inapplicable = "xform")))
+  # The finite value, not merely finiteness: HSJ scores this 1.5, and every
+  # concrete reading of the unobserved secondary gives the same 2.
+  expect_equal(TreeLength(tree, ds, hierarchy = h, inapplicable = "xform"), 2)
 })
 
 
 test_that("Search over a subset with an unobserved secondary completes", {
-  ds <- xss_dat(xss_degenerate())
+  ds <- XssDat(XssDegenerate())
   h <- CharacterHierarchy("1" = 2:3)
   tree <- ape::read.tree(text = "((s1,s2),(s3,s4));")
 
   res <- suppressWarnings(
     MaximizeParsimony(ds, tree = tree, hierarchy = h, inapplicable = "xform",
                       maxReplicates = 2L, verbosity = 0L))
-  expect_true(is.finite(attr(res, "score")))
+  # Not merely finite: the reported score must still be the length of a tree
+  # returned, which is what the taxon-subset path threatens.
+  expect_equal(
+    attr(res, "score"),
+    min(TreeLength(res, TreeSearch:::.Recompress(ds[res[[1]][["tip.label"]]]),
+                   hierarchy = h, inapplicable = "xform")))
 })
 
 
@@ -215,8 +267,20 @@ test_that("Contracted trees are not scored as though binary", {
     error = function(e) NA_real_)
   expect_false(isTRUE(all.equal(direct, reference)))
 
-  # The reporting path scores the binary pool the tree was contracted from.
-  expect_equal(
-    TreeSearch:::.XformPoolScore(list(polytomous), list(binary), ds, h, -1),
-    reference)
+  # The reporting path scores the binary pool the tree was contracted from,
+  # and says so -- that length is not one `TreeLength()` of a returned tree
+  # reproduces, which is the discrepancy T-385 was filed for.
+  expect_warning(
+    substituted <- TreeSearch:::.XformPoolScore(list(polytomous), list(binary),
+                                               ds, h, -1),
+    "Returned trees contain polytomies")
+  expect_equal(substituted, reference)
+
+  # With nothing binary to fall back on, the search's own score is reported
+  # rather than a length read off a contracted tree.
+  expect_warning(
+    fellBack <- TreeSearch:::.XformPoolScore(list(polytomous),
+                                             list(polytomous), ds, h, -1),
+    "Returned trees contain polytomies")
+  expect_equal(fellBack, -1)
 })
diff --git a/vignettes/inapplicable.Rmd b/vignettes/inapplicable.Rmd
index 403c95230..d516f99ba 100644
--- a/vignettes/inapplicable.Rmd
+++ b/vignettes/inapplicable.Rmd
@@ -79,6 +79,16 @@ This asymmetry captures the idea that independently evolving a complex
 structure (and its associated secondary characters) is less parsimonious
 than losing it.
 
+A secondary character's states are read from the dataset's contrast matrix, so
+an ambiguity token such as `{01}` denotes the set of states it stands for
+rather than a state of its own: a tree's length is the smallest that any
+resolution of the ambiguity attains.
+A secondary character with no observed state -- which scoring a tree that bears
+only some of the taxa can produce -- contributes nothing to any tree's length,
+but still counts towards its block's gain cost, since that cost is fixed by the
+hierarchy rather than by the taxa sampled.
+A single secondary character may take at most 31 states.
+
 ```{r xform, eval = FALSE}
 hierarchy <- CharacterHierarchy("1" = 2:5)
 MaximizeParsimony(dataset, hierarchy = hierarchy,
diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd
index a87d2b292..cc3358559 100644
--- a/vignettes/search-algorithm.Rmd
+++ b/vignettes/search-algorithm.Rmd
@@ -308,11 +308,17 @@ output stage has since canonicalised.
 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.
+Where `collapse = TRUE` has contracted a branch, the binary trees it was
+contracted from are scored instead, with a warning: the length of a contracted
+tree is the length of its best resolution rather than of the edge matrix as
+given, which the Sankoff kernel does not compute.
 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.
+It likewise warns, and falls back to the score the search itself compared,
+should a returned tree have no finite length at all.
 
 ### Zero-length edge skipping
 

From df3bd73d0ca8cce23c9b224f602b2d98fce916dd Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Tue, 4 Aug 2026 17:38:09 +0100
Subject: [PATCH 4/6] test: stop driving a polytomy through the scoring kernel

The T-401 regression test scored a contracted tree directly, to show that
doing so does not reproduce the binary length.  That call is the very defect
the test exists to guard against: until the kernel bounds-checks a non-binary
edge matrix (T-400) it writes past the end of its Fitch word vector, which
tryCatch cannot contain.  It read as a wrong length on x86 and aborted
ubuntu-arm64 with 'double free or corruption'.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 tests/testthat/test-ts-xform-statespace.R | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/tests/testthat/test-ts-xform-statespace.R b/tests/testthat/test-ts-xform-statespace.R
index 07bdff242..b0973cb14 100644
--- a/tests/testthat/test-ts-xform-statespace.R
+++ b/tests/testthat/test-ts-xform-statespace.R
@@ -257,15 +257,12 @@ test_that("Contracted trees are not scored as though binary", {
 
   reference <- TreeLength(binary, ds, hierarchy = h, inapplicable = "xform")
 
-  # Precondition: handing the contracted tree straight to `TreeLength()` does
-  # not reproduce that length.  Tolerant of an error, because the kernel-level
-  # boundary check on non-binary edge matrices (T-400) turns the silent form of
-  # this into a loud one.
-  direct <- tryCatch(
-    TreeLength(structure(list(polytomous), class = "multiPhylo"), ds,
-               hierarchy = h, inapplicable = "xform"),
-    error = function(e) NA_real_)
-  expect_false(isTRUE(all.equal(direct, reference)))
+  # What handing `polytomous` straight to `TreeLength()` returns is deliberately
+  # NOT asserted here: until the kernel bounds-checks a non-binary edge matrix
+  # (T-400) that call writes past the end of its Fitch word vector, which no
+  # `tryCatch()` can contain -- measured as a silently wrong length on x86 and a
+  # `double free or corruption` abort on arm64.  That is the behaviour this call
+  # site must avoid, so the test must not perform it either.
 
   # The reporting path scores the binary pool the tree was contracted from,
   # and says so -- that length is not one `TreeLength()` of a returned tree

From fe79c502609e34647fe93d0ada727240c25d58e7 Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Tue, 4 Aug 2026 18:06:38 +0100
Subject: [PATCH 5/6] fix: keep the polytomy warning truthful; drop the
 redundant token-string matrix

The warning fired before the substituted pool was known to be non-empty, so
it could promise a length recomputed from binary trees and then report the
search's own score instead.

`RecodeHierarchy()` built a full taxon-by-character character matrix that,
once the state space came from the contrast matrix, was read only to fetch
each tip's primary token; the integer token matrix already carries it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 R/MaximizeParsimony.R | 10 +++++++---
 R/recode_hierarchy.R  | 18 +++++++-----------
 2 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R
index ba9f82d7b..405a7b86a 100644
--- a/R/MaximizeParsimony.R
+++ b/R/MaximizeParsimony.R
@@ -1794,9 +1794,13 @@ MaximizeParsimony <- function(
   if (!all(.Binary(pool))) {
     pool <- binaryPool[.Binary(binaryPool)]
     warning("Returned trees contain polytomies, whose x-transformation length ",
-            "is that of their best resolution; reporting the length of the ",
-            "binary trees they were contracted from, which `TreeLength()` of a ",
-            "returned tree need not reproduce.", call. = FALSE)
+            "is that of their best resolution; reporting ",
+            if (length(pool) > 0L) {
+              paste0("the length of the binary trees they were contracted ",
+                     "from, which `TreeLength()` of a returned tree need not ",
+                     "reproduce.")
+            } else "the search's own score.",
+            call. = FALSE)
   }
   if (length(pool) == 0L) {
     return(fallback)
diff --git a/R/recode_hierarchy.R b/R/recode_hierarchy.R
index 538fd9eb4..75288800e 100644
--- a/R/recode_hierarchy.R
+++ b/R/recode_hierarchy.R
@@ -101,16 +101,12 @@ RecodeHierarchy <- function(dataset, hierarchy) {
   nChar <- length(idx)
   nTip <- length(dataset)
 
-  # Original character matrix (taxon × char), as token strings
-  origMat <- do.call(rbind, lapply(dataset, function(x) {
-    allLevels[x[idx]]
-  }))
-  # ... and as `contrast` row indices, which is what a secondary's state space
-  # must be read from.  Deriving it from the token strings instead made an
-  # ambiguity token such as "{01}" a level of its own -- one Hamming step from
-  # both "0" and "1" rather than matching either, and one more factor in the
-  # combination count (T-393).  `tokenLevels` is the R-side counterpart of
-  # `DataSet::token_states`, which the HSJ path already reads.
+  # Original character matrix (taxon × char), as `contrast` row indices -- which
+  # is what a secondary's state space must be read from.  Reading it off the
+  # token strings instead made an ambiguity token such as "{01}" a level of its
+  # own: one Hamming step from both "0" and "1" rather than matching either, and
+  # one more factor in the combination count (T-393).  `tokenLevels` is the
+  # R-side counterpart of `DataSet::token_states`, which the HSJ path reads.
   tokenMat <- do.call(rbind, lapply(dataset, function(x) x[idx]))
   applicable <- which(levels != "-")
   tokenLevels <- lapply(seq_len(nrow(contrast)), function(tk) {
@@ -200,7 +196,7 @@ RecodeHierarchy <- function(dataset, hierarchy) {
     tipStates <- integer(nTip)
     tipSecKnown <- matrix(0L, nrow = nTip, ncol = nSec)
     for (t in seq_len(nTip)) {
-      pri <- origMat[t, ctrl]
+      pri <- allLevels[tokenMat[t, ctrl]]
 
       if (pri == "?") {
         tipStates[t] <- -1L  # fully ambiguous

From 21de433f84d2eb9766d5631b9d6a7f25b7643e3d Mon Sep 17 00:00:00 2001
From: R script <1695515+ms609@users.noreply.github.com>
Date: Wed, 5 Aug 2026 04:07:29 +0100
Subject: [PATCH 6/6] test: do not claim the two mask readings are shown to
 agree

MaximizeParsimony() derives its reported score by calling TreeLength() on the
pool (T-385), so asserting the two agree is true by construction.  Compare the
search's optimum against a brute force over all six-taxon topologies instead,
and state plainly that what keeps the two mask consumers honest is being
edited together, not this test.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 tests/testthat/test-ts-xform-statespace.R | 34 +++++++++++++++++------
 1 file changed, 25 insertions(+), 9 deletions(-)

diff --git a/tests/testthat/test-ts-xform-statespace.R b/tests/testthat/test-ts-xform-statespace.R
index b0973cb14..92b11d324 100644
--- a/tests/testthat/test-ts-xform-statespace.R
+++ b/tests/testthat/test-ts-xform-statespace.R
@@ -122,12 +122,20 @@ test_that("Polymorphism narrows a multistate secondary without freeing it", {
 })
 
 
-test_that("Search and TreeLength agree on a multi-bit secondary mask", {
+test_that("A multi-bit mask drives the search without mis-scoring", {
   # The mask is read in two places -- `unpack_xform()` for the search and
-  # `ts_sankoff_test()` for `TreeLength()` -- and a divergence between them
-  # would mis-guide the search while the reported score stayed self-consistent.
-  # Only a mask of more than one bit tells the two encodings apart, so drive
-  # both over data that produces one.
+  # `ts_sankoff_test()` for `TreeLength()`.  This drives the first of them over
+  # data that produces a mask of more than one bit, which is the only shape that
+  # tells the mask encoding apart from the single level index it replaced.
+  #
+  # It does NOT establish that the two readings agree, and no test here does.
+  # `MaximizeParsimony()` derives its reported score by calling `TreeLength()`
+  # on the pool (T-385), so comparing the two is true by construction; and on
+  # this data every reading of the mask yields the same optimum (3, brute-forced
+  # over all 105 six-taxon topologies for `{01}` and for each of `0`, `1`, `2`
+  # and `?`), so the search cannot be misled into a measurable difference
+  # either.  What keeps the two sites honest is that they are edited together --
+  # recorded in `.AGENTS/memory/feature-inapplicable.md`.
   mat <- matrix(c(
     "1", "2",
     "1", "{01}",
@@ -143,8 +151,14 @@ test_that("Search and TreeLength agree on a multi-bit secondary mask", {
   res <- MaximizeParsimony(ds, tree = XssTree(), hierarchy = h,
                            inapplicable = "xform", maxReplicates = 3L,
                            verbosity = 0L)
-  expect_equal(attr(res, "score"),
-               min(TreeLength(res, ds, hierarchy = h, inapplicable = "xform")))
+  # The optimum, computed without reference to anything the search reports.
+  topologies <- lapply(seq_len(NUnrooted(6)), function(i) {
+    RootTree(as.phylo(i - 1L, 6, tipLabels = names(ds)), 1)
+  })
+  expect_equal(
+    attr(res, "score"),
+    min(TreeLength(structure(topologies, class = "multiPhylo"), ds,
+                   hierarchy = h, inapplicable = "xform")))
 })
 
 
@@ -189,8 +203,10 @@ test_that("Search over a subset with an unobserved secondary completes", {
   res <- suppressWarnings(
     MaximizeParsimony(ds, tree = tree, hierarchy = h, inapplicable = "xform",
                       maxReplicates = 2L, verbosity = 0L))
-  # Not merely finite: the reported score must still be the length of a tree
-  # returned, which is what the taxon-subset path threatens.
+  # Not merely finite: this also pins `.XformPoolScore()`'s edge-count test
+  # against the SUBSET dataset.  Were it measured against the dataset as
+  # supplied, no returned tree would look binary, the pool would empty and the
+  # search's own mid-search score would be reported here instead.
   expect_equal(
     attr(res, "score"),
     min(TreeLength(res, TreeSearch:::.Recompress(ds[res[[1]][["tip.label"]]]),