fix(hsj): translate token indices to state indices in HSJ scoring (T-375, T-376)

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

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.

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

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,

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(

src/ts_data.h

@@ -68,7 +68,12 @@ struct HierarchyBlock {
   int primary_char;               // original character index (0-based)
   std::vector<int> 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<uint32_t> 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<HierarchyBlock> 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<int> tip_labels;
   int n_orig_chars = 0;

src/ts_hsj.cpp

@@ -32,29 +32,36 @@ std::vector<int> 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<int>& tip_labels,
     int char_idx,
     int n_orig_chars,
-    int inapp_state,
+    const std::vector<uint32_t>& token_states,
+    int n_levels,
     std::vector<uint32_t>& 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<int>::max();
   std::vector<int> tb_cnt(static_cast<size_t>(n_node) * K, 0);
   std::vector<int> tb_mintip(static_cast<size_t>(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<size_t>(t) * K + label] = 1;
-      tb_mintip[static_cast<size_t>(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<size_t>(t) * K + s] = 1;
+      tb_mintip[static_cast<size_t>(t) * K + s] = t;
     }
   }
   for (int i = 0; i < static_cast<int>(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<int>& tip_labels,
     int n_orig_chars,
-    int inapp_state)
+    int inapp_state,
+    const std::vector<uint32_t>& 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<uint32_t> 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<bool> 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;

src/ts_hsj.h

@@ -51,7 +51,10 @@ std::vector<int> 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)