diff --git a/DESCRIPTION b/DESCRIPTION index 5c3f5f71b..435557dfe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -47,6 +47,7 @@ Imports: fastmatch (>= 1.1.3), graphics, grDevices, + parallel, Rcpp, Rdpack (>= 0.7), stats, @@ -60,6 +61,7 @@ Suggests: knitr, MaxMin, phangorn (>= 2.2.1), + pkgload, PlotTools, promises, protoclust, @@ -95,6 +97,8 @@ Config/testthat/edition: 3 Collate: 'AdditionTree.R' 'Bootstrap.R' + 'Bremer.R' + 'BremerParallel.R' 'CharacterHierarchy.R' 'ClusterStrings.R' 'Concordance.R' @@ -122,6 +126,7 @@ Collate: 'recode_hierarchy.R' 'SPR.R' 'ScoreSpectrum.R' + 'SuboptimalTrees.R' 'SuccessiveApproximations.R' 'TBR.R' 'TaxonInfluence.R' diff --git a/NAMESPACE b/NAMESPACE index a7039a09b..df08cee37 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -25,6 +25,7 @@ export(.NonDuplicateRoot) export(.UniqueExceptHits) export(AdditionTree) export(BootstrapTree) +export(Bremer) export(Carter1) export(CharacterHierarchy) export(CharacterLength) @@ -111,6 +112,7 @@ export(SingleCharData) export(SingleCharMorphy) export(StepInformation) export(StopUnlessBifurcating) +export(SuboptimalTrees) export(Suboptimality) export(SuccessiveApproximations) export(SuccessiveWeights) diff --git a/NEWS.md b/NEWS.md index 8fa93f015..d904bf541 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,22 @@ # 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 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 + 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..8bb0ce8ac --- /dev/null +++ b/R/Bremer.R @@ -0,0 +1,497 @@ +#' 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`. +#' @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"`. 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`, +#' `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{ +#' # `set.seed()` makes the heuristic search reproducible +#' set.seed(0) +#' trees <- MaximizeParsimony(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, +#' # whose node numbers key `decay`) +#' 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{} +#' @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", 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, + 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)) + } + + # 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) + + 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) +} + +# 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) +} + +# 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) { + # 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")) { + s <- attr(tree, "score") + if (!is.null(s) && is.finite(s)) s else NULL + } else { + NULL + } + if (is.null(suppliedLstar)) { + return(invisible(NULL)) + } + + # 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") + } else { + resolve(tree) + } + refLen <- min(suppressWarnings( + do.call(TreeLength, c(list(tree = resolved, dataset = dataset), scoringArgs)))) + if (!is.finite(refLen)) { + return(invisible(NULL)) + } + + # 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) +} + +# 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) + # 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. 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 + } + } + 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) +} + +# 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. +# +# 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, + cl = NULL, .runConverse = NULL) { + disabled <- list( + driftCycles = 0L, sectorGoDrift = 0L, sectorDriftCycles = 0L, + annealCycles = 0L + ) + # 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"))] + + # 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") + } + + # 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) + 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) + } + # 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) { + .BremerProcessResult(runConverse(ref$splits[[i]]), ref$reference, + splitNames[i]) + } + + # 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)) { + # 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)) + 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/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/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index cb12c07a5..f7636fdff 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. @@ -568,12 +600,25 @@ #' 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 #' 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 +#' 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.} @@ -641,6 +686,7 @@ MaximizeParsimony <- function( inapplicable = "bgs", hsj_alpha = 1.0, constraint, + .negativeConstraint = NULL, strategy = "auto", maxReplicates = 96L, targetHits = NULL, @@ -664,6 +710,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)) @@ -950,6 +1006,27 @@ 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)) { + # 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 + if (verbosity > 0L) { + cli_alert_info("Forbidding {nrow(negMatrix)} clade{?s}") + } + } + } + # --- Profile parsimony: extract info_amounts --- profileArgs <- list() if (useProfile) { @@ -1041,6 +1118,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) @@ -1089,9 +1171,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 --- @@ -1114,6 +1213,8 @@ MaximizeParsimony <- function( structure( outTrees, score = result$best_score, + scoring = scoringSignature, + scores = perTreeScores, replicates = result$replicates, hits_to_best = result$hits_to_best, n_topologies = nTopologies, @@ -1129,6 +1230,57 @@ 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"]], "\"") +} + + #' Launch tree search graphical user interface #' #' Opens a "shiny" app for interactive parsimony tree search and results diff --git a/R/SuboptimalTrees.R b/R/SuboptimalTrees.R new file mode 100644 index 000000000..bf5fbaf16 --- /dev/null +++ b/R/SuboptimalTrees.R @@ -0,0 +1,108 @@ +#' 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 +#' @family tree scoring +#' @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.") + } + + # 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 + } + # `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) + + do.call(MaximizeParsimony, + c(list(dataset = dataset, tree = tree, control = control, + collapse = FALSE), dots)) +} 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/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 new file mode 100644 index 000000000..bbcb001a0 --- /dev/null +++ b/man/Bremer.Rd @@ -0,0 +1,225 @@ +% 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", + cl = NULL, + ... +) +} +\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{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"}. 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}, +\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{ +# `set.seed()` makes the heuristic search reproducible +set.seed(0) +trees <- MaximizeParsimony(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, +# whose node numbers key `decay`) +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{} +} +\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[=MostContradictedFreq]{MostContradictedFreq()}}, +\code{\link[=PaintCharacters]{PaintCharacters()}}, +\code{\link[=PresCont]{PresCont()}}, +\code{\link[=Resample]{Resample()}}, +\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/CharacterHierarchy.Rd b/man/CharacterHierarchy.Rd index e24157dba..949175e7f 100644 --- a/man/CharacterHierarchy.Rd +++ b/man/CharacterHierarchy.Rd @@ -56,6 +56,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/CharacterLength.Rd b/man/CharacterLength.Rd index 2c4e07d0e..6c1753a7e 100644 --- a/man/CharacterLength.Rd +++ b/man/CharacterLength.Rd @@ -54,6 +54,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 9d6a0831b..e1e8b1234 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[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, \code{\link[=Resample]{Resample()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \concept{split support functions} diff --git a/man/ExpectedLength.Rd b/man/ExpectedLength.Rd index 02b7c5201..b7e5a2d88 100644 --- a/man/ExpectedLength.Rd +++ b/man/ExpectedLength.Rd @@ -46,6 +46,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/HierarchyFromNames.Rd b/man/HierarchyFromNames.Rd index 436e0f622..547baa762 100644 --- a/man/HierarchyFromNames.Rd +++ b/man/HierarchyFromNames.Rd @@ -39,6 +39,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/JackLabels.Rd b/man/JackLabels.Rd index dc7475f5a..d40a774ca 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[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, \code{\link[=Resample]{Resample()}}, -\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 c91a958ee..ca59ec122 100644 --- a/man/Jackknife.Rd +++ b/man/Jackknife.Rd @@ -140,13 +140,15 @@ engine. } Other split support functions: +\code{\link[=Bremer]{Bremer()}}, \code{\link[=ConcordanceTable]{ConcordanceTable()}}, \code{\link[=JackLabels]{JackLabels()}}, \code{\link[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, \code{\link[=Resample]{Resample()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} Other custom search functions: \code{\link[=BootstrapTree]{BootstrapTree()}}, diff --git a/man/LengthAdded.Rd b/man/LengthAdded.Rd index 8da377732..ab31f7fab 100644 --- a/man/LengthAdded.Rd +++ b/man/LengthAdded.Rd @@ -76,6 +76,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 02e6a38b2..3221a2f00 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -16,6 +16,7 @@ MaximizeParsimony( inapplicable = "bgs", hsj_alpha = 1, constraint, + .negativeConstraint = NULL, strategy = "auto", maxReplicates = 96L, targetHits = NULL, @@ -115,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{ @@ -234,6 +239,16 @@ A \code{multiPhylo} object containing the best tree(s) found, with attributes: \describe{ \item{\code{score}}{Best parsimony score.} +\item{\code{scoring}}{A list recording the scoring conditions the score is +optimal under (\code{concavity}, \code{inapplicable}, ...), so a saved score +remains interpretable. \code{\link[=Bremer]{Bremer()}} checks it against its own +scoring arguments and warns on a mismatch.} +\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.} @@ -362,6 +377,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/MinimumLength.Rd b/man/MinimumLength.Rd index ebfdce34d..fe849fd35 100644 --- a/man/MinimumLength.Rd +++ b/man/MinimumLength.Rd @@ -90,6 +90,7 @@ Other tree scoring: \code{\link[=MaximizeParsimony]{MaximizeParsimony()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/MostContradictedFreq.Rd b/man/MostContradictedFreq.Rd index 0ccbe0d70..e35ae1023 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[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, \code{\link[=Resample]{Resample()}}, -\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 49500254f..d21856b78 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[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PresCont]{PresCont()}}, \code{\link[=Resample]{Resample()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \concept{split support functions} diff --git a/man/ParsSim.Rd b/man/ParsSim.Rd index 6c8fc8d27..82f365d2a 100644 --- a/man/ParsSim.Rd +++ b/man/ParsSim.Rd @@ -124,6 +124,7 @@ Other tree scoring: \code{\link[=MaximizeParsimony]{MaximizeParsimony()}}, \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/PresCont.Rd b/man/PresCont.Rd index d1a1be7a3..e4ae9a180 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[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=Resample]{Resample()}}, -\code{\link{SiteConcordance}} +\code{\link{SiteConcordance}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}} } \author{ \href{https://smithlabdurham.github.io/}{Martin R. Smith} diff --git a/man/RecodeHierarchy.Rd b/man/RecodeHierarchy.Rd index afe971a06..f53392aa3 100644 --- a/man/RecodeHierarchy.Rd +++ b/man/RecodeHierarchy.Rd @@ -74,6 +74,7 @@ Other tree scoring: \code{\link[=MaximizeParsimony]{MaximizeParsimony()}}, \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/Resample.Rd b/man/Resample.Rd index e1ad9fea4..48fbcec3f 100644 --- a/man/Resample.Rd +++ b/man/Resample.Rd @@ -191,12 +191,14 @@ JackLabels(result[[1]], jackTrees) } \seealso{ 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()}} } \concept{split support functions} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 7dc208b54..c36a84cde 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[=MostContradictedFreq]{MostContradictedFreq()}}, \code{\link[=PaintCharacters]{PaintCharacters()}}, \code{\link[=PresCont]{PresCont()}}, -\code{\link[=Resample]{Resample()}} +\code{\link[=Resample]{Resample()}}, +\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..3bb1c9405 --- /dev/null +++ b/man/SuboptimalTrees.Rd @@ -0,0 +1,102 @@ +% 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[=MostContradictedFreq]{MostContradictedFreq()}}, +\code{\link[=PaintCharacters]{PaintCharacters()}}, +\code{\link[=PresCont]{PresCont()}}, +\code{\link[=Resample]{Resample()}}, +\code{\link{SiteConcordance}} + +Other tree scoring: +\code{\link[=CharacterHierarchy]{CharacterHierarchy()}}, +\code{\link[=CharacterLength]{CharacterLength()}}, +\code{\link[=ExpectedLength]{ExpectedLength()}}, +\code{\link[=HierarchyFromNames]{HierarchyFromNames()}}, +\code{\link[=IWScore]{IWScore()}}, +\code{\link[=LengthAdded]{LengthAdded()}}, +\code{\link[=MaximizeParsimony]{MaximizeParsimony()}}, +\code{\link[=MinimumLength]{MinimumLength()}}, +\code{\link[=ParsSim]{ParsSim()}}, +\code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=TaxonInfluence]{TaxonInfluence()}}, +\code{\link[=TreeScore]{TreeScore()}}, +\code{\link[=WideSample]{WideSample()}} +} +\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} +\concept{tree scoring} diff --git a/man/TaxonInfluence.Rd b/man/TaxonInfluence.Rd index 1cb91151f..282ab2e55 100644 --- a/man/TaxonInfluence.Rd +++ b/man/TaxonInfluence.Rd @@ -149,6 +149,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} } diff --git a/man/TreeLength.Rd b/man/TreeLength.Rd index 4f76a73fb..2717c4b2d 100644 --- a/man/TreeLength.Rd +++ b/man/TreeLength.Rd @@ -203,6 +203,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}}, \code{\link[=WideSample]{WideSample()}} diff --git a/man/TreeScore.Rd b/man/TreeScore.Rd index 8cb083f8c..02298d1fc 100644 --- a/man/TreeScore.Rd +++ b/man/TreeScore.Rd @@ -66,6 +66,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=WideSample]{WideSample()}} } diff --git a/man/WideSample.Rd b/man/WideSample.Rd index 045ef6214..013771e51 100644 --- a/man/WideSample.Rd +++ b/man/WideSample.Rd @@ -166,6 +166,7 @@ Other tree scoring: \code{\link[=MinimumLength]{MinimumLength()}}, \code{\link[=ParsSim]{ParsSim()}}, \code{\link[=RecodeHierarchy]{RecodeHierarchy()}}, +\code{\link[=SuboptimalTrees]{SuboptimalTrees()}}, \code{\link[=TaxonInfluence]{TaxonInfluence()}}, \code{\link[=TreeScore]{TreeScore()}} } 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..f8eb74fc5 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,27 @@ 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. 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. // 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..91eaa1a4e 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(); @@ -1073,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_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_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 4d742b375..924df6b0c 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,25 @@ 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) { + // 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); + } + } + return cd; } @@ -1691,6 +1712,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 +1722,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_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/src/ts_tbr.cpp b/src/ts_tbr.cpp index a00978a7b..2ac169db5 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); } @@ -2318,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) { @@ -2337,6 +2345,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 && cd->n_splits > 0) { + 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()) { @@ -2498,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 new file mode 100644 index 000000000..c1c2c78ee --- /dev/null +++ b/tests/testthat/test-Bremer.R @@ -0,0 +1,535 @@ +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)) +} + +# 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. + +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 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) + 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(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 (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 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", { + # 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 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) + 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. 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), + "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", + 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") + 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(mpts, dat, driftCycles = 5L, annealCycles = 3L, + maxReplicates = 12L, verbosity = 0L), + "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)) +}) + +# --- 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 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") + 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"))) +}) + +# --- 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))) +}) diff --git a/tests/testthat/test-BremerParallel.R b/tests/testthat/test-BremerParallel.R new file mode 100644 index 000000000..fb4a4367c --- /dev/null +++ b/tests/testthat/test-BremerParallel.R @@ -0,0 +1,192 @@ +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) +}) + +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)) +}) diff --git a/tests/testthat/test-SuboptimalTrees.R b/tests/testthat/test-SuboptimalTrees.R new file mode 100644 index 000000000..741e33b46 --- /dev/null +++ b/tests/testthat/test-SuboptimalTrees.R @@ -0,0 +1,92 @@ +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") +}) + +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. + 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") +})