fix: derive XFORM secondary state spaces from contrasts, not token strings

.AGENTS/memory/feature-inapplicable.md

@@ -94,3 +94,26 @@ 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()` 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.

NEWS.md

@@ -1,5 +1,34 @@
 # 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 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
+  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.
 - `TreeLength()`, `CharacterLength()`, `TreeScore()` and `EdgeListScore()` -- and
   so `Consistency()`, `ExpectedLength()`, `ConcordantInformation()`,
   `LengthAdded()` and `SuccessiveApproximations()`, which score trees through

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,84 @@ 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.  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) {
+    vapply(trees, function(tr) dim(tr[["edge"]])[[1]], integer(1)) == nEdgeBinary
+  }
+  if (!all(.Binary(pool))) {
+    pool <- binaryPool[.Binary(binaryPool)]
+    warning("Returned trees contain polytomies, whose x-transformation length ",
+            "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)
+  }
+  .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).  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)) {
+    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.", call. = FALSE)
+  }
+  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.
+    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.", call. = FALSE)
+  }
+  min(canonicalScores)
+}
+
 #' Launch tree search graphical user interface
 #'
 #' Opens a "shiny" app for interactive parsimony tree search and results

R/recode_hierarchy.R

@@ -4,15 +4,28 @@
 #' \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
+#' 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.  A single secondary may take at most 31 levels.
 #'
 #' ## Cost matrix
 #'
@@ -61,9 +74,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,13 +96,26 @@ 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)
 
-  # Original character matrix (taxon × char), as token strings
-  origMat <- do.call(rbind, lapply(dataset, function(x) {
-    allLevels[x[idx]]
-  }))
+  # 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) {
+    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 +126,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 +160,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,17 +185,18 @@ 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)) {
-      pri <- origMat[t, ctrl]
+      pri <- allLevels[tokenMat[t, ctrl]]
 
       if (pri == "?") {
         tipStates[t] <- -1L  # fully ambiguous
@@ -171,28 +212,39 @@ 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.
+        # `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]])
+        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 +253,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
     }