From 4e6adae593b6dd4821f63df7285cb44b39aa51e6 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 17:24:26 +0100 Subject: [PATCH 1/3] fix: enforce user constraints at the start-tree, pool-capture and collapse boundaries Three independent holes let a `constraint` stop binding the trees the caller is handed. A start tree supplied through `tree =` was never checked against the constraint (T-402). Constrained rearrangement cannot climb out of a violating tree -- an unmapped split makes every candidate regraft illegal -- so the replicate froze on it and reported its unconstrained score, which then evicted the compliant trees the other replicates found. The repair has to happen before anything takes the start's score as a baseline: a violating tree is drawn from a wider set of topologies and so scores better, making the legal repair look like a regression to any later accept test. `run_single_replicate` now repairs any violating start, whatever its source, and falls back to a constrained Wagner build where the heuristic repair does not take. `MaximizeParsimony()` warns when it was the caller's `tree` that conflicted. The per-replicate pool capture had no constraint gate, asymmetrically to the fuse capture beside it (T-324). All three capture sites -- interrupted and normal in the serial driver, and the parallel worker's -- now go through `capture_satisfies_constraint()`, which tests with `violates_constraint_posthoc` rather than the fuse's `constraint_node < 0`: a tree can map every constraint node and still fail the full Fitch check. Discards are counted and reported from the main thread. With the pool gated, an empty pool under a constraint is an error rather than a fall back to the unvalidated start tree. The collapse pass protected only a node whose tip set was the 1 group exactly (T-403). Tips coded `?` for a constraint character are free to sit on either side, so the split is often realised by a wider node -- left collapsible, and contracted away under the default `collapse = TRUE`. `consZero` is now plumbed through `.PrepareConstraint()` to the kernel, which additionally protects the MRCA of either group when it excludes the other. Co-Authored-By: Claude Opus 5 --- R/AdditionTree.R | 4 +- R/MaximizeParsimony.R | 94 ++++++++- R/RcppExports.R | 4 +- R/Resample.R | 4 +- R/SuccessiveApproximations.R | 4 +- man/MaximizeParsimony.Rd | 6 + src/RcppExports.cpp | 7 +- src/TreeSearch-init.c | 4 +- src/ts_driven.cpp | 61 +++++- src/ts_driven.h | 13 ++ src/ts_parallel.cpp | 22 ++- src/ts_rcpp.cpp | 90 ++++++++- .../test-MaximizeParsimony-features.R | 4 +- tests/testthat/test-ts-constraint-holes.R | 178 ++++++++++++++++++ vignettes/search-algorithm.Rmd | 21 +++ 15 files changed, 491 insertions(+), 25 deletions(-) create mode 100644 tests/testthat/test-ts-constraint-holes.R diff --git a/R/AdditionTree.R b/R/AdditionTree.R index 6b9ce7877..498657b80 100644 --- a/R/AdditionTree.R +++ b/R/AdditionTree.R @@ -142,7 +142,9 @@ AdditionTree <- function(dataset, concavity = Inf, constraint, sequence) { addition_order = addition_order, 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..752b5a48c 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -190,6 +190,7 @@ list( consSplitMatrix = consSplits, + consZero = consZero, consContrast = consContrast, consTipData = consTipData, consWeight = as.integer(consWeight), @@ -198,6 +199,57 @@ ) } +# 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. +.kernelConstraintArgs <- c("consSplitMatrix", "consContrast", "consTipData", + "consWeight", "consLevels", "consExpectedScore") + +.KernelConstraintArgs <- function(consArgs) { + consArgs[intersect(names(consArgs), .kernelConstraintArgs)] +} + +# 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. +.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: columns 1..nRow are the + # "1" groups, the rest the "0" groups. + counts <- matrix(0L, nrow = max(edge), ncol = 2L * nRow) + counts[seq_len(nTip), ] <- t(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 +667,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 @@ -1520,6 +1578,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.", + call. = FALSE) + } + } + # --- Profile parsimony: extract info_amounts --- profileArgs <- list() if (useProfile) { @@ -1654,12 +1729,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 +1763,15 @@ MaximizeParsimony <- function( }) } if (length(outTrees) == 0L) { + # `treeTpl` is a starting tree, which under a constraint is exactly what + # may not be handed back: an empty pool means no replicate produced a tree + # the constraint gate accepted (or none finished at all), and returning an + # unvalidated tree would break the guarantee `constraint` makes. + if (!is.null(constraintConfig)) { + 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..58b89ea6d 100644 --- a/R/Resample.R +++ b/R/Resample.R @@ -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/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 917751868..40ec6bca1 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 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..5f5a3f8fd 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -51,6 +51,25 @@ ProgressInfo make_progress(int rep, const DrivenParams& params, } // 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 (!violates_constraint_posthoc(tree, *cd)) return true; + + impose_constraint(tree, *cd); + tree.build_postorder(); + tree.reset_states(ds); + score = score_tree(tree, ds); + return !violates_constraint_posthoc(tree, *cd); +} + // --- Single-replicate pipeline --- ReplicateResult run_single_replicate( @@ -151,6 +170,30 @@ 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 && + violates_constraint_posthoc(result.tree, *cd)) { + impose_constraint(result.tree, *cd); + result.tree.build_postorder(); + result.tree.reset_states(ds); + if (violates_constraint_posthoc(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. + random_wagner_tree(result.tree, ds, 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,14 +1098,24 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, result.timings += rep_result.timings; + // 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; + // 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_ok) { + compute_collapsed_flags(rep_result.tree, ds, rep_collapsed); + } if (rep_result.interrupted) { - if (rep_result.score < 1e18) { + if (rep_ok && rep_result.score < 1e18) { pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); } result.timed_out = true; @@ -1071,7 +1124,9 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, // Add to pool with collapsed-topology dedup double prev_best = pool.best_score(); - pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); + if (rep_ok) { + pool.add_collapsed(rep_result.tree, rep_result.score, rep_collapsed); + } bool score_improved = pool.best_score() < prev_best; if (score_improved) { result.last_improved_rep = rep1; 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..f9b39b644 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,11 +233,17 @@ 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); + // 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); + } else { + ++ctx.thread_constraint_discards[ctx.thread_id]; + } // Record per-replicate score for Chao1 coverage estimation ctx.thread_scores[ctx.thread_id].push_back(rep_result.score); @@ -361,6 +370,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 +378,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 +590,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 fcb6678c0..a4199ade4 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2022,6 +2022,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, @@ -2158,7 +2168,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); @@ -2173,16 +2184,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 — and that node, being + // unmatched, was left collapsible, contracting the enforced grouping away. + // `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; @@ -2190,6 +2216,13 @@ 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)); } std::vector reps; // representative collapsed edges @@ -2270,6 +2303,53 @@ List ts_collapse_pool( if (eq) { flags[v] = 0; break; } } } + + // Protect the node that realises each split under the looser, promised + // reading: the MRCA of one group, when it holds none of the other. The + // postorder visits every node before its parent, so the first node to + // hold a whole group is its MRCA; keeping that one edge is enough, + // because contracting an edge below it leaves its descendant set — and so + // the split it displays — unchanged. Which of the two groups is the + // clade depends on the rooting alone, so try each in turn. Groups of + // fewer than two tips are skipped: such a split is realised by a terminal + // edge, which is never a collapse candidate. + for (size_t r = 0; r < cons_one.size(); ++r) { + const std::vector* grp[2] = { &cons_one[r], &cons_zero[r] }; + int n_in_group[2] = {0, 0}; + for (int side = 0; side < 2; ++side) { + for (int w = 0; w < wps; ++w) { + n_in_group[side] += ts::popcount64((*grp[side])[w]); + } + } + if (n_in_group[0] < 2 || n_in_group[1] < 2) continue; + for (int side = 0; side < 2; ++side) { + const std::vector& in = *grp[side]; + const std::vector& out = *grp[1 - side]; + // postorder holds internal nodes only, and the MRCA of two or more + // tips is internal, so the first match is that MRCA. + int mrca = -1; + for (size_t pi = 0; pi < tree.postorder.size() && mrca < 0; ++pi) { + const int node = tree.postorder[pi]; + const uint64_t* nb = &tb[static_cast(node) * wps]; + bool holds = true; + for (int w = 0; w < wps; ++w) { + if ((nb[w] & in[w]) != in[w]) { holds = false; break; } + } + if (holds) mrca = node; + } + if (mrca < 0) continue; + const uint64_t* mb = &tb[static_cast(mrca) * wps]; + bool clean = true; + for (int w = 0; w < wps; ++w) { + if (mb[w] & out[w]) { clean = false; break; } + } + if (!clean) continue; // this side is not the clade + if (mrca > n_tip && mrca < static_cast(flags.size())) { + flags[mrca] = 0; + } + break; + } + } } // 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..2cec1f36f --- /dev/null +++ b/tests/testthat/test-ts-constraint-holes.R @@ -0,0 +1,178 @@ +# Tier 2: skipped on CRAN; see tests/testing-strategy.md +skip_on_cran() + +## Three holes through which a `constraint` stopped binding the trees the user +## is handed (T-402, T-324, T-403). +## +## T-402: a start tree supplied via `tree =` was never checked against the +## constraint. Constrained rearrangement cannot repair such a tree — an +## unmapped split makes regraft_violates_constraint() reject every move — so +## the replicate froze on it and reported its unconstrained score, which then +## evicted the compliant trees other replicates found. +## +## T-324: the per-replicate pool capture had no constraint gate at all, +## asymmetrically to the fuse capture beside it, so any violating tree that +## reached it was handed straight back. +## +## T-403: the collapse pass protected only a node whose tip set was the "1" +## group EXACTLY. Tips ambiguous for a constraint character are free to sit on +## either side, so the split is often realised by a wider node — left +## unprotected, and contracted away under the default `collapse = TRUE`. + +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: the violating tree's illegal score used to evict every + # compliant tree the other replicates found. + 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("a violating tree never enters the pool (T-324)", { + # The Wagner retry-exhaustion route that motivated T-324 is not constructible + # on demand, so the shared downstream half -- the ungated pool capture -- is + # driven through T-402's start instead: without the gate the replicate's + # frozen, violating tree is captured verbatim. `poolSuboptimal` keeps + # non-best trees too, so a violating tree would be visible even if a better + # compliant one existed. + 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) +}) + + +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}, which is why an exact-match protection missed 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..2346b56ed 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -217,6 +217,27 @@ 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. + ## The driven search pipeline From 8cfbf3b7923a2640a22f15b4c70553cb0fc6606c Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Tue, 4 Aug 2026 18:18:37 +0100 Subject: [PATCH 2/3] fix: address review findings on the constraint gates `Resample(constraint =, nReplicates > 1)` errored: `R/Resample.R:372` splats `.PrepareConstraint()`'s whole list into `ts_parallel_resample()`, which has no `consZero` formal. It was the one splat site of four not filtered, and no test covered a constrained `Resample()` at all; a test now asserts the filter against every flat kernel's formals. The collapse protection over-reached. It protected the MRCA of a group unconditionally, so an unsupported branch that merely happened to be the narrowest node containing the 0 group was returned resolved -- the "unsupported non-constraint branches still collapse" half of the promise, and a shift in `n_topologies`. It now protects only when no realising edge survives the contraction on its own. A discarded replicate no longer feeds the strategy bandit, the Chao1 coverage scores or the replicate report: a violating tree scores better than any legal one, so those would be credited to whatever produced it. Its stopping rules still run -- skipping them would outlive the deadline and swallow an interrupt. An interrupted replicate is now tested but not repaired, since `impose_constraint()` has no interrupt check of its own. The constraint check short-circuits on the locked-node mapping, which is cheaper than the post-hoc Fitch check and strictly stronger, so only an unmapped split pays for the latter. Constrained wall-clock on Vinther2008 (30 paired seeds) goes from a 1.12 median ratio, 11 seeds >10% slower, to 1.003 with none >10%. Co-Authored-By: Claude Opus 5 --- NEWS.md | 19 +++++ R/MaximizeParsimony.R | 23 +++--- R/Resample.R | 2 +- src/ts_driven.cpp | 84 +++++++++++++++------- src/ts_parallel.cpp | 7 +- src/ts_rcpp.cpp | 69 +++++++++--------- tests/testthat/test-ts-constraint-holes.R | 88 ++++++++++++++++------- vignettes/search-algorithm.Rmd | 11 +++ 8 files changed, 204 insertions(+), 99 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3492cf2d4..7c8a12de7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,24 @@ # 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. + - `inapplicable = "xform"` scores are now reported at a canonical rooting, so a reported score is reproducible. The x-transformation's step matrix is asymmetric -- a gain costs one more than the number of secondary characters it diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 752b5a48c..a6740e3ca 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -203,11 +203,11 @@ # 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. -.kernelConstraintArgs <- c("consSplitMatrix", "consContrast", "consTipData", - "consWeight", "consLevels", "consExpectedScore") +.kernelConsFields <- c("consSplitMatrix", "consContrast", "consTipData", + "consWeight", "consLevels", "consExpectedScore") .KernelConstraintArgs <- function(consArgs) { - consArgs[intersect(names(consArgs), .kernelConstraintArgs)] + consArgs[intersect(names(consArgs), .kernelConsFields)] } # Does `tree` display a split separating a constraint character's "1" group @@ -1590,8 +1590,8 @@ MaximizeParsimony <- function( 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.", - call. = FALSE) + "they will be rearranged to comply before the search starts, ", + "or replaced if that fails.", call. = FALSE) } } @@ -1763,11 +1763,14 @@ MaximizeParsimony <- function( }) } if (length(outTrees) == 0L) { - # `treeTpl` is a starting tree, which under a constraint is exactly what - # may not be handed back: an empty pool means no replicate produced a tree - # the constraint gate accepted (or none finished at all), and returning an - # unvalidated tree would break the guarantee `constraint` makes. - if (!is.null(constraintConfig)) { + # `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`.") diff --git a/R/Resample.R b/R/Resample.R index 58b89ea6d..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) diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index 5f5a3f8fd..5e6561dd3 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -49,6 +49,26 @@ 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. 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. +bool constraint_satisfied(TreeState& tree, ConstraintData& cd) { + map_constraint_nodes(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, @@ -61,13 +81,13 @@ bool capture_satisfies_constraint(TreeState& tree, ConstraintData* cd, // 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 (!violates_constraint_posthoc(tree, *cd)) 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 !violates_constraint_posthoc(tree, *cd); + return constraint_satisfied(tree, *cd); } // --- Single-replicate pipeline --- @@ -179,11 +199,11 @@ ReplicateResult run_single_replicate( // 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 && - violates_constraint_posthoc(result.tree, *cd)) { + !constraint_satisfied(result.tree, *cd)) { impose_constraint(result.tree, *cd); result.tree.build_postorder(); result.tree.reset_states(ds); - if (violates_constraint_posthoc(result.tree, *cd)) { + 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. @@ -1098,6 +1118,21 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, result.timings += rep_result.timings; + if (rep_result.interrupted) { + // 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; + } + // 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 @@ -1106,28 +1141,25 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, rep_result.score); if (!rep_ok) ++result.constraint_discards; - // 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; + // 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); - } - if (rep_result.interrupted) { - if (rep_ok && rep_result.score < 1e18) { - 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(); - if (rep_ok) { + // 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; } - bool score_improved = pool.best_score() < prev_best; if (score_improved) { result.last_improved_rep = rep1; unsuccessful_reps = 0; @@ -1140,7 +1172,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); } @@ -1153,10 +1185,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_parallel.cpp b/src/ts_parallel.cpp index f9b39b644..4035d0f30 100644 --- a/src/ts_parallel.cpp +++ b/src/ts_parallel.cpp @@ -241,13 +241,14 @@ void worker_thread(WorkerContext ctx) { 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]; } - // Record per-replicate score for Chao1 coverage estimation - ctx.thread_scores[ctx.thread_id].push_back(rep_result.score); - ctx.replicates_done->fetch_add(1, std::memory_order_relaxed); // Check convergence diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index a4199ade4..923af55d9 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2189,10 +2189,10 @@ List ts_collapse_pool( // 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 — and that node, being - // unmatched, was left collapsible, contracting the enforced grouping away. - // `cons_one` / `cons_zero` are the raw (uncanonicalised) groups, from which - // the realising node is found per tree below. + // 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; @@ -2304,15 +2304,22 @@ List ts_collapse_pool( } } - // Protect the node that realises each split under the looser, promised - // reading: the MRCA of one group, when it holds none of the other. The + // 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; keeping that one edge is enough, - // because contracting an edge below it leaves its descendant set — and so - // the split it displays — unchanged. Which of the two groups is the - // clade depends on the rooting alone, so try each in turn. Groups of - // fewer than two tips are skipped: such a split is realised by a terminal - // edge, which is never a collapse candidate. + // 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. Groups of fewer than two tips are + // skipped: such a split is realised by a terminal edge, never a collapse + // candidate. for (size_t r = 0; r < cons_one.size(); ++r) { const std::vector* grp[2] = { &cons_one[r], &cons_zero[r] }; int n_in_group[2] = {0, 0}; @@ -2322,33 +2329,29 @@ List ts_collapse_pool( } } if (n_in_group[0] < 2 || n_in_group[1] < 2) continue; - for (int side = 0; side < 2; ++side) { + + 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]; - // postorder holds internal nodes only, and the MRCA of two or more - // tips is internal, so the first match is that MRCA. - int mrca = -1; - for (size_t pi = 0; pi < tree.postorder.size() && mrca < 0; ++pi) { - const int node = tree.postorder[pi]; - const uint64_t* nb = &tb[static_cast(node) * wps]; - bool holds = true; + 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]) { holds = false; break; } + if ((nb[w] & in[w]) != in[w] || (nb[w] & out[w])) { + realises = false; + break; + } } - if (holds) mrca = node; - } - if (mrca < 0) continue; - const uint64_t* mb = &tb[static_cast(mrca) * wps]; - bool clean = true; - for (int w = 0; w < wps; ++w) { - if (mb[w] & out[w]) { clean = false; break; } - } - if (!clean) continue; // this side is not the clade - if (mrca > n_tip && mrca < static_cast(flags.size())) { - flags[mrca] = 0; + if (!realises) continue; + if (!flags[v]) { survives = true; break; } + if (to_protect < 0) to_protect = v; // the MRCA, in postorder } - break; } + if (!survives && to_protect >= 0) flags[to_protect] = 0; } } diff --git a/tests/testthat/test-ts-constraint-holes.R b/tests/testthat/test-ts-constraint-holes.R index 2cec1f36f..b81ecb039 100644 --- a/tests/testthat/test-ts-constraint-holes.R +++ b/tests/testthat/test-ts-constraint-holes.R @@ -1,23 +1,14 @@ # Tier 2: skipped on CRAN; see tests/testing-strategy.md skip_on_cran() -## Three holes through which a `constraint` stopped binding the trees the user -## is handed (T-402, T-324, T-403). +## 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). ## -## T-402: a start tree supplied via `tree =` was never checked against the -## constraint. Constrained rearrangement cannot repair such a tree — an -## unmapped split makes regraft_violates_constraint() reject every move — so -## the replicate froze on it and reported its unconstrained score, which then -## evicted the compliant trees other replicates found. -## -## T-324: the per-replicate pool capture had no constraint gate at all, -## asymmetrically to the fuse capture beside it, so any violating tree that -## reached it was handed straight back. -## -## T-403: the collapse pass protected only a node whose tip set was the "1" -## group EXACTLY. Tips ambiguous for a constraint character are free to sit on -## either side, so the split is often realised by a wider node — left -## unprotected, and contracted away under the default `collapse = TRUE`. +## 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) @@ -101,8 +92,9 @@ test_that("a violating `tree` cannot beat the constrained optimum (T-402)", { expect_equal(AllShown(one, c("a", "b"), setdiff(taxa, c("a", "b"))), length(one)) - # Several replicates: the violating tree's illegal score used to evict every - # compliant tree the other replicates found. + # 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, @@ -116,13 +108,13 @@ test_that("a violating `tree` cannot beat the constrained optimum (T-402)", { }) -test_that("a violating tree never enters the pool (T-324)", { - # The Wagner retry-exhaustion route that motivated T-324 is not constructible - # on demand, so the shared downstream half -- the ungated pool capture -- is - # driven through T-402's start instead: without the gate the replicate's - # frozen, violating tree is captured verbatim. `poolSuboptimal` keeps - # non-best trees too, so a violating tree would be visible even if a better - # compliant one existed. +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, @@ -134,6 +126,49 @@ test_that("a violating tree never enters the pool (T-324)", { 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) + filtered <- names(TreeSearch:::.KernelConstraintArgs( + TreeSearch:::.PrepareConstraint(abConstraint, abDataset) + )) + for (k in kernels) { + expect_true(all(filtered %in% names(formals(k)))) + } + + 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") }) @@ -141,7 +176,8 @@ 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}, which is why an exact-match protection missed it. + # {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), diff --git a/vignettes/search-algorithm.Rmd b/vignettes/search-algorithm.Rmd index 2346b56ed..605633638 100644 --- a/vignettes/search-algorithm.Rmd +++ b/vignettes/search-algorithm.Rmd @@ -238,6 +238,17 @@ 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 From 239762070cded342dcb4229ee8843bfa3223c9c3 Mon Sep 17 00:00:00 2001 From: R script <1695515+ms609@users.noreply.github.com> Date: Wed, 5 Aug 2026 03:59:58 +0100 Subject: [PATCH 3/3] fix: make the constraint standard one thing, and say what it is Self-review of the two commits before it. A constraint is enforced as a split throughout -- the locked-node filter screens rearrangements on it, and impose_constraint() repairs to it and nothing else -- so a character with a third state has taxa nothing constrains. The capture gate was judging that same tree by the stricter full-Fitch reading, which is not a standard the search can reach: a probe on a three-state constraint discarded all four replicates and errored with an empty pool, where a partial answer existed. The gate now follows the mapping, and `.PrepareConstraint()` warns at input that an intermediate state is unconstrained, rather than leaving the caller to infer from `@param constraint` that it is not. `constraint_satisfied()` refreshes the DFS timestamps alongside the node ids. map_constraint_nodes() alone left the two out of step, and spr_search() reads both without re-mapping, so the `sprFirst = TRUE` warmup could classify a regraft against this tree's nodes and another tree's timestamps. A start whose repair fails falls back to a Wagner build, which is exactly the constructor that can exhaust its reshuffles and return a violating tree; it is now repaired rather than trusted. Also: the collapse protection's group sizes are counted once instead of per tree; the kernel-formals test now asserts that no constraint field is dropped, not just that none is unknown; `.ConstraintViolated()` indexes nodes by column so its accumulation runs down a column-major matrix rather than across it. Constrained wall-clock on Vinther2008 (20 paired seeds) is 0.95 median against `cpp-search`, scores identical. Co-Authored-By: Claude Opus 5 --- R/MaximizeParsimony.R | 34 ++++++++++--- man/AdditionTree.Rd | 5 +- man/MaximizeParsimony.Rd | 5 +- man/Resample.Rd | 5 +- man/SuccessiveApproximations.Rd | 5 +- src/ts_driven.cpp | 37 +++++++++++--- src/ts_rcpp.cpp | 23 +++++---- tests/testthat/test-ts-constraint-holes.R | 61 +++++++++++++++++++++-- 8 files changed, 142 insertions(+), 33 deletions(-) diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index a6740e3ca..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" @@ -217,26 +228,32 @@ # 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: columns 1..nRow are the - # "1" groups, the rest the "0" groups. - counts <- matrix(0L, nrow = max(edge), ncol = 2L * nRow) - counts[seq_len(nTip), ] <- t(rbind(consOne, consZero)) + # 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], ] + 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] + one <- counts[r, ] + zero <- counts[nRow + r, ] nOne <- sum(consOne[r, ]) nZero <- sum(consZero[r, ]) mrcaOne <- nodes[one[nodes] == nOne][1] @@ -754,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. diff --git a/man/AdditionTree.Rd b/man/AdditionTree.Rd index 13fbc9e53..e4acb2403 100644 --- a/man/AdditionTree.Rd +++ b/man/AdditionTree.Rd @@ -36,7 +36,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 40ec6bca1..56b99d993 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -144,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/ts_driven.cpp b/src/ts_driven.cpp index 5e6561dd3..fae5ddfe3 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -53,14 +53,28 @@ ProgressInfo make_progress(int rep, const DrivenParams& params, // coded 1 for each constraint character from those coded 0? // // violates_constraint_posthoc() answers that directly, but builds a whole -// TreeState and scores it. 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. +// 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) { - map_constraint_nodes(tree, 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); @@ -206,10 +220,17 @@ ReplicateResult run_single_replicate( 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. + // 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); } diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 923af55d9..2daef0078 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -2224,6 +2224,17 @@ List ts_collapse_pool( } 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 std::vector rep_hash; // collapsed-split hash per rep @@ -2317,18 +2328,10 @@ List ts_collapse_pool( // 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. Groups of fewer than two tips are - // skipped: such a split is realised by a terminal edge, never a collapse - // candidate. + // 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] }; - int n_in_group[2] = {0, 0}; - for (int side = 0; side < 2; ++side) { - for (int w = 0; w < wps; ++w) { - n_in_group[side] += ts::popcount64((*grp[side])[w]); - } - } - if (n_in_group[0] < 2 || n_in_group[1] < 2) continue; bool survives = false; int to_protect = -1; diff --git a/tests/testthat/test-ts-constraint-holes.R b/tests/testthat/test-ts-constraint-holes.R index b81ecb039..797ad468a 100644 --- a/tests/testthat/test-ts-constraint-holes.R +++ b/tests/testthat/test-ts-constraint-holes.R @@ -154,11 +154,16 @@ test_that("every flat kernel takes .PrepareConstraint()'s output", { TreeSearch:::ts_resample_search, TreeSearch:::ts_parallel_resample, TreeSearch:::ts_successive_approx) - filtered <- names(TreeSearch:::.KernelConstraintArgs( - TreeSearch:::.PrepareConstraint(abConstraint, abDataset) - )) + consArgs <- TreeSearch:::.PrepareConstraint(abConstraint, abDataset) + filtered <- names(TreeSearch:::.KernelConstraintArgs(consArgs)) for (k in kernels) { - expect_true(all(filtered %in% names(formals(k)))) + 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) @@ -172,6 +177,54 @@ test_that("every flat kernel takes .PrepareConstraint()'s output", { }) +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