diff --git a/R/SearchControl.R b/R/SearchControl.R index 82ea4d31f..c27e3a560 100644 --- a/R/SearchControl.R +++ b/R/SearchControl.R @@ -206,6 +206,22 @@ #' uncertain regions. Constraints are cleared whenever a new best score #' is found. Only active when no user-supplied `constraint` is #' present. Default `FALSE`. +#' @param transientAutoconst Logical; use TNT-style *transient* +#' auto-constraint? Distinct from the rigid `consensusConstrain`: each +#' replicate builds a per-replicate constraint from the strict consensus +#' of the previous replicates' pool trees intersected with this +#' replicate's own Wagner start, enforces it only over the early pipeline +#' stages (initial TBR through the first sectorial pass), then *releases* +#' it before the ratchet, drift, later sectorial and final TBR stages. +#' The first replicate of each hit (replicate 0, or any replicate that +#' immediately follows a new best score) is left unconstrained. Because +#' the constraint refreshes every replicate and is never permanently +#' locked, it can guide early search without entrenching a basin. Only +#' active when neither a user `constraint` nor `consensusConstrain` is in +#' effect. Default `FALSE`. +#' @param transientAutoconstMinReps Integer; number of completed replicates +#' (with at least two best-score pool trees) required before +#' `transientAutoconst` engages. Default `2`. #' @param wagnerBias Integer; criterion for biasing taxon addition order #' during Wagner tree construction. 0 = random (default), #' 1 = Goloboff (2014) non-ambiguous-character priority, @@ -342,6 +358,8 @@ SearchControl <- function( perturbStopFactor = 2L, adaptiveLevel = FALSE, consensusConstrain = FALSE, + transientAutoconst = FALSE, + transientAutoconstMinReps = 2L, # Taxon pruning-reinsertion (T-266) pruneReinsertCycles = 0L, pruneReinsertDrop = 0.10, @@ -441,6 +459,8 @@ SearchControl <- function( perturbStopFactor = as.integer(perturbStopFactor), adaptiveLevel = as.logical(adaptiveLevel), consensusConstrain = as.logical(consensusConstrain), + transientAutoconst = as.logical(transientAutoconst), + transientAutoconstMinReps = as.integer(transientAutoconstMinReps), pruneReinsertCycles = as.integer(pruneReinsertCycles), pruneReinsertDrop = as.double(pruneReinsertDrop), pruneReinsertSelection = as.integer(pruneReinsertSelection), @@ -488,7 +508,8 @@ print.SearchControl <- function(x, ...) { "poolMaxSize", "poolSuboptimal"), "Stopping" = c("consensusStableReps", "perturbStopFactor", "adaptiveLevel", - "consensusConstrain", "adaptiveStart", + "consensusConstrain", "transientAutoconst", + "transientAutoconstMinReps", "adaptiveStart", "enumTimeFraction") ) cat("SearchControl object\n") diff --git a/R/ts-driven-compat.R b/R/ts-driven-compat.R index fea8e6127..b16f9ed34 100644 --- a/R/ts-driven-compat.R +++ b/R/ts-driven-compat.R @@ -73,6 +73,8 @@ ts_driven_search <- function( perturbStopFactor = 2L, adaptiveLevel = FALSE, consensusConstrain = FALSE, + transientAutoconst = FALSE, + transientAutoconstMinReps = 2L, nniPerturbCycles = 0L, nniPerturbFraction = 0.5, wagnerBias = 0L, @@ -134,6 +136,8 @@ ts_driven_search <- function( perturbStopFactor = as.integer(perturbStopFactor), adaptiveLevel = as.logical(adaptiveLevel), consensusConstrain = as.logical(consensusConstrain), + transientAutoconst = as.logical(transientAutoconst), + transientAutoconstMinReps = as.integer(transientAutoconstMinReps), pruneReinsertCycles = as.integer(pruneReinsertCycles), pruneReinsertDrop = as.double(pruneReinsertDrop), pruneReinsertSelection = as.integer(pruneReinsertSelection), diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 0c16caba7..22e2babfb 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -1,25 +1,3 @@ -@article{LapointeCucumel1997, - title = {The average consensus procedure: combination of weighted trees containing identical or overlapping sets of taxa}, - author = {Lapointe, Fran{\c c}ois-Joseph and Cucumel, Guy}, - year = {1997}, - journal = {Systematic Biology}, - volume = {46}, - number = {2}, - pages = {306--312}, - doi = {10.1093/sysbio/46.2.306} -} - -@article{Arias2004, - title = {Profile parsimony ({{PP}}): an analysis under implied weights ({{IW}})}, - author = {Arias, J. Salvador and {Miranda-Esquivel}, Daniel Rafael}, - year = {2004}, - volume = {20}, - pages = {56--63}, - doi = {10.1111/j.1096-0031.2003.00001.x}, - journal = {Cladistics}, - number = {1} -} - @article{Aberer2013, title = {Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice}, author = {Aberer, Andre J. and Krompass, Desmithconsnis and Stamatakis, Alexandros}, @@ -42,6 +20,17 @@ @article{Archie1989 doi = {10.2307/2992286} } +@article{Arias2004, + title = {Profile parsimony ({{PP}}): an analysis under implied weights ({{IW}})}, + author = {Arias, J. Salvador and {Miranda-Esquivel}, Daniel Rafael}, + year = {2004}, + volume = {20}, + pages = {56--63}, + doi = {10.1111/j.1096-0031.2003.00001.x}, + journal = {Cladistics}, + number = {1} +} + @article{Bien2011, title = {Hierarchical Clustering with Prototypes via Minimax Linkage}, author = {Bien, Jacob and Tibshirani, Robert}, @@ -93,6 +82,27 @@ @article{Carter1990 doi = {10.1137/0403005} } +@article{Chao1984, + author = {Chao, Anne}, + title = {Nonparametric estimation of the number of classes in a population}, + journal = {Scandinavian Journal of Statistics}, + year = {1984}, + volume = {11}, + number = {4}, + pages = {265--270} +} + +@article{Chao2012, + author = {Chao, Anne and Jost, Lou}, + title = {Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size}, + journal = {Ecology}, + year = {2012}, + volume = {93}, + number = {12}, + pages = {2533--2547}, + doi = {10.1890/11-1952.1} +} + @article{Congreve2016, author = {Congreve, Curtis R. and Lamsdell, James C.}, doi = {10.1111/pala.12236}, @@ -357,6 +367,17 @@ @article{Gonzalez1985 doi = {10.1016/0304-3975(85)90224-5} } +@article{Good1953, + author = {Good, I. J.}, + title = {The population frequencies of species and the estimation of population parameters}, + journal = {Biometrika}, + year = {1953}, + volume = {40}, + number = {3/4}, + pages = {237--264}, + doi = {10.2307/2333344} +} + @article{Hartigan1979, title = {Algorithm {{AS}} 136: A {{{\emph{K}}}}-Means Clustering Algorithm}, shorttitle = {Algorithm {{AS}} 136}, @@ -400,6 +421,17 @@ @article{Kumar2012 number = {2} } +@article{LapointeCucumel1997, + title = {The average consensus procedure: combination of weighted trees containing identical or overlapping sets of taxa}, + author = {Lapointe, Fran{\c c}ois-Joseph and Cucumel, Guy}, + year = {1997}, + journal = {Systematic Biology}, + volume = {46}, + number = {2}, + pages = {306--312}, + doi = {10.1093/sysbio/46.2.306} +} + @article{Maddison1991, title = {Null models for the number of evolutionary steps in a character on a phylogenetic tree}, author = {Maddison, Wayne P. and Slatkin, M}, @@ -474,6 +506,17 @@ @article{Muller2005 number = {1} } +@article{Nguyen2015, + title = {{IQ-TREE}: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies}, + author = {Nguyen, Lam-Tung and Schmidt, Heiko A. and von Haeseler, Arndt and Minh, Bui Quang}, + year = {2015}, + volume = {32}, + pages = {268--274}, + doi = {10.1093/molbev/msu300}, + journal = {Molecular Biology and Evolution}, + number = {1} +} + @article{Nixon1999, author = {Nixon, Kevin C.}, doi = {10.1111/j.1096-0031.1999.tb00277.x}, @@ -705,6 +748,17 @@ @article{Sun2018 doi = {10.1098/rspb.2018.1780} } +@article{Thompson1933, + title = {On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples}, + author = {Thompson, William R.}, + year = {1933}, + volume = {25}, + pages = {285--294}, + doi = {10.2307/2332286}, + journal = {Biometrika}, + number = {3/4} +} + @article{Vinther2008, author = {Vinther, Jakob and {Van Roy}, Peter and Briggs, Derek E. G.}, doi = {10.1038/nature06474}, @@ -753,57 +807,3 @@ @article{Wilkinson2017 number = {4} } -@article{Nguyen2015, - title = {{IQ-TREE}: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies}, - author = {Nguyen, Lam-Tung and Schmidt, Heiko A. and von Haeseler, Arndt and Minh, Bui Quang}, - year = {2015}, - volume = {32}, - pages = {268--274}, - doi = {10.1093/molbev/msu300}, - journal = {Molecular Biology and Evolution}, - number = {1} -} - -@article{Thompson1933, - title = {On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples}, - author = {Thompson, William R.}, - year = {1933}, - volume = {25}, - pages = {285--294}, - doi = {10.2307/2332286}, - journal = {Biometrika}, - number = {3/4} -} - -@article{Chao1984, - author = {Chao, Anne}, - title = {Nonparametric estimation of the number of classes in a population}, - journal = {Scandinavian Journal of Statistics}, - year = {1984}, - volume = {11}, - number = {4}, - pages = {265--270} -} - -@article{Good1953, - author = {Good, I. J.}, - title = {The population frequencies of species and the estimation of population parameters}, - journal = {Biometrika}, - year = {1953}, - volume = {40}, - number = {3/4}, - pages = {237--264}, - doi = {10.2307/2333344} -} - -@article{Chao2012, - author = {Chao, Anne and Jost, Lou}, - title = {Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size}, - journal = {Ecology}, - year = {2012}, - volume = {93}, - number = {12}, - pages = {2533--2547}, - doi = {10.1890/11-1952.1} -} - diff --git a/man/SearchControl.Rd b/man/SearchControl.Rd index 3c49cda79..b26b5e89d 100644 --- a/man/SearchControl.Rd +++ b/man/SearchControl.Rd @@ -56,6 +56,8 @@ SearchControl( perturbStopFactor = 2L, adaptiveLevel = FALSE, consensusConstrain = FALSE, + transientAutoconst = FALSE, + transientAutoconstMinReps = 2L, pruneReinsertCycles = 0L, pruneReinsertDrop = 0.1, pruneReinsertSelection = 0L, @@ -299,6 +301,24 @@ uncertain regions. Constraints are cleared whenever a new best score is found. Only active when no user-supplied \code{constraint} is present. Default \code{FALSE}.} +\item{transientAutoconst}{Logical; use TNT-style \emph{transient} +auto-constraint? Distinct from the rigid \code{consensusConstrain}: each +replicate builds a per-replicate constraint from the strict consensus +of the previous replicates' pool trees intersected with this +replicate's own Wagner start, enforces it only over the early pipeline +stages (initial TBR through the first sectorial pass), then \emph{releases} +it before the ratchet, drift, later sectorial and final TBR stages. +The first replicate of each hit (replicate 0, or any replicate that +immediately follows a new best score) is left unconstrained. Because +the constraint refreshes every replicate and is never permanently +locked, it can guide early search without entrenching a basin. Only +active when neither a user \code{constraint} nor \code{consensusConstrain} is in +effect. Default \code{FALSE}.} + +\item{transientAutoconstMinReps}{Integer; number of completed replicates +(with at least two best-score pool trees) required before +\code{transientAutoconst} engages. Default \code{2}.} + \item{pruneReinsertCycles}{Integer; number of taxon pruning-reinsertion perturbation cycles per replicate. Each cycle drops a fraction of leaves, runs TBR on the reduced tree to let the backbone restructure, then diff --git a/src/ts_driven.cpp b/src/ts_driven.cpp index bc7272073..9c33ed932 100644 --- a/src/ts_driven.cpp +++ b/src/ts_driven.cpp @@ -24,6 +24,7 @@ #include #include #include +#include namespace ts { @@ -60,11 +61,21 @@ ReplicateResult run_single_replicate( TreeState* starting_tree, const SplitFrequencyTable* split_freq, StartStrategy strategy, - const TreePool* pool) + const TreePool* pool, + bool transient_build) { ReplicateResult result; result.interrupted = false; + // TNT-style transient autoconstraint state (see DrivenParams docs). + // `cd` is a by-value pointer parameter, so we reassign it: nullptr for the + // (unconstrained) Wagner start, then a per-rep constraint over the early + // stages, then nullptr again before ratchet. transient_cd must outlive + // every stage that reads `cd`, so it lives at function scope. + ConstraintData transient_cd; + bool transient_active = false; + int transient_n_mapped = 0; + bool tree_large_enough_for_sectors = ds.n_tips >= 2 * params.sector_min_size; @@ -158,6 +169,83 @@ ReplicateResult run_single_replicate( } } + // 1b. TNT-style transient autoconstraint (only when requested by the + // driver for this replicate). Build a per-rep positive constraint = + // strict consensus of the PREVIOUS replicates' best-score pool trees, + // INTERSECTED with this replicate's own (unconstrained) Wagner start — + // i.e. consensus({previous trees} ∪ {this Wagner}). Because every kept + // split is present in the Wagner start, the start satisfies the + // constraint by construction (no repair needed). The constraint is + // enforced only over the early stages (initial TBR + first sectorial + // pass) and released before ratchet (see the release below). + if (transient_build && pool && pool->size() >= 2 && !starting_tree) { + int n_unan = 0, wps = 0; + std::vector cons_bits = pool->extract_consensus_splits(n_unan, wps); + int n_kept = 0; + if (n_unan > 0 && wps > 0) { + // Hash the Wagner start's canonical splits (same canonicalisation as + // the pool's stored splits, so bitset hashes are directly comparable). + SplitSet wag = compute_splits(result.tree); + std::unordered_set wag_hashes; + wag_hashes.reserve(static_cast(wag.n_splits) * 2 + 1); + for (int i = 0; i < wag.n_splits; ++i) { + wag_hashes.insert(hash_single_split(wag.split(i), wag.words_per_split)); + } + // Keep only consensus splits that the Wagner start also displays. + std::vector kept; + kept.reserve(cons_bits.size()); + for (int s = 0; s < n_unan; ++s) { + const uint64_t* sb = cons_bits.data() + static_cast(s) * wps; + if (wag_hashes.count(hash_single_split(sb, wps))) { + kept.insert(kept.end(), sb, sb + wps); + ++n_kept; + } + } + int n_kept_mapped = 0; + if (n_kept > 0) { + transient_cd = + build_constraint_from_bitsets(kept.data(), n_kept, wps, ds.n_tips); + // Map constraint nodes on the start. Every kept split is displayed + // by the Wagner start, but map_constraint_nodes matches a ROOTED node + // subtree EXACTLY, so a split whose tip-0-free clade is not a rooted + // clade of this tree (tip 0 sits inside the clade) will not map + // (constraint_node == -1) and would then spuriously reject moves. + // Empirically ~20-30% of kept splits fail to map on a 62-tip tree. + // Drop the unmapped splits and rebuild so the constraint is + // valid-by-construction (every enforced split maps to a real node), + // rather than enforcing a garbled over-restrictive constraint. + // (Mid-search, TBR re-maps after each accepted move — ts_tbr.cpp — so + // a split can still flip to -1 as tip 0 moves; that reverts to the + // conservative reject-its-violations semantics shared with the rigid + // and user-constraint paths, and only touches the early stages since + // the constraint is released before ratchet.) + update_constraint(result.tree, transient_cd); + std::vector mapped_bits; + mapped_bits.reserve(kept.size()); + for (int s = 0; s < transient_cd.n_splits; ++s) { + if (transient_cd.constraint_node[s] >= 0) { + const uint64_t* sb = + &transient_cd.split_tips[static_cast(s) * wps]; + mapped_bits.insert(mapped_bits.end(), sb, sb + wps); + ++n_kept_mapped; + } + } + if (n_kept_mapped > 0) { + transient_cd = build_constraint_from_bitsets( + mapped_bits.data(), n_kept_mapped, wps, ds.n_tips); + update_constraint(result.tree, transient_cd); // all map now + transient_n_mapped = n_kept_mapped; + cd = &transient_cd; + transient_active = true; + } + } + } + if (verbosity >= 2) { + Rprintf(" Transient autoconstraint: consensus=%d, kept(∩wagner)=%d, " + "enforced(mapped)=%d\n", n_unan, n_kept, transient_n_mapped); + } + } + // 2. Hill-climbing to local optimum. // When NNI is active and unconstrained, each Wagner start was already // NNI-optimized, and SPR is skipped (NNI→TBR outperforms NNI→SPR→TBR). @@ -316,6 +404,21 @@ ReplicateResult run_single_replicate( return result; } + // Transient autoconstraint RELEASE (TNT autoconst). The constraint has + // now guided the early stages (initial TBR + the first sectorial pass); + // release it before ratchet/drift/later-sectorial/final-TBR so late + // search is free to cross into the (unconstrained) global optimum. This + // is the whole point of the *transient* variant vs the rigid one: it + // never permanently locks a split. Fires once (outer == 0 boundary); + // transient_active is false thereafter. + if (transient_active) { + cd = nullptr; + transient_active = false; + if (verbosity >= 2) { + Rprintf(" Transient autoconstraint released (pre-ratchet)\n"); + } + } + // 4. Ratchet perturbation to escape local optima. // Skipped entirely when ratchet is disabled (ratchet_per == 0), mirroring // the drift / nni-perturb guards below. This matters because ratchet_search() @@ -741,6 +844,23 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, ConstraintData auto_cd; // built from pool consensus; reused across reps double auto_cd_best_score = 1e18; // score when auto_cd was last built + // TNT-style TRANSIENT autoconstraint (per-rep, released before ratchet; + // see run_single_replicate). Mutually exclusive with the rigid path and + // with a user constraint. Env kill-switch TS_NO_TRANSIENT_AUTOCONST + // (read once here, not per replicate: std::getenv is comparatively slow + // on some platforms, so it is hoisted out of the replicate loop below). + bool use_transient_autoconst = + params.transient_autoconst && !use_auto_constraint && (!cd || !cd->active); + { + const char* kill = std::getenv("TS_NO_TRANSIENT_AUTOCONST"); + if (kill && kill[0] != '\0' && kill[0] != '0') use_transient_autoconst = false; + } + // First-rep-of-hit tracking: a replicate is left UNCONSTRAINED when the + // immediately-preceding replicate improved the best score (a "new hit"), + // keeping hits independent (TNT: "hits are totally independent"). + // Initialised true so rep 0 is also unconstrained. + bool prev_rep_improved = true; + for (int rep = 0; rep < params.max_replicates; ++rep) { int rep1 = rep + 1; @@ -904,6 +1024,15 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, } } + // Transient autoconstraint: request a per-rep constraint build inside + // run_single_replicate for this replicate. Skip the first rep of each + // hit (prev_rep_improved) so hits stay independent; need >= 2 best-score + // pool trees for a consensus. rep_cd stays the (null) user constraint — + // the transient constraint is built AFTER the Wagner start, internally. + bool transient_now = use_transient_autoconst && !prev_rep_improved && + result.replicates_completed >= params.transient_autoconst_min_reps && + pool.size() >= 2 && !start_ptr; + // Select starting-tree strategy for this replicate. StartStrategy rep_strategy = StartStrategy::WAGNER_RANDOM; if (start_ptr) { @@ -922,7 +1051,7 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, // Run the single-replicate pipeline ReplicateResult rep_result = run_single_replicate( ds, rep_params, rep_cd, check_timeout, params.verbosity, start_ptr, - sft_ptr, rep_strategy, &pool); + sft_ptr, rep_strategy, &pool, transient_now); result.timings += rep_result.timings; @@ -950,6 +1079,9 @@ DrivenResult driven_search(TreePool& pool, DataSet& ds, } else { ++unsuccessful_reps; } + // Feed the transient-autoconstraint first-rep-of-hit rule: the NEXT rep + // is left unconstrained iff THIS rep just established a new best. + prev_rep_improved = score_improved; // Update strategy bandit (T-190) if (params.adaptive_start) { diff --git a/src/ts_driven.h b/src/ts_driven.h index 4febe1cab..ac01bfd20 100644 --- a/src/ts_driven.h +++ b/src/ts_driven.h @@ -263,6 +263,22 @@ struct DrivenParams { bool consensus_constrain = false; int consensus_constrain_min_reps = 5; // minimum replicates before engaging + // TNT-style TRANSIENT autoconstraint (distinct from the RIGID + // consensus_constrain above). Each replicate builds a per-rep constraint + // from the strict consensus of the PREVIOUS replicates' pool trees, + // intersected with THIS rep's own (unconstrained) Wagner start — i.e. + // consensus({previous trees} ∪ {this Wagner}) — then enforces it only over + // the EARLY pipeline stages (initial TBR → first sectorial pass) and + // RELEASES it before ratchet/drift/later-sectorial/final-TBR. The first + // replicate of each "hit" (rep 0, or any rep immediately after a new best) + // is left unconstrained so hits stay independent (TNT: "hits are totally + // independent"). Unlike the rigid path this refreshes EVERY rep and never + // permanently locks a split, so it can guide early search without + // entrenching a basin. Only active when no user constraint is present and + // consensus_constrain is off. Env kill-switch: TS_NO_TRANSIENT_AUTOCONST. + bool transient_autoconst = false; + int transient_autoconst_min_reps = 2; // completed reps before engaging + // Fraction of the time budget reserved for MPT enumeration (T-202). // The main search loop exits at budget × (1 - enum_time_fraction), // leaving the remainder for the plateau walk. Default 0.1 = 10%. @@ -366,7 +382,8 @@ ReplicateResult run_single_replicate( TreeState* starting_tree = nullptr, const SplitFrequencyTable* split_freq = nullptr, StartStrategy strategy = StartStrategy::WAGNER_RANDOM, - const TreePool* pool = nullptr); + const TreePool* pool = nullptr, + bool transient_build = false); // Run the full driven search. Returns search statistics. // The pool contents (all retained trees) are accessible via the pool diff --git a/src/ts_rcpp.cpp b/src/ts_rcpp.cpp index 4d742b375..21f2d2cfd 100644 --- a/src/ts_rcpp.cpp +++ b/src/ts_rcpp.cpp @@ -1571,6 +1571,9 @@ static void unpack_search_control(List ctrl, ts::DrivenParams& params) { params.perturb_stop_factor = as(ctrl["perturbStopFactor"]); params.adaptive_level = as(ctrl["adaptiveLevel"]); params.consensus_constrain = as(ctrl["consensusConstrain"]); + params.transient_autoconst = as(ctrl["transientAutoconst"]); + params.transient_autoconst_min_reps = + as(ctrl["transientAutoconstMinReps"]); params.adaptive_start = as(ctrl["adaptiveStart"]); params.enum_time_fraction = as(ctrl["enumTimeFraction"]);