From e2d3e0b1c9ca0590d1abdba13e927f67aa3fc75c Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 12:24:14 +0100 Subject: [PATCH 01/11] fix(hsj): translate token indices to state indices before comparison (T-375, T-376) HSJ's primary_present test and fitch_label_char() both compared tip_labels' token (allLevels/contrast-row) indices directly against state (levels) indices -- absent_state, inapp_state, and bit positions in state_sets are all in state space, but tip_labels is in token space. The two spaces coincide often enough that earlier fixes and their test batteries (which permuted levels, not allLevels) stayed green, but on real reductively-coded data (e.g. Vinther2008.nex) absent/inapplicable primaries could read as present while "?" read as absent, and the score was provably not a function of the data under contrast-row permutation alone. DataSet now stores token_states (the token -> state-set bitmask already computed transiently in build_dataset()) and n_levels, so the HSJ kernel can translate a token label into its state set before testing set membership against absent_state/inapp_state, and before building Fitch state_sets for secondary characters (fixing "?" scoring as a concrete, conflicting state instead of a wildcard). Also corrects the stale "token index" docstrings on .HSJAbsentState() and its caller, which described a levels index. --- R/CharacterHierarchy.R | 8 ++- R/MaximizeParsimony.R | 5 +- src/ts_data.cpp | 5 ++ src/ts_data.h | 25 ++++++- src/ts_hsj.cpp | 128 +++++++++++++++++++++-------------- src/ts_hsj.h | 5 +- tests/testthat/test-ts-hsj.R | 2 +- 7 files changed, 118 insertions(+), 60 deletions(-) diff --git a/R/CharacterHierarchy.R b/R/CharacterHierarchy.R index ab708f5a4..385acf3e2 100644 --- a/R/CharacterHierarchy.R +++ b/R/CharacterHierarchy.R @@ -415,8 +415,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/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..53390c1d7 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; }; 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..5f0892256 100644 --- a/src/ts_hsj.cpp +++ b/src/ts_hsj.cpp @@ -32,29 +32,34 @@ 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; 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]; } // --- Downpass --- @@ -77,27 +82,36 @@ 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 + const int K = n_levels; 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) { + for (int s = 0; s < K; ++s) { + if (set & (1u << s)) { + tb_cnt[static_cast(t) * K + s] = 1; + tb_mintip[static_cast(t) * K + s] = t; + break; + } + } } } for (int i = 0; i < static_cast(tree.postorder.size()); ++i) { @@ -115,9 +129,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 +200,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 +216,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 +253,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 +340,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/tests/testthat/test-ts-hsj.R b/tests/testthat/test-ts-hsj.R index 179d81935..c3bb2fd58 100644 --- a/tests/testthat/test-ts-hsj.R +++ b/tests/testthat/test-ts-hsj.R @@ -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) From 7812a5363d3a073889ec321e46bfff86464438b0 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 12:41:05 +0100 Subject: [PATCH 02/11] test(hsj): add regression coverage that permutes contrast rows, not levels (T-375, T-376) make_hsj_dat()'s construction keeps allLevels in the same relative order as levels, so none of the existing absolute-value tests could have caught the token/state index confusion -- confirmed by running this file's original tests unmodified against a pre-fix build (all still pass, unchanged). Adds: (1) contrast-row permutation invariance, lifted from dev/red-team/heavy-tests/hsj-token-permutation.R, for zero and one secondaries; (2) a MatrixToPhyDat()-emergent token/state misalignment (levels = "- 0 1", allLevels = "1 0 -", the same shape as the shipped Vinther2008.nex dataset) with a hand-derived absolute score, plus the paper's alpha=0-equals-Fitch identity on the same data. Verified both new tests fail against a throwaway pre-fix build of cpp-search and pass against the fix. --- tests/testthat/test-ts-hsj.R | 122 +++++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) diff --git a/tests/testthat/test-ts-hsj.R b/tests/testthat/test-ts-hsj.R index c3bb2fd58..0e4669e50 100644 --- a/tests/testthat/test-ts-hsj.R +++ b/tests/testthat/test-ts-hsj.R @@ -703,6 +703,128 @@ 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 +} + +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_present set-membership test + # (ts_hsj.cpp:220) alone, with fitch_label_char() (T-375) never entered. + b1 <- rbind(t1 = c("-", "1"), t2 = c("0", "0"), + t3 = c("1", "?"), t4 = c("?", "-")) + d1 <- make_hsj_dat(b1) + ref1 <- PhyDatToMatrix(d1)[tips, , drop = FALSE] + h1 <- CharacterHierarchy(`2` = integer(0)) + nTok1 <- length(attr(d1, "allLevels")) + perms1 <- as.matrix(expand.grid(rep(list(seq_len(nTok1)), nTok1))) + perms1 <- perms1[apply(perms1, 1, function(r) !anyDuplicated(r)), , drop = FALSE] + scores1 <- apply(perms1, 1, function(perm) { + d <- .PermuteTokens(d1, perm) + stopifnot(identical(PhyDatToMatrix(d)[tips, , drop = FALSE], ref1)) + hsj_score(tree, d, h1, alpha = 1) + }) + expect_equal(scores1, rep(scores1[[1]], length(scores1)), + info = "zero secondaries: isolates primary_present (ts_hsj.cpp:220)") + + # (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", "?")) + d2 <- make_hsj_dat(b2) + ref2 <- PhyDatToMatrix(d2)[tips, , drop = FALSE] + h2 <- CharacterHierarchy(`2` = 3L) + nTok2 <- length(attr(d2, "allLevels")) + perms2 <- as.matrix(expand.grid(rep(list(seq_len(nTok2)), nTok2))) + perms2 <- perms2[apply(perms2, 1, function(r) !anyDuplicated(r)), , drop = FALSE] + scores2 <- apply(perms2, 1, function(perm) { + d <- .PermuteTokens(d2, perm) + stopifnot(identical(PhyDatToMatrix(d)[tips, , drop = FALSE], ref2)) + hsj_score(tree, d, h2, alpha = 1) + }) + 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: HSJ + sectorial search (T-303 guard) # ========================================================================= From 317a9d562971bd5ec209b578ea740b78eb621553 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 12:41:20 +0100 Subject: [PATCH 03/11] docs: regenerate MaximizeParsimony.Rd (roxygen sync, pre-existing drift) Unrelated to the HSJ fix: R/MaximizeParsimony.R's @examples block was updated by an earlier commit (18397db2, "regenerate a MaximizeParsimony example left stale by the effort rename") but man/MaximizeParsimony.Rd wasn't regenerated to match. Caught by running devtools::check_man() before dispatch, as this branch's fix required no doc changes of its own. --- man/MaximizeParsimony.Rd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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{ From a29ff1d9be239cd517b78d37d6af7db1ee0c4e96 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 12:42:57 +0100 Subject: [PATCH 04/11] docs(hsj-paper-oracle): caveat that check [3] passing here isn't a T-374 fix Verified while validating the T-375/T-376 fix: check [3] now runs live (the alpha term is no longer inert) and passes on this small, Fig.1-derived matrix -- but a fresh sweep over larger/messier matrices shows T-374's rooting-dependence still reproduces readily post-fix, exactly as expected since T-374 needs the paper's two-state DP, which this fix does not touch. Add a comment so a future session doesn't read an all-green run here as T-374 being closed. --- dev/red-team/heavy-tests/hsj-paper-oracle.R | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/dev/red-team/heavy-tests/hsj-paper-oracle.R b/dev/red-team/heavy-tests/hsj-paper-oracle.R index 596f6e4c1..f08eaf2ca 100644 --- a/dev/red-team/heavy-tests/hsj-paper-oracle.R +++ b/dev/red-team/heavy-tests/hsj-paper-oracle.R @@ -154,6 +154,15 @@ 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 { + # A PASS on this small, Fig.1-derived matrix is NOT evidence that T-374 is + # fixed -- it only means this particular matrix/topology doesn't happen to + # trigger fitch_label_char()'s rooting-sensitivity. T-374's rooting bug is + # data/topology-dependent (measured 165-197/300 random cases) and remains + # OPEN and readily reproducible on other matrices after the T-375/T-376 fix + # (verified: a 10-tip matrix with more secondary ambiguity gives distinct + # scores across rootings at alpha=1, same mechanism as before). T-374's fix + # is the paper's two-state DP (or marginal-MPR resolution), not anything in + # this commit -- do not close T-374 on the strength of this check alone. for (alpha in c(0, 0.5, 1)) { scores <- vapply(rownames(extended), function(tip) { TreeLength(RootTree(base, tip), extDs, hierarchy = hierarchy, From 4c96ec36a50270e4f6e69eef9063b2a9a1d6b140 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 14:35:51 +0100 Subject: [PATCH 05/11] test(hsj): add a value-level T-375 regression test (secondary "?" resolution invariant) The permutation and misalignment tests added earlier both use hsj_alpha=0, which never calls count_mismatches() -- they exercise T-376's primary_present term exclusively and can't see whether fitch_label_char() resolves an ambiguous secondary correctly. This isolates T-375 with all primaries "1" (so the primary DP is loss-free and the whole score reduces to the secondary's Fitch step count) and checks T-375's own stated acceptance criterion: score("?") <= min over concrete resolutions. Verified to fail against a throwaway pre-fix build (scored 1, violating the invariant against min(0,1,1)=0) and pass post-fix. --- tests/testthat/test-ts-hsj.R | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/tests/testthat/test-ts-hsj.R b/tests/testthat/test-ts-hsj.R index 0e4669e50..eae29778a 100644 --- a/tests/testthat/test-ts-hsj.R +++ b/tests/testthat/test-ts-hsj.R @@ -824,6 +824,41 @@ test_that("HSJ handles a MatrixToPhyDat token/state misalignment (T-376)", { } }) +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) From 33702c9d745a21b9a940ba81e343fa728508af48 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 14:36:03 +0100 Subject: [PATCH 06/11] fix(ts_hsj_score): validate tip_labels/absent_state bounds Sharpened by the T-375/T-376 fix: a bad value used to be a harmless wrong scalar comparison, but now indexes DataSet::token_states directly (out-of-bounds read) or shifts by absent_state (undefined for >= 32). Only the internal ts_hsj_score() test bridge is reachable with hand-crafted values -- public callers always derive tip_labels/absent_state from a validated phyDat -- so this mirrors the existing validate_tip_data_values guard rather than changing any public behaviour. --- src/ts_rcpp.cpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 914314b1d..f4063d590 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -54,6 +54,28 @@ void validate_tip_data_values(const int* tip_data_r, int n_tips, } } +// Validate `tip_labels`/`absent_state` VALUES at the Rcpp boundary, for the +// `ts_hsj_score()` test bridge only: 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 (`1u << absent_state`), 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. +void validate_hsj_tip_labels(const IntegerMatrix& tip_labels_r, + int absent_state, int n_tokens) { + if (absent_state < 0 || absent_state >= 32) { + Rcpp::stop("`absent_state` must be in [0, 31]; found %d", 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 @@ -3122,6 +3144,8 @@ double ts_hsj_score( IntegerMatrix tip_labels_r, int absent_state) { + validate_hsj_tip_labels(tip_labels_r, absent_state, contrast.nrow()); + // Build DataSet for non-hierarchy characters (weight already adjusted) ts::DataSet ds = make_dataset(contrast, tip_data, weight, levels); From 71319a809d22d370291d86576e3eb7d44b35ef52 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 14:36:16 +0100 Subject: [PATCH 07/11] docs: promote T-374 caveat to printed output; log HSJ scoring change in NEWS hsj-paper-oracle.R's earlier caveat was a source comment on the branch that never executes when check [3] passes -- a future session reading stdout would miss it. Moved into a cat() so it prints on every passing run. NEWS.md: HSJ scores under inapplicable = "hsj" can change from previous versions now that the controlling primary's presence/absence and an ambiguous secondary are read correctly -- log it next to the existing x-transformation rooting entry, same "To integrate into 2.0.0" convention. --- NEWS.md | 20 ++++++++++++++++++++ dev/red-team/heavy-tests/hsj-paper-oracle.R | 16 +++++++--------- 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5b14f97a9..8c7c0d2c3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -28,6 +28,26 @@ 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** (generally + increasing, since a genuinely present or absent controlling primary is now + always counted). This is a correctness fix; the method's design -- + including its intentional rooting-sensitivity, not addressed here -- is + unchanged. + - `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/dev/red-team/heavy-tests/hsj-paper-oracle.R b/dev/red-team/heavy-tests/hsj-paper-oracle.R index f08eaf2ca..32e65a532 100644 --- a/dev/red-team/heavy-tests/hsj-paper-oracle.R +++ b/dev/red-team/heavy-tests/hsj-paper-oracle.R @@ -154,15 +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 { - # A PASS on this small, Fig.1-derived matrix is NOT evidence that T-374 is - # fixed -- it only means this particular matrix/topology doesn't happen to - # trigger fitch_label_char()'s rooting-sensitivity. T-374's rooting bug is - # data/topology-dependent (measured 165-197/300 random cases) and remains - # OPEN and readily reproducible on other matrices after the T-375/T-376 fix - # (verified: a 10-tip matrix with more secondary ambiguity gives distinct - # scores across rootings at alpha=1, same mechanism as before). T-374's fix - # is the paper's two-state DP (or marginal-MPR resolution), not anything in - # this commit -- do not close T-374 on the strength of this check alone. + 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, From d2a2ebca2f92447b53374b986d0851f6af930c10 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 18:21:30 +0100 Subject: [PATCH 08/11] fix(hsj): validate tip_labels/absent_state on the production search path too /code-review (xhigh) flagged this independently across 6 of 10 review angles: validate_hsj_tip_labels() only guarded the ts_hsj_score() test bridge, leaving unpack_hsj() -- which populates the same fields for the real ts_driven_search()/ts_collapse_pool() search path -- completely unvalidated. This is concretely reachable, not just hypothetical: R/ts-driven-compat.R's legacy flat-argument wrapper exposes hsjTipLabels/ hsjAbsentState as independent, hand-settable arguments feeding straight into unpack_hsj(). Since this PR removed fitch_label_char()'s old defensive clamp in favour of direct token_states[label] indexing, a bad value here is now an out-of-bounds heap read instead of a harmless wrong comparison. Also tightens the absent_state bound from a hardcoded [0,31] (shift safety only) to [0, n_levels) -- the dataset's actual state count -- since a value in between passed silently but reproduced the exact T-375/T-376 symptom class (the primary is never classified absent). Adds a default initializer to HierarchyBlock::absent_state for defensive symmetry with DataSet::inapp_state's existing default. Verified: a hand-crafted out-of-range hsjTipLabels/hsjAbsentState now raises a clean Rcpp::stop() from all three entry points (ts_hsj_score, and now ts_driven_search via the compat wrapper), where it previously either mis-scored silently or would have read out of bounds. --- src/ts_data.h | 2 +- src/ts_rcpp.cpp | 27 ++++++++++++++++++--------- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/src/ts_data.h b/src/ts_data.h index 53390c1d7..5c078fda2 100644 --- a/src/ts_data.h +++ b/src/ts_data.h @@ -73,7 +73,7 @@ struct HierarchyBlock { // 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; + int absent_state = -1; }; struct CharBlock { diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index f4063d590..fe0a17639 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -54,19 +54,24 @@ void validate_tip_data_values(const int* tip_data_r, int n_tips, } } -// Validate `tip_labels`/`absent_state` VALUES at the Rcpp boundary, for the -// `ts_hsj_score()` test bridge only: 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 (`1u << absent_state`), with no further +// 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. +// 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) { - if (absent_state < 0 || absent_state >= 32) { - Rcpp::stop("`absent_state` must be in [0, 31]; found %d", absent_state); + 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) { @@ -1844,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; @@ -3144,7 +3152,8 @@ double ts_hsj_score( IntegerMatrix tip_labels_r, int absent_state) { - validate_hsj_tip_labels(tip_labels_r, absent_state, contrast.nrow()); + 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); From bd2c74b85932f4d1b4a970df828fc6b9bb3894c3 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 18:21:47 +0100 Subject: [PATCH 09/11] perf(hsj): restore per-character tie-break array sizing, reuse ctz64() /code-review flagged this independently via two angles: sizing fitch_label_char()'s tb_cnt/tb_mintip arrays by the dataset-global n_levels (instead of the old per-character state count) meant a dataset with a high global state count but binary hierarchy secondaries paid up to ~16x more allocation/scan cost per secondary character per HSJ full-rescore call than before this PR. The downpass only ever intersects/unions existing tip bits, so no bit outside the union of tip state-sets can appear at any node -- accumulate that union while initializing state_sets (one extra OR per tip, already being read) and size K from its highest set bit instead of the global n_levels. Restores the old tightness with no extra tree traversal. Also replaces the tie-break credit loop's manual "scan for the single set bit" with ctz64() (already declared in ts_data.h, already used elsewhere in this codebase for the identical operation), removing a 7-line loop-with-break in favour of one call. No test changes needed: this is a pure sizing/lookup optimisation, not a behaviour change -- full test-ts-hsj.R (85 assertions) and the hsj-token-permutation.R/hsj-paper-oracle.R heavy tests still pass unchanged. --- src/ts_hsj.cpp | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/src/ts_hsj.cpp b/src/ts_hsj.cpp index 5f0892256..5d70e7f30 100644 --- a/src/ts_hsj.cpp +++ b/src/ts_hsj.cpp @@ -57,9 +57,11 @@ static int fitch_label_char( int n_tip = tree.n_tip; int n_node = tree.n_node; + uint32_t used_mask = 0; for (int t = 0; t < n_tip; ++t) { int label = tip_labels[t * n_orig_chars + char_idx]; state_sets[t] = token_states[label]; + used_mask |= state_sets[t]; } // --- Downpass --- @@ -95,7 +97,18 @@ static int fitch_label_char( // // 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 - const int K = n_levels; + // + // 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); @@ -105,13 +118,9 @@ static int fitch_label_char( // 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) { - for (int s = 0; s < K; ++s) { - if (set & (1u << s)) { - tb_cnt[static_cast(t) * K + s] = 1; - tb_mintip[static_cast(t) * K + s] = t; - break; - } - } + 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) { From a76338f8200bfd6d70cae7ba746c83da77042ffe Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 18:22:08 +0100 Subject: [PATCH 10/11] docs(hsj): fix review-flagged doc gaps -- NEWS accuracy, stale terminology, vignette /code-review (xhigh) findings addressed: - NEWS.md's new HSJ entry called T-374's rooting-dependence "intentional" (contradicts this repo's own settled finding that it's a tracked bug per the Hopkins & St John paper) and claimed scores are "generally increasing" (contradicted by this PR's own T-375 test, which shows a decrease case). Corrected both. - .BuildTipLabels()'s own docstring (R/CharacterHierarchy.R) still called its output "state labels" -- the exact token/state conflation this PR fixed everywhere else except at the one place that actually produces the token-space data. - A pre-existing test helper comment (test-ts-hsj.R) called absent_state a "token index"; it is a state index, per this PR's own terminology fix elsewhere. - A new test's comment + expect_equal info= string cited a hardcoded ts_hsj.cpp line number that had already drifted from the code it described (and would surface the wrong line in CI failure output). - ts_sector.cpp's two T-303 guard-rationale comments enumerated the DataSet fields build_reduced_dataset() doesn't copy, but didn't mention the two fields (token_states, n_levels) this PR added. - vignettes/search-algorithm.Rmd wasn't updated despite AGENTS.md's mandatory-checks rule for any scoring-behaviour change (this PR's own NEWS.md entry documents exactly that); added a short note that HSJ scoring is correct regardless of a dataset's internal token ordering. Also refactors the new contrast-row-permutation test's two near-identical permutation-generation blocks into one shared .AllTokenOrderingScores() helper, called twice -- flagged as duplicated logic that would otherwise need fixing in two places. --- NEWS.md | 11 +++---- R/CharacterHierarchy.R | 4 ++- src/ts_sector.cpp | 10 ++++--- tests/testthat/test-ts-hsj.R | 53 +++++++++++++++++----------------- vignettes/search-algorithm.Rmd | 6 ++++ 5 files changed, 47 insertions(+), 37 deletions(-) diff --git a/NEWS.md b/NEWS.md index 8c7c0d2c3..1f91a204a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -42,11 +42,12 @@ 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** (generally - increasing, since a genuinely present or absent controlling primary is now - always counted). This is a correctness fix; the method's design -- - including its intentional rooting-sensitivity, not addressed here -- is - unchanged. + **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 HSJ's + known rooting-sensitivity, 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 diff --git a/R/CharacterHierarchy.R b/R/CharacterHierarchy.R index 385acf3e2..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") 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 eae29778a..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, @@ -733,44 +733,43 @@ test_that("HSJ secondary dissimilarity is level-order invariant (multistate)", { 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_present set-membership test - # (ts_hsj.cpp:220) alone, with fitch_label_char() (T-375) never entered. + # (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("?", "-")) - d1 <- make_hsj_dat(b1) - ref1 <- PhyDatToMatrix(d1)[tips, , drop = FALSE] - h1 <- CharacterHierarchy(`2` = integer(0)) - nTok1 <- length(attr(d1, "allLevels")) - perms1 <- as.matrix(expand.grid(rep(list(seq_len(nTok1)), nTok1))) - perms1 <- perms1[apply(perms1, 1, function(r) !anyDuplicated(r)), , drop = FALSE] - scores1 <- apply(perms1, 1, function(perm) { - d <- .PermuteTokens(d1, perm) - stopifnot(identical(PhyDatToMatrix(d)[tips, , drop = FALSE], ref1)) - hsj_score(tree, d, h1, alpha = 1) - }) + 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_present (ts_hsj.cpp:220)") + 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", "?")) - d2 <- make_hsj_dat(b2) - ref2 <- PhyDatToMatrix(d2)[tips, , drop = FALSE] - h2 <- CharacterHierarchy(`2` = 3L) - nTok2 <- length(attr(d2, "allLevels")) - perms2 <- as.matrix(expand.grid(rep(list(seq_len(nTok2)), nTok2))) - perms2 <- perms2[apply(perms2, 1, function(r) !anyDuplicated(r)), , drop = FALSE] - scores2 <- apply(perms2, 1, function(perm) { - d <- .PermuteTokens(d2, perm) - stopifnot(identical(PhyDatToMatrix(d)[tips, , drop = FALSE], ref2)) - hsj_score(tree, d, h2, alpha = 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)") }) 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 From 465be292311d323acc06116b9d536fac0ff289c5 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Sat, 1 Aug 2026 18:39:35 +0100 Subject: [PATCH 11/11] fix(news): reword to avoid an unlisted possessive-acronym spelling flag GHA's R-CMD-check caught what my local spell-check missed: "HSJ's" isn't in inst/WORDLIST (only bare "HSJ" is -- spelling::spell_check_test() doesn't stem across the apostrophe, same reason "TNT's" needed its own separate WORDLIST entry elsewhere in this file). Reworded rather than added a wordlist entry, per AGENTS.md's stated preference. My own verification gap: I'd been running spelling::spell_check_package(vignettes = FALSE) all session, which apparently doesn't scan NEWS.md the same way tests/spelling.R's actual spell_check_test(vignettes = TRUE) does. Verified against that exact invocation this time; both spell_check_package() (default args) and spell_check_test(vignettes = TRUE) now pass clean. --- NEWS.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 1f91a204a..4a8994670 100644 --- a/NEWS.md +++ b/NEWS.md @@ -46,8 +46,9 @@ 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 HSJ's - known rooting-sensitivity, which remains a separate, open issue. + 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