diff --git a/NEWS.md b/NEWS.md index 3492cf2d4..47724a8cb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,51 @@ # To integrate into 2.0.0 notes +- `inapplicable = "hsj"` no longer charges for an inapplicable secondary + character. Where a secondary was coded `"-"` but its controlling primary did + not certainly code the structure absent -- because the primary was `"?"`, or + because the matrix codes the structure present and the secondary inapplicable + anyway -- the `"-"` was treated as an ordinary state of that secondary. Being + disjoint from every other state, it could then be propagated inwards and used + to label an internal node in the middle of a region where the structure *is* + present, mismatching every secondary of the block at once and charging that + branch the full weight of the scaling parameter. + + This was the behaviour the method exists to avoid. Hopkins & St John (2021) + count only the secondary characters that apply, and note that treating + inapplicable cells as a separate state "increases the dissimilarity of all + pairwise comparisons", overweighting the controlling primary and favouring + clades that separate taxa possessing the structure from those lacking it. An + inapplicable secondary now contributes nothing to the dissimilarity, whatever + its controlling primary codes. + + **HSJ scores may therefore fall on matrices that code a secondary + inapplicable where its primary does not code absence**, by the scaling + parameter divided by the number of secondary characters in the block, for each + affected branch. Scores are unchanged wherever the matrix is coded + consistently, and remain independent of where the tree is rooted. + +- Nested hierarchies now validate, so a `CharacterHierarchy()` describing + tertiary characters can be scored under `inapplicable = "hsj"`. A + sub-controlling character is deliberately recorded both as a dependent of the + character above it and as the controlling character of the one below -- that + dual role is what nesting means -- but validation counted the second + occurrence as one character appearing in two blocks, and so rejected every + nested hierarchy that could be written, including the one documented in + `?CharacterHierarchy`. Genuine double claims are still rejected. + + `HierarchyFromNames()` now detects nesting as its documentation describes. A + controlling character whose tag extends another's, as `sup_tail_tip` extends + `sup_tail`, is nested beneath it, and each dependent attaches to the longest + tag it extends, so `sub_tail_tip_gloss` belongs to `sup_tail_tip` rather than + to `sup_tail`. Previously a tag was matched only as far as its first + underscore, which collapsed every depth onto the outermost tag; nesting was + silently dropped and could not be expressed at all. A shared prefix without + an underscore boundary does not nest, so `sup_tailfin` remains independent of + `sup_tail`. + + The x-transformation still does not implement nesting, and now says so + directly instead of failing validation first. + - `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 diff --git a/R/CharacterHierarchy.R b/R/CharacterHierarchy.R index 9a6f07e9a..47fd16ee1 100644 --- a/R/CharacterHierarchy.R +++ b/R/CharacterHierarchy.R @@ -183,7 +183,15 @@ ValidateHierarchy <- function(hierarchy, dataset) { claimed <- integer(0) - .ValidateBlock <- function(node, depth = 1L) { + # `ctrlClaimedByParent` is TRUE for a nested block, whose controlling + # character is *deliberately* also a dependent of its parent: .ParseOneBlock() + # records a sub-controller in both places, because it is simultaneously a + # secondary of the character above it and the primary of the one below. That + # dual role is the whole content of "nested", so counting it as a second claim + # rejected every nested hierarchy that could be written -- including this + # file's own documented example. The double-claim check still applies in full + # to the block's own dependents, and to a controlling character at top level. + .ValidateBlock <- function(node, depth = 1L, ctrlClaimedByParent = FALSE) { ctrl <- node$controlling deps <- node$dependents @@ -198,14 +206,15 @@ ValidateHierarchy <- function(hierarchy, dataset) { } # Check no double-claiming - overlap <- intersect(allIdx, claimed) + newIdx <- if (ctrlClaimedByParent) deps else allIdx + overlap <- intersect(newIdx, claimed) if (length(overlap) > 0L) { stop(sprintf( "Character(s) %s appear in multiple hierarchy blocks.", paste(overlap, collapse = ", ") )) } - claimed <<- c(claimed, allIdx) + claimed <<- c(claimed, newIdx) # Check controlling character is binary (has exactly states "0" and "1", # possibly with inapplicable/missing) @@ -236,9 +245,10 @@ ValidateHierarchy <- function(hierarchy, dataset) { } } - # Recurse into children + # Recurse into children. A child's controlling character was already + # claimed just above, as one of this block's dependents. for (child in node$children) { - .ValidateBlock(child, depth + 1L) + .ValidateBlock(child, depth + 1L, ctrlClaimedByParent = TRUE) } } @@ -254,9 +264,16 @@ ValidateHierarchy <- function(hierarchy, dataset) { #' #' Parse character names following the TNT convention where controlling #' characters are named `sup_` and their dependent characters are -#' named `sub_[_suffix]`. Tags must match between a controlling -#' character and its dependents. Nested hierarchies are detected when a -#' `sub_` character is also a `sup_` for further characters. +#' named `sub_[_suffix]`. Each dependent is attached to the longest +#' `sup_` tag that its own tag extends, at an underscore boundary. +#' +#' Nested hierarchies are written by giving a controlling character a tag that +#' itself extends another controlling character's tag: `sup_tail_tip` controls +#' the tag `tail_tip`, and because `tail_tip` extends `tail` it is also a +#' dependent of `sup_tail`. Longest-match is what keeps `sub_tail_tip_gloss` +#' with `sup_tail_tip` rather than with `sup_tail`. A shared prefix alone does +#' not nest: `sup_tailfin` is independent of `sup_tail`, since `tailfin` does +#' not extend `tail` at an underscore boundary. #' #' @param charNames Character vector of names, one per original character. #' @@ -268,6 +285,11 @@ ValidateHierarchy <- function(hierarchy, dataset) { #' "sup_wing", "sub_wing_venation", "eyes") #' HierarchyFromNames(names) #' +#' # Nesting: `tail_tip` extends `tail`, so character 3 is both a dependent of +#' # character 1 and the controlling character of character 4. +#' HierarchyFromNames(c("sup_tail", "sub_tail_colour", +#' "sup_tail_tip", "sub_tail_tip_gloss")) +#' #' @family tree scoring #' @seealso [CharacterHierarchy()] #' @export @@ -284,72 +306,79 @@ HierarchyFromNames <- function(charNames) { return(NULL) } - # Extract tags + # Extract tags. A tag may itself contain underscores, and that is what makes + # nesting expressible: `sup_tail_tip` controls the tag "tail_tip", and because + # "tail_tip" extends the existing tag "tail" it is simultaneously a dependent + # of `sup_tail`. This is the only way one character can play both roles, since + # a single name cannot carry both the `sup_` and `sub_` prefix -- which is why + # the previous `intersect(subIdx, supIdx)` test could never fire, leaving the + # documented nesting support unreachable. + # + # Every sub_/sup_ character attaches to the LONGEST sup_ tag its own tag + # extends. Longest-wins is what keeps `sub_tail_tip_gloss` with `sup_tail_tip` + # instead of with `sup_tail`; matching only the first underscore-delimited + # component (as this function used to) collapsed every depth onto the outermost + # tag, so no dependent could ever reach a nested controller. supTags <- sub("^sup_", "", charNames[supIdx]) subTagsFull <- sub("^sub_", "", charNames[subIdx]) - # The tag is the first component before any additional underscore-suffix - # e.g. "sub_tail_colour" → tag = "tail" - subTags <- sub("_.*", "", subTagsFull) - # Build mapping: tag → controlling index, tag → dependent indices + # Build mapping: tag → controlling index tagToSup <- setNames(supIdx, supTags) - # Group sub characters by tag - tagToSubs <- split(subIdx, subTags) + # The longest sup_ tag that `tag` sits under: an exact match, or an extension + # at an underscore boundary (so "tailfin" does NOT sit under "tail"). + # `exclude` stops a sup_ tag being its own parent. + .ParentTag <- function(tag, exclude = "") { + eligible <- supTags[supTags != exclude] + if (!length(eligible)) return(NA_character_) + under <- tag == eligible | startsWith(tag, paste0(eligible, "_")) + if (!any(under)) return(NA_character_) + cand <- eligible[under] + cand[[which.max(nchar(cand))]] + } - # Check for sub_ characters referencing nonexistent sup_ tags - orphanTags <- setdiff(names(tagToSubs), supTags) - if (length(orphanTags) > 0L) { + subParent <- vapply(subTagsFull, .ParentTag, character(1), USE.NAMES = FALSE) + orphan <- is.na(subParent) + if (any(orphan)) { warning(sprintf( "sub_ characters reference tags with no corresponding sup_: %s", - paste(orphanTags, collapse = ", ") + paste(unique(subTagsFull[orphan]), collapse = ", ") )) } - # Detect nested hierarchies: a sub_ character that is also a sup_ - # Find sub_ chars that are also in supIdx - subAlsoSup <- intersect(subIdx, supIdx) - - # Build hierarchy - # First pass: create flat blocks for all sup_ tags - args <- list() - for (tag in supTags) { - ctrl <- tagToSup[[tag]] - subs <- tagToSubs[[tag]] - if (is.null(subs)) subs <- integer(0) - - # Check which subs are themselves controlling (nested hierarchy) - nestedSubs <- intersect(subs, supIdx) - flatSubs <- setdiff(subs, supIdx) - - if (length(nestedSubs) == 0L) { - # Simple block - args[[as.character(ctrl)]] <- as.integer(subs) - } else { - # Nested: build list with named sub-hierarchies - block <- as.list(as.integer(flatSubs)) - for (ns in nestedSubs) { - nsTag <- supTags[supIdx == ns] - nsSubs <- tagToSubs[[nsTag]] - if (is.null(nsSubs)) nsSubs <- integer(0) - block[[as.character(ns)]] <- as.integer(nsSubs) - } - args[[as.character(ctrl)]] <- block + # A sup_ tag that extends another sup_ tag is a nested controller. + supParent <- vapply(supTags, function(s) .ParentTag(s, exclude = s), + character(1), USE.NAMES = FALSE) + + # Dependents of one tag: its own sub_ characters, plus a named sub-hierarchy + # for each sup_ tag nested directly beneath it. Recursive, so nesting works + # to arbitrary depth. A nested controller is NOT added to the flat dependents + # here -- .ParseOneBlock() records a named sub-controller as a dependent of the + # enclosing block itself. + .BuildTag <- function(tag) { + flat <- as.integer(subIdx[!orphan & subParent == tag]) + kidTags <- supTags[!is.na(supParent) & supParent == tag] + if (!length(kidTags)) { + return(flat) + } + block <- as.list(flat) + for (k in kidTags) { + block[[as.character(tagToSup[[k]])]] <- .BuildTag(k) } + block } - # Filter out sup_ chars whose index also appears in subIdx - # (they'll be included as children of their parent) - topLevelSup <- setdiff(supIdx, subIdx) - if (length(topLevelSup) == 0L) { - # All sup_ characters are also sub_ — circular or all nested. - # Fall back to treating all as top-level with a warning. - warning("All sup_ characters are also sub_ characters. ", + topTags <- supTags[is.na(supParent)] + if (!length(topTags)) { + # Every sup_ tag extends another, which needs a cycle and so cannot arise + # from prefix matching; kept as a guard rather than a reachable branch. + warning("Every sup_ tag extends another sup_ tag. ", "Treating all as top-level.") - topLevelSup <- supIdx + topTags <- supTags } - topLevelCtrls <- as.character(topLevelSup) - args <- args[topLevelCtrls] + + args <- lapply(topTags, .BuildTag) + names(args) <- as.character(tagToSup[topTags]) do.call(CharacterHierarchy, args) } diff --git a/man/HierarchyFromNames.Rd b/man/HierarchyFromNames.Rd index 436e0f622..06456751a 100644 --- a/man/HierarchyFromNames.Rd +++ b/man/HierarchyFromNames.Rd @@ -16,15 +16,28 @@ detected. \description{ Parse character names following the TNT convention where controlling characters are named \verb{sup_} and their dependent characters are -named \verb{sub_[_suffix]}. Tags must match between a controlling -character and its dependents. Nested hierarchies are detected when a -\code{sub_} character is also a \code{sup_} for further characters. +named \verb{sub_[_suffix]}. Each dependent is attached to the longest +\code{sup_} tag that its own tag extends, at an underscore boundary. +} +\details{ +Nested hierarchies are written by giving a controlling character a tag that +itself extends another controlling character's tag: \code{sup_tail_tip} controls +the tag \code{tail_tip}, and because \code{tail_tip} extends \code{tail} it is also a +dependent of \code{sup_tail}. Longest-match is what keeps \code{sub_tail_tip_gloss} +with \code{sup_tail_tip} rather than with \code{sup_tail}. A shared prefix alone does +not nest: \code{sup_tailfin} is independent of \code{sup_tail}, since \code{tailfin} does +not extend \code{tail} at an underscore boundary. } \examples{ names <- c("sup_tail", "sub_tail_colour", "sub_tail_shape", "sup_wing", "sub_wing_venation", "eyes") HierarchyFromNames(names) +# Nesting: `tail_tip` extends `tail`, so character 3 is both a dependent of +# character 1 and the controlling character of character 4. +HierarchyFromNames(c("sup_tail", "sub_tail_colour", + "sup_tail_tip", "sub_tail_tip_gloss")) + } \seealso{ \code{\link[=CharacterHierarchy]{CharacterHierarchy()}} diff --git a/src/ts_hsj.cpp b/src/ts_hsj.cpp index 9c68a193f..436c2252e 100644 --- a/src/ts_hsj.cpp +++ b/src/ts_hsj.cpp @@ -128,21 +128,44 @@ static CanonOrder build_canon_order(const TreeState& tree) { // which is what state_sets must hold to make the Fitch downpass/uppass below // correct for ambiguous tokens. // -// `pri_free[t]` marks the tips at which this secondary carries no constraint: -// those whose controlling primary CANNOT code the structure present, so the -// character does not exist there and its "-" is not a state it takes (T-374). -// score_hierarchy_block() computes it and documents why the test is that -// strict one rather than "may be absent". Admitting "-" as an ordinary -// concrete state -- as this function formerly did, and as the comment here -// formerly asserted was deliberate -- let the uppass propagate it INWARDS and -// resolve a node in the middle of the PRESENT region to it, where it is -// disjoint from every present neighbour in every secondary at once, and -// score_hierarchy_block() charged that branch d = m, the full alpha, for a -// node that by construction has no inapplicable secondaries. That over-charge -// is wrong under any rooting (the paper's d counts "nonmatching secondary -// characters", p.5, among characters that APPLY), and because whether it fired -// depended on the DELTRAN direction it was also the dominant source of -// T-374's rooting-dependence. +// THE INAPPLICABLE TOKEN IS NOT A STATE OF A SECONDARY CHARACTER. This is the +// paper's central claim, not a convention we are free to pick: Hopkins & St John +// (2021) define d as "the number of nonmatching secondary characters dependent +// on that primary" (p.5) among the characters that APPLY, and state plainly that +// where secondaries are inapplicable to a taxon "they have no influence on the +// estimated dissimilarity" (p.5). The paper introduces HSJ precisely to avoid +// the alternative: "Treating inapplicable characters as a new, separate state +// will ... skew the analysis, because having a new separate state increases the +// dissimilarity of all pairwise comparisons ... This results in overweighting +// the [controlling] primary character and favors clades that separate taxa with +// secondary characters from those without" (p.5). Admitting "-" as a concrete +// state is therefore not a stricter reading of HSJ; it is the FitchS behaviour +// HSJ exists to replace. +// +// So `inapp_bit` is stripped from every secondary's observed set, at every tip, +// regardless of what the controlling primary codes. Two distinct routes reach a +// non-constraining tip and both must be handled here: +// * `pri_free[t]` -- the controlling primary CANNOT code the structure present, +// so the character does not exist at t (T-374). score_hierarchy_block() +// computes this and documents why the test is that strict one rather than +// the laxer "may be absent". +// * the secondary's own token carries no applicable state once "-" is removed +// -- i.e. the cell is coded "-" (or an ambiguity resolving only to "-") +// while the primary does NOT certainly code absence. That combination is +// contradictory coding -- ValidateHierarchy() enforces the converse +// direction (no applicable secondary where the primary codes absence) but +// not this one -- and whatever a validator decides to do about it, the score +// must not invent a state for it (T-396). Before this, such a tip +// kept a single concrete inapp state that contributed tie-break support and +// that the uppass could propagate INWARDS, resolving a node in the middle of +// the PRESENT region to it -- disjoint from every present neighbour in every +// secondary at once, so score_hierarchy_block() charged that branch d = m, +// the full alpha, for a node that by construction has no inapplicable +// secondaries. Measured over-charge was exactly alpha/m per affected +// branch, and it exceeded every concrete resolution as well as the missing +// treatment, so no reading of the data made it right. +// Because whether the old over-charge fired depended on the DELTRAN direction, +// it was also the dominant source of T-374's rooting-dependence. static int fitch_label_char( const TreeState& tree, const std::vector& tip_labels, @@ -150,6 +173,7 @@ static int fitch_label_char( int n_orig_chars, const std::vector& token_states, int n_levels, + uint32_t inapp_bit, const std::vector& pri_free, const CanonOrder& co, std::vector& state_sets) @@ -157,13 +181,20 @@ static int fitch_label_char( int n_tip = tree.n_tip; int n_node = tree.n_node; + // Applicable state set of this secondary at tip t: the states its token + // denotes, minus the inapplicable state, which is not one of them (above). + // Zero means the cell constrains nothing. + auto applicable_at = [&](int t) -> uint32_t { + return token_states[tip_labels[t * n_orig_chars + char_idx]] & ~inapp_bit; + }; + // The applicable domain: the states this character is observed in at tips // where it actually applies. Wildcarding to this rather than to all // n_levels bits keeps the tie-break arrays below as small as they were. uint32_t domain = 0; for (int t = 0; t < n_tip; ++t) { if (!pri_free[t]) { - domain |= token_states[tip_labels[t * n_orig_chars + char_idx]]; + domain |= applicable_at(t); } } // The character applies nowhere: it constrains nothing. Give every node one @@ -173,8 +204,8 @@ static int fitch_label_char( uint32_t used_mask = 0; std::vector observed(n_tip); for (int t = 0; t < n_tip; ++t) { - int label = tip_labels[t * n_orig_chars + char_idx]; - observed[t] = pri_free[t] ? domain : token_states[label]; + const uint32_t applicable = applicable_at(t); + observed[t] = (pri_free[t] || applicable == 0) ? domain : applicable; state_sets[t] = observed[t]; used_mask |= observed[t]; } @@ -401,7 +432,8 @@ static double score_hierarchy_block( std::vector buf(n_node); for (int j = 0; j < m; ++j) { fitch_label_char(tree, tip_labels, block.secondary_chars[j], - n_orig_chars, token_states, n_levels, pri_free, co, buf); + n_orig_chars, token_states, n_levels, inapp_bit, + pri_free, co, buf); for (int nd = 0; nd < n_node; ++nd) { sec_states[j * n_node + nd] = buf[nd]; } diff --git a/tests/testthat/test-CharacterHierarchy.R b/tests/testthat/test-CharacterHierarchy.R index 96a483232..17cd66d63 100644 --- a/tests/testthat/test-CharacterHierarchy.R +++ b/tests/testthat/test-CharacterHierarchy.R @@ -152,7 +152,7 @@ test_that(".NonHierarchyWeights subtracts hierarchy chars", { h <- CharacterHierarchy("1" = 2L) w_orig <- attr(ds, "weight") - w_adj <- .NonHierarchyWeights(ds, h) + w_adj <- TreeSearch:::.NonHierarchyWeights(ds, h) # Adjusted weights should be non-negative @@ -171,7 +171,7 @@ test_that(".BuildTipLabels creates correct matrix", { ds <- phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"), ambiguity = "?") - tl <- .BuildTipLabels(ds) + tl <- TreeSearch:::.BuildTipLabels(ds) expect_equal(nrow(tl), 3L) expect_equal(ncol(tl), 3L) # Values should be 0-based token indices @@ -180,7 +180,7 @@ test_that(".BuildTipLabels creates correct matrix", { test_that(".HierarchyToBlocks converts to 0-based flat list", { h <- CharacterHierarchy("1" = 2:4, "5" = 6:7) - blocks <- .HierarchyToBlocks(h) + blocks <- TreeSearch:::.HierarchyToBlocks(h) expect_length(blocks, 2) expect_equal(blocks[[1]]$primary, 0L) expect_equal(blocks[[1]]$secondaries, 1:3) @@ -190,7 +190,7 @@ test_that(".HierarchyToBlocks converts to 0-based flat list", { test_that(".HierarchyToBlocks flattens nested hierarchies", { h <- CharacterHierarchy("1" = list(2, 4, "3" = 9:10)) - blocks <- .HierarchyToBlocks(h) + blocks <- TreeSearch:::.HierarchyToBlocks(h) expect_gte(length(blocks), 2) # First block: primary=0, secondaries should include 1 and 3 (chars 2 and 4) expect_equal(blocks[[1]]$primary, 0L) @@ -213,7 +213,7 @@ test_that(".NonHierarchyWeights preserves non-hierarchy patterns", { h <- CharacterHierarchy("1" = 2L) idx <- attr(ds, "index") w_orig <- attr(ds, "weight") - w_adj <- .NonHierarchyWeights(ds, h) + w_adj <- TreeSearch:::.NonHierarchyWeights(ds, h) # Character 3 is not in the hierarchy; its pattern should keep its weight # unless it shares a pattern with a hierarchy character @@ -225,3 +225,127 @@ test_that(".NonHierarchyWeights preserves non-hierarchy patterns", { expect_gte(w_adj[pat], 0L) } }) + + +# ========================================================================= +# Nested hierarchies: the documented example must actually validate (T-395) +# ========================================================================= +# .ParseOneBlock() deliberately records a sub-controller BOTH as a dependent of +# its parent and as the controlling character of its own block -- that dual role +# is the entire content of "nested". .ValidateBlock() used to count the second +# occurrence as a double claim, so every nested hierarchy that could be written +# was rejected, including this file's own roxygen example. R CMD check could not +# see it because CharacterHierarchy() itself never validates: the example +# constructs fine and is never scored. + +test_that("a nested hierarchy validates, and the documented example works", { + # Char 1 controls {2, 3, 4, 5}; char 3 additionally controls {9, 10}. + nested <- CharacterHierarchy("1" = list(2, 3, 4, 5, "3" = 9:10)) + + expect_equal(nested[[1]]$controlling, 1L) + expect_equal(sort(nested[[1]]$dependents), c(2L, 3L, 4L, 5L)) + expect_equal(length(nested[[1]]$children), 1L) + expect_equal(nested[[1]]$children[[1]]$controlling, 3L) + expect_equal(nested[[1]]$children[[1]]$dependents, 9:10) + + # Coding invariants: secondaries "-" where their controller codes absence; + # char 3 binary where it applies; chars 9-10 "-" wherever char 3 is not "1". + mat <- rbind( + t1 = c("1", "0", "1", "0", "1", "0", "1", "0", "0", "1"), + t2 = c("1", "1", "1", "1", "0", "1", "0", "1", "1", "1"), + t3 = c("1", "0", "0", "1", "1", "0", "1", "1", "-", "-"), + t4 = c("0", "-", "-", "-", "-", "1", "0", "0", "-", "-"), + t5 = c("0", "-", "-", "-", "-", "0", "1", "1", "-", "-"), + t6 = c("0", "-", "-", "-", "-", "1", "1", "0", "-", "-")) + ds <- phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"), + ambiguity = "?") + + expect_silent(ValidateHierarchy(nested, ds)) + + # A genuine double claim must STILL be rejected -- the fix must not have + # simply disabled the check. Char 2 is claimed by both blocks here. + expect_error( + ValidateHierarchy(CharacterHierarchy("1" = 2:3, "5" = c(2L, 6L)), ds), + "multiple hierarchy blocks" + ) + # And a nested block's own dependents are still checked against other blocks. + expect_error( + ValidateHierarchy(CharacterHierarchy("1" = list(2, "3" = 9:10), "5" = 9L), + ds), + "multiple hierarchy blocks" + ) +}) + +test_that("a nested hierarchy scores, and stays rooting-invariant", { + nested <- CharacterHierarchy("1" = list(2, 3, 4, 5, "3" = 9:10)) + mat <- rbind( + t1 = c("1", "0", "1", "0", "1", "0", "1", "0", "0", "1"), + t2 = c("1", "1", "1", "1", "0", "1", "0", "1", "1", "1"), + t3 = c("1", "0", "0", "1", "1", "0", "1", "1", "-", "-"), + t4 = c("0", "-", "-", "-", "-", "1", "0", "0", "-", "-"), + t5 = c("0", "-", "-", "-", "-", "0", "1", "1", "-", "-"), + t6 = c("0", "-", "-", "-", "-", "1", "1", "0", "-", "-")) + ds <- phangorn::phyDat(mat, type = "USER", levels = c("-", "0", "1"), + ambiguity = "?") + tree <- Renumber(RenumberTips( + ape::read.tree(text = "(((t1,t2),t3),(t4,(t5,t6)));"), rownames(mat))) + + score <- TreeLength(tree, ds, hierarchy = nested, inapplicable = "hsj", + hsj_alpha = 1) + expect_true(is.finite(score)) + expect_gt(score, 0) + + # The HSJ score is a minimum over labellings of a sum of SYMMETRIC + # dissimilarities on an unrooted tree, so it cannot depend on the rooting. + rooted <- vapply(rownames(mat), function(tip) { + TreeLength(Renumber(RootTree(tree, tip)), ds, hierarchy = nested, + inapplicable = "hsj", hsj_alpha = 1) + }, double(1)) + expect_equal(diff(range(rooted)), 0) + + # The x-transformation genuinely does not implement nesting; its error must + # be the informative one, not a double-claim complaint from the validator. + expect_error(RecodeHierarchy(ds, nested), "[Nn]ested") +}) + +test_that("HierarchyFromNames detects nesting by tag extension", { + # Flat, documented case: unchanged. + flat <- HierarchyFromNames(c("sup_tail", "sub_tail_colour", "sub_tail_shape", + "sup_wing", "sub_wing_venation", "eyes")) + expect_equal(length(flat), 2L) + expect_equal(sort(flat[[1]]$dependents), c(2L, 3L)) + expect_equal(length(flat[[1]]$children), 0L) + + # Nested: `tail_tip` extends `tail`, so char 3 is both a dependent of char 1 + # and the controller of char 4. The old first-component tag match collapsed + # every depth onto the outermost tag, so this returned two flat blocks -- one + # of them empty -- and the nesting was silently lost. + nested <- HierarchyFromNames(c("sup_tail", "sub_tail_colour", + "sup_tail_tip", "sub_tail_tip_gloss")) + expect_equal(length(nested), 1L) + expect_equal(nested[[1]]$controlling, 1L) + expect_equal(sort(nested[[1]]$dependents), c(2L, 3L)) + expect_equal(length(nested[[1]]$children), 1L) + expect_equal(nested[[1]]$children[[1]]$controlling, 3L) + expect_equal(nested[[1]]$children[[1]]$dependents, 4L) + + # Longest-match, not first-match: the deeper dependent must not be captured + # by the shallower tag. + expect_false(4L %in% nested[[1]]$dependents) + + # A shared prefix without an underscore boundary is NOT nesting. + fin <- HierarchyFromNames(c("sup_tail", "sup_tailfin", "sub_tailfin_x")) + expect_equal(length(fin), 2L) + + # Three levels deep. + deep <- HierarchyFromNames(c("sup_a", "sub_a_x", "sup_a_b", "sub_a_b_y", + "sup_a_b_c", "sub_a_b_c_z")) + expect_equal(length(deep), 1L) + expect_equal(deep[[1]]$children[[1]]$controlling, 3L) + expect_equal(deep[[1]]$children[[1]]$children[[1]]$controlling, 5L) + expect_equal(deep[[1]]$children[[1]]$children[[1]]$dependents, 6L) + + # An orphan sub_ still warns. + expect_warning(HierarchyFromNames(c("sup_tail", "sub_nose_shape")), + "no corresponding sup_") +}) diff --git a/tests/testthat/test-ts-hsj.R b/tests/testthat/test-ts-hsj.R index c25c60470..a4d31a483 100644 --- a/tests/testthat/test-ts-hsj.R +++ b/tests/testthat/test-ts-hsj.R @@ -834,14 +834,32 @@ test_that("HSJ secondary '?' obeys the resolution invariant (T-375)", { # collapses the whole score to the secondary's ordinary Fitch step count. # # Tree ((t1,t2),(t3,t4)); primaries all "1"; secondary t1=t2=t3="0", t4 - # varies. Hand-derived: t4="0" ties all four -> 0 steps. t4="1" or t4="-" - # each disagree with the (t3,t4) clade's neighbour -> 1 step (the downpass - # intersect((t3=0),(t4=1 or -)) is empty, forcing a union). t4="?" must - # resolve to whichever concrete state is compatible AND cheapest -- here - # that's "0" (matching t1/t2/t3), giving 0 steps, so score("?") == 0 == - # min(0, 1, 1). Before the fix, fitch_label_char() bit-encoded the "?" - # TOKEN index as its own concrete state bit, indistinguishable from a - # genuine mismatch, and scored 1 -- violating the invariant (1 > 0). + # varies. Hand-derived: t4="0" ties all four -> 0 steps. t4="1" disagrees with + # the (t3,t4) clade's neighbour -> 1 step (the downpass intersect((t3=0), + # (t4=1)) is empty, forcing a union). t4="?" must resolve to whichever concrete + # state is compatible AND cheapest -- here that's "0" (matching t1/t2/t3), + # giving 0 steps. Before the fix, fitch_label_char() bit-encoded the "?" TOKEN + # index as its own concrete state bit, indistinguishable from a genuine + # mismatch, and scored 1 -- violating the invariant (1 > 0). + # + # t4="-" MUST ALSO SCORE 0, and this expectation CHANGED on 2026-08-03 (T-396): + # it previously asserted 1, reasoning that "-" disagrees with its neighbour + # just as "1" does. That reasoning treats the inapplicable token as an ordinary + # third state, which is precisely the behaviour Hopkins & St John (2021) + # introduce HSJ to avoid -- "Treating inapplicable characters as a new, + # separate state[] will ... skew the analysis, because having a new separate + # state increases the dissimilarity of all pairwise comparisons" (p.5). The + # paper's d counts "nonmatching secondary characters" among those that APPLY, + # and where secondaries are inapplicable to a taxon "they have no influence on + # the estimated dissimilarity" (p.5). An inapplicable secondary therefore + # contributes nothing to d whatever its controlling primary codes. + # + # Note every primary here is "1" (present), so this is the CONTRADICTORY coding + # case: "-" in a secondary whose controlling primary says the structure exists. + # The old value was not a defensible alternative reading -- it exceeded every + # concrete resolution AND the missing-data treatment, so no resolution of the + # cell produced it. Whether such a matrix should be rejected by + # ValidateHierarchy() is a separate question; the score must not invent a state. h <- CharacterHierarchy("1" = 2L) tree <- Renumber(RenumberTips( ape::read.tree(text = "((t1,t2),(t3,t4));"), paste0("t", 1:4))) @@ -854,8 +872,13 @@ test_that("HSJ secondary '?' obeys the resolution invariant (T-375)", { } scores <- vapply(c("0", "1", "-", "?"), score_for, double(1)) - expect_equal(unname(scores), c(0, 1, 1, 0)) + expect_equal(unname(scores), c(0, 1, 0, 0)) expect_lte(scores[["?"]], min(scores[c("0", "1", "-")])) + # An inapplicable secondary must cost no more than the cheapest concrete + # resolution of that cell -- the paper's "no influence" property. Asserted + # separately from the exact vector so a future change to the matrix still + # checks the property rather than only the number. + expect_lte(scores[["-"]], min(scores[c("0", "1")])) })