fix: reject over-segmented skyline grids

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
 

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
 

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,95 +26,159 @@ 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<GraphTimetree, Report> {
+    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, &params)?;
 
-    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, &params)?;
 
-    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]
       );
     }
 
     Ok(())
   }
 
   #[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, &params)?;
+    assert_error!(
+      optimize_skyline(&graph, &params),
+      "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, &params),
+      "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, &params)?;
+
+    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, &params)?;
 
-    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()
     };