Refactor/unified coalescent

kb/decisions/coalescent-analytic-tc-optimization.md

@@ -13,10 +13,13 @@ log-likelihood with Brent's method.
 candidate $T_c$.
 
 **v1 location**: `optimize_tc()` at
-[`packages/treetime/src/coalescent/optimize_tc.rs`](../../packages/treetime/src/coalescent/optimize_tc.rs). The
-$T_c$-independent integral $I$ reuses `compute_integral_merger_rate()` at
-[`packages/treetime/src/coalescent/integration.rs`](../../packages/treetime/src/coalescent/integration.rs)
-with a constant $T_c = 1$, rather than a dedicated integrator.
+[`packages/treetime/src/coalescent/optimize_tc.rs`](../../packages/treetime/src/coalescent/optimize_tc.rs).
+It is now a thin facade over the one-segment skyline — a constant $T_c$ is exactly
+an `n_points = 1` skyline — so the $I$/$M$ accumulation and likelihood reporting are
+the shared skyline machinery (see
+[coalescent-skyline-convex-log-tc.md](coalescent-skyline-convex-log-tc.md)), not a
+separate integrator. For one segment the skyline's decoupled optimum is precisely
+$T_c^{*}=I/M$.
 
 ## Background
 
@@ -75,8 +78,11 @@ skyline case.
   where Brent's bracket bit.
 - The skyline path (piecewise-constant $T_c$) was subsequently given the same
   analytic convex treatment; see
