From 7c8c3ab0676df9d00570a2913a22ed06cb6e5cf9 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:31:23 +0100 Subject: [PATCH] fix: reject non-binary trees at the TreeState boundary TreeState::init_from_edge derives n_tip, n_internal and n_node from the edge count alone, which identifies a tree only when it is binary. A multifurcating edge list broke that derivation in both parities: an odd edge count left the topology arrays one element short, so parent[] and left[]/right[] were written out of bounds, while an even one kept the indices in range but rooted the tree on a real tip, leaving a one-element postorder whose downpass read prelim.data() - total_words. Either way the caller got a plausible number instead of an error, and the number varied between identical calls; a polytomous startEdge segfaulted. Refuse the edge list at that boundary rather than at each of the R entry points that reach it: edge_list_is_binary() checks the shape from the edge arrays alone, before anything is written, and init_from_edge throws for Rcpp to forward. build_topology_tree() clones the same derivation for the least-squares path and gets the same check. ts_driven_search() screens start trees with the predicate on the main thread, since an uncaught throw on a parallel worker would terminate the session. TreeLength.list(), .CheckTreeCharLen(), TreeScore() and EdgeListScore() gain the R-level check so the message matches the one TreeLength.phylo() has always given. FastCharacterLength() is left unchecked, as documented; the kernel now gives it the same message. The Shiny app scores every tree it displays, so it now searches with collapse = FALSE. Also extends the T-261 zeroing audit in reset_states() to name the collapse kernels, whose whole-row memcmps read words no pass writes, and to name every path that zero-fills the state arrays rather than only init_from_edge. The T-382 one-sidedness comment keeps its original reasoning for prelim -- a tip sibling always carries real states, so a ratchet-zeroed block makes equality harder -- and gains the down2 / subtree_actives case, where the words really are always zero. Fixes #16 Fixes #24 Co-Authored-By: Claude Opus 5 --- NEWS.md | 19 ++ R/tree_length.R | 18 ++ inst/Parsimony/server/mod_search.R | 6 + src/ts_collapsed.cpp | 20 +- src/ts_rcpp.cpp | 16 ++ src/ts_tree.cpp | 49 +++++ src/ts_tree.h | 12 ++ .../test-ts-t400-multifurcating-guard.R | 179 ++++++++++++++++++ 8 files changed, 312 insertions(+), 7 deletions(-) create mode 100644 tests/testthat/test-ts-t400-multifurcating-guard.R diff --git a/NEWS.md b/NEWS.md index 3492cf2d4..0deb3057e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,24 @@ # To integrate into 2.0.0 notes +- `TreeLength()`, `CharacterLength()`, `TreeScore()` and `EdgeListScore()` -- and + so `Consistency()`, `ExpectedLength()`, `ConcordantInformation()`, + `LengthAdded()` and `SuccessiveApproximations()`, which score trees through + them -- now reject a + tree that contains a polytomy, with the "`tree` must be binary" error that + `TreeLength()` already gave for a single `phylo` tree. Such a tree + previously returned a number. The scoring engine derives its node counts from + the number of edges, which identifies a tree only if that tree is binary: a + polytomous tree with an odd number of edges wrote past the end of the arrays + holding its topology, and one with an even number of edges was rooted on a + leaf and then scored from memory outside its own state buffer, so repeating + the same call could return a different answer each time. `MaximizeParsimony()` + collapses the trees it returns unless `collapse = FALSE`, so scoring its output + reached this path; search with `collapse = FALSE` to obtain trees that can be + scored, whose lengths are the score the search reports. Resolving a collapsed + tree instead, with `TreeTools::MakeTreeBinary()`, does not recover that score: + an arbitrary resolution of a polytomy need not be one of the most parsimonious + ones. + - `inapplicable = "xform"` scores are now reported at a canonical rooting, so a reported score is reproducible. The x-transformation's step matrix is asymmetric -- a gain costs one more than the number of secondary characters it diff --git a/R/tree_length.R b/R/tree_length.R index 54191f270..7d558d3ed 100644 --- a/R/tree_length.R +++ b/R/tree_length.R @@ -328,6 +328,9 @@ TreeLength.list <- function(tree, dataset, concavity = Inf, paste0(nEdge, collapse = ", "), "); try collapsing polytomies?)") } + if (nEdge != nTip + nTip - 2) { + stop("`tree` must be binary") + } if (is.null(attr(dataset, "levels")) || ncol(attr(dataset, "contrast")) == 0L) { return(rep(0L, length(tree))) @@ -448,6 +451,10 @@ Fitch <- function(tree, dataset) { if (!TreeIsRooted(tree)) { stop("`tree` must be rooted; try RootTree(tree)") } + nTip <- length(TipLabels(tree)) + if (dim(tree[["edge"]])[1] != nTip + nTip - 2) { + stop("`tree` must be binary") + } } #' @importFrom cli cli_alert @@ -576,6 +583,9 @@ TreeScore <- function(tree, dataset) { stop("Number of taxa in dataset (", nTaxa, ") not equal to number of tips in tree") } + if (dim(tree[["edge"]])[1] != nTaxa + nTaxa - 2) { + stop("`tree` must be binary") + } tree <- RenumberTips(tree, dataset[["tip.label"]]) el <- RenumberEdges(tree[["edge"]][, 1], tree[["edge"]][, 2]) # Return: @@ -598,6 +608,14 @@ EdgeListScore <- function(parent, child, dataset, inPostorder = FALSE, ...) { stop("`dataset` must be a `ParsimonyData` object; prepare it first with ", "`PrepareData()`, or supply your own `TreeScorer`.") } + # Every internal node of a rooted binary tree parents exactly two children; + # the scoring kernel derives its node counts from the edge count alone, so a + # polytomy makes it index out of bounds. This catches that case, to give the + # same message as the other entry points; the kernel checks the rest. + nChild <- tabulate(parent) + if (any(nChild != 0L & nChild != 2L)) { + stop("`tree` must be binary") + } if (!inPostorder) { edgeList <- Preorder(cbind(parent, child)) edgeList <- edgeList[PostorderOrder(edgeList), , drop = FALSE] diff --git a/inst/Parsimony/server/mod_search.R b/inst/Parsimony/server/mod_search.R index d091a7e47..c5dd1b8cd 100644 --- a/inst/Parsimony/server/mod_search.R +++ b/inst/Parsimony/server/mod_search.R @@ -627,6 +627,11 @@ search_server <- function(id, r, AnyTrees, HaveData, UpdateAllTrees, log_fns) { targetHits = targetHits, maxSeconds = maxSeconds, nThreads = nThreads, + # The app scores every returned tree, to display its length and to + # apply the suboptimality tolerance; only a binary tree can be + # scored, and an arbitrary resolution of a collapsed one need not + # be most parsimonious. + collapse = FALSE, verbosity = 0L ) # Only pass control when non-default, so the effort rung applies @@ -833,6 +838,7 @@ search_server <- function(id, r, AnyTrees, HaveData, UpdateAllTrees, log_fns) { if (identical(searchInapplicable, "hsj") && !is.null(searchHierarchy) && searchHsjAlpha != 1.0) paste0(" hsj_alpha = ", searchHsjAlpha, ","), + " collapse = FALSE,", " verbosity = 0", ")")) diff --git a/src/ts_collapsed.cpp b/src/ts_collapsed.cpp index 7d41e98d0..ec855aa84 100644 --- a/src/ts_collapsed.cpp +++ b/src/ts_collapsed.cpp @@ -113,19 +113,25 @@ void compute_collapsed_flags( } // --- Condition 3: prelim[sibling] == prelim[parent] --- - // This full-row memcmp also reads words belonging to ratchet-zeroed - // blocks (active_mask == 0), which fitch_downpass leaves stale rather - // than updating. That staleness only ever makes equality *harder* to - // reach (a stale word is unlikely to coincidentally match), so its only - // effect is to under-flag collapsible edges — a lost optimisation, never - // a false collapse. One-sided safe; not worth a masked per-word compare. - // See red-team T-382. + // This full-row memcmp also spans words fitch_downpass does not write: the + // SIMD pad word, and the words of ratchet-zeroed blocks (active_mask == 0), + // which it skips. The pad word reads zero on both sides and so contributes + // nothing. A zeroed block's words are whatever the node last held, while a + // tip sibling always carries its real states — load_tip_states copies every + // word, active or not — so the two rows usually differ. That only makes + // equality *harder* to reach: it costs a collapse flag, never invents one. + // One-sided safe; not worth a masked per-word compare. See red-team T-382. size_t sb = static_cast(s) * tw; size_t pb = static_cast(p) * tw; if (std::memcmp(&tree.prelim[sb], &tree.prelim[pb], word_bytes) != 0) continue; // --- Conditions 4–5 (NA only): down2 and subtree_actives preservation --- + // Unlike condition 3, these rows are written only for NA blocks — even + // load_tip_states skips subtree_actives for the rest — so a non-NA block's + // words stay at the zeros the arrays were sized with and always compare + // equal, dropping out of a test that has nothing to say about them. Here + // the unwritten words make equality *easier*, not harder (T-411). if (has_na) { if (std::memcmp(&tree.down2[sb], &tree.down2[pb], word_bytes) != 0) continue; diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index fcb6678c0..a190244b3 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -270,6 +270,15 @@ IntegerMatrix tree_to_collapsed_edge(const ts::TreeState& tree, // first-encountered child of each node goes left. ts::TreeState build_topology_tree(const IntegerMatrix& edge) { int n_edge = edge.nrow(); + // Same derivation, and so the same out-of-bounds writes, as init_from_edge. + // ncol is checked first: the child column is read as edge(i, 1), which on an + // n x 1 matrix indexes past the end of the underlying vector. + if (edge.ncol() != 2) { + stop("`tree` edge matrix must have exactly 2 columns."); + } + if (n_edge < 2 || !ts::edge_list_is_binary(&edge(0, 0), &edge(0, 1), n_edge)) { + stop("`tree` must be binary"); + } int n_tip = n_edge / 2 + 1; ts::TreeState tree; @@ -1733,6 +1742,13 @@ static int unpack_runtime(List rt, ts::DrivenParams& params) { flat[i] = se(i, 0); flat[n_edge + i] = se(i, 1); } + // init_from_edge refuses a non-binary tree by throwing, but under + // nThreads > 1 it runs on a worker thread, where an uncaught throw + // terminates the session. Reject here, on the main thread. + if (!ts::edge_list_is_binary(flat.data(), flat.data() + n_edge, + n_edge)) { + stop("Each `startEdge` matrix must describe a binary tree."); + } params.start_edges.push_back(std::move(flat)); } } diff --git a/src/ts_tree.cpp b/src/ts_tree.cpp index a6a14aae4..353b265b5 100644 --- a/src/ts_tree.cpp +++ b/src/ts_tree.cpp @@ -1,13 +1,53 @@ #include "ts_tree.h" #include #include +#include namespace ts { +bool edge_list_is_binary(const int* edge_parent, const int* edge_child, + int n_edge) { + // A rooted binary tree on n tips has 2 * (n - 1) edges, so an odd count + // cannot describe one; n_edge < 2 leaves no root to attach. + if (n_edge < 2 || (n_edge & 1)) return false; + const int n_tip = (n_edge / 2) + 1; + const int n_internal = n_tip - 1; + const int n_node = n_tip + n_internal; + + // Every non-root node must appear exactly once as a child and every internal + // node at most twice as a parent. n_edge == 2 * n_internal then forces + // "at most twice" to "exactly twice", which is binarity. A multifurcating + // edge list has more real tips than the n_tip derived above, so its extra + // tips are counted as internal nodes and parent no children at all. + std::vector child_of_an_edge(n_node, 0); + std::vector n_child(n_internal, 0); + for (int i = 0; i < n_edge; ++i) { + const int p = edge_parent[i] - 1; + const int c = edge_child[i] - 1; + if (p < n_tip || p >= n_node) return false; + if (c < 0 || c >= n_node || c == n_tip) return false; + if (child_of_an_edge[c]) return false; + child_of_an_edge[c] = 1; + if (++n_child[p - n_tip] > 2) return false; + } + return true; +} + void TreeState::init_from_edge( const int* edge_parent, const int* edge_child, int n_edge, const DataSet& ds) { + // Every count below is derived from n_edge on the assumption that the edge + // list is binary, and nothing downstream rechecks it. On a multifurcating + // list the derived n_tip falls short of the real tip count, so the loop + // writes past the end of parent[]/left[]/right[] (odd n_edge) or roots the + // tree on a real tip, leaving a one-element postorder whose downpass reads + // prelim[-total_words] (even n_edge). Refuse the tree instead. Rcpp + // forwards this to R as an error at every export boundary. + if (!edge_list_is_binary(edge_parent, edge_child, n_edge)) { + throw std::invalid_argument("`tree` must be binary"); + } + n_tip = (n_edge / 2) + 1; n_internal = n_tip - 1; n_node = n_tip + n_internal; @@ -294,6 +334,15 @@ void TreeState::reset_states(const DataSet& ds) { // subtree_a — only NA blocks; tips: load_tip_states + pass 2 update; // internals: pass 1 + pass 3 // local_cost— only standard blocks; written in pass 1 + // + // T-411: the collapse kernels (ts_collapsed.cpp) are a THIRD consumer the + // audit above does not cover. They compare whole rows by memcmp, so they + // also read words no pass ever writes: the SIMD pad word, and — for + // down2 / subtree_actives — the non-NA blocks of an NA dataset. Those read + // as zero only because every path that sizes a TreeState's state arrays + // zero-fills them (assign in init_from_edge, ts_sector.cpp and + // ts_constraint.cpp; resize on a fresh TreeState's empty vectors in + // ts_fuse.cpp). Re-audit the collapse kernels before relaxing that. load_tip_states(ds); } diff --git a/src/ts_tree.h b/src/ts_tree.h index 1a903f48f..236affbeb 100644 --- a/src/ts_tree.h +++ b/src/ts_tree.h @@ -196,6 +196,18 @@ struct TreeState { void reset_states(const DataSet& ds); }; +// True iff the 1-based edge list of `n_edge` rows has the degree spectrum of a +// rooted binary tree under the node convention above: every parent internal, +// every non-root node claimed as a child exactly once, each internal claiming +// two. That is what `init_from_edge` needs — it derives every node count from +// `n_edge` alone, and anything else makes it index past the end of +// parent[]/left[]/right[]. It is NOT full tree validation: a list satisfying +// it can still hold a cycle unreachable from the root, which `build_postorder` +// catches instead. Callers that can report an error more helpfully than the +// throw in `init_from_edge` should test with this first. +bool edge_list_is_binary(const int* edge_parent, const int* edge_child, + int n_edge); + } // namespace ts #endif // TS_TREE_H diff --git a/tests/testthat/test-ts-t400-multifurcating-guard.R b/tests/testthat/test-ts-t400-multifurcating-guard.R new file mode 100644 index 000000000..60284edc8 --- /dev/null +++ b/tests/testthat/test-ts-t400-multifurcating-guard.R @@ -0,0 +1,179 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +# T-400: `TreeState::init_from_edge()` derives every node count from the edge +# count alone, which holds only for a binary tree. A multifurcating tree wrote +# past the end of the topology arrays when the edge count was odd, and rooted +# the tree on a real tip when it was even -- leaving a one-element postorder +# whose downpass read the words immediately *before* the state buffer. Either +# way the caller got a plausible number rather than an error, so both parities +# are exercised below. + +# `(a,(b,((e,f),(g,h),(c,d))));` 8 tips, 6 internal nodes: 13 edges (odd) +.OddPolytomy <- function() { + ape::read.tree(text = "(a,(b,((e,f),(g,h),(c,d))));") +} + +# `(a,(b,(c,(d,(e,f,g,h)))));` 8 tips, 5 internal nodes: 12 edges (even) +.EvenPolytomy <- function() { + ape::read.tree(text = "(a,(b,(c,(d,(e,f,g,h)))));") +} + +.EightTaxonData <- function() { + MatrixToPhyDat(matrix( + c("0", "0", "1", "1", "0", "1", "0", "1", + "0", "1", "0", "1", "1", "1", "0", "0", + "1", "1", "1", "0", "0", "0", "1", "0", + "0", "0", "0", "1", "1", "0", "1", "1"), + nrow = 8, dimnames = list(letters[1:8], NULL))) +} + +test_that("The test polytomies have the edge counts the guard must handle", { + # Both parities must be covered: they failed by different mechanisms. + expect_equal(dim(.OddPolytomy()[["edge"]])[1], 13L) + expect_equal(dim(.EvenPolytomy()[["edge"]])[1], 12L) +}) + +test_that("TreeLength() rejects a multifurcating tree in a list", { + dat <- .EightTaxonData() + for (tr in list(.OddPolytomy(), .EvenPolytomy())) { + # Length-1 sets are the dangerous case: a heterogeneous set was already + # caught, accidentally, by the differing-edge-count check. + expect_error(TreeLength(structure(list(tr), class = "multiPhylo"), dat), + "must be binary") + expect_error(TreeLength(list(tr, tr), dat), "must be binary") + } +}) + +test_that("TreeLength() rejects a multifurcating tree through `[` and `[[`", { + dat <- .EightTaxonData() + tr <- .OddPolytomy() + trees <- structure(list(tr, tr), class = "multiPhylo") + # `[[` dispatches to the phylo method, which was already guarded; `[` keeps + # the multiPhylo class and reached the kernel. + expect_error(TreeLength(trees[[1]], dat), "must be binary") + expect_error(TreeLength(trees[1], dat), "must be binary") +}) + +test_that("Repeated scoring of one tree gives one answer", { + # The score was read from memory before the state buffer, so identical calls + # could disagree. Whatever the answer is, it must not vary between calls. + dat <- .EightTaxonData() + trees <- structure(list(.OddPolytomy()), class = "multiPhylo") + outcomes <- vapply(seq_len(5), function(i) { + tryCatch(paste(TreeLength(trees, dat), collapse = ","), + error = function(e) conditionMessage(e)) + }, character(1)) + expect_length(unique(outcomes), 1L) +}) + +test_that("CharacterLength() rejects a multifurcating tree", { + dat <- .EightTaxonData() + expect_error(CharacterLength(.OddPolytomy(), dat), "must be binary") + expect_error(CharacterLength(.EvenPolytomy(), dat), "must be binary") +}) + +test_that("TreeScore() and EdgeListScore() reject a multifurcating tree", { + dat <- PrepareData(.EightTaxonData()) + tr <- .OddPolytomy() + expect_error(TreeScore(tr, dat), "must be binary") + expect_error(EdgeListScore(tr[["edge"]][, 1], tr[["edge"]][, 2], dat), + "must be binary") + ev <- .EvenPolytomy() + expect_error(TreeScore(ev, dat), "must be binary") + expect_error(EdgeListScore(ev[["edge"]][, 1], ev[["edge"]][, 2], dat), + "must be binary") +}) + +test_that("The scoring kernel itself refuses a multifurcating edge matrix", { + # The R-level guards above are convenience; this is the boundary that every + # other kernel entry point sits behind. + dat <- .EightTaxonData() + tr <- RenumberTips(Renumber(.OddPolytomy()), names(dat)) + at <- attributes(dat) + tipData <- matrix(unlist(dat, use.names = FALSE), nrow = length(dat), + byrow = TRUE) + expect_error( + TreeSearch:::ts_fitch_score(tr[["edge"]], at[["contrast"]], tipData, + TreeSearch:::.ScaleWeight(at[["weight"]]), + at[["levels"]]), + "must be binary") +}) + +test_that("The kernel refuses malformed edge lists of binary length", { + # A polytomy is not the only edge list that would send the kernel out of + # bounds; these have the edge count of a four-tip binary tree (root = 5) but + # a shape it cannot index. + dat <- MatrixToPhyDat(matrix(c("0", "0", "1", "1", + "0", "1", "0", "1", + "1", "1", "0", "0"), + nrow = 4, dimnames = list(letters[1:4], NULL))) + at <- attributes(dat) + tipData <- matrix(unlist(dat, use.names = FALSE), nrow = 4, byrow = TRUE) + Score <- function(edge) { + TreeSearch:::ts_fitch_score(edge, at[["contrast"]], tipData, + TreeSearch:::.ScaleWeight(at[["weight"]]), + at[["levels"]]) + } + Edge <- function(...) matrix(c(...), ncol = 2, byrow = TRUE) + + # The same topology, accepted and scored as `(a,(b,(c,d)));` + expect_equal(Score(Edge(5, 1, 5, 6, 6, 2, 6, 7, 7, 3, 7, 4)), + TreeLength(ape::read.tree(text = "(a,(b,(c,d)));"), dat)) + # One node claimed as a child twice, leaving another with no parent + expect_error(Score(Edge(5, 1, 5, 6, 6, 2, 6, 7, 7, 3, 7, 3)), "must be binary") + # The root claimed as a child + expect_error(Score(Edge(5, 1, 5, 6, 6, 2, 6, 5, 7, 3, 7, 4)), "must be binary") + # A tip used as a parent + expect_error(Score(Edge(5, 1, 5, 6, 6, 2, 6, 7, 1, 3, 1, 4)), "must be binary") +}) + +test_that("A non-binary `startEdge` is refused on the main thread", { + # init_from_edge also runs on a search worker, where a throw would terminate + # the session rather than raise an R error, so the driven search screens + # start trees before dispatching. nThreads = 2 covers the threaded path. + dat <- MatrixToPhyDat(matrix( + c("0", "0", "0", "0", "0", "1", "1", "1", + "0", "0", "1", "1", "1", "0", "0", "1", + "0", "1", "0", "1", "1", "0", "1", "0"), + nrow = 8, dimnames = list(letters[1:8], NULL))) + ds <- make_ts_data(dat) + Driven <- function(edge, nThreads = 1L) { + TreeSearch:::ts_driven_search( + ds$contrast, ds$tip_data, ds$weight, ds$levels, + maxReplicates = 2L, ratchetCycles = 1L, verbosity = 0L, + nThreads = nThreads, startEdge = edge) + } + binary <- RenumberTips(Preorder(BalancedTree(letters[1:8])), + names(dat))[["edge"]] + expect_error(Driven(.OddPolytomy()[["edge"]]), "binary") + expect_error(Driven(.EvenPolytomy()[["edge"]]), "binary") + expect_error(Driven(.EvenPolytomy()[["edge"]], nThreads = 2L), "binary") + # A binary start must still be accepted, on both paths: a search has to have + # run and scored something, which an empty `scores` would not show. + for (nThreads in c(1L, 2L)) { + scores <- Driven(binary, nThreads = nThreads)[["scores"]] + expect_gt(length(scores), 0L) + expect_true(all(is.finite(scores))) + } +}) + +test_that("Binary trees are unaffected by the guard", { + dat <- .EightTaxonData() + pd <- PrepareData(dat) + for (tr in lapply(list(MakeTreeBinary(.OddPolytomy()), + MakeTreeBinary(.EvenPolytomy()), + BalancedTree(letters[1:8]), + PectinateTree(letters[1:8])), + Preorder)) { + score <- TreeLength(tr, dat) + expect_true(is.finite(score)) + expect_equal(unname(TreeLength(structure(list(tr), class = "multiPhylo"), + dat)), + score) + expect_equal(sum(CharacterLength(tr, dat, compress = TRUE) * + attr(dat, "weight")), + score) + expect_equal(TreeScore(RenumberTips(tr, names(dat)), pd), score) + } +})