diff --git a/NEWS.md b/NEWS.md index 5b14f97a9..4a8994670 100644 --- a/NEWS.md +++ b/NEWS.md @@ -28,6 +28,28 @@ rooting-invariance is a property the method requires rather than a convention to pick. +- `inapplicable = "hsj"` scoring fixed an index-space confusion that could + under- or over-count the controlling primary's gains and losses, and could + score an ambiguous (`"?"`) secondary character as though it were a specific, + conflicting state. Internally, a tip's data value was read as an index into + the dataset's character *states*, but it is actually an index into the + dataset's observed *tokens* (a distinct, dataset-specific ordering that only + sometimes coincides with state order) -- so, depending on a dataset's + internal token ordering, an absent or inapplicable primary could be scored + as present, a present primary as absent, and a genuinely ambiguous secondary + character as a forced, arbitrary state. This was independent of tree + topology, so no search or comparison using `inapplicable = "hsj"` was + reliable: the same dataset and tree could report different scores merely by + virtue of the order characters happened to appear in. + + **HSJ scores may therefore differ from previous versions**, in either + direction: scores typically rise where a genuinely present or absent + controlling primary is now always counted, but can also fall where an + ambiguous secondary is no longer forced into a spurious mismatch. This is + a correctness fix to how a tip's data is looked up; it does not touch the + known rooting-sensitivity of HSJ scoring, which remains a separate, open + issue. + - `MaximizeParsimony(effort = )` replaces `strategy = `, which is removed (it was never released). `effort` is a **relative** offset, not an absolute level: `0` (the default) accepts the amount of search the dataset's size and diff --git a/R/CharacterHierarchy.R b/R/CharacterHierarchy.R index ab708f5a4..219fbed0c 100644 --- a/R/CharacterHierarchy.R +++ b/R/CharacterHierarchy.R @@ -395,7 +395,9 @@ HierarchyControlling <- function(hierarchy) { # Build the tip-labels matrix for HSJ scoring. # # Converts a phyDat dataset into an integer matrix of per-tip, per-character -# state labels (0-based) for the C++ HSJ scorer: length(dataset) rows (tips) by +# TOKEN labels (0-based indices into `attr(dataset, "allLevels")`, NOT state +# indices into `attr(dataset, "levels")`; the two only sometimes coincide -- +# see T-375/T-376) for the C++ HSJ scorer: length(dataset) rows (tips) by # length(attr(dataset, "index")) columns (original characters). .BuildTipLabels <- function(dataset) { idx <- attr(dataset, "index") @@ -415,8 +417,12 @@ HierarchyControlling <- function(hierarchy) { # Identify the primary "absent" state for HSJ scoring. # -# Returns the 0-based token index of the controlling primary character's -# *absent* state, for the C++ HSJ scorer's `absent_state` argument. +# Returns the 0-based STATE (`levels`) index of the controlling primary +# character's *absent* state, for the C++ HSJ scorer's `absent_state` +# argument -- NOT a token/`allLevels` index (T-375/T-376: `tip_labels`, built +# by .BuildTipLabels() below, holds token indices, a different index space; +# the C++ kernel translates a tip's token into this state space via +# `DataSet::token_states` before comparing it to this value). # # Under reductive coding (Hopkins & St John 2021) the primary codes a # structure's presence/absence, conventionally "0" = absent, "1" = present. diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 685938561..ba3cc4bfa 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -1536,8 +1536,9 @@ MaximizeParsimony <- function( hsjArgs$hierarchyBlocks <- .HierarchyToBlocks(hierarchy) hsjArgs$hsjTipLabels <- .BuildTipLabels(dataset) hsjArgs$hsjAlpha <- as.double(hsj_alpha) - # 0-based token index of the primary's "absent" state (depends on level - # ordering, so computed from the data rather than hard-coded). + # 0-based STATE (levels) index of the primary's "absent" state (depends on + # level ordering, so computed from the data rather than hard-coded; the + # C++ kernel translates tip token labels into this same index space). hsjArgs$hsjAbsentState <- .HSJAbsentState(dataset) # Adjust weights: subtract hierarchy characters so Fitch scores non-hierarchy diff --git a/dev/red-team/heavy-tests/hsj-paper-oracle.R b/dev/red-team/heavy-tests/hsj-paper-oracle.R index 596f6e4c1..32e65a532 100644 --- a/dev/red-team/heavy-tests/hsj-paper-oracle.R +++ b/dev/red-team/heavy-tests/hsj-paper-oracle.R @@ -154,6 +154,13 @@ if (isTRUE(all.equal(alphaLive, 0))) { cat(" across its 10 tip-rootings while alpha=0 gives 20 ten times, and\n") cat(" XFORM is rooting-dependent on 165/300 random topologies (T-374).\n") } else { + cat(" NOTE a PASS below is NOT evidence that T-374 is fixed -- it only means\n") + cat(" this small, Fig.1-derived matrix/topology doesn't happen to trigger\n") + cat(" fitch_label_char()'s rooting-sensitivity. T-374's rooting bug is\n") + cat(" data/topology-dependent (measured 165-197/300 random cases) and can\n") + cat(" remain OPEN and readily reproducible on OTHER matrices even when every\n") + cat(" check below is green. T-374's fix is the paper's two-state DP (or\n") + cat(" marginal-MPR resolution) -- do not close T-374 on this check alone.\n") for (alpha in c(0, 0.5, 1)) { scores <- vapply(rownames(extended), function(tip) { TreeLength(RootTree(base, tip), extDs, hierarchy = hierarchy, diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index b68f00fc5..fdd7d6775 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -498,9 +498,9 @@ result <- MaximizeParsimony( result attr(result, "score") -# Ask for one notch more search than this dataset would get by default, +# Ask for one notch less search than this dataset would get by default, # whatever its size: -harder <- MaximizeParsimony(dataset, effort = 1L, maxReplicates = 12L) +sprint <- MaximizeParsimony(dataset, effort = -1, maxReplicates = 12) } \references{ diff --git a/src/ts_data.cpp b/src/ts_data.cpp index cdd45d72e..50bf960c4 100644 --- a/src/ts_data.cpp +++ b/src/ts_data.cpp @@ -52,6 +52,11 @@ DataSet build_dataset( } } } + // Stored verbatim (pre-simplification, original global state space) for + // HSJ scoring, which must translate tip_labels' token indices into this + // same state space (T-375/T-376) -- do not apply state_remap to it. + ds.token_states = token_states; + ds.n_levels = n_states; // --- Character simplification --- SimplificationResult simpl = simplify_patterns( diff --git a/src/ts_data.h b/src/ts_data.h index 8c048dc45..5c078fda2 100644 --- a/src/ts_data.h +++ b/src/ts_data.h @@ -68,7 +68,12 @@ struct HierarchyBlock { int primary_char; // original character index (0-based) std::vector secondary_chars; // original character indices (0-based) int n_secondaries; // = secondary_chars.size() - int absent_state; // state index meaning "absent" in primary + // State (levels) index meaning "absent" in primary -- NOT a token/allLevels + // index. tip_labels (see DataSet::tip_labels) holds token indices, so a + // tip's label must be translated via DataSet::token_states before it can be + // compared to this field or to DataSet::inapp_state (T-375/T-376: the two + // index spaces were previously compared directly via `==`). + int absent_state = -1; }; struct CharBlock { @@ -166,10 +171,26 @@ struct DataSet { // Populated by build_dataset(); used by HSJ scoring. int inapp_state = -1; + // Per-token (contrast-row / allLevels) state-set bitmask: bit s of + // token_states[t] is set iff the phyDat contrast matrix has contrast[t, s] + // > 0.5 (token t is compatible with state s). Populated by build_dataset() + // directly from the ORIGINAL, pre-simplification contrast matrix (T-375/ + // T-376) -- do not apply state_remap to it. This is the translation HSJ + // scoring needs: tip_labels (below) holds TOKEN indices, but absent_state/ + // inapp_state are STATE indices, so a tip's label must be looked up here + // before it can be compared to either. n_levels is the number of columns + // (states) these bitmasks range over -- deliberately not named n_states, + // which on CharBlock means something different (per-block, post- + // simplification). + std::vector token_states; + int n_levels = 0; + // HSJ scoring data (populated when scoring_mode == HSJ). // These are set by the Rcpp bridge after build_dataset(). std::vector hierarchy_blocks; - // tip_labels: per-tip per-original-char state labels (0-based). + // tip_labels: per-tip per-original-char state labels (0-based TOKEN index, + // i.e. an index into token_states above -- NOT directly comparable to + // absent_state/inapp_state; see token_states' comment). // Layout: tip_labels[tip * n_orig_chars + char] std::vector tip_labels; int n_orig_chars = 0; diff --git a/src/ts_hsj.cpp b/src/ts_hsj.cpp index f9a5025db..5d70e7f30 100644 --- a/src/ts_hsj.cpp +++ b/src/ts_hsj.cpp @@ -32,29 +32,36 @@ std::vector partition_weights( // single state so that parent–child mismatches can be detected for HSJ // secondary dissimilarity. // Returns number of Fitch steps (union operations in the downpass). +// +// tip_labels holds 0-based TOKEN (allLevels/contrast-row) indices, not state +// indices (T-375): a token like "?" is its own row in the contrast matrix, +// generally with several columns set, so treating the token index itself as +// a bit position (as this function formerly did) is a category error -- +// `label > 30` can never fire on a valid token index, so "?" silently scored +// as one concrete, arbitrary state instead of the wildcard it denotes. +// token_states[label] gives the actual bitmask of states the token is +// compatible with (populated by build_dataset() from the contrast matrix), +// which is what state_sets must hold to make the Fitch downpass/uppass below +// correct for ambiguous tokens. inapp_state needs no special handling here -- +// per this module's design (see feature-inapplicable.md), the inapplicable +// state is deliberately just another concrete state in the secondary pass. static int fitch_label_char( const TreeState& tree, const std::vector& tip_labels, int char_idx, int n_orig_chars, - int inapp_state, + const std::vector& token_states, + int n_levels, std::vector& state_sets) { int n_tip = tree.n_tip; int n_node = tree.n_node; - // Initialize tips; track the number of distinct states (highest concrete - // token index + 1) so the order-invariant tie-break arrays can be sized. - int n_states = 1; + uint32_t used_mask = 0; for (int t = 0; t < n_tip; ++t) { int label = tip_labels[t * n_orig_chars + char_idx]; - if (label < 0 || label > 30) { - // Ambiguous: all states - state_sets[t] = 0xFFFFFFFFu; - } else { - state_sets[t] = 1u << label; - if (label + 1 > n_states) n_states = label + 1; - } + state_sets[t] = token_states[label]; + used_mask |= state_sets[t]; } // --- Downpass --- @@ -77,27 +84,43 @@ static int fitch_label_char( // --- Order-invariant tie-break support ------------------------------- // The uppass below must resolve ambiguous nodes to a single state so that // parent-child mismatches (the HSJ secondary dissimilarity) can be counted. - // Resolving by bit index (the old "lowest set bit") makes the result depend - // on the arbitrary phyDat `levels` ordering, because which token maps to the - // lowest bit is determined by `levels`. Instead we resolve toward the token - // with the most support in the node's own subtree, breaking ties by the - // smallest supporting tip index. Both keys are properties of the *tokens* - // and the tree, not of the bit encoding, so the resolution — and hence the - // mismatch count — is invariant to level ordering. This still yields a valid - // most-parsimonious reconstruction, so the dissimilarity stays non-zero (the - // concern that motivated adding the uppass in the first place). + // Resolving by bit index (the old "lowest set bit") would make the result + // depend on the arbitrary phyDat `levels` ordering, since which state + // occupies the lowest bit is determined by `levels`. Instead we resolve + // toward the state with the most support in the node's own subtree, + // breaking ties by the smallest supporting tip index. Both keys are + // properties of the *states* and the tree, not of the bit encoding, so the + // resolution — and hence the mismatch count — is invariant to level + // ordering. This still yields a valid most-parsimonious reconstruction, so + // the dissimilarity stays non-zero (the concern that motivated adding the + // uppass in the first place). // - // tb_cnt[node * K + s] = # tips in subtree(node) carrying concrete token s - // tb_mintip[node * K + s] = smallest tip index in subtree(node) with token s - const int K = n_states; + // tb_cnt[node * K + s] = # tips in subtree(node) carrying concrete state s + // tb_mintip[node * K + s] = smallest tip index in subtree(node) with state s + // + // K only needs to cover states this CHARACTER actually uses -- the downpass + // above only ever intersects/unions existing tip bits, so no bit outside + // `used_mask` (accumulated while initializing state_sets) can appear at any + // node. Bounding K by `used_mask`'s highest set bit (rather than the + // dataset-global n_levels) keeps these arrays as small as the old + // per-character sizing did, even when other characters in the dataset use + // many more states than this one. + int K = 0; + for (uint32_t m = used_mask; m; m >>= 1) ++K; + if (K == 0) K = 1; + if (K > n_levels) K = n_levels; // defensive: never exceed the real state space const int INF_TIP = std::numeric_limits::max(); std::vector tb_cnt(static_cast(n_node) * K, 0); std::vector tb_mintip(static_cast(n_node) * K, INF_TIP); for (int t = 0; t < n_tip; ++t) { - int label = tip_labels[t * n_orig_chars + char_idx]; - if (label >= 0 && label < K) { // concrete (ambiguous tips favour nothing) - tb_cnt[static_cast(t) * K + label] = 1; - tb_mintip[static_cast(t) * K + label] = t; + uint32_t set = state_sets[t]; + // Only a tip resolved to exactly one concrete state contributes tie-break + // support -- an ambiguous tip (e.g. "?", now correctly multi-bit) must + // not bias which state the uppass prefers, same as before this fix. + if (set != 0 && (set & (set - 1)) == 0) { + int s = ctz64(set); + tb_cnt[static_cast(t) * K + s] = 1; + tb_mintip[static_cast(t) * K + s] = t; } } for (int i = 0; i < static_cast(tree.postorder.size()); ++i) { @@ -115,9 +138,9 @@ static int fitch_label_char( } // Resolve `state_sets[node]` to a single state, preferring the best-supported - // token (max subtree count; ties broken by smallest supporting tip index — - // a strict order, since distinct tokens never share a supporting tip). When - // no token in the set has concrete support (the whole subtree is ambiguous + // state (max subtree count; ties broken by smallest supporting tip index — + // a strict order, since distinct states never share a supporting tip). When + // no state in the set has concrete support (the whole subtree is ambiguous // for this character), fall back to the lowest set bit: every node then // inherits it and no mismatch is affected, so the choice is score-neutral. auto pick_state = [&](int node) -> uint32_t { @@ -186,7 +209,9 @@ static double score_hierarchy_block( double alpha, const std::vector& tip_labels, int n_orig_chars, - int inapp_state) + int inapp_state, + const std::vector& token_states, + int n_levels) { const int n_tip = tree.n_tip; const int n_node = tree.n_node; @@ -200,25 +225,35 @@ 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, inapp_state, buf); + n_orig_chars, token_states, n_levels, buf); for (int nd = 0; nd < n_node; ++nd) { sec_states[j * n_node + nd] = buf[nd]; } } } - // Step 2: Determine primary state at each tip - // primary_present[tip] = true unless the primary codes the structure as - // absent. The structure is absent when the primary token is the explicit - // "absent" state (block.absent_state, e.g. "0") OR the inapplicable token - // ("-"). This mirrors the x-transform recoding (recode_hierarchy.R), which - // treats `pri == "0" || pri == "-"` as absent, and is required for nested - // hierarchies where a controlling primary may itself be inapplicable. - std::vector primary_present(n_tip, false); - for (int t = 0; t < n_tip; ++t) { - int label = tip_labels[t * n_orig_chars + block.primary_char]; - primary_present[t] = (label != block.absent_state) && (label != inapp_state); - } + // Step 2: Determine primary feasibility at each tip via SET MEMBERSHIP. + // `label` is a TOKEN (allLevels/contrast-row) index, but block.absent_state + // and inapp_state are STATE (levels) indices -- comparing them directly via + // `==` (the former code) was T-375/T-376's bug: it happened to work only + // when `levels` and `allLevels` coincide, which is common but not + // guaranteed, and it always mis-scored "?" (whose token index is never + // itself a valid comparison target for either state index). token_states + // translates the token into the state-space bitmask it actually denotes, so + // both sides of the test are in the same index space. + // + // The structure is absent when the primary's state set includes the + // explicit "absent" state (block.absent_state, e.g. "0") OR the + // inapplicable state ("-"). This mirrors the x-transform recoding + // (recode_hierarchy.R), which treats `pri == "0" || pri == "-"` as absent, + // and is required for nested hierarchies where a controlling primary may + // itself be inapplicable. a(leaf) = 0 iff the token's state set meets the + // absent-coding states; p(leaf) = 0 iff it meets the present-coding states + // (anything else). A fully ambiguous token (e.g. "?") meets both, so + // a = p = 0 there: genuinely unconstrained, not "present" (the old code's + // effective classification of "?"). + const uint32_t inapp_bit = (inapp_state >= 0) ? (1u << inapp_state) : 0u; + const uint32_t absent_bits = (1u << block.absent_state) | inapp_bit; // Step 3: a(n)/p(n) DP // a[node], p[node] @@ -227,13 +262,10 @@ static double score_hierarchy_block( // Initialize leaves for (int t = 0; t < n_tip; ++t) { - if (primary_present[t]) { - a[t] = INF; - p[t] = 0.0; - } else { - a[t] = 0.0; - p[t] = INF; - } + int label = tip_labels[t * n_orig_chars + block.primary_char]; + uint32_t set = token_states[label]; + a[t] = (set & absent_bits) ? 0.0 : INF; + p[t] = (set & ~absent_bits) ? 0.0 : INF; } // Helper: count secondary mismatches between two nodes @@ -317,7 +349,8 @@ double hsj_score( double hsj_total = 0.0; for (const auto& block : hierarchy_blocks) { hsj_total += score_hierarchy_block( - tree, block, alpha, tip_labels, n_orig_chars, ds.inapp_state); + tree, block, alpha, tip_labels, n_orig_chars, ds.inapp_state, + ds.token_states, ds.n_levels); } return fitch_total + hsj_total; diff --git a/src/ts_hsj.h b/src/ts_hsj.h index 737fe49fe..3c9c2baac 100644 --- a/src/ts_hsj.h +++ b/src/ts_hsj.h @@ -51,7 +51,10 @@ std::vector partition_weights( // ds: dataset (used for non-hierarchy characters AND hierarchy data) // hierarchy_blocks: hierarchy specification // alpha: HSJ scaling parameter in [0, 1] -// tip_labels: per-tip, per-original-char state labels (0-based token index). +// tip_labels: per-tip, per-original-char labels (0-based TOKEN/allLevels +// index -- NOT a state index; translate via ds.token_states before +// comparing to a HierarchyBlock's absent_state or to ds.inapp_state, both +// of which are state indices; see T-375/T-376 and DataSet::token_states). // Layout: tip_labels[tip * n_orig_chars + char]. This is the full // (uncompressed) original matrix needed for secondary character matching. // n_orig_chars: number of original characters (before compression) diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 914314b1d..fe0a17639 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -54,6 +54,33 @@ void validate_tip_data_values(const int* tip_data_r, int n_tips, } } +// Validate `tip_labels`/`absent_state` VALUES at the Rcpp boundary: the HSJ +// kernel treats each tip_labels entry as a 0-based index into +// `DataSet::token_states` (size n_tokens) and `absent_state` as a bit +// position within an `n_levels`-bit state-space mask, with no further +// checking (T-375/T-376: before that fix, a bad value was merely a wrong +// scalar comparison; now it is an out-of-bounds read / undefined shift). +// Public wrappers always derive both from a validated phyDat via +// `.BuildTipLabels()`/`.HSJAbsentState()`, so this only guards a direct +// internal call (`TreeSearch:::`) with hand-crafted values -- but that +// includes both `ts_hsj_score()` (the test bridge, below) and +// `ts_driven_search()`/`ts_collapse_pool()` via their `hsjConfig` argument +// (see `unpack_hsj()`), since `R/ts-driven-compat.R`'s legacy wrapper exposes +// `hsjTipLabels`/`hsjAbsentState` as independent, hand-settable arguments. +void validate_hsj_tip_labels(const IntegerMatrix& tip_labels_r, + int absent_state, int n_tokens, int n_levels) { + if (absent_state < 0 || absent_state >= n_levels) { + Rcpp::stop("`absent_state` must be in [0, %d); found %d", + n_levels, absent_state); + } + for (int v : tip_labels_r) { + if (v < 0 || v >= n_tokens) { + Rcpp::stop("`tip_labels` values must be in [0, nrow(contrast)) (%d); " + "found %d", n_tokens, v); + } + } +} + // Validate an `addition_order` VALUE vector at the Rcpp boundary: // wagner_tree() treats a non-empty `order` as a length-n_tips permutation of // 0..n_tips-1 and reads order[0..2] / order[i] for i in [3, n_tips) with no @@ -1822,6 +1849,9 @@ static void unpack_hsj(Nullable hsjConfig, ts::DataSet& ds) { if (hc.containsElementNamed("hsjTipLabels") && !Rf_isNull(hc["hsjTipLabels"])) { IntegerMatrix tl = as(hc["hsjTipLabels"]); + validate_hsj_tip_labels(tl, hsjAbsentState, + static_cast(ds.token_states.size()), + ds.n_levels); int n_t = tl.nrow(); int n_c = tl.ncol(); ds.n_orig_chars = n_c; @@ -3122,6 +3152,9 @@ double ts_hsj_score( IntegerMatrix tip_labels_r, int absent_state) { + validate_hsj_tip_labels(tip_labels_r, absent_state, contrast.nrow(), + contrast.ncol()); + // Build DataSet for non-hierarchy characters (weight already adjusted) ts::DataSet ds = make_dataset(contrast, tip_data, weight, levels); diff --git a/src/ts_sector.cpp b/src/ts_sector.cpp index 2c438b771..25452eec2 100644 --- a/src/ts_sector.cpp +++ b/src/ts_sector.cpp @@ -1540,8 +1540,9 @@ SectorResult rss_search(TreeState& tree, DataSet& ds, result.total_steps_saved = 0; // build_reduced_dataset() does not copy hierarchy_blocks, tip_labels, - // n_orig_chars, hsj_alpha, or sankoff_* fields (T-303). Sector-internal - // scoring would silently degrade to Fitch-only. Same class as T-275 guard. + // n_orig_chars, hsj_alpha, token_states, n_levels, or sankoff_* fields + // (T-303). Sector-internal scoring would silently degrade to Fitch-only. + // Same class as T-275 guard. if (ds.scoring_mode == ScoringMode::HSJ || ds.scoring_mode == ScoringMode::XFORM) { return result; @@ -1908,8 +1909,9 @@ SectorResult xss_search(TreeState& tree, DataSet& ds, result.total_steps_saved = 0; // build_reduced_dataset() does not copy hierarchy_blocks, tip_labels, - // n_orig_chars, hsj_alpha, or sankoff_* fields (T-303). Sector-internal - // scoring would silently degrade to Fitch-only. Same class as T-275 guard. + // n_orig_chars, hsj_alpha, token_states, n_levels, or sankoff_* fields + // (T-303). Sector-internal scoring would silently degrade to Fitch-only. + // Same class as T-275 guard. if (ds.scoring_mode == ScoringMode::HSJ || ds.scoring_mode == ScoringMode::XFORM) { return result; diff --git a/tests/testthat/test-ts-hsj.R b/tests/testthat/test-ts-hsj.R index 179d81935..d6ef318c0 100644 --- a/tests/testthat/test-ts-hsj.R +++ b/tests/testthat/test-ts-hsj.R @@ -26,8 +26,8 @@ hsj_score <- function(tree, dataset, hierarchy, alpha = 1.0) { nrow = length(dataset), byrow = TRUE) blocks <- .HierarchyToBlocks(hierarchy) tl <- .BuildTipLabels(dataset) - # absent_state = 0-based token index of "0" (= 1 for levels c("-","0","1")), - # computed the same way the driven pipeline does. + # absent_state = 0-based STATE (levels) index of "0" (= 1 for levels + # c("-","0","1")), computed the same way the driven pipeline does. ts_hsj_score( edge = tree$edge, contrast = at$contrast, @@ -577,7 +577,7 @@ test_that("HSJ handles extreme absent/present ratios", { # (Driven pipeline previously hard-coded 0L = index of "-", so primaries # coded "0" were treated as present and gain/loss was never counted.) # ========================================================================= -test_that(".HSJAbsentState() tracks the '0' token across level orderings", { +test_that(".HSJAbsentState() tracks the '0' state (levels index) across level orderings", { expect_equal(.HSJAbsentState(make_hsj_dat( matrix(c("0", "1", "0", "1"), 2, dimnames = list(c("a", "b"), NULL)), levels = c("-", "0", "1"))), 1L) @@ -703,6 +703,162 @@ test_that("HSJ secondary dissimilarity is level-order invariant (multistate)", { }) +# ========================================================================= +# Regression: T-375/T-376 -- tip_labels holds TOKEN (allLevels/contrast-row) +# indices, but the primary's absent_state/inapp_state (and the bitmask built +# by fitch_label_char() for secondaries) are STATE (levels) indices. +# `make_hsj_dat()`'s construction (phyDat(type = "USER", levels =, ambiguity +# = "?")) keeps allLevels and levels in the SAME relative order, so it cannot +# expose this bug -- every absolute-value test above using it would pass +# whether or not the two index spaces were confused. These tests instead +# permute the contrast-ROW order directly (holding levels/taxa/tree fixed), +# which is the only thing that moves tip_labels, and separately use +# MatrixToPhyDat(), whose allLevels is ordered by first appearance and so +# routinely disagrees with levels -- the same shape of misalignment as the +# shipped Vinther2008.nex dataset (see dev/red-team/findings.md T-376). +# ========================================================================= + +# Relabel the (arbitrary) contrast row order ONLY, exactly as in +# dev/red-team/heavy-tests/hsj-token-permutation.R. Taxa, tip numbering, tree, +# levels and every token's state set are untouched, so the dataset is +# identical -- any score change is not explicable by anything but the T-376 +# index-space bug. +.PermuteTokens <- function(d, perm) { + at <- attributes(d) + inv <- order(perm) + out <- lapply(unclass(d), function(x) inv[x]) + at$allLevels <- at$allLevels[perm] + at$contrast <- at$contrast[perm, , drop = FALSE] + attributes(out) <- at + out +} + +# HSJ score of `d` under every contrast-row permutation of its token +# alphabet, asserting each permuted dataset is byte-identical to `d` via +# PhyDatToMatrix() (the load-bearing confound-free check: only the arbitrary +# token order moves, nothing else). +.AllTokenOrderingScores <- function(d, tips, tree, h, alpha = 1) { + ref <- PhyDatToMatrix(d)[tips, , drop = FALSE] + nTok <- length(attr(d, "allLevels")) + perms <- as.matrix(expand.grid(rep(list(seq_len(nTok)), nTok))) + perms <- perms[apply(perms, 1, function(r) !anyDuplicated(r)), , drop = FALSE] + apply(perms, 1, function(perm) { + dp <- .PermuteTokens(d, perm) + stopifnot(identical(PhyDatToMatrix(dp)[tips, , drop = FALSE], ref)) + hsj_score(tree, dp, h, alpha = alpha) + }) +} + +test_that("HSJ score is invariant to contrast-row (token) order", { + tips <- paste0("t", 1:4) + tree <- Renumber(RenumberTips( + ape::read.tree(text = "((t1,t2),(t3,t4));"), tips)) + + # (1) Zero secondaries: isolates the primary-feasibility set-membership + # test in score_hierarchy_block() alone, with fitch_label_char() (T-375) + # never entered. + b1 <- rbind(t1 = c("-", "1"), t2 = c("0", "0"), + t3 = c("1", "?"), t4 = c("?", "-")) + scores1 <- .AllTokenOrderingScores( + make_hsj_dat(b1), tips, tree, CharacterHierarchy(`2` = integer(0))) + expect_equal(scores1, rep(scores1[[1]], length(scores1)), + info = "zero secondaries: isolates primary feasibility in score_hierarchy_block()") + + # (2) One secondary: adds fitch_label_char()'s token-to-state translation + # (T-375) on top of (1). + b2 <- rbind(t1 = c("-", "1", "0"), t2 = c("0", "1", "1"), + t3 = c("1", "0", "-"), t4 = c("?", "1", "?")) + scores2 <- .AllTokenOrderingScores( + make_hsj_dat(b2), tips, tree, CharacterHierarchy(`2` = 3L)) + expect_equal(scores2, rep(scores2[[1]], length(scores2)), + info = "one secondary: adds fitch_label_char() (T-375)") +}) + +test_that("HSJ handles a MatrixToPhyDat token/state misalignment (T-376)", { + # MatrixToPhyDat() orders allLevels by first appearance, which routinely + # disagrees with `levels` -- unlike make_hsj_dat() above. This specific + # matrix reproduces exactly the shape of misalignment found on the + # package's own shipped Vinther2008.nex dataset: levels = "- 0 1", but + # allLevels = "1 0 -" (confirmed by inspection below), so token("0") = 1 + # coincides with the STATE index of "0" only by luck of levels' order, and + # more importantly token("1") = 0 sits where inapp_state would be checked + # against under the old buggy code. + mat <- rbind( + t1 = c("1", "0"), t2 = c("1", "1"), t3 = c("1", "1"), + t4 = c("0", "-"), t5 = c("1", "0") + ) + ds <- MatrixToPhyDat(mat) + expect_equal(attr(ds, "levels"), c("-", "0", "1")) + expect_equal(attr(ds, "allLevels"), c("1", "0", "-")) + + h <- CharacterHierarchy("1" = 2L) + + # Hand-derived expected value: exactly one tip (t4) is absent among five + # present tips. Under the absent/present DP's symmetric branch costs + # (absent<->present = 1, absent<->absent = present<->present = 0), a binary + # character with a single differing tip always costs exactly 1 step, + # regardless of tree topology (there is always one branch, somewhere, that + # can carry the single change) -- so this is topology-invariant, checked + # across three unrelated topologies. Before the T-375/T-376 fix, the + # token/state confusion misclassified every one of these five tips as + # ABSENT (all read the same, wrongly), so the block cost 0 instead of 1: + # the controlling primary silently contributed nothing, exactly the T-376 + # escalation ("HSJ(alpha=0) == Fitch(non-controlling primaries only)"). + for (topo_txt in c("((t1,t2),(t3,(t4,t5)));", "(((t1,t4),t2),(t3,t5));", + "(t4,(t1,(t2,(t3,t5))));")) { + tree <- Renumber(RenumberTips(ape::read.tree(text = topo_txt), names(ds))) + expect_equal( + TreeLength(tree, ds, hierarchy = h, inapplicable = "hsj", hsj_alpha = 0), + 1, + info = topo_txt + ) + # The paper's alpha=0 identity (Hopkins & St John 2021, p.6): HSJ(alpha=0) + # must equal plain Fitch scoring of the primary character alone, + # controlling primary included. + expect_equal( + TreeLength(tree, ds, hierarchy = h, inapplicable = "hsj", hsj_alpha = 0), + TreeLength(tree, MatrixToPhyDat(mat[, 1, drop = FALSE])), + info = topo_txt + ) + } +}) + +test_that("HSJ secondary '?' obeys the resolution invariant (T-375)", { + # T-375's own acceptance criterion: score("?") <= min over concrete + # resolutions. The tests above (contrast-row permutation, alpha=0 on a + # misaligned layout) only exercise the T-376 primary_present term -- alpha=0 + # never calls count_mismatches(), so none of them can see whether + # fitch_label_char() resolves a "?" secondary correctly. This one isolates + # T-375 by keeping every primary "1" (present, so the primary DP is + # loss-free throughout and contributes nothing but a floor of 0), which + # 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). + h <- CharacterHierarchy("1" = 2L) + tree <- Renumber(RenumberTips( + ape::read.tree(text = "((t1,t2),(t3,t4));"), paste0("t", 1:4))) + + score_for <- function(t4sec) { + mat <- matrix(c("1", "0", "1", "0", "1", "0", "1", t4sec), + nrow = 4, byrow = TRUE, + dimnames = list(paste0("t", 1:4), NULL)) + hsj_score(tree, make_hsj_dat(mat), h, alpha = 1) + } + + scores <- vapply(c("0", "1", "-", "?"), score_for, double(1)) + expect_equal(unname(scores), c(0, 1, 1, 0)) + expect_lte(scores[["?"]], min(scores[c("0", "1", "-")])) +}) + + # ========================================================================= # Test: HSJ + sectorial search (T-303 guard) # ========================================================================= diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index d77b371f5..7d635ba84 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -80,6 +80,12 @@ The HSJ and x-transformation methods require a character hierarchy specification via the `CharacterHierarchy` class; see `vignette("inapplicable")` for details. +HSJ scoring correctly identifies the controlling primary's absence, +presence, and any partial ambiguity (e.g. `?`), regardless of the arbitrary +internal order in which a dataset's observed character-state tokens happen +to be stored -- that internal order need not match, and often does not +match, the order of the character states themselves. + ## Starting trees ### Wagner tree construction