diff --git a/NEWS.md b/NEWS.md index a08d94daf..19ad92e64 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,23 @@ # To integrate into 2.0.0 notes +- `constraint` now binds the trees `MaximizeParsimony()` returns, at three + boundaries where it did not. A starting tree supplied through `tree` was + never checked against the constraint; because a constrained search rejects + every rearrangement away from a violating tree, the replicate froze on it and + reported a score no constraint-satisfying tree could reach, which then evicted + the compliant trees other replicates had found. A violating start is now + rearranged until it complies before the search begins, **with a warning**. + Separately, a replicate's own tree entered the pool unchecked, and the final + collapse of unsupported branches could contract the very branch that displayed + an enforced grouping -- so under the default `collapse = TRUE` a returned tree + could break the constraint outright. Both paths are now checked. + + **Constrained results may therefore differ from previous versions**: scores + can rise to the true constrained optimum, and returned trees will display the + constrained groupings. `MaximizeParsimony()` also warns if any replicate + 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. - `TreeLength()`, `CharacterLength()`, `TreeScore()` and `EdgeListScore()` -- and so `Consistency()`, `ExpectedLength()`, `ConcordantInformation()`, `LengthAdded()` and `SuccessiveApproximations()`, which score trees through diff --git a/R/AdditionTree.R b/R/AdditionTree.R index cbdaac288..ee1ce13cc 100644 --- a/R/AdditionTree.R +++ b/R/AdditionTree.R @@ -154,7 +154,9 @@ AdditionTree <- function(dataset, concavity = Inf, constraint, sequence) { min_steps = minSteps, concavity = as.double(concavity) ) - result <- do.call(ts_wagner_tree, c(searchArgs, consArgs, profileArgs)) + result <- do.call(ts_wagner_tree, + c(searchArgs, .KernelConstraintArgs(consArgs), + profileArgs)) # Reconstruct phylo from edge matrix tree <- list( diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 3ccce1f44..fc01ab357 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -123,6 +123,17 @@ nConsStates <- ncol(consContrast) if (nConsStates < 2L) 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 + # group and go unconstrained. Say so rather than let the caller infer, from + # `@param constraint`'s "compatible with each character", that a third state + # groups its taxa too. + if (nConsStates > 2L) { + warning("`constraint` characters with more than two states are enforced ", + "as the split between their first and last state only; taxa in ", + "any intermediate state are left unconstrained.", call. = FALSE) + } + consMat <- matrix(unlist(constraint, use.names = FALSE), nrow = length(constraint), byrow = TRUE) # For each constraint character, record the tips unambiguously in the "1" @@ -190,6 +201,7 @@ list( consSplitMatrix = consSplits, + consZero = consZero, consContrast = consContrast, consTipData = consTipData, consWeight = as.integer(consWeight), @@ -198,6 +210,63 @@ ) } +# Constraint fields the flat `ts_*` kernels declare as formals, in contrast to +# the list-config entry points, which ignore anything they do not name. A +# `do.call()` onto a flat kernel has to be filtered through this, or a field +# added for the list-config path becomes an unused-argument error there. +.kernelConsFields <- c("consSplitMatrix", "consContrast", "consTipData", + "consWeight", "consLevels", "consExpectedScore") + +.KernelConstraintArgs <- function(consArgs) { + consArgs[intersect(names(consArgs), .kernelConsFields)] +} + +# Does `tree` display a split separating a constraint character's "1" group +# from its "0" group? This is the phyDat reading `constraint` is documented +# in: tips ambiguous for the character sit on either side, so the test is +# "some edge separates the two groups", not the stricter "the 1 group is +# exactly a clade" that the search's locked-node machinery enforces +# internally. `consOne` / `consZero` are .PrepareConstraint()'s matrices, in +# `tip_data` column order; `tree`'s tips must already be renumbered to match. +# +# The two groups are the character's extreme states, so this answers for +# exactly what the engine enforces -- an intermediate state's taxa are in +# neither group here and are unconstrained there too (.PrepareConstraint() +# warns about that at input). +.ConstraintViolated <- function(tree, consOne, consZero) { + edge <- Postorder(tree)[["edge"]] + parent <- edge[, 1L] + child <- edge[, 2L] + nTip <- ncol(consOne) + nRow <- nrow(consOne) + # One accumulation pass carries every group at once: rows 1..nRow are the + # "1" groups, the rest the "0" groups. Nodes index the COLUMNS, so each + # accumulation touches one contiguous stretch of a column-major matrix. + counts <- matrix(0L, nrow = 2L * nRow, ncol = max(edge)) + counts[, seq_len(nTip)] <- rbind(consOne, consZero) + for (i in seq_along(parent)) { + counts[, parent[i]] <- counts[, parent[i]] + counts[, child[i]] + } + # Postorder lists every node before its parent, so the first node holding a + # whole group is that group's MRCA; the groups are separated iff one MRCA + # holds none of the other group. + nodes <- c(child, parent[length(parent)]) + for (r in seq_len(nRow)) { + one <- counts[r, ] + zero <- counts[nRow + r, ] + nOne <- sum(consOne[r, ]) + nZero <- sum(consZero[r, ]) + mrcaOne <- nodes[one[nodes] == nOne][1] + mrcaZero <- nodes[zero[nodes] == nZero][1] + displayed <- (!is.na(mrcaOne) && zero[mrcaOne] == 0L) || + (!is.na(mrcaZero) && one[mrcaZero] == 0L) + if (!displayed) { + return(TRUE) + } + } + FALSE +} + # Ratchet depth for implied weights under `thorough`/`large`, applied after the # strategy preset (see MaximizeParsimony()). Kept out of `.StrategyPresets()` so # the preset table stays scorer-agnostic: this depth is calibrated for implied @@ -615,6 +684,12 @@ #' says so in a warning. Raise `targetHits` as well as `maxReplicates` to #' use more of it. #' If unspecified, all replicates start from random Wagner trees. +#' A start tree that does not satisfy `constraint` is rearranged until it +#' does before the search begins, with a warning: `constraint` is a +#' guarantee about the trees returned, whereas `tree` only says where to +#' begin, so when the two conflict the guarantee wins. A taxon coded `?` +#' for a constraint character is unconstrained by it and may start on +#' either side of that split. #' Edge lengths are not supported and will be deleted. #' Rooted and unrooted trees are both accepted; an unrooted tree is rooted #' arbitrarily (on its first tip) before the search begins, which may @@ -696,6 +771,9 @@ #' in any output tree. #' 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 +#' than two states is read as the split between its first and last state +#' alone: taxa in an intermediate state are left unconstrained, with a warning. #' @param effort Integer: how much search effort to spend, **relative to the #' amount chosen automatically** for this dataset. `0` (the default) accepts #' the automatic choice; `1` asks for one notch more, `-1` one notch less. @@ -1520,6 +1598,23 @@ MaximizeParsimony <- function( cli_alert_info("Constraint: {nrow(consArgs$consSplitMatrix)} split{?s}") } + # A start tree that breaks the constraint is not something the search can + # rearrange its way out of -- every constrained move from it is rejected, so + # it would freeze the replicate on a tree scoring better than any legal one. + # The engine repairs such a start before scoring it, but the conflict is the + # caller's to know about: either `tree` or `constraint` is not what they + # meant, and the tree they get back will not be the one they supplied. + if (userTree && length(consArgs) > 0L) { + violating <- vapply(startTrees, .ConstraintViolated, logical(1), + consArgs[["consSplitMatrix"]], consArgs[["consZero"]]) + if (any(violating)) { + warning(sum(violating), " of the ", length(startTrees), + " tree(s) supplied to `tree` do not satisfy `constraint`; ", + "they will be rearranged to comply before the search starts, ", + "or replaced if that fails.", call. = FALSE) + } + } + # --- Profile parsimony: extract info_amounts --- profileArgs <- list() if (useProfile) { @@ -1654,12 +1749,20 @@ MaximizeParsimony <- function( # 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. 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 + scoringConfig, hsjConfig, xformConfig, consSplits, consZero ) outTrees <- lapply(collapsed$trees, function(edgeMat) { tr <- list( @@ -1680,6 +1783,18 @@ MaximizeParsimony <- function( }) } if (length(outTrees) == 0L) { + # `treeTpl` is a starting tree, so under a constraint it is exactly what may + # not be handed back unchecked: an empty pool means no replicate produced a + # tree the constraint gate accepted -- or, benignly, that the time limit + # expired before the first one finished. Check rather than assume, so a + # short budget still returns a tree when the fallback happens to comply. + if (!is.null(constraintConfig) && + .ConstraintViolated(treeTpl, constraintConfig[["consSplitMatrix"]], + constraintConfig[["consZero"]])) { + stop("The search returned no tree satisfying `constraint`. Check that ", + "the constraint is compatible with the data, and allow more search ", + "with `maxReplicates` or `maxSeconds`.") + } outTrees <- list(treeTpl) } diff --git a/R/RcppExports.R b/R/RcppExports.R index dc9d764e1..46704fd61 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) { - .Call(`_TreeSearch_ts_collapse_pool`, edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix) +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_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/R/Resample.R b/R/Resample.R index 820e23e2b..22d316204 100644 --- a/R/Resample.R +++ b/R/Resample.R @@ -369,7 +369,7 @@ Resample <- function(dataset, tree, method = "jack", proportion = 2 / 3, # Batch mode: run all replicates at once (optionally in parallel) batchArgs <- c(searchArgs, list(nReplicates = nReplicates, nThreads = nThreads), - consArgs, profileArgs) + .KernelConstraintArgs(consArgs), profileArgs) result <- do.call(ts_parallel_resample, batchArgs) trees <- vector("list", nReplicates) @@ -393,7 +393,9 @@ Resample <- function(dataset, tree, method = "jack", proportion = 2 / 3, } # Single-replicate path (original behavior) - result <- do.call(ts_resample_search, c(searchArgs, consArgs, profileArgs)) + result <- do.call(ts_resample_search, + c(searchArgs, .KernelConstraintArgs(consArgs), + profileArgs)) if (nrow(result$edge) == 0L) { tr <- if (!missing(tree) && inherits(tree, "phylo")) tree diff --git a/R/SuccessiveApproximations.R b/R/SuccessiveApproximations.R index d5a414645..52361f9c5 100644 --- a/R/SuccessiveApproximations.R +++ b/R/SuccessiveApproximations.R @@ -104,7 +104,9 @@ SuccessiveApproximations <- function (tree, dataset, outgroup = NULL, k = 3, xpiwe_max_f = as.double(xpiwe_max_f), obs_count = if (useXpiwe) obsCount else integer(0) ) - result <- do.call(ts_successive_approx, c(searchArgs, consArgs, profileArgs)) + result <- do.call(ts_successive_approx, + c(searchArgs, .KernelConstraintArgs(consArgs), + profileArgs)) if (result$converged && verbosity > 0) { message("Successive approximations converged after ", diff --git a/man/AdditionTree.Rd b/man/AdditionTree.Rd index c810ac7d3..105b21191 100644 --- a/man/AdditionTree.Rd +++ b/man/AdditionTree.Rd @@ -38,7 +38,10 @@ 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. Constraint searches are supported natively: all tree rearrangements -are filtered to respect the constraint topology.} +are filtered to respect the constraint topology. +Each constraint character is enforced as a single split, so one with more +than two states is read as the split between its first and last state +alone: taxa in an intermediate state are left unconstrained, with a warning.} \item{sequence}{Character or numeric vector listing sequence in which to add taxa. Randomized if not provided.} diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 917751868..56b99d993 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -50,6 +50,12 @@ limit, whichever fires first — draws on only part of a large pool, and says so in a warning. Raise \code{targetHits} as well as \code{maxReplicates} to use more of it. If unspecified, all replicates start from random Wagner trees. +A start tree that does not satisfy \code{constraint} is rearranged until it +does before the search begins, with a warning: \code{constraint} is a +guarantee about the trees returned, whereas \code{tree} only says where to +begin, so when the two conflict the guarantee wins. A taxon coded \verb{?} +for a constraint character is unconstrained by it and may start on +either side of that split. Edge lengths are not supported and will be deleted. Rooted and unrooted trees are both accepted; an unrooted tree is rooted arbitrarily (on its first tip) before the search begins, which may @@ -138,7 +144,10 @@ 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. Constraint searches are supported natively: all tree rearrangements -are filtered to respect the constraint topology.} +are filtered to respect the constraint topology. +Each constraint character is enforced as a single split, so one with more +than two states is read as the split between its first and last state +alone: taxa in an intermediate state are left unconstrained, with a warning.} \item{effort}{Integer: how much search effort to spend, \strong{relative to the amount chosen automatically} for this dataset. \code{0} (the default) accepts diff --git a/man/Resample.Rd b/man/Resample.Rd index e99630af0..530a527af 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -72,7 +72,10 @@ 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. Constraint searches are supported natively: all tree rearrangements -are filtered to respect the constraint topology.} +are filtered to respect the constraint topology. +Each constraint character is enforced as a single split, so one with more +than two states is read as the split between its first and last state +alone: taxa in an intermediate state are left unconstrained, with a warning.} \item{verbosity}{Integer specifying level of messaging; higher values give more detail. Set to \code{0} to run silently. diff --git a/man/SuccessiveApproximations.Rd b/man/SuccessiveApproximations.Rd index cdd4e3e8c..6dfbce7c4 100644 --- a/man/SuccessiveApproximations.Rd +++ b/man/SuccessiveApproximations.Rd @@ -84,7 +84,10 @@ 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. Constraint searches are supported natively: all tree rearrangements -are filtered to respect the constraint topology.} +are filtered to respect the constraint topology. +Each constraint character is enforced as a single split, so one with more +than two states is read as the split between its first and last state +alone: taxa in an intermediate state are left unconstrained, with a warning.} \item{extended_iw}{Logical: if \code{TRUE} (default) and \code{concavity} is finite, apply the missing-entries correction of diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dc52e9b59..cc19df6e8 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); -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) { +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) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -604,7 +604,8 @@ 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_result_gen = Rcpp::wrap(ts_collapse_pool(edges, contrast, tip_data, weight, levels, scoringConfig, hsjConfig, xformConfig, consSplitMatrix)); + 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)); return rcpp_result_gen; END_RCPP } diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c index 7b166a616..ceabefc71 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); +extern SEXP _TreeSearch_ts_collapse_pool(SEXP, 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, 9}, + {"_TreeSearch_ts_collapse_pool", (DL_FUNC) &_TreeSearch_ts_collapse_pool, 10}, {NULL, NULL, 0} }; diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index a2ad975b8..fae5ddfe3 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -49,8 +49,61 @@ ProgressInfo make_progress(int rep, const DrivenParams& params, return pi; } +// Does the tree satisfy the user constraint -- some edge separating the taxa +// coded 1 for each constraint character from those coded 0? +// +// 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. +// +// 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 +// 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 +// cannot produce anything better. The R layer warns at input that an +// intermediate state goes unconstrained. +// +// update_constraint(), not map_constraint_nodes(): the DFS timestamps have to +// move with the node ids, or a consumer that reads both without re-mapping +// (spr_search) sees this tree's nodes against another tree's timestamps. +bool constraint_satisfied(TreeState& tree, ConstraintData& cd) { + update_constraint(tree, cd); + for (int s = 0; s < cd.n_splits; ++s) { + if (cd.constraint_node[s] < 0) { + return !violates_constraint_posthoc(tree, cd); + } + } + return true; +} + } // anonymous namespace +bool capture_satisfies_constraint(TreeState& tree, ConstraintData* cd, + const DataSet& ds, double& score) +{ + // Gate on the post-hoc DataSet, which only a *user* constraint carries. The + // cross-replicate consensus constraint is a search heuristic, not a promise + // about the answer, so a tree that breaks it is not a wrong result and must + // not be thrown away. The post-hoc check is also the right test even for a + // user constraint: a tree can map every constraint node and still fail the + // full Fitch check, which is the case the post-hoc DataSet exists for. + if (!cd || !cd->active || !cd->has_posthoc) return true; + if (constraint_satisfied(tree, *cd)) return true; + + impose_constraint(tree, *cd); + tree.build_postorder(); + tree.reset_states(ds); + score = score_tree(tree, ds); + return constraint_satisfied(tree, *cd); +} + // --- Single-replicate pipeline --- ReplicateResult run_single_replicate( @@ -151,6 +204,37 @@ ReplicateResult run_single_replicate( } } + // A start that breaks the constraint has to be repaired here, before anything + // takes its score as a baseline. Constrained rearrangement cannot undo it: + // regraft_violates_constraint() reads an unmapped split as "already + // violating" and rejects every move, so the search freezes on the start and + // reports its unconstrained — and therefore unbeatably low — score. Nor can + // a later verify-and-revert gate help, for the same reason: the repaired tree + // is legal and so necessarily scores worse than the violation it replaces. + // The R layer warns when a caller's `tree` is what arrived here violating. + if (cd && cd->active && cd->has_posthoc && + !constraint_satisfied(result.tree, *cd)) { + impose_constraint(result.tree, *cd); + result.tree.build_postorder(); + result.tree.reset_states(ds); + if (!constraint_satisfied(result.tree, *cd)) { + // impose_constraint() is heuristic. Discard the start rather than search + // from a tree the constraint machinery cannot move: a constrained Wagner + // build, with its own post-hoc reshuffles, is the better bet. It is not + // a guarantee either -- exhausting those reshuffles returns a violating + // tree -- so repair whatever it hands back rather than trusting it. + random_wagner_tree(result.tree, ds, cd); + result.tree.build_postorder(); + result.tree.reset_states(ds); + if (!constraint_satisfied(result.tree, *cd)) { + impose_constraint(result.tree, *cd); + result.tree.build_postorder(); + result.tree.reset_states(ds); + } + } + best_wag = score_tree(result.tree, ds); + } + result.timings.wagner_ms = ph_lap(); if (verbosity >= 2) { if (starting_tree) { @@ -1055,24 +1139,48 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, result.timings += rep_result.timings; - // Compute collapsed flags for collapsed-topology pool dedup. - // Trees that differ only in zero-length resolutions are treated - // as duplicates, improving pool diversity (Goloboff & Farris 2001). - std::vector rep_collapsed; - compute_collapsed_flags(rep_result.tree, ds, rep_collapsed); - if (rep_result.interrupted) { - if (rep_result.score < 1e18) { + // Tested but not repaired: the deadline has already passed, and + // impose_constraint() is an unbounded SPR loop with no interrupt check + // of its own, so repairing here would extend an overrun. + const bool keep = !cd || !cd->active || !cd->has_posthoc || + constraint_satisfied(rep_result.tree, *cd); + if (keep && rep_result.score < 1e18) { + std::vector rep_collapsed; + compute_collapsed_flags(rep_result.tree, ds, rep_collapsed); pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); } result.timed_out = true; goto finish; } - // Add to pool with collapsed-topology dedup - double prev_best = pool.best_score(); - pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); - bool score_improved = pool.best_score() < prev_best; + // A replicate can still finish on a constraint-violating tree: a Wagner + // start whose reshuffles all failed, or a phase that accepts on a looser + // check than the pool promises. The pool is what the caller is handed, so + // gate it here, as the fuse capture below already does. + const bool rep_ok = capture_satisfies_constraint(rep_result.tree, cd, ds, + rep_result.score); + if (!rep_ok) ++result.constraint_discards; + + // A discarded replicate contributes its count and nothing else. Its score + // is that of a violating tree, which beats any legal one, so letting it + // through would credit the strategy arm that produced it, bias the coverage + // estimate downwards and report a figure no returned tree attains. The + // stopping rules at the foot of the loop still run: skipping them would + // outlive the deadline and swallow an interrupt. + bool score_improved = false; + if (rep_ok) { + // Compute collapsed flags for collapsed-topology pool dedup. + // Trees that differ only in zero-length resolutions are treated + // as duplicates, improving pool diversity (Goloboff & Farris 2001). + std::vector rep_collapsed; + compute_collapsed_flags(rep_result.tree, ds, rep_collapsed); + + // Add to pool with collapsed-topology dedup + double prev_best = pool.best_score(); + pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); + score_improved = pool.best_score() < prev_best; + } if (score_improved) { result.last_improved_rep = rep1; unsuccessful_reps = 0; @@ -1085,7 +1193,7 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, // not use a fresh-start arm, so crediting/blaming one would corrupt the // bandit. Together these two flags mean exactly `start_ptr == nullptr`; // any future warm-start source must be excluded here too. - if (params.adaptive_start && !pr_reseeded && !user_started) { + if (params.adaptive_start && !pr_reseeded && !user_started && rep_ok) { bool hit_best = (rep_result.score <= pool.best_score()); strategy_tracker.update(rep_strategy, hit_best); } @@ -1098,10 +1206,10 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, } ++result.replicates_completed; - result.replicate_scores.push_back(rep_result.score); - - // Report end of replicate - report("replicate", 1, rep_result.score, rep1); + if (rep_ok) { + result.replicate_scores.push_back(rep_result.score); + report("replicate", 1, rep_result.score, rep1); + } // Periodic tree fusing if (params.fuse_interval > 0 && diff --git a/src/ts_driven.h b/src/ts_driven.h index 44f3ff93e..3a54c0e80 100644 --- a/src/ts_driven.h +++ b/src/ts_driven.h @@ -359,6 +359,11 @@ struct DrivenResult { // search (TNT "Total rearrangements examined" analogue). Serial path only; // 0 when run in parallel. See DataSet::n_candidates_evaluated. long long candidates_evaluated = 0; + + // Replicates whose finished tree still violated the user constraint after + // repair and so never entered the pool (see capture_satisfies_constraint). + // Reported by the caller: Rf_warning() is not safe from a worker thread. + int constraint_discards = 0; }; // Result of a single replicate (tree + score, no pool interaction). @@ -391,6 +396,14 @@ ReplicateResult run_single_replicate( StartStrategy strategy = StartStrategy::WAGNER_RANDOM, const TreePool* pool = nullptr); +// Gate a finished replicate's tree on its way into the pool. Mirrors the fuse +// capture: repair a constraint violation, verify the repair took, and return +// false when it did not, so nothing the caller is handed can break the +// constraint the caller asked for. `score` is refreshed when a repair moves +// the tree. Inert unless a *user* constraint is active. +bool capture_satisfies_constraint(TreeState& tree, ConstraintData* cd, + const DataSet& ds, double& score); + // Run the full driven search. Returns search statistics. // The pool contents (all retained trees) are accessible via the pool // reference stored in `pool_out`. Caller should extract edge matrices. diff --git a/src/ts_parallel.cpp b/src/ts_parallel.cpp index 1b6cf2e60..4035d0f30 100644 --- a/src/ts_parallel.cpp +++ b/src/ts_parallel.cpp @@ -153,6 +153,9 @@ struct WorkerContext { // Per-thread score accumulator (index = thread_id) std::vector* thread_scores; + // Per-thread count of replicates dropped by the constraint capture gate + int* thread_constraint_discards; + // Wall-clock deadline for sector phases (only meaningful when use_timeout) bool use_timeout; std::chrono::steady_clock::time_point deadline; @@ -230,14 +233,21 @@ void worker_thread(WorkerContext ctx) { // Accumulate phase timings for this thread ctx.thread_timings[ctx.thread_id] += rep_result.timings; - // Add to shared pool with collapsed-topology dedup - std::vector rep_collapsed; - compute_collapsed_flags(rep_result.tree, ds_local, rep_collapsed); - ctx.shared_pool->add_collapsed(rep_result.tree, rep_result.score, - rep_collapsed); - - // Record per-replicate score for Chao1 coverage estimation - ctx.thread_scores[ctx.thread_id].push_back(rep_result.score); + // Add to shared pool with collapsed-topology dedup, gated on the + // constraint exactly as the serial driver's capture is. + if (capture_satisfies_constraint(rep_result.tree, cd_ptr, ds_local, + rep_result.score)) { + std::vector rep_collapsed; + compute_collapsed_flags(rep_result.tree, ds_local, rep_collapsed); + ctx.shared_pool->add_collapsed(rep_result.tree, rep_result.score, + rep_collapsed); + // Record per-replicate score for Chao1 coverage estimation. A discarded + // replicate is left out: its score is a violating tree's, which no + // returned tree attains. + ctx.thread_scores[ctx.thread_id].push_back(rep_result.score); + } else { + ++ctx.thread_constraint_discards[ctx.thread_id]; + } ctx.replicates_done->fetch_add(1, std::memory_order_relaxed); @@ -361,6 +371,7 @@ DrivenResult parallel_driven_search( // Per-thread timing and score accumulators std::vector thread_timings(n_threads); std::vector> thread_scores(n_threads); + std::vector thread_constraint_discards(n_threads, 0); // Spawn worker threads std::vector workers; @@ -368,6 +379,7 @@ DrivenResult parallel_driven_search( for (int t = 0; t < n_threads; ++t) { ctx.thread_timings = thread_timings.data(); ctx.thread_scores = thread_scores.data(); + ctx.thread_constraint_discards = thread_constraint_discards.data(); ctx.thread_id = t; workers.emplace_back(worker_thread, ctx); } @@ -579,6 +591,7 @@ DrivenResult parallel_driven_search( // Sum per-thread timings; merge per-thread replicate scores for (int t = 0; t < n_threads; ++t) { result.timings += thread_timings[t]; + result.constraint_discards += thread_constraint_discards[t]; for (double s : thread_scores[t]) { result.replicate_scores.push_back(s); } diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index fe9abee38..57b7959c3 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2105,6 +2105,16 @@ List ts_driven_search( result = ts::driven_search(pool, ds, params, cd_ptr); } + // Reported here rather than where it is detected: the count accumulates on + // worker threads, and Rf_warning() is a main-thread-only call. + if (result.constraint_discards > 0) { + Rf_warning( + "%d replicate(s) ended on a tree that could not be made to satisfy " + "`constraint`, and were discarded. The remaining trees do satisfy it; " + "raise `maxReplicates` if too few trees were found.", + result.constraint_discards); + } + // Build timings as a NumericVector (lighter than List) NumericVector timings = NumericVector::create( Named("wagner_ms") = result.timings.wagner_ms, @@ -2241,7 +2251,8 @@ List ts_collapse_pool( List scoringConfig, Nullable hsjConfig = R_NilValue, Nullable xformConfig = R_NilValue, - Nullable consSplitMatrix = R_NilValue) + Nullable consSplitMatrix = R_NilValue, + Nullable consZero = R_NilValue) { ts::DataSet ds = unpack_scoring(contrast, tip_data, weight, levels, scoringConfig); @@ -2256,16 +2267,31 @@ List ts_collapse_pool( // 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. + // + // 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. 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(wps, 0); - for (int c = 0; c < n_tip && c < cs.ncol(); ++c) { - if (cs(r, c)) b[c >> 6] |= (1ULL << (c & 63)); - } + 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; @@ -2273,6 +2299,24 @@ List ts_collapse_pool( } cons_canon.push_back(std::move(b)); } + 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: + // such a split is realised by a terminal edge, never a collapse candidate. + std::vector n_one_tips(cons_one.size(), 0); + std::vector n_zero_tips(cons_one.size(), 0); + for (size_t r = 0; r < cons_one.size(); ++r) { + for (int w = 0; w < wps; ++w) { + n_one_tips[r] += ts::popcount64(cons_one[r][w]); + n_zero_tips[r] += ts::popcount64(cons_zero[r][w]); + } } std::vector reps; // representative collapsed edges @@ -2353,6 +2397,48 @@ List ts_collapse_pool( 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. + // 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. + // + // Where none survives, the MRCA of a group is the node protected: the + // postorder visits every node before its parent, so the first node to + // hold a whole group is its MRCA, and keeping that one edge suffices, + // since contracting an edge below it leaves its descendant set — and so + // the split it displays — unchanged. + for (size_t r = 0; r < cons_one.size(); ++r) { + if (n_one_tips[r] < 2 || n_zero_tips[r] < 2) continue; + const std::vector* grp[2] = { &cons_one[r], &cons_zero[r] }; + + 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]; + 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 (!flags[v]) { survives = true; break; } + if (to_protect < 0) to_protect = v; // the MRCA, in postorder + } + } + if (!survives && to_protect >= 0) flags[to_protect] = 0; + } } // Dedup on the collapsed split set (skips the flagged zero-length edges). diff --git a/tests/testthat/test-MaximizeParsimony-features.R b/tests/testthat/test-MaximizeParsimony-features.R index cfd494f6f..a273eabca 100644 --- a/tests/testthat/test-MaximizeParsimony-features.R +++ b/tests/testthat/test-MaximizeParsimony-features.R @@ -535,7 +535,9 @@ test_that("Constrained Wagner tree works with multiple seeds", { list(contrast = at$contrast, tip_data = matrix(unlist(ds5, use.names = FALSE), nrow = 5, byrow = TRUE), weight = at$weight, levels = at$levels), - consArgs)) + # The flat kernels declare the constraint arguments as formals, so they + # take only the subset .PrepareConstraint() builds for them. + TreeSearch:::.KernelConstraintArgs(consArgs))) expect_true(is.finite(result$score), info = paste("seed", s)) expect_equal(nrow(result$edge), 8L, info = paste("seed", s)) } diff --git a/tests/testthat/test-ts-constraint-holes.R b/tests/testthat/test-ts-constraint-holes.R new file mode 100644 index 000000000..797ad468a --- /dev/null +++ b/tests/testthat/test-ts-constraint-holes.R @@ -0,0 +1,267 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +## A `constraint` must bind the trees the caller is handed, at each of the three +## boundaries where it can be lost: the starting tree, the pool capture, and the +## final collapse (T-402, T-324, T-403). +## +## Each test asserts COMPLIANCE of the returned trees, not the score alone. A +## constraint-violating tree is drawn from a wider set of topologies than a legal +## one, so it scores better; a score assertion alone would pass on exactly the +## tree that breaks the contract. + +library("TreeTools", quietly = TRUE) + +taxa <- letters[1:8] + +# Does `tr` display a split with all of `one` on one side and all of `zero` on +# the other? Spelled out from the edge matrix rather than via `Splits`: `%in%` +# on a Splits object dispatches differently under test_check() than under +# load_all(), and this has to answer the same way in both. +ConstraintShown <- function(tr, one, zero) { + tr <- Postorder(tr) + edge <- tr[["edge"]] + label <- tr[["tip.label"]] + nOne <- integer(max(edge)) + nZero <- integer(max(edge)) + nOne[match(one, label)] <- 1L + nZero[match(zero, label)] <- 1L + for (i in seq_len(nrow(edge))) { + nOne[edge[i, 1]] <- nOne[edge[i, 1]] + nOne[edge[i, 2]] + nZero[edge[i, 1]] <- nZero[edge[i, 1]] + nZero[edge[i, 2]] + } + # Postorder lists every node before its parent, so the first node holding a + # whole group is that group's MRCA. + nodes <- c(edge[, 2], edge[nrow(edge), 1]) + mrcaOne <- nodes[nOne[nodes] == length(one)][1] + mrcaZero <- nodes[nZero[nodes] == length(zero)][1] + (!is.na(mrcaOne) && nZero[mrcaOne] == 0L) || + (!is.na(mrcaZero) && nOne[mrcaZero] == 0L) +} + +AllShown <- function(trees, one, zero) { + sum(vapply(trees, ConstraintShown, logical(1), one, zero)) +} + +# Two characters supporting (a, e) and two supporting (b, f): the unconstrained +# optimum groups a with e and b with f, which no tree holding {a, b} together +# can do. Optimum 6 unconstrained, 10 under an {a, b} constraint. +abDataset <- local({ + m <- rbind( + c(1, 0, 0, 0, 1, 0, 0, 0), + c(1, 0, 0, 0, 1, 0, 0, 0), + c(0, 1, 0, 0, 0, 1, 0, 0), + c(0, 1, 0, 0, 0, 1, 0, 0), + c(0, 0, 0, 0, 0, 0, 1, 1), + c(0, 0, 0, 0, 0, 0, 1, 1) + ) + colnames(m) <- taxa + MatrixToPhyDat(t(m)) +}) + +# {a, b} against every other taxon: no ambiguous tip, so the constraint the +# user states and the stricter one the search enforces internally coincide. +abConstraint <- MatrixToPhyDat(matrix( + c(1, 1, 0, 0, 0, 0, 0, 0), ncol = 1, dimnames = list(taxa, NULL) +)) + +# Scores 6 -- better than any {a, b}-compliant tree -- and violates {a, b}. +abViolatingStart <- ape::read.tree(text = "(((a,e),(b,f)),((c,d),(g,h)));") + + +test_that("a violating `tree` cannot beat the constrained optimum (T-402)", { + # Control: the same constraint from a cold start reaches 10 and complies. + set.seed(1) + cold <- MaximizeParsimony(abDataset, constraint = abConstraint, + maxReplicates = 4L, verbosity = 0L) + expect_equal(as.numeric(attr(cold, "score")), 10) + expect_equal(AllShown(cold, c("a", "b"), setdiff(taxa, c("a", "b"))), + length(cold)) + + # maxReplicates = 1: gating the pool capture alone would leave the pool + # empty here, and MaximizeParsimony() would fall back to returning the + # supplied start. The start must be dealt with at the boundary. + set.seed(1) + expect_warning( + one <- MaximizeParsimony(abDataset, tree = abViolatingStart, + constraint = abConstraint, maxReplicates = 1L, + verbosity = 0L), + "do not satisfy `constraint`" + ) + expect_equal(as.numeric(attr(one, "score")), 10) + expect_equal(AllShown(one, c("a", "b"), setdiff(taxa, c("a", "b"))), + length(one)) + + # Several replicates: an illegal score is better than any legal one, so it + # evicts every compliant tree the other replicates find. One bad start must + # not cost the whole search. + set.seed(1) + expect_warning( + many <- MaximizeParsimony(abDataset, tree = abViolatingStart, + constraint = abConstraint, maxReplicates = 8L, + verbosity = 0L), + "do not satisfy `constraint`" + ) + expect_equal(as.numeric(attr(many, "score")), 10) + expect_equal(AllShown(many, c("a", "b"), setdiff(taxa, c("a", "b"))), + length(many)) +}) + + +test_that("no tree in the returned pool breaks the constraint (T-324)", { + # What this asserts is the outcome -- every tree handed back complies -- over + # the whole pool, not just the best-score trees: `poolSuboptimal` retains the + # near-misses, which is where an ungated capture shows up. It does NOT prove + # the capture gate itself fires; the route that motivated T-324 is Wagner + # retry-exhaustion, whose reachability is unconfirmed and which cannot be + # forced from R. Treat this as a contract test, not a gate test. + set.seed(2) + expect_warning( + result <- MaximizeParsimony(abDataset, tree = abViolatingStart, + constraint = abConstraint, maxReplicates = 3L, + verbosity = 0L, collapse = FALSE, + poolSuboptimal = 4), + "do not satisfy `constraint`" + ) + expect_equal(AllShown(result, c("a", "b"), setdiff(taxa, c("a", "b"))), + length(result)) + expect_gte(as.numeric(attr(result, "score")), 10) + + # The parallel driver has its own copy of the capture, on a per-thread + # constraint and pool; two threads is the project's per-agent core limit. + set.seed(2) + expect_warning( + parallel <- MaximizeParsimony(abDataset, tree = abViolatingStart, + constraint = abConstraint, maxReplicates = 4L, + nThreads = 2L, verbosity = 0L), + "do not satisfy `constraint`" + ) + expect_equal(AllShown(parallel, c("a", "b"), setdiff(taxa, c("a", "b"))), + length(parallel)) + expect_equal(as.numeric(attr(parallel, "score")), 10) +}) + + +test_that("every flat kernel takes .PrepareConstraint()'s output", { + # The flat `ts_*` kernels declare their constraint arguments as formals, so a + # field .PrepareConstraint() adds for the list-config entry points is an + # unused-argument error at any site that splats the whole list into one. + # Assert the filter covers every formal each kernel actually declares, and + # exercise the entry points that splat -- `Resample(nReplicates > 1)` had no + # constrained coverage at all, so an unfiltered splat there stayed green. + kernels <- list(TreeSearch:::ts_wagner_tree, + TreeSearch:::ts_random_wagner_tree, + TreeSearch:::ts_resample_search, + TreeSearch:::ts_parallel_resample, + TreeSearch:::ts_successive_approx) + consArgs <- TreeSearch:::.PrepareConstraint(abConstraint, abDataset) + filtered <- names(TreeSearch:::.KernelConstraintArgs(consArgs)) + for (k in kernels) { + kernelFormals <- names(formals(k)) + # Nothing the kernel does not declare -- an unused-argument error... + expect_equal(setdiff(filtered, kernelFormals), character(0)) + # ...and nothing it declares left behind, which would silently fall back to + # the kernel's own default instead of the constraint the caller gave. + expect_equal(setdiff(intersect(names(consArgs), kernelFormals), filtered), + character(0)) + } + + set.seed(4) + expect_s3_class( + Resample(abDataset, constraint = abConstraint, nReplicates = 2L, + maxReplicates = 2L), + "multiPhylo" + ) + set.seed(4) + expect_s3_class(AdditionTree(abDataset, constraint = abConstraint), "phylo") +}) + + +test_that("a three-state constraint says what it does and does not enforce", { + # Constraints are enforced as bipartitions throughout -- locked-node filter, + # capture gate, and impose_constraint(), which can repair to nothing else. A + # third state's taxa are therefore unconstrained, and a character can sit + # above its minimum length with the enforced split intact. What must not + # happen is that silently: judging captures by the stricter full-Fitch reading + # instead would discard every replicate of a search that cannot do better, + # turning a partial answer into no answer at all. + taxa6 <- letters[1:6] + constraint <- MatrixToPhyDat(matrix( + c("2", "2", "0", "0", "1", "1"), ncol = 1, + dimnames = list(taxa6, NULL) + )) + expect_equal(as.numeric(MinimumLength(constraint)), 2) + # {a, b} is an exact clade here -- the enforced split holds -- yet the + # character costs 3, which is the gap the warning is about. + expect_equal(as.numeric(TreeLength( + ape::read.tree(text = "((a,b),(e,(c,(f,d))));"), constraint)), 3) + + expect_warning(TreeSearch:::.PrepareConstraint(constraint, constraint), + "more than two states") + + # Data pulling against the constraint: it supports (a,b) but also (c,e) and + # (d,f), which splits the intermediate state apart. + m <- rbind( + c(1, 1, 0, 0, 0, 0), c(1, 1, 0, 0, 0, 0), + c(0, 0, 1, 0, 1, 0), c(0, 0, 1, 0, 1, 0), + c(0, 0, 0, 1, 0, 1), c(0, 0, 0, 1, 0, 1) + ) + colnames(m) <- taxa6 + dataset <- MatrixToPhyDat(t(m)) + + # collapse = FALSE keeps the trees binary for TreeLength(). + set.seed(7) + expect_warning( + result <- MaximizeParsimony(dataset, constraint = constraint, + maxReplicates = 4L, verbosity = 0L, + collapse = FALSE), + "more than two states" + ) + # The enforced split still binds on every returned tree... + expect_equal(AllShown(result, c("a", "b"), c("c", "d")), length(result)) + # ...and a search that can only partially satisfy the constraint returns + # trees rather than erroring with an empty pool. + expect_gt(length(result), 0) +}) + + +test_that("collapse keeps the constraint visible (T-403)", { + # Only (a, e) and (b, f) are supported, so the branch that separates + # {a, b} from {c, d} is unsupported and collapses -- taking the constraint + # with it. The node realising the split is {a, e, b, f}, not the "1" group + # {a, b}, so protection keyed on an exact match with the "1" group does not + # reach it. + m <- rbind( + c(1, 0, 0, 0, 1, 0, 0, 0), + c(1, 0, 0, 0, 1, 0, 0, 0), + c(0, 1, 0, 0, 0, 1, 0, 0), + c(0, 1, 0, 0, 0, 1, 0, 0) + ) + colnames(m) <- taxa + dataset <- MatrixToPhyDat(t(m)) + + # e--h ambiguous: the constraint asks only that {a, b} be separated from + # {c, d}, which the start below already does. + constraint <- MatrixToPhyDat(matrix( + c("1", "1", "0", "0", "?", "?", "?", "?"), + ncol = 1, dimnames = list(taxa, NULL) + )) + start <- ape::read.tree(text = "(((a,e),(b,f)),(c,(d,(g,h))));") + + set.seed(1) + collapsed <- MaximizeParsimony(dataset, tree = start, + constraint = constraint, maxReplicates = 2L, + verbosity = 0L) + expect_equal(AllShown(collapsed, c("a", "b"), c("c", "d")), + length(collapsed)) + + # Built-in control: without collapsing, the split was never at risk. + set.seed(1) + resolved <- MaximizeParsimony(dataset, tree = start, + constraint = constraint, maxReplicates = 2L, + verbosity = 0L, collapse = FALSE) + expect_equal(AllShown(resolved, c("a", "b"), c("c", "d")), + length(resolved)) + expect_equal(as.numeric(attr(collapsed, "score")), + as.numeric(attr(resolved, "score"))) +}) diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index a87d2b292..605633638 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -217,6 +217,38 @@ Per-strategy attempt and success counts are returned in the `strategy_diagnostics` attribute of the search result for post-hoc inspection. +### Starting trees under a constraint + +Whatever its source -- a Wagner build, a random topology, or a tree supplied +through `tree` -- a start that does not satisfy the constraint is repaired +before the replicate scores it, by the same minimal-SPR routine that repairs +fused trees. +This has to happen at the start rather than by rejecting the tree afterwards. +Constrained rearrangement cannot climb out of a violating tree: an unmapped +constraint split makes every candidate regraft illegal, so the replicate +freezes on the tree it began with. +Nor can the violation be caught by comparison later, because a violating tree +is drawn from a wider set of topologies than a legal one and so tends to score +*better*; taken as a baseline, it makes the legal repair look like a +regression, and no subsequent phase can accept it. +Where the repair does not succeed, the start is discarded in favour of a fresh +constrained Wagner build. +The same check gates each replicate's finished tree on its way into the pool, +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. + ## The driven search pipeline