-  [coalescent-skyline-convex-inverse-tc.md](coalescent-skyline-convex-inverse-tc.md).
-- On a degenerate tree ($M \le 0$ or non-finite $I$) there is no hardcoded initial
-  $T_c$: the caller falls back to the previous round's $T_c$, then the user value,
-  then no prior (see
+  [coalescent-skyline-convex-log-tc.md](coalescent-skyline-convex-log-tc.md). The
+  two were then unified: constant $T_c$ is implemented as the one-segment skyline,
+  leaving a single optimizer.
+- On a degenerate tree ($M \le 0$ or non-finite/non-positive $I$) there is no
+  hardcoded initial $T_c$ and no fallback: the shared solve returns an error, which
+  the pipeline propagates to stop the run rather than substituting an invented
+  timescale (see
   [issues/N-coalescent-initial-tc-hardcoded-fallback.md](../issues/N-coalescent-initial-tc-hardcoded-fallback.md)).

kb/decisions/coalescent-skyline-convex-log-tc.md

@@ -0,0 +1,138 @@
+# Skyline Tc solved by exact convex Newton in log Tc
+
+v1 estimates the piecewise-constant skyline $T_c(t)$ by an exact convex
+optimization in the log time scale $z_i = \ln T_{c,i}$. This replaced two earlier
+approaches: the original v1 derivative-free search over $\log T_c$ knots, and an
+intermediate exact convex solve in the inverse time scale $x_i = 1/T_{c,i}$.
+
+**Type**: Implementation change (parametrization and solver; also diverges from
+v0's optimizer, extending the constant-$T_c$ decision to the skyline).
+
+**Original v1**: `optimize_skyline()` used `argmin`'s `NelderMead` over a
+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.
+
+**Intermediate v1**: exact convex Newton in $x_i = 1/T_{c,i}$ minimizing
+$\sum_i (I_i x_i - M_i \ln x_i) + \frac{\gamma}{2}\sum_i (x_{i+1}-x_i)^2$, with
+positivity-preserving damped steps.
+
+**v0**: `MergerModel.optimize_skyline()` at
+[`packages/legacy/treetime/treetime/merger_models.py#L281-L318`](../../packages/legacy/treetime/treetime/merger_models.py) uses
+`scipy.optimize.minimize(method='SLSQP')`.
+
+**v1 now**: `optimize_skyline()` at
+[`packages/treetime/src/coalescent/skyline.rs`](../../packages/treetime/src/coalescent/skyline.rs).
+
+## Background
+
+Splitting the tree into segments with constant $T_c$, the constant-$T_c$ likelihood
+(see [coalescent-analytic-tc-optimization.md](coalescent-analytic-tc-optimization.md))
+applies per segment. Writing $z_i = \ln T_{c,i}$ (so the coalescent rate is
+$1/T_{c,i} = e^{-z_i}$) and penalizing squared log-fold-changes of $T_c$, the cost
+to minimize is
+
+$$
+C(z) = \sum_i \big(I_i\,e^{-z_i} + M_i z_i\big) + \frac{\gamma}{2}\sum_i (z_{i+1}-z_i)^2,
+$$
+
+with $I_i = \int_{\text{seg }i} k(k-1)/2\,dt$ and merger count $M_i$. Each data term
+is convex ($e^{-z_i}$ convex, $M_i z_i$ linear) and the penalty is a PSD quadratic,
+so $C$ has a unique minimizer, found by Newton's method on the symmetric tridiagonal
+Hessian (Thomas algorithm, $O(n)$ per step), warm-started from the decoupled optimum
+$z_i = \ln(I_i/M_i)$ and globalized with an Armijo line search. This is a
+Poisson-likelihood / Gaussian-smoothing MAP estimate — the frequentist analog of the
+Bayesian GMRF skygrid, which is likewise formulated on $\log$ effective population
+size.
+
+## Rationale
+
+- **Convex with a unique optimum.** In $z = \ln T_c$ the data term is convex and the
+  penalty is a PSD quadratic, so there are no local optima and no simplex or bracket
+  to tune. Optimizing $T_c$ directly is not convex.
+- **Scale-independent smoothing.** The penalty charges $(z_{i+1}-z_i)^2 =
+  \ln(T_{c,i+1}/T_{c,i})^2$, i.e. squared log-fold-changes, so the stiffness is
+  dimensionless and its effect does not depend on the absolute scale of $T_c$ (or the
+  time unit). In the earlier $x = 1/T_c$ parametrization the penalty
+  $(x_{i+1}-x_i)^2$ carried units of time$^{-2}$, so the same stiffness meant
+  different smoothing at different scales.
+- **Positivity for free.** $T_c = e^z > 0$ for any real $z$, so the Newton step needs
+  no positivity cap — only an Armijo line search for global convergence (the data
+  term is no longer quadratic in the optimization variable, so Newton is iterative).
+- **Exact and cheap.** A few $O(n)$ Newton steps versus many derivative-free
+  objective evaluations, each rebuilding the merger-rate integral.
+- **Self-consistent likelihood.** $I_i$ and $M_i$ use the same interval-midpoint and
+  node-time attribution as `CoalescentModel`, so the analytic optimum is the true
+  maximizer of the model-evaluated likelihood, and the reported LH matches
+  `compute_coalescent_total_lh`.
+- **Simpler.** Removes the Nelder-Mead executor, the initial-simplex construction,
+  the boundary penalty, and the piecewise-linear knot interpolation.
+
+## Equal-width segment boundaries
+
+Segments are equally spaced over the tree's time span (`equal_width_boundaries`, a
+`linspace`). An earlier design placed boundaries at merger-time quantiles so every
+segment would own mergers ($M_i>0$), keeping the linear $M_i z_i$ term active and
+each per-segment optimum finite even without smoothing. That guarantee turned out to
+be illusory and costly:
+
+- **It is conditional and unenforced.** By pigeonhole, quantile bins can give every
+  segment a merger only when `n_points` does not exceed the number of merger times —
+  which depends on topology (polytomies) and bad-branch exclusion, is unknown when
+  the user picks `n_points`, and is already violated in the test suite (5 segments on
+  a tree with 2 mergers). Outside that regime quantile bins also leave empty
+  segments, so well-posedness rests on `stiffness > 0` regardless.
+- **Ties produce zero-width segments.** Coincident merger times (polytomies) collapse
+  adjacent quantile boundaries into degenerate zero-width segments.
+- **Variable widths make the stiffness scale-ambiguous.** On a non-uniform grid
+  $\sum (z_{i+1}-z_i)^2$ penalizes fold-change per segment *index*, ignoring the time
+  gap. On a uniform grid it is a consistent discretization of
+  $\int (d\ln T_c/dt)^2\,dt$, so the stiffness has a clean, grid-independent meaning.
+
+Equal width wins on tie-robustness, on a well-defined stiffness, and on simplicity
+(the sort/quantile/clamp/tie logic is gone), and it matches the skygrid we emulate,
+which uses a regular grid and lets the smoothing carry low-data cells. The cost is
+that merger-sparse regions (typically near the root and the tips) can yield empty
+segments ($M_i=0$); their $z_i$ is then pinned by the smoothing prior rather than the
+data, so **`stiffness > 0` is required for `n_points > 1`** (enforced in
+`solve_log_tc`) to keep the Hessian positive-definite and every $z_i$ finite — a lone
+data term $I_i e^{-z_i}$ would otherwise drive $T_c\to\infty$. With positive stiffness
+the PD tridiagonal Hessian keeps the solve well-posed and an empty segment's $T_c$
+interpolates its neighbours instead of diverging. (The historical ebola/20 blow-up,
+$T_c\approx 3\times10^{150}$, occurred with equal-*time* boundaries under the earlier
+$x=1/T_c$ parametrization before reliable smoothing, and does not recur.)
+
+## Options considered
+
+1. **Port v0's SLSQP over log Tc.** Rejected: retains a general-purpose numerical
+   search with optimizer-dependent local optima; the constant-$T_c$ case already
+   showed the numerical search landing on wrong values.
+2. **Newton in $x=1/T_c$ (intermediate v1).** Convex and exact with a tridiagonal
+   Hessian, but the smoothing penalty $(x_{i+1}-x_i)^2$ is scale-dependent and the
+   step needs a positivity cap to keep $x_i>0$.
+3. **Newton in $z=\ln T_c$, equal-width boundaries (chosen).** Convex and exact with
+   the same tridiagonal Hessian, plus scale-independent smoothing and unconstrained
+   positivity. The data term becomes non-quadratic (exponential), so Newton iterates
+   under an Armijo line search rather than converging in one step. Boundaries are
+   uniform rather than merger-quantile (see [above](#equal-width-segment-boundaries));
+   empty segments are handled by requiring `stiffness > 0`.
+
+## Impact
+
+- **Constant $T_c$ is unified as the one-segment case of this solve.** `optimize_tc`
+  is a facade calling `optimize_skyline` with `n_points = 1`, where the smoothing is
+  inert and the decoupled optimum is the closed form $T_c^{*}=I/M$. There is now a
+  single coalescent optimizer. A whole-tree degeneracy guard (`m_tot > 0` and
+  positive finite `i_tot`) lives in `optimize_skyline` and is shared by both cases:
+  it errors on a tree with no time span or no mergers rather than letting the
+  per-segment pooled fallback silently return $T_c = 1$.
+- Skyline estimates no longer match v0 numerically; they are the exact optimum of
+  the (regularized) piecewise-constant Kingman likelihood.
+- The stiffness parameter penalizes adjacent $\ln T_c$ differences
+  (log-fold-changes), so it is dimensionless and scale-independent. Its numeric
+  scale differs from the intermediate $1/T_c$ implementation; the CLI default is
+  retained at `2.0` and reinterpreted on the log scale.
+- Supersedes the Nelder-Mead optimizer and simplex-initialization concerns in
+  [issues/M-timetree-skyline-nelder-mead-optimizer.md](../issues/M-timetree-skyline-nelder-mead-optimizer.md)
+  and [issues/N-coalescent-skyline-simplex-initialization-undecided.md](../issues/N-coalescent-skyline-simplex-initialization-undecided.md).
+  Boundary out-of-domain evaluation still uses constant extrapolation
+  ([issues/N-coalescent-skyline-extrapolation-policy-undecided.md](../issues/N-coalescent-skyline-extrapolation-policy-undecided.md)).