diff --git a/NEWS.md b/NEWS.md index eb2e71ad6..fb7534252 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,25 @@ # To integrate into 2.0.0 notes +- New `SearchControl()` parameter `enumMaxTrees`: a retention ceiling applied to + the post-search MPT-enumeration phase alone. `0` (the default) + keeps `poolMaxSize` throughout, so behaviour is unchanged unless you set it. + From `effort` rung 5 the ladder now doubles it each notch alongside the + replicate budget and the hit target, so asking for more effort also asks for a + more complete tree set — previously a search could be given eight times the + budget and still return only the default 100 trees. + + `poolMaxSize` is deliberately **not** scaled, and the split is the point. + During the replicate loop the pool cap is not a ceiling on what is returned but + the size of the working set the search reads: fusing draws its donors from the + whole pool, conflict-guided sectorial search reads the pool's split + frequencies once per replicate, and `consensusConstrain` reads its consensus + splits. Raising it therefore changes which trees the search *visits*, so the + anytime-dominance argument that licenses raising `maxReplicates` — a higher cap + only appends later replicates and can never delay an earlier improvement — does + not transfer to it. Once the loop is over the pool is pure output, and there + the same argument does hold, which is why the ceiling is raised at that point + instead. Raising `poolMaxSize` yourself still works and still governs both + phases; `enumMaxTrees` is the side-effect-free way to keep more trees. - Profile parsimony computes exactly for more multi-state characters, where it previously approximated nearly all of them. The exact Maddison & Slatkin solver caches into fixed-capacity memo tables and bails out when one fills -- diff --git a/R/Concordance.R b/R/Concordance.R index 0f3077352..9cb70455d 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -177,20 +177,10 @@ ClusteringConcordance <- function( dataset <- dataset[keep] # Prepare data - # `tree` may carry tips absent from `dataset` (already dropped from `keep`). - # Restrict `splits` to the shared taxon set via `Subsplit()` rather than - # pruning `tree` itself: retaining extra tips computes bipartitions over a - # different taxon set than `dataset` describes (and column-indexing the - # unpruned splits matrix by `keep` cannot recover the correct, smaller set - # of splits), but `KeepTip()` would renumber nodes and break any caller - # (e.g. `ConcordanceTable()`, `PaintCharacters()`) that matches these split - # names against `tree$edge`. + # `tree` may carry tips absent from `dataset`; Subsplit restricts `splits` to the + # shared taxa without renumbering nodes. splits <- as.logical(Subsplit(as.Splits(tree), keep)) - # `Subsplit()` drops row names entirely when exactly one split survives - # restriction (its own version of the drop-to-a-vector bug this file works - # around elsewhere) -- recover each surviving split's original node number - # by matching it (or its complement) against `tree`'s own splits, likewise - # restricted to `keep`'s columns. + # Recover each surviving split's original node number if (is.null(rownames(splits))) { fullRestricted <- as.logical(as.Splits(tree))[, TipLabels(tree) %in% keep, drop = FALSE] diff --git a/R/Consistency.R b/R/Consistency.R index 639e085fb..dc59d113d 100644 --- a/R/Consistency.R +++ b/R/Consistency.R @@ -156,19 +156,9 @@ ExpectedLength <- function(dataset, tree, nRelabel = 1000, compress = FALSE) { as.integer(intToBits(x)[1:nLevels]) }, integer(nLevels))) - # Key on the unlabelled rooted shape, which is what the sampled distribution - # is a function of: leaf states are permuted uniformly, and relabelling - # composes with a uniform permutation to leave it uniform, so any two trees - # of the same shape are sampling the same distribution. Keying on the - # labelled topology instead would be sound but strictly weaker -- identical - # topologies are a subset of identical shapes, so it would miss every reuse - # this catches and none of its own. Rooting is part of the shape, as these - # characters may contain inapplicable tokens, whose lengths are not - # rooting-invariant. + # Key on the unlabelled rooted shape treeKey <- .ShapeKey(tree) - # Cache per shape, and within that per character, rather than pasting both - # into one key: that keeps the shape key out of every character's entry, and - # leaves no ambiguity about where the shape key ends and the counts begin. + # Cache per shape, and within that per character treeCache <- .CharLengthCache[[treeKey]] if (is.null(treeCache)) { treeCache <- new.env(hash = TRUE, parent = emptyenv()) @@ -218,11 +208,7 @@ ExpectedLength <- function(dataset, tree, nRelabel = 1000, compress = FALSE) { # Canonical identifier of a rooted tree's unlabelled shape, after # Aho, Hopcroft & Ullman: a leaf encodes as `01`, and an internal node wraps -# its children's codes, sorted into a fixed order, in `0`...`1`. Sorting is -# what makes the code canonical, so it is already invariant to edge order and -# to node rotation, and two rooted shapes are isomorphic exactly if their codes -# agree. Unlike `TreeTools::RootedTreeShape()`, which enumerates shapes into -# an integer and so stops at 55 leaves, this is bounded only by string length. +# its children's codes, sorted into a fixed order, in `0`...`1`. # @param tree A rooted, binary tree of class `phylo`. # @return A string identifying the shape of `tree`. #' @importFrom TreeTools NTip Postorder diff --git a/R/MaximizeParsimony.R b/R/MaximizeParsimony.R index 4c76daa60..d2f39209d 100644 --- a/R/MaximizeParsimony.R +++ b/R/MaximizeParsimony.R @@ -122,15 +122,7 @@ consContrast <- attr(constraint, "contrast") nConsStates <- ncol(consContrast) if (nConsStates < 2L) { - # One state means no taxon is coded `0`, so this is the extreme case of the - # inert character warned about below -- and the loudest one, because it is - # what `MatrixToPhyDat(c(a = 1, b = 1, c = 1))` produces: a user asking for - # a clade and getting no constraint at all. Warn here rather than returning - # silently; the group-size test below never sees these characters. - warning("Constraint constrains nothing, and is ignored: no taxon is coded ", - "`0`, so every tree separates the `1` taxa from the (empty) `0` ", - "group. Code the taxa that must fall outside the group as `0`.", - call. = FALSE) + warning("Ignoring empty constraint", call. = FALSE) return(list()) } @@ -150,8 +142,7 @@ # For each constraint character, record the tips unambiguously in the "1" # group (derived state present, ancestral absent) and, separately, those in # the "0" group (ancestral present, derived absent). Tips ambiguous for the - # character ("?", or unconstrained taxa) belong to neither group and are free - # to plot on either side of the split. + # character ("?", or unconstrained taxa) may plot anywhere. consSplits <- matrix(0L, nrow = ncol(consMat), ncol = length(constraint)) consZero <- matrix(0L, nrow = ncol(consMat), ncol = length(constraint)) for (ch in seq_len(ncol(consMat))) { @@ -167,28 +158,13 @@ } } - # Every tree separates a group of fewer than two taxa from anything: the edge - # above a lone tip already does it, and an empty group needs no edge at all. - # Such a character constrains nothing under the documented contract, so - # enforcing its group as a clade would restrict the search for a guarantee it - # already has -- which is the over-strict reading agent-issues/TreeSearch#54 - # is about. The test is symmetric in the two groups because they are - # interchangeable: which one a user calls "1" is arbitrary, and - # build_constraint() swaps them freely to canonicalise. - # - # Warn rather than drop silently: a character coding only "1" and "?" almost - # certainly means "group these taxa", which is not what it says. + # Ignore trivial constraints nOne <- rowSums(consSplits) nZero <- rowSums(consZero) inert <- nOne < 2 | nZero < 2 if (any(inert)) { - warning("Constraint character", if (sum(inert) > 1) "s" else "", " ", - paste(which(inert), collapse = ", "), - if (sum(inert) > 1) " constrain" else " constrains", - " nothing, and", if (sum(inert) > 1) " are" else " is", - " ignored: every tree separates a group of fewer than two taxa ", - "from the rest. Taxa coded `?` join neither group; code those ", - "that must fall outside the group as `0`.", call. = FALSE) + warning("Ignoring trivial constraint character", if (sum(inert) > 1) "s" else "", " ", + paste(which(inert), collapse = ", "), call. = FALSE) } keep <- !inert consSplits <- consSplits[keep, , drop = FALSE] @@ -196,14 +172,6 @@ if (nrow(consSplits) == 0L) return(list()) # Every returned tree must display all constraint splits simultaneously. - # Two splits are jointly displayable iff they are compatible in the - # four-gamete sense: treating each as a bipartition of the tips it - # constrains (its "1" group vs its "0" group, ambiguous tips excluded), the - # pair is compatible iff at least one of the four group intersections is - # empty. A laminar (nested-or-disjoint) test alone is too strict: it rejects - # the case where the two "0" groups are disjoint -- i.e. the splits' "1" - # sides jointly cover the constrained tips -- which is perfectly displayable, - # e.g. ab | cef and abcd | ef coexist on ((a,b),(d,(c,(e,f)))). nSplits <- nrow(consSplits) if (nSplits > 1L) { for (i in seq_len(nSplits - 1L)) { @@ -216,7 +184,7 @@ !any(aZero & bOne) || !any(aZero & bZero) if (!compatible) { stop("Constraint is impossible to satisfy: splits ", i, " and ", j, - " are incompatible (all four taxon groupings co-occur)") + " are incompatible") } } } @@ -602,7 +570,32 @@ # through .IwRatchetDepth()'s targetHits/defaultHits escalation (capped at # .iwRatchetMaxCycles), which is a genuine reach lever the equal-weights # measurement above cannot see. - hitMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L)) + hitMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L)), + # `enumMaxTrees` multiplier: the size of the returned MPT set, relative to + # `poolMaxSize`. Doubling in step with the other two knobs, from rung 5, so + # that one notch keeps meaning "roughly twice the work" on this axis too. + # + # This scales the ENUMERATION ceiling only, never `poolMaxSize` itself, and + # the distinction is the whole point. During the replicate loop the pool cap + # is the size of the working set the search reads -- fuse donors are the + # entire pool (uncapped, and taken under the pool mutex on the parallel + # path), conflict-guided sector selection reads the pool's split frequencies + # once per replicate, and `consensusConstrain` reads its consensus splits -- + # so scaling it would change which trees the search VISITS. The + # anytime-dominance argument that licenses raising `maxReplicates` above + # therefore does NOT transfer to `poolMaxSize`: a bigger pool can delay + # every later improvement rather than merely appending to the result. + # After the loop, the pool is pure output and a bigger ceiling can only + # append equal-score topologies, so the same argument DOES hold there. + # + # It is bounded in practice without needing a cap: enumeration shares the + # `maxSeconds * enumTimeFraction` reserve, and its loop exits as soon as the + # pool fills, so an over-generous ceiling costs enumeration time, never a + # worse tree. As with `hitMultiplier` the doubling shape is an OPERATING + # POINT rather than a measurement -- what is measured is that the July 2026 + # 182-tip runs returned exactly `poolMaxSize` trees in all four analyses, + # i.e. the ceiling bound the answer rather than the MPT count doing so. + enumMultiplier = if (rung <= 4L) 1L else as.integer(2^(rung - 4L)) ) } @@ -677,9 +670,13 @@ #' most-parsimonious tree (\acronym{MPT}) is recovered. #' The size of the returned set is bounded by, in order: #' \enumerate{ -#' \item **`poolMaxSize`** (default `100`) — a hard ceiling on the number of -#' trees retained. Raise it (via [`SearchControl()`]) to keep more MPTs; -#' with the default you will never see more than 100. +#' \item **`enumMaxTrees`**, falling back to **`poolMaxSize`** (default `100`) +#' when `enumMaxTrees` is `0` — a hard ceiling on the number of trees +#' retained; with the default you will never see more than 100. Prefer +#' raising `enumMaxTrees` (via [`SearchControl()`]): it applies only once the +#' search is over, so it cannot alter which trees are visited, whereas +#' `poolMaxSize` also sizes the working set that fusing and sectorial search +#' read. From `effort` rung 5 the ladder raises `enumMaxTrees` for you. #' \item **MPT-enumeration time.** After the main search, a TBR plateau walk #' enumerates equal-score neighbours of each pool tree, within a time #' reserve of `maxSeconds * enumTimeFraction`. If this phase times out it @@ -851,9 +848,13 @@ #' TBR-disconnected islands that random restarts alone miss.} #' \item{4, `large`}{`thorough`'s provisioning with `maxReplicates` raised #' to 500, to suit the higher per-replicate cost of big trees.} -#' \item{5 and up}{`thorough`'s provisioning, with both the replicate budget -#' and the hit target doubling each notch (1000, 2000, 4000 ... -#' replicates), so that one notch always means roughly twice the work. +#' \item{5 and up}{`thorough`'s provisioning, with the replicate budget, the +#' hit target and the \acronym{MPT}-enumeration ceiling (`enumMaxTrees`) +#' all doubling each notch (1000, 2000, 4000 ... replicates), so that one +#' notch always means roughly twice the work. `poolMaxSize` is +#' deliberately *not* scaled: it sizes the working set that fusing and +#' sectorial search read during the run, so raising it would change which +#' trees are visited rather than only how many are returned. #' There is no policy ceiling: extra replicates cannot cost reach, only #' wall, which is what you asked to spend. The ladder stops only at rung #' 26, where the replicate budget outgrows R's integer type.} @@ -1283,6 +1284,18 @@ MaximizeParsimony <- function( targetHits <- as.integer(targetHits * spec[["hitMultiplier"]]) } + # Rung-scaled MPT-enumeration ceiling (rung 5 and up). Keyed off the + # POST-merge `poolMaxSize`, so a user who raised the pool gets a + # proportionally larger returned set rather than having their value + # ignored. Skipped when the user named `enumMaxTrees` themselves. + # `poolMaxSize` is deliberately not touched -- see .RungSpec(). + if (!("enumMaxTrees" %in% union(names(controlDots), + attr(control, "explicit"))) && + spec[["enumMultiplier"]] > 1L) { + control[["enumMaxTrees"]] <- + as.integer(control[["poolMaxSize"]] * spec[["enumMultiplier"]]) + } + # Implied-weights ratchet depth. Under implied weights the optimum often # sits in a small basin at fine score resolution, separated from an # easy-to-find near-optimum by a fraction of a step; character reweighting diff --git a/R/ParsSim.R b/R/ParsSim.R index 51f1bb3b6..11c8acd38 100644 --- a/R/ParsSim.R +++ b/R/ParsSim.R @@ -132,11 +132,7 @@ ParsSim <- function(tree, # --- Determine state counts per character ---------------------------------- n_states_vec <- rep(seq_along(nChar) + 1L, times = nChar) if (any(n_states_vec > 31L)) { - stop("ParsSim() supports at most 31 states per character (state codes ", - "0:30): the internal Fitch bit-set representation packs states ", - "into a 32-bit integer via bitwShiftL(), which silently overflows ", - "to NA beyond that. Requested up to ", max(n_states_vec), - " states via `nChar`.") + stop("ParsSim() supports at most 31 states per character.") } # --- Validate and expand rootState ------------------------------------------ diff --git a/R/SearchControl.R b/R/SearchControl.R index a528ac9c4..c3607baac 100644 --- a/R/SearchControl.R +++ b/R/SearchControl.R @@ -123,8 +123,22 @@ #' within each replicate, after TBR polish. This approximates TNT's #' within-replicate fusing pattern. Default: `FALSE`. #' @param poolMaxSize Integer; maximum trees retained in the pool. +#' This governs the pool throughout the search, so it is not only a ceiling on +#' the trees returned: fuse draws its donors from the pool, and conflict-guided +#' sector selection and `consensusConstrain` both read it. Raising it +#' therefore changes the trajectory as well as the output; to keep more +#' most-parsimonious trees without that side effect, use `enumMaxTrees`. #' @param poolSuboptimal Numeric; retain trees that are this many steps #' worse than the best tree. 0 (default) keeps only optimal trees. +#' @param enumMaxTrees Integer; retention ceiling applied to the +#' \acronym{MPT}-enumeration phase alone, which runs after the last replicate. +#' `0` (default) keeps `poolMaxSize` throughout, reproducing the behaviour +#' before this argument existed. Because enumeration happens once the search +#' is over, a larger ceiling here only appends further equal-score topologies: +#' it cannot change which trees the search visits. This is the knob +#' [`MaximizeParsimony()`]'s `effort` scales; `poolMaxSize` is deliberately +#' left alone. Values below `poolMaxSize` are ignored (the ceiling is only +#' ever raised). #' @param consensusStableReps Integer; stop when the strict consensus of #' best-score pool trees has been unchanged for this many consecutive #' replicates. @@ -377,7 +391,8 @@ SearchControl <- function( # sampling from {Wagner-random, Wagner-Goloboff, Wagner-entropy, # random-tree, pool-ratchet, pool-NNI-perturb}. Overrides wagnerBias. adaptiveStart = FALSE, - enumTimeFraction = 0.1 + enumTimeFraction = 0.1, + enumMaxTrees = 0L ) { # Record which fields the caller set explicitly (by name or position; # `match.call()` normalises positional args to their names). This lets @@ -398,6 +413,14 @@ SearchControl <- function( stop("`", .p, "` must be a single positive integer") } } + # `enumMaxTrees` takes 0 ("follow poolMaxSize") but never a negative: the + # kernel only ever RAISES the ceiling, so a negative would be silently inert + # rather than reported, hiding a sign typo. + .emt <- as.integer(enumMaxTrees) + if (length(.emt) != 1L || is.na(.emt) || .emt < 0L) { + stop("`enumMaxTrees` must be a single non-negative integer ", + "(0 follows `poolMaxSize`)") + } # `stopPatience` is a replicate count, so a negative value is meaningless; the # kernel treats anything <= 0 as "off", which would silently ignore a typo # such as -20 rather than honouring the obvious intent. @@ -498,7 +521,8 @@ SearchControl <- function( annealTEnd = as.double(annealTEnd), annealMovesPerPhase = as.integer(annealMovesPerPhase), adaptiveStart = as.logical(adaptiveStart), - enumTimeFraction = as.double(enumTimeFraction) + enumTimeFraction = as.double(enumTimeFraction), + enumMaxTrees = .emt ), class = "SearchControl", explicit = .explicit @@ -530,7 +554,7 @@ print.SearchControl <- function(x, ...) { "sectorCombStarts", "sectorFuseRounds", "postRatchetSectorial"), "Fuse/Pool" = c("fuseInterval", "fuseAcceptEqual", "intraFuse", - "poolMaxSize", "poolSuboptimal"), + "poolMaxSize", "poolSuboptimal", "enumMaxTrees"), "Stopping" = c("consensusStableReps", "perturbStopFactor", "stopPatience", "adaptiveLevel", "consensusConstrain", "adaptiveStart", diff --git a/R/ts-driven-compat.R b/R/ts-driven-compat.R index 8b6b141b0..5e9aad82b 100644 --- a/R/ts-driven-compat.R +++ b/R/ts-driven-compat.R @@ -42,6 +42,7 @@ ts_driven_search <- function( fuseAcceptEqual = FALSE, poolMaxSize = 100L, poolSuboptimal = 0.0, + enumMaxTrees = 0L, maxSeconds = 0.0, verbosity = 0L, min_steps = integer(0), @@ -140,7 +141,8 @@ ts_driven_search <- function( pruneReinsertDrop = as.double(pruneReinsertDrop), pruneReinsertSelection = as.integer(pruneReinsertSelection), adaptiveStart = as.logical(adaptiveStart), - enumTimeFraction = as.double(enumTimeFraction) + enumTimeFraction = as.double(enumTimeFraction), + enumMaxTrees = as.integer(enumMaxTrees) ) # Anneal config: fold into SearchControl if provided diff --git a/man/MaximizeParsimony.Rd b/man/MaximizeParsimony.Rd index 11f33cfe2..d8a8a8fde 100644 --- a/man/MaximizeParsimony.Rd +++ b/man/MaximizeParsimony.Rd @@ -181,9 +181,13 @@ cycle loop. The drift cycles also recover equal-score trees on TBR-disconnected islands that random restarts alone miss.} \item{4, \code{large}}{\code{thorough}'s provisioning with \code{maxReplicates} raised to 500, to suit the higher per-replicate cost of big trees.} -\item{5 and up}{\code{thorough}'s provisioning, with both the replicate budget -and the hit target doubling each notch (1000, 2000, 4000 ... -replicates), so that one notch always means roughly twice the work. +\item{5 and up}{\code{thorough}'s provisioning, with the replicate budget, the +hit target and the \acronym{MPT}-enumeration ceiling (\code{enumMaxTrees}) +all doubling each notch (1000, 2000, 4000 ... replicates), so that one +notch always means roughly twice the work. \code{poolMaxSize} is +deliberately \emph{not} scaled: it sizes the working set that fusing and +sectorial search read during the run, so raising it would change which +trees are visited rather than only how many are returned. There is no policy ceiling: extra replicates cannot cost reach, only wall, which is what you asked to spend. The ladder stops only at rung 26, where the replicate budget outgrows R's integer type.} @@ -439,9 +443,13 @@ topologies held in its tree pool; it does not guarantee that every most-parsimonious tree (\acronym{MPT}) is recovered. The size of the returned set is bounded by, in order: \enumerate{ -\item \strong{\code{poolMaxSize}} (default \code{100}) — a hard ceiling on the number of -trees retained. Raise it (via \code{\link[=SearchControl]{SearchControl()}}) to keep more MPTs; -with the default you will never see more than 100. +\item \strong{\code{enumMaxTrees}}, falling back to \strong{\code{poolMaxSize}} (default \code{100}) +when \code{enumMaxTrees} is \code{0} — a hard ceiling on the number of trees +retained; with the default you will never see more than 100. Prefer +raising \code{enumMaxTrees} (via \code{\link[=SearchControl]{SearchControl()}}): it applies only once the +search is over, so it cannot alter which trees are visited, whereas +\code{poolMaxSize} also sizes the working set that fusing and sectorial search +read. From \code{effort} rung 5 the ladder raises \code{enumMaxTrees} for you. \item \strong{MPT-enumeration time.} After the main search, a TBR plateau walk enumerates equal-score neighbours of each pool tree, within a time reserve of \code{maxSeconds * enumTimeFraction}. If this phase times out it diff --git a/man/SearchControl.Rd b/man/SearchControl.Rd index d46f7c037..9666063ba 100644 --- a/man/SearchControl.Rd +++ b/man/SearchControl.Rd @@ -69,7 +69,8 @@ SearchControl( annealTEnd = 0, annealMovesPerPhase = 0L, adaptiveStart = FALSE, - enumTimeFraction = 0.1 + enumTimeFraction = 0.1, + enumMaxTrees = 0L ) } \arguments{ @@ -258,7 +259,12 @@ TNT's interleaved sectorial pattern. Default: \code{FALSE}.} within each replicate, after TBR polish. This approximates TNT's within-replicate fusing pattern. Default: \code{FALSE}.} -\item{poolMaxSize}{Integer; maximum trees retained in the pool.} +\item{poolMaxSize}{Integer; maximum trees retained in the pool. +This governs the pool throughout the search, so it is not only a ceiling on +the trees returned: fuse draws its donors from the pool, and conflict-guided +sector selection and \code{consensusConstrain} both read it. Raising it +therefore changes the trajectory as well as the output; to keep more +most-parsimonious trees without that side effect, use \code{enumMaxTrees}.} \item{poolSuboptimal}{Numeric; retain trees that are this many steps worse than the best tree. 0 (default) keeps only optimal trees.} @@ -393,6 +399,16 @@ equal-score topologies). The main search loop exits at \code{maxSeconds * (1 - enumTimeFraction)}. Set to 0 to disable the reserve (pre-v1.6 behaviour: enumeration skipped if the main loop times out). Default: \code{0.1} (10\%).} + +\item{enumMaxTrees}{Integer; retention ceiling applied to the +\acronym{MPT}-enumeration phase alone, which runs after the last replicate. +\code{0} (default) keeps \code{poolMaxSize} throughout, reproducing the behaviour +before this argument existed. Because enumeration happens once the search +is over, a larger ceiling here only appends further equal-score topologies: +it cannot change which trees the search visits. This is the knob +\code{\link[=MaximizeParsimony]{MaximizeParsimony()}}'s \code{effort} scales; \code{poolMaxSize} is deliberately +left alone. Values below \code{poolMaxSize} are ignored (the ceiling is only +ever raised).} } \value{ A named list of class \code{"SearchControl"}. diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index f5b33c99a..3a9708eb1 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -1436,6 +1436,10 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, // one tree, so different TBR-connected islands are only discovered if // different replicates landed on them. We enumerate from each seed // tree to explore its island, stopping when the pool is full. + // + // The retention ceiling is raised HERE, not before the loop: past this + // point the pool is pure output, so a larger cap only appends topologies. + pool.raise_max_size(params.enum_pool_max_size); if (pool.size() > 0 && pool.size() < pool.max_size) { TBRParams tp; tp.accept_equal = true; diff --git a/src/ts_driven.h b/src/ts_driven.h index 3a54c0e80..c51a511b6 100644 --- a/src/ts_driven.h +++ b/src/ts_driven.h @@ -142,6 +142,17 @@ struct DrivenParams { int pool_max_size = 100; double pool_suboptimal = 0.0; // 0 = keep only optimal + // Retention ceiling for the MPT-enumeration phase alone; 0 means "no separate + // ceiling", i.e. keep `pool_max_size` throughout, which is the default and is + // byte-identical to the behaviour before this field existed. + // + // Split from `pool_max_size` because only this half is safe to scale with + // search effort: during the replicate loop the cap is the size of the working + // set that fuse, sector selection and consensusConstrain all read, so raising + // it changes the search trajectory; after the loop it is purely how many + // equal-score topologies get returned. See TreePool::raise_max_size(). + int enum_pool_max_size = 0; + // Timeout (seconds). 0 or negative = no timeout. double max_seconds = 0.0; diff --git a/src/ts_parallel.cpp b/src/ts_parallel.cpp index 4035d0f30..b85269104 100644 --- a/src/ts_parallel.cpp +++ b/src/ts_parallel.cpp @@ -621,6 +621,11 @@ DrivenResult parallel_driven_search( return false; }; + // Raise the retention ceiling for enumeration only. `pool_out` is already + // the OUTPUT pool -- the shared search pool keeps its own `pool_max_size`, so + // fuse donors and sector selection are untouched by this. + pool_out.raise_max_size(params.enum_pool_max_size); + if (pool_out.size() > 0 && pool_out.size() < pool_out.max_size) { TBRParams tp; tp.accept_equal = true; diff --git a/src/ts_pool.h b/src/ts_pool.h index 84327301d..8a12e7c70 100644 --- a/src/ts_pool.h +++ b/src/ts_pool.h @@ -62,6 +62,21 @@ class TreePool { // Get all entries. const std::vector& all() const { return entries_; } + // Raise the retention ceiling; never lowers it, so entries already held can + // not be orphaned above the cap. + // + // Call this ONLY at the MPT-enumeration phase, where the pool is pure OUTPUT. + // During the replicate loop `max_size` is not a ceiling on what is *returned* + // but the size of the working set the search reads: fuse donors are + // `all().size()` (uncapped, and held under the pool mutex on the parallel + // path), conflict-guided sector selection calls compute_split_frequencies() + // over best-score entries once per replicate, and consensusConstrain calls + // extract_consensus_splits(). Raising it mid-loop would therefore change what + // the search DOES, not merely how much it keeps -- which is why the effort + // ladder raises it here, after the loop has finished, rather than scaling + // `poolMaxSize` itself. + void raise_max_size(int n) { if (n > max_size) max_size = n; } + // Evict entries worse than best_score + suboptimal. void evict(); diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index d970e2cc9..9fcc11bf5 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -1655,6 +1655,10 @@ static void unpack_search_control(List ctrl, ts::DrivenParams& params) { params.intra_fuse = as(ctrl["intraFuse"]); params.pool_max_size = as(ctrl["poolMaxSize"]); params.pool_suboptimal = as(ctrl["poolSuboptimal"]); + // Absent in a control list built by an older caller: treat as 0 ("no separate + // enumeration ceiling") rather than letting as() throw on R_NilValue. + params.enum_pool_max_size = ctrl.containsElementNamed("enumMaxTrees") + ? as(ctrl["enumMaxTrees"]) : 0; // Stopping / adaptive params.consensus_stable_reps = as(ctrl["consensusStableReps"]); diff --git a/tests/testthat/test-MaximizeParsimony-features.R b/tests/testthat/test-MaximizeParsimony-features.R index a273eabca..91daaadbe 100644 --- a/tests/testthat/test-MaximizeParsimony-features.R +++ b/tests/testthat/test-MaximizeParsimony-features.R @@ -214,6 +214,59 @@ test_that("both budget knobs double, so a notch is the same size either way", { for (r in 1:4) expect_equal(RS(r)[["hitMultiplier"]], 1L) }) +test_that("the ladder scales the MPT ceiling at the same rate", { + RS <- TreeSearch:::.RungSpec + for (r in 5:9) { + expect_equal(RS(r)[["enumMultiplier"]] / RS(r - 1L)[["enumMultiplier"]], 2) + } + for (r in 1:4) expect_equal(RS(r)[["enumMultiplier"]], 1L) +}) + +test_that("no preset touches enumMaxTrees, and an explicit value survives", { + # The ladder must be the only thing that sets this field automatically, so no + # preset may carry a value of its own. + for (p in TreeSearch:::.StrategyPresets()) { + expect_equal(p[["enumMaxTrees"]], 0L) + } + # A user-supplied `enumMaxTrees` outranks the ladder, exactly as a + # user-supplied `targetHits` or `maxReplicates` does; the merge must record it + # as explicit for the rung block to skip it. + ctrl <- TreeSearch:::.ApplyStrategyPreset( + SearchControl(enumMaxTrees = 7L), + TreeSearch:::.StrategyPresets()[["thorough"]], + explicitDots = character(0) + ) + expect_equal(ctrl[["enumMaxTrees"]], 7L) + expect_true("enumMaxTrees" %in% attr(ctrl, "explicit")) +}) + +test_that("enumMaxTrees defaults to 0 and rejects a negative", { + expect_equal(SearchControl()[["enumMaxTrees"]], 0L) + expect_equal(SearchControl(enumMaxTrees = 250L)[["enumMaxTrees"]], 250L) + expect_error(SearchControl(enumMaxTrees = -1L), "non-negative") + # 0 is legal (it means "follow poolMaxSize"), unlike poolMaxSize itself. + expect_error(SearchControl(poolMaxSize = 0L), "positive integer") +}) + +test_that("a raised enumMaxTrees can return more trees than poolMaxSize", { + # End-to-end: the ceiling has to reach the C++ enumeration phase, so this + # fails if the field is dropped anywhere along SearchControl -> ctrl list -> + # DrivenParams -> TreePool::raise_max_size(). + skip_on_cran() + set.seed(90210) + small <- MaximizeParsimony(ds, maxReplicates = 3L, targetHits = 99L, + poolMaxSize = 4L, verbosity = 0L) + set.seed(90210) + big <- MaximizeParsimony(ds, maxReplicates = 3L, targetHits = 99L, + poolMaxSize = 4L, enumMaxTrees = 40L, + verbosity = 0L) + expect_lte(length(small), 4L) + expect_gt(length(big), length(small)) + # Raising the ceiling must not change the score: enumeration only ever adds + # EQUAL-score topologies. + expect_equal(attr(big, "score"), attr(small, "score")) +}) + test_that("effort = 0 reproduces the automatic choice on every size band", { # The whole point of an offset rather than an absolute level: the default # must be byte-identical to what the package chose before `effort` existed. diff --git a/tests/testthat/test-ts-constraint-free-taxa.R b/tests/testthat/test-ts-constraint-free-taxa.R index 677246a24..805eaf2e3 100644 --- a/tests/testthat/test-ts-constraint-free-taxa.R +++ b/tests/testthat/test-ts-constraint-free-taxa.R @@ -222,7 +222,7 @@ test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { expect_warning( dropped <- TreeSearch:::.PrepareConstraint( Inert("1", "1", "?", "?", "?", "?", "?", "?"), dataset), - "constrains nothing") + "trivial constraint") expect_equal(dropped, list()) # A `0` group of one. The two groups are interchangeable, so this must be # treated exactly like its mirror image below -- which the old @@ -230,11 +230,11 @@ test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { expect_warning( TreeSearch:::.PrepareConstraint( Inert("1", "1", "0", "?", "?", "?", "?", "?"), dataset), - "constrains nothing") + "trivial constraint") expect_warning( TreeSearch:::.PrepareConstraint( Inert("0", "0", "1", "?", "?", "?", "?", "?"), dataset), - "constrains nothing") + "trivial constraint") # Two and two: kept, and kept silently. expect_silent(TreeSearch:::.PrepareConstraint( Inert("1", "1", "0", "0", "?", "?", "?", "?"), dataset)) @@ -246,7 +246,7 @@ test_that(".PrepareConstraint codes free taxa as NA and drops vacuous rows", { expect_warning( TreeSearch:::.PrepareConstraint( TreeTools::MatrixToPhyDat(c(a = "1", b = "1", c = "1")), dataset), - "constrains nothing") + "empty constraint") }) test_that("the Wagner build places free taxa freely", {