diff --git a/kb/decisions/coalescent-skyline-convex-inverse-tc.md b/kb/decisions/coalescent-skyline-convex-inverse-tc.md index 9fd6a15c..cc8f7018 100644 --- a/kb/decisions/coalescent-skyline-convex-inverse-tc.md +++ b/kb/decisions/coalescent-skyline-convex-inverse-tc.md @@ -9,7 +9,7 @@ over $\log T_c$ knots. from v0's optimizer, extending the constant-$T_c$ decision to the skyline). **Previous v1**: `optimize_skyline()` used `argmin`'s `NelderMead` over a -piecewise-*linear* $\log T_c$ grid with a stiffness penalty on adjacent $\log T_c$ +piecewise-_linear_ $\log T_c$ grid with a stiffness penalty on adjacent $\log T_c$ and a boundary penalty, plus a hand-built initial simplex. **v0**: `MergerModel.optimize_skyline()` at @@ -59,10 +59,18 @@ to the tree's time span). This is deliberate: with equal-time boundaries the boundary segment between the youngest coalescence and the tips has lineages ($I_i>0$) but no mergers ($M_i=0$). With $M_i=0$ the $-M_i\ln x_i$ barrier vanishes and the segment's optimum sits at $x_i\to 0$, i.e. $T_c\to\infty$ (only -weakly restrained by smoothing). Quantile boundaries guarantee every segment — -including both boundary segments — owns mergers, so $M_i>0$ keeps every $x_i$ -bounded away from zero. On ebola/20 this changed a divergent most-recent segment -($T_c\approx 3\times10^{150}$) into a finite, smoothly varying trajectory. +weakly restrained by smoothing). Quantile boundaries let every segment own mergers +so $M_i>0$ keeps every $x_i$ bounded away from zero. On ebola/20 this changed a +divergent most-recent segment ($T_c\approx 3\times10^{150}$) into a finite, smoothly +varying trajectory. + +This holds only when the requested segment count does not exceed the number of +distinct merger times. When it does, the clamped quantile construction emits +duplicate (zero-width) or zero-exposure segments, for which $I_i x - M_i \ln x$ is +unbounded below. `optimize_skyline` now rejects such an over-segmented grid up front +(and validates strictly increasing boundaries and positive per-segment exposure as +backstops) rather than masking the degeneracy with the pooled warm-start fallback, +so the guarantee is enforced instead of assumed. ## Options considered diff --git a/kb/issues/N-coalescent-skyline-grid-validation-incomplete.md b/kb/issues/N-coalescent-skyline-grid-validation-incomplete.md index 104052c1..7ea7c6de 100644 --- a/kb/issues/N-coalescent-skyline-grid-validation-incomplete.md +++ b/kb/issues/N-coalescent-skyline-grid-validation-incomplete.md @@ -1,23 +1,23 @@ # Skyline grid construction accepts invalid endpoint arrays -> **Partially obsolete.** The `SkylineCostFunction` and the externally supplied -> `log_tc` array are gone: the segment grid is now derived internally from merger -> quantiles and the tree's time span (`merger_quantile_boundaries()`), and the -> optimizer is a convex Newton solve. The remaining validation gap is narrower — -> internal `boundaries` indexing still lacks a typed nonempty/ordering guard. See -> [decisions/coalescent-skyline-convex-inverse-tc.md](../decisions/coalescent-skyline-convex-inverse-tc.md). +> **Validation added; only the typed-grid refactor remains.** `optimize_skyline` +> now rejects an over-segmented grid up front (requested segments exceeding the +> number of distinct merger times) and validates strictly increasing boundaries and +> positive per-segment exposure as backstops, erroring with an actionable message +> instead of masking the degeneracy with the pooled warm-start fallback +> ([packages/treetime/src/coalescent/skyline.rs](../../packages/treetime/src/coalescent/skyline.rs)). +> The skyline tests exercise the production `optimize_skyline` objective directly and +> assert the invariant and the over-segmentation error. The core defect (silent +> acceptance of invalid segmentations) is resolved; what remains is the optional O1 +> refactor that moves these inline guards into a single validated skyline-grid type. +> See [decisions/coalescent-skyline-convex-inverse-tc.md](../decisions/coalescent-skyline-convex-inverse-tc.md). > The original text below predates the rewrite. Skyline construction indexes the first and last time-grid values before a typed boundary proves that the time and log-$T_c$ arrays are nonempty, equal in length, finite, and strictly ordered. Objective tests also exercise reconstructed formula fragments instead of the production cost function. -## Potential solutions +## Remaining work -- O1. Parse both arrays into one validated skyline-grid type whose constructor establishes every length, finiteness, and ordering invariant. -- O2. Repeat guards in each constructor and objective entry point. This can reject invalid arrays but permits the guards to drift and leaves invalid intermediate states representable. - -## Recommendation - -Use O1: parse both arrays into a validated skyline-grid type before endpoint indexing, and test the production objective directly. This validation is independent of the open extrapolation, quadrature, and simplex policies. +- O1. Parse the segment grid into one validated skyline-grid type whose constructor establishes every length, finiteness, and ordering invariant, replacing the current inline guards. This is a structural cleanup; the invariants themselves are already enforced. ## Related issues diff --git a/packages/treetime/src/coalescent/__tests__/test_skyline.rs b/packages/treetime/src/coalescent/__tests__/test_skyline.rs index 4afe2da0..c62f5288 100644 --- a/packages/treetime/src/coalescent/__tests__/test_skyline.rs +++ b/packages/treetime/src/coalescent/__tests__/test_skyline.rs @@ -9,7 +9,11 @@ mod tests { use maplit::btreemap; use treetime_io::dates_csv::{DateConstraint, DatesMap}; use treetime_io::nwk::nwk_read_str; + use treetime_utils::assert_error; + // Two internal nodes at distinct times (root, internal1), so two distinct merger + // times. The two lineage-count intervals are coarser than a two-way split, so this + // tree supports only a single skyline segment; it is used for the degenerate cases. const SMALL_TREE_NWK: &str = "((leaf1:1.0,leaf2:1.0)internal1:1.0,leaf3:1.0)root:1.0;"; fn small_tree_dates() -> DatesMap { @@ -22,37 +26,77 @@ mod tests { } } + // A balanced tree whose internal nodes sit at three distinct times (2000, 2002, + // 2004), giving three distinct merger times and enough lineage-count resolution to + // support a genuine multi-segment skyline. + const MEDIUM_TREE_NWK: &str = "(((a:1,b:1)ab:1,(c:1,d:1)cd:1)abcd:1,((e:1,f:1)ef:1,(g:1,h:1)gh:1)efgh:1)root:1;"; + const MEDIUM_TREE_DISTINCT_MERGERS: usize = 3; + + fn medium_tree_dates() -> DatesMap { + btreemap! { + "root".to_owned() => Some(DateConstraint::exact(2000.0)), + "abcd".to_owned() => Some(DateConstraint::exact(2002.0)), + "efgh".to_owned() => Some(DateConstraint::exact(2002.0)), + "ab".to_owned() => Some(DateConstraint::exact(2004.0)), + "cd".to_owned() => Some(DateConstraint::exact(2004.0)), + "ef".to_owned() => Some(DateConstraint::exact(2004.0)), + "gh".to_owned() => Some(DateConstraint::exact(2004.0)), + "a".to_owned() => Some(DateConstraint::exact(2006.0)), + "b".to_owned() => Some(DateConstraint::exact(2007.0)), + "c".to_owned() => Some(DateConstraint::exact(2008.0)), + "d".to_owned() => Some(DateConstraint::exact(2009.0)), + "e".to_owned() => Some(DateConstraint::exact(2010.0)), + "f".to_owned() => Some(DateConstraint::exact(2011.0)), + "g".to_owned() => Some(DateConstraint::exact(2012.0)), + "h".to_owned() => Some(DateConstraint::exact(2013.0)), + } + } + + fn medium_graph() -> Result { + helpers::create_graph_with_dates(MEDIUM_TREE_NWK, &medium_tree_dates()) + } + #[test] fn test_optimize_skyline_returns_result() -> Result<(), Report> { - let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; + let graph = medium_graph()?; let params = SkylineParams { - n_points: 5, + n_points: 3, ..SkylineParams::default() }; let result = optimize_skyline(&graph, ¶ms)?; - assert_eq!(result.segment_boundaries.len(), 6); - assert_eq!(result.tc_values.len(), 5); + assert_eq!(result.segment_boundaries.len(), 4); + assert_eq!(result.tc_values.len(), 3); assert!(result.log_likelihood.is_finite()); Ok(()) } #[test] - fn test_optimize_skyline_tc_values_positive() -> Result<(), Report> { - let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; + fn test_optimize_skyline_boundaries_strictly_increasing() -> Result<(), Report> { + // The reported grid must be strictly ordered, and the segment count must not + // exceed the number of distinct merger times (the every-segment-owns-a-merger + // invariant the optimizer relies on to keep each Tc finite). + let graph = medium_graph()?; let params = SkylineParams { - n_points: 5, + n_points: MEDIUM_TREE_DISTINCT_MERGERS, ..SkylineParams::default() }; let result = optimize_skyline(&graph, ¶ms)?; - for &tc in &result.tc_values { + assert!( + result.tc_values.len() <= MEDIUM_TREE_DISTINCT_MERGERS, + "segment count {} must not exceed the {MEDIUM_TREE_DISTINCT_MERGERS} distinct merger times", + result.tc_values.len() + ); + for pair in result.segment_boundaries.windows(2) { assert!( - tc > 0.0 && tc.is_finite(), - "Tc segment value must be positive and finite, got {tc}" + pair[1] > pair[0], + "segment boundaries must be strictly increasing, got {} then {}", + pair[0], + pair[1] ); } @@ -60,57 +104,81 @@ mod tests { } #[test] - fn test_optimize_skyline_tc_distribution_evaluates() -> Result<(), Report> { + fn test_optimize_skyline_rejects_oversegmentation() -> Result<(), Report> { + // Requesting more segments than distinct merger times cannot give every segment + // its own merger, so the optimizer rejects it up front with an actionable message + // instead of emitting zero-exposure segments whose Tc is set by the smoothing prior. let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; let params = SkylineParams { n_points: 5, ..SkylineParams::default() }; - let result = optimize_skyline(&graph, ¶ms)?; + assert_error!( + optimize_skyline(&graph, ¶ms), + "Skyline requested 5 segments but the tree has only 2 distinct merger times; each segment must own at least one merger. Reduce --n-skyline to at most 2." + ); - let t_min = result.segment_boundaries[0]; - let t_max = result.segment_boundaries[result.segment_boundaries.len() - 1]; - let t_mid = f64::midpoint(t_min, t_max); + Ok(()) + } - assert!(result.tc_distribution.eval(t_min)? > 0.0); - assert!(result.tc_distribution.eval(t_mid)? > 0.0); - assert!(result.tc_distribution.eval(t_max)? > 0.0); + #[test] + fn test_optimize_skyline_rejects_zero_exposure_segment() -> Result<(), Report> { + // The distinct-merger count alone does not guarantee positive exposure: on this + // coarse tree a two-way split leaves the oldest segment with no lineage-count + // interval, which the exposure backstop rejects rather than masking with the + // pooled warm-start fallback. + let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; + let params = SkylineParams { + n_points: 2, + ..SkylineParams::default() + }; + + assert_error!( + optimize_skyline(&graph, ¶ms), + "Skyline segment 0 has non-positive coalescent exposure (I = 0.000000e0); the requested 2 segments exceed what the tree's merger structure supports. Reduce --n-skyline." + ); Ok(()) } #[test] - fn test_optimize_skyline_larger_tree() -> Result<(), Report> { - const TREE_NWK: &str = "(((a:1,b:1)ab:1,(c:1,d:1)cd:1)abcd:1,((e:1,f:1)ef:1,(g:1,h:1)gh:1)efgh:1)root:1;"; - let dates = btreemap! { - "root".to_owned() => Some(DateConstraint::exact(2000.0)), - "abcd".to_owned() => Some(DateConstraint::exact(2002.0)), - "efgh".to_owned() => Some(DateConstraint::exact(2002.0)), - "ab".to_owned() => Some(DateConstraint::exact(2004.0)), - "cd".to_owned() => Some(DateConstraint::exact(2004.0)), - "ef".to_owned() => Some(DateConstraint::exact(2004.0)), - "gh".to_owned() => Some(DateConstraint::exact(2004.0)), - "a".to_owned() => Some(DateConstraint::exact(2006.0)), - "b".to_owned() => Some(DateConstraint::exact(2007.0)), - "c".to_owned() => Some(DateConstraint::exact(2008.0)), - "d".to_owned() => Some(DateConstraint::exact(2009.0)), - "e".to_owned() => Some(DateConstraint::exact(2010.0)), - "f".to_owned() => Some(DateConstraint::exact(2011.0)), - "g".to_owned() => Some(DateConstraint::exact(2012.0)), - "h".to_owned() => Some(DateConstraint::exact(2013.0)), + fn test_optimize_skyline_tc_values_positive() -> Result<(), Report> { + let graph = medium_graph()?; + let params = SkylineParams { + n_points: 3, + ..SkylineParams::default() }; - let graph = helpers::create_graph_with_dates(TREE_NWK, &dates)?; + let result = optimize_skyline(&graph, ¶ms)?; + + for &tc in &result.tc_values { + assert!( + tc > 0.0 && tc.is_finite(), + "Tc segment value must be positive and finite, got {tc}" + ); + } + + Ok(()) + } + + #[test] + fn test_optimize_skyline_tc_distribution_evaluates() -> Result<(), Report> { + let graph = medium_graph()?; let params = SkylineParams { - n_points: 10, + n_points: 3, ..SkylineParams::default() }; let result = optimize_skyline(&graph, ¶ms)?; - assert_eq!(result.tc_values.len(), 10); - assert!(result.log_likelihood.is_finite()); + let t_min = result.segment_boundaries[0]; + let t_max = result.segment_boundaries[result.segment_boundaries.len() - 1]; + let t_mid = f64::midpoint(t_min, t_max); + + assert!(result.tc_distribution.eval(t_min)? > 0.0); + assert!(result.tc_distribution.eval(t_mid)? > 0.0); + assert!(result.tc_distribution.eval(t_max)? > 0.0); Ok(()) } @@ -119,9 +187,9 @@ mod tests { fn test_skyline_reported_likelihood_matches_model_evaluation() -> Result<(), Report> { // The reported log-likelihood must equal the shared per-edge model cost // evaluated on the returned piecewise-constant Tc(t). - let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; + let graph = medium_graph()?; let params = SkylineParams { - n_points: 4, + n_points: 3, ..SkylineParams::default() }; @@ -134,6 +202,7 @@ mod tests { #[test] fn test_skyline_reported_likelihood_matches_per_edge_cost_for_polytomy() -> Result<(), Report> { + // A single-time polytomy has one distinct merger time, so exactly one segment. const TREE_NWK: &str = "(a:1,b:1,c:1,d:1)root:1;"; let dates = btreemap! { "root".to_owned() => Some(DateConstraint::exact(2000.0)), @@ -144,7 +213,7 @@ mod tests { }; let graph = helpers::create_graph_with_dates(TREE_NWK, &dates)?; let params = SkylineParams { - n_points: 2, + n_points: 1, ..SkylineParams::default() }; @@ -161,9 +230,9 @@ mod tests { fn test_skyline_beats_or_matches_constant_tc() -> Result<(), Report> { // The regularized skyline optimum should not have lower likelihood than the // best constant Tc (a skyline with all segments equal is a feasible point). - let graph = helpers::create_graph_with_dates(SMALL_TREE_NWK, &small_tree_dates())?; + let graph = medium_graph()?; let params = SkylineParams { - n_points: 4, + n_points: 3, stiffness: 1e-6, // near-zero smoothing: skyline free to fit each segment ..SkylineParams::default() }; diff --git a/packages/treetime/src/coalescent/skyline.rs b/packages/treetime/src/coalescent/skyline.rs index 254e6e2e..43d0ecee 100644 --- a/packages/treetime/src/coalescent/skyline.rs +++ b/packages/treetime/src/coalescent/skyline.rs @@ -6,6 +6,8 @@ use crate::payload::traits::TimetreeNode; use eyre::Report; use log::info; use ndarray::Array1; +use ordered_float::OrderedFloat; +use std::collections::BTreeSet; use treetime_distribution::{Distribution, DistributionFormula}; use treetime_graph::edge::GraphEdge; use treetime_graph::graph::Graph; @@ -107,9 +109,25 @@ where let t_min = breakpoints[0]; let t_max = breakpoints[breakpoints.len() - 1]; + // Every segment must own at least one merger for the `-Mᵢ ln xᵢ` barrier to keep + // its `Tc` finite, so the requested segment count cannot exceed the number of + // distinct merger times. Reject an over-segmented grid up front with a clear, + // actionable message rather than emitting duplicate boundaries or zero-exposure + // segments (which the structural guards below would otherwise reject). + let n_distinct_mergers = count_distinct_merger_times(&edges); + if params.n_points > n_distinct_mergers { + return make_error!( + "Skyline requested {} segments but the tree has only {n_distinct_mergers} distinct merger times; \ + each segment must own at least one merger. Reduce --n-skyline to at most {n_distinct_mergers}.", + params.n_points + ); + } + let boundaries = merger_quantile_boundaries(&edges, t_min, t_max, params.n_points); + validate_segment_boundaries(&boundaries, params.n_points)?; let (i_seg, m_seg) = accumulate_segment_terms(precomputed.lineage_counts(), &edges, &boundaries); + validate_segment_exposure(&i_seg, params.n_points)?; let x = solve_inverse_tc(&i_seg, &m_seg, params.stiffness, params.tolerance, params.max_iter)?; let tc_values = Array1::from_iter(x.iter().map(|&xi| 1.0 / xi)); @@ -182,6 +200,59 @@ fn merger_quantile_boundaries(edges: &[CoalescentEdgeData], t_min: f64, t_max: f boundaries } +/// Counts the number of distinct merger (parent) times among the collected edges. +/// +/// Coincident mergers (a polytomy, or ties) count once: boundaries cannot separate +/// events at the same time, so they cap the number of data-supported segments. +fn count_distinct_merger_times(edges: &[CoalescentEdgeData]) -> usize { + edges + .iter() + .map(|edge| OrderedFloat(edge.parent_time().value())) + .collect::>() + .len() +} + +/// Rejects a degenerate segmentation whose boundaries are not strictly increasing. +/// +/// `merger_quantile_boundaries` clamps interior candidates to be non-decreasing, so +/// when the requested segment count exceeds the number of distinct merger times it +/// emits duplicate (zero-width) boundaries. A zero-width segment carries no data, so +/// its `Tc` would be set only by the smoothing prior. Non-finite boundaries (e.g. a +/// NaN merger time) also fail the `>` comparison and are rejected here. +fn validate_segment_boundaries(boundaries: &[f64], n_points: usize) -> Result<(), Report> { + for pair in boundaries.windows(2) { + if !(pair[1] > pair[0]) { + return make_error!( + "Skyline produced a degenerate segmentation: boundary {:.6e} is not strictly less than the next \ + boundary {:.6e}. The requested {n_points} segments exceed the tree's distinct merger times; \ + reduce --n-skyline.", + pair[0], + pair[1] + ); + } + } + Ok(()) +} + +/// Rejects a segmentation with a zero-exposure segment. +/// +/// A segment with a merger but no coalescent exposure has `M_i > 0` and `I_i = 0`, +/// for which `I_i x - M_i ln x` is unbounded below. The guarantee that every segment +/// owns mergers holds only when the requested segment count does not exceed the +/// tree's distinct merger times, so a non-positive `I_i` signals an over-segmented +/// grid the data cannot support. +fn validate_segment_exposure(i_seg: &[f64], n_points: usize) -> Result<(), Report> { + for (i, &exposure) in i_seg.iter().enumerate() { + if !(exposure > 0.0) { + return make_error!( + "Skyline segment {i} has non-positive coalescent exposure (I = {exposure:.6e}); the requested \ + {n_points} segments exceed what the tree's merger structure supports. Reduce --n-skyline.", + ); + } + } + Ok(()) +} + /// Accumulates the per-segment pairwise-rate integral `Iᵢ` and merger count `Mᵢ`. /// /// `Iᵢ` sums, over lineage-count intervals whose midpoint falls in segment `i`, the diff --git a/packages/treetime/src/commands/timetree/args.rs b/packages/treetime/src/commands/timetree/args.rs index f448d346..aaccd2b0 100644 --- a/packages/treetime/src/commands/timetree/args.rs +++ b/packages/treetime/src/commands/timetree/args.rs @@ -122,7 +122,10 @@ pub struct TreetimeTimetreeArgs { /// /// When set, TreeTime will find the optimal Tc using Brent's method. This is similar /// to Python v0's `--coalescent=opt`, but uses a closed-form solution. - #[cfg_attr(feature = "clap", clap(long, conflicts_with = "coalescent", conflicts_with = "coalescent_skyline"))] + #[cfg_attr( + feature = "clap", + clap(long, conflicts_with = "coalescent", conflicts_with = "coalescent_skyline") + )] pub coalescent_opt: bool, /// Use skyline coalescent model instead of constant Tc.