diff --git a/NEWS.md b/NEWS.md index 3492cf2d4..e104fe8c3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -88,6 +88,15 @@ known rooting-sensitivity of HSJ scoring, which remains a separate, open issue. +- Zero-length-branch collapse (`collapse = TRUE`) no longer disables itself + for an `inapplicable = "hsj"`/`"xform"` search whenever *no* hierarchy + block actually exists in that replicate -- previously it keyed on the + scoring mode alone. This only affects `Resample()`, whose bootstrap and + jackknife replicates can drop every hierarchy block from a unit while + still passing a (now-empty) hierarchy config through; those replicates are + ordinary Fitch data and now collapse like any other. A replicate that + retains any hierarchy block is unaffected. + - `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/src/ts_collapsed.cpp b/src/ts_collapsed.cpp index 7d41e98d0..ce33dfe61 100644 --- a/src/ts_collapsed.cpp +++ b/src/ts_collapsed.cpp @@ -21,8 +21,14 @@ void compute_collapsed_flags( // (all-zero flags == nothing collapses — the safe conservative outcome). // Falling back to the conservative flags is NOT sufficient: it is equally // blind to hierarchy/Sankoff support. See red-team T-330. - if (ds.scoring_mode == ScoringMode::HSJ || - ds.scoring_mode == ScoringMode::XFORM) return; + // + // Gate on whether hierarchy data actually EXISTS (T-408), not merely on + // scoring_mode: an HSJ/XFORM config with no hierarchy_blocks / sankoff_n_chars + // carries no topology-dependent support this kernel is blind to, so collapse + // is safe. Matches DataSet::topology_independent()'s predicate. + if ((ds.scoring_mode == ScoringMode::HSJ || + ds.scoring_mode == ScoringMode::XFORM) && + (!ds.hierarchy_blocks.empty() || ds.sankoff_n_chars > 0)) return; // If all characters were simplified away (total_words == 0), every binary // resolution ties at the same score: no internal branch carries support, @@ -150,8 +156,12 @@ void compute_collapsed_flags_aggressive( // modes (all-zero flags). Guarded independently of compute_collapsed_flags: // the has_na delegation at the bottom of this block only reaches it on NA // data, not the general HSJ/XFORM case. See red-team T-330. - if (ds.scoring_mode == ScoringMode::HSJ || - ds.scoring_mode == ScoringMode::XFORM) { + // + // Gate on hierarchy data presence, not scoring_mode alone (T-408); see + // compute_collapsed_flags() above for the rationale. + if ((ds.scoring_mode == ScoringMode::HSJ || + ds.scoring_mode == ScoringMode::XFORM) && + (!ds.hierarchy_blocks.empty() || ds.sankoff_n_chars > 0)) { collapsed.assign(tree.n_node, 0); return; } diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index fcb6678c0..8ff590920 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -1850,14 +1850,44 @@ static void unpack_hsj(Nullable hsjConfig, ts::DataSet& ds) { ds.hsj_alpha = as(hc["hsjAlpha"]); ds.scoring_mode = ts::ScoringMode::HSJ; - if (hc.containsElementNamed("hsjTipLabels") && - !Rf_isNull(hc["hsjTipLabels"])) { + // hsjTipLabels must be present and non-NULL whenever HSJ is enabled: + // score_hierarchy_block() reads ds.tip_labels unconditionally once + // scoring_mode == HSJ, and that field is only populated inside this + // branch. `list(hsjTipLabels = NULL)` keeps the element name, so + // containsElementNamed() alone does not catch an empty tip_labels (T-398). + if (!hc.containsElementNamed("hsjTipLabels") || + Rf_isNull(hc["hsjTipLabels"])) { + Rcpp::stop("hsjConfig$hsjTipLabels must be provided (non-NULL) whenever " + "hsjConfig is supplied (which enables HSJ scoring)."); + } + + { 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(); + // hsjTipLabels must cover every block's primary/secondary character + // index: score_hierarchy_block() reads + // tip_labels[t * n_orig_chars + block.primary_char] (and likewise for + // secondaries) unconditionally once scoring_mode == HSJ, so a + // non-NULL but too-narrow matrix reads past ds.tip_labels the same + // way a NULL one did (T-398). + for (const ts::HierarchyBlock& block : ds.hierarchy_blocks) { + if (block.primary_char < 0 || block.primary_char >= n_c) { + Rcpp::stop("hsjConfig$hsjTipLabels has %d columns, but a hierarchy " + "block's primary character index is %d", + n_c, block.primary_char); + } + for (int sec : block.secondary_chars) { + if (sec < 0 || sec >= n_c) { + Rcpp::stop("hsjConfig$hsjTipLabels has %d columns, but a " + "hierarchy block's secondary character index is %d", + n_c, sec); + } + } + } ds.n_orig_chars = n_c; ds.tip_labels.resize(n_t * n_c); for (int t = 0; t < n_t; ++t) { @@ -1900,6 +1930,16 @@ static void unpack_xform(Nullable xformConfig, List rc = xf_list[ch]; NumericMatrix cm = as(rc["cost_matrix"]); int ns = ns_vec[ch]; + // Validate cost matrix dimensions match the character's state count + // (mirrors the check ts_sankoff_test() already performs; T-397 — + // Rcpp's Matrix indexing never bounds-checks a mis-shaped-but- + // same-length matrix, so an unguarded read here silently scores + // garbage instead of erroring). + if (cm.nrow() != ns || cm.ncol() != ns) { + Rcpp::stop("xformChars[[%d]]$cost_matrix has dimensions %d x %d, but " + "character %d has %d states (expected %d x %d)", + ch + 1, cm.nrow(), cm.ncol(), ch + 1, ns, ns, ns); + } double* dst = ds.sankoff_cost_matrices.data() + static_cast(ch) * max_ns * max_ns; for (int r = 0; r < ns; ++r) @@ -1926,6 +1966,22 @@ static void unpack_xform(Nullable xformConfig, IntegerMatrix combo_grid = as(rc["combo_grid"]); IntegerMatrix tip_sec = as(rc["tip_sec_known"]); int n_sec = combo_grid.ncol(); + // combo_grid must carry one row per present state (states 1..ns-1); + // state == -2 below indexes it at (s - 1) for s up to ns - 1, so + // fewer rows than that reads out of bounds (T-397). + if (combo_grid.nrow() != ns - 1) { + Rcpp::stop("xformChars[[%d]]$combo_grid has %d rows, but character " + "%d has %d states (expected %d rows)", + ch + 1, combo_grid.nrow(), ch + 1, ns, ns - 1); + } + // tip_sec_known is read at (t, d) for t in [0, n_t), d in + // [0, n_sec) in the state == -2 branch below; a truncated matrix + // reads past the SEXP the same way an unguarded combo_grid would. + if (tip_sec.nrow() != n_t || tip_sec.ncol() != n_sec) { + Rcpp::stop("xformChars[[%d]]$tip_sec_known has dimensions %d x %d, " + "but expected %d x %d (n_tips x n_secondaries)", + ch + 1, tip_sec.nrow(), tip_sec.ncol(), n_t, n_sec); + } for (int t = 0; t < n_t; ++t) { int state = ts_r[t]; double* tip_ptr = ds.sankoff_tip_costs.data() + @@ -1950,6 +2006,15 @@ static void unpack_xform(Nullable xformConfig, } } else if (state >= 0 && state < ns) { tip_ptr[state] = 0.0; + } else { + // Any other value (e.g. state >= ns) falls through every branch + // above, leaving tip_ptr all-INF; that INF then propagates through + // the pool sentinel (1e18) rather than a true Inf and passes + // is.finite(), silently corrupting the score instead of erroring + // (T-397). + Rcpp::stop("xformChars[[%d]]$tip_states[%d] = %d is out of range; " + "must be -1, -2, or in [0, %d)", + ch + 1, t + 1, state, ns); } } } diff --git a/tests/testthat/test-ts-hsj-xform-guards.R b/tests/testthat/test-ts-hsj-xform-guards.R new file mode 100644 index 000000000..9bb3e803a --- /dev/null +++ b/tests/testthat/test-ts-hsj-xform-guards.R @@ -0,0 +1,328 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +# Regression tests for three guard clauses at the HSJ/XFORM bridge +# (src/ts_rcpp.cpp: unpack_hsj(), unpack_xform()) and the collapse gate +# (src/ts_collapsed.cpp). +# +# T-398 (#14): unpack_hsj() enabled HSJ scoring (scoring_mode = HSJ) even when +# hsjTipLabels was present-but-NULL, leaving ds.tip_labels empty and +# segfaulting score_hierarchy_block(). Fixed with an explicit presence/NULL +# check that Rcpp::stop()s instead of silently skipping population. Review +# also surfaced a second entrance to the same crash class -- a non-NULL but +# too-narrow hsjTipLabels not covering every block's primary/secondary index +# -- closed with the same guard. +# +# T-397 (#13): unpack_xform() never validated cost_matrix dimensions, +# combo_grid row count, tip_sec_known dimensions, or tip_states range +# against n_states, unlike its sibling ts_sankoff_test() (guarded since +# 0856748f). A same-length-but-wrong-shape cost matrix reads garbage with no +# warning at all (Rcpp's Matrix::operator() only bounds-checks the linear +# offset, not (row, col)). +# +# T-408 (#22): the collapse guards in ts_collapsed.cpp keyed on +# ds.scoring_mode alone, disabling collapse for an HSJ/XFORM config with NO +# hierarchy data -- a case collapse is provably safe for. Fixed to gate on +# hierarchy-data presence (matching DataSet::topology_independent()). + +library("TreeTools") + +.HierarchyToBlocks <- TreeSearch:::.HierarchyToBlocks +.BuildTipLabels <- TreeSearch:::.BuildTipLabels +.HSJAbsentState <- TreeSearch:::.HSJAbsentState +.NonHierarchyWeights <- TreeSearch:::.NonHierarchyWeights +ts_driven_search <- TreeSearch:::ts_driven_search +ts_collapse_pool <- TreeSearch:::ts_collapse_pool + +make_dat <- function(mat, levels = c("-", "0", "1")) { + phangorn::phyDat(mat, type = "USER", levels = levels, ambiguity = "?") +} + + +# ========================================================================= +# T-398 / #14: hsjTipLabels omitted (left at compat-wrapper default NULL) +# must error cleanly, not segfault. +# +# A segfault kills the test process outright, so it cannot be asserted from +# inside testthat. The pre-fix segfault was reproduced separately with a +# standalone Rscript (exit code 139) -- see the PR body for that output. +# ========================================================================= +test_that("HSJ scoring with omitted hsjTipLabels errors instead of segfaulting", { + mat <- matrix(c( + "0", "0", + "0", "0", + "1", "0", + "1", "1", + "1", "1", + "1", "-" + ), nrow = 6, byrow = TRUE, + dimnames = list(paste0("t", 1:6), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("2" = integer(0)) + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + blocks <- .HierarchyToBlocks(h) + + # hsjTipLabels intentionally omitted -- exercises the compat wrapper's own + # default (R/ts-driven-compat.R), which previously reached unpack_hsj() + # with hsjConfig$hsjTipLabels == NULL. + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + hierarchyBlocks = blocks, + hsjAlpha = 1.0, + hsjAbsentState = .HSJAbsentState(ds), + maxReplicates = 1L + ), + "hsjTipLabels" + ) +}) + +test_that("HSJ scoring with too-narrow hsjTipLabels errors", { + mat <- matrix(c( + "0", "0", + "0", "0", + "1", "0", + "1", "1", + "1", "1", + "1", "-" + ), nrow = 6, byrow = TRUE, + dimnames = list(paste0("t", 1:6), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("2" = integer(0)) + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + blocks <- .HierarchyToBlocks(h) + tl <- .BuildTipLabels(ds) + + # A non-NULL but too-narrow hsjTipLabels (missing the primary's column) + # must be caught the same way an omitted one is: score_hierarchy_block() + # indexes tip_labels at [t * n_orig_chars + block$primary], so dropping + # the column covering that index reads past the vector. + narrow_tl <- tl[, -ncol(tl), drop = FALSE] + + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + hierarchyBlocks = blocks, + hsjTipLabels = narrow_tl, + hsjAlpha = 1.0, + hsjAbsentState = .HSJAbsentState(ds), + maxReplicates = 1L + ), + "hsjTipLabels" + ) +}) + + +# ========================================================================= +# T-397 / #13: mis-shaped cost matrix through unpack_xform() must error +# with the expected/actual dimensions, matching ts_sankoff_test()'s style. +# Covers the same-length-but-wrong-shape case (1x9 for a 3x3) that Rcpp's +# own indexing does not warn about. +# ========================================================================= +test_that("Xform bridge errors on mis-shaped cost matrix", { + # Two informative secondary states (sec = 0, 1) -> n_states = 1(absent) + + # 2(present combos) = 3, giving a 3x3 cost matrix. + mat <- matrix(c( + "0", "-", + "1", "0", + "1", "1", + "1", "0" + ), nrow = 4, byrow = TRUE, + dimnames = list(paste0("t", 1:4), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("1" = 2L) + + recoded <- RecodeHierarchy(ds, h) + blk <- recoded$sankoff_chars[[1]] + expect_equal(dim(blk$cost_matrix), c(3, 3)) + + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + + # Same-length-but-wrong-shape: 1x9 carries the same 9 values as the 3x3 + # matrix, so Rcpp's linear-offset indexing reads it with no warning. + bad_blk <- blk + bad_blk$cost_matrix <- matrix(as.vector(blk$cost_matrix), nrow = 1, ncol = 9) + + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + xformChars = list(bad_blk), + maxReplicates = 1L + ), + "cost_matrix has dimensions 1 x 9.*3 states" + ) +}) + +test_that("Xform bridge errors on undersized combo_grid", { + mat <- matrix(c( + "1", "0", "0", + "1", "0", "0", + "1", "1", "1", + "1", "1", "1" + ), nrow = 4, byrow = TRUE, + dimnames = list(paste0("t", 1:4), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("1" = 2:3) + + recoded <- RecodeHierarchy(ds, h) + blk <- recoded$sankoff_chars[[1]] + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + + bad_blk <- blk + bad_blk$combo_grid <- blk$combo_grid[-1, , drop = FALSE] + + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + xformChars = list(bad_blk), + maxReplicates = 1L + ), + "combo_grid has" + ) +}) + +test_that("Xform bridge errors on mis-shaped tip_sec_known", { + mat <- matrix(c( + "1", "0", "0", + "1", "0", "0", + "1", "1", "1", + "1", "1", "1" + ), nrow = 4, byrow = TRUE, + dimnames = list(paste0("t", 1:4), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("1" = 2:3) + + recoded <- RecodeHierarchy(ds, h) + blk <- recoded$sankoff_chars[[1]] + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + + bad_blk <- blk + bad_blk$tip_sec_known <- blk$tip_sec_known[-1, , drop = FALSE] + + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + xformChars = list(bad_blk), + maxReplicates = 1L + ), + "tip_sec_known has" + ) +}) + +test_that("Xform bridge errors on out-of-range tip_states", { + mat <- matrix(c( + "0", "-", "0", + "1", "0", "1", + "1", "0", "0", + "1", "0", "1" + ), nrow = 4, byrow = TRUE, + dimnames = list(paste0("t", 1:4), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("1" = 2L) + + recoded <- RecodeHierarchy(ds, h) + blk <- recoded$sankoff_chars[[1]] + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + + bad_blk <- blk + bad_blk$tip_states[1] <- blk$n_states + 5L # out of [0, n_states) + + expect_error( + ts_driven_search( + contrast = at$contrast, + tip_data = tip_data, + weight = adj_w, + levels = at$levels, + xformChars = list(bad_blk), + maxReplicates = 1L + ), + "out of range" + ) +}) + + +# ========================================================================= +# T-408 / #22: an HSJ config with empty hierarchyBlocks must collapse +# zero-length branches exactly as the no-hsjConfig (EW) case does -- the +# guard must not disable collapse when no hierarchy data is actually +# present, only when it is. +# ========================================================================= +test_that("Collapse fires for an HSJ config with no hierarchy blocks", { + # Reuses the T-330 reproducing configuration (char1 zero-length-supports + # the (t4,t5) clade once char2's weight is zeroed) but with an HSJ config + # whose hierarchy_blocks is empty -- scoring_mode == HSJ, yet no hierarchy + # data exists for the collapse kernel to be blind to. + mat <- matrix(c( + "0", "0", + "0", "0", + "1", "0", + "1", "1", + "1", "1", + "1", "-" + ), nrow = 6, byrow = TRUE, + dimnames = list(paste0("t", 1:6), NULL)) + ds <- make_dat(mat) + h <- CharacterHierarchy("2" = integer(0)) + tr <- Preorder(RenumberTips( + ape::read.tree(text = "(((t1,t2),(t3,(t4,t5))),t6);"), names(ds))) + at <- attributes(ds) + adj_w <- as.integer(.NonHierarchyWeights(ds, h)) + tip_data <- matrix(unlist(ds, use.names = FALSE), + nrow = length(ds), byrow = TRUE) + scoringConfig <- list(min_steps = integer(0), concavity = Inf, + xpiwe = FALSE, xpiwe_r = 0.5, xpiwe_max_f = 5.0, + obs_count = integer(0), infoAmounts = NULL) + + hsjConfig <- list( + hierarchyBlocks = list(), + hsjAlpha = 1.0, + hsjTipLabels = matrix(integer(0), nrow = length(ds), ncol = 0), + hsjAbsentState = 0L) + + n_in <- nrow(tr$edge) + cp_ew <- ts_collapse_pool( + list(tr$edge), at$contrast, tip_data, adj_w, at$levels, + scoringConfig, NULL, NULL, NULL) + cp_hsj <- ts_collapse_pool( + list(tr$edge), at$contrast, tip_data, adj_w, at$levels, + scoringConfig, hsjConfig, NULL, NULL) + + # Sanity: the (t4,t5) clade genuinely is collapsible under plain EW. + expect_lt(nrow(cp_ew$trees[[1]]), n_in) + # The empty-hierarchy HSJ config must collapse identically -- not be + # blocked by the scoring_mode-only guard. + expect_equal(nrow(cp_hsj$trees[[1]]), nrow(cp_ew$trees[[1]])) +})