From 773eec1587a15f38ee9c8dc18f543c88770480d8 Mon Sep 17 00:00:00 2001 From: R script Date: Thu, 9 Jul 2026 15:11:54 +0100 Subject: [PATCH 01/12] feat(bremer): Bremer (decay) support via converse-constraint + pool Add standalone Bremer(tree, dataset, method=c("constraint","pool")) computing the per-clade decay index L(-C) - L*, plus SuboptimalTrees() and per-tree pool scores on MaximizeParsimony(collapse=FALSE) for landscape analysis. Rigorous default engine forces each clade absent via new C++ negative-constraint machinery (ConstraintData.neg_*, add_negative_constraint, displays_forbidden_clade) plumbed through consNegSplitMatrix + an internal .negativeConstraint arg. Soundness + optimality via three layers: a TBR regraft move-guard, a TreePool backstop (set_forbidden rejects any clade-displaying tree), and a not-C start repair. Drift and annealing accept worse moves through an unguarded path, so they are disabled in converse searches; ratchet/sectorial/nni-perturb re-optimise through the guarded TBR and stay on for search power. Oracle-validated against exhaustive enumeration of all binary trees, including the adversarial case (optimum displays C, not-C 3 steps away) and implied weights. Docs regenerated: adds Bremer/SuboptimalTrees to the split-support family and syncs the previously-stale SearchControl.Rd in-sector-drift params. Co-Authored-By: Claude Opus 4.8 --- NAMESPACE | 2 + NEWS.md | 14 ++ R/Bremer.R | 285 ++++++++++++++++++++++++++ R/MaximizeParsimony.R | 74 +++++++ R/SuboptimalTrees.R | 73 +++++++ inst/REFERENCES.bib | 22 ++ man/Bremer.Rd | 198 ++++++++++++++++++ man/ConcordanceTable.Rd | 4 +- man/JackLabels.Rd | 4 +- man/Jackknife.Rd | 4 +- man/MaximizeParsimony.Rd | 17 ++ man/Morphy.Rd | 4 +- man/MostContradictedFreq.Rd | 4 +- man/PaintCharacters.Rd | 4 +- man/PresCont.Rd | 4 +- man/SearchControl.Rd | 35 ++++ man/SiteConcordance.Rd | 4 +- man/SuboptimalTrees.Rd | 86 ++++++++ src/ts_constraint.cpp | 85 ++++++++ src/ts_constraint.h | 27 +++ src/ts_driven.cpp | 21 ++ src/ts_pool.cpp | 11 + src/ts_pool.h | 12 ++ src/ts_rcpp.cpp | 23 ++- src/ts_tbr.cpp | 19 ++ tests/testthat/test-Bremer.R | 213 +++++++++++++++++++ tests/testthat/test-SuboptimalTrees.R | 65 ++++++ 27 files changed, 1303 insertions(+), 11 deletions(-) create mode 100644 R/Bremer.R create mode 100644 R/SuboptimalTrees.R create mode 100644 man/Bremer.Rd create mode 100644 man/SuboptimalTrees.Rd create mode 100644 tests/testthat/test-Bremer.R create mode 100644 tests/testthat/test-SuboptimalTrees.R diff --git a/NAMESPACE b/NAMESPACE index 218352864..b46ae5f2a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -25,6 +25,7 @@ S3method(summary,morphyPtr) export(.NonDuplicateRoot) export(.UniqueExceptHits) export(AdditionTree) +export(Bremer) export(C_MorphyLength) export(Carter1) export(CharacterHierarchy) @@ -113,6 +114,7 @@ export(SharedPhylogeneticConcordance) export(SingleCharMorphy) export(StepInformation) export(StopUnlessBifurcating) +export(SuboptimalTrees) export(Suboptimality) export(SuccessiveApproximations) export(SuccessiveWeights) diff --git a/NEWS.md b/NEWS.md index 1ae2c3e99..5b34d2dc3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,19 @@ # To integrate into 2.0.0 notes +- New `Bremer()` calculates Bremer (decay) support for each clade of a reference + tree. The default `method = "constraint"` runs a rigorous converse-constraint + search per clade -- the shortest tree forced to *lack* the clade -- giving + trustworthy decay indices with bounded memory. A fast, approximate + `method = "pool"` instead reads a pool of suboptimal trees (an upper bound, + right-censored at the sampling depth). Both engines match an exhaustive + brute-force oracle on small datasets, including under implied weights. + +- New `SuboptimalTrees()` returns every tree the search retained within a given + number of steps of the optimum, each annotated with its parsimony score, for + landscape analysis. `MaximizeParsimony(collapse = FALSE)` now surfaces the + per-tree scores of the retained pool via a `scores` attribute (and a `score` + attribute on each tree), so `Suboptimality()` works on the result directly. + - `MaximizeParsimony()` now contracts zero-length (unsupported) branches into polytomies by default (`collapse = TRUE`), deduplicating the returned trees on the resulting collapsed topologies, à la TNT's "collapse zero-length diff --git a/R/Bremer.R b/R/Bremer.R new file mode 100644 index 000000000..22ca465ce --- /dev/null +++ b/R/Bremer.R @@ -0,0 +1,285 @@ +#' Bremer (decay) support +#' +#' `Bremer()` calculates the Bremer support (decay index) +#' \insertCite{Bremer1988,Bremer1994}{TreeSearch} of each clade in a reference +#' tree: the number of extra steps required before the clade is no longer +#' present in an optimal tree. Formally, for a clade _C_, +#' \deqn{Bremer(C) = L(\neg C) - L^*}{Bremer(C) = L(not C) - L*} +#' where \eqn{L(\neg C)}{L(not C)} is the length of the shortest tree that does +#' *not* contain _C_ and \eqn{L^*}{L*} is the length of the most-parsimonious +#' tree. Larger values indicate better-supported clades. +#' +#' Two engines are available: +#' +#' \describe{ +#' \item{`method = "constraint"` (default, rigorous)}{Runs one *converse- +#' constraint* search per clade, forcing that clade to be absent and taking +#' the length of the shortest resulting tree. This directly targets +#' \eqn{L(\neg C)}{L(not C)} and needs only bounded memory (one tree per +#' clade). It is the reliable choice for reported support values.} +#' \item{`method = "pool"` (fast, approximate)}{Collects a pool of suboptimal +#' trees with [`SuboptimalTrees()`] and, for each clade, takes the shortest +#' retained tree that lacks it. This is quick but **over-estimates** +#' support (the minimum over a sampled subset can only exceed the true +#' minimum) and is **right-censored** at the sampling depth `maxBremer`: a +#' clade broken by no retained tree is reported as `Inf` and flagged in the +#' `censored` attribute. Use it as a preview, not for publication values.} +#' } +#' +#' The reference `tree` supplies the clades to be evaluated. Pass a single +#' `phylo` (e.g. one most-parsimonious tree) or a `multiPhylo` search result, in +#' which case the strict consensus is used and support is calculated only for +#' its resolved bipartitions. Scoring options (`concavity`, `inapplicable`, +#' ...) **must match those used to find the trees**, or the extra-step counts +#' will be meaningless; they default to equal-weights Fitch parsimony. +#' +#' @param tree A most-parsimonious tree (`phylo`) whose clades are to be +#' evaluated, or a `multiPhylo` whose strict consensus provides them. A +#' `multiPhylo` returned by [`MaximizeParsimony()`] additionally supplies the +#' optimal score \eqn{L^*}{L*} via its `score` attribute. +#' @param dataset A phylogenetic data matrix of class `phyDat`. +#' @param method Character: `"constraint"` (default) for rigorous converse- +#' constraint searches, or `"pool"` for the fast suboptimal-pool approximation. +#' @param maxBremer Numeric: the largest decay value to resolve. Under +#' `method = "pool"` this is the pool's suboptimality depth (defaults to `10`); +#' clades not broken within it are censored. Under `method = "constraint"` it +#' is ignored (exact values are computed). +#' @param optimalScore Numeric: the optimal tree length \eqn{L^*}{L*}. If +#' `NULL` (default) it is taken from `attr(tree, "score")` when available, or +#' computed by search. Supplying it (or a scored `multiPhylo`) avoids a +#' redundant search. +#' @param format Character specifying return format, as in [`JackLabels()`]: +#' `"numeric"` (default) returns named numeric values for further analysis; +#' `"character"` returns a vector shaped for `phylo$node.label`. +#' @inheritParams MaximizeParsimony +#' @param \dots Further arguments passed to [`MaximizeParsimony()`] / +#' [`SuboptimalTrees()`], e.g. `maxReplicates`, `maxSeconds`, `strategy`, +#' `nThreads`, `verbosity`. +#' +#' @return A numeric vector (or, if `format = "character"`, a +#' `phylo$node.label`-shaped character vector) giving the Bremer support of each +#' resolved clade in `tree`, named by node number (the row names of +#' [`TreeTools::as.Splits()`]). Annotate a plot with +#' [`TreeTools::LabelSplits()`], or assign to `tree$node.label`. Under +#' `method = "pool"` the numeric result carries a logical `censored` attribute +#' marking clades whose support exceeds `maxBremer`. +#' +#' @examples +#' data("inapplicable.phyData", package = "TreeSearch") +#' dataset <- inapplicable.phyData[["Vinther2008"]] +#' \donttest{ +#' trees <- MaximizeParsimony(dataset, maxReplicates = 8, verbosity = 0) +#' decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) +#' decay +#' +#' # Annotate a tree +#' TreeTools::LabelSplits(trees[[1]], decay) +#' } +#' @references \insertAllCited{} +#' @template MRS +#' @seealso +#' Other clade support measures: [`JackLabels()`], [`SiteConcordance`]; +#' [`SuboptimalTrees()`] collects the pool used by `method = "pool"`. +#' @family split support functions +#' @importFrom TreeTools as.Splits NTip SplitFrequency TipLabels +#' @export +Bremer <- function(tree, dataset, + concavity = Inf, extended_iw = TRUE, xpiwe_r = 0.5, + xpiwe_max_f = 5, hierarchy = NULL, inapplicable = "bgs", + hsj_alpha = 1.0, + method = c("constraint", "pool"), + maxBremer = Inf, optimalScore = NULL, + format = "numeric", ...) { + method <- match.arg(method) + + scoringArgs <- list(concavity = concavity, extended_iw = extended_iw, + xpiwe_r = xpiwe_r, xpiwe_max_f = xpiwe_max_f, + hierarchy = hierarchy, inapplicable = inapplicable, + hsj_alpha = hsj_alpha) + + # --- Resolve the reference tree, its clades, and L* --- + ref <- .BremerReference(tree, dataset, optimalScore) + + if (length(ref$splitNames) == 0L) { + warning("Reference tree has no resolved internal clades; ", + "nothing to calculate.") + return(.BremerFormat(setNames(numeric(0), character(0)), + logical(0), ref$reference, format)) + } + + engine <- if (method == "pool") .BremerPool else .BremerConstraint + res <- engine(ref, dataset, scoringArgs, maxBremer, list(...)) + + .BremerFormat(res$bremer, res$censored, ref$reference, format) +} + +# Resolve reference tree + clades (Splits) + optimal score L*. +#' @importFrom TreeTools as.Splits NTip +.BremerReference <- function(tree, dataset, optimalScore) { + Lstar <- optimalScore + if (inherits(tree, "multiPhylo")) { + if (is.null(Lstar)) { + s <- attr(tree, "score") + if (!is.null(s) && is.finite(s)) Lstar <- s + } + reference <- if (length(tree) == 1L) tree[[1L]] else ape::consensus(tree, p = 1) + } else if (inherits(tree, "phylo")) { + reference <- tree + } else { + stop("`tree` must be a `phylo` or `multiPhylo` object.") + } + + splits <- as.Splits(reference, tipLabels = names(dataset)) + splitNames <- rownames(as.matrix(splits)) + if (is.null(splitNames)) splitNames <- character(0) + + list(reference = reference, splits = splits, splitNames = splitNames, + Lstar = Lstar) +} + +# Format a node-indexed numeric result, mirroring JackLabels(). +#' @importFrom TreeTools NTip +.BremerFormat <- function(values, censored, reference, format) { + numericFmt <- c("numeric", "number", "double") + characterFmt <- c("character", "text") + returnMode <- c(rep("numeric", length(numericFmt)), + rep("character", length(characterFmt)))[ + pmatch(tolower(format), c(numericFmt, characterFmt))] + if (is.na(returnMode)) returnMode <- "numeric" + + switch(returnMode, + "character" = { + ret <- character(reference[["Nnode"]]) + if (length(values)) { + idx <- as.integer(names(values)) - NTip(reference) + ret[idx] <- as.character(values) + } + ret + }, + { + if (length(censored)) attr(values, "censored") <- censored + values + }) +} + +# Approximate Bremer from a pool of suboptimal trees. +# Returns list(bremer = named numeric, censored = logical). +#' @importFrom TreeTools SplitFrequency +.BremerPool <- function(ref, dataset, scoringArgs, maxBremer, dots) { + poolK <- if (is.finite(maxBremer)) maxBremer else 10 + if (!is.finite(maxBremer)) { + message("method = \"pool\": using maxBremer = ", poolK, + " for the pool depth (set `maxBremer` to control it).") + } + + pool <- do.call(SuboptimalTrees, + c(list(dataset = dataset, maxSuboptimal = poolK), + scoringArgs, dots)) + poolScores <- attr(pool, "scores") + if (is.null(poolScores)) { + poolScores <- do.call(TreeLength, + c(list(tree = pool, dataset = dataset), scoringArgs)) + } + poolBest <- min(poolScores) + + Lstar <- ref$Lstar + tol <- 1e-8 + if (is.null(Lstar)) { + Lstar <- poolBest + } else if (poolBest < Lstar - tol) { + warning("Pool search found a tree (length ", signif(poolBest, 7), + ") shorter than the supplied optimalScore (", signif(Lstar, 7), + "); adopting the shorter length as L*.") + Lstar <- poolBest + } + + # Per-clade displaying matrix: refFreq gives the node names; then test each + # pool tree individually (SplitFrequency of a length-1 forest is 0/1). + refFreq <- SplitFrequency(ref$reference, pool) + nSplits <- length(refFreq) + displayMat <- vapply(seq_along(pool), function(j) { + SplitFrequency(ref$reference, pool[j]) + }, double(nSplits)) + dim(displayMat) <- c(nSplits, length(pool)) + + bremer <- numeric(nSplits) + censored <- logical(nSplits) + for (i in seq_len(nSplits)) { + nonDisplaying <- displayMat[i, ] < 0.5 + if (any(nonDisplaying)) { + bremer[i] <- min(poolScores[nonDisplaying]) - Lstar + } else { + bremer[i] <- Inf + censored[i] <- TRUE + } + } + names(bremer) <- names(refFreq) + names(censored) <- names(refFreq) + list(bremer = bremer, censored = censored) +} + +# Rigorous Bremer by converse-constraint search: for each clade, the shortest +# tree forced to LACK it. Uses the negative-constraint engine wired into +# MaximizeParsimony() via the internal `.negativeConstraint` argument. +# +# The converse search keeps the full search machinery -- Wagner starts, TBR, +# ratchet, sectorial search and NNI perturbation -- which all re-optimize +# through the negative-constraint-guarded TBR and so stay in the space of trees +# lacking the clade; the tree pool additionally rejects any tree that displays +# it. Only the phases that accept score-worsening moves through an unguarded +# path (drift, in-sector drift, simulated annealing) are disabled -- they would +# otherwise wander onto the clade and, being worse-accepting, report it as +# unsupported. Runs serially: the pool guard is on the serial search path. +.BremerConstraint <- function(ref, dataset, scoringArgs, maxBremer, dots) { + disabled <- list( + driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, + annealCycles = 0L + ) + overrides <- dots[intersect(names(dots), names(disabled))] + converseFixed <- modifyList(disabled, overrides) + # nThreads is forced to 1 (the negative-constraint pool guard is on the + # serial search path only); everything else the user passes flows through. + passthrough <- dots[setdiff(names(dots), c(names(disabled), "nThreads"))] + + # L* (the unconstrained optimum) is found by a full-strength search -- the + # worse-accepting phases are safe here because there is no forbidden clade. + Lstar <- ref$Lstar + if (is.null(Lstar)) { + res0 <- do.call(MaximizeParsimony, + c(list(dataset = dataset, collapse = TRUE, nThreads = 1L), + scoringArgs, passthrough)) + Lstar <- attr(res0, "score") + } + + runConverse <- function(negSplit) { + args <- c(list(dataset = dataset, collapse = TRUE, nThreads = 1L, + .negativeConstraint = negSplit), + scoringArgs, converseFixed, passthrough) + attr(do.call(MaximizeParsimony, args), "score") + } + + splitNames <- ref$splitNames + scores <- vapply(seq_along(splitNames), function(i) { + s <- runConverse(ref$splits[[i]]) + # A negative score is the engine's "no tree found" sentinel (empty pool): + # the converse search failed to locate any tree lacking the clade. + if (!is.finite(s) || s < 0) NA_real_ else s + }, double(1)) + + if (anyNA(scores)) { + warning(sum(is.na(scores)), " converse-constraint search(es) found no tree ", + "lacking the clade; reported as NA. Increase `maxReplicates`.") + } + + trueBest <- suppressWarnings(min(c(Lstar, scores), na.rm = TRUE)) + if (is.finite(trueBest) && trueBest < Lstar - 1e-8) { + warning("A converse-constraint search found a tree (length ", + signif(trueBest, 7), ") shorter than the optimal score (", + signif(Lstar, 7), "); the original search was suboptimal. ", + "Adopting the shorter length as L*.") + Lstar <- trueBest + } + + bremer <- setNames(scores - Lstar, splitNames) + list(bremer = bremer, censored = logical(length(bremer))) +} diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index e1d454bf8..bf6628f37 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -57,6 +57,38 @@ } # Internal helper: prepare constraint data for C++ engine. +# Build the forbidden-clade (negative / converse constraint) matrix for the +# C++ engine. A tree that displays any of these bipartitions will be rejected. +# @param negConstraint A `Splits` object, or anything `as.Splits()` accepts +# (e.g. a `phylo`), giving the clade(s) to forbid. +# @param dataset A phyDat whose names define the tip ordering. +# @return An integer n_neg x n_tips membership matrix (columns in dataset tip +# order), or NULL if there is nothing to forbid. +# @keywords internal +#' @importFrom TreeTools as.Splits +.PrepareNegativeConstraint <- function(negConstraint, dataset) { + if (is.null(negConstraint)) return(NULL) + tipLabels <- names(dataset) + if (!inherits(negConstraint, "Splits")) { + negConstraint <- as.Splits(negConstraint, tipLabels = tipLabels) + } + if (length(negConstraint) == 0L) return(NULL) + + # as.logical() yields an n_splits x n_tips logical membership matrix. + membership <- as.logical(negConstraint) + if (is.null(dim(membership))) { + membership <- matrix(membership, nrow = 1L) + } + splitLabels <- attr(negConstraint, "tip.label") + colOrder <- match(tipLabels, splitLabels) + if (anyNA(colOrder)) { + stop("Forbidden-clade taxa do not match the dataset taxa.") + } + membership <- membership[, colOrder, drop = FALSE] + storage.mode(membership) <- "integer" + membership +} + # Returns a named list of constraint arguments (empty list if no constraint). # @param constraint A phyDat, phylo, or NULL. # @param dataset A phyDat whose names define the tip ordering. @@ -585,6 +617,12 @@ #' attributes: #' \describe{ #' \item{`score`}{Best parsimony score.} +#' \item{`scores`}{Present only when `collapse = FALSE`: a numeric vector of +#' the parsimony score of each returned tree, aligned with the returned +#' `multiPhylo` (the same values are attached as a `score` attribute on +#' each individual tree). With `poolSuboptimal > 0` this exposes the +#' retained suboptimal pool for landscape analysis (see +#' [`SuboptimalTrees()`], [`Suboptimality()`]).} #' \item{`replicates`}{Number of replicates completed.} #' \item{`hits_to_best`}{Number of independent discoveries of the best #' score.} @@ -648,6 +686,7 @@ MaximizeParsimony <- function( inapplicable = "bgs", hsj_alpha = 1.0, constraint, + .negativeConstraint = NULL, strategy = "auto", maxReplicates = 96L, targetHits = NULL, @@ -973,6 +1012,18 @@ MaximizeParsimony <- function( cli_alert_info("Constraint: {nrow(consArgs$consSplitMatrix)} split{?s}") } + # Negative (converse) constraints: forbidden clades. Internal argument used + # by Bremer() -- the returned trees must NOT display any supplied split. + if (!is.null(.negativeConstraint)) { + negMatrix <- .PrepareNegativeConstraint(.negativeConstraint, dataset) + if (!is.null(negMatrix)) { + consArgs[["consNegSplitMatrix"]] <- negMatrix + if (verbosity > 0L) { + cli_alert_info("Forbidding {nrow(negMatrix)} clade{?s}") + } + } + } + # --- Profile parsimony: extract info_amounts --- profileArgs <- list() if (useProfile) { @@ -1064,6 +1115,11 @@ MaximizeParsimony <- function( resultTrees <- list() } nTopologies <- result$n_topologies + # Per-tree scores aligned with result$trees; surfaced only in the collapse = + # FALSE path (the collapse = TRUE path keeps only best-score trees, so their + # scores are all result$best_score and the collapsed/deduped list no longer + # aligns with result$scores). NULL here -> no "scores" attribute is attached. + perTreeScores <- NULL if (isTRUE(collapse) && length(resultTrees) > 0L) { # Contract zero-length (unsupported) branches into polytomies, à la TNT's # "collapse zero-length branches" -- done entirely in C++ (ts_collapse_pool) @@ -1112,9 +1168,26 @@ MaximizeParsimony <- function( # C++ edge order may differ from template; renumber to valid preorder Renumber(tr) }) + # Surface the per-tree parsimony scores the engine returns aligned with + # result$trees. Attaching a "score" attribute to each tree lets + # Suboptimality(), SuboptimalTrees() and Bremer(method = "pool") read the + # suboptimal pool's landscape directly without re-scoring. + perTreeScores <- result$scores + if (length(perTreeScores) == length(outTrees)) { + outTrees <- Map(function(tr, sc) { + attr(tr, "score") <- sc + tr + }, outTrees, perTreeScores) + } else { + perTreeScores <- NULL + } } if (length(outTrees) == 0L) { outTrees <- list(treeTpl) + if (!isTRUE(collapse)) { + attr(outTrees[[1]], "score") <- result$best_score + perTreeScores <- result$best_score + } } # --- Output --- @@ -1137,6 +1210,7 @@ MaximizeParsimony <- function( structure( outTrees, score = result$best_score, + scores = perTreeScores, replicates = result$replicates, hits_to_best = result$hits_to_best, n_topologies = nTopologies, diff --git a/R/SuboptimalTrees.R b/R/SuboptimalTrees.R new file mode 100644 index 000000000..a14e8d492 --- /dev/null +++ b/R/SuboptimalTrees.R @@ -0,0 +1,73 @@ +#' Collect suboptimal trees for landscape analysis +#' +#' `SuboptimalTrees()` performs a parsimony search with [`MaximizeParsimony()`] +#' and returns *every* tree retained within a specified number of steps of the +#' optimum, each annotated with its parsimony score. This exposes the shape of +#' the parsimony landscape near the optimum -- for example to visualise the +#' distribution of near-optimal scores, to measure tree-to-tree distances among +#' competing resolutions, or as input to a fast approximate +#' [`Bremer()`][Bremer] support calculation (`method = "pool"`). +#' +#' The retained trees are the search engine's internal tree pool. Its size is +#' bounded by `maxPool`, so memory use is capped even when many trees fall +#' within `maxSuboptimal` steps; once the pool is full the engine evicts trees +#' to preserve topological diversity. The pool therefore reflects the islands +#' the search actually visited, and is not an exhaustive enumeration of every +#' tree within `maxSuboptimal` steps: raise `maxReplicates` / `maxSeconds` +#' (passed via `...`) for a denser sample. +#' +#' @param dataset A phylogenetic data matrix of class `phyDat`, as accepted by +#' [`MaximizeParsimony()`]. +#' @param tree Optional starting tree (of class `phylo`) or `multiPhylo`; if +#' `NULL` (default) the search begins from random addition sequence trees. +#' @param maxSuboptimal Numeric: retain trees scoring up to this many steps +#' worse than the best tree found (in the search's optimality units -- integer +#' steps under equal weights, fractional under implied weights or profile +#' parsimony). Sets `poolSuboptimal` in [`SearchControl()`]. +#' @param maxPool Integer: maximum number of trees to retain in the pool. Sets +#' `poolMaxSize` in [`SearchControl()`]; raised from the search default (100) +#' so that a suboptimal sample is not prematurely truncated. +#' @param \dots Further arguments passed to [`MaximizeParsimony()`], including +#' scoring options (`concavity`, `inapplicable`, ...) and search effort +#' (`maxReplicates`, `maxSeconds`, `strategy`, `nThreads`, `verbosity`). +#' Named [`SearchControl()`] fields may also be passed here to override the +#' constructed control. +#' +#' @return A `multiPhylo` object listing the retained trees, best tree(s) +#' first. Each tree carries a `score` attribute giving its parsimony length, +#' and the object carries a `scores` attribute: a numeric vector of those +#' lengths aligned with the returned trees. [`Suboptimality()`] reports each +#' tree's excess over the optimum. +#' +#' @examples +#' data("inapplicable.phyData", package = "TreeSearch") +#' dataset <- inapplicable.phyData[["Vinther2008"]] +#' \donttest{ +#' trees <- SuboptimalTrees(dataset, maxSuboptimal = 3, maxReplicates = 8, +#' verbosity = 0) +#' attr(trees, "scores") # parsimony length of each retained tree +#' Suboptimality(trees) # excess over the optimum +#' } +#' @seealso +#' [`MaximizeParsimony()`] performs the underlying search; +#' [`Bremer()`][Bremer] uses this pool for approximate decay indices; +#' [`Suboptimality()`] summarises the excess length of each tree. +#' @template MRS +#' @family split support functions +#' @export +SuboptimalTrees <- function(dataset, tree = NULL, + maxSuboptimal = 5, maxPool = 1000L, ...) { + if (!is.numeric(maxSuboptimal) || length(maxSuboptimal) != 1L || + maxSuboptimal < 0) { + stop("`maxSuboptimal` must be a single non-negative number.") + } + maxPool <- as.integer(maxPool) + if (is.na(maxPool) || maxPool < 1L) { + stop("`maxPool` must be a single positive integer.") + } + + control <- SearchControl(poolSuboptimal = maxSuboptimal, poolMaxSize = maxPool) + + MaximizeParsimony(dataset, tree = tree, control = control, + collapse = FALSE, ...) +} diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 0c16caba7..c5cf21ae0 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -462,6 +462,28 @@ @article{Mir2013 number = {1} } +@article{Bremer1988, + title = {The limits of amino acid sequence data in angiosperm phylogenetic reconstruction}, + author = {Bremer, Kare}, + year = {1988}, + volume = {42}, + pages = {795--803}, + doi = {10.1111/j.1558-5646.1988.tb02497.x}, + journal = {Evolution}, + number = {4} +} + +@article{Bremer1994, + title = {Branch support and tree stability}, + author = {Bremer, Kare}, + year = {1994}, + volume = {10}, + pages = {295--304}, + doi = {10.1111/j.1096-0031.1994.tb00179.x}, + journal = {Cladistics}, + number = {3} +} + @article{Muller2005, title = {The efficiency of different search strategies in estimating parsimony jackknife, bootstrap, and {{Bremer}} support}, todo_author = {M{\"u}ller, Kai F.}, diff --git a/man/Bremer.Rd b/man/Bremer.Rd new file mode 100644 index 000000000..b098a734f --- /dev/null +++ b/man/Bremer.Rd @@ -0,0 +1,198 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Bremer.R +\name{Bremer} +\alias{Bremer} +\title{Bremer (decay) support} +\usage{ +Bremer( + tree, + dataset, + concavity = Inf, + extended_iw = TRUE, + xpiwe_r = 0.5, + xpiwe_max_f = 5, + hierarchy = NULL, + inapplicable = "bgs", + hsj_alpha = 1, + method = c("constraint", "pool"), + maxBremer = Inf, + optimalScore = NULL, + format = "numeric", + ... +) +} +\arguments{ +\item{tree}{A most-parsimonious tree (\code{phylo}) whose clades are to be +evaluated, or a \code{multiPhylo} whose strict consensus provides them. A +\code{multiPhylo} returned by \code{\link[=MaximizeParsimony]{MaximizeParsimony()}} additionally supplies the +optimal score \eqn{L^*}{L*} via its \code{score} attribute.} + +\item{dataset}{A phylogenetic data matrix of class \code{phyDat}.} + +\item{concavity}{Determines the degree to which extra steps beyond the first +are penalized. Specify a numeric value to use implied weighting +\insertCite{Goloboff1993}{TreeSearch}; \code{concavity} specifies \emph{k} in +\emph{k} / \emph{e} + \emph{k}. A value of 10 is recommended; +TNT sets a default of 3, but this is too low in some circumstances +\insertCite{Goloboff2018,Smith2019}{TreeSearch}. +Better still explore the sensitivity of results under a range of +concavity values, e.g. \code{k = 2 ^ (1:7)}. +Specify \code{Inf} to weight each additional step equally. +Specify \code{"profile"} to employ profile parsimony +\insertCite{Faith2001}{TreeSearch}.} + +\item{extended_iw}{Logical: if \code{TRUE} (default) and \code{concavity} is finite, +apply the missing-entries correction of +\insertCite{Goloboff2014;textual}{TreeSearch}. +Characters with missing data receive a reduced effective concavity +\emph{k_c} = \emph{k} / \emph{f_c}, making their weights drop off faster. +This compensates for the artificially low homoplasy of poorly sampled +characters. Set \code{FALSE} for legacy Goloboff (1993) behaviour. +Ignored when \code{concavity = Inf} (equal weights) or \code{"profile"}.} + +\item{xpiwe_r}{Numeric in (0, 1]: proportion of observed homoplasy +expected in unobserved (missing) entries. Default 0.5 (following TNT). +Only used when \code{extended_iw = TRUE}.} + +\item{xpiwe_max_f}{Numeric >= 1: maximum extrapolation factor. +Characters with very few observed entries are clamped so that the +extrapolation factor does not exceed this value. Default 5 (following +TNT). Only used when \code{extended_iw = TRUE}.} + +\item{hierarchy}{A \code{\link{CharacterHierarchy}} object specifying which +characters are controlling primaries and which are their dependent +secondaries. Required when \code{inapplicable} is \code{"hsj"} or \code{"xform"}; +ignored when \code{inapplicable = "bgs"} (the default). +See \code{\link[=CharacterHierarchy]{CharacterHierarchy()}} for how to construct one, and +\code{\link[=HierarchyFromNames]{HierarchyFromNames()}} for automated construction from +TNT-style character names.} + +\item{inapplicable}{Character: method for handling inapplicable characters. +Case-insensitive. +See \code{vignette("inapplicable", package = "TreeSearch")} for details. +\describe{ +\item{\code{"bgs"} (default)}{Three-pass algorithm of +\insertCite{Brazeau2019;textual}{TreeSearch}, inferring applicability +regions from the \code{"-"} token. No hierarchy required.} +\item{\code{"missing"}}{Pure Fitch parsimony +\insertCite{Fitch1971}{TreeSearch}: the inapplicable (\code{"-"}) state is +treated as missing data, so any token that includes a gap is recoded +as fully ambiguous (\code{"?"}) and contributes no steps -- including +polymorphisms such as \verb{\{0,-\}}, which become \verb{?}. Reproduces standard +Fitch analyses (e.g. PAUP*, or TNT with gaps read as missing) that do +not use the Brazeau-Gardner-Smith inapplicable algorithm. No +hierarchy required.} +\item{\code{"hsj"}}{Dissimilarity-metric scoring of +\insertCite{Hopkins2021;textual}{TreeSearch}. Requires a +\code{hierarchy}; controlled by \code{hsj_alpha}.} +\item{\code{"xform"}}{Step-matrix recoding approximating maximum homology +via x-transformations +\insertCite{Goloboff2021;textual}{TreeSearch}. Requires a +\code{hierarchy}.} +}} + +\item{hsj_alpha}{Numeric in [0, 1]: scaling parameter for secondary- +character contributions under the HSJ method. 0 = secondaries ignored; +1 (default) = secondaries contribute up to 1 per branch per hierarchy +block. Only used when \code{inapplicable = "hsj"}.} + +\item{method}{Character: \code{"constraint"} (default) for rigorous converse- +constraint searches, or \code{"pool"} for the fast suboptimal-pool approximation.} + +\item{maxBremer}{Numeric: the largest decay value to resolve. Under +\code{method = "pool"} this is the pool's suboptimality depth (defaults to \code{10}); +clades not broken within it are censored. Under \code{method = "constraint"} it +is ignored (exact values are computed).} + +\item{optimalScore}{Numeric: the optimal tree length \eqn{L^*}{L*}. If +\code{NULL} (default) it is taken from \code{attr(tree, "score")} when available, or +computed by search. Supplying it (or a scored \code{multiPhylo}) avoids a +redundant search.} + +\item{format}{Character specifying return format, as in \code{\link[=JackLabels]{JackLabels()}}: +\code{"numeric"} (default) returns named numeric values for further analysis; +\code{"character"} returns a vector shaped for \code{phylo$node.label}.} + +\item{\dots}{Further arguments passed to \code{\link[=MaximizeParsimony]{MaximizeParsimony()}} / +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, e.g. \code{maxReplicates}, \code{maxSeconds}, \code{strategy}, +\code{nThreads}, \code{verbosity}.} +} +\value{ +A numeric vector (or, if \code{format = "character"}, a +\code{phylo$node.label}-shaped character vector) giving the Bremer support of each +resolved clade in \code{tree}, named by node number (the row names of +\code{\link[TreeTools:as.Splits]{TreeTools::as.Splits()}}). Annotate a plot with +\code{\link[TreeTools:LabelSplits]{TreeTools::LabelSplits()}}, or assign to \code{tree$node.label}. Under +\code{method = "pool"} the numeric result carries a logical \code{censored} attribute +marking clades whose support exceeds \code{maxBremer}. +} +\description{ +\code{Bremer()} calculates the Bremer support (decay index) +\insertCite{Bremer1988,Bremer1994}{TreeSearch} of each clade in a reference +tree: the number of extra steps required before the clade is no longer +present in an optimal tree. Formally, for a clade \emph{C}, +\deqn{Bremer(C) = L(\neg C) - L^*}{Bremer(C) = L(not C) - L*} +where \eqn{L(\neg C)}{L(not C)} is the length of the shortest tree that does +\emph{not} contain \emph{C} and \eqn{L^*}{L*} is the length of the most-parsimonious +tree. Larger values indicate better-supported clades. +} +\details{ +Two engines are available: + +\describe{ +\item{\code{method = "constraint"} (default, rigorous)}{Runs one \emph{converse- +constraint} search per clade, forcing that clade to be absent and taking +the length of the shortest resulting tree. This directly targets +\eqn{L(\neg C)}{L(not C)} and needs only bounded memory (one tree per +clade). It is the reliable choice for reported support values.} +\item{\code{method = "pool"} (fast, approximate)}{Collects a pool of suboptimal +trees with \code{\link[=SuboptimalTrees]{SuboptimalTrees()}} and, for each clade, takes the shortest +retained tree that lacks it. This is quick but \strong{over-estimates} +support (the minimum over a sampled subset can only exceed the true +minimum) and is \strong{right-censored} at the sampling depth \code{maxBremer}: a +clade broken by no retained tree is reported as \code{Inf} and flagged in the +\code{censored} attribute. Use it as a preview, not for publication values.} +} + +The reference \code{tree} supplies the clades to be evaluated. Pass a single +\code{phylo} (e.g. one most-parsimonious tree) or a \code{multiPhylo} search result, in +which case the strict consensus is used and support is calculated only for +its resolved bipartitions. Scoring options (\code{concavity}, \code{inapplicable}, +...) \strong{must match those used to find the trees}, or the extra-step counts +will be meaningless; they default to equal-weights Fitch parsimony. +} +\examples{ +data("inapplicable.phyData", package = "TreeSearch") +dataset <- inapplicable.phyData[["Vinther2008"]] +\donttest{ +trees <- MaximizeParsimony(dataset, maxReplicates = 8, verbosity = 0) +decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) +decay + +# Annotate a tree +TreeTools::LabelSplits(trees[[1]], decay) +} +} +\references{ +\insertAllCited{} +} +\seealso{ +Other clade support measures: \code{\link[=JackLabels]{JackLabels()}}, \code{\link{SiteConcordance}}; +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} collects the pool used by \code{method = "pool"}. + +Other split support functions: +\code{\link[=ConcordanceTable]{ConcordanceTable()}}, +\code{\link[=JackLabels]{JackLabels()}}, +\code{\link[=Jackknife]{Jackknife()}}, +\code{\link[=Morphy]{Morphy()}}, +\code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, +\code{\link[=PaintCharacters]{PaintCharacters()}}, +\code{\link[=PresCont]{PresCont()}}, +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} +} +\author{ +\href{https://smithlabdurham.github.io/}{Martin R. Smith} +(\href{mailto:martin.smith@durham.ac.uk}{martin.smith@durham.ac.uk}) +} +\concept{split support functions} diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 3cae0d982..de62a5923 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -123,12 +123,14 @@ image(t(`mode<-`(PhyDatToMatrix(dataset), "numeric")), axes = FALSE, } Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \concept{split support functions} diff --git a/man/JackLabels.Rd b/man/JackLabels.Rd index e11dfa48e..c48686781 100644 --- a/man/JackLabels.Rd +++ b/man/JackLabels.Rd @@ -85,13 +85,15 @@ Generate trees by jackknife resampling using \code{\link[=Resample]{Resample()}} parsimony searches, or \code{\link[=Jackknife]{Jackknife()}} for custom search criteria. Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \author{ \href{https://smithlabdurham.github.io/}{Martin R. Smith} diff --git a/man/Jackknife.Rd b/man/Jackknife.Rd index 888a5b037..2ff9e07d3 100644 --- a/man/Jackknife.Rd +++ b/man/Jackknife.Rd @@ -77,13 +77,15 @@ engine. } Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} Other custom search functions: \code{\link[=EdgeListSearch]{EdgeListSearch()}}, diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 5940a0b41..a41efe208 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -17,6 +17,7 @@ MaximizeParsimony( inapplicable = "bgs", hsj_alpha = 1, constraint, + .negativeConstraint = NULL, strategy = "auto", maxReplicates = 96L, targetHits = NULL, @@ -158,6 +159,16 @@ replicates (default: 96). The default is a multiple of 48 (= LCM(12, 16)) so that replicates divide evenly across common 12- or 16-core machines when running in parallel. +When \code{strategy} resolves to \code{"large"} (automatically selected for +datasets of \eqn{\ge}{>=} 120 tips and \eqn{\ge}{>=} 100 characters) and +\code{maxReplicates} is left at its default, the +cap is raised to 500: a 120--180-tip sweep showed the fraction of runs +reaching the best-known score climbing from 0.68 at 96 replicates to 0.79 +at 250 and still rising at 500, with no plateau. Raising the cap only +appends later replicates, so it never delays an earlier improvement, and +easy datasets still stop early once \code{targetHits} is met; only genuinely +hard datasets run the extra replicates. An explicit \code{maxReplicates} +is always respected. For large or complex datasets a higher value improves the chance of finding all MPTs. A rough minimum is \code{max(10, ceiling(NTip * NChar / 5000))}, where \code{NChar = sum(weight)}. @@ -234,6 +245,12 @@ A \code{multiPhylo} object containing the best tree(s) found, with attributes: \describe{ \item{\code{score}}{Best parsimony score.} +\item{\code{scores}}{Present only when \code{collapse = FALSE}: a numeric vector of +the parsimony score of each returned tree, aligned with the returned +\code{multiPhylo} (the same values are attached as a \code{score} attribute on +each individual tree). With \code{poolSuboptimal > 0} this exposes the +retained suboptimal pool for landscape analysis (see +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=Suboptimality]{Suboptimality()}}).} \item{\code{replicates}}{Number of replicates completed.} \item{\code{hits_to_best}}{Number of independent discoveries of the best score.} diff --git a/man/Morphy.Rd b/man/Morphy.Rd index 0bbd579de..acb7f8702 100644 --- a/man/Morphy.Rd +++ b/man/Morphy.Rd @@ -396,13 +396,15 @@ Morphy(characters, constraint = constraint, verbosity = 0) Tree search \emph{via} graphical user interface: \code{\link[=EasyTrees]{EasyTrees()}} Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \author{ \href{https://smithlabdurham.github.io/}{Martin R. Smith} diff --git a/man/MostContradictedFreq.Rd b/man/MostContradictedFreq.Rd index 31404191c..c5b135259 100644 --- a/man/MostContradictedFreq.Rd +++ b/man/MostContradictedFreq.Rd @@ -36,12 +36,14 @@ contradicted" score. \code{PresCont()} calculates the "groups present / contradicted" score. Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \concept{split support functions} diff --git a/man/PaintCharacters.Rd b/man/PaintCharacters.Rd index bef233339..7e5cdba84 100644 --- a/man/PaintCharacters.Rd +++ b/man/PaintCharacters.Rd @@ -65,12 +65,14 @@ plot(tree, edge.color = paint$edgeCol, edge.width = 2) \code{\link[TreeTools:PaintTree]{TreeTools::PaintTree()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}} Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \concept{split support functions} diff --git a/man/PresCont.Rd b/man/PresCont.Rd index 5e1b38546..a7f285297 100644 --- a/man/PresCont.Rd +++ b/man/PresCont.Rd @@ -103,13 +103,15 @@ of trees that contain the split, and the frequency of the most contradicted split, respectively. Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \author{ \href{https://smithlabdurham.github.io/}{Martin R. Smith} diff --git a/man/SearchControl.Rd b/man/SearchControl.Rd index f492ca563..3c49cda79 100644 --- a/man/SearchControl.Rd +++ b/man/SearchControl.Rd @@ -39,6 +39,13 @@ SearchControl( sectorMaxHits = 1L, sectorCollapseTarget = 0L, rssPicks = 0L, + sectorGoDrift = 0L, + sectorGoComb = 0L, + sectorDriftCycles = 5L, + sectorDriftAfd = 3L, + sectorDriftRfd = 0.1, + sectorCombStarts = 3L, + sectorFuseRounds = 3L, postRatchetSectorial = FALSE, fuseInterval = 3L, fuseAcceptEqual = FALSE, @@ -210,6 +217,34 @@ auto-sizes to \code{2 * nTip / meanSectorSize} (~5). Higher values chain more sequential sector replacements before the next global cleanup, matching the ~20-25 of \insertCite{Goloboff1999;textual}{TreeSearch} / \acronym{TNT}.} +\item{sectorGoDrift, sectorGoComb}{Integer; size thresholds (in real sector +tips) above which a sector is solved by tree-drifting (\code{sectorGoDrift}) or by +combined analysis (\code{sectorGoComb}) instead of plain +\acronym{RAS}+\acronym{TBR}, following \acronym{TNT}'s \code{sectsch} +\code{godrift}/\code{gocomb}. Small sectors are cheap to optimise exhaustively; large +sectors have more reach but need drift to escape their own local optima. +\code{sectorGoComb} solves the sector with \code{sectorCombStarts} \acronym{RAS}+drift +starts and then \emph{fuses} them (recombines shared clades across the starts over +\code{sectorFuseRounds} rounds), keeping the fused tree only if it beats the best +start. \code{sectorGoComb} takes precedence when both trigger. \code{0} (default) +disables each, leaving all sectors on \acronym{RAS}+\acronym{TBR}. To +exercise these, also raise \code{sectorMaxSize} so large sectors are selected.} + +\item{sectorDriftCycles}{Integer; drift cycles used when a sector is solved by +drift or combined analysis (\acronym{TNT} \code{drift}).} + +\item{sectorDriftAfd}{Integer; absolute-fit-difference limit (steps) for +in-sector drift.} + +\item{sectorDriftRfd}{Numeric; relative-fit-difference limit for in-sector +drift.} + +\item{sectorCombStarts}{Integer; \acronym{RAS}+drift starts per combined +sector (\acronym{TNT} \code{combstarts}).} + +\item{sectorFuseRounds}{Integer; max fuse rounds for the combined-analysis +recombination sub-step (\acronym{TNT} \code{fuse}).} + \item{postRatchetSectorial}{Logical; when \code{TRUE}, run XSS+RSS+CSS again after ratchet perturbation using the same round counts. Approximates TNT's interleaved sectorial pattern. Default: \code{FALSE}.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 0139ea63f..e550d237d 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -245,13 +245,15 @@ ClusteringConcordance(TreeTools::NJTree(myMatrix), myMatrix) } Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=Jackknife]{Jackknife()}}, \code{\link[=Morphy]{Morphy()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, -\code{\link[=PresCont]{PresCont()}} +\code{\link[=PresCont]{PresCont()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \author{ \href{https://smithlabdurham.github.io/}{Martin R. Smith} diff --git a/man/SuboptimalTrees.Rd b/man/SuboptimalTrees.Rd new file mode 100644 index 000000000..1225acea8 --- /dev/null +++ b/man/SuboptimalTrees.Rd @@ -0,0 +1,86 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/SuboptimalTrees.R +\name{SuboptimalTrees} +\alias{SuboptimalTrees} +\title{Collect suboptimal trees for landscape analysis} +\usage{ +SuboptimalTrees(dataset, tree = NULL, maxSuboptimal = 5, maxPool = 1000L, ...) +} +\arguments{ +\item{dataset}{A phylogenetic data matrix of class \code{phyDat}, as accepted by +\code{\link[=MaximizeParsimony]{MaximizeParsimony()}}.} + +\item{tree}{Optional starting tree (of class \code{phylo}) or \code{multiPhylo}; if +\code{NULL} (default) the search begins from random addition sequence trees.} + +\item{maxSuboptimal}{Numeric: retain trees scoring up to this many steps +worse than the best tree found (in the search's optimality units -- integer +steps under equal weights, fractional under implied weights or profile +parsimony). Sets \code{poolSuboptimal} in \code{\link[=SearchControl]{SearchControl()}}.} + +\item{maxPool}{Integer: maximum number of trees to retain in the pool. Sets +\code{poolMaxSize} in \code{\link[=SearchControl]{SearchControl()}}; raised from the search default (100) +so that a suboptimal sample is not prematurely truncated.} + +\item{\dots}{Further arguments passed to \code{\link[=MaximizeParsimony]{MaximizeParsimony()}}, including +scoring options (\code{concavity}, \code{inapplicable}, ...) and search effort +(\code{maxReplicates}, \code{maxSeconds}, \code{strategy}, \code{nThreads}, \code{verbosity}). +Named \code{\link[=SearchControl]{SearchControl()}} fields may also be passed here to override the +constructed control.} +} +\value{ +A \code{multiPhylo} object listing the retained trees, best tree(s) +first. Each tree carries a \code{score} attribute giving its parsimony length, +and the object carries a \code{scores} attribute: a numeric vector of those +lengths aligned with the returned trees. \code{\link[=Suboptimality]{Suboptimality()}} reports each +tree's excess over the optimum. +} +\description{ +\code{SuboptimalTrees()} performs a parsimony search with \code{\link[=MaximizeParsimony]{MaximizeParsimony()}} +and returns \emph{every} tree retained within a specified number of steps of the +optimum, each annotated with its parsimony score. This exposes the shape of +the parsimony landscape near the optimum -- for example to visualise the +distribution of near-optimal scores, to measure tree-to-tree distances among +competing resolutions, or as input to a fast approximate +\code{\link[=Bremer]{Bremer()}} support calculation (\code{method = "pool"}). +} +\details{ +The retained trees are the search engine's internal tree pool. Its size is +bounded by \code{maxPool}, so memory use is capped even when many trees fall +within \code{maxSuboptimal} steps; once the pool is full the engine evicts trees +to preserve topological diversity. The pool therefore reflects the islands +the search actually visited, and is not an exhaustive enumeration of every +tree within \code{maxSuboptimal} steps: raise \code{maxReplicates} / \code{maxSeconds} +(passed via \code{...}) for a denser sample. +} +\examples{ +data("inapplicable.phyData", package = "TreeSearch") +dataset <- inapplicable.phyData[["Vinther2008"]] +\donttest{ +trees <- SuboptimalTrees(dataset, maxSuboptimal = 3, maxReplicates = 8, + verbosity = 0) +attr(trees, "scores") # parsimony length of each retained tree +Suboptimality(trees) # excess over the optimum +} +} +\seealso{ +\code{\link[=MaximizeParsimony]{MaximizeParsimony()}} performs the underlying search; +\code{\link[=Bremer]{Bremer()}} uses this pool for approximate decay indices; +\code{\link[=Suboptimality]{Suboptimality()}} summarises the excess length of each tree. + +Other split support functions: +\code{\link[=Bremer]{Bremer()}}, +\code{\link[=ConcordanceTable]{ConcordanceTable()}}, +\code{\link[=JackLabels]{JackLabels()}}, +\code{\link[=Jackknife]{Jackknife()}}, +\code{\link[=Morphy]{Morphy()}}, +\code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, +\code{\link[=PaintCharacters]{PaintCharacters()}}, +\code{\link[=PresCont]{PresCont()}}, +\code{\link{SiteConcordance}} +} +\author{ +\href{https://smithlabdurham.github.io/}{Martin R. Smith} +(\href{mailto:martin.smith@durham.ac.uk}{martin.smith@durham.ac.uk}) +} +\concept{split support functions} diff --git a/src/ts_constraint.cpp b/src/ts_constraint.cpp index 34a400d26..4fc66486d 100644 --- a/src/ts_constraint.cpp +++ b/src/ts_constraint.cpp @@ -149,6 +149,91 @@ std::vector compute_node_tips(const TreeState& tree, int n_words) return node_tips; } +// ========================================================================= +// Negative (converse) constraints +// ========================================================================= + +void add_negative_constraint( + ConstraintData& cd, const int* neg_matrix, int n_neg, int n_tips) +{ + if (n_neg <= 0) return; + + int n_words = (n_tips + 63) / 64; + if (!cd.active) { + // Initialize an otherwise-empty ConstraintData so the search's + // `constrained` path runs (it will be a no-op over zero positive splits). + cd.active = true; + cd.n_words = n_words; + int n_node = 2 * n_tips - 1; + cd.dfs_entry.assign(n_node, 0); + cd.dfs_exit.assign(n_node, 0); + cd.clip_tip_mask.resize(n_words, 0ULL); + } + + cd.neg_active = true; + cd.n_neg_splits = n_neg; + cd.neg_split_tips.assign( + static_cast(n_neg) * cd.n_words, 0ULL); + + // Pack rows into bitmasks (column-major from R: [s, t] at s + n_neg * t), + // canonicalizing so tip 0 is always outside (bit 0 = 0), exactly as + // build_constraint() does for positive splits. + for (int s = 0; s < n_neg; ++s) { + uint64_t* mask = &cd.neg_split_tips[static_cast(s) * cd.n_words]; + for (int t = 0; t < n_tips; ++t) { + if (neg_matrix[s + n_neg * t]) { + mask[t / 64] |= (1ULL << (t % 64)); + } + } + if (mask[0] & 1ULL) { + for (int w = 0; w < cd.n_words; ++w) { + mask[w] = ~mask[w]; + } + int remainder = n_tips % 64; + if (remainder > 0) { + mask[cd.n_words - 1] &= (1ULL << remainder) - 1; + } + } + } +} + +bool displays_forbidden_clade(const TreeState& tree, const ConstraintData& cd) +{ + if (!cd.neg_active || cd.n_neg_splits == 0) return false; + + auto node_tips = compute_node_tips(tree, cd.n_words); + + // The internal tree is not guaranteed to be rooted on tip 0, so a bipartition + // can be displayed as EITHER side. Canonicalize each node's subtree mask to + // the tip-0-outside form (flipping if tip 0 is inside) before comparing to the + // canonical forbidden split; this matches a clade regardless of which side is + // the rooted subtree. + int n_tips = tree.n_tip; + int rem = n_tips % 64; + uint64_t last_mask = rem ? ((1ULL << rem) - 1) : ~0ULL; + + std::vector canon(cd.n_words); + for (int node : tree.postorder) { + const uint64_t* nd = &node_tips[static_cast(node) * cd.n_words]; + bool flip = (nd[0] & 1ULL) != 0; + for (int w = 0; w < cd.n_words; ++w) { + canon[w] = flip ? ~nd[w] : nd[w]; + } + if (flip) canon[cd.n_words - 1] &= last_mask; + + for (int s = 0; s < cd.n_neg_splits; ++s) { + const uint64_t* split = + &cd.neg_split_tips[static_cast(s) * cd.n_words]; + bool match = true; + for (int w = 0; w < cd.n_words; ++w) { + if (canon[w] != split[w]) { match = false; break; } + } + if (match) return true; + } + } + return false; +} + // ========================================================================= // Map constraint nodes: find which internal node holds each split // ========================================================================= diff --git a/src/ts_constraint.h b/src/ts_constraint.h index cb86016dc..1bc494c35 100644 --- a/src/ts_constraint.h +++ b/src/ts_constraint.h @@ -58,6 +58,16 @@ struct ConstraintData { int expected_score = 0; bool has_posthoc = false; + // --- Negative (converse) constraints: forbidden clades --- + // A tree is REJECTED if it displays ANY of these bipartitions. Used for + // Bremer support (the shortest tree that does NOT contain a given clade). + // Stored as canonicalized tip bitmasks (tip 0 outside), exactly like + // split_tips. Independent of the positive splits above: a ConstraintData + // may carry positive splits, negative splits, or (for Bremer) only negative. + bool neg_active = false; + int n_neg_splits = 0; + std::vector neg_split_tips; // n_neg_splits * n_words + // Per-clip workspace (reused across clips, sized at init) std::vector clip_zones; // [n_splits] std::vector clip_tip_mask; // [n_words] @@ -132,6 +142,23 @@ bool violates_constraint_posthoc(const TreeState& tree, // For tips: bit[t] = 1. For internal nodes: OR of children. std::vector compute_node_tips(const TreeState& tree, int n_words); +// --- Negative (converse) constraints --- + +// Add negative (converse) constraints from an R-side split matrix +// (n_neg rows x n_tips cols, column-major, 0/1 membership). A tree that +// displays any of these bipartitions is rejected. Initializes `cd` (active, +// n_words, DFS arrays, clip workspace) if it is not already built, so a +// negative-only ConstraintData is safe to pass through the positive-constraint +// machinery (which then runs as a no-op over zero positive splits). +void add_negative_constraint( + ConstraintData& cd, const int* neg_matrix, int n_neg, int n_tips); + +// Returns true if the tree currently displays any forbidden clade. The tree +// is rooted on tip 0, so a canonical (tip-0-outside) bipartition is displayed +// iff some internal node's subtree tip mask equals it -- the same test +// map_constraint_nodes() uses to locate a positive split. +bool displays_forbidden_clade(const TreeState& tree, const ConstraintData& cd); + // Repair constraint violations by minimal SPR moves. // After return, all constraint splits are displayed and // update_constraint() has been called. Caller must rescore. diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index bc7272073..ae033f1aa 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -143,6 +143,19 @@ ReplicateResult run_single_replicate( } } + // Converse-constraint (Bremer) search: guarantee a start tree that does NOT + // display any forbidden clade. A Wagner tree often builds the very clade we + // must avoid, and the guarded TBR below cannot climb out of it downhill, so + // the pool would stay empty. Random topologies almost never display an + // arbitrary clade, so fall back to them until one is clean. + if (cd && cd->neg_active && displays_forbidden_clade(result.tree, *cd)) { + int tries = 0; + do { + random_topology_tree(result.tree, ds); + } while (displays_forbidden_clade(result.tree, *cd) && ++tries < 100); + best_wag = score_tree(result.tree, ds); + } + result.timings.wagner_ms = ph_lap(); if (verbosity >= 2) { if (starting_tree) { @@ -661,6 +674,14 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, // Reset the candidate-evaluation counter for this search (serial path). ds.n_candidates_evaluated = 0; + // Converse-constraint (Bremer) search: make the pool reject any tree that + // displays a forbidden clade, so the returned best tree/score can never + // contain it regardless of which phase produced the candidate. The TBR + // hill-climber (above) additionally steers the search into that space. + if (cd && cd->neg_active) { + pool.set_forbidden(cd); + } + bool use_timeout = params.max_seconds > 0.0; auto start_time = std::chrono::steady_clock::now(); diff --git a/src/ts_pool.cpp b/src/ts_pool.cpp index 397c0f6b6..0419618f6 100644 --- a/src/ts_pool.cpp +++ b/src/ts_pool.cpp @@ -1,4 +1,5 @@ #include "ts_pool.h" +#include "ts_constraint.h" // displays_forbidden_clade (negative-constraint filter) #include #include #include @@ -22,6 +23,11 @@ bool TreePool::add(const TreeState& tree, double score) { return false; } + // Reject any tree displaying a forbidden clade (converse-constraint search). + if (forbidden_ && displays_forbidden_clade(tree, *forbidden_)) { + return false; + } + SplitSet ss = compute_splits(tree); uint64_t hash = hash_splits(ss); @@ -105,6 +111,11 @@ bool TreePool::add_collapsed(const TreeState& tree, double score, return false; } + // Reject any tree displaying a forbidden clade (converse-constraint search). + if (forbidden_ && displays_forbidden_clade(tree, *forbidden_)) { + return false; + } + // Compute collapsed splits (skip zero-length edges) for dedup SplitSet css = compute_collapsed_splits(tree, collapsed); uint64_t chash = hash_splits(css); diff --git a/src/ts_pool.h b/src/ts_pool.h index 84327301d..9a9666a18 100644 --- a/src/ts_pool.h +++ b/src/ts_pool.h @@ -18,6 +18,8 @@ namespace ts { +struct ConstraintData; // fwd decl (ts_constraint.h); used only by pointer + // Per-split frequency table for conflict-guided sector selection. // Maps per-split hash → count across best-score pool trees. struct SplitFrequencyTable { @@ -46,6 +48,13 @@ class TreePool { best_score_(1e18), hits_to_best_(0), consensus_hash_(0), consensus_unchanged_(0) {} + // Register negative (converse) constraints: any tree displaying a forbidden + // clade is rejected by add()/add_collapsed(). Guarantees the pool -- and + // hence the returned best tree/score -- never contains a forbidden clade, + // regardless of which search phase produced the candidate. Pass nullptr + // (default) to disable. Used for Bremer converse-constraint searches. + void set_forbidden(const ConstraintData* f) { forbidden_ = f; } + // Add a tree if it's not a duplicate and meets score threshold. // Returns true if the tree was actually added. bool add(const TreeState& tree, double score); @@ -118,6 +127,9 @@ class TreePool { double best_score_; int hits_to_best_; + // Negative-constraint filter (not owned); nullptr disables. + const ConstraintData* forbidden_ = nullptr; + // Consensus stability tracking uint64_t consensus_hash_; int consensus_unchanged_; diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 4d742b375..b182afaa4 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -276,7 +276,8 @@ static ts::ConstraintData build_constraint_from_r( Nullable consTipData, Nullable consWeight, Nullable consLevels, - int consExpectedScore); + int consExpectedScore, + Nullable consNegSplitMatrix = R_NilValue); // [[Rcpp::export]] double ts_fitch_score( @@ -1432,7 +1433,8 @@ static ts::ConstraintData build_constraint_from_r( Nullable consTipData, Nullable consWeight, Nullable consLevels, - int consExpectedScore) + int consExpectedScore, + Nullable consNegSplitMatrix) { ts::ConstraintData cd; if (consSplitMatrix.isNotNull()) { @@ -1469,6 +1471,17 @@ static ts::ConstraintData build_constraint_from_r( } } } + + // Negative (converse) constraints: forbidden clades (Bremer support). May be + // supplied alone (cd otherwise empty) or alongside positive splits. + if (consNegSplitMatrix.isNotNull()) { + IntegerMatrix nsm(consNegSplitMatrix.get()); + int n_neg = nsm.nrow(); + if (n_neg > 0 && nsm.ncol() == n_tips) { + ts::add_negative_constraint(cd, INTEGER(nsm), n_neg, n_tips); + } + } + return cd; } @@ -1691,6 +1704,9 @@ static ts::ConstraintData unpack_constraint(int n_tips, if (cc.containsElementNamed("consExpectedScore")) consExpectedScore = as(cc["consExpectedScore"]); + SEXP nsm = cc.containsElementNamed("consNegSplitMatrix") + ? SEXP(cc["consNegSplitMatrix"]) : R_NilValue; + return build_constraint_from_r( n_tips, Nullable(csm), @@ -1698,7 +1714,8 @@ static ts::ConstraintData unpack_constraint(int n_tips, Nullable(ctd), Nullable(cw), Nullable(cl), - consExpectedScore); + consExpectedScore, + Nullable(nsm)); } return ts::ConstraintData{}; } diff --git a/src/ts_tbr.cpp b/src/ts_tbr.cpp index a00978a7b..06df78b3f 100644 --- a/src/ts_tbr.cpp +++ b/src/ts_tbr.cpp @@ -1291,6 +1291,10 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, // Initialize constraint mapping if active bool constrained = cd && cd->active; + // Negative (converse) constraint: reject any accepted move whose result + // displays a forbidden clade, directing the hill-climb into the space of + // trees that lack it (used for Bremer support). + bool neg_constrained = cd && cd->neg_active; if (constrained) { update_constraint(tree, *cd); } @@ -2337,6 +2341,21 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, } } + // Negative-constraint guard: reject any move whose resulting tree + // displays a forbidden clade. Mirrors the positive post-hoc block + // above (revert the applied move and skip), keeping every accepted + // tree -- and hence best_score and the pool -- free of the clade. + if (neg_constrained && displays_forbidden_clade(tree, *cd)) { + restore_topology(tree, snap); + state_snap.restore(tree); + score_fresh = true; + if (constrained) { + map_constraint_nodes(tree, *cd); + compute_dfs_timestamps(tree, *cd); + } + continue; + } + // Compute topology hash for tabu checking uint64_t tree_hash = 0; if (tabu.active()) { diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R new file mode 100644 index 000000000..d2398386d --- /dev/null +++ b/tests/testthat/test-Bremer.R @@ -0,0 +1,213 @@ +library("TreeTools", quietly = TRUE) + +data("inapplicable.phyData", package = "TreeSearch") +ds <- inapplicable.phyData[["Vinther2008"]] + +# Exact Bremer by exhaustive enumeration of every binary tree on the taxa +# (feasible for <= 7 tips): the minimum length among trees LACKING each +# reference clade, minus the global minimum length. The gold-standard oracle +# for the converse-constraint engine. +allBinaryTrees <- function(tips) { + trees <- list(ape::read.tree( + text = paste0("(", tips[1], ",", tips[2], ",", tips[3], ");"))) + for (k in 4:length(tips)) { + trees <- unlist(lapply(trees, function(tr) { + AddTipEverywhere(tr, tips[k], includeRoot = FALSE) + }), recursive = FALSE) + } + trees +} + +oracleBremer <- function(reference, dataset, concavity = Inf) { + tips <- names(dataset) + trees <- allBinaryTrees(tips) + trees <- structure(lapply(trees, RootTree, tips[1]), class = "multiPhylo") + lengths <- TreeLength(trees, dataset, concavity = concavity) + Lstar <- min(lengths) + refFreq <- SplitFrequency(reference, trees) + nSplits <- length(refFreq) + disp <- vapply(seq_along(trees), function(j) { + SplitFrequency(reference, trees[j]) + }, double(nSplits)) + dim(disp) <- c(nSplits, length(trees)) + out <- vapply(seq_len(nSplits), function(i) { + lacking <- disp[i, ] < 0.5 + if (any(lacking)) min(lengths[lacking]) - Lstar else Inf + }, double(1)) + setNames(out, names(refFreq)) +} + +# --- method = "pool" (approximate engine) --- +# A multiPhylo of several MPTs is summarised by its strict consensus; support +# is calculated for the consensus's resolved clades only. + +test_that("Bremer(method = 'pool') returns node-keyed non-negative support", { + set.seed(3418) + trees <- MaximizeParsimony(ds, maxReplicates = 6L, targetHits = 2L, + verbosity = 0L) + refCons <- if (length(trees) == 1L) trees[[1]] else ape::consensus(trees, p = 1) + + decay <- Bremer(trees, ds, method = "pool", maxBremer = 6, + maxReplicates = 6L, targetHits = 2L, verbosity = 0L) + + expect_type(decay, "double") + # One value per resolved clade of the consensus reference. + expect_length(decay, NSplits(refCons)) + # Names are internal node numbers. + expect_true(all(as.integer(names(decay)) > NTip(ds))) + + # Finite (uncensored) support values are non-negative. + expect_true(all(decay[is.finite(decay)] >= 0)) + + # Censored attribute marks exactly the infinite entries. + cens <- attr(decay, "censored") + expect_type(cens, "logical") + expect_length(cens, length(decay)) + expect_equal(unname(cens), unname(!is.finite(decay))) +}) + +test_that("Bremer(format = 'character') yields a node.label-shaped vector", { + set.seed(3418) + trees <- MaximizeParsimony(ds, maxReplicates = 6L, targetHits = 2L, + verbosity = 0L) + refCons <- if (length(trees) == 1L) trees[[1]] else ape::consensus(trees, p = 1) + + lab <- Bremer(trees, ds, method = "pool", maxBremer = 6, format = "character", + maxReplicates = 6L, targetHits = 2L, verbosity = 0L) + expect_type(lab, "character") + expect_length(lab, refCons[["Nnode"]]) + # Assignable straight onto the reference tree. + refCons[["node.label"]] <- lab + expect_length(refCons[["node.label"]], refCons[["Nnode"]]) +}) + +test_that("Bremer accepts a single phylo reference with an explicit L*", { + set.seed(3418) + trees <- MaximizeParsimony(ds, maxReplicates = 6L, targetHits = 2L, + verbosity = 0L) + Lstar <- attr(trees, "score") + decay <- Bremer(trees[[1]], ds, method = "pool", maxBremer = 6, + optimalScore = Lstar, + maxReplicates = 6L, targetHits = 2L, verbosity = 0L) + expect_type(decay, "double") + expect_length(decay, NSplits(trees[[1]])) + expect_true(all(decay[is.finite(decay)] >= 0)) + expect_setequal(names(decay), rownames(as.matrix(as.Splits(trees[[1]])))) +}) + +test_that("Bremer errors on a bad reference type", { + expect_error(Bremer("not a tree", ds, method = "pool"), + "must be a `phylo` or `multiPhylo`") +}) + +# --- method = "constraint" (rigorous engine): validated against brute force --- + +test_that("the oracle enumerates every binary tree (guards the gate itself)", { + # (2n - 5)!! unrooted binary trees; a gappy enumerator would make the oracle + # itself an over-estimate that could rubber-stamp an over-estimating search. + doubleFactorial <- function(m) prod(seq(m, 1, by = -2)) + expect_equal(length(allBinaryTrees(LETTERS[1:6])), doubleFactorial(2 * 6 - 5)) + expect_equal(length(allBinaryTrees(c(LETTERS[1:6], "out"))), + doubleFactorial(2 * 7 - 5)) +}) + +test_that("constraint Bremer matches enumeration when L* is computed internally", { + # Exercises the full-strength L* search path (no optimalScore supplied). + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + con <- Bremer(ref, dat, method = "constraint", # no optimalScore -> internal L* + maxReplicates = 25L, verbosity = 0L) + orc <- oracleBremer(ref, dat) + expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) +}) + +test_that("constraint Bremer matches exact enumeration (graded pectinate)", { + dat <- StringToPhyDat( + "1100000 1110000 1111000 1111100 1100000 1110000 1111000 1111100 1001000", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + con <- Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score"), + maxReplicates = 25L, verbosity = 0L) + orc <- oracleBremer(ref, dat) + expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) +}) + +test_that("constraint Bremer matches enumeration in the adversarial case", { + # {A,B,C} is supported by three characters (decay 3); the optimum DISPLAYS it + # and the nearest tree lacking it is three steps longer -- the case a stalling + # post-hoc filter would over-estimate. + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + con <- Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score"), + maxReplicates = 25L, verbosity = 0L) + orc <- oracleBremer(ref, dat) + expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) + expect_gte(max(con), 3) +}) + +test_that("constraint Bremer matches enumeration under implied weights", { + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, concavity = 3, maxReplicates = 8L, + verbosity = 0L) + ref <- mpts[[1]] + con <- Bremer(ref, dat, method = "constraint", concavity = 3, + optimalScore = attr(mpts, "score"), + maxReplicates = 25L, verbosity = 0L) + orc <- oracleBremer(ref, dat, concavity = 3) + expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) + # Fractional decay under implied weights. + expect_false(all(con == round(con))) +}) + +test_that("constraint is the default and dominates the pool estimate", { + set.seed(42) + mpts <- MaximizeParsimony(ds, maxReplicates = 8L, verbosity = 0L) + con <- Bremer(mpts, ds, maxReplicates = 12L, verbosity = 0L) # default method + expect_type(con, "double") + expect_true(all(con >= 0)) + expect_true(all(is.finite(con))) + + pool <- Bremer(mpts, ds, method = "pool", maxBremer = 8, + maxReplicates = 8L, verbosity = 0L) + # Both are upper bounds on the true decay; the directed constraint search is + # never looser than the incidental pool where the pool is uncensored. + shared <- intersect(names(con), names(pool)) + uncensored <- shared[is.finite(pool[shared])] + expect_true(all(con[uncensored] <= pool[uncensored] + 1e-6)) +}) + +test_that("Bremer warns and stays non-negative when optimalScore is too high", { + dat <- StringToPhyDat("1100000 1110000 1111000 1111100 1001000", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + # Claim an optimum 5 steps worse than reality; converse searches find shorter + # trees, so the guard must warn and clamp rather than report negative decay. + expect_warning( + con <- Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score") + 5, + maxReplicates = 20L, verbosity = 0L), + "shorter than the optimal score") + expect_true(all(con >= 0)) +}) diff --git a/tests/testthat/test-SuboptimalTrees.R b/tests/testthat/test-SuboptimalTrees.R new file mode 100644 index 000000000..e4940d788 --- /dev/null +++ b/tests/testthat/test-SuboptimalTrees.R @@ -0,0 +1,65 @@ +library("TreeTools", quietly = TRUE) + +data("inapplicable.phyData", package = "TreeSearch") +ds <- inapplicable.phyData[["Vinther2008"]] + +# --- Scores surfaced only when collapse = FALSE --- + +test_that("MaximizeParsimony(collapse = TRUE) attaches no per-tree scores", { + set.seed(3418) + result <- MaximizeParsimony(ds, maxReplicates = 2L, targetHits = 1L, + verbosity = 0L) # collapse = TRUE by default + expect_null(attr(result, "scores")) + expect_true(is.finite(attr(result, "score"))) +}) + +test_that("MaximizeParsimony(collapse = FALSE) surfaces aligned pool scores", { + set.seed(3418) + result <- MaximizeParsimony(ds, maxReplicates = 3L, targetHits = 1L, + verbosity = 0L, collapse = FALSE, + control = SearchControl(poolSuboptimal = 5, + poolMaxSize = 500L)) + scores <- attr(result, "scores") + expect_type(scores, "double") + expect_length(scores, length(result)) + + # Each tree carries a matching `score` attribute (so Suboptimality() works). + perTree <- vapply(result, attr, double(1), "score") + expect_equal(perTree, scores) + + # Scores match an independent re-scoring of the trees. + expect_equal(scores, TreeLength(result, ds), tolerance = 1e-8) + + # The best of the pool equals the reported best score, and every retained + # tree is within the requested suboptimality window. + expect_equal(min(scores), attr(result, "score")) + expect_true(all(scores <= attr(result, "score") + 5 + 1e-8)) +}) + +# --- SuboptimalTrees() accessor --- + +test_that("SuboptimalTrees() returns a scored, bounded pool", { + set.seed(3418) + trees <- SuboptimalTrees(ds, maxSuboptimal = 3, maxPool = 200L, + maxReplicates = 3L, targetHits = 1L, + verbosity = 0L) + expect_s3_class(trees, "multiPhylo") + scores <- attr(trees, "scores") + expect_false(is.null(scores)) + expect_length(scores, length(trees)) + expect_true(length(trees) <= 200L) + expect_equal(min(scores), attr(trees, "score")) + expect_true(all(scores <= attr(trees, "score") + 3 + 1e-8)) + + # Suboptimality() reads the per-tree `score` attribute we attached. + subopt <- Suboptimality(trees) + expect_equal(subopt, scores - min(scores)) + expect_true(all(subopt >= 0)) +}) + +test_that("SuboptimalTrees() validates its arguments", { + expect_error(SuboptimalTrees(ds, maxSuboptimal = -1), + "must be a single non-negative number") + expect_error(SuboptimalTrees(ds, maxPool = 0L), + "must be a single positive integer") +}) From 02a7f7e20dd3dbc10380320f23f08104adc84cd9 Mon Sep 17 00:00:00 2001 From: R script Date: Thu, 9 Jul 2026 17:12:00 +0100 Subject: [PATCH 02/12] fix(bremer): address red-team findings (scoring guard, reach, soundness) Adversarial review (two background red-team agents) of the Bremer implementation surfaced ten findings; this lands the fixes. Parallel fan-out follows separately. R (Bremer.R / SuboptimalTrees.R / MaximizeParsimony.R): - BR-1 (P1): guard against a scoring-units mismatch. A multiPhylo whose stored optimal score was measured under different scoring args (e.g. trees found under implied weights, then Bremer() called with the default equal weights) silently subtracted an IW optimum from EW lengths. Re-score the reference under the active scoring args (resolving collapsed polytomies) and error on a mismatch, or on a supplied optimalScore longer than the reference tree. - BR-2 (P2): drift / in-sector drift / annealing are unsound in a converse search; a user value no longer re-enables them (was modifyList, silently emptying the pool -> all NA). Ignore with a warning. - BR-6 (P3): verify in R that each converse-search tree really lacks the clade, surfacing an engine regression as a visible NA rather than a deflated decay. - BR-3 (P3): SuboptimalTrees() strips a colliding control=/collapse= from ... (was a hard "matched by multiple actual arguments" crash on the pool path). - BR-5 (P3): carry the censored flag into format="character" output. - BRM-2 (P2, latent): MaximizeParsimony() forces nThreads=1 when a negative constraint is set (the parallel pool has no soundness backstop). C++ (ts_tbr.cpp / ts_driven.cpp), all gated on neg_active so the mainline path is byte-identical (confirmed: 180 constraint/driven/parallel/sector/fuse/ratchet/ wagner test blocks unchanged): - BRM-1 (P2): the outer-reroot sweep (try_root_edge_moves / exact_verify_sweep) was unguarded, so a root-edge TBR move could carry the search into forbidden space and inflate the reported decay (or spuriously NA). Snapshot the clade-free tree and restore+stop if the improved tree displays the clade. - BRM-3 (P3): gate the periodic-fuse improvement bookkeeping on the actual pool insertion, so a backstop-rejected fused tree no longer resets the stopping rules. - BRM-4 (perf): skip the positive post-hoc block when there are no positive splits (a negative-only constraint), removing a redundant compute_node_tips sweep. Tests: NA/bgs enumeration oracle for the DEFAULT scoring mode (previously unvalidated); BR-1 and BR-2 regressions; de-vacuumed the cross-engine test (it had been comparing empty vectors). All oracle checks (EW / IW / bgs / adversarial) pass; broad regression 180 blocks / 0 failures. Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 122 +++++++++++++++++++++++++++++++++-- R/MaximizeParsimony.R | 9 +++ R/SuboptimalTrees.R | 25 ++++++- src/ts_driven.cpp | 10 ++- src/ts_tbr.cpp | 26 +++++++- tests/testthat/test-Bremer.R | 111 ++++++++++++++++++++++++++----- 6 files changed, 275 insertions(+), 28 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index 22ca465ce..b2c59899c 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -107,6 +107,11 @@ Bremer <- function(tree, dataset, logical(0), ref$reference, format)) } + # Scoring-units guard: L* must be measured with the same scoring arguments as + # the converse searches / pool re-scores, or the decay mixes two optimality + # criteria (e.g. an equal-weights length minus an implied-weights optimum). + .BremerCheckScoring(tree, dataset, scoringArgs, optimalScore) + engine <- if (method == "pool") .BremerPool else .BremerConstraint res <- engine(ref, dataset, scoringArgs, maxBremer, list(...)) @@ -137,6 +142,55 @@ Bremer <- function(tree, dataset, Lstar = Lstar) } +# Error if the reference's optimal score L* was measured under different scoring +# arguments than those now in effect -- a silent-units-mismatch guard. A decay +# of "extra steps" is only meaningful when L* and L(not C) use the same ruler. +.BremerCheckScoring <- function(tree, dataset, scoringArgs, optimalScore) { + refLen <- function() { + # Collapse-derived MPTs may be multifurcating, but TreeLength needs binary + # trees. Resolving the zero-length polytomies preserves the parsimony length + # (every resolution of an equally-optimal collapsed polytomy is itself + # optimal), and taking the minimum guards against an unlucky resolution. + resolved <- if (inherits(tree, "multiPhylo")) { + structure(lapply(tree, ape::multi2di, random = FALSE), class = "multiPhylo") + } else { + ape::multi2di(tree, random = FALSE) + } + min(suppressWarnings( + do.call(TreeLength, c(list(tree = resolved, dataset = dataset), scoringArgs)))) + } + tol <- 1e-6 + if (!is.null(optimalScore)) { + # A user-asserted optimum cannot exceed the reference's own length under + # these arguments; if it does, the score or the arguments are wrong. + Lmin <- refLen() + if (is.finite(Lmin) && optimalScore > Lmin + tol) { + stop("The supplied `optimalScore` (", signif(optimalScore, 7), + ") exceeds the reference tree's length under the supplied scoring ", + "arguments (", signif(Lmin, 7), "). Check that `optimalScore` and ", + "the scoring arguments (`concavity`, `inapplicable`, ...) match ", + "those used to find the trees.", call. = FALSE) + } + } else if (inherits(tree, "multiPhylo")) { + # L* will come from the trees' stored optimal score. Because the supplied + # trees are optimal, their minimum length UNDER scoringArgs equals L* under + # scoringArgs; a disagreement means the scoring arguments differ from those + # used to find them -- fatal, as the decay would be meaningless. + attrScore <- attr(tree, "score") + if (!is.null(attrScore) && is.finite(attrScore)) { + Lmin <- refLen() + if (is.finite(Lmin) && abs(attrScore - Lmin) > tol) { + stop("The reference trees' length under the supplied scoring arguments ", + "(", signif(Lmin, 7), ") does not equal their stored optimal score ", + "(", signif(attrScore, 7), "). Pass the same `concavity`, ", + "`inapplicable` and other scoring arguments used to find the trees.", + call. = FALSE) + } + } + } + invisible(NULL) +} + # Format a node-indexed numeric result, mirroring JackLabels(). #' @importFrom TreeTools NTip .BremerFormat <- function(values, censored, reference, format) { @@ -153,6 +207,14 @@ Bremer <- function(tree, dataset, if (length(values)) { idx <- as.integer(names(values)) - NTip(reference) ret[idx] <- as.character(values) + # Inf renders as the literal "Inf"; preserve the censoring flag as an + # attribute so pool-method callers can still distinguish "> maxBremer" + # (censored) from a genuine value (the numeric format carries it too). + if (length(censored) && any(censored)) { + censVec <- logical(reference[["Nnode"]]) + censVec[idx] <- as.logical(censored) + attr(ret, "censored") <- censVec + } } ret }, @@ -235,8 +297,18 @@ Bremer <- function(tree, dataset, driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, annealCycles = 0L ) - overrides <- dots[intersect(names(dots), names(disabled))] - converseFixed <- modifyList(disabled, overrides) + # The worse-accepting phases (drift, in-sector drift, annealing) reach trees + # through a path the negative-constraint hill-climb does not guard, so they + # are ALWAYS disabled here. A user value would only reintroduce the + # unsoundness the converse search exists to avoid, so warn rather than honour + # it (previously `modifyList` let the user win, silently emptying the pool). + reenabled <- intersect(names(dots), names(disabled)) + if (length(reenabled)) { + warning("Ignoring ", paste(reenabled, collapse = ", "), + " in the converse-constraint search: drift, in-sector drift and ", + "annealing are always disabled here for soundness.") + } + converseFixed <- disabled # nThreads is forced to 1 (the negative-constraint pool guard is on the # serial search path only); everything else the user passes flows through. passthrough <- dots[setdiff(names(dots), c(names(disabled), "nThreads"))] @@ -251,20 +323,56 @@ Bremer <- function(tree, dataset, Lstar <- attr(res0, "score") } + # One converse-constraint search: the shortest tree forced to lack `negSplit`. + # Returns the score AND the tree, so the caller can verify the clade really is + # absent (defence in depth around the C++ move-guard + pool backstop). runConverse <- function(negSplit) { args <- c(list(dataset = dataset, collapse = TRUE, nThreads = 1L, .negativeConstraint = negSplit), scoringArgs, converseFixed, passthrough) - attr(do.call(MaximizeParsimony, args), "score") + res <- do.call(MaximizeParsimony, args) + # MaximizeParsimony() may return a single `phylo` or a `multiPhylo`; take the + # first tree either way (never `res[[1L]]` on a bare phylo, which would grab + # its `edge` matrix). + tree <- if (inherits(res, "phylo")) { + res + } else if (inherits(res, "multiPhylo") && length(res) >= 1L) { + res[[1L]] + } else { + NULL + } + list(score = attr(res, "score"), tree = tree) } splitNames <- ref$splitNames - scores <- vapply(seq_along(splitNames), function(i) { - s <- runConverse(ref$splits[[i]]) + + # Per-clade work unit (also the unit of parallelism for the cluster fan-out). + processConverse <- function(i) { + out <- runConverse(ref$splits[[i]]) + s <- out$score # A negative score is the engine's "no tree found" sentinel (empty pool): # the converse search failed to locate any tree lacking the clade. - if (!is.finite(s) || s < 0) NA_real_ else s - }, double(1)) + if (!is.finite(s) || s < 0) { + return(NA_real_) + } + # Verify the returned tree really lacks the clade. The C++ guard + pool + # backstop already guarantee this; checking here means a future engine + # regression surfaces as a visible NA rather than a silently deflated decay + # for a well-supported clade (the worst outcome for a published statistic). + if (!is.null(out$tree)) { + disp <- SplitFrequency(ref$reference, + structure(list(out$tree), class = "multiPhylo")) + # Index by the reference node number (robust to split ordering). + if (isTRUE(unname(disp[splitNames[i]]) >= 0.5)) { + warning("The converse search for clade ", splitNames[i], + " returned a tree that still displays it; reported as NA.") + return(NA_real_) + } + } + s + } + + scores <- vapply(seq_along(splitNames), processConverse, double(1)) if (anyNA(scores)) { warning(sum(is.na(scores)), " converse-constraint search(es) found no tree ", diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index bf6628f37..ec1b04770 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -1015,6 +1015,15 @@ MaximizeParsimony <- function( # Negative (converse) constraints: forbidden clades. Internal argument used # by Bremer() -- the returned trees must NOT display any supplied split. if (!is.null(.negativeConstraint)) { + # The negative-constraint soundness backstop (TreePool::set_forbidden) is + # wired only into the serial search path; the parallel pool has no such + # guard, so a forbidden clade could slip into the result under nThreads > 1. + # Force serial (Bremer() already does; this protects any other caller). + if (!identical(as.integer(nThreads), 1L)) { + warning("Negative (converse) constraints are supported only in serial ", + "search; forcing `nThreads = 1`.") + nThreads <- 1L + } negMatrix <- .PrepareNegativeConstraint(.negativeConstraint, dataset) if (!is.null(negMatrix)) { consArgs[["consNegSplitMatrix"]] <- negMatrix diff --git a/R/SuboptimalTrees.R b/R/SuboptimalTrees.R index a14e8d492..89d267617 100644 --- a/R/SuboptimalTrees.R +++ b/R/SuboptimalTrees.R @@ -66,8 +66,29 @@ SuboptimalTrees <- function(dataset, tree = NULL, stop("`maxPool` must be a single positive integer.") } + # SuboptimalTrees() manages `control` (to size the pool) and `collapse` (which + # must stay FALSE to retain the full suboptimal sample). A whole `control=` + # or `collapse=` argument arriving through `...` would otherwise collide with + # these explicit arguments ("matched by multiple actual arguments"). Strip + # them; individual SearchControl() fields (e.g. `ratchetCycles = `) still pass + # through `...` and are merged by MaximizeParsimony(). + dots <- list(...) + if ("control" %in% names(dots)) { + warning("`control` is managed by SuboptimalTrees(); pass individual ", + "SearchControl() fields (e.g. `ratchetCycles = `) via `...` instead.") + dots[["control"]] <- NULL + } + if ("collapse" %in% names(dots)) { + if (!identical(dots[["collapse"]], FALSE)) { + warning("`collapse` is forced to FALSE by SuboptimalTrees() so the full ", + "suboptimal pool is returned.") + } + dots[["collapse"]] <- NULL + } + control <- SearchControl(poolSuboptimal = maxSuboptimal, poolMaxSize = maxPool) - MaximizeParsimony(dataset, tree = tree, control = control, - collapse = FALSE, ...) + do.call(MaximizeParsimony, + c(list(dataset = dataset, tree = tree, control = control, + collapse = FALSE), dots)) } diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index ae033f1aa..91eaa1a4e 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -1094,13 +1094,19 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, } } } + bool fused_added = false; if (fused_ok) { std::vector fused_collapsed; compute_collapsed_flags(fused, ds, fused_collapsed); - pool.add_collapsed(fused, fused_score, fused_collapsed); + // The pool can reject the fused tree even when it scores better -- e.g. + // a converse-constraint (Bremer) search's backstop rejects any tree + // that displays the forbidden clade. Only treat it as an improvement + // (resetting the perturb-stop / target-hits counters) if it actually + // entered the pool; otherwise the stopping rules act on a phantom best. + fused_added = pool.add_collapsed(fused, fused_score, fused_collapsed); } - if (fused_ok && fused_score < best_before) { + if (fused_added && fused_score < best_before) { pool.set_hits_to_best(0); result.last_improved_rep = rep1; unsuccessful_reps = 0; // fuse found a better score; reset perturb-stop counter diff --git a/src/ts_tbr.cpp b/src/ts_tbr.cpp index 06df78b3f..2ac169db5 100644 --- a/src/ts_tbr.cpp +++ b/src/ts_tbr.cpp @@ -2322,7 +2322,11 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, // clip phase (the rerooting changes which constraint tips // end up on which side of the attachment edge). Reject // any move that introduces a constraint violation. - if (constrained) { + // `cd->n_splits > 0`: a negative-only constraint (Bremer converse + // search) has `active == true` but zero positive splits, so this whole + // block is a no-op -- skip it to avoid a redundant compute_node_tips + // sweep on top of the displays_forbidden_clade sweep just below. + if (constrained && cd->n_splits > 0) { map_constraint_nodes(tree, *cd); bool violation = false; for (int _s = 0; _s < cd->n_splits; ++_s) { @@ -2349,7 +2353,7 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, restore_topology(tree, snap); state_snap.restore(tree); score_fresh = true; - if (constrained) { + if (constrained && cd->n_splits > 0) { map_constraint_nodes(tree, *cd); compute_dfs_timestamps(tree, *cd); } @@ -2517,9 +2521,27 @@ TBRResult tbr_search(TreeState& tree, const DataSet& ds, // fast additive for EW / apply+rescore for IW). NA: the indirect scan is // only approximate, so an EXACT full-neighbourhood sweep is required to // certify a true unrooted-TBR optimum (see exact_verify_sweep). + // Negative (converse) constraint: unlike the inner clip loop, this + // root-edge / full-neighbourhood sweep is NOT constraint-guarded, so its + // best move can introduce a forbidden clade (a root-edge TBR move changes + // the unrooted topology). Snapshot the clade-free tree first; if the + // improved tree displays the clade, keep the snapshot and stop. We are + // then at a local optimum within the space of trees lacking the clade -- + // conservative (an improving root-edge move that stays clade-free is + // forgone, costing reach) but never unsound (the pool backstop already + // guarantees soundness; this keeps best_score consistent with a ¬C tree). + TopoSnapshot neg_pre_snap; + if (neg_constrained) save_topology(tree, neg_pre_snap); bool improved = has_na ? exact_verify_sweep(tree, ds, best_score) : try_root_edge_moves(tree, ds, best_score, ew_directional); + if (neg_constrained && improved && displays_forbidden_clade(tree, *cd)) { + restore_topology(tree, neg_pre_snap); + tree.build_postorder(); + best_score = full_rescore(tree, ds); // pre-sweep score; refreshes states + score_fresh = true; + break; + } if (!improved) break; score_fresh = true; if (!collapsed.empty()) { diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index d2398386d..0f91e33d9 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -37,6 +37,27 @@ oracleBremer <- function(reference, dataset, concavity = Inf) { setNames(out, names(refFreq)) } +# Enumeration oracle scoring under the inapplicable "bgs" kernel (the DEFAULT +# scoring mode). Validates the converse-constraint engine on data with +# inapplicable (`-`) tokens -- the path that runs exact_verify_sweep. +oracleBremerNA <- function(reference, dataset) { + tips <- names(dataset) + trees <- structure(lapply(allBinaryTrees(tips), RootTree, tips[1]), + class = "multiPhylo") + lengths <- TreeLength(trees, dataset, inapplicable = "bgs") + Lstar <- min(lengths) + refFreq <- SplitFrequency(reference, trees) + nSplits <- length(refFreq) + disp <- vapply(seq_along(trees), function(j) { + SplitFrequency(reference, trees[j]) + }, double(nSplits)) + dim(disp) <- c(nSplits, length(trees)) + setNames(vapply(seq_len(nSplits), function(i) { + lacking <- disp[i, ] < 0.5 + if (any(lacking)) min(lengths[lacking]) - Lstar else Inf + }, double(1)), names(refFreq)) +} + # --- method = "pool" (approximate engine) --- # A multiPhylo of several MPTs is summarised by its strict consensus; support # is calculated for the consensus's resolved clades only. @@ -178,36 +199,96 @@ test_that("constraint Bremer matches enumeration under implied weights", { expect_false(all(con == round(con))) }) -test_that("constraint is the default and dominates the pool estimate", { +test_that("constraint is the default and never looser than an uncensored pool", { + # The matrix leaves parts of the tree unresolved, so a fully-resolved reference + # has decay-0 clades that OTHER equally-optimal trees lack. A best-score pool + # therefore breaks them, making the cross-engine comparison NON-vacuous (guards + # the earlier bug where a consensus reference left every shared clade censored + # and the key assertion compared empty vectors). + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") set.seed(42) - mpts <- MaximizeParsimony(ds, maxReplicates = 8L, verbosity = 0L) - con <- Bremer(mpts, ds, maxReplicates = 12L, verbosity = 0L) # default method + mptsFull <- MaximizeParsimony(dat, collapse = FALSE, maxReplicates = 12L, + verbosity = 0L) + ref <- mptsFull[[1]] # a fully-resolved (binary) MPT + Lstar <- attr(mptsFull, "score") + + con <- Bremer(ref, dat, optimalScore = Lstar, maxReplicates = 20L, + verbosity = 0L) # default method (constraint) expect_type(con, "double") expect_true(all(con >= 0)) expect_true(all(is.finite(con))) - pool <- Bremer(mpts, ds, method = "pool", maxBremer = 8, - maxReplicates = 8L, verbosity = 0L) - # Both are upper bounds on the true decay; the directed constraint search is - # never looser than the incidental pool where the pool is uncensored. + pool <- Bremer(ref, dat, method = "pool", optimalScore = Lstar, maxBremer = 10, + maxReplicates = 50L, verbosity = 0L) shared <- intersect(names(con), names(pool)) uncensored <- shared[is.finite(pool[shared])] + expect_gt(length(uncensored), 0L) # the pool broke at least one clade + # Both are upper bounds on the true decay; the directed constraint search is + # never looser than the incidental pool where the pool is uncensored. expect_true(all(con[uncensored] <= pool[uncensored] + 1e-6)) }) -test_that("Bremer warns and stays non-negative when optimalScore is too high", { +test_that("Bremer errors when optimalScore exceeds the reference tree length", { dat <- StringToPhyDat("1100000 1110000 1111000 1111100 1001000", 1:7, byTaxon = FALSE) names(dat) <- c(LETTERS[1:6], "out") set.seed(1) mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) ref <- mpts[[1]] - # Claim an optimum 5 steps worse than reality; converse searches find shorter - # trees, so the guard must warn and clamp rather than report negative decay. + # A claimed optimum longer than the reference tree's own length is impossible + # under these scoring arguments; the scoring guard catches it before searching + # (a stale or mis-scaled L*) rather than reporting a bogus decay. + expect_error( + Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score") + 5, + maxReplicates = 20L, verbosity = 0L), + "exceeds the reference tree") +}) + +test_that("constraint Bremer matches enumeration on inapplicable (bgs) data", { + # The default scoring mode (inapplicable = "bgs") previously had no ground- + # truth oracle. Validate the negative-constraint machinery + exact_verify_sweep + # against exact enumeration on a 6-taxon matrix with inapplicable (`-`) tokens. + dat <- StringToPhyDat("111000 111000 111000 110000 000110 -00-11 -00-11", + 1:6, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:5], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 10L, verbosity = 0L) + ref <- mpts[[1]] + con <- Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score"), + maxReplicates = 30L, verbosity = 0L) + orc <- oracleBremerNA(ref, dat) + expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) +}) + +test_that("Bremer errors on a scoring-units mismatch (IW trees vs EW default)", { + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, concavity = 3, maxReplicates = 8L, verbosity = 0L) + # Trees found under implied weights; the default equal-weights Bremer would + # subtract an IW optimum from EW lengths -- caught before searching (BR-1). + expect_error(Bremer(mpts, dat, maxReplicates = 8L, verbosity = 0L), + "does not equal their stored optimal score") + # Passing the matching concavity is accepted. + ok <- Bremer(mpts, dat, concavity = 3, maxReplicates = 12L, verbosity = 0L) + expect_type(ok, "double") +}) + +test_that("Bremer ignores (with a warning) drift/anneal args in a converse search", { + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + # driftCycles is a documented forwarded arg but unsound in a converse search; + # it must be ignored (with a warning), not silently empty the pool (BR-2). expect_warning( - con <- Bremer(ref, dat, method = "constraint", - optimalScore = attr(mpts, "score") + 5, - maxReplicates = 20L, verbosity = 0L), - "shorter than the optimal score") - expect_true(all(con >= 0)) + con <- Bremer(mpts, dat, driftCycles = 5L, annealCycles = 3L, + maxReplicates = 12L, verbosity = 0L), + "always disabled") + expect_true(all(is.finite(con))) }) From 952120ef20f72c323f1a8faa72890746c2352a1b Mon Sep 17 00:00:00 2001 From: R script Date: Thu, 9 Jul 2026 19:55:31 +0100 Subject: [PATCH 03/12] feat(bremer): parallelise per-clade converse searches over a cluster `Bremer(method = "constraint")` runs one independent converse search per clade; its dominant cost is that N-search sequential fan-out. Add a `cl =` argument distributing the per-clade searches over a `parallel` cluster (each worker runs one serial search -- the negative-constraint pool guard is serial-only, so parallelism must be across clades, never within a search). The fan-out primitive `.BremerConverseScores()` seeds each clade independently (seeds drawn once from the caller's stream), so a parallel run reproduces a serial one under the same `set.seed()`, regardless of task scheduling. The default `cl = NULL` keeps the original serial path bit-for-bit; `method = "pool"` ignores `cl` (with a warning). Validated against the exhaustive enumeration oracle both serially and over a live 2-worker PSOCK cluster (identical results); the public `Bremer(cl =)` matches the serial path exactly. Adds `parallel` to Imports; completes the (previously stale) Collate with Bremer.R / BremerParallel.R / SuboptimalTrees.R. Co-Authored-By: Claude Opus 4.8 --- DESCRIPTION | 4 + NEWS.md | 5 +- R/Bremer.R | 29 ++++- R/BremerParallel.R | 59 ++++++++++ man/Bremer.Rd | 8 ++ tests/testthat/test-BremerParallel.R | 155 +++++++++++++++++++++++++++ 6 files changed, 254 insertions(+), 6 deletions(-) create mode 100644 R/BremerParallel.R create mode 100644 tests/testthat/test-BremerParallel.R diff --git a/DESCRIPTION b/DESCRIPTION index 2662de991..22d8c1cf3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -43,6 +43,7 @@ Imports: fastmatch (>= 1.1.3), graphics, grDevices, + parallel, Rcpp, Rdpack (>= 0.7), stats, @@ -91,6 +92,8 @@ Config/testthat/edition: 3 Collate: 'AdditionTree.R' 'Bootstrap.R' + 'Bremer.R' + 'BremerParallel.R' 'CharacterHierarchy.R' 'ClusterStrings.R' 'Concordance.R' @@ -119,6 +122,7 @@ Collate: 'SPR.R' 'ScoreSpectrum.R' 'Sectorial.R' + 'SuboptimalTrees.R' 'SuccessiveApproximations.R' 'TBR.R' 'TaxonInfluence.R' diff --git a/NEWS.md b/NEWS.md index 5b34d2dc3..f0dc29995 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,7 +6,10 @@ trustworthy decay indices with bounded memory. A fast, approximate `method = "pool"` instead reads a pool of suboptimal trees (an upper bound, right-censored at the sampling depth). Both engines match an exhaustive - brute-force oracle on small datasets, including under implied weights. + brute-force oracle on small datasets, including under implied weights and + inapplicable ("bgs") data. Because each clade's constraint search is + independent, they can be distributed over a `parallel` cluster with `cl =` + for a near-linear speed-up on trees with many clades. - New `SuboptimalTrees()` returns every tree the search retained within a given number of steps of the optimum, each annotated with its parsimony score, for diff --git a/R/Bremer.R b/R/Bremer.R index b2c59899c..fd3f2f0da 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -51,6 +51,12 @@ #' @param format Character specifying return format, as in [`JackLabels()`]: #' `"numeric"` (default) returns named numeric values for further analysis; #' `"character"` returns a vector shaped for `phylo$node.label`. +#' @param cl Optional \pkg{parallel} cluster (e.g. from +#' [`parallel::makeCluster()`]) over which to distribute the per-clade converse +#' searches of `method = "constraint"`. Each clade is an independent search, so +#' this gives a near-linear speed-up for trees with many clades; workers must +#' have \pkg{TreeSearch} loaded. `NULL` (default) runs serially. Ignored by +#' `method = "pool"`. #' @inheritParams MaximizeParsimony #' @param \dots Further arguments passed to [`MaximizeParsimony()`] / #' [`SuboptimalTrees()`], e.g. `maxReplicates`, `maxSeconds`, `strategy`, @@ -89,7 +95,7 @@ Bremer <- function(tree, dataset, hsj_alpha = 1.0, method = c("constraint", "pool"), maxBremer = Inf, optimalScore = NULL, - format = "numeric", ...) { + format = "numeric", cl = NULL, ...) { method <- match.arg(method) scoringArgs <- list(concavity = concavity, extended_iw = extended_iw, @@ -112,8 +118,14 @@ Bremer <- function(tree, dataset, # criteria (e.g. an equal-weights length minus an implied-weights optimum). .BremerCheckScoring(tree, dataset, scoringArgs, optimalScore) - engine <- if (method == "pool") .BremerPool else .BremerConstraint - res <- engine(ref, dataset, scoringArgs, maxBremer, list(...)) + res <- if (method == "pool") { + if (!is.null(cl)) { + warning("`cl` is ignored for method = \"pool\" (a single pooled search).") + } + .BremerPool(ref, dataset, scoringArgs, maxBremer, list(...)) + } else { + .BremerConstraint(ref, dataset, scoringArgs, maxBremer, list(...), cl = cl) + } .BremerFormat(res$bremer, res$censored, ref$reference, format) } @@ -292,7 +304,8 @@ Bremer <- function(tree, dataset, # path (drift, in-sector drift, simulated annealing) are disabled -- they would # otherwise wander onto the clade and, being worse-accepting, report it as # unsupported. Runs serially: the pool guard is on the serial search path. -.BremerConstraint <- function(ref, dataset, scoringArgs, maxBremer, dots) { +.BremerConstraint <- function(ref, dataset, scoringArgs, maxBremer, dots, + cl = NULL) { disabled <- list( driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, annealCycles = 0L @@ -372,7 +385,13 @@ Bremer <- function(tree, dataset, s } - scores <- vapply(seq_along(splitNames), processConverse, double(1)) + # Per-clade searches are independent; fan them out over `cl` when supplied. + # The default (cl = NULL) keeps the original serial vapply, bit-for-bit. + scores <- if (is.null(cl)) { + vapply(seq_along(splitNames), processConverse, double(1)) + } else { + .BremerConverseScores(length(splitNames), processConverse, cl = cl) + } if (anyNA(scores)) { warning(sum(is.na(scores)), " converse-constraint search(es) found no tree ", diff --git a/R/BremerParallel.R b/R/BremerParallel.R new file mode 100644 index 000000000..b1df4986f --- /dev/null +++ b/R/BremerParallel.R @@ -0,0 +1,59 @@ +# Parallel fan-out for the converse-constraint Bremer engine. +# +# Each clade's converse-constraint search is fully independent of the others, so +# the dominant cost of `Bremer(method = "constraint")` -- N sequential +# `MaximizeParsimony()` searches -- parallelises cleanly ACROSS clades. It must +# parallelise across clades and NOT within a single search: every converse +# search runs serially (`nThreads = 1`) because the negative-constraint pool +# guard (`TreePool::set_forbidden`) lives on the serial search path only. +# +# `.BremerConverseScores()` runs the N per-clade tasks -- optionally over a +# user-supplied `parallel` cluster -- and collects their results. The per-clade +# task (`processFn`, i.e. `.BremerConstraint()`'s `processConverse`) runs one +# converse search, maps the engine's "no tree found" sentinel to NA and verifies +# the returned tree lacks the clade, so this helper stays agnostic to that logic +# and can be validated in isolation against the enumeration oracle. +# +# Determinism: each clade is seeded independently (seeds drawn once, in order, +# from the caller's RNG stream). Seeding per CLADE rather than per WORKER makes +# the result independent of how tasks are scheduled onto workers, so a parallel +# run reproduces a serial one under the same `set.seed()` -- serial +# `MaximizeParsimony()` is itself reproducible under `set.seed()`. +# +# @param n Number of clades (tasks). +# @param processFn `function(i)` returning the (scalar) Bremer contribution of +# clade `i` -- a length, or `NA_real_`. The closure carries the reference, +# dataset and scoring/search arguments. +# @param cl A `parallel` cluster whose workers have TreeSearch loaded, or `NULL` +# for a seeded serial fan-out. +# @param seeds Optional integer per-clade RNG seeds; drawn reproducibly from the +# caller's stream when `NULL`. +# @return A numeric vector of length `n`, one entry per clade. +.BremerConverseScores <- function(n, processFn, cl = NULL, seeds = NULL) { + if (n == 0L) { + return(numeric(0)) + } + if (is.null(seeds)) { + seeds <- sample.int(.Machine$integer.max, n) + } + + # One task per clade. Seed inside the task so the outcome depends only on the + # clade index, never on which worker (or in what order) runs it. + worker <- function(i) { + set.seed(seeds[[i]]) + processFn(i) + } + + scores <- if (is.null(cl)) { + lapply(seq_len(n), worker) + } else { + if (!inherits(cl, "cluster")) { + stop("`cl` must be a `parallel` cluster or NULL.") + } + parallel::parLapply(cl, seq_len(n), worker) + } + + vapply(scores, function(s) { + if (length(s) == 0L) NA_real_ else as.numeric(s)[[1L]] + }, double(1)) +} diff --git a/man/Bremer.Rd b/man/Bremer.Rd index b098a734f..653bd92cd 100644 --- a/man/Bremer.Rd +++ b/man/Bremer.Rd @@ -18,6 +18,7 @@ Bremer( maxBremer = Inf, optimalScore = NULL, format = "numeric", + cl = NULL, ... ) } @@ -113,6 +114,13 @@ redundant search.} \code{"numeric"} (default) returns named numeric values for further analysis; \code{"character"} returns a vector shaped for \code{phylo$node.label}.} +\item{cl}{Optional \pkg{parallel} cluster (e.g. from +\code{\link[parallel:makeCluster]{parallel::makeCluster()}}) over which to distribute the per-clade converse +searches of \code{method = "constraint"}. Each clade is an independent search, so +this gives a near-linear speed-up for trees with many clades; workers must +have \pkg{TreeSearch} loaded. \code{NULL} (default) runs serially. Ignored by +\code{method = "pool"}.} + \item{\dots}{Further arguments passed to \code{\link[=MaximizeParsimony]{MaximizeParsimony()}} / \code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, e.g. \code{maxReplicates}, \code{maxSeconds}, \code{strategy}, \code{nThreads}, \code{verbosity}.} diff --git a/tests/testthat/test-BremerParallel.R b/tests/testthat/test-BremerParallel.R new file mode 100644 index 000000000..2cba07b56 --- /dev/null +++ b/tests/testthat/test-BremerParallel.R @@ -0,0 +1,155 @@ +library("TreeTools", quietly = TRUE) + +# Exact Bremer by exhaustive enumeration of every binary tree (<= 7 tips) -- the +# same gold-standard oracle used by test-Bremer.R. The parallel fan-out is +# validated to hit exactly this, seed- and schedule-independently. +allBinaryTrees <- function(tips) { + trees <- list(ape::read.tree( + text = paste0("(", tips[1], ",", tips[2], ",", tips[3], ");"))) + for (k in 4:length(tips)) { + trees <- unlist(lapply(trees, function(tr) { + AddTipEverywhere(tr, tips[k], includeRoot = FALSE) + }), recursive = FALSE) + } + trees +} + +oracleBremer <- function(reference, dataset, concavity = Inf) { + tips <- names(dataset) + trees <- structure(lapply(allBinaryTrees(tips), RootTree, tips[1]), + class = "multiPhylo") + lengths <- TreeLength(trees, dataset, concavity = concavity) + Lstar <- min(lengths) + refFreq <- SplitFrequency(reference, trees) + nSplits <- length(refFreq) + disp <- vapply(seq_along(trees), function(j) { + SplitFrequency(reference, trees[j]) + }, double(nSplits)) + dim(disp) <- c(nSplits, length(trees)) + setNames(vapply(seq_len(nSplits), function(i) { + lacking <- disp[i, ] < 0.5 + if (any(lacking)) min(lengths[lacking]) - Lstar else Inf + }, double(1)), names(refFreq)) +} + +# A per-clade task mirroring .BremerConstraint()'s processConverse: an EW, +# serial converse search (worse-accepting phases disabled) returning the raw +# optimal score of the tree forced to lack clade `i`, with the engine's +# "no tree found" sentinel mapped to NA. `TreeSearch::` qualification so the +# closure resolves the same way inside a cluster worker. +makeProcessFn <- function(dataset, splitList) { + function(i) { + s <- attr(TreeSearch::MaximizeParsimony( + dataset, collapse = TRUE, nThreads = 1L, + .negativeConstraint = splitList[[i]], + driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, + annealCycles = 0L, + maxReplicates = 25L, verbosity = 0L), "score") + if (!is.finite(s) || s < 0) NA_real_ else s + } +} + +# Shared fixture: a graded pectinate 7-taxon matrix with clades of varying +# support, so the fan-out is exercised over several distinct decay values. +bremerFixture <- function() { + dat <- StringToPhyDat( + "1100000 1110000 1111000 1111100 1100000 1110000 1111000 1111100 1001000", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + splits <- as.Splits(ref, tipLabels = names(dat)) + splitNames <- rownames(as.matrix(splits)) + list(dat = dat, ref = ref, Lstar = attr(mpts, "score"), + splits = splits, splitNames = splitNames, + splitList = lapply(seq_along(splitNames), function(i) splits[[i]])) +} + +test_that(".BremerConverseScores (serial fan-out) matches the enumeration oracle", { + skip_on_cran() + fx <- bremerFixture() + n <- length(fx$splitList) + + set.seed(7) + raw <- TreeSearch:::.BremerConverseScores( + n, makeProcessFn(fx$dat, fx$splitList), cl = NULL) + bremer <- setNames(raw - fx$Lstar, fx$splitNames) + + orc <- oracleBremer(fx$ref, fx$dat) + expect_length(bremer, length(orc)) + expect_equal(as.numeric(bremer[names(orc)]), as.numeric(orc), tolerance = 1e-6) +}) + +test_that(".BremerConverseScores is reproducible under set.seed (per-clade seeding)", { + skip_on_cran() + fx <- bremerFixture() + n <- length(fx$splitList) + pf <- makeProcessFn(fx$dat, fx$splitList) + + set.seed(7) + a <- TreeSearch:::.BremerConverseScores(n, pf, cl = NULL) + set.seed(7) + b <- TreeSearch:::.BremerConverseScores(n, pf, cl = NULL) + # Serial MaximizeParsimony is reproducible under set.seed, and the fan-out + # draws the per-clade seeds in a fixed order, so two runs are bit-identical. + expect_identical(a, b) + + # Explicit per-clade seeds bypass the RNG draw entirely. + seeds <- c(11L, 22L, 33L, 44L, 55L)[seq_len(n)] + c1 <- TreeSearch:::.BremerConverseScores(n, pf, cl = NULL, seeds = seeds) + c2 <- TreeSearch:::.BremerConverseScores(n, pf, cl = NULL, seeds = seeds) + expect_identical(c1, c2) +}) + +test_that(".BremerConverseScores handles an empty clade list", { + expect_identical( + TreeSearch:::.BremerConverseScores(0L, function(i) 0, cl = NULL), + numeric(0)) +}) + +test_that(".BremerConverseScores rejects a non-cluster cl", { + expect_error( + TreeSearch:::.BremerConverseScores(3L, function(i) 0, cl = 4L), + "must be a `parallel` cluster") +}) + +test_that(".BremerConverseScores (PSOCK cluster) matches serial and the oracle", { + skip_on_cran() + skip_on_ci() # PSOCK worker startup + load_all is slow/flaky in CI + if (!requireNamespace("parallel", quietly = TRUE) || + !requireNamespace("pkgload", quietly = TRUE)) { + skip("parallel / pkgload unavailable") + } + pkgPath <- tryCatch(pkgload::pkg_path(), error = function(e) NULL) + if (is.null(pkgPath)) skip("cannot locate package source for worker load_all") + + cl <- tryCatch(parallel::makePSOCKcluster(2L), error = function(e) NULL) + if (is.null(cl)) skip("cannot create PSOCK cluster") + on.exit(parallel::stopCluster(cl), add = TRUE) + + loaded <- tryCatch({ + parallel::clusterCall(cl, function(p) { + suppressMessages(pkgload::load_all(p, quiet = TRUE, helpers = FALSE, + attach_testthat = FALSE)) + TRUE + }, pkgPath) + }, error = function(e) NULL) + if (is.null(loaded) || !all(vapply(loaded, isTRUE, logical(1)))) { + skip("cannot load TreeSearch on cluster workers") + } + + fx <- bremerFixture() + n <- length(fx$splitList) + pf <- makeProcessFn(fx$dat, fx$splitList) + seeds <- seq_len(n) + 100L # fixed seeds -> deterministic regardless of worker + + par <- TreeSearch:::.BremerConverseScores(n, pf, cl = cl, seeds = seeds) + ser <- TreeSearch:::.BremerConverseScores(n, pf, cl = NULL, seeds = seeds) + # Same per-clade seeds -> identical searches whether run in parallel or serial. + expect_equal(par, ser, tolerance = 1e-6) + + bremer <- setNames(par - fx$Lstar, fx$splitNames) + orc <- oracleBremer(fx$ref, fx$dat) + expect_equal(as.numeric(bremer[names(orc)]), as.numeric(orc), tolerance = 1e-6) +}) From 2845b220a767f2ee44802008e8b8a9e2581b7662 Mon Sep 17 00:00:00 2001 From: R script Date: Thu, 9 Jul 2026 21:09:41 +0100 Subject: [PATCH 04/12] fix(bremer): address self-review findings on the red-team fixes A 3-agent adversarial pass over the red-team fix diff surfaced four issues (the C++ changes came back clean): - VF-1 (P2): the BR-1 scoring guard could raise a SPURIOUS fatal error on the canonical Bremer(trees, ds) call. It compared the stored L* to the reference's length under the current scoring, resolving collapsed MPTs with multi2di(random=FALSE); but an arbitrary resolution can be suboptimal, so min(refLen) > L* whenever poolMaxSize truncation leaves no cleanly-resolving member. Since TreeLength rejects non-binary trees, refLen can only be an UPPER bound on L* (exact L* is NP-hard). Reworked to a SOUND one-sided error (supplied L* exceeds an achievable length -> definitely wrong; multi2di inflation can never trigger it) plus a heuristic WARNING for the opposite direction (an implied-weights optimum scored under equal weights), which no sound score comparison can catch. - VF-2/VF-3: documented that a parallel run is reproducible under set.seed() but, because its per-clade seeding differs from the serial stream, may differ from a serial run where a search does not saturate (both valid upper bounds); added a public Bremer(cl=) vs serial equality test on saturating data (the primitive tests only covered .BremerConverseScores). - VF-4: the aggregate NA warning no longer misattributes a "returned tree still displays the clade" NA (whose per-clade warning is lost on cluster workers) to a search-budget shortfall. Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 110 +++++++++++++++++---------- man/Bremer.Rd | 7 +- tests/testthat/test-Bremer.R | 18 +++-- tests/testthat/test-BremerParallel.R | 37 +++++++++ 4 files changed, 123 insertions(+), 49 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index fd3f2f0da..38be646b1 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -56,7 +56,12 @@ #' searches of `method = "constraint"`. Each clade is an independent search, so #' this gives a near-linear speed-up for trees with many clades; workers must #' have \pkg{TreeSearch} loaded. `NULL` (default) runs serially. Ignored by -#' `method = "pool"`. +#' `method = "pool"`. Parallel runs are reproducible under `set.seed()` (each +#' clade is seeded independently); because that seeding differs from the serial +#' run's stream, a parallel result may differ from a serial one wherever a +#' converse search does not reach its true optimum -- exactly as re-running +#' serially with a different seed would -- with both remaining valid upper bounds +#' on the decay. #' @inheritParams MaximizeParsimony #' @param \dots Further arguments passed to [`MaximizeParsimony()`] / #' [`SuboptimalTrees()`], e.g. `maxReplicates`, `maxSeconds`, `strategy`, @@ -158,46 +163,63 @@ Bremer <- function(tree, dataset, # arguments than those now in effect -- a silent-units-mismatch guard. A decay # of "extra steps" is only meaningful when L* and L(not C) use the same ruler. .BremerCheckScoring <- function(tree, dataset, scoringArgs, optimalScore) { - refLen <- function() { - # Collapse-derived MPTs may be multifurcating, but TreeLength needs binary - # trees. Resolving the zero-length polytomies preserves the parsimony length - # (every resolution of an equally-optimal collapsed polytomy is itself - # optimal), and taking the minimum guards against an unlucky resolution. - resolved <- if (inherits(tree, "multiPhylo")) { - structure(lapply(tree, ape::multi2di, random = FALSE), class = "multiPhylo") - } else { - ape::multi2di(tree, random = FALSE) - } - min(suppressWarnings( - do.call(TreeLength, c(list(tree = resolved, dataset = dataset), scoringArgs)))) + suppliedLstar <- if (!is.null(optimalScore)) { + optimalScore + } else if (inherits(tree, "multiPhylo")) { + s <- attr(tree, "score") + if (!is.null(s) && is.finite(s)) s else NULL + } else { + NULL + } + if (is.null(suppliedLstar)) { + return(invisible(NULL)) + } + + # Reference length under the CURRENT scoring arguments. TreeLength requires + # binary trees and MaximizeParsimony(collapse = TRUE) can return multifurcating + # trees, so resolve the zero-length polytomies first. An arbitrary resolution + # can be SUBOPTIMAL, so `refLen` is only an UPPER bound on L* under scoringArgs + # (refLen >= L*); the checks below are chosen to stay correct under that. + resolve <- function(t) ape::multi2di(t, random = FALSE) + resolved <- if (inherits(tree, "multiPhylo")) { + structure(lapply(tree, resolve), class = "multiPhylo") + } else { + resolve(tree) + } + refLen <- min(suppressWarnings( + do.call(TreeLength, c(list(tree = resolved, dataset = dataset), scoringArgs)))) + if (!is.finite(refLen)) { + return(invisible(NULL)) } tol <- 1e-6 - if (!is.null(optimalScore)) { - # A user-asserted optimum cannot exceed the reference's own length under - # these arguments; if it does, the score or the arguments are wrong. - Lmin <- refLen() - if (is.finite(Lmin) && optimalScore > Lmin + tol) { - stop("The supplied `optimalScore` (", signif(optimalScore, 7), - ") exceeds the reference tree's length under the supplied scoring ", - "arguments (", signif(Lmin, 7), "). Check that `optimalScore` and ", - "the scoring arguments (`concavity`, `inapplicable`, ...) match ", - "those used to find the trees.", call. = FALSE) - } - } else if (inherits(tree, "multiPhylo")) { - # L* will come from the trees' stored optimal score. Because the supplied - # trees are optimal, their minimum length UNDER scoringArgs equals L* under - # scoringArgs; a disagreement means the scoring arguments differ from those - # used to find them -- fatal, as the decay would be meaningless. - attrScore <- attr(tree, "score") - if (!is.null(attrScore) && is.finite(attrScore)) { - Lmin <- refLen() - if (is.finite(Lmin) && abs(attrScore - Lmin) > tol) { - stop("The reference trees' length under the supplied scoring arguments ", - "(", signif(Lmin, 7), ") does not equal their stored optimal score ", - "(", signif(attrScore, 7), "). Pass the same `concavity`, ", - "`inapplicable` and other scoring arguments used to find the trees.", - call. = FALSE) - } + + # SOUND error: refLen is an achievable length under scoringArgs, hence an upper + # bound on the optimum. A supplied optimum that EXCEEDS it cannot be the + # optimum under these arguments -- the scoring arguments differ from those used + # to find the trees, or the score is stale. (An unlucky multi2di resolution + # only inflates refLen, so it can never trigger this -- unlike an exact-equality + # check, which spuriously errored on a collapsed reference.) + if (suppliedLstar > refLen + tol) { + stop("The reference's optimal score (", signif(suppliedLstar, 7), + ") exceeds an achievable length under the supplied scoring arguments (", + signif(refLen, 7), "). Pass the same scoring arguments (`concavity`, ", + "`inapplicable`, ...) used to find the trees.", call. = FALSE) + } + + # Heuristic warning for the opposite direction (supplied optimum well BELOW the + # reference length -- e.g. an implied-weights optimum scored under equal + # weights). Exact detection is infeasible (computing L* is NP-hard and refLen + # is only an upper bound), so WARN rather than error, and only for a multiPhylo + # (whose members are optimal, so a large gap signals a scoring-MODE mismatch, + # not a merely-suboptimal single reference tree). + if (inherits(tree, "multiPhylo")) { + gap <- refLen - suppliedLstar + if (gap > tol && gap > 0.05 * max(abs(refLen), abs(suppliedLstar), 1)) { + warning("The reference trees' optimal score (", signif(suppliedLstar, 7), + ") differs substantially from their length under the supplied ", + "scoring arguments (~", signif(refLen, 7), "). If you did not pass ", + "the same scoring arguments (`concavity`, `inapplicable`, ...) used ", + "to find the trees, the decay values will be meaningless.") } } invisible(NULL) @@ -394,8 +416,14 @@ Bremer <- function(tree, dataset, } if (anyNA(scores)) { - warning(sum(is.na(scores)), " converse-constraint search(es) found no tree ", - "lacking the clade; reported as NA. Increase `maxReplicates`.") + # Two causes map to NA: no clade-free tree was found (a budget shortfall), or + # a returned tree unexpectedly displayed the clade (an engine regression the + # BR-6 check caught). The per-clade warning that distinguishes them is lost + # on cluster workers, so keep the aggregate message honest about both. + warning(sum(is.na(scores)), " converse-constraint search(es) returned NA: ", + "either no tree lacking the clade was found (increase `maxReplicates`)", + ", or a returned tree still displayed the clade (an engine bug worth ", + "reporting).") } trueBest <- suppressWarnings(min(c(Lstar, scores), na.rm = TRUE)) diff --git a/man/Bremer.Rd b/man/Bremer.Rd index 653bd92cd..5ee617d6f 100644 --- a/man/Bremer.Rd +++ b/man/Bremer.Rd @@ -119,7 +119,12 @@ redundant search.} searches of \code{method = "constraint"}. Each clade is an independent search, so this gives a near-linear speed-up for trees with many clades; workers must have \pkg{TreeSearch} loaded. \code{NULL} (default) runs serially. Ignored by -\code{method = "pool"}.} +\code{method = "pool"}. Parallel runs are reproducible under \code{set.seed()} (each +clade is seeded independently); because that seeding differs from the serial +run's stream, a parallel result may differ from a serial one wherever a +converse search does not reach its true optimum -- exactly as re-running +serially with a different seed would -- with both remaining valid upper bounds +on the decay.} \item{\dots}{Further arguments passed to \code{\link[=MaximizeParsimony]{MaximizeParsimony()}} / \code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, e.g. \code{maxReplicates}, \code{maxSeconds}, \code{strategy}, diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 0f91e33d9..11efbe564 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -243,7 +243,7 @@ test_that("Bremer errors when optimalScore exceeds the reference tree length", { Bremer(ref, dat, method = "constraint", optimalScore = attr(mpts, "score") + 5, maxReplicates = 20L, verbosity = 0L), - "exceeds the reference tree") + "exceeds an achievable length") }) test_that("constraint Bremer matches enumeration on inapplicable (bgs) data", { @@ -263,7 +263,7 @@ test_that("constraint Bremer matches enumeration on inapplicable (bgs) data", { expect_equal(as.numeric(con[names(orc)]), as.numeric(orc), tolerance = 1e-6) }) -test_that("Bremer errors on a scoring-units mismatch (IW trees vs EW default)", { +test_that("Bremer warns on a likely scoring-units mismatch (IW trees vs EW default)", { dat <- StringToPhyDat( "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", 1:7, byTaxon = FALSE) @@ -271,11 +271,15 @@ test_that("Bremer errors on a scoring-units mismatch (IW trees vs EW default)", set.seed(1) mpts <- MaximizeParsimony(dat, concavity = 3, maxReplicates = 8L, verbosity = 0L) # Trees found under implied weights; the default equal-weights Bremer would - # subtract an IW optimum from EW lengths -- caught before searching (BR-1). - expect_error(Bremer(mpts, dat, maxReplicates = 8L, verbosity = 0L), - "does not equal their stored optimal score") - # Passing the matching concavity is accepted. - ok <- Bremer(mpts, dat, concavity = 3, maxReplicates = 12L, verbosity = 0L) + # subtract an IW optimum from EW lengths. Exact detection is infeasible (an + # arbitrary resolution only upper-bounds L*, and computing L* is NP-hard), so + # the guard WARNS on the large gap rather than erroring or silently mis-scoring + # (BR-1). + expect_warning(Bremer(mpts, dat, maxReplicates = 8L, verbosity = 0L), + "differs substantially") + # Passing the matching concavity is accepted without the mismatch warning. + ok <- suppressWarnings( + Bremer(mpts, dat, concavity = 3, maxReplicates = 12L, verbosity = 0L)) expect_type(ok, "double") }) diff --git a/tests/testthat/test-BremerParallel.R b/tests/testthat/test-BremerParallel.R index 2cba07b56..fb4a4367c 100644 --- a/tests/testthat/test-BremerParallel.R +++ b/tests/testthat/test-BremerParallel.R @@ -153,3 +153,40 @@ test_that(".BremerConverseScores (PSOCK cluster) matches serial and the oracle", orc <- oracleBremer(fx$ref, fx$dat) expect_equal(as.numeric(bremer[names(orc)]), as.numeric(orc), tolerance = 1e-6) }) + +test_that("public Bremer(cl=) matches the serial path on saturating data", { + # Exercises the WHOLE public dispatch (Bremer -> .BremerConstraint -> fan-out), + # which the .BremerConverseScores tests above do not. On 7 tips the converse + # searches saturate, so the serial and parallel paths -- despite their + # different per-clade RNG streams -- reach the same optima and must agree. + skip_on_cran() + skip_on_ci() # PSOCK worker startup + load_all is slow/flaky in CI + if (!requireNamespace("parallel", quietly = TRUE) || + !requireNamespace("pkgload", quietly = TRUE)) { + skip("parallel / pkgload unavailable") + } + pkgPath <- tryCatch(pkgload::pkg_path(), error = function(e) NULL) + if (is.null(pkgPath)) skip("cannot locate package source for worker load_all") + cl <- tryCatch(parallel::makePSOCKcluster(2L), error = function(e) NULL) + if (is.null(cl)) skip("cannot create PSOCK cluster") + on.exit(parallel::stopCluster(cl), add = TRUE) + loaded <- tryCatch( + parallel::clusterCall(cl, function(p) { + suppressMessages(pkgload::load_all(p, quiet = TRUE, helpers = FALSE, + attach_testthat = FALSE)) + TRUE + }, pkgPath), error = function(e) NULL) + if (is.null(loaded) || !all(vapply(loaded, isTRUE, logical(1)))) { + skip("cannot load TreeSearch on cluster workers") + } + + fx <- bremerFixture() + set.seed(1) + mpts <- MaximizeParsimony(fx$dat, maxReplicates = 8L, verbosity = 0L) + set.seed(5) + ser <- Bremer(mpts, fx$dat, maxReplicates = 20L, verbosity = 0L) + set.seed(5) + par <- Bremer(mpts, fx$dat, cl = cl, maxReplicates = 20L, verbosity = 0L) + expect_equal(unname(par), unname(ser), tolerance = 1e-6) + expect_identical(names(par), names(ser)) +}) From c5f700b5a4f7c5b69ef85ce3b716d55b75847afb Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 06:05:03 +0100 Subject: [PATCH 05/12] fix(bremer): guard converse sector/prune-reinsert against rebuilding the forbidden clade The reduced-dataset sector solve (rss_search/xss_search) and the prune-reinsert backbone TBR run constraint-blind, so a reinsertion could rebuild a forbidden clade even though the pre-move tree lacked it. The pool backstop still guaranteed soundness (a clade-displaying tree is rejected at pool insert), but the LIVE search tree could strand on the clade with no improving escape move -- wasting the converse replicate and inflating/NA-ing the reported decay on datasets with >= 12 tips (where sectors engage). This is structurally invisible to the <= 7-tip enumeration oracle, so the gap was never caught. Mirror the positive post-hoc constraint check with a negative displays_forbidden_clade revert in rss_search/xss_search and the prune-reinsert accept path, gated on neg_active so mainline (positive/unconstrained) searches are byte-identical. Also fail loudly on a consNegSplitMatrix column-count mismatch (previously a silent unconstrained converse search -> decay reported as 0) and correct a misleading displays_forbidden_clade doc comment. Red-team round-2 findings A-11/A-12/A-15/A-17. Validated: mainline sector/prune-reinsert/positive-constraint/driven suites unchanged (613 assertions), <= 7-tip enumeration oracle exact, and a new >= 12-tip sector test asserts the converse trees genuinely lack the clade. Reach-benefit magnitude is a Hamilton follow-up (not measurable on the exhaustive oracle). Co-Authored-By: Claude Opus 4.8 --- src/ts_constraint.h | 12 ++++-- src/ts_prune_reinsert.cpp | 11 +++++- src/ts_rcpp.cpp | 10 ++++- src/ts_sector.cpp | 32 +++++++++++++++ tests/testthat/test-Bremer.R | 75 ++++++++++++++++++++++++++++++++++++ 5 files changed, 133 insertions(+), 7 deletions(-) diff --git a/src/ts_constraint.h b/src/ts_constraint.h index 1bc494c35..f8eb74fc5 100644 --- a/src/ts_constraint.h +++ b/src/ts_constraint.h @@ -153,10 +153,14 @@ std::vector compute_node_tips(const TreeState& tree, int n_words); void add_negative_constraint( ConstraintData& cd, const int* neg_matrix, int n_neg, int n_tips); -// Returns true if the tree currently displays any forbidden clade. The tree -// is rooted on tip 0, so a canonical (tip-0-outside) bipartition is displayed -// iff some internal node's subtree tip mask equals it -- the same test -// map_constraint_nodes() uses to locate a positive split. +// Returns true if the tree currently displays any forbidden clade. A +// bipartition is displayed iff some internal node's subtree tip mask equals it +// on EITHER side, so each node mask is canonicalized to tip-0-outside form +// (flipped when tip 0 is inside) before comparison. This bilateral +// canonicalization is REQUIRED because the internal search tree is not +// guaranteed to be rooted on tip 0 -- unlike map_constraint_nodes(), which does +// a one-sided (already-tip-0-outside) compare for a positive split and so must +// not be assumed bilateral-safe. bool displays_forbidden_clade(const TreeState& tree, const ConstraintData& cd); // Repair constraint violations by minimal SPR moves. diff --git a/src/ts_prune_reinsert.cpp b/src/ts_prune_reinsert.cpp index 68a3ef092..22f75d9f2 100644 --- a/src/ts_prune_reinsert.cpp +++ b/src/ts_prune_reinsert.cpp @@ -654,12 +654,19 @@ PruneReinsertResult prune_reinsert_search( // 7. Accept or revert double new_score = score_tree(tree, ds); - if (new_score < current_score - 1e-10) { + // Negative (converse/Bremer) constraint: the reduced-backbone TBR (step 4) + // runs cd-blind, so expand_and_reinsert can rebuild a forbidden clade. + // Reject such a tree even if it scores better, reverting to the (clade-free) + // pre-prune backup -- otherwise the replicate can strand on the clade. + // Soundness is already guaranteed by the pool backstop; this preserves reach. + bool neg_violated = + cd && cd->neg_active && displays_forbidden_clade(tree, *cd); + if (new_score < current_score - 1e-10 && !neg_violated) { current_score = new_score; result.best_score = new_score; ++result.n_improvements; } else { - tree = backup; // revert + tree = backup; // revert (worse, or would display the forbidden clade) // Re-sync constraint metadata after topology revert. // Same bug class as F-015 (ratchet), F-016 (NNI-perturb). if (cd) update_constraint(tree, *cd); diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index b182afaa4..924df6b0c 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -1477,7 +1477,15 @@ static ts::ConstraintData build_constraint_from_r( if (consNegSplitMatrix.isNotNull()) { IntegerMatrix nsm(consNegSplitMatrix.get()); int n_neg = nsm.nrow(); - if (n_neg > 0 && nsm.ncol() == n_tips) { + if (n_neg > 0) { + // A column-count mismatch would otherwise silently drop the negative + // constraint, running an UNCONSTRAINED converse search and reporting a + // clade's Bremer support as 0 -- a silent wrong answer. Fail loudly. + if (nsm.ncol() != n_tips) { + Rcpp::stop("consNegSplitMatrix has %d column(s) but the dataset has %d " + "tip(s); a negative-constraint matrix needs one column per " + "tip.", nsm.ncol(), n_tips); + } ts::add_negative_constraint(cd, INTEGER(nsm), n_neg, n_tips); } } diff --git a/src/ts_sector.cpp b/src/ts_sector.cpp index a55134b42..7c2c78156 100644 --- a/src/ts_sector.cpp +++ b/src/ts_sector.cpp @@ -1504,6 +1504,16 @@ SectorResult rss_search(TreeState& tree, DataSet& ds, // (T-S6c micro-bank). static const bool _sect_debug = std::getenv("TS_SECT_DEBUG") != nullptr; bool constrained = cd && cd->active && cd->has_posthoc; + // Negative (converse/Bremer) constraint: the sector-internal reduced-dataset + // solve (search_sector/reinsert_sector) is constraint-blind, so a reinsertion + // can rebuild a forbidden clade even though the pre-sector tree lacked it. + // Reject such a reinsertion post-hoc (mirrors `constrained` above), keeping the + // live tree -- and hence every downstream phase -- in the space of trees + // lacking the clade. Soundness is already guaranteed by the pool backstop; + // this preserves REACH (an un-reverted reinsertion can strand the replicate on + // the clade with no improving escape move). Not covered by the <=7-tip oracle + // (sectors need >=12 tips). + bool neg_constrained = cd && cd->neg_active; // Seed RNG (from R in serial mode, from thread-local in parallel mode) std::mt19937 rng = ts::make_rng(); @@ -1747,6 +1757,16 @@ SectorResult rss_search(TreeState& tree, DataSet& ds, score_tree(tree, ds); continue; } + // Post-hoc negative (converse) constraint check: reject a reinsertion that + // rebuilt a forbidden clade. restore_clade undoes only the sector clade, + // which is the sole change since the pre-sector tree, so the reverted tree + // is clade-free. + if (neg_constrained && displays_forbidden_clade(tree, *cd)) { + restore_clade(tree, snap); + tree.build_postorder(); + score_tree(tree, ds); + continue; + } bool kept; if (new_score < result.best_score) { @@ -1891,6 +1911,11 @@ SectorResult xss_search(TreeState& tree, DataSet& ds, } bool constrained = cd && cd->active && cd->has_posthoc; + // Negative (converse/Bremer) constraint: the reduced-dataset sector solve is + // constraint-blind, so reject post-hoc any reinsertion that rebuilds a + // forbidden clade (see rss_search for the full rationale). REACH-preserving; + // soundness is already guaranteed by the pool backstop. + bool neg_constrained = cd && cd->neg_active; for (int round = 0; round < params.xss_rounds; ++round) { double score_before_round = result.best_score; @@ -1949,6 +1974,13 @@ SectorResult xss_search(TreeState& tree, DataSet& ds, score_tree(tree, ds); continue; } + // Post-hoc negative (converse) constraint check (see rss_search). + if (neg_constrained && displays_forbidden_clade(tree, *cd)) { + restore_clade(tree, snap); + tree.build_postorder(); + score_tree(tree, ds); + continue; + } if (new_score < result.best_score) { result.total_steps_saved += diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 11efbe564..c4249deb3 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -296,3 +296,78 @@ test_that("Bremer ignores (with a warning) drift/anneal args in a converse searc "always disabled") expect_true(all(is.finite(con))) }) + +# --- method = "constraint" on >= 12 tips: exercises the SECTORIAL search path --- +# rss_search / xss_search engage once a tree has >= 12 tips. Their reduced- +# dataset sector solve is constraint-blind, so a reinsertion could rebuild a +# forbidden clade; the fix rejects such a reinsertion post-hoc (a reach guard -- +# the pool backstop already guaranteed soundness). This whole path is +# structurally invisible to the <= 7-tip enumeration oracle above (sectors need +# >= 12 tips), so it needs its own coverage. + +# A graded pectinate matrix on `nIn` ingroup taxa + one outgroup: nested clades +# {t1..tk} (k = 2..nIn-1), each duplicated for graded support. +sectorFixture <- function(nIn = 13) { + tips <- c(paste0("t", seq_len(nIn)), "out") + nTip <- length(tips) + mkChar <- function(k) paste0(c(rep("1", k), rep("0", nTip - k)), collapse = "") + chars <- unlist(lapply(rep(2:(nIn - 1L), each = 2L), mkChar)) + dat <- StringToPhyDat(paste(chars, collapse = " "), seq_len(nTip), + byTaxon = FALSE) + names(dat) <- tips + dat +} + +test_that("converse search on >= 12 tips (sectors engaged) returns clade-free trees", { + skip_on_cran() + dat <- sectorFixture(13) # 14 tips -> sectorial search engages + expect_gte(length(dat), 12L) + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 6L, verbosity = 0L) + ref <- mpts[[1]] + splits <- as.Splits(ref, tipLabels = names(dat)) + splitNames <- rownames(as.matrix(splits)) + expect_gt(length(splitNames), 0L) + + # Directly drive the negative-constraint sector path for a spread of clades and + # confirm each returned tree genuinely LACKS its forbidden clade. The pool + # backstop guarantees this regardless of the sector fix, so it is a robust + # soundness assertion; the fix's REACH benefit (fewer stranded replicates) is a + # Hamilton-scale measurement, not asserted here. + probe <- unique(round(seq(1, length(splitNames), + length.out = min(4L, length(splitNames))))) + for (i in probe) { + set.seed(100L + i) + res <- MaximizeParsimony(dat, collapse = TRUE, nThreads = 1L, + .negativeConstraint = splits[[i]], + driftCycles = 0L, sectorGoDrift = 0L, + sectorDriftCycles = 0L, annealCycles = 0L, + maxReplicates = 12L, verbosity = 0L) + tr <- if (inherits(res, "phylo")) res else res[[1L]] + disp <- SplitFrequency(ref, structure(list(tr), class = "multiPhylo")) + expect_lt(unname(disp[splitNames[i]]), 0.5) # returned tree lacks the clade + } +}) + +test_that("constraint Bremer on >= 12 tips is non-negative and <= pool", { + skip_on_cran() + dat <- sectorFixture(13) + set.seed(2) + mpts <- MaximizeParsimony(dat, maxReplicates = 6L, verbosity = 0L) + ref <- mpts[[1]] + Lstar <- attr(mpts, "score") + + con <- Bremer(ref, dat, method = "constraint", optimalScore = Lstar, + maxReplicates = 10L, verbosity = 0L) + expect_type(con, "double") + expect_true(any(is.finite(con))) # the sector path did not strand + expect_true(all(con[is.finite(con)] >= 0)) # correct L* (no negative decay) + + pool <- Bremer(ref, dat, method = "pool", optimalScore = Lstar, maxBremer = 12, + maxReplicates = 30L, verbosity = 0L) + shared <- intersect(names(con), names(pool)) + uncensored <- shared[is.finite(pool[shared]) & is.finite(con[shared])] + # Both are upper bounds on the true decay; the directed constraint search is + # never looser than the incidental pool where both are finite. + expect_true(all(con[uncensored] <= pool[uncensored] + 1e-6)) +}) From 8a0f56b126f914f2a226c24e4da52e491f79d042 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 06:09:51 +0100 Subject: [PATCH 06/12] fix(bremer): validate optimalScore, make censored attr consistent, strip pool-size overrides Three red-team round-2 R findings: - optimalScore = NA_real_ slipped past the is.null() "not supplied" sentinel and crashed the one-sided scoring guard with "missing value where TRUE/FALSE needed"; now rejected up front with a clear message (B-8). - Bremer(format = "character") dropped the `censored` attribute entirely when no clade was censored, while the numeric format kept an all-FALSE vector -- so downstream code assuming its presence broke on the character path whenever nothing happened to be censored. The character format now always attaches a node.label-shaped censored vector (B-9). - SuboptimalTrees() stripped colliding `control`/`collapse` but not the raw `poolSuboptimal`/`poolMaxSize` SearchControl fields, which would silently win over the values set from `maxSuboptimal`/`maxPool` via MaximizeParsimony()'s dots-override-control merge; now warned and ignored (B-10). Tests added for each; Bremer / SuboptimalTrees / BremerParallel suites green. Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 18 +++++++++-- R/SuboptimalTrees.R | 13 ++++++++ tests/testthat/test-Bremer.R | 43 +++++++++++++++++++++++++++ tests/testthat/test-SuboptimalTrees.R | 14 +++++++++ 4 files changed, 86 insertions(+), 2 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index 38be646b1..a69282e23 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -103,6 +103,16 @@ Bremer <- function(tree, dataset, format = "numeric", cl = NULL, ...) { method <- match.arg(method) + # `optimalScore = NULL` is the documented "not supplied" sentinel. A non-NULL + # value must be a single finite number: NA_real_ in particular would otherwise + # slip past the is.null() checks and crash the one-sided scoring guard with + # "missing value where TRUE/FALSE needed". + if (!is.null(optimalScore) && + (length(optimalScore) != 1L || !is.numeric(optimalScore) || + !is.finite(optimalScore))) { + stop("`optimalScore` must be a single finite number, or NULL.") + } + scoringArgs <- list(concavity = concavity, extended_iw = extended_iw, xpiwe_r = xpiwe_r, xpiwe_max_f = xpiwe_max_f, hierarchy = hierarchy, inapplicable = inapplicable, @@ -243,8 +253,12 @@ Bremer <- function(tree, dataset, ret[idx] <- as.character(values) # Inf renders as the literal "Inf"; preserve the censoring flag as an # attribute so pool-method callers can still distinguish "> maxBremer" - # (censored) from a genuine value (the numeric format carries it too). - if (length(censored) && any(censored)) { + # (censored) from a genuine value. Attach it whenever there ARE entries + # (not only when something is censored), so the character format carries + # the same fixed-length `censored` attribute the numeric format does -- + # otherwise downstream code assuming its presence breaks on the character + # path whenever no clade happens to be censored. + if (length(censored)) { censVec <- logical(reference[["Nnode"]]) censVec[idx] <- as.logical(censored) attr(ret, "censored") <- censVec diff --git a/R/SuboptimalTrees.R b/R/SuboptimalTrees.R index 89d267617..79ca0cc0c 100644 --- a/R/SuboptimalTrees.R +++ b/R/SuboptimalTrees.R @@ -85,6 +85,19 @@ SuboptimalTrees <- function(dataset, tree = NULL, } dots[["collapse"]] <- NULL } + # `poolSuboptimal`/`poolMaxSize` are the raw SearchControl() fields that + # `maxSuboptimal`/`maxPool` set below. Passed via `...` they would flow into + # MaximizeParsimony()'s dots-override-control merge and SILENTLY win over the + # constructed control (last-writer-wins), defeating `maxSuboptimal`/`maxPool`. + # Strip them with a warning, directing the caller to the blessed arguments. + managed <- intersect(c("poolSuboptimal", "poolMaxSize"), names(dots)) + if (length(managed)) { + warning("Set the suboptimal-pool size via `maxSuboptimal` / `maxPool`, not ", + "`", paste(managed, collapse = "` / `"), + "`; the latter would silently override the former and ", + if (length(managed) > 1L) "are" else "is", " ignored.") + dots[managed] <- NULL + } control <- SearchControl(poolSuboptimal = maxSuboptimal, poolMaxSize = maxPool) diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index c4249deb3..8475d48b6 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -371,3 +371,46 @@ test_that("constraint Bremer on >= 12 tips is non-negative and <= pool", { # never looser than the incidental pool where both are finite. expect_true(all(con[uncensored] <= pool[uncensored] + 1e-6)) }) + +# --- input validation / format consistency (red-team round 2) --- + +test_that("Bremer rejects a non-finite or non-scalar optimalScore (B-8)", { + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + # NA_real_ slips past is.null() and previously crashed the one-sided scoring + # guard with "missing value where TRUE/FALSE needed"; now a clear error. + expect_error( + Bremer(mpts, dat, optimalScore = NA_real_, maxReplicates = 6L, + verbosity = 0L), + "single finite number") + expect_error( + Bremer(mpts, dat, optimalScore = c(1, 2), maxReplicates = 6L, + verbosity = 0L), + "single finite number") +}) + +test_that("Bremer(format='character') always carries a censored attribute (B-9)", { + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, collapse = FALSE, maxReplicates = 8L, + verbosity = 0L) + ref <- mpts[[1]] + Lstar <- attr(mpts, "score") + # method = "constraint" never censors, so this is the "nothing censored" path + # where the character format previously dropped the attribute entirely while + # the numeric format kept an all-FALSE vector. + lab <- Bremer(ref, dat, optimalScore = Lstar, format = "character", + maxReplicates = 12L, verbosity = 0L) + cens <- attr(lab, "censored") + expect_false(is.null(cens)) # attribute is present... + expect_type(cens, "logical") + expect_length(cens, ref[["Nnode"]]) # ...node.label-shaped... + expect_false(any(cens)) # ...and correctly all-FALSE. + # The numeric format for the same call also carries a censored attribute. + num <- Bremer(ref, dat, optimalScore = Lstar, + maxReplicates = 12L, verbosity = 0L) + expect_false(is.null(attr(num, "censored"))) +}) diff --git a/tests/testthat/test-SuboptimalTrees.R b/tests/testthat/test-SuboptimalTrees.R index e4940d788..4d88e7f03 100644 --- a/tests/testthat/test-SuboptimalTrees.R +++ b/tests/testthat/test-SuboptimalTrees.R @@ -63,3 +63,17 @@ test_that("SuboptimalTrees() validates its arguments", { expect_error(SuboptimalTrees(ds, maxPool = 0L), "must be a single positive integer") }) + +test_that("SuboptimalTrees() warns and ignores raw poolSuboptimal/poolMaxSize (B-10)", { + # These raw SearchControl fields would silently override the values set from + # maxSuboptimal/maxPool via MaximizeParsimony()'s dots-override-control merge. + set.seed(3418) + expect_warning( + SuboptimalTrees(ds, maxSuboptimal = 2, poolSuboptimal = 9, + maxReplicates = 2L, targetHits = 1L, verbosity = 0L), + "maxSuboptimal") + expect_warning( + SuboptimalTrees(ds, maxPool = 100L, poolMaxSize = 5L, + maxReplicates = 2L, targetHits = 1L, verbosity = 0L), + "maxPool") +}) From 22e6646f9efe313817cd1c3705ff940fa4a633b8 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 06:29:02 +0100 Subject: [PATCH 07/12] fix(bremer): make the scoring-mismatch guard a reliable warning, not a fuzzy error The scoring-units guard now checks a supplied optimalScore against the reference's length under the scoring arguments in effect and WARNS (in either direction) on any material difference, then accepts the score and proceeds -- per the package author's decision to trust the user's scoring choice and only alert them in case they forgot to pass the scoring arguments used to find the trees. This replaces the previous one-sided 5%-of-length heuristic, which had a P1 blind spot: two criteria that place L* close together (a finder reproduced inapplicable="bgs" vs "missing" giving an identical length) slipped under the threshold and produced a silently mis-scored decay vector. The tight two-sided tolerance now catches a small but genuine mismatch; taking the minimum over the supplied MPTs absorbs multi2di resolution slop, so a consistent call does not false-warn (verified across the EW / IW / bgs test corpus). The former hard error on an "impossible" L* is downgraded to the same accept-and-proceed warning. Red-team round-2 finding B-7 / TA-6. Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 68 ++++++++++++++++-------------------- tests/testthat/test-Bremer.R | 50 ++++++++++++++++++-------- 2 files changed, 67 insertions(+), 51 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index a69282e23..34a6589ea 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -169,9 +169,15 @@ Bremer <- function(tree, dataset, Lstar = Lstar) } -# Error if the reference's optimal score L* was measured under different scoring -# arguments than those now in effect -- a silent-units-mismatch guard. A decay -# of "extra steps" is only meaningful when L* and L(not C) use the same ruler. +# Sanity-check a SUPPLIED optimal score L* against the reference's length under +# the scoring arguments now in effect. A decay of "extra steps" is only +# meaningful when L* and L(not C) are measured with the same ruler, so a +# difference usually means the reference's L* was computed under different +# scoring arguments (e.g. `Bremer()` left at the default equal weights for trees +# found under implied weights). We trust the user to keep the scoring mode +# consistent, so this only WARNS (never errors) and Bremer proceeds with the +# supplied score -- the warning is a safety net for a forgotten scoring argument, +# not a veto on a deliberately different analysis. .BremerCheckScoring <- function(tree, dataset, scoringArgs, optimalScore) { suppliedLstar <- if (!is.null(optimalScore)) { optimalScore @@ -185,11 +191,13 @@ Bremer <- function(tree, dataset, return(invisible(NULL)) } - # Reference length under the CURRENT scoring arguments. TreeLength requires - # binary trees and MaximizeParsimony(collapse = TRUE) can return multifurcating - # trees, so resolve the zero-length polytomies first. An arbitrary resolution - # can be SUBOPTIMAL, so `refLen` is only an UPPER bound on L* under scoringArgs - # (refLen >= L*); the checks below are chosen to stay correct under that. + # Length of the reference under the CURRENT scoring arguments -- cheap (scoring, + # not searching). TreeLength requires binary trees and collapse = TRUE can + # return multifurcating ones, so resolve the zero-length polytomies first; an + # arbitrary resolution only ever LENGTHENS a tree, and we take the minimum over + # the supplied trees, so refLen is a tight estimate of the reference's + # achievable length (== L* whenever the reference is optimal under these + # arguments -- as an MPT set is, when the scoring modes match). resolve <- function(t) ape::multi2di(t, random = FALSE) resolved <- if (inherits(tree, "multiPhylo")) { structure(lapply(tree, resolve), class = "multiPhylo") @@ -201,36 +209,22 @@ Bremer <- function(tree, dataset, if (!is.finite(refLen)) { return(invisible(NULL)) } - tol <- 1e-6 - # SOUND error: refLen is an achievable length under scoringArgs, hence an upper - # bound on the optimum. A supplied optimum that EXCEEDS it cannot be the - # optimum under these arguments -- the scoring arguments differ from those used - # to find the trees, or the score is stale. (An unlucky multi2di resolution - # only inflates refLen, so it can never trigger this -- unlike an exact-equality - # check, which spuriously errored on a collapsed reference.) - if (suppliedLstar > refLen + tol) { - stop("The reference's optimal score (", signif(suppliedLstar, 7), - ") exceeds an achievable length under the supplied scoring arguments (", - signif(refLen, 7), "). Pass the same scoring arguments (`concavity`, ", - "`inapplicable`, ...) used to find the trees.", call. = FALSE) - } - - # Heuristic warning for the opposite direction (supplied optimum well BELOW the - # reference length -- e.g. an implied-weights optimum scored under equal - # weights). Exact detection is infeasible (computing L* is NP-hard and refLen - # is only an upper bound), so WARN rather than error, and only for a multiPhylo - # (whose members are optimal, so a large gap signals a scoring-MODE mismatch, - # not a merely-suboptimal single reference tree). - if (inherits(tree, "multiPhylo")) { - gap <- refLen - suppliedLstar - if (gap > tol && gap > 0.05 * max(abs(refLen), abs(suppliedLstar), 1)) { - warning("The reference trees' optimal score (", signif(suppliedLstar, 7), - ") differs substantially from their length under the supplied ", - "scoring arguments (~", signif(refLen, 7), "). If you did not pass ", - "the same scoring arguments (`concavity`, `inapplicable`, ...) used ", - "to find the trees, the decay values will be meaningless.") - } + # A material difference in EITHER direction signals a probable scoring-mode + # mismatch: a supplied L* ABOVE an achievable length cannot be the optimum + # here, and one well BELOW it means the reference is far from optimal under + # these arguments (e.g. an implied-weights optimum scored under equal weights). + # Unlike the previous one-sided 5%-of-length heuristic, a tight tolerance also + # catches a small but genuine mismatch that two criteria happen to place close + # together. Accept and proceed regardless (trusting the user's choice). + tol <- 1e-6 * max(1, abs(refLen), abs(suppliedLstar)) + if (abs(suppliedLstar - refLen) > tol) { + warning("The supplied optimal score (", signif(suppliedLstar, 7), + ") differs from the reference tree length under the supplied scoring ", + "arguments (", signif(refLen, 7), "). If you did not pass the same ", + "scoring arguments (`concavity`, `inapplicable`, ...) used to find ", + "the trees, the decay values mix two optimality criteria and are ", + "meaningless. Proceeding with the supplied score.") } invisible(NULL) } diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 8475d48b6..48a44cd7e 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -229,21 +229,24 @@ test_that("constraint is the default and never looser than an uncensored pool", expect_true(all(con[uncensored] <= pool[uncensored] + 1e-6)) }) -test_that("Bremer errors when optimalScore exceeds the reference tree length", { +test_that("Bremer warns (does not error) on an inconsistent optimalScore, then proceeds", { dat <- StringToPhyDat("1100000 1110000 1111000 1111100 1001000", 1:7, byTaxon = FALSE) names(dat) <- c(LETTERS[1:6], "out") set.seed(1) mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) ref <- mpts[[1]] - # A claimed optimum longer than the reference tree's own length is impossible - # under these scoring arguments; the scoring guard catches it before searching - # (a stale or mis-scaled L*) rather than reporting a bogus decay. - expect_error( - Bremer(ref, dat, method = "constraint", - optimalScore = attr(mpts, "score") + 5, - maxReplicates = 20L, verbosity = 0L), - "exceeds an achievable length") + # A supplied L* inconsistent with the reference length under these scoring + # arguments is WARNED, not errored: Bremer trusts the user's scoring choice and + # proceeds (they may deliberately want a different analysis). A value BELOW the + # reference length exercises the warning without triggering the adopt-shorter-L* + # path (which would add a second, unrelated warning). + expect_warning( + con <- Bremer(ref, dat, method = "constraint", + optimalScore = attr(mpts, "score") - 3, + maxReplicates = 20L, verbosity = 0L), + "differs from the reference") + expect_type(con, "double") }) test_that("constraint Bremer matches enumeration on inapplicable (bgs) data", { @@ -271,18 +274,37 @@ test_that("Bremer warns on a likely scoring-units mismatch (IW trees vs EW defau set.seed(1) mpts <- MaximizeParsimony(dat, concavity = 3, maxReplicates = 8L, verbosity = 0L) # Trees found under implied weights; the default equal-weights Bremer would - # subtract an IW optimum from EW lengths. Exact detection is infeasible (an - # arbitrary resolution only upper-bounds L*, and computing L* is NP-hard), so - # the guard WARNS on the large gap rather than erroring or silently mis-scoring - # (BR-1). + # subtract an IW optimum from EW lengths. The scoring guard checks the supplied + # L* against the reference's length under the current arguments and WARNS on the + # difference (accepting and proceeding), rather than erroring or silently + # mis-scoring (BR-1/B-7). expect_warning(Bremer(mpts, dat, maxReplicates = 8L, verbosity = 0L), - "differs substantially") + "differs from the reference") # Passing the matching concavity is accepted without the mismatch warning. ok <- suppressWarnings( Bremer(mpts, dat, concavity = 3, maxReplicates = 12L, verbosity = 0L)) expect_type(ok, "double") }) +test_that("the scoring-mismatch warning catches a small gap the 5% band missed (B-7)", { + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) + ref <- mpts[[1]] + Lstar <- attr(mpts, "score") + # A gap of 0.3 is < 5% of the reference length, so the old one-sided + # 5%-of-length heuristic would NOT have warned (the P1 blind spot two criteria + # placing L* close together); the tight two-sided tolerance now does. Below + # the reference length, so the adopt-shorter-L* path does not add a 2nd warning. + expect_warning( + Bremer(ref, dat, method = "constraint", optimalScore = Lstar - 0.3, + maxReplicates = 15L, verbosity = 0L), + "differs from the reference") +}) + test_that("Bremer ignores (with a warning) drift/anneal args in a converse search", { dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) names(dat) <- c(LETTERS[1:6], "out") From 9ee1d3cf41f46d52682d07c38c750a75a4f1bd88 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 06:44:40 +0100 Subject: [PATCH 08/12] feat(bremer): record scoring provenance on MaximizeParsimony output; check it exactly in Bremer A saved optimal score is only meaningful alongside the conditions it was optimal under. MaximizeParsimony() now stamps a `scoring` attribute -- a canonical signature of the scoring arguments (concavity, extended_iw, xpiwe_r, xpiwe_max_f, hierarchy presence, inapplicable, hsj_alpha) -- on its returned multiPhylo. Bremer() compares that signature EXACTLY against the scoring arguments in effect: a mismatch (e.g. implied-weights trees passed to a default equal-weights Bremer) warns and proceeds (trusting the user's choice); a match is silent. This closes the B-7/TA-6 blind spot at the root -- it needs no re-scoring, so it has neither the resolution slop nor the near-equal-length hole of the length heuristic, which is retained only as a fallback for a bare optimalScore or a reference built outside MaximizeParsimony. Additive attribute; validated across EW / IW / profile / hsj / xform / hierarchy modes (537 assertions in MaximizeParsimony / Jackknife / scoring-mode suites unchanged) and the Bremer suite (provenance match = no warning, mismatch = warning, length fallback intact). Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 23 +++++++++++++ R/MaximizeParsimony.R | 65 ++++++++++++++++++++++++++++++++++++ tests/testthat/test-Bremer.R | 36 +++++++++++++++++--- 3 files changed, 119 insertions(+), 5 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index 34a6589ea..4db64cc39 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -179,6 +179,29 @@ Bremer <- function(tree, dataset, # supplied score -- the warning is a safety net for a forgotten scoring argument, # not a veto on a deliberately different analysis. .BremerCheckScoring <- function(tree, dataset, scoringArgs, optimalScore) { + # Exact check: a MaximizeParsimony() result records the scoring conditions it + # is optimal under (attr "scoring"). When present, compare that signature + # DIRECTLY to the arguments now in effect -- the meaningful test, since a saved + # optimal score is only interpretable alongside the conditions it was found + # under. It needs no re-scoring, so it has neither the resolution-slop nor the + # near-equal-length blind spot of the length fallback below. + recorded <- attr(tree, "scoring", exact = TRUE) + if (!is.null(recorded)) { + current <- do.call(.ScoringSignature, scoringArgs) + if (!.ScoringSignatureMatch(recorded, current)) { + warning("The reference trees were found under ", .DescribeScoring(recorded), + " but Bremer() is scoring with ", .DescribeScoring(current), + ". The decay values mix two optimality criteria and are ", + "meaningless unless you pass the same scoring arguments ", + "(`concavity`, `inapplicable`, ...) used to find the trees. ", + "Proceeding with the supplied score.") + } + return(invisible(NULL)) + } + + # Fallback (no recorded signature: a bare optimalScore, a single-tree reference, + # or a tree built outside MaximizeParsimony). Compare the supplied optimal + # score to the reference's length under the current scoring arguments. suppliedLstar <- if (!is.null(optimalScore)) { optimalScore } else if (inherits(tree, "multiPhylo")) { diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index ec1b04770..14ac93a64 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -617,6 +617,10 @@ #' attributes: #' \describe{ #' \item{`score`}{Best parsimony score.} +#' \item{`scoring`}{A list recording the scoring conditions the score is +#' optimal under (`concavity`, `inapplicable`, ...), so a saved score +#' remains interpretable. [`Bremer()`][Bremer] checks it against its own +#' scoring arguments and warns on a mismatch.} #' \item{`scores`}{Present only when `collapse = FALSE`: a numeric vector of #' the parsimony score of each returned tree, aligned with the returned #' `multiPhylo` (the same values are attached as a `score` attribute on @@ -710,6 +714,16 @@ MaximizeParsimony <- function( # the only reliable read. userSetReps <- !missing(maxReplicates) + # Capture the scoring conditions the result is optimal under, BEFORE any + # normalization (e.g. inapplicable aliasing, profile -> Inf). Attached to the + # returned trees so a saved optimal score is interpretable -- a bare score is + # meaningless without knowing the criterion it was found under, and Bremer() + # checks this signature against its own scoring arguments. + scoringSignature <- .ScoringSignature( + concavity = concavity, extended_iw = extended_iw, xpiwe_r = xpiwe_r, + xpiwe_max_f = xpiwe_max_f, hierarchy = hierarchy, + inapplicable = inapplicable, hsj_alpha = hsj_alpha) + # --- Set targetHits default if not provided --- if (is.null(targetHits)) { targetHits <- max(10L, as.integer(NTip(dataset) / 5)) @@ -1219,6 +1233,7 @@ MaximizeParsimony <- function( structure( outTrees, score = result$best_score, + scoring = scoringSignature, scores = perTreeScores, replicates = result$replicates, hits_to_best = result$hits_to_best, @@ -1235,6 +1250,56 @@ MaximizeParsimony <- function( ) } +# Canonical scoring identity: the arguments that define the optimality criterion, +# normalized so aliases compare equal. Recorded on a MaximizeParsimony() result +# (attr "scoring") and checked by Bremer() -- a saved optimal score is only +# meaningful alongside the conditions it was optimal under. +.ScoringSignature <- function(concavity = Inf, extended_iw = TRUE, xpiwe_r = 0.5, + xpiwe_max_f = 5, hierarchy = NULL, + inapplicable = "bgs", hsj_alpha = 1.0) { + inap <- tolower(as.character(inapplicable)[[1]]) + if (inap == "brazeau") inap <- "bgs" # documented alias for the bgs kernel + list( + concavity = if (is.numeric(concavity)) concavity else + tolower(as.character(concavity)[[1]]), + extended_iw = isTRUE(extended_iw), + xpiwe_r = as.double(xpiwe_r), + xpiwe_max_f = as.double(xpiwe_max_f), + # Presence only: comparing full hierarchy objects is heavy and brittle, and a + # present-vs-absent difference is the practically important mismatch. + hierarchy = !is.null(hierarchy), + inapplicable = inap, + hsj_alpha = as.double(hsj_alpha) + ) +} + +# TRUE if two scoring signatures denote the same optimality criterion. +.ScoringSignatureMatch <- function(a, b, tol = 1e-8) { + numMatch <- function(x, y) { + if (is.numeric(x) && is.numeric(y)) { + (is.infinite(x) && is.infinite(y) && sign(x) == sign(y)) || + (is.finite(x) && is.finite(y) && abs(x - y) <= tol) + } else { + identical(as.character(x), as.character(y)) + } + } + numMatch(a[["concavity"]], b[["concavity"]]) && + identical(a[["extended_iw"]], b[["extended_iw"]]) && + isTRUE(abs(a[["xpiwe_r"]] - b[["xpiwe_r"]]) <= tol) && + isTRUE(abs(a[["xpiwe_max_f"]] - b[["xpiwe_max_f"]]) <= tol) && + identical(a[["hierarchy"]], b[["hierarchy"]]) && + identical(a[["inapplicable"]], b[["inapplicable"]]) && + isTRUE(abs(a[["hsj_alpha"]] - b[["hsj_alpha"]]) <= tol) +} + +# Compact human-readable rendering of a scoring signature for diagnostics. +.DescribeScoring <- function(sig) { + cn <- sig[["concavity"]] + weight <- if (is.numeric(cn) && is.infinite(cn)) "equal weights" else + paste0("concavity = ", if (is.numeric(cn)) signif(cn, 4) else cn) + paste0(weight, ", inapplicable = \"", sig[["inapplicable"]], "\"") +} + #' @rdname MaximizeParsimony #' @usage MaximizeParsimony2(...) diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 48a44cd7e..478226222 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -274,18 +274,44 @@ test_that("Bremer warns on a likely scoring-units mismatch (IW trees vs EW defau set.seed(1) mpts <- MaximizeParsimony(dat, concavity = 3, maxReplicates = 8L, verbosity = 0L) # Trees found under implied weights; the default equal-weights Bremer would - # subtract an IW optimum from EW lengths. The scoring guard checks the supplied - # L* against the reference's length under the current arguments and WARNS on the - # difference (accepting and proceeding), rather than erroring or silently - # mis-scoring (BR-1/B-7). + # subtract an IW optimum from EW lengths. The recorded scoring signature makes + # this an EXACT criterion mismatch (concavity 3 vs Inf), so it WARNS (accepting + # and proceeding) rather than erroring or silently mis-scoring (BR-1/B-7). expect_warning(Bremer(mpts, dat, maxReplicates = 8L, verbosity = 0L), - "differs from the reference") + "found under") # Passing the matching concavity is accepted without the mismatch warning. ok <- suppressWarnings( Bremer(mpts, dat, concavity = 3, maxReplicates = 12L, verbosity = 0L)) expect_type(ok, "double") }) +test_that("MaximizeParsimony records a scoring signature Bremer trusts (B-7 provenance)", { + dat <- StringToPhyDat( + "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", + 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, maxReplicates = 8L, verbosity = 0L) # EW default + sig <- attr(mpts, "scoring") + expect_false(is.null(sig)) + expect_true(is.infinite(sig[["concavity"]])) # equal weights + expect_identical(sig[["inapplicable"]], "bgs") # default kernel + + # Matching scoring arguments -> the exact provenance check passes -> no warning + # (and, being all defaults here, this is the canonical Bremer() call). + expect_warning( + con <- Bremer(mpts, dat, maxReplicates = 12L, verbosity = 0L), NA) + expect_true(all(is.finite(con))) + + # An IW analysis under matching concavity also matches its own signature. + set.seed(1) + mptsIW <- MaximizeParsimony(dat, concavity = 5, maxReplicates = 8L, + verbosity = 0L) + expect_equal(attr(mptsIW, "scoring")[["concavity"]], 5) + expect_warning( + Bremer(mptsIW, dat, concavity = 5, maxReplicates = 8L, verbosity = 0L), NA) +}) + test_that("the scoring-mismatch warning catches a small gap the 5% band missed (B-7)", { dat <- StringToPhyDat( "1110000 1110000 1110000 1100000 0001100 0001100 0000110 1111111", From 2cc63117682018c79b7c6f95d3d59b232480985c Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 06:47:13 +0100 Subject: [PATCH 09/12] test(bremer): cover zero-clade early return and pool-eviction score alignment Two of the deferred round-2 test-adequacy gaps, added with fast searches: - TA-9: Bremer() on a fully-unresolved (star) reference warns and returns an empty numeric (the "nothing to calculate" early return, previously untested). - TA-5: SuboptimalTrees() under a tiny maxPool exercises the diversity-eviction path and asserts the `scores` attribute (and each tree's `score`) stays aligned with the evicted-down set -- the exact truncation mechanism that hid VF-1. Co-Authored-By: Claude Opus 4.8 --- tests/testthat/test-Bremer.R | 12 ++++++++++++ tests/testthat/test-SuboptimalTrees.R | 13 +++++++++++++ 2 files changed, 25 insertions(+) diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 478226222..713b53e36 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -439,6 +439,18 @@ test_that("Bremer rejects a non-finite or non-scalar optimalScore (B-8)", { "single finite number") }) +test_that("Bremer on a fully-unresolved reference warns and returns nothing (TA-9)", { + dat <- StringToPhyDat("111000 110000 001100", 1:6, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:5], "out") + star <- ape::stree(length(dat), type = "star", tip.label = names(dat)) + # No resolved internal clades -> nothing to calculate; early return, no search. + expect_warning( + res <- Bremer(star, dat, verbosity = 0L), + "no resolved internal clades") + expect_type(res, "double") + expect_length(res, 0L) +}) + test_that("Bremer(format='character') always carries a censored attribute (B-9)", { dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) names(dat) <- c(LETTERS[1:6], "out") diff --git a/tests/testthat/test-SuboptimalTrees.R b/tests/testthat/test-SuboptimalTrees.R index 4d88e7f03..741e33b46 100644 --- a/tests/testthat/test-SuboptimalTrees.R +++ b/tests/testthat/test-SuboptimalTrees.R @@ -64,6 +64,19 @@ test_that("SuboptimalTrees() validates its arguments", { "must be a single positive integer") }) +test_that("SuboptimalTrees() keeps scores aligned under pool eviction (TA-5)", { + # A tiny maxPool forces the diversity-eviction path; the scores attribute -- and + # each tree's own score attribute -- must stay aligned with the evicted-down set. + set.seed(3418) + trees <- SuboptimalTrees(ds, maxSuboptimal = 8, maxPool = 5L, + maxReplicates = 4L, targetHits = 1L, verbosity = 0L) + expect_true(length(trees) <= 5L) + scores <- attr(trees, "scores") + expect_length(scores, length(trees)) + expect_equal(vapply(trees, attr, double(1), "score"), scores) # alignment held + expect_equal(scores, TreeLength(trees, ds), tolerance = 1e-8) # independent re-score +}) + test_that("SuboptimalTrees() warns and ignores raw poolSuboptimal/poolMaxSize (B-10)", { # These raw SearchControl fields would silently override the values set from # maxSuboptimal/maxPool via MaximizeParsimony()'s dots-override-control merge. From 41556cdc79d47b52c762fe25a99d4c54828037b8 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 07:06:28 +0100 Subject: [PATCH 10/12] test(bremer): close TA-3/TA-4 via a search-injection seam The aggregate-NA warning (VF-4) and the per-clade "returned tree still displays the clade" check (BR-6) guard against engine states the C++ move-guard + pool backstop are designed never to reach, so an ordinary search cannot provoke them. Add a proper test seam instead of leaving them untested: - Extract .BremerProcessResult() -- the per-clade result handler -- and unit-test its NA sentinels (non-finite / negative score) and its displays-clade warning path directly (TA-4 / BR-6). - Add an internal .runConverse injection argument to .BremerConstraint() so a stub engine can drive every clade to NA, exercising the aggregate NA warning and the NA assembly end-to-end (TA-3 / VF-4). Pure refactor of the real path (processConverse now delegates to the extracted helper); the full Bremer + BremerParallel suites are unchanged. Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 58 +++++++++++++++++++++-------------- tests/testthat/test-Bremer.R | 59 ++++++++++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+), 23 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index 4db64cc39..2622a311d 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -345,6 +345,32 @@ Bremer <- function(tree, dataset, list(bremer = bremer, censored = censored) } +# Convert one converse-search result `out` = list(score, tree) to a decay- +# eligible score, or NA. Two ways to NA: a non-finite/negative score is the +# engine's "no tree lacking the clade found" sentinel (empty pool -> budget +# shortfall); a returned tree that STILL displays the clade is an engine +# regression the C++ guard + pool backstop should preclude -- surface it as NA + +# a warning rather than a silently deflated decay (the worst outcome for a +# published statistic). +#' @importFrom TreeTools SplitFrequency +.BremerProcessResult <- function(out, refReference, splitName) { + s <- out$score + if (!is.finite(s) || s < 0) { + return(NA_real_) + } + if (!is.null(out$tree)) { + disp <- SplitFrequency(refReference, + structure(list(out$tree), class = "multiPhylo")) + # Index by the reference node number (robust to split ordering). + if (isTRUE(unname(disp[splitName]) >= 0.5)) { + warning("The converse search for clade ", splitName, + " returned a tree that still displays it; reported as NA.") + return(NA_real_) + } + } + s +} + # Rigorous Bremer by converse-constraint search: for each clade, the shortest # tree forced to LACK it. Uses the negative-constraint engine wired into # MaximizeParsimony() via the internal `.negativeConstraint` argument. @@ -358,7 +384,7 @@ Bremer <- function(tree, dataset, # otherwise wander onto the clade and, being worse-accepting, report it as # unsupported. Runs serially: the pool guard is on the serial search path. .BremerConstraint <- function(ref, dataset, scoringArgs, maxBremer, dots, - cl = NULL) { + cl = NULL, .runConverse = NULL) { disabled <- list( driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, annealCycles = 0L @@ -409,33 +435,19 @@ Bremer <- function(tree, dataset, } list(score = attr(res, "score"), tree = tree) } + # Test seam: an injected search function (same `(negSplit) -> list(score, tree)` + # contract) bypasses the real engine, so the NA / displays-clade handling and + # the aggregate NA warning can be exercised deterministically. + if (!is.null(.runConverse)) { + runConverse <- .runConverse + } splitNames <- ref$splitNames # Per-clade work unit (also the unit of parallelism for the cluster fan-out). processConverse <- function(i) { - out <- runConverse(ref$splits[[i]]) - s <- out$score - # A negative score is the engine's "no tree found" sentinel (empty pool): - # the converse search failed to locate any tree lacking the clade. - if (!is.finite(s) || s < 0) { - return(NA_real_) - } - # Verify the returned tree really lacks the clade. The C++ guard + pool - # backstop already guarantee this; checking here means a future engine - # regression surfaces as a visible NA rather than a silently deflated decay - # for a well-supported clade (the worst outcome for a published statistic). - if (!is.null(out$tree)) { - disp <- SplitFrequency(ref$reference, - structure(list(out$tree), class = "multiPhylo")) - # Index by the reference node number (robust to split ordering). - if (isTRUE(unname(disp[splitNames[i]]) >= 0.5)) { - warning("The converse search for clade ", splitNames[i], - " returned a tree that still displays it; reported as NA.") - return(NA_real_) - } - } - s + .BremerProcessResult(runConverse(ref$splits[[i]]), ref$reference, + splitNames[i]) } # Per-clade searches are independent; fan them out over `cl` when supplied. diff --git a/tests/testthat/test-Bremer.R b/tests/testthat/test-Bremer.R index 713b53e36..c1c2c78ee 100644 --- a/tests/testthat/test-Bremer.R +++ b/tests/testthat/test-Bremer.R @@ -474,3 +474,62 @@ test_that("Bremer(format='character') always carries a censored attribute (B-9)" maxReplicates = 12L, verbosity = 0L) expect_false(is.null(attr(num, "censored"))) }) + +# --- engine-failure handling, via a test seam (no real search needed) --- +# These deterministically drive the NA / clade-displaying paths that the real +# engine (C++ guard + pool backstop) is designed never to reach, so they cannot +# be provoked by an ordinary search. + +test_that(".BremerProcessResult maps engine sentinels to NA (BR-6/TA-4)", { + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, collapse = FALSE, maxReplicates = 8L, + verbosity = 0L) + ref <- mpts[[1]] + sn <- rownames(as.matrix(as.Splits(ref, tipLabels = names(dat))))[[1]] + + # "No tree lacking the clade found" sentinel (negative / non-finite) -> NA. + expect_identical( + TreeSearch:::.BremerProcessResult(list(score = -1, tree = NULL), ref, sn), + NA_real_) + expect_identical( + TreeSearch:::.BremerProcessResult(list(score = NA_real_, tree = NULL), ref, sn), + NA_real_) + + # A returned tree that STILL displays the clade (the reference does) -> the + # BR-6 defence-in-depth check warns and returns NA rather than a deflated decay. + expect_warning( + r <- TreeSearch:::.BremerProcessResult(list(score = 5, tree = ref), ref, sn), + "still displays it") + expect_identical(r, NA_real_) + + # A genuinely clade-free tree (a star displays no clade) -> score passes through. + star <- ape::stree(length(dat), type = "star", tip.label = names(dat)) + expect_identical( + TreeSearch:::.BremerProcessResult(list(score = 5, tree = star), ref, sn), 5) +}) + +test_that(".BremerConstraint emits the aggregate NA warning (VF-4/TA-3)", { + dat <- StringToPhyDat("1111000 1111000 1100000", 1:7, byTaxon = FALSE) + names(dat) <- c(LETTERS[1:6], "out") + set.seed(1) + mpts <- MaximizeParsimony(dat, collapse = FALSE, maxReplicates = 8L, + verbosity = 0L) + ref <- TreeSearch:::.BremerReference(mpts, dat, + optimalScore = attr(mpts, "score")) + expect_gte(length(ref$splitNames), 2L) + sa <- list(concavity = Inf, extended_iw = TRUE, xpiwe_r = 0.5, xpiwe_max_f = 5, + hierarchy = NULL, inapplicable = "bgs", hsj_alpha = 1.0) + + # Inject an engine stub returning the "no tree found" sentinel for every clade, + # so every converse search NAs and the aggregate warning fires. Lstar is + # supplied via the reference, so no real search runs. + stubNA <- function(negSplit) list(score = -1, tree = NULL) + expect_warning( + res <- TreeSearch:::.BremerConstraint(ref, dat, sa, Inf, list(), + cl = NULL, .runConverse = stubNA), + "returned NA") + expect_length(res$bremer, length(ref$splitNames)) + expect_true(all(is.na(res$bremer))) +}) From 855c19940d9e0f2e4d49e0cc10f3f76d02dfff41 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 14:02:59 +0100 Subject: [PATCH 11/12] fix(bremer): green the R CMD check (examples, spelling, docs, deps) The sense-check run on the PR surfaced 2 ERRORs + 2 WARNINGs; all four are doc/hygiene, not code: - examples ERROR: the Bremer example called LabelSplits() with no active device (edgelabels -> `last_plot.phylo` not found), and keyed decay onto trees[[1]] rather than the consensus it is computed against. Plot and label the consensus reference so node numbers align and a device exists. - tests ERROR (spelling): add bgs, bipartitions, clade's, run's, search's to inst/WORDLIST (new terms in the Bremer/SuboptimalTrees docs). - Rd WARNING: document the internal .negativeConstraint argument of MaximizeParsimony() (roxygen requires every \usage arg documented). - tests WARNING: declare pkgload in Suggests (test-BremerParallel.R uses it for the PSOCK worker load_all()). Co-Authored-By: Claude Opus 4.8 --- DESCRIPTION | 1 + R/Bremer.R | 7 +++++-- R/MaximizeParsimony.R | 3 +++ inst/WORDLIST | 5 +++++ man/Bremer.Rd | 7 +++++-- man/MaximizeParsimony.Rd | 4 ++++ 6 files changed, 23 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8d70dbab4..435557dfe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -61,6 +61,7 @@ Suggests: knitr, MaxMin, phangorn (>= 2.2.1), + pkgload, PlotTools, promises, protoclust, diff --git a/R/Bremer.R b/R/Bremer.R index 2622a311d..29b047552 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -83,8 +83,11 @@ #' decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) #' decay #' -#' # Annotate a tree -#' TreeTools::LabelSplits(trees[[1]], decay) +#' # Annotate the reference tree (the strict consensus of the optimal trees, +#' # whose node numbers key `decay`) +#' reference <- ape::consensus(trees, p = 1) +#' plot(reference) +#' TreeTools::LabelSplits(reference, decay) #' } #' @references \insertAllCited{} #' @template MRS diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 305940aa4..f7636fdff 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -600,6 +600,9 @@ #' collapse, so an enforced-but-unsupported clade (a zero-length branch) stays #' visible (a constraint encodes external evidence the matrix does not #' capture); unsupported non-constraint branches still collapse. +#' @param .negativeConstraint Internal. A splits object (or `NULL`) naming +#' clades that returned trees must *not* display; used by [`Bremer()`] to run +#' converse-constraint searches. Not intended for direct use. #' @param ... Backward compatibility. #' #' @return A `multiPhylo` object containing the best tree(s) found, with diff --git a/inst/WORDLIST b/inst/WORDLIST index 16b7e17b2..2dd985f81 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -196,12 +196,15 @@ abcd ac aculiferan archaeopriapulid +bgs +bipartitions bristleworm cdef cdot cf characterwise cla +clade's cleanup codecov colourblind @@ -272,6 +275,8 @@ rerootings rescored rescoring rightarrow +run's +search's secondaries sim softmax diff --git a/man/Bremer.Rd b/man/Bremer.Rd index 3fc73cc7f..268193df9 100644 --- a/man/Bremer.Rd +++ b/man/Bremer.Rd @@ -182,8 +182,11 @@ trees <- MaximizeParsimony(dataset, maxReplicates = 8, verbosity = 0) decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) decay -# Annotate a tree -TreeTools::LabelSplits(trees[[1]], decay) +# Annotate the reference tree (the strict consensus of the optimal trees, +# whose node numbers key `decay`) +reference <- ape::consensus(trees, p = 1) +plot(reference) +TreeTools::LabelSplits(reference, decay) } } \references{ diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 8b1532ec2..3221a2f00 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -116,6 +116,10 @@ in any output tree. Constraint searches are supported natively: all tree rearrangements are filtered to respect the constraint topology.} +\item{.negativeConstraint}{Internal. A splits object (or \code{NULL}) naming +clades that returned trees must \emph{not} display; used by \code{\link[=Bremer]{Bremer()}} to run +converse-constraint searches. Not intended for direct use.} + \item{strategy}{Character: named strategy preset controlling the search heuristic parameters. Presets: \describe{ From 9a67c04970f63d0b1e0d25bc36b7c25f69b78689 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 16:09:00 +0100 Subject: [PATCH 12/12] fix(bremer): make the runnable example cheap, seeded, and hang-proof The R 4.1 CI leg hung 81 min (-> job cancelled) in the Bremer `\donttest` example, while the same example completes in ~40s on the release/devel legs. Root cause is not slowness (4x would be ~3 min): the example had no `set.seed()`, so each leg's RNG stream yields a different consensus/clade set, and on R 4.1 one clade's converse-constraint search (method="constraint" fans out to one search per clade, ~15 here) fails to terminate under the default unbounded budget. Make the executed example structurally unable to hang: - `set.seed(0)` for a reproducible search. - Run the fast `method = "pool"` path (a single bounded SuboptimalTrees search, no converse fan-out) as the `\donttest` demo. - Demote the rigorous default `method = "constraint"` to `\dontrun`, and model `maxSeconds` there as the recommended bound. Bremer's correctness is covered by the oracle + 126 testthat assertions, so the example is illustrative, not test surface. The underlying converse-search non-termination under an unbounded budget is tracked separately (non-blocking). Co-Authored-By: Claude Opus 4.8 --- R/Bremer.R | 13 ++++++++++++- man/Bremer.Rd | 13 ++++++++++++- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/R/Bremer.R b/R/Bremer.R index 29b047552..8bb0ce8ac 100644 --- a/R/Bremer.R +++ b/R/Bremer.R @@ -79,8 +79,12 @@ #' data("inapplicable.phyData", package = "TreeSearch") #' dataset <- inapplicable.phyData[["Vinther2008"]] #' \donttest{ +#' # `set.seed()` makes the heuristic search reproducible +#' set.seed(0) #' trees <- MaximizeParsimony(dataset, maxReplicates = 8, verbosity = 0) -#' decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) +#' +#' # A fast, approximate decay index read from the suboptimal-tree pool +#' decay <- Bremer(trees, dataset, method = "pool", verbosity = 0) #' decay #' #' # Annotate the reference tree (the strict consensus of the optimal trees, @@ -89,6 +93,13 @@ #' plot(reference) #' TreeTools::LabelSplits(reference, decay) #' } +#' +#' \dontrun{ +#' # The default `method = "constraint"` is rigorous but much slower: it runs +#' # one converse-constraint search per clade. Bound each search with +#' # `maxSeconds`, and/or distribute the searches over a cluster with `cl`. +#' Bremer(trees, dataset, maxReplicates = 8, maxSeconds = 10, verbosity = 0) +#' } #' @references \insertAllCited{} #' @template MRS #' @seealso diff --git a/man/Bremer.Rd b/man/Bremer.Rd index 268193df9..bbcb001a0 100644 --- a/man/Bremer.Rd +++ b/man/Bremer.Rd @@ -178,8 +178,12 @@ will be meaningless; they default to equal-weights Fitch parsimony. data("inapplicable.phyData", package = "TreeSearch") dataset <- inapplicable.phyData[["Vinther2008"]] \donttest{ +# `set.seed()` makes the heuristic search reproducible +set.seed(0) trees <- MaximizeParsimony(dataset, maxReplicates = 8, verbosity = 0) -decay <- Bremer(trees, dataset, maxReplicates = 8, verbosity = 0) + +# A fast, approximate decay index read from the suboptimal-tree pool +decay <- Bremer(trees, dataset, method = "pool", verbosity = 0) decay # Annotate the reference tree (the strict consensus of the optimal trees, @@ -188,6 +192,13 @@ reference <- ape::consensus(trees, p = 1) plot(reference) TreeTools::LabelSplits(reference, decay) } + +\dontrun{ +# The default `method = "constraint"` is rigorous but much slower: it runs +# one converse-constraint search per clade. Bound each search with +# `maxSeconds`, and/or distribute the searches over a cluster with `cl`. +Bremer(trees, dataset, maxReplicates = 8, maxSeconds = 10, verbosity = 0) +} } \references{ \insertAllCited{}