From 25031ccafebbf890a4f7fa2c6f6b1e0a9532bd1e Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 04:24:26 +0100 Subject: [PATCH 1/4] fix: loosen constraint matching to the documented free-taxa contract MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `?MaximizeParsimony`'s `constraint` documents the phyDat reading: a tree is compliant when some edge separates the taxa coded `1` from those coded `0`, with `?`-coded taxa free on either side. The locked-node machinery enforced a strictly stronger one — some node's tip set had to EQUAL the `1` group (or its exact complement), free taxa excluded. Strict implies loose, so no wrong answer was ever returned. What broke was movement: a tree satisfying the documented contract without making either group an exact clade mapped to no node, which `regraft_violates_constraint()` reads as "already violating" and answers by rejecting every rearrangement. The replicate froze on its start. One reading, applied at every entry point: * `.PrepareConstraint()` now folds both groups into `consSplitMatrix` as 1 / 0 / NA, NA marking a free taxon. `build_constraint()` reads any value that is neither 1 nor 0 as free, so a hand-built 0/1 matrix (tests, `build_constraint_from_bitsets()`'s pool splits) still means "no free tips" and takes exactly the old path. * `map_constraint_nodes()` maps a split to the chain of nodes that DISPLAY it — covering one group, holding none of the other — instead of matching one exactly, and records both ends. Tips are candidates now too: a single-taxon group's "clade" is the tip itself, and scanning `postorder` alone (internal nodes only) left those splits unmapped and froze the replicate the same way. * `regraft_violates_constraint()` uses the end of that chain that permits most per question: the highest displaying node for a clip that must land inside, the tightest for one that must land outside. * `classify_clip_constraints()` reads "outside" as the apart-group rather than as `~split_tips`, so a clip of purely free taxa is UNCONSTRAINED and may be regrafted anywhere. * `wagner_tree_displays_constraint()` and `wagner_collect_active_splits()` get the same reading, so the Wagner build path cannot diverge from the search path again; `violates_constraint_posthoc()` already used it (a Fitch score against the constraint phyDat), so it is unchanged. * `impose_constraint()` no longer counts a documented-compliant tree as violating, and never moves a free taxon when it does repair one. * `ts_collapse_pool()` protects the tightest node realising each constraint rather than one matching it exactly, which with free taxa protected nothing. * A constraint character with no `0` taxa is dropped: it is vacuous under the documented contract, so enforcing its `1` group as a clade would restrict the search for nothing. `vignettes/search-algorithm.Rmd` gains a section stating the single contract. Fixes #54 Co-Authored-By: Claude Opus 5 --- R/MaximizeParsimony.R | 25 +- inst/WORDLIST | 1 + man/AdditionTree.Rd | 4 + man/MaximizeParsimony.Rd | 4 + man/Resample.Rd | 4 + man/SuccessiveApproximations.Rd | 4 + src/ts_constraint.cpp | 312 +++++++++++------- src/ts_constraint.h | 62 +++- src/ts_rcpp.cpp | 70 ++-- src/ts_wagner.cpp | 84 +++-- tests/testthat/test-ts-constraint-free-taxa.R | 192 +++++++++++ vignettes/search-algorithm.Rmd | 30 ++ 12 files changed, 604 insertions(+), 188 deletions(-) create mode 100644 tests/testthat/test-ts-constraint-free-taxa.R diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 3ccce1f44..d0f8842dd 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -145,10 +145,14 @@ } } - keep <- apply(consSplits, 1, function(row) { - s <- sum(row) - s >= 1 && s < length(constraint) - 1 - }) + # A character only constrains anything when both groups are occupied: with no + # "0" tips there is no edge for the "1" tips to be separated *from*, so the + # documented contract ("some edge separates the 1 taxa from the 0 taxa") is + # vacuously true and enforcing the group as a clade would be a restriction the + # user never asked for. + nOne <- rowSums(consSplits) + nZero <- rowSums(consZero) + keep <- nOne >= 1 & nZero >= 1 & nOne < length(constraint) - 1 consSplits <- consSplits[keep, , drop = FALSE] consZero <- consZero[keep, , drop = FALSE] if (nrow(consSplits) == 0L) return(list()) @@ -188,6 +192,15 @@ consTipData <- matrix(unlist(constraint, use.names = FALSE), nrow = length(constraint), byrow = TRUE) + # Fold the two groups into the single membership matrix the C++ engine reads: + # 1 = "together", 0 = "apart", NA = free to fall on either side. A tip that + # is in neither group must not be coded 0, or the engine would enforce the + # stricter "the 1 group is an exact clade" reading and refuse to move a start + # tree that already satisfies the documented one (agent-issues/TreeSearch#54). + # build_constraint() (src/ts_constraint.cpp) treats any value that is neither + # 1 nor 0 as free, so a plain 0/1 matrix still means "no free tips". + consSplits[consSplits == 0L & consZero == 0L] <- NA_integer_ + list( consSplitMatrix = consSplits, consContrast = consContrast, @@ -694,6 +707,10 @@ #' returned trees will be perfectly compatible with each character in #' `constraint`; or a tree of class `phylo`, all of whose nodes will occur #' in any output tree. +#' A returned tree is compatible with a constraint character when some edge +#' separates the taxa coded `1` from those coded `0`. Taxa coded `?`, and taxa +#' that `constraint` does not mention, are unconstrained: they may fall on +#' either side of that edge, and are not required to join either group. #' Constraint searches are supported natively: all tree rearrangements #' are filtered to respect the constraint topology. #' @param effort Integer: how much search effort to spend, **relative to the diff --git a/inst/WORDLIST b/inst/WORDLIST index 940b38b7e..0765f408e 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -266,6 +266,7 @@ rearranger reconverged reconverges regraft +regrafted regrafting regrafts reoptimisation diff --git a/man/AdditionTree.Rd b/man/AdditionTree.Rd index 13fbc9e53..380bea854 100644 --- a/man/AdditionTree.Rd +++ b/man/AdditionTree.Rd @@ -35,6 +35,10 @@ construction begins.} returned trees will be perfectly compatible with each character in \code{constraint}; or a tree of class \code{phylo}, all of whose nodes will occur in any output tree. +A returned tree is compatible with a constraint character when some edge +separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa +that \code{constraint} does not mention, are unconstrained: they may fall on +either side of that edge, and are not required to join either group. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 917751868..28c82714b 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -137,6 +137,10 @@ block. Only used when \code{inapplicable = "hsj"}.} returned trees will be perfectly compatible with each character in \code{constraint}; or a tree of class \code{phylo}, all of whose nodes will occur in any output tree. +A returned tree is compatible with a constraint character when some edge +separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa +that \code{constraint} does not mention, are unconstrained: they may fall on +either side of that edge, and are not required to join either group. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/Resample.Rd b/man/Resample.Rd index e99630af0..329b065bf 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -71,6 +71,10 @@ Specify \code{"profile"} to employ profile parsimony returned trees will be perfectly compatible with each character in \code{constraint}; or a tree of class \code{phylo}, all of whose nodes will occur in any output tree. +A returned tree is compatible with a constraint character when some edge +separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa +that \code{constraint} does not mention, are unconstrained: they may fall on +either side of that edge, and are not required to join either group. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/SuccessiveApproximations.Rd b/man/SuccessiveApproximations.Rd index cdd4e3e8c..864bd5481 100644 --- a/man/SuccessiveApproximations.Rd +++ b/man/SuccessiveApproximations.Rd @@ -83,6 +83,10 @@ Specify \code{"profile"} to employ profile parsimony returned trees will be perfectly compatible with each character in \code{constraint}; or a tree of class \code{phylo}, all of whose nodes will occur in any output tree. +A returned tree is compatible with a constraint character when some edge +separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa +that \code{constraint} does not mention, are unconstrained: they may fall on +either side of that edge, and are not required to join either group. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/src/ts_constraint.cpp b/src/ts_constraint.cpp index 114dc543b..908f7d0f6 100644 --- a/src/ts_constraint.cpp +++ b/src/ts_constraint.cpp @@ -25,31 +25,33 @@ ConstraintData build_constraint( cd.split_tips.resize( static_cast(n_splits) * cd.n_words, 0ULL); + cd.split_zeros.resize( + static_cast(n_splits) * cd.n_words, 0ULL); cd.constraint_node.assign(n_splits, -1); + cd.constraint_node_hi.assign(n_splits, -1); cd.constraint_complement.assign(n_splits, 0); - // Pack split_matrix rows into bitmasks. + // Pack split_matrix rows into a pair of bitmasks. // split_matrix is column-major (from R): element [s, t] is at - // index s + n_splits * t. + // index s + n_splits * t. 1 -> "together" group, 0 -> "apart" group, + // anything else (NA_INTEGER) -> free, in neither mask (T-386). for (int s = 0; s < n_splits; ++s) { - uint64_t* mask = &cd.split_tips[static_cast(s) * cd.n_words]; + uint64_t* ones = &cd.split_tips[static_cast(s) * cd.n_words]; + uint64_t* zeros = &cd.split_zeros[static_cast(s) * cd.n_words]; for (int t = 0; t < n_tips; ++t) { - if (split_matrix[s + n_splits * t]) { - int w = t / 64; - int b = t % 64; - mask[w] |= (1ULL << b); - } + const int v = split_matrix[s + n_splits * t]; + if (v != 1 && v != 0) continue; // free tip + int w = t / 64; + int b = t % 64; + (v == 1 ? ones : zeros)[w] |= (1ULL << b); } - // Canonicalize: tip 0 must be on the "outside" (bit 0 = 0). - // If bit 0 is set, flip the entire mask. - if (mask[0] & 1ULL) { + // Canonicalize: tip 0 must be outside split_tips (bit 0 = 0). The two + // groups of a bipartition are interchangeable, so SWAP them rather than + // complementing either — complementing would swallow the free tips into + // the "apart" group and reinstate the exact-clade reading. + if (ones[0] & 1ULL) { for (int w = 0; w < cd.n_words; ++w) { - mask[w] = ~mask[w]; - } - // Clear bits beyond n_tips - int remainder = n_tips % 64; - if (remainder > 0) { - mask[cd.n_words - 1] &= (1ULL << remainder) - 1; + std::swap(ones[w], zeros[w]); } } } @@ -83,7 +85,23 @@ ConstraintData build_constraint_from_bitsets( // Copy split data size_t total = static_cast(n_splits) * words_per_split; cd.split_tips.assign(split_bits, split_bits + total); + // These splits come from pool bipartitions, which partition every tip: there + // are no free tips, so the "apart" group is exactly the complement and the + // free-taxa machinery collapses back to the exact-clade test (T-386). + cd.split_zeros.assign(total, 0ULL); + { + const int rem = n_tips % 64; + const uint64_t top = rem ? ((1ULL << rem) - 1ULL) : ~0ULL; + for (int s = 0; s < n_splits; ++s) { + const size_t off = static_cast(s) * words_per_split; + for (int w = 0; w < words_per_split; ++w) { + cd.split_zeros[off + w] = ~cd.split_tips[off + w]; + if (w == words_per_split - 1) cd.split_zeros[off + w] &= top; + } + } + } cd.constraint_node.assign(n_splits, -1); + cd.constraint_node_hi.assign(n_splits, -1); cd.constraint_complement.assign(n_splits, 0); int n_node = 2 * n_tips - 1; @@ -155,26 +173,86 @@ std::vector compute_node_tips(const TreeState& tree, int n_words) // Map constraint nodes: find which internal node holds each split // ========================================================================= -// Width mask for the highest word of a tip bitmask: node tip sets carry zeros -// above tip n_tip - 1, so a *complemented* split mask has to be trimmed to the -// same width before it can be compared with one. -static inline uint64_t tip_mask_top_word(int n_tip) { - const int rem = n_tip % 64; - return rem ? ((1ULL << rem) - 1ULL) : ~0ULL; +// Does the edge above `node` separate `together` from `apart`? It does when +// the node's descendant tip set covers every tip of `together` and holds none +// of `apart`; the tips in neither group are free and are not looked at (T-386). +// +// With `apart` the exact complement of `together` — a constraint with no free +// tips, and every split built by build_constraint_from_bitsets() — the two +// conditions together force set equality, which is the exact-clade test this +// replaced. +static inline bool node_displays_split( + const uint64_t* nd, const uint64_t* together, const uint64_t* apart, + int n_words) +{ + for (int w = 0; w < n_words; ++w) { + if ((together[w] & ~nd[w]) != 0ULL) return false; // a required tip missing + if ((apart[w] & nd[w]) != 0ULL) return false; // an excluded tip present + } + return true; } -// Does node `node`'s descendant tip set equal `split` (complement = false) or -// the complement of `split` over tips 0..n_tip-1 (complement = true)? -static inline bool node_matches_split( - const uint64_t* nd, const uint64_t* split, int n_words, - uint64_t top_word, bool complement) +// Tightest and highest node displaying `together` | `apart`, or {-1, -1}. +// +// Every node that displays the split covers `together`, so all of them are +// ancestors of LCA(together) and they form one unbroken upward chain: each +// step up adds tips, and the moment a step adds a tip of `apart` the chain +// ends (higher nodes keep it). So the tight end is the first match in an +// order that visits descendants before ancestors, and the high end is found by +// walking parents from there. Tips are candidates too, for the single-taxon +// group whose "clade" is the tip itself — tree.postorder holds only internal +// nodes, so scanning it alone left those splits unmapped, which +// regraft_violates_constraint() reads as "already violating". +static void find_displaying_chain( + const TreeState& tree, const std::vector& node_tips, + const uint64_t* together, const uint64_t* apart, int n_words, + int& lo, int& hi) { + lo = -1; + hi = -1; + + // Tip candidates without scanning the tips: a tip's set is the singleton + // {t}, so it can only cover `together` when `together` is {t} itself (or, + // degenerately, empty — then the lowest tip outside `apart` wins). + int n_together = 0, lone_together = -1; for (int w = 0; w < n_words; ++w) { - uint64_t want = complement ? ~split[w] : split[w]; - if (complement && w == n_words - 1) want &= top_word; - if (nd[w] != want) return false; + if (together[w]) { + n_together += popcount64(together[w]); + lone_together = w * 64 + ctz64(together[w]); + } + } + if (n_together == 1) { + const uint64_t* nd = &node_tips[static_cast(lone_together) * n_words]; + if (node_displays_split(nd, together, apart, n_words)) lo = lone_together; + } else if (n_together == 0) { + for (int w = 0; w < n_words && lo < 0; ++w) { + uint64_t free_here = ~apart[w]; + const int lim = tree.n_tip - w * 64; + if (lim < 64) free_here &= (1ULL << lim) - 1ULL; + if (free_here) lo = w * 64 + ctz64(free_here); + } + } + if (lo < 0) { + for (int node : tree.postorder) { + const uint64_t* nd = &node_tips[static_cast(node) * n_words]; + if (node_displays_split(nd, together, apart, n_words)) { lo = node; break; } + } + } + if (lo < 0) return; + + // Walk to the top of the chain. Bounded by n_node rather than trusting the + // root to be reachable: impose_one_pass() calls this on trees it is midway + // through repairing, and a parent-ascending loop over a corrupt parent[] is + // exactly the hang T-327/T-333 had to be defended against elsewhere. + hi = lo; + const int root = tree.n_tip; + for (int guard = 0; guard < tree.n_node && hi != root; ++guard) { + const int up = tree.parent[hi]; + if (up < 0 || up >= tree.n_node || up == hi) break; + const uint64_t* nd = &node_tips[static_cast(up) * n_words]; + if (!node_displays_split(nd, together, apart, n_words)) break; + hi = up; } - return true; } void map_constraint_nodes(const TreeState& tree, ConstraintData& cd) @@ -182,7 +260,6 @@ void map_constraint_nodes(const TreeState& tree, ConstraintData& cd) if (!cd.active) return; auto node_tips = compute_node_tips(tree, cd.n_words); - const uint64_t top_word = tip_mask_top_word(tree.n_tip); // For each constraint split, find the node that displays it. // @@ -190,38 +267,35 @@ void map_constraint_nodes(const TreeState& tree, ConstraintData& cd) // a rooted subtree, so the split is displayed whenever EITHER side is a // clade. Exactly one of the two is, except when the split is the root's own // bipartition (then both are): for an edge (parent(v), v) with v != root the - // two sides are desc(v) and its complement, so a tree displays A|B iff some - // node's tip set equals A or equals B. build_constraint() canonicalises A so - // that tip 0 is outside it, which makes A the clade side only when tip 0 sits - // on the root's own edge -- true of a tip-0-rooted tree and of nothing else. - // Testing the complement as well is what makes this mapping rooting-agnostic, - // and it costs one extra scan only for splits that used to map to -1 (which - // regraft_violates_constraint reads as "tree already violates", rejecting - // every move). Phase 1 is run to completion first so that every tree which - // mapped successfully before maps to exactly the same node now. + // two sides are desc(v) and its complement. build_constraint() canonicalises + // A so that tip 0 is outside it, which makes A the clade side only when tip 0 + // sits on the root's own edge -- true of a tip-0-rooted tree and of nothing + // else. Testing the complement as well is what makes this mapping + // rooting-agnostic, and it costs one extra scan only for splits that used to + // map to -1 (which regraft_violates_constraint reads as "tree already + // violates", rejecting every move). Phase 1 is run to completion first so + // that every tree which mapped successfully before maps to exactly the same + // node now. + // + // T-386: "is a clade" is the free-taxa reading, not set equality -- see + // node_displays_split(). The chain of displaying nodes is recorded at both + // ends, because a regraft that must land INSIDE the constrained group may use + // the whole chain while one that must land outside may not; see + // regraft_violates_constraint(). for (int s = 0; s < cd.n_splits; ++s) { - const uint64_t* split = &cd.split_tips[static_cast(s) * cd.n_words]; - cd.constraint_node[s] = -1; + const uint64_t* ones = &cd.split_tips[static_cast(s) * cd.n_words]; + const uint64_t* zeros = &cd.split_zeros[static_cast(s) * cd.n_words]; cd.constraint_complement[s] = 0; - for (int node : tree.postorder) { - const uint64_t* nd = &node_tips[static_cast(node) * cd.n_words]; - if (node_matches_split(nd, split, cd.n_words, top_word, false)) { - cd.constraint_node[s] = node; - break; - } - } - if (cd.constraint_node[s] >= 0) continue; - - // Phase 2: the tip-0 side is the clade in this rooting. - for (int node : tree.postorder) { - const uint64_t* nd = &node_tips[static_cast(node) * cd.n_words]; - if (node_matches_split(nd, split, cd.n_words, top_word, true)) { - cd.constraint_node[s] = node; - cd.constraint_complement[s] = 1; - break; - } + int lo = -1, hi = -1; + find_displaying_chain(tree, node_tips, ones, zeros, cd.n_words, lo, hi); + if (lo < 0) { + // Phase 2: the tip-0 side is the clade in this rooting. + find_displaying_chain(tree, node_tips, zeros, ones, cd.n_words, lo, hi); + if (lo >= 0) cd.constraint_complement[s] = 1; } + cd.constraint_node[s] = lo; + cd.constraint_node_hi[s] = hi; } } @@ -318,15 +392,24 @@ void classify_clip_constraints(const TreeState& tree, int clip_node, compute_clip_tip_mask(tree, clip_node, cd.clip_tip_mask); + // "Inside" means the clip holds a tip of the group that must stay together; + // "outside", a tip of the group that must stay apart from it. A clip made + // only of FREE tips is in neither, and lands here as UNCONSTRAINED — the + // whole point of the free-taxa reading, and what lets such a clip be + // regrafted anywhere without breaking the separating edge (T-386). Reading + // "outside" as ~split, which is what this did before free tips existed, + // pinned every free tip to the far side of the constraint. for (int s = 0; s < cd.n_splits; ++s) { - const uint64_t* split = + const uint64_t* ones = &cd.split_tips[static_cast(s) * cd.n_words]; + const uint64_t* zeros = + &cd.split_zeros[static_cast(s) * cd.n_words]; bool any_inside = false; bool any_outside = false; for (int w = 0; w < cd.n_words; ++w) { - if (cd.clip_tip_mask[w] & split[w]) any_inside = true; - if (cd.clip_tip_mask[w] & ~split[w]) any_outside = true; + if (cd.clip_tip_mask[w] & ones[w]) any_inside = true; + if (cd.clip_tip_mask[w] & zeros[w]) any_outside = true; if (any_inside && any_outside) break; } @@ -339,23 +422,10 @@ void classify_clip_constraints(const TreeState& tree, int clip_node, bool rest_has_in = false; bool rest_has_out = false; for (int w = 0; w < cd.n_words; ++w) { - uint64_t rest = ~cd.clip_tip_mask[w]; - // Mask out bits beyond n_tips in the last word - if (w == cd.n_words - 1) { - int remainder = tree.n_tip % 64; - if (remainder > 0) - rest &= (1ULL << remainder) - 1; - } - if (rest & split[w]) rest_has_in = true; - if (rest & ~split[w]) { - uint64_t out_bits = ~split[w]; - if (w == cd.n_words - 1) { - int remainder = tree.n_tip % 64; - if (remainder > 0) - out_bits &= (1ULL << remainder) - 1; - } - if (rest & out_bits) rest_has_out = true; - } + const uint64_t rest = ~cd.clip_tip_mask[w]; + // No width mask needed: ones/zeros carry zeros above tip n_tip - 1. + if (rest & ones[w]) rest_has_in = true; + if (rest & zeros[w]) rest_has_out = true; } if (rest_has_in && rest_has_out) { cd.clip_zones[s] = ClipZone::FORBIDDEN; @@ -396,13 +466,14 @@ bool regraft_violates_constraint(int below, // can preserve this split — reject unconditionally. if (cd.clip_zones[s] == ClipZone::FORBIDDEN) return true; - int cn = cd.constraint_node[s]; - if (cn < 0) { + const int cn_lo = cd.constraint_node[s]; + if (cn_lo < 0) { // Constraint genuinely not displayed by the current tree (both sides // tested — see map_constraint_nodes). Reject all moves to avoid // entrenching a bad state. return true; } + const int cn_hi = cd.constraint_node_hi[s]; // Which side of the split does cn's subtree hold? Under the canonical // orientation it is the split itself; in a rooting where only the tip-0 @@ -414,18 +485,27 @@ bool regraft_violates_constraint(int below, const ClipZone zone_out = cd.constraint_complement[s] ? ClipZone::MUST_INSIDE : ClipZone::MUST_OUTSIDE; - // Is `below` a descendant of cn (= inside the mapped clade)? - bool inside = is_ancestor_or_equal(cn, below, - cd.dfs_entry, cd.dfs_exit); - - if (cd.clip_zones[s] == zone_in && !inside) { + // The two ends of the displaying chain answer two different questions, and + // each wants the end that permits most (T-386; with no free tips the chain + // is one node long and both reduce to the pre-T-386 test): + // + // * a clip that must land INSIDE carries tips of the together-group but + // none of the apart-group, so anywhere within the HIGHEST displaying + // node keeps that node covering the group and free of the other. Only + // above cn_hi does the enclosing node pick up an apart-group tip. + // * a clip that must land OUTSIDE carries apart-group tips, so it may go + // anywhere that leaves some displaying node intact — and the TIGHTEST + // is the one hardest to contaminate, so it forbids least. + if (cd.clip_zones[s] == zone_in && + !is_ancestor_or_equal(cn_hi, below, cd.dfs_entry, cd.dfs_exit)) { return true; } // Exclude the boundary edge (above_cn, cn): regrafting an outside-only // clade just above the constraint clade makes it a sibling of that clade, // preserving monophyly. Only reject if the clade would land *strictly // inside* the constraint clade. - if (cd.clip_zones[s] == zone_out && inside && below != cn) { + if (cd.clip_zones[s] == zone_out && below != cn_lo && + is_ancestor_or_equal(cn_lo, below, cd.dfs_entry, cd.dfs_exit)) { return true; } } @@ -715,24 +795,21 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, // already displayed every constraint, spending up to n_tip / 4 + 2 arbitrary // SPR moves on it. That mattered most at ts_nni_perturb.cpp's unconditional // impose_constraint() call, which runs after every perturbation cycle. - // The repair below still aims at making the canonical side the clade, which - // displays the split either way. - const uint64_t top_word = tip_mask_top_word(tree.n_tip); + // "Is a clade" is the free-taxa reading (T-386), so a tree that satisfies + // what `constraint` documents is likewise left alone. The repair below aims + // at making the canonical side a clade, which displays the split either way. std::vector violated; for (int s = 0; s < cd.n_splits; ++s) { - const uint64_t* split = + const uint64_t* ones = &cd.split_tips[static_cast(s) * n_words]; - bool found = false; - for (int node : tree.postorder) { - const uint64_t* nd = - &node_tips[static_cast(node) * n_words]; - if (node_matches_split(nd, split, n_words, top_word, false) || - node_matches_split(nd, split, n_words, top_word, true)) { - found = true; - break; - } + const uint64_t* zeros = + &cd.split_zeros[static_cast(s) * n_words]; + int lo = -1, hi = -1; + find_displaying_chain(tree, node_tips, ones, zeros, n_words, lo, hi); + if (lo < 0) { + find_displaying_chain(tree, node_tips, zeros, ones, n_words, lo, hi); } - if (!found) violated.push_back(s); + if (lo < 0) violated.push_back(s); } if (violated.empty()) return 0; @@ -754,10 +831,25 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, int total_moves = 0; + // Tips that keep node `nd` from displaying the split: those of the + // together-group it is missing, plus those of the apart-group it holds. + // Free tips appear in neither, so the repair never moves one (T-386) — they + // may sit on whichever side they already do. + auto repair_cost = [&](const uint64_t* nd, const uint64_t* ones, + const uint64_t* zeros) { + int cost = 0; + for (int w = 0; w < n_words; ++w) { + cost += popcount64(ones[w] & ~nd[w]) + popcount64(zeros[w] & nd[w]); + } + return cost; + }; + for (size_t vi = 0; vi < violated.size(); ++vi) { int s = violated[vi]; const uint64_t* split = &cd.split_tips[static_cast(s) * n_words]; + const uint64_t* split_out = + &cd.split_zeros[static_cast(s) * n_words]; // Rebuild bitmasks after previous split's moves if (vi > 0) { @@ -765,16 +857,13 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, node_tips = compute_node_tips(tree, n_words); } - // --- Find best candidate node (min symmetric difference) --- + // --- Find best candidate node (fewest misplaced tips) --- int best_node = -1; int best_cost = tree.n_tip + 1; for (int node : tree.postorder) { const uint64_t* nd = &node_tips[static_cast(node) * n_words]; - int cost = 0; - for (int w = 0; w < n_words; ++w) { - cost += popcount64(nd[w] ^ split[w]); - } + int cost = repair_cost(nd, split, split_out); if (cost < best_cost) { best_cost = cost; best_node = node; @@ -789,7 +878,7 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, const uint64_t* best_nt = &node_tips[static_cast(best_node) * n_words]; for (int w = 0; w < n_words; ++w) { - move_out_mask[w] = best_nt[w] & ~split[w]; + move_out_mask[w] = best_nt[w] & split_out[w]; move_in_mask[w] = split[w] & ~best_nt[w]; } @@ -839,10 +928,7 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, int bc = tree.n_tip + 1; for (int node : tree.postorder) { const uint64_t* nd = &nt[static_cast(node) * n_words]; - int cost = 0; - for (int w = 0; w < n_words; ++w) { - cost += popcount64(nd[w] ^ split[w]); - } + int cost = repair_cost(nd, split, split_out); if (cost < bc) { bc = cost; bn = node; } } return bn; diff --git a/src/ts_constraint.h b/src/ts_constraint.h index 44e345103..ee4a16061 100644 --- a/src/ts_constraint.h +++ b/src/ts_constraint.h @@ -3,10 +3,20 @@ // Topological constraint enforcement (TNT-style locked nodes). // -// A constraint is a set of splits (bipartitions). A tree satisfies the -// constraint iff every constraint split is displayed — i.e., for each -// split there is an internal node whose subtree tip set matches (after -// accounting for unconstrained taxa that may sit on either side). +// A constraint is a set of splits. Each split names two disjoint groups of +// tips — the "1" group and the "0" group of one constraint character — and any +// remaining tips are FREE: coded `?`, or absent from the constraint phyDat +// altogether. A tree satisfies the split iff some edge separates the 1 group +// from the 0 group, free tips falling on either side. That is the contract +// `?MaximizeParsimony`'s `constraint` argument documents, and since T-386 it is +// the one enforced here: a node DISPLAYS the split when its descendant tip set +// is a superset of one group and disjoint from the other. +// +// Requiring the tip set to EQUAL a group (the pre-T-386 test) is strictly +// stronger. It never returned a wrong answer, but a start tree that satisfied +// the documented contract without making either group an exact clade mapped to +// no node at all, which regraft_violates_constraint() reads as "the tree +// already violates" — freezing the replicate on its start. // // Implementation: // 1. At init: store constraint splits as tip bitmasks. @@ -36,18 +46,37 @@ struct ConstraintData { int n_splits = 0; int n_words = 0; // ceil(n_tips / 64) - // Tip bitmasks: split_tips[i * n_words .. (i+1) * n_words - 1] + // Tip bitmasks: split_tips[i * n_words .. (i+1) * n_words - 1]. + // The tips that must end up TOGETHER, on one side of some edge. // Canonical: bit 0 (tip 0) is always on the "outside" (= 0). std::vector split_tips; - // Current mapping: constraint_node[i] = the internal node whose - // subtree tips match split i in the current tree. + // The tips that must end up on the OTHER side of that edge, same layout. + // Disjoint from split_tips; the two need NOT be complements — a tip in + // neither mask is free to fall on either side (T-386). When the caller + // supplies no free tips this is exactly ~split_tips, and every check below + // reduces to the pre-T-386 exact-clade test. + std::vector split_zeros; + + // Current mapping: constraint_node[i] = the TIGHTEST node (tip or internal) + // that displays split i in the current tree — its descendant tip set covers + // one of the two groups and avoids the other. // -1 if not yet mapped (or the tree does not display split i). std::vector constraint_node; + // The HIGHEST node that displays split i, in the same polarity as + // constraint_node[i]; equal to it when no free tip sits directly above. + // The displaying nodes form an unbroken chain from constraint_node[i] up to + // this one (each step adds only free tips), so the two ends are all a + // regraft test needs — see regraft_violates_constraint(), which uses this + // end for "must land inside" and the tight end for "must land outside". + // -1 exactly when constraint_node[i] is. + std::vector constraint_node_hi; + // Polarity of constraint_node[i] (T-384). 0: the node's descendant tip set - // is split_tips[i] itself. 1: it is the *complement* of split_tips[i], i.e. - // the tip-0 side of the bipartition. A constraint split is an UNROOTED + // covers split_tips[i] and avoids split_zeros[i]. 1: the other way round — + // it covers split_zeros[i], the tip-0 side of the bipartition, and avoids + // split_tips[i]. A constraint split is an UNROOTED // bipartition, so a tree displays it whenever EITHER side is a rooted clade, // and which side that is depends on the rooting alone -- see // map_constraint_nodes(). Consumers that treat constraint_node[i] as "the @@ -73,9 +102,18 @@ struct ConstraintData { std::vector clip_tip_mask; // [n_words] }; -// Build ConstraintData from R-side split bitmask matrix. -// split_matrix: n_splits x n_tips, each row is 0/1 indicating split membership. -// The matrix is canonicalized so tip 0 is always "outside" (= 0). +// Build ConstraintData from R-side split membership matrix. +// split_matrix: n_splits x n_tips, column-major. Element [s, t] is +// 1 tip t is in split s's "together" group; +// 0 tip t is in split s's "apart" group; +// anything else (NA_INTEGER, as .PrepareConstraint() writes for a `?`-coded +// or unconstrained taxon) — tip t is FREE, and may fall on either side. +// A pure 0/1 matrix therefore means "no free tips", i.e. the exact-clade +// reading that predates T-386; callers that build one by hand keep it. +// The two groups are swapped where needed so that tip 0 is never in +// split_tips ("outside" the canonical side) — the same invariant the Wagner +// and pool paths have always relied on, and harmless because a split is an +// unrooted bipartition whose two groups are interchangeable. ConstraintData build_constraint( const int* split_matrix, int n_splits, int n_tips); diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 8ff590920..0dfae2af7 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2235,25 +2235,29 @@ List ts_collapse_pool( // must stay visible even when its branch is zero-length (it would otherwise be // contracted, leaving the result looking unconstrained). We protect every // constraint split from collapse — the unsupported NON-constraint branches - // still collapse. Store each constraint split as a canonical (tip-0-excluded) - // bitset: trees are re-rooted on tip 0 below, so every internal node's - // descendant set excludes tip 0 and is directly comparable to these. + // still collapse. + // + // Each split is stored as the pair of groups build_constraint() reads + // (1 = together, 0 = apart, anything else = free; see ts_constraint.cpp), not + // as one canonical bitset: with free tips the enforced grouping is generally + // NOT any node's exact tip set, and the exact-match test this replaced then + // protected nothing at all (agent-issues/TreeSearch#54). A pure 0/1 matrix + // still gives apart == the complement, and the test below still fires on + // exactly the node the equality test used to find. const int n_tip = tip_data.nrow(); const int wps = (n_tip + 63) / 64; - std::vector> cons_canon; + std::vector> cons_one, cons_zero; if (consSplitMatrix.isNotNull()) { IntegerMatrix cs(consSplitMatrix.get()); for (int r = 0; r < cs.nrow(); ++r) { - std::vector b(wps, 0); + std::vector one(wps, 0), zero(wps, 0); for (int c = 0; c < n_tip && c < cs.ncol(); ++c) { - if (cs(r, c)) b[c >> 6] |= (1ULL << (c & 63)); + const int v = cs(r, c); + if (v != 1 && v != 0) continue; // free tip + (v == 1 ? one : zero)[c >> 6] |= (1ULL << (c & 63)); } - if (b[0] & 1ULL) { // canonicalize: exclude tip 0 - for (int w = 0; w < wps; ++w) b[w] = ~b[w]; - int rem = n_tip & 63; - if (rem) b[wps - 1] &= ((1ULL << rem) - 1); // clear padding bits - } - cons_canon.push_back(std::move(b)); + cons_one.push_back(std::move(one)); + cons_zero.push_back(std::move(zero)); } } @@ -2306,11 +2310,14 @@ List ts_collapse_pool( ts::compute_collapsed_flags_aggressive(tree, ds, flags); - // Protect constraint splits: clear the collapse flag of any internal edge - // whose bipartition realises a constraint (keeps the enforced clade - // visible). Per-node descendant tip sets via a postorder OR; rooted on - // tip 0, so every internal set excludes tip 0 == the canonical form above. - if (!cons_canon.empty()) { + // Protect constraint splits: clear the collapse flag of the internal edge + // that realises each constraint (keeps the enforced clade visible). That + // is the TIGHTEST node displaying the split — collapsing it is what would + // hide the grouping, whereas the looser nodes above it (which differ only + // by free tips) show nothing the tight one does not. Per-node descendant + // tip sets via a postorder OR; rooted on tip 0, so every internal set + // excludes tip 0 and only one of the two groups can be the clade side. + if (!cons_one.empty()) { std::vector tb(static_cast(tree.n_node) * wps, 0); for (int tp = 0; tp < n_tip; ++tp) { tb[static_cast(tp) * wps + (tp >> 6)] = 1ULL << (tp & 63); @@ -2324,15 +2331,26 @@ List ts_collapse_pool( const uint64_t* R = &tb[static_cast(tree.right[ni]) * wps]; for (int w = 0; w < wps; ++w) dst[w] = L[w] | R[w]; } - for (int v = n_tip + 1; v < tree.n_node; ++v) { - if (v >= static_cast(flags.size()) || !flags[v]) continue; - const uint64_t* nb = &tb[static_cast(v) * wps]; - for (const auto& cb : cons_canon) { - bool eq = true; - for (int w = 0; w < wps; ++w) { - if (nb[w] != cb[w]) { eq = false; break; } - } - if (eq) { flags[v] = 0; break; } + auto displays = [&](const uint64_t* nb, const std::vector& in, + const std::vector& out) { + for (int w = 0; w < wps; ++w) { + if (in[w] & ~nb[w]) return false; // a required tip missing + if (out[w] & nb[w]) return false; // an excluded tip present + } + return true; + }; + for (size_t ci = 0; ci < cons_one.size(); ++ci) { + int tight = -1, tight_size = n_tip + 1; + for (int v = n_tip + 1; v < tree.n_node; ++v) { + const uint64_t* nb = &tb[static_cast(v) * wps]; + if (!displays(nb, cons_one[ci], cons_zero[ci]) && + !displays(nb, cons_zero[ci], cons_one[ci])) continue; + int sz = 0; + for (int w = 0; w < wps; ++w) sz += ts::popcount64(nb[w]); + if (sz < tight_size) { tight_size = sz; tight = v; } + } + if (tight >= 0 && tight < static_cast(flags.size())) { + flags[tight] = 0; } } } diff --git a/src/ts_wagner.cpp b/src/ts_wagner.cpp index 6607a36dd..db59a118a 100644 --- a/src/ts_wagner.cpp +++ b/src/ts_wagner.cpp @@ -347,20 +347,19 @@ static int wagner_smallest_containing_node( // further insertion can make it. static int wagner_map_complement( const TreeState& tree, int n_tip, int nw, - const std::vector& node_tips, const uint64_t* split, + const std::vector& node_tips, const uint64_t* split_out, const std::vector& added_tips, WagnerConstraintScratch& scratch) { + (void)n_tip; uint64_t* needed_out = scratch.needed_out.data(); int n_out = 0; int lone_out = -1; for (int w = 0; w < nw; ++w) { - uint64_t outside_mask = ~split[w]; - if (w == nw - 1) { - const int rem = n_tip % 64; - if (rem > 0) outside_mask &= (1ULL << rem) - 1; - } - const uint64_t owd = outside_mask & added_tips[w]; + // T-386: the outside group is cd.split_zeros, not ~split_tips. A `?`-coded + // tip is in neither, so it never pulls this LCA about — which is what lets + // it be placed on either side of the constraint, as documented. + const uint64_t owd = split_out[w] & added_tips[w]; needed_out[w] = owd; if (owd) { n_out += popcount64(owd); @@ -426,6 +425,8 @@ static void wagner_map_constraint_nodes( for (int s = 0; s < cd.n_splits; ++s) { const uint64_t* split = &cd.split_tips[static_cast(s) * nw]; + const uint64_t* split_out = + &cd.split_zeros[static_cast(s) * nw]; // Once the inside LCA has reached the root it can never come back down — // grafting a leaf preserves ancestor relations among existing nodes, so the @@ -441,9 +442,13 @@ static void wagner_map_constraint_nodes( // keep ConstraintData's T-384 flag at its "names the split itself" // default so a value left over from an earlier map_constraint_nodes() // cannot reach regraft_violates_constraint() before the next full remap. + // constraint_node_hi is pinned to constraint_node for the same reason: + // Wagner's LCA mapping has no displaying-chain, so the tight end is the + // only anchor it can honestly offer (T-386). cd.constraint_complement[s] = 0; + cd.constraint_node_hi[s] = cd.constraint_node[s]; scratch.outside_node[s] = wagner_map_complement( - tree, n_tip, nw, node_tips, split, added_tips, scratch); + tree, n_tip, nw, node_tips, split_out, added_tips, scratch); continue; } @@ -466,6 +471,7 @@ static void wagner_map_constraint_nodes( tree, nw, node_tips, needed, n_needed, lone_needed); cd.constraint_node[s] = inside_node; cd.constraint_complement[s] = 0; // see the latched branch above + cd.constraint_node_hi[s] = inside_node; // A split is an *unrooted* bipartition, but a clade is a rooted subtree, so // "inside is monophyletic" is only one of the two ways this tree can display @@ -485,7 +491,7 @@ static void wagner_map_constraint_nodes( if (inside_node == n_tip) { scratch.use_complement[s] = 1; scratch.outside_node[s] = wagner_map_complement( - tree, n_tip, nw, node_tips, split, added_tips, scratch); + tree, n_tip, nw, node_tips, split_out, added_tips, scratch); } } } @@ -534,8 +540,18 @@ static void wagner_collect_active_splits( for (int s = 0; s < cd.n_splits; ++s) { const uint64_t* split = &cd.split_tips[static_cast(s) * cd.n_words]; + const uint64_t* split_out = + &cd.split_zeros[static_cast(s) * cd.n_words]; + + const bool tip_inside = (split[tw] >> tb) & 1; + const bool tip_outside = (split_out[tw] >> tb) & 1; - bool tip_inside = (split[tw] >> tb) & 1; + // T-386: a tip in neither group is FREE — the constraint says nothing + // about which side of the separating edge it belongs on, so no edge is + // barred to it. Reading "outside" as ~split_tips, as this did before free + // tips were represented, forced every `?`-coded taxon out of the + // constrained group and could leave the filter with no legal edge at all. + if (!tip_inside && !tip_outside) continue; // Split constrains placement when the opposite side of the new tip // has at least one previously-added tip. An inside tip is only @@ -544,15 +560,10 @@ static void wagner_collect_active_splits( for (int w = 0; w < cd.n_words; ++w) { uint64_t prev = added_tips[w]; if (prev & split[w]) has_prev_inside = true; - uint64_t outside_mask = ~split[w]; - if (w == cd.n_words - 1) { - int rem = tree.n_tip % 64; - if (rem > 0) outside_mask &= (1ULL << rem) - 1; - } - if (prev & outside_mask) has_prev_outside = true; + if (prev & split_out[w]) has_prev_outside = true; } if (tip_inside && !has_prev_outside) continue; - if (!tip_inside && !has_prev_inside) continue; + if (tip_outside && !has_prev_inside) continue; // The constraint is active. constraint_node[s] is the LCA of // added inside tips (set by wagner_map_constraint_nodes). @@ -586,12 +597,20 @@ static void wagner_collect_active_splits( // Does the finished tree display every constraint split? // -// A bipartition is displayed iff some edge separates its two sides, i.e. iff -// some node's subtree tip set equals one side exactly. Only non-root nodes are -// candidates: the root subtends every tip, and its two children already cover -// the single edge the degree-two root sits on. Tips are included so trivial -// (single-taxon) splits are recognised. Orientation-agnostic by construction, -// so it stays correct however the tree happens to be rooted. +// A split is displayed iff some edge separates its two groups, i.e. iff some +// node's subtree tip set covers one group and holds none of the other — the +// same free-taxa reading map_constraint_nodes() applies (T-386), spelled out +// again here because Wagner runs before any of that machinery is valid. It +// MUST stay in step with node_displays_split() in ts_constraint.cpp: this is +// the gate that decides whether the caller reshuffles and rebuilds, and the +// stricter of the two entry points would reject trees the other then searches +// happily (or, worse, accept ones it will not move from). +// +// Only non-root nodes are candidates: the root subtends every tip, and its two +// children already cover the single edge the degree-two root sits on. Tips are +// included so trivial (single-taxon) groups are recognised. +// Orientation-agnostic by construction, so it stays correct however the tree +// happens to be rooted. static bool wagner_tree_displays_constraint(const TreeState& tree, const ConstraintData& cd) { const int n_tip = tree.n_tip; @@ -612,22 +631,21 @@ static bool wagner_tree_displays_constraint(const TreeState& tree, for (int s = 0; s < cd.n_splits; ++s) { const uint64_t* split = &cd.split_tips[static_cast(s) * nw]; + const uint64_t* split_out = &cd.split_zeros[static_cast(s) * nw]; bool found = false; for (int node = 0; node < tree.n_node && !found; ++node) { if (node == n_tip) continue; // root subtends everything const uint64_t* nd = &node_tips[static_cast(node) * nw]; - bool eq = true, eqCompl = true; + // Either group may be the clade side; free tips (in neither mask) are + // not looked at. No width mask is needed — both masks carry zeros above + // tip n_tip - 1, so padding bits can never make a test fail. + bool holds = true, holdsCompl = true; for (int w = 0; w < nw; ++w) { - uint64_t tip_mask = ~0ULL; - if (w == nw - 1) { - int rem = n_tip % 64; - if (rem > 0) tip_mask = (1ULL << rem) - 1; - } - if (nd[w] != (split[w] & tip_mask)) eq = false; - if (nd[w] != (~split[w] & tip_mask)) eqCompl = false; - if (!eq && !eqCompl) break; + if ((split[w] & ~nd[w]) || (split_out[w] & nd[w])) holds = false; + if ((split_out[w] & ~nd[w]) || (split[w] & nd[w])) holdsCompl = false; + if (!holds && !holdsCompl) break; } - if (eq || eqCompl) found = true; + if (holds || holdsCompl) found = true; } if (!found) return false; } diff --git a/tests/testthat/test-ts-constraint-free-taxa.R b/tests/testthat/test-ts-constraint-free-taxa.R new file mode 100644 index 000000000..8a33068dd --- /dev/null +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -0,0 +1,192 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +## T-386 / agent-issues/TreeSearch#54: the constraint the search enforces must +## be the one `?MaximizeParsimony`'s `constraint` argument documents. +## +## The documented contract is the phyDat reading: a returned tree is compliant +## when some edge separates the taxa coded `1` from those coded `0`, with +## `?`-coded taxa free to fall on either side. The locked-node machinery used +## to enforce a strictly stronger one — some node's tip set had to EQUAL the 1 +## group (or its exact complement), free taxa excluded. +## +## Strict implies loose, so no wrong answer was ever returned. What broke was +## movement: a tree that satisfies the documented contract without making +## either group an exact clade mapped to no node at all, which +## regraft_violates_constraint() reads as "the tree already violates" and +## answers by rejecting EVERY rearrangement. The replicate froze on its start. +## +## So these tests assert two things together, and neither alone is the fix: +## that the search MOVES from such a start, and that what it returns is still +## compliant. + +library("TreeTools") + +# 8 taxa; three characters agree on {a,b,e,f} | {c,d,g,h} and one cuts across +# it, so the start tree below is a local optimum only for a search that cannot +# move. +freeTaxaData <- function() { + taxa <- letters[1:8] + phangorn::phyDat( + matrix(c("0", "0", "1", "1", "0", "0", "1", "1", + "0", "0", "1", "1", "0", "0", "1", "1", + "1", "1", "0", "0", "1", "1", "0", "0", + "0", "1", "0", "1", "0", "1", "0", "1", + "0", "1", "0", "1", "0", "1", "0", "1"), + nrow = 8, dimnames = list(taxa, NULL)), + type = "USER", levels = c("0", "1") + ) +} + +# c(a = 1, b = 1, c = 0, d = 0, e:h = "?") +freeTaxaConstraint <- function() { + phangorn::phyDat( + matrix(c("1", "1", "0", "0", "?", "?", "?", "?"), + nrow = 8, dimnames = list(letters[1:8], NULL)), + type = "USER", levels = c("0", "1") + ) +} + +# `{a,e,b,f}` | `{c,d,g,h}` separates {a,b} from {c,d}, so this satisfies the +# documented constraint — but {a,b} is not a clade, and neither is {c,d}. +freeTaxaStart <- function() { + ape::read.tree(text = "(((a,e),(b,f)),(c,(d,(g,h))));") +} + +# Does `tree` display a split with every tip of `inGroup` on one side and every +# tip of `outGroup` on the other? Tips in neither group are ignored — this is +# the documented contract, spelled out independently of the engine. +SeparatesGroups <- function(tree, labels, inGroup, outGroup) { + splits <- as.logical(as.Splits(tree, tipLabels = labels)) + if (!is.matrix(splits)) splits <- matrix(splits, nrow = 1) + isIn <- labels %in% inGroup + isOut <- labels %in% outGroup + any(apply(splits, 1, function(row) { + (all(row[isIn]) && !any(row[isOut])) || + (!any(row[isIn]) && all(row[isOut])) + })) +} + +# The engine hands back bare edge matrices for rooted binary trees. +EdgeSeparatesGroups <- function(edge, labels, inGroup, outGroup) { + tree <- structure( + list(edge = edge, Nnode = length(labels) - 1L, tip.label = labels), + class = "phylo") + SeparatesGroups(tree, labels, inGroup, outGroup) +} + +# TBR only: every phase that could rescue a frozen replicate by starting +# somewhere else is switched off, so the reported score is what rearranging the +# supplied start achieved. Mirrors tbrOnlyRun() in test-ts-constraint-rooting.R. +tbrOnly <- function(ds, startEdge, splitMatrix) { + TreeSearch:::ts_driven_search( + ds$contrast, ds$tip_data, ds$weight, ds$levels, + maxReplicates = 1L, targetHits = 99L, tbrMaxHits = 1L, + ratchetCycles = 0L, driftCycles = 0L, nniPerturbCycles = 0L, + xssRounds = 0L, rssRounds = 0L, cssRounds = 0L, + pruneReinsertCycles = 0L, fuseInterval = 0L, + outerCycles = 1L, maxOuterResets = 0L, + nniFirst = FALSE, sprFirst = FALSE, + poolMaxSize = 100L, poolSuboptimal = 0, maxSeconds = 0, verbosity = 0L, + nThreads = 1L, startEdge = startEdge, consSplitMatrix = splitMatrix + ) +} + +test_that("free `?` taxa do not freeze a compliant start tree", { + dataset <- freeTaxaData() + labels <- names(dataset) + ds <- make_ts_data(dataset) + start <- Preorder(RenumberTips(freeTaxaStart(), labels)) + startScore <- TreeLength(start, dataset) + + # Premise: the start satisfies the documented constraint but neither group is + # a clade, so the pre-T-386 exact-clade test could not map it. + expect_true(SeparatesGroups(start, labels, c("a", "b"), c("c", "d"))) + expect_false(SeparatesGroups(start, labels, c("a", "b"), + setdiff(labels, c("a", "b")))) + + # Take the split matrix from .PrepareConstraint rather than writing it out, + # so this exercises the same R -> C++ contract the user's `constraint =` + # phyDat travels along: pre-T-386 it coded the free taxa 0 (making {a,b} an + # exact clade), now it codes them NA. + free <- TreeSearch:::.PrepareConstraint( + freeTaxaConstraint(), dataset)[["consSplitMatrix"]] + expect_equal(as.vector(free), c(1L, 1L, 0L, 0L, NA, NA, NA, NA)) + + set.seed(386) + result <- tbrOnly(ds, start[["edge"]], free) + + # The defect: every regraft was rejected and the start came back unimproved. + expect_lt(result$best_score, startScore) + + # ... and what comes back still honours the documented constraint. + expect_true(all(vapply(result$trees, EdgeSeparatesGroups, logical(1), + labels = labels, inGroup = c("a", "b"), + outGroup = c("c", "d")))) +}) + +test_that("a 0/1 constraint matrix still enforces the exact clade", { + # Guard against over-loosening: with no free tips the two groups are + # complements, and every check must collapse back to the exact-clade test. + dataset <- freeTaxaData() + labels <- names(dataset) + ds <- make_ts_data(dataset) + start <- Preorder(RenumberTips( + ape::read.tree(text = "(((a,b),(e,f)),(c,(d,(g,h))));"), labels)) + + strict <- matrix(c(1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), nrow = 1) + set.seed(386) + result <- tbrOnly(ds, start[["edge"]], strict) + + expect_true(all(vapply(result$trees, EdgeSeparatesGroups, logical(1), + labels = labels, inGroup = c("a", "b"), + outGroup = setdiff(labels, c("a", "b"))))) +}) + +test_that("MaximizeParsimony honours and searches under a `?` constraint", { + # End to end. Unlike the TBR-only test above this cannot isolate the freeze + # — ratchet, drift and nni-perturb all get a replicate moving again by other + # means — so what it adds is that the whole pipeline still returns compliant + # trees once every entry point reads the constraint the same way. + dataset <- freeTaxaData() + labels <- names(dataset) + start <- freeTaxaStart() + + set.seed(386) + result <- suppressWarnings( + MaximizeParsimony(dataset, tree = start, + constraint = freeTaxaConstraint(), + maxReplicates = 4L, verbosity = 0L)) + expect_lt(attr(result, "score"), TreeLength(start, dataset)) + + expect_true(all(vapply(result, SeparatesGroups, logical(1), + labels = labels, inGroup = c("a", "b"), + outGroup = c("c", "d")))) +}) + +test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { + dataset <- freeTaxaData() + + consArgs <- TreeSearch:::.PrepareConstraint(freeTaxaConstraint(), dataset) + expect_equal(nrow(consArgs[["consSplitMatrix"]]), 1L) + expect_equal(as.vector(consArgs[["consSplitMatrix"]]), + c(1L, 1L, 0L, 0L, NA, NA, NA, NA)) + + # A taxon on the tree but absent from the constraint is free too. + partial <- phangorn::phyDat( + matrix(c("1", "1", "0", "0"), nrow = 4, + dimnames = list(letters[1:4], NULL)), + type = "USER", levels = c("0", "1")) + consArgs <- TreeSearch:::.PrepareConstraint(partial, dataset) + expect_equal(as.vector(consArgs[["consSplitMatrix"]]), + c(1L, 1L, 0L, 0L, NA, NA, NA, NA)) + + # A character with no `0` taxa constrains nothing under the documented + # contract — there is no group for the `1` taxa to be separated FROM — so it + # must not be enforced as a clade. + vacuous <- phangorn::phyDat( + matrix(c("1", "1", "?", "?", "?", "?", "?", "?"), nrow = 8, + dimnames = list(letters[1:8], NULL)), + type = "USER", levels = c("0", "1")) + expect_equal(TreeSearch:::.PrepareConstraint(vacuous, dataset), list()) +}) diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index a87d2b292..b561127a6 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -217,6 +217,36 @@ Per-strategy attempt and success counts are returned in the `strategy_diagnostics` attribute of the search result for post-hoc inspection. +### What a topological constraint requires + +A tree satisfies a constraint character when some edge separates the taxa +coded `1` from those coded `0`. +Taxa coded `?`, and taxa that the constraint does not mention at all, are +*free*: the constraint says nothing about which side of that edge they belong +on, and a tree is compliant however they fall. +This is what `constraint` promises the user, and every part of the search reads +it the same way -- the locked-node filter that screens individual +rearrangements, the check that gates a finished replicate on its way into the +pool, the constrained Wagner build, and the minimal-SPR repair that fixes a +violating start. + +Reading a constraint more strictly -- as "the `1` group is a clade exactly", +excluding the free taxa -- never returns a wrong answer, because every +strictly-compliant tree also satisfies the user's constraint. +It does, however, cost search: a tree that satisfies the constraint the user +wrote, without making either group an exact clade, matches no node under the +stricter reading, and an unmapped constraint split makes every candidate +regraft illegal, so the replicate freezes on the tree it began with. +Enforcing one reading everywhere is therefore a search-quality matter rather +than a correctness one. + +The free-taxa reading also determines which rearrangements are legal, not just +which trees are. +A subtree made only of free taxa may be regrafted anywhere, since moving it +cannot disturb the separating edge; and a subtree carrying `1`-group taxa may +be regrafted anywhere within the largest clade that still covers the `1` group +and excludes the `0` group, not merely within the smallest one. + ## The driven search pipeline From bf11273dcd359d27e5c20509d98f1949680313e8 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 04:32:12 +0100 Subject: [PATCH 2/4] test: pin the collapse pass's constraint protection `ts_collapse_pool()` identified the branch realising a constraint by matching a node's tip set to the `1` group exactly, so with free taxa it protected nothing: the separating edge was contracted and the RETURNED tree broke the documented constraint, even though every tree the search visited satisfied it. This is the one place the strict reading did produce a wrong answer, so it gets a deterministic test of its own rather than riding on the search tests. Co-Authored-By: Claude Opus 5 --- tests/testthat/test-ts-constraint-free-taxa.R | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/testthat/test-ts-constraint-free-taxa.R b/tests/testthat/test-ts-constraint-free-taxa.R index 8a33068dd..6aed14d14 100644 --- a/tests/testthat/test-ts-constraint-free-taxa.R +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -164,6 +164,47 @@ test_that("MaximizeParsimony honours and searches under a `?` constraint", { outGroup = c("c", "d")))) }) +test_that("the collapse pass keeps the enforced grouping visible", { + # The one place the strict reading did return a wrong answer. A constraint + # is external evidence for a grouping, so ts_collapse_pool() protects the + # branch that realises it from contraction — but it identified that branch by + # matching a node's tip set to the `1` group EXACTLY. With free taxa the + # realising node is generally not that set, so nothing was protected and the + # separating edge was contracted away: the returned tree broke the documented + # constraint even though every tree the search visited satisfied it. + labels <- letters[1:6] + # Two characters support (a,e), two support (b,f); none supports (c,d), so + # the branch that separates {a,b} from {c,d} is unsupported and collapses + # unless it is protected. + charDat <- StringToPhyDat( + c("100010", "100010", "010001", "010001", "000000"), labels) + at <- attributes(charDat) + 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 + ) + # {a,e,b,f} | {c,d} separates {a,b} from {c,d}; neither group is a clade. + tree <- Preorder(RenumberTips( + ape::read.tree(text = "(((a,e),(b,f)),(c,d));"), labels)) + cons <- phangorn::phyDat( + matrix(c("1", "1", "0", "0", "?", "?"), nrow = 6, + dimnames = list(labels, NULL)), + type = "USER", levels = c("0", "1")) + + collapsed <- TreeSearch:::ts_collapse_pool( + list(tree[["edge"]]), at$contrast, + matrix(unlist(charDat, use.names = FALSE), nrow = 6, byrow = TRUE), + at$weight, at$levels, scoringConfig, NULL, NULL, + TreeSearch:::.PrepareConstraint(cons, charDat)[["consSplitMatrix"]]) + + out <- structure( + list(edge = collapsed$trees[[1]], tip.label = labels, + Nnode = max(collapsed$trees[[1]]) - 6L), + class = "phylo") + expect_true(SeparatesGroups(Renumber(out), labels, c("a", "b"), c("c", "d"))) +}) + test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { dataset <- freeTaxaData() From 57ee3802b4c79bb2628d58777b38c56ae2362454 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 08:47:31 +0100 Subject: [PATCH 3/4] fix: address review of the loosened constraint contract Three independent reviews of the first commit. Substantive changes: * A clip whose OWN tip set displays the split is now UNCONSTRAINED (`classify_clip_constraints`). It carries the constraint with it: the node at whatever attachment point it lands on has exactly the clip's tip set, so the split stays displayed, and TBR's rerooting cannot change that. Without this, a clip containing the whole displaying chain left no surviving `below` that could be the anchor's descendant, and every regraft of a subtree that was in fact free to go anywhere was rejected. * `node_displays_split()` moves to `ts_constraint.h`, so the search's mapping, the Wagner build's check and the collapse pass's branch protection call ONE predicate instead of three lookalikes. The comment demanding they stay in step is now enforced by construction. * `wagner_tree_displays_constraint()` no longer skips the root, which made it stricter than `find_displaying_chain()` (which scans `tree.postorder`, and that includes the root): a split with an empty apart-group mapped there and nowhere else, so the two entry points disagreed about a tree every constraint is satisfied by. Including the root can never accept a violation. * The inert-character filter is symmetric. Since `build_constraint()` now swaps the two groups to canonicalise, they are interchangeable, and a test on the `1` group alone was incoherent: `c(a = 1, b = 1, c = 0)` and `c(a = 0, b = 0, c = 1)` state the same constraint and were treated differently. A group of fewer than two taxa is separated from the rest by every tree, so such a character is ignored -- with a warning, because coding only `1` and `?` almost always means "group these taxa", which is not what it says. Plus: NEWS entries for both behaviour changes; the `@param constraint` doc states the inert-character rule; `.AGENTS/memory/architecture.md` records the 1/0/NA encoding and the shared predicate; `wagner_map_complement()` loses its dead `n_tip` parameter; the TBR-only test harness moves to `helper-ts.R` instead of being copied; `wagner_tree()`'s comment about `has_posthoc` is corrected (AdditionTree DOES build the posthoc DataSet -- it reaches `ts_wagner_tree`, which bypasses the reshuffle loop, which is the real reason); and `random_constrained_tree()`'s narrower sampling is documented rather than changed, since widening it is a search-quality change to measure on its own. Co-Authored-By: Claude Opus 5 --- .AGENTS/memory/architecture.md | 19 +++- NEWS.md | 23 +++++ R/MaximizeParsimony.R | 32 +++++-- man/AdditionTree.Rd | 4 + man/MaximizeParsimony.Rd | 4 + man/Resample.Rd | 4 + man/SuccessiveApproximations.Rd | 4 + src/ts_constraint.cpp | 54 ++++++----- src/ts_constraint.h | 38 ++++++-- src/ts_rcpp.cpp | 17 ++-- src/ts_wagner.cpp | 79 ++++++++------- tests/testthat/helper-ts.R | 21 ++++ tests/testthat/test-ts-constraint-free-taxa.R | 95 +++++++++++++------ tests/testthat/test-ts-constraint-rooting.R | 19 +--- vignettes/search-algorithm.Rmd | 2 +- 15 files changed, 282 insertions(+), 133 deletions(-) diff --git a/.AGENTS/memory/architecture.md b/.AGENTS/memory/architecture.md index 702ad3492..650abcb97 100644 --- a/.AGENTS/memory/architecture.md +++ b/.AGENTS/memory/architecture.md @@ -94,8 +94,25 @@ Profile mode sets `ds.concavity = 1.0` (finite sentinel) so existing ## Constraint enforcement +- A constraint split names **two disjoint groups** plus a FREE remainder. A tree + satisfies it iff some edge separates the groups; free tips may fall either + side. This is what `?MaximizeParsimony`'s `constraint` documents and, since + agent-issues/TreeSearch#54, what every entry point enforces. `ts::node_displays_split()` + (`ts_constraint.h`) is THE shared predicate — `map_constraint_nodes()`, + `wagner_tree_displays_constraint()` and `ts_collapse_pool()` all call it. + Reintroducing an exact-clade test at any one of them freezes replicates. - `build_constraint()` reads R split matrix with **column-major** indexing: - `split_matrix[s + n_splits * t]`. + `split_matrix[s + n_splits * t]`. Values: `1` = together-group, `0` = + apart-group, **anything else (`NA_INTEGER`) = free**. A hand-built 0/1 matrix + therefore means "no free tips" and reduces to the exact-clade behaviour, which + is what `build_constraint_from_bitsets()` (consensus constraints) relies on. +- `ConstraintData` carries `split_zeros` (the apart-group) alongside + `split_tips`, and both ends of the displaying-node chain: + `constraint_node` (tightest, used for "must land outside") and + `constraint_node_hi` (highest, "must land inside"). Any writer of one must + write the other — `ts_wagner.cpp` pins hi to the tight anchor. +- `.PrepareConstraint()` drops (and warns about) a character with no `0` taxa: + vacuous under the documented contract. - Wagner uses LCA-based constraint mapping (`wagner_map_constraint_nodes`) since splits aren't fully present during incremental construction. - Wagner has a posthoc retry loop (up to 100 random addition orders) as a diff --git a/NEWS.md b/NEWS.md index e104fe8c3..6a209c179 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,28 @@ # To integrate into 2.0.0 notes +- `constraint` now enforces exactly what it documents: a returned tree is + compatible with a constraint character when some edge separates the taxa + coded `1` from those coded `0`, with `?`-coded and unmentioned taxa free to + fall on either side. The enforcement machinery previously required the `1` + group to be a clade *exactly*, free taxa excluded. That is strictly + stronger, so the search never accepted a tree that broke the documented + constraint; but a start tree that satisfied the documented constraint without + making either group an exact clade matched no node, every rearrangement was + rejected, and the replicate returned its start unimproved. Constrained + searches with `?`-coded taxa therefore reach better scores. + The collapse pass is fixed with it: it identified the branch realising a + constraint by exact match too, so with free taxa it protected nothing and the + separating edge could be contracted away -- the one route by which a + *returned* tree could break the constraint. +- A constraint character whose `1` or `0` group holds fewer than two taxa now + warns and is ignored, rather than being enforced as a clade. Every tree + separates such a group from the rest, so the character constrains nothing + under the documented reading. The test is symmetric in the two groups, which + the old one was not: `c(a = 1, b = 1, c = 0)` and `c(a = 0, b = 0, c = 1)` + state the same constraint and are now treated the same way. Code the taxa + that must fall outside a group as `0`, rather than leaving them `?`, to keep + it enforced. + - `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/MaximizeParsimony.R b/R/MaximizeParsimony.R index d0f8842dd..637f7337b 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -145,14 +145,30 @@ } } - # A character only constrains anything when both groups are occupied: with no - # "0" tips there is no edge for the "1" tips to be separated *from*, so the - # documented contract ("some edge separates the 1 taxa from the 0 taxa") is - # vacuously true and enforcing the group as a clade would be a restriction the - # user never asked for. + # Every tree separates a group of fewer than two taxa from anything: the edge + # above a lone tip already does it, and an empty group needs no edge at all. + # Such a character constrains nothing under the documented contract, so + # enforcing its group as a clade would restrict the search for a guarantee it + # already has -- which is the over-strict reading agent-issues/TreeSearch#54 + # is about. The test is symmetric in the two groups because they are + # interchangeable: which one a user calls "1" is arbitrary, and + # build_constraint() swaps them freely to canonicalise. + # + # Warn rather than drop silently: a character coding only "1" and "?" almost + # certainly means "group these taxa", which is not what it says. nOne <- rowSums(consSplits) nZero <- rowSums(consZero) - keep <- nOne >= 1 & nZero >= 1 & nOne < length(constraint) - 1 + inert <- nOne < 2 | nZero < 2 + if (any(inert)) { + warning("Constraint character", if (sum(inert) > 1) "s" else "", " ", + paste(which(inert), collapse = ", "), + if (sum(inert) > 1) " constrain" else " constrains", + " nothing, and", if (sum(inert) > 1) " are" else " is", + " ignored: every tree separates a group of fewer than two taxa ", + "from the rest. Taxa coded `?` join neither group; code those ", + "that must fall outside the group as `0`.", call. = FALSE) + } + keep <- !inert consSplits <- consSplits[keep, , drop = FALSE] consZero <- consZero[keep, , drop = FALSE] if (nrow(consSplits) == 0L) return(list()) @@ -711,6 +727,10 @@ #' separates the taxa coded `1` from those coded `0`. Taxa coded `?`, and taxa #' that `constraint` does not mention, are unconstrained: they may fall on #' either side of that edge, and are not required to join either group. +#' A character whose `1` or `0` group contains fewer than two taxa therefore +#' constrains nothing -- every tree separates such a group from the rest -- and +#' is ignored with a warning. To group taxa, code the taxa they must be +#' separated from as `0` rather than leaving them `?`. #' Constraint searches are supported natively: all tree rearrangements #' are filtered to respect the constraint topology. #' @param effort Integer: how much search effort to spend, **relative to the diff --git a/man/AdditionTree.Rd b/man/AdditionTree.Rd index 380bea854..ef6d60801 100644 --- a/man/AdditionTree.Rd +++ b/man/AdditionTree.Rd @@ -39,6 +39,10 @@ A returned tree is compatible with a constraint character when some edge separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa that \code{constraint} does not mention, are unconstrained: they may fall on either side of that edge, and are not required to join either group. +A character whose \code{1} or \code{0} group contains fewer than two taxa therefore +constrains nothing -- every tree separates such a group from the rest -- and +is ignored with a warning. To group taxa, code the taxa they must be +separated from as \code{0} rather than leaving them \verb{?}. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 28c82714b..91a8d126c 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -141,6 +141,10 @@ A returned tree is compatible with a constraint character when some edge separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa that \code{constraint} does not mention, are unconstrained: they may fall on either side of that edge, and are not required to join either group. +A character whose \code{1} or \code{0} group contains fewer than two taxa therefore +constrains nothing -- every tree separates such a group from the rest -- and +is ignored with a warning. To group taxa, code the taxa they must be +separated from as \code{0} rather than leaving them \verb{?}. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/Resample.Rd b/man/Resample.Rd index 329b065bf..cc4bed912 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -75,6 +75,10 @@ A returned tree is compatible with a constraint character when some edge separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa that \code{constraint} does not mention, are unconstrained: they may fall on either side of that edge, and are not required to join either group. +A character whose \code{1} or \code{0} group contains fewer than two taxa therefore +constrains nothing -- every tree separates such a group from the rest -- and +is ignored with a warning. To group taxa, code the taxa they must be +separated from as \code{0} rather than leaving them \verb{?}. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/man/SuccessiveApproximations.Rd b/man/SuccessiveApproximations.Rd index 864bd5481..35ecb7f5b 100644 --- a/man/SuccessiveApproximations.Rd +++ b/man/SuccessiveApproximations.Rd @@ -87,6 +87,10 @@ A returned tree is compatible with a constraint character when some edge separates the taxa coded \code{1} from those coded \code{0}. Taxa coded \verb{?}, and taxa that \code{constraint} does not mention, are unconstrained: they may fall on either side of that edge, and are not required to join either group. +A character whose \code{1} or \code{0} group contains fewer than two taxa therefore +constrains nothing -- every tree separates such a group from the rest -- and +is ignored with a warning. To group taxa, code the taxa they must be +separated from as \code{0} rather than leaving them \verb{?}. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} diff --git a/src/ts_constraint.cpp b/src/ts_constraint.cpp index 908f7d0f6..881e7e492 100644 --- a/src/ts_constraint.cpp +++ b/src/ts_constraint.cpp @@ -34,7 +34,7 @@ ConstraintData build_constraint( // Pack split_matrix rows into a pair of bitmasks. // split_matrix is column-major (from R): element [s, t] is at // index s + n_splits * t. 1 -> "together" group, 0 -> "apart" group, - // anything else (NA_INTEGER) -> free, in neither mask (T-386). + // anything else (NA_INTEGER) -> free, in neither mask (#54). for (int s = 0; s < n_splits; ++s) { uint64_t* ones = &cd.split_tips[static_cast(s) * cd.n_words]; uint64_t* zeros = &cd.split_zeros[static_cast(s) * cd.n_words]; @@ -87,7 +87,7 @@ ConstraintData build_constraint_from_bitsets( cd.split_tips.assign(split_bits, split_bits + total); // These splits come from pool bipartitions, which partition every tip: there // are no free tips, so the "apart" group is exactly the complement and the - // free-taxa machinery collapses back to the exact-clade test (T-386). + // free-taxa machinery collapses back to the exact-clade test (#54). cd.split_zeros.assign(total, 0ULL); { const int rem = n_tips % 64; @@ -173,24 +173,9 @@ std::vector compute_node_tips(const TreeState& tree, int n_words) // Map constraint nodes: find which internal node holds each split // ========================================================================= -// Does the edge above `node` separate `together` from `apart`? It does when -// the node's descendant tip set covers every tip of `together` and holds none -// of `apart`; the tips in neither group are free and are not looked at (T-386). -// -// With `apart` the exact complement of `together` — a constraint with no free -// tips, and every split built by build_constraint_from_bitsets() — the two -// conditions together force set equality, which is the exact-clade test this -// replaced. -static inline bool node_displays_split( - const uint64_t* nd, const uint64_t* together, const uint64_t* apart, - int n_words) -{ - for (int w = 0; w < n_words; ++w) { - if ((together[w] & ~nd[w]) != 0ULL) return false; // a required tip missing - if ((apart[w] & nd[w]) != 0ULL) return false; // an excluded tip present - } - return true; -} +// node_displays_split() — the shared "does this node display the split" +// predicate — lives in ts_constraint.h, so the Wagner build and the collapse +// pass answer the question with the same code rather than a lookalike. // Tightest and highest node displaying `together` | `apart`, or {-1, -1}. // @@ -277,7 +262,7 @@ void map_constraint_nodes(const TreeState& tree, ConstraintData& cd) // that every tree which mapped successfully before maps to exactly the same // node now. // - // T-386: "is a clade" is the free-taxa reading, not set equality -- see + // #54: "is a clade" is the free-taxa reading, not set equality -- see // node_displays_split(). The chain of displaying nodes is recorded at both // ends, because a regraft that must land INSIDE the constrained group may use // the whole chain while one that must land outside may not; see @@ -396,7 +381,7 @@ void classify_clip_constraints(const TreeState& tree, int clip_node, // "outside", a tip of the group that must stay apart from it. A clip made // only of FREE tips is in neither, and lands here as UNCONSTRAINED — the // whole point of the free-taxa reading, and what lets such a clip be - // regrafted anywhere without breaking the separating edge (T-386). Reading + // regrafted anywhere without breaking the separating edge (#54). Reading // "outside" as ~split, which is what this did before free tips existed, // pinned every free tip to the far side of the constraint. for (int s = 0; s < cd.n_splits; ++s) { @@ -405,6 +390,23 @@ void classify_clip_constraints(const TreeState& tree, int clip_node, const uint64_t* zeros = &cd.split_zeros[static_cast(s) * cd.n_words]; + // The clip carries the constraint with it: its own tip set covers one + // group and holds none of the other. Wherever it is regrafted, the node + // at the attachment point has exactly the clip's tip set, so the split + // stays displayed — and TBR's rerooting of the clip cannot change that, + // since the set is the same however the subtree hangs. Testing this + // first is what unpins a clip that contains the whole displaying chain: + // the anchor is then inside the clipped subtree, no surviving `below` can + // be its descendant, and the MUST_INSIDE test below would reject every + // regraft of a subtree that is in fact free to go anywhere. + if (node_displays_split(cd.clip_tip_mask.data(), ones, zeros, + cd.n_words) || + node_displays_split(cd.clip_tip_mask.data(), zeros, ones, + cd.n_words)) { + cd.clip_zones[s] = ClipZone::UNCONSTRAINED; + continue; + } + bool any_inside = false; bool any_outside = false; for (int w = 0; w < cd.n_words; ++w) { @@ -486,8 +488,8 @@ bool regraft_violates_constraint(int below, ? ClipZone::MUST_INSIDE : ClipZone::MUST_OUTSIDE; // The two ends of the displaying chain answer two different questions, and - // each wants the end that permits most (T-386; with no free tips the chain - // is one node long and both reduce to the pre-T-386 test): + // each wants the end that permits most (#54; with no free tips the chain + // is one node long and both reduce to the pre-#54 test): // // * a clip that must land INSIDE carries tips of the together-group but // none of the apart-group, so anywhere within the HIGHEST displaying @@ -795,7 +797,7 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, // already displayed every constraint, spending up to n_tip / 4 + 2 arbitrary // SPR moves on it. That mattered most at ts_nni_perturb.cpp's unconditional // impose_constraint() call, which runs after every perturbation cycle. - // "Is a clade" is the free-taxa reading (T-386), so a tree that satisfies + // "Is a clade" is the free-taxa reading (#54), so a tree that satisfies // what `constraint` documents is likewise left alone. The repair below aims // at making the canonical side a clade, which displays the split either way. std::vector violated; @@ -833,7 +835,7 @@ static int impose_one_pass(TreeState& tree, ConstraintData& cd, // Tips that keep node `nd` from displaying the split: those of the // together-group it is missing, plus those of the apart-group it holds. - // Free tips appear in neither, so the repair never moves one (T-386) — they + // Free tips appear in neither, so the repair never moves one (#54) — they // may sit on whichever side they already do. auto repair_cost = [&](const uint64_t* nd, const uint64_t* ones, const uint64_t* zeros) { diff --git a/src/ts_constraint.h b/src/ts_constraint.h index ee4a16061..e929489e6 100644 --- a/src/ts_constraint.h +++ b/src/ts_constraint.h @@ -8,11 +8,12 @@ // remaining tips are FREE: coded `?`, or absent from the constraint phyDat // altogether. A tree satisfies the split iff some edge separates the 1 group // from the 0 group, free tips falling on either side. That is the contract -// `?MaximizeParsimony`'s `constraint` argument documents, and since T-386 it is -// the one enforced here: a node DISPLAYS the split when its descendant tip set -// is a superset of one group and disjoint from the other. +// `?MaximizeParsimony`'s `constraint` argument documents, and since +// agent-issues/TreeSearch#54 it is the one enforced here: a node DISPLAYS the +// split when its descendant tip set is a superset of one group and disjoint +// from the other. // -// Requiring the tip set to EQUAL a group (the pre-T-386 test) is strictly +// Requiring the tip set to EQUAL a group (the pre-#54 test) is strictly // stronger. It never returned a wrong answer, but a start tree that satisfied // the documented contract without making either group an exact clade mapped to // no node at all, which regraft_violates_constraint() reads as "the tree @@ -53,9 +54,9 @@ struct ConstraintData { // The tips that must end up on the OTHER side of that edge, same layout. // Disjoint from split_tips; the two need NOT be complements — a tip in - // neither mask is free to fall on either side (T-386). When the caller + // neither mask is free to fall on either side (#54). When the caller // supplies no free tips this is exactly ~split_tips, and every check below - // reduces to the pre-T-386 exact-clade test. + // reduces to the pre-#54 exact-clade test. std::vector split_zeros; // Current mapping: constraint_node[i] = the TIGHTEST node (tip or internal) @@ -109,7 +110,7 @@ struct ConstraintData { // anything else (NA_INTEGER, as .PrepareConstraint() writes for a `?`-coded // or unconstrained taxon) — tip t is FREE, and may fall on either side. // A pure 0/1 matrix therefore means "no free tips", i.e. the exact-clade -// reading that predates T-386; callers that build one by hand keep it. +// reading that predates #54; callers that build one by hand keep it. // The two groups are swapped where needed so that tip 0 is never in // split_tips ("outside" the canonical side) — the same invariant the Wagner // and pool paths have always relied on, and harmless because a split is an @@ -128,6 +129,29 @@ void build_constraint_posthoc( // --- Node mapping and DFS timestamps --- +// Does the edge above a node whose descendant tip set is `nd` separate +// `together` from `apart`? It does when the set covers every tip of +// `together` and holds none of `apart`; the tips in neither group are free and +// are not looked at. With `apart` the exact complement of `together` the two +// conditions force set equality, which is the exact-clade test this replaced. +// +// THE definition of "displays a constraint split", shared by every entry point +// that has to decide it: the search/TBR mapping (map_constraint_nodes), the +// Wagner build's own check (wagner_tree_displays_constraint, ts_wagner.cpp), +// and the collapse pass's branch protection (ts_collapse_pool, ts_rcpp.cpp). +// They must not drift apart: the stricter of any two would reject trees +// another searches happily, or accept ones it will not move from. +inline bool node_displays_split( + const uint64_t* nd, const uint64_t* together, const uint64_t* apart, + int n_words) +{ + for (int w = 0; w < n_words; ++w) { + if ((together[w] & ~nd[w]) != 0ULL) return false; // a required tip missing + if ((apart[w] & nd[w]) != 0ULL) return false; // an excluded tip present + } + return true; +} + // Find which internal node holds each constraint split in the current tree. // Must be called after each accepted move and at search init. void map_constraint_nodes(const TreeState& tree, ConstraintData& cd); diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 0dfae2af7..6529c4a19 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2331,20 +2331,17 @@ List ts_collapse_pool( const uint64_t* R = &tb[static_cast(tree.right[ni]) * wps]; for (int w = 0; w < wps; ++w) dst[w] = L[w] | R[w]; } - auto displays = [&](const uint64_t* nb, const std::vector& in, - const std::vector& out) { - for (int w = 0; w < wps; ++w) { - if (in[w] & ~nb[w]) return false; // a required tip missing - if (out[w] & nb[w]) return false; // an excluded tip present - } - return true; - }; + // ts::node_displays_split() (ts_constraint.h) is the shared definition — + // the search's mapping and the Wagner build's check use the same one, so + // the branch protected here is the branch they enforce. for (size_t ci = 0; ci < cons_one.size(); ++ci) { + const uint64_t* one = cons_one[ci].data(); + const uint64_t* zero = cons_zero[ci].data(); int tight = -1, tight_size = n_tip + 1; for (int v = n_tip + 1; v < tree.n_node; ++v) { const uint64_t* nb = &tb[static_cast(v) * wps]; - if (!displays(nb, cons_one[ci], cons_zero[ci]) && - !displays(nb, cons_zero[ci], cons_one[ci])) continue; + if (!ts::node_displays_split(nb, one, zero, wps) && + !ts::node_displays_split(nb, zero, one, wps)) continue; int sz = 0; for (int w = 0; w < wps; ++w) sz += ts::popcount64(nb[w]); if (sz < tight_size) { tight_size = sz; tight = v; } diff --git a/src/ts_wagner.cpp b/src/ts_wagner.cpp index db59a118a..eb5a17c25 100644 --- a/src/ts_wagner.cpp +++ b/src/ts_wagner.cpp @@ -346,17 +346,16 @@ static int wagner_smallest_containing_node( // neither side is a clade, the partial tree does not display the split, and no // further insertion can make it. static int wagner_map_complement( - const TreeState& tree, int n_tip, int nw, + const TreeState& tree, int nw, const std::vector& node_tips, const uint64_t* split_out, const std::vector& added_tips, WagnerConstraintScratch& scratch) { - (void)n_tip; uint64_t* needed_out = scratch.needed_out.data(); int n_out = 0; int lone_out = -1; for (int w = 0; w < nw; ++w) { - // T-386: the outside group is cd.split_zeros, not ~split_tips. A `?`-coded + // #54: the outside group is cd.split_zeros, not ~split_tips. A `?`-coded // tip is in neither, so it never pulls this LCA about — which is what lets // it be placed on either side of the constraint, as documented. const uint64_t owd = split_out[w] & added_tips[w]; @@ -368,7 +367,7 @@ static int wagner_map_complement( } const int on = wagner_smallest_containing_node( tree, nw, node_tips, needed_out, n_out, lone_out); - return (on == n_tip) ? -1 : on; + return (on == tree.n_tip) ? -1 : on; } // Wagner-specific constraint node mapping. @@ -444,11 +443,11 @@ static void wagner_map_constraint_nodes( // cannot reach regraft_violates_constraint() before the next full remap. // constraint_node_hi is pinned to constraint_node for the same reason: // Wagner's LCA mapping has no displaying-chain, so the tight end is the - // only anchor it can honestly offer (T-386). + // only anchor it can honestly offer (#54). cd.constraint_complement[s] = 0; cd.constraint_node_hi[s] = cd.constraint_node[s]; scratch.outside_node[s] = wagner_map_complement( - tree, n_tip, nw, node_tips, split_out, added_tips, scratch); + tree, nw, node_tips, split_out, added_tips, scratch); continue; } @@ -491,7 +490,7 @@ static void wagner_map_constraint_nodes( if (inside_node == n_tip) { scratch.use_complement[s] = 1; scratch.outside_node[s] = wagner_map_complement( - tree, n_tip, nw, node_tips, split_out, added_tips, scratch); + tree, nw, node_tips, split_out, added_tips, scratch); } } } @@ -546,7 +545,7 @@ static void wagner_collect_active_splits( const bool tip_inside = (split[tw] >> tb) & 1; const bool tip_outside = (split_out[tw] >> tb) & 1; - // T-386: a tip in neither group is FREE — the constraint says nothing + // #54: a tip in neither group is FREE — the constraint says nothing // about which side of the separating edge it belongs on, so no edge is // barred to it. Reading "outside" as ~split_tips, as this did before free // tips were represented, forced every `?`-coded taxon out of the @@ -597,18 +596,24 @@ static void wagner_collect_active_splits( // Does the finished tree display every constraint split? // -// A split is displayed iff some edge separates its two groups, i.e. iff some -// node's subtree tip set covers one group and holds none of the other — the -// same free-taxa reading map_constraint_nodes() applies (T-386), spelled out -// again here because Wagner runs before any of that machinery is valid. It -// MUST stay in step with node_displays_split() in ts_constraint.cpp: this is -// the gate that decides whether the caller reshuffles and rebuilds, and the -// stricter of the two entry points would reject trees the other then searches -// happily (or, worse, accept ones it will not move from). +// A split is displayed iff some edge separates its two groups, which is +// node_displays_split() (ts_constraint.h) in one polarity or the other. This +// runs before any of map_constraint_nodes()'s machinery is valid, but it calls +// the same predicate rather than restating it: this is the gate that decides +// whether the caller reshuffles and rebuilds, and the stricter of the two entry +// points would reject trees the other then searches happily (or, worse, accept +// ones it will not move from). // -// Only non-root nodes are candidates: the root subtends every tip, and its two -// children already cover the single edge the degree-two root sits on. Tips are -// included so trivial (single-taxon) groups are recognised. +// Every node is a candidate, tips and root alike. Tips matter for a +// single-taxon group, whose "clade" is the tip itself. The root looks +// redundant — it subtends every tip, and its two children already cover the +// single edge the degree-two root sits on — but excluding it made this test +// STRICTER than find_displaying_chain(), which scans tree.postorder and so +// does reach the root: a split whose apart-group is empty maps there and +// nowhere else, and the two entry points then disagreed about a tree that +// every constraint is satisfied by. Including it costs one comparison and +// can never accept a violation, since the root displays a split only when the +// apart-group is empty, and then so does every tree. // Orientation-agnostic by construction, so it stays correct however the tree // happens to be rooted. static bool wagner_tree_displays_constraint(const TreeState& tree, @@ -634,18 +639,12 @@ static bool wagner_tree_displays_constraint(const TreeState& tree, const uint64_t* split_out = &cd.split_zeros[static_cast(s) * nw]; bool found = false; for (int node = 0; node < tree.n_node && !found; ++node) { - if (node == n_tip) continue; // root subtends everything const uint64_t* nd = &node_tips[static_cast(node) * nw]; - // Either group may be the clade side; free tips (in neither mask) are - // not looked at. No width mask is needed — both masks carry zeros above - // tip n_tip - 1, so padding bits can never make a test fail. - bool holds = true, holdsCompl = true; - for (int w = 0; w < nw; ++w) { - if ((split[w] & ~nd[w]) || (split_out[w] & nd[w])) holds = false; - if ((split_out[w] & ~nd[w]) || (split[w] & nd[w])) holdsCompl = false; - if (!holds && !holdsCompl) break; + // Either group may be the clade side. + if (node_displays_split(nd, split, split_out, nw) || + node_displays_split(nd, split_out, split, nw)) { + found = true; } - if (holds || holdsCompl) found = true; } if (!found) return false; } @@ -886,11 +885,14 @@ WagnerResult wagner_tree(TreeState& tree, const DataSet& ds, // been exhaustive. `constraint_fallback` alone is not enough: it only fires // when the filter rejected *every* edge, which the T-364/T-370 leak never did // -- it returned violating trees mutely. Nor is the caller's post-hoc check - // enough, because AdditionTree() never sets `has_posthoc` (it is built only - // at the search entry, ts_rcpp.cpp), so on that path there is no reshuffle to - // fall back on -- including for the both-sides-straddle case above. This - // check holds for every cause, known or not. The caller reports it: - // Rf_warning() is not safe from a search worker thread. + // enough: AdditionTree() reaches ts_wagner_tree(), which calls this function + // directly and so never runs random_wagner_tree()'s reshuffle loop. (It does + // build the posthoc DataSet -- R/AdditionTree.R splices the whole + // .PrepareConstraint() list through -- but nothing on that path consults it.) + // So on that path there is no retry to fall back on, including for the + // both-sides-straddle case above. This check holds for every cause, known or + // not. The caller reports it: Rf_warning() is not safe from a search worker + // thread. if (constrained) { result.constraint_violated = constraint_fallback || !wagner_tree_displays_constraint(tree, *cd); @@ -1202,6 +1204,15 @@ void random_topology_tree(TreeState& tree, const DataSet& ds) { // all constraint splits. (Uniform conditional on the split nesting // structure, which determines the partition of items across polytomy // resolution steps.) +// +// With free tips the "among those that satisfy" is narrower than the +// documented contract: this reads cd.split_tips only, so every free tip is a +// root-level item and the together-group comes out as an EXACT clade. That is +// strictly compliant, hence always legal (agent-issues/TreeSearch#54) — but it +// samples a strict subset of the legal topologies, so a free tip never starts +// inside the constrained group. Widening it would change which start trees +// the search sees, which is a search-quality change to measure on its own +// rather than a correctness fix to make here. namespace { diff --git a/tests/testthat/helper-ts.R b/tests/testthat/helper-ts.R index d121a8ded..2a714a1d3 100644 --- a/tests/testthat/helper-ts.R +++ b/tests/testthat/helper-ts.R @@ -100,3 +100,24 @@ validate_result <- function(result, n_tip) { tips <- sort(children[children <= n_tip]) testthat::expect_equal(tips, seq_len(n_tip)) } + +#' Run a constrained driven search with NOTHING but TBR enabled. +#' +#' Every phase that could rescue a replicate that cannot rearrange its start is +#' switched off, so `best_score` is what TBR alone achieved from `startEdge`: +#' a Wagner start is re-rooted on tip 0 (796a29d3), fuse re-roots its recipient, +#' and nni-perturb calls impose_constraint(), which repairs a start as a +#' side-effect. Used by the constraint tests that assert the search MOVES. +tbrOnlyRun <- function(ds, startEdge, splitMatrix) { + TreeSearch:::ts_driven_search( + ds$contrast, ds$tip_data, ds$weight, ds$levels, + maxReplicates = 1L, targetHits = 99L, tbrMaxHits = 1L, + ratchetCycles = 0L, driftCycles = 0L, nniPerturbCycles = 0L, + xssRounds = 0L, rssRounds = 0L, cssRounds = 0L, + pruneReinsertCycles = 0L, fuseInterval = 0L, + outerCycles = 1L, maxOuterResets = 0L, + nniFirst = FALSE, sprFirst = FALSE, + poolMaxSize = 100L, poolSuboptimal = 0, maxSeconds = 0, verbosity = 0L, + nThreads = 1L, startEdge = startEdge, consSplitMatrix = splitMatrix + ) +} diff --git a/tests/testthat/test-ts-constraint-free-taxa.R b/tests/testthat/test-ts-constraint-free-taxa.R index 6aed14d14..ca536949f 100644 --- a/tests/testthat/test-ts-constraint-free-taxa.R +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -1,7 +1,7 @@ # Tier 2: skipped on CRAN; see tests/testing-strategy.md skip_on_cran() -## T-386 / agent-issues/TreeSearch#54: the constraint the search enforces must +## agent-issues/TreeSearch#54: the constraint the search enforces must ## be the one `?MaximizeParsimony`'s `constraint` argument documents. ## ## The documented contract is the phyDat reading: a returned tree is compliant @@ -70,28 +70,11 @@ SeparatesGroups <- function(tree, labels, inGroup, outGroup) { # The engine hands back bare edge matrices for rooted binary trees. EdgeSeparatesGroups <- function(edge, labels, inGroup, outGroup) { tree <- structure( - list(edge = edge, Nnode = length(labels) - 1L, tip.label = labels), + list(edge = edge, Nnode = max(edge) - length(labels), tip.label = labels), class = "phylo") SeparatesGroups(tree, labels, inGroup, outGroup) } -# TBR only: every phase that could rescue a frozen replicate by starting -# somewhere else is switched off, so the reported score is what rearranging the -# supplied start achieved. Mirrors tbrOnlyRun() in test-ts-constraint-rooting.R. -tbrOnly <- function(ds, startEdge, splitMatrix) { - TreeSearch:::ts_driven_search( - ds$contrast, ds$tip_data, ds$weight, ds$levels, - maxReplicates = 1L, targetHits = 99L, tbrMaxHits = 1L, - ratchetCycles = 0L, driftCycles = 0L, nniPerturbCycles = 0L, - xssRounds = 0L, rssRounds = 0L, cssRounds = 0L, - pruneReinsertCycles = 0L, fuseInterval = 0L, - outerCycles = 1L, maxOuterResets = 0L, - nniFirst = FALSE, sprFirst = FALSE, - poolMaxSize = 100L, poolSuboptimal = 0, maxSeconds = 0, verbosity = 0L, - nThreads = 1L, startEdge = startEdge, consSplitMatrix = splitMatrix - ) -} - test_that("free `?` taxa do not freeze a compliant start tree", { dataset <- freeTaxaData() labels <- names(dataset) @@ -100,21 +83,21 @@ test_that("free `?` taxa do not freeze a compliant start tree", { startScore <- TreeLength(start, dataset) # Premise: the start satisfies the documented constraint but neither group is - # a clade, so the pre-T-386 exact-clade test could not map it. + # a clade, so the pre-#54 exact-clade test could not map it. expect_true(SeparatesGroups(start, labels, c("a", "b"), c("c", "d"))) expect_false(SeparatesGroups(start, labels, c("a", "b"), setdiff(labels, c("a", "b")))) # Take the split matrix from .PrepareConstraint rather than writing it out, # so this exercises the same R -> C++ contract the user's `constraint =` - # phyDat travels along: pre-T-386 it coded the free taxa 0 (making {a,b} an + # phyDat travels along: pre-#54 it coded the free taxa 0 (making {a,b} an # exact clade), now it codes them NA. free <- TreeSearch:::.PrepareConstraint( freeTaxaConstraint(), dataset)[["consSplitMatrix"]] expect_equal(as.vector(free), c(1L, 1L, 0L, 0L, NA, NA, NA, NA)) set.seed(386) - result <- tbrOnly(ds, start[["edge"]], free) + result <- tbrOnlyRun(ds, start[["edge"]], free) # The defect: every regraft was rejected and the start came back unimproved. expect_lt(result$best_score, startScore) @@ -136,7 +119,7 @@ test_that("a 0/1 constraint matrix still enforces the exact clade", { strict <- matrix(c(1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), nrow = 1) set.seed(386) - result <- tbrOnly(ds, start[["edge"]], strict) + result <- tbrOnlyRun(ds, start[["edge"]], strict) expect_true(all(vapply(result$trees, EdgeSeparatesGroups, logical(1), labels = labels, inGroup = c("a", "b"), @@ -222,12 +205,62 @@ test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { expect_equal(as.vector(consArgs[["consSplitMatrix"]]), c(1L, 1L, 0L, 0L, NA, NA, NA, NA)) - # A character with no `0` taxa constrains nothing under the documented - # contract — there is no group for the `1` taxa to be separated FROM — so it - # must not be enforced as a clade. - vacuous <- phangorn::phyDat( - matrix(c("1", "1", "?", "?", "?", "?", "?", "?"), nrow = 8, - dimnames = list(letters[1:8], NULL)), - type = "USER", levels = c("0", "1")) - expect_equal(TreeSearch:::.PrepareConstraint(vacuous, dataset), list()) + # A group of fewer than two taxa is separated from the rest by every tree, so + # such a character constrains nothing under the documented contract and must + # not be enforced as a clade. It is dropped, but not silently: coding only + # `1` and `?` almost always means "group these taxa", which is not what it + # says, and the alternative reading is the one that froze replicates. + Inert <- function(...) { + phangorn::phyDat(matrix(c(...), nrow = 8, + dimnames = list(letters[1:8], NULL)), + type = "USER", levels = c("0", "1")) + } + # No `0` group at all. + expect_warning( + dropped <- TreeSearch:::.PrepareConstraint( + Inert("1", "1", "?", "?", "?", "?", "?", "?"), dataset), + "constrains nothing") + expect_equal(dropped, list()) + # A `0` group of one. The two groups are interchangeable, so this must be + # treated exactly like its mirror image below -- which the old + # `1`-group-only test did not do. + expect_warning( + TreeSearch:::.PrepareConstraint( + Inert("1", "1", "0", "?", "?", "?", "?", "?"), dataset), + "constrains nothing") + expect_warning( + TreeSearch:::.PrepareConstraint( + Inert("0", "0", "1", "?", "?", "?", "?", "?"), dataset), + "constrains nothing") + # Two and two: kept, and kept silently. + expect_silent(TreeSearch:::.PrepareConstraint( + Inert("1", "1", "0", "0", "?", "?", "?", "?"), dataset)) +}) + +test_that("the Wagner build places free taxa freely", { + # wagner_tree_displays_constraint() and wagner_collect_active_splits() are a + # second, independent implementation of the same reading. AdditionTree() is + # the path with no post-hoc retry to fall back on (has_posthoc is set only at + # the search entry), so a Wagner build that read the constraint strictly + # would warn here — and, before the fix, was forced to place every `?` taxon + # outside the constrained group. + dataset <- freeTaxaData() + labels <- names(dataset) + cons <- freeTaxaConstraint() + + for (seed in 1:8) { + set.seed(seed) + tree <- expect_silent(AdditionTree(dataset, constraint = cons)) + expect_true(SeparatesGroups(tree, labels, c("a", "b"), c("c", "d")), + info = paste("seed", seed)) + } + + # A free taxon is genuinely free: over several addition orders the Wagner + # build is not forced to keep `e` out of the {a,b} group. + inGroup <- vapply(1:12, function(seed) { + set.seed(seed) + tr <- AdditionTree(dataset, constraint = cons) + SeparatesGroups(tr, labels, c("a", "b", "e"), c("c", "d")) + }, logical(1)) + expect_true(any(inGroup)) }) diff --git a/tests/testthat/test-ts-constraint-rooting.R b/tests/testthat/test-ts-constraint-rooting.R index 2ac2da24f..66cc9ebed 100644 --- a/tests/testthat/test-ts-constraint-rooting.R +++ b/tests/testthat/test-ts-constraint-rooting.R @@ -24,23 +24,8 @@ skip_on_cran() library("TreeTools") -# Everything that could rescue a rooting-dependent TBR is switched off: a Wagner -# start is re-rooted on tip 0 (796a29d3), fuse re-roots its recipient, and -# nni-perturb calls impose_constraint(), which repairs the rooting as a -# side-effect of repairing the split. -tbrOnlyRun <- function(ds, startEdge, splitMatrix) { - TreeSearch:::ts_driven_search( - ds$contrast, ds$tip_data, ds$weight, ds$levels, - maxReplicates = 1L, targetHits = 99L, tbrMaxHits = 1L, - ratchetCycles = 0L, driftCycles = 0L, nniPerturbCycles = 0L, - xssRounds = 0L, rssRounds = 0L, cssRounds = 0L, - pruneReinsertCycles = 0L, fuseInterval = 0L, - outerCycles = 1L, maxOuterResets = 0L, - nniFirst = FALSE, sprFirst = FALSE, - poolMaxSize = 100L, poolSuboptimal = 0, maxSeconds = 0, verbosity = 0L, - nThreads = 1L, startEdge = startEdge, consSplitMatrix = splitMatrix - ) -} +# tbrOnlyRun() lives in helper-ts.R: everything that could rescue a +# rooting-dependent TBR is switched off there, so the scores below are TBR's. # Phases whose timing must be zero, so a passing test cannot be one that quietly # searched its way around the mapping. Each is guarded by an explicit `> 0` diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index b561127a6..ef59db293 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -217,7 +217,7 @@ Per-strategy attempt and success counts are returned in the `strategy_diagnostics` attribute of the search result for post-hoc inspection. -### What a topological constraint requires +## Topological constraints A tree satisfies a constraint character when some edge separates the taxa coded `1` from those coded `0`. From 99ca3f049a6179fe8e41f9ff3255b9b1ec8b93ae Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 09:25:29 +0100 Subject: [PATCH 4/4] fix: close the review's remaining gaps in the constraint contract * The single-state constraint now warns. `.PrepareConstraint()` returned `list()` at the `nConsStates < 2` guard before the inert-character test could see it -- and that guard is exactly what `MatrixToPhyDat(c(a = "1", b = "1", c = "1"))` hits, the "make these a clade" idiom. It is the case the warning's own rationale names, and it was the one case that stayed silent. * The Wagner test asserted nothing. It checked that {a,b,e} ends up separated from {c,d}, which an exact {a,b} clade satisfies too, so it passed against the pre-fix build. It now measures the tightest node covering {a,b} and avoiding {c,d}: pre-fix that node is EXACTLY {a,b} in 25 of 25 seeds -- every `?` taxon forced out of the constrained clade -- and now holds a free taxon in all 25. * The over-loosening guard needed a guard: its start already has {a,b} as a clade, so a frozen search would have satisfied it for the wrong reason. It now asserts the score improved as well. * `random_constrained_tree()`'s new comment claimed its exact-clade sampling is "always legal". True only when the together-groups are laminar, which `.PrepareConstraint`'s four-gamete gate does not guarantee; the non-laminar case is handled by the T-329 collapse path and the post-hoc check, not by the claim. Corrected. Co-Authored-By: Claude Opus 5 --- R/MaximizeParsimony.R | 13 +++++- src/ts_wagner.cpp | 20 ++++++--- tests/testthat/test-ts-constraint-free-taxa.R | 44 ++++++++++++++++--- 3 files changed, 63 insertions(+), 14 deletions(-) diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 637f7337b..00b4f54a6 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -121,7 +121,18 @@ consContrast <- attr(constraint, "contrast") nConsStates <- ncol(consContrast) - if (nConsStates < 2L) return(list()) + if (nConsStates < 2L) { + # One state means no taxon is coded `0`, so this is the extreme case of the + # inert character warned about below -- and the loudest one, because it is + # what `MatrixToPhyDat(c(a = 1, b = 1, c = 1))` produces: a user asking for + # a clade and getting no constraint at all. Warn here rather than returning + # silently; the group-size test below never sees these characters. + warning("Constraint constrains nothing, and is ignored: no taxon is coded ", + "`0`, so every tree separates the `1` taxa from the (empty) `0` ", + "group. Code the taxa that must fall outside the group as `0`.", + call. = FALSE) + return(list()) + } consMat <- matrix(unlist(constraint, use.names = FALSE), nrow = length(constraint), byrow = TRUE) diff --git a/src/ts_wagner.cpp b/src/ts_wagner.cpp index eb5a17c25..da522c1e6 100644 --- a/src/ts_wagner.cpp +++ b/src/ts_wagner.cpp @@ -1207,12 +1207,20 @@ void random_topology_tree(TreeState& tree, const DataSet& ds) { // // With free tips the "among those that satisfy" is narrower than the // documented contract: this reads cd.split_tips only, so every free tip is a -// root-level item and the together-group comes out as an EXACT clade. That is -// strictly compliant, hence always legal (agent-issues/TreeSearch#54) — but it -// samples a strict subset of the legal topologies, so a free tip never starts -// inside the constrained group. Widening it would change which start trees -// the search sees, which is a search-quality change to measure on its own -// rather than a correctness fix to make here. +// root-level item and the together-group comes out as an EXACT clade, which +// displays the split under the free-taxa reading too +// (agent-issues/TreeSearch#54). So it samples a strict subset of the legal +// topologies, and a free tip never starts inside the constrained group. +// Widening it would change which start trees the search sees, which is a +// search-quality change to measure on its own rather than a correctness fix to +// make here. +// +// Making each together-group an exact clade is not always *possible*: the +// R-side gate (.PrepareConstraint) admits four-gamete-compatible splits that +// are not laminar, and those cannot all be clades at once. That case is +// handled by the collapse path below (a split that loses every tip to tighter +// non-laminar splits gets split_root == -1 and is skipped, T-329) and caught +// afterwards by the caller's post-hoc check, not by this comment's claim. namespace { diff --git a/tests/testthat/test-ts-constraint-free-taxa.R b/tests/testthat/test-ts-constraint-free-taxa.R index ca536949f..677246a24 100644 --- a/tests/testthat/test-ts-constraint-free-taxa.R +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -121,6 +121,9 @@ test_that("a 0/1 constraint matrix still enforces the exact clade", { set.seed(386) result <- tbrOnlyRun(ds, start[["edge"]], strict) + # Guard the guard: a search that froze would satisfy the compliance test + # below for the wrong reason, since the start already has {a,b} as a clade. + expect_lt(result$best_score, TreeLength(start, dataset)) expect_true(all(vapply(result$trees, EdgeSeparatesGroups, logical(1), labels = labels, inGroup = c("a", "b"), outGroup = setdiff(labels, c("a", "b"))))) @@ -235,6 +238,15 @@ test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { # Two and two: kept, and kept silently. expect_silent(TreeSearch:::.PrepareConstraint( Inert("1", "1", "0", "0", "?", "?", "?", "?"), dataset)) + + # The loudest case of all, and the one the group-size test never sees: a + # constraint with a single state, which is what MatrixToPhyDat() returns for + # the `c(a = 1, b = 1, c = 1)` "make these a clade" idiom. It reaches the + # `nConsStates < 2` early return, so it must warn there. + expect_warning( + TreeSearch:::.PrepareConstraint( + TreeTools::MatrixToPhyDat(c(a = "1", b = "1", c = "1")), dataset), + "constrains nothing") }) test_that("the Wagner build places free taxa freely", { @@ -255,12 +267,30 @@ test_that("the Wagner build places free taxa freely", { info = paste("seed", seed)) } - # A free taxon is genuinely free: over several addition orders the Wagner - # build is not forced to keep `e` out of the {a,b} group. - inGroup <- vapply(1:12, function(seed) { + # A free taxon is genuinely free. `wagner_collect_active_splits()` used to + # read "outside the split" as ~split_tips, which put every `?` taxon in the + # apart group and forced it out of the constrained clade: the tightest node + # covering {a,b} and avoiding {c,d} was EXACTLY {a,b} in 25 of 25 seeds. It + # now holds at least one free taxon in all 25. Asserting only that {a,b} and + # {c,d} end up separated would not detect this -- an exact {a,b} clade + # separates them too. + tightest <- vapply(1:12, function(seed) { set.seed(seed) - tr <- AdditionTree(dataset, constraint = cons) - SeparatesGroups(tr, labels, c("a", "b", "e"), c("c", "d")) - }, logical(1)) - expect_true(any(inGroup)) + splits <- as.logical(as.Splits(AdditionTree(dataset, constraint = cons), + tipLabels = labels)) + isOne <- labels %in% c("a", "b") + isZero <- labels %in% c("c", "d") + sizes <- c( + rowSums(splits)[apply(splits, 1, function(r) { + all(r[isOne]) && !any(r[isZero]) + })], + (length(labels) - rowSums(splits))[apply(splits, 1, function(r) { + !any(r[isOne]) && all(r[isZero]) + })] + ) + if (length(sizes)) min(sizes) else NA_integer_ + }, numeric(1)) + # Compliant in every seed (no NA), and never pinned to the bare `1` group. + expect_false(anyNA(tightest)) + expect_true(all(tightest > 2)) })