diff --git a/.AGENTS/memory/architecture.md b/.AGENTS/memory/architecture.md index 702ad3492..526caacfd 100644 --- a/.AGENTS/memory/architecture.md +++ b/.AGENTS/memory/architecture.md @@ -94,8 +94,34 @@ 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 whose `1` **or** + `0` group holds fewer than two taxa: vacuous under the documented contract, + since every tree separates such a group from the rest. +- A user constraint binds at three boundaries besides the per-move filter + (agent-issues/TreeSearch#59): the start tree is repaired by `impose_constraint()` + before it is scored, each replicate's finished tree is gated by + `capture_satisfies_constraint()` on its way into the pool, and the enforced + split is kept out of the final collapse. `impose_constraint()` is heuristic and + can fail, so every caller re-verifies. The collapse protects a realising node + only when no other realising node already survives — protecting + unconditionally would resolve a branch the constraint never asked for. - 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 c4372e120..61e82b9f2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,6 +18,29 @@ ended on a tree that could not be made to satisfy the constraint, and now raises an error rather than returning an unverified tree if no constraint-satisfying tree was found at all. +- Every part of the search now reads `constraint` the way it is documented: a + 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 locked-node filter that screens individual + rearrangements, the constrained Wagner build and the collapse pass previously + required the `1` group to be a clade *exactly*, free taxa excluded. That is + strictly stronger, so the search never accepted a rearrangement 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 exact match also blunted the collapse protection described above: with + free taxa it matched no branch, so the separating edge could still 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. - `TreeLength()`, `CharacterLength()`, `TreeScore()` and `EdgeListScore()` -- and so `Consistency()`, `ExpectedLength()`, `ConcordantInformation()`, `LengthAdded()` and `SuccessiveApproximations()`, which score trees through diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index fc01ab357..4c76daa60 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()) + } # Constraints are enforced as bipartitions, so only the two extreme states of # a character are read: taxa carrying an intermediate state are in neither @@ -156,10 +167,30 @@ } } - keep <- apply(consSplits, 1, function(row) { - s <- sum(row) - s >= 1 && s < length(constraint) - 1 - }) + # 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) + 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()) @@ -199,6 +230,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, consZero = consZero, @@ -234,6 +274,11 @@ # neither group here and are unconstrained there too (.PrepareConstraint() # warns about that at input). .ConstraintViolated <- function(tree, consOne, consZero) { + # `consOne` is the membership matrix the C++ kernels read, so it uses their + # coding: 1 = in the group, anything else -- including the NA that marks a + # free tip -- out of it. Reduce it to 0/1 here rather than let an NA + # propagate through the accumulation below and turn every comparison NA. + consOne <- (!is.na(consOne) & consOne == 1L) * 1L edge <- Postorder(tree)[["edge"]] parent <- edge[, 1L] child <- edge[, 2L] @@ -769,6 +814,14 @@ #' 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. +#' 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. #' Each constraint character is enforced as a single split, so one with more @@ -1748,21 +1801,17 @@ MaximizeParsimony <- function( # enforced clade"): a constraint is external evidence for a grouping the # matrix doesn't capture, so it stays visible even at zero length, while the # unsupported non-constraint branches still collapse. consSplitMatrix rows - # are the enforced bipartitions in tip_data order (see .PrepareConstraint). - # `consZero` names the tips the constraint places on the far side of the - # split; tips ambiguous for the character are in neither group. Without it - # the kernel can only recognise a node whose tip set is the 1 group exactly, - # and a split realised by any wider node goes unprotected -- collapsing the - # enforced grouping out of the returned tree. + # are the enforced bipartitions in tip_data order, carrying both groups + # (1 = together, 0 = apart, NA = free; see .PrepareConstraint). The kernel + # needs both: a tree with free tips generally realises the split at a node + # whose tip set is wider than the 1 group, which no exact match reaches, so + # the enforced grouping would collapse out of the returned tree. consSplits <- if (!is.null(constraintConfig)) { constraintConfig[["consSplitMatrix"]] } - consZero <- if (!is.null(constraintConfig)) { - constraintConfig[["consZero"]] - } collapsed <- ts_collapse_pool( bestTrees, contrast, tip_data, weight, levels, - scoringConfig, hsjConfig, xformConfig, consSplits, consZero + scoringConfig, hsjConfig, xformConfig, consSplits ) outTrees <- lapply(collapsed$trees, function(edgeMat) { tr <- list( diff --git a/R/RcppExports.R b/R/RcppExports.R index 46704fd61..dc9d764e1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -197,8 +197,8 @@ ts_driven_search <- function(contrast, tip_data, weight, levels, searchControl, .Call(`_TreeSearch_ts_driven_search`, contrast, tip_data, weight, levels, searchControl, runtimeConfig, scoringConfig, constraintConfig, hsjConfig, xformConfig) } -ts_collapse_pool <- function(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig = NULL, xformConfig = NULL, consSplitMatrix = NULL, consZero = NULL) { - .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix, consZero) +ts_collapse_pool <- function(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig = NULL, xformConfig = NULL, consSplitMatrix = NULL) { + .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix) } ts_resample_search <- function(contrast, tip_data, weight, levels, bootstrap = FALSE, jackProportion = 2.0 / 3.0, maxReplicates = 5L, targetHits = 2L, tbrMaxHits = 1L, ratchetCycles = 3L, ratchetPerturbProb = 0.04, driftCycles = 0L, min_steps = integer(), concavity = -1.0, consSplitMatrix = NULL, consContrast = NULL, consTipData = NULL, consWeight = NULL, consLevels = NULL, consExpectedScore = 0L, infoAmounts = NULL, xpiwe = FALSE, xpiwe_r = 0.5, xpiwe_max_f = 5.0, obs_count = integer()) { 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 105b21191..dcc28b8ad 100644 --- a/man/AdditionTree.Rd +++ b/man/AdditionTree.Rd @@ -37,6 +37,14 @@ 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. +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. Each constraint character is enforced as a single split, so one with more diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 56b99d993..11f33cfe2 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -143,6 +143,14 @@ 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. +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. Each constraint character is enforced as a single split, so one with more diff --git a/man/Resample.Rd b/man/Resample.Rd index 530a527af..98a744665 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -71,6 +71,14 @@ 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. +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. Each constraint character is enforced as a single split, so one with more diff --git a/man/SuccessiveApproximations.Rd b/man/SuccessiveApproximations.Rd index 6dfbce7c4..5d602dd61 100644 --- a/man/SuccessiveApproximations.Rd +++ b/man/SuccessiveApproximations.Rd @@ -83,6 +83,14 @@ 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. +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. Each constraint character is enforced as a single split, so one with more diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index cc19df6e8..dc52e9b59 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -590,8 +590,8 @@ BEGIN_RCPP END_RCPP } // ts_collapse_pool -List ts_collapse_pool(List edges, NumericMatrix contrast, IntegerMatrix tip_data, IntegerVector weight, CharacterVector levels, List scoringConfig, Nullable hsjConfig, Nullable xformConfig, Nullable consSplitMatrix, Nullable consZero); -RcppExport SEXP _TreeSearch_ts_collapse_pool(SEXP edgesSEXP, SEXP contrastSEXP, SEXP tip_dataSEXP, SEXP weightSEXP, SEXP levelsSEXP, SEXP scoringConfigSEXP, SEXP hsjConfigSEXP, SEXP xformConfigSEXP, SEXP consSplitMatrixSEXP, SEXP consZeroSEXP) { +List ts_collapse_pool(List edges, NumericMatrix contrast, IntegerMatrix tip_data, IntegerVector weight, CharacterVector levels, List scoringConfig, Nullable hsjConfig, Nullable xformConfig, Nullable consSplitMatrix); +RcppExport SEXP _TreeSearch_ts_collapse_pool(SEXP edgesSEXP, SEXP contrastSEXP, SEXP tip_dataSEXP, SEXP weightSEXP, SEXP levelsSEXP, SEXP scoringConfigSEXP, SEXP hsjConfigSEXP, SEXP xformConfigSEXP, SEXP consSplitMatrixSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -604,8 +604,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< Nullable >::type hsjConfig(hsjConfigSEXP); Rcpp::traits::input_parameter< Nullable >::type xformConfig(xformConfigSEXP); Rcpp::traits::input_parameter< Nullable >::type consSplitMatrix(consSplitMatrixSEXP); - Rcpp::traits::input_parameter< Nullable >::type consZero(consZeroSEXP); - rcpp_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix, consZero)); + rcpp_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix)); return rcpp_result_gen; END_RCPP } diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c index ceabefc71..7b166a616 100644 --- a/src/TreeSearch-init.c +++ b/src/TreeSearch-init.c @@ -60,7 +60,7 @@ extern SEXP _TreeSearch_ts_ev_cache_key_probe(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP extern SEXP _TreeSearch_ts_ls_fit(SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_ls_search(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_collapsed_flags_debug(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _TreeSearch_ts_collapse_pool(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _TreeSearch_ts_collapse_pool(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); static const R_CallMethodDef callMethods[] = { {"_TreeSearch_nni", (DL_FUNC) &_TreeSearch_nni, 3}, @@ -116,7 +116,7 @@ static const R_CallMethodDef callMethods[] = { {"_TreeSearch_ts_ls_fit", (DL_FUNC) &_TreeSearch_ts_ls_fit, 4}, {"_TreeSearch_ts_ls_search", (DL_FUNC) &_TreeSearch_ts_ls_search, 6}, {"_TreeSearch_ts_collapsed_flags_debug", (DL_FUNC) &_TreeSearch_ts_collapsed_flags_debug, 6}, - {"_TreeSearch_ts_collapse_pool", (DL_FUNC) &_TreeSearch_ts_collapse_pool, 10}, + {"_TreeSearch_ts_collapse_pool", (DL_FUNC) &_TreeSearch_ts_collapse_pool, 9}, {NULL, NULL, 0} }; diff --git a/src/ts_constraint.cpp b/src/ts_constraint.cpp index 114dc543b..881e7e492 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 (#54). 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 (#54). + 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,71 @@ 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 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) +// 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}. +// +// 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 +245,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 +252,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. + // + // #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 + // 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 +377,41 @@ 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 (#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) { - 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]; + + // 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) { - 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 +424,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 +468,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 +487,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 (#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 + // 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 +797,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 (#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; 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 +833,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 (#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) { + 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 +859,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 +880,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 +930,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 341c42d94..6ffd68e7b 100644 --- a/src/ts_constraint.h +++ b/src/ts_constraint.h @@ -3,10 +3,21 @@ // 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 +// 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-#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 +// already violates" — freezing the replicate on its start. // // Implementation: // 1. At init: store constraint splits as tip bitmasks. @@ -36,18 +47,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 (#54). When the caller + // supplies no free tips this is exactly ~split_tips, and every check below + // reduces to the pre-#54 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 +103,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 #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 +// unrooted bipartition whose two groups are interchangeable. ConstraintData build_constraint( const int* split_matrix, int n_splits, int n_tips); @@ -90,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_driven.cpp b/src/ts_driven.cpp index fae5ddfe3..f5b33c99a 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -54,16 +54,16 @@ ProgressInfo make_progress(int rep, const DrivenParams& params, // // violates_constraint_posthoc() answers that directly, but builds a whole // TreeState and scores it. For a BINARY constraint the locked-node mapping is -// much cheaper and is strictly the stronger test: it asks for the 1 group to be -// a clade exactly, excluding the taxa coded `?`, and a tree that manages that -// necessarily separates the two coded groups. So a full mapping settles the -// case the search puts us in almost every time -- every rearrangement it -// accepts is filtered on that same mapping -- and only an unmapped split pays -// for Fitch. +// much cheaper and asks exactly the same question: since #54 it maps a split to +// any node holding one whole group and none of the other, with the taxa coded +// `?` free to fall on either side, which is the separating edge itself. So a +// full mapping settles the case the search puts us in almost every time -- +// every rearrangement it accepts is filtered on that same mapping -- and only +// an unmapped split pays for Fitch. // -// With a third state the two tests diverge -- its taxa belong to no split_tips -// entry, so the character can sit above its minimum length with every split -// mapped -- and the mapping is the one to follow. It is the standard the rest +// With a third state the two tests diverge -- its taxa belong to neither group, +// so the character can sit above its minimum length with every split mapped -- +// and the mapping is the one to follow. It is the standard the rest // of the engine enforces: the locked-node filter screens rearrangements on it, // and impose_constraint() repairs to it and nothing more, so judging a capture // by the stricter Fitch check would discard every replicate of a search that diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 57b7959c3..a0da31485 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2251,8 +2251,7 @@ List ts_collapse_pool( List scoringConfig, Nullable hsjConfig = R_NilValue, Nullable xformConfig = R_NilValue, - Nullable consSplitMatrix = R_NilValue, - Nullable consZero = R_NilValue) + Nullable consSplitMatrix = R_NilValue) { ts::DataSet ds = unpack_scoring(contrast, tip_data, weight, levels, scoringConfig); @@ -2264,48 +2263,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. // - // The canonical bitsets alone protect only a node whose descendant set is the - // 1 group EXACTLY, which is what the search's locked-node machinery enforces. - // The constraint the user is promised is looser: tips ambiguous for the - // constraint character are free to sit on either side, so the split can be - // realised by a node that is not exactly the 1 group, which no exact match - // reaches — and contracting that node's edge takes the enforced grouping with - // it. `cons_one` / `cons_zero` are the raw (uncanonicalised) groups, from - // which the realising node is found per tree below. + // Each split is stored as the pair of groups build_constraint() reads + // (1 = together, 0 = apart, anything else = free; see ts_constraint.cpp). + // Both come out of the one membership matrix, so the protection here cannot + // drift from the constraint the search enforced. A pure 0/1 matrix gives + // apart == the complement, and the test below then fires on exactly the node + // an exact-match test would have found. const int n_tip = tip_data.nrow(); const int wps = (n_tip + 63) / 64; - std::vector> cons_canon; std::vector> cons_one, cons_zero; - auto row_bits = [&](const IntegerMatrix& m, int r) { - std::vector b(wps, 0); - for (int c = 0; c < n_tip && c < m.ncol(); ++c) { - if (m(r, c)) b[c >> 6] |= (1ULL << (c & 63)); - } - return b; - }; if (consSplitMatrix.isNotNull()) { IntegerMatrix cs(consSplitMatrix.get()); for (int r = 0; r < cs.nrow(); ++r) { - std::vector b = row_bits(cs, r); - cons_one.push_back(b); - 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 + std::vector one(wps, 0), zero(wps, 0); + for (int c = 0; c < n_tip && c < cs.ncol(); ++c) { + const int v = cs(r, c); + if (v != 1 && v != 0) continue; // free tip + (v == 1 ? one : zero)[c >> 6] |= (1ULL << (c & 63)); } - cons_canon.push_back(std::move(b)); + cons_one.push_back(std::move(one)); + cons_zero.push_back(std::move(zero)); } - if (consZero.isNotNull()) { - IntegerMatrix cz(consZero.get()); - for (int r = 0; r < cz.nrow() && r < cs.nrow(); ++r) { - cons_zero.push_back(row_bits(cz, r)); - } - } - cons_zero.resize(cons_one.size(), std::vector(wps, 0)); } // Group sizes depend only on the constraint, so they are counted once here // rather than per tree. A group of fewer than two taxa is skipped below: @@ -2368,11 +2348,11 @@ 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: keep an internal edge that realises each + // constraint out of the contraction, so the enforced grouping stays + // visible. Per-node descendant tip sets via a postorder OR; rooted on + // tip 0, so every internal set excludes tip 0. + 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); @@ -2386,23 +2366,16 @@ 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; } - } - } - - // A split can also be realised by a node that is not the 1 group exactly, - // and that node needs protecting too — but only when nothing else keeps - // the split visible. A node realises the split when it holds one whole - // group and none of the other; every such node's own edge displays it, so - // if any of them already survives the contraction there is nothing to do. + // A node realises the split when it holds one whole group and none of the + // other -- ts::node_displays_split() (ts_constraint.h), the same predicate + // the search's mapping and the Wagner build read, so the branch protected + // here is the branch they enforce. With free tips that node is generally + // NOT the 1 group exactly, and the exact-match test this replaced then + // protected nothing at all (agent-issues/TreeSearch#54). + // + // Protect one such node, and only when nothing else keeps the split + // visible: every realising node's own edge displays the split, so if any + // of them already survives the contraction there is nothing to do. // Protecting unconditionally would instead force the resolution of a // branch the constraint does not ask for, which is the "unsupported // non-constraint branches still collapse" half of the promise. @@ -2419,20 +2392,13 @@ List ts_collapse_pool( bool survives = false; int to_protect = -1; for (int side = 0; side < 2 && !survives; ++side) { - const std::vector& in = *grp[side]; - const std::vector& out = *grp[1 - side]; + const uint64_t* in = grp[side]->data(); + const uint64_t* out = grp[1 - side]->data(); for (size_t pi = 0; pi < tree.postorder.size(); ++pi) { const int v = tree.postorder[pi]; if (v <= n_tip || v >= static_cast(flags.size())) continue; const uint64_t* nb = &tb[static_cast(v) * wps]; - bool realises = true; - for (int w = 0; w < wps; ++w) { - if ((nb[w] & in[w]) != in[w] || (nb[w] & out[w])) { - realises = false; - break; - } - } - if (!realises) continue; + if (!ts::node_displays_split(nb, in, out, wps)) continue; if (!flags[v]) { survives = true; break; } if (to_protect < 0) to_protect = v; // the MRCA, in postorder } diff --git a/src/ts_wagner.cpp b/src/ts_wagner.cpp index ecd2ddcaa..d48e54f54 100644 --- a/src/ts_wagner.cpp +++ b/src/ts_wagner.cpp @@ -346,8 +346,8 @@ 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 std::vector& node_tips, const uint64_t* split, + const TreeState& tree, int nw, + const std::vector& node_tips, const uint64_t* split_out, const std::vector& added_tips, WagnerConstraintScratch& scratch) { @@ -355,12 +355,10 @@ static int wagner_map_complement( 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]; + // #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]; needed_out[w] = owd; if (owd) { n_out += popcount64(owd); @@ -369,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. @@ -426,6 +424,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 +441,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 (#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, added_tips, scratch); + tree, nw, node_tips, split_out, added_tips, scratch); continue; } @@ -466,6 +470,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 +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, added_tips, scratch); + tree, nw, node_tips, split_out, added_tips, scratch); } } } @@ -534,8 +539,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; + // #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 + // 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 +559,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 +596,26 @@ 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, 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). +// +// 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, const ConstraintData& cd) { const int n_tip = tree.n_tip; @@ -612,22 +636,15 @@ 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; - 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; + // 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 (eq || eqCompl) found = true; } if (!found) return false; } @@ -868,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); @@ -1190,6 +1210,23 @@ 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, 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/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 new file mode 100644 index 000000000..677246a24 --- /dev/null +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -0,0 +1,296 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +## 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 = max(edge) - length(labels), tip.label = labels), + class = "phylo") + SeparatesGroups(tree, labels, inGroup, outGroup) +} + +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-#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-#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 <- tbrOnlyRun(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 <- 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"))))) +}) + +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("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() + + 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 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)) + + # 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", { + # 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. `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) + 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)) +}) 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 605633638..bb06dc501 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. +## Topological constraints + +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 is then treated exactly like a violating start (below) -- +frozen, or repaired to a topology it did not need to be moved to. +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. + ### Starting trees under a constraint Whatever its source -- a Wagner build, a random topology, or a tree supplied @@ -238,16 +268,9 @@ so a tree that breaks the constraint is never returned; and the enforced splits are protected from the final collapse pass, whichever branch happens to realise them. -"Satisfies the constraint" here means what `constraint` promises the user: some -edge separates the taxa coded `1` for a constraint character from those coded -`0`, with `?` taxa free to sit on either side. -Note that the locked-node filter used to screen individual rearrangements reads -the constraint more strictly, as "the `1` group is a clade exactly", excluding -the free taxa. -Every strictly-compliant tree satisfies the user's constraint, so the search -never returns a tree that breaks it; but a start that satisfies the user's -constraint without satisfying the stricter form maps to no node, and the -replicate makes no moves from it. +A start that already separates the two groups is left alone, however its free +taxa happen to fall; the repair fires only on a start that genuinely violates +the constraint. ## The driven search pipeline