diff --git a/.AGENTS/memory/feature-inapplicable.md b/.AGENTS/memory/feature-inapplicable.md index 7c776b65c..f322ee1a1 100644 --- a/.AGENTS/memory/feature-inapplicable.md +++ b/.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. diff --git a/NEWS.md b/NEWS.md index a08d94daf..83aaf4aa2 100644 --- a/NEWS.md +++ b/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 diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 3ccce1f44..405a7b86a 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,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 diff --git a/R/recode_hierarchy.R b/R/recode_hierarchy.R index 2cc435a50..75288800e 100644 --- a/R/recode_hierarchy.R +++ b/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 } diff --git a/man/RecodeHierarchy.Rd b/man/RecodeHierarchy.Rd index 91e24135a..81d5d62df 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.} } } @@ -45,15 +46,28 @@ 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 +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. A single secondary may take at most 31 levels. } \subsection{Cost matrix}{ diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index fe9abee38..dd7415ed5 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -1979,8 +1979,9 @@ static void unpack_xform(Nullable 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(rc["combo_grid"]); IntegerMatrix tip_sec = as(rc["tip_sec_known"]); int n_sec = combo_grid.ncol(); @@ -2007,15 +2008,16 @@ static void unpack_xform(Nullable 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; } @@ -3377,10 +3379,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) { @@ -3411,7 +3414,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..92b11d324 --- /dev/null +++ b/tests/testthat/test-ts-xform-statespace.R @@ -0,0 +1,299 @@ +# 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. +XssDat <- function(mat) MatrixToPhyDat(mat) + +XssTree <- 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 <- XssDat(mat) + h <- CharacterHierarchy("1" = 2L) + 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, XssDat(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 <- XssDat(mat) + h <- CharacterHierarchy("1" = 2:5) + + expect_silent(recoded <- RecodeHierarchy(ds, h)) + expect_equal(recoded$sankoff_chars[[1]]$n_states, 17) +}) + + +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. + 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) + tree <- XssTree() + + resolved <- vapply(c("0", "1", "2"), function(state) { + m <- mat + m[2, 2] <- state + 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. + expect_lt(resolved[["2"]], min(resolved[c("0", "1")])) + + expect_equal( + 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("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()`. 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}", + "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) + # 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"))) +}) + + +# ===== T-394: a secondary with no observed state ============================ +# 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`. + +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( + "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 <- 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) + + # 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 <- 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)) + # 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"]]]), + hierarchy = h, inapplicable = "xform"))) +}) + + +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") + + # 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 + # 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