fix: reject non-binary trees at the TreeState boundary

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

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]

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",
         ")"))
 

src/ts_collapsed.cpp

@@ -119,19 +119,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<size_t>(s) * tw;
     size_t pb = static_cast<size_t>(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;

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));
       }
     }

src/ts_tree.cpp

@@ -1,13 +1,53 @@
 #include "ts_tree.h"
 #include <algorithm>
 #include <cstring>
+#include <stdexcept>
 
 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<uint8_t> child_of_an_edge(n_node, 0);
+  std::vector<uint8_t> 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);
 }
 

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