From bbae2bdbf68a97ee837fd55fdf96c0c6d216b8ca Mon Sep 17 00:00:00 2001 From: Vince Buffalo Date: Sat, 11 Jul 2026 21:11:10 -0700 Subject: [PATCH] =?UTF-8?q?docs(notes):=20Stan=20re-implementation=20of=20?= =?UTF-8?q?Garki=20ctl=5Fbb=20=E2=80=94=20ODE-fit=20benchmark?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A faithful Stan port of the ctl_bb ODE model (70 states, beta-binomial, full 60.9-yr burn-in), validated at camdl's MAP. Measured timing: 0.055 s forward solve, ~1.2 s per NUTS gradient, full fit extrapolating to ~7 hr– few days vs camdl mh-ode's ~28 min. NUTS-through-the-burn-in is impractical (reproduces garki's F21 independently). Records a standalone lever: Stan's compiled RK4 solve (0.055 s) beats camdl's resolved-expression evaluator per solve — worth chasing for camdl's own ODE fits. --- .../notes/2026-07-12-stan-ctl-bb-benchmark.md | 208 ++++++++++++++++++ 1 file changed, 208 insertions(+) create mode 100644 docs/dev/notes/2026-07-12-stan-ctl-bb-benchmark.md diff --git a/docs/dev/notes/2026-07-12-stan-ctl-bb-benchmark.md b/docs/dev/notes/2026-07-12-stan-ctl-bb-benchmark.md new file mode 100644 index 00000000..fe3ddaf1 --- /dev/null +++ b/docs/dev/notes/2026-07-12-stan-ctl-bb-benchmark.md @@ -0,0 +1,208 @@ +# Stan re-implementation of Garki `ctl_bb` — ODE-fit benchmark vs camdl + +Date: 2026-07-12 Project: camdl Tags: ode, stan, benchmark, performance, garki, +burn-in + +## Headline + +**NUTS-through-a-60-year ODE is impractical for `ctl_bb`.** A faithful Stan +model compiles and samples in the right region, but a full fit extrapolates to +**~7 hours (optimistic floor) to ~3–4 days (realistic) per chain**, vs camdl's +**~28 minutes** for the whole 6-chain fit. The cost is the ~60.9-year burn-in: +every NUTS gradient differentiates through a ~23,355-step RK4 tape, and NUTS +pays tens-to-thousands of gradients per iteration where camdl's Metropolis– +Hastings pays one forward solve. This independently reproduces garki's own +finding (friction-log F21). + +A **separately actionable sub-finding** is in §5: Stan's tight compiled RK4 +forward solve (0.055 s) is faster than camdl's resolved-expression evaluator per +solve. The wall-clock loss is entirely NUTS × burn-in, not Stan's engine — and +the per-solve gap is a standalone lever for camdl's own ODE fits. + +## 1. Model and faithfulness + +Translated verbatim from `models/ctl_bb.camdl`: + +- **70 ODE states** = 2 villages (kwaru, ajura) × 5 age bands × 7 compartments + (S_naive, E_naive, I_infectious, I_recovering, I_immune, S_immune, E_immune). +- Saturating FOI `h_v = g·(1 − exp(−Cscale·C_v(t)·(I_infectious_vil/N_vil)))`. +- Superinfection-corrected recovery `R1 = h/(exp(h/r1)−1)`, `R2` with + `r2 = 10·r1` (implemented with `expm1`; algebraically identical). +- Incubation `nu`, loss-of-infectivity `a1`, immunity acquisition `a2`. +- Aging across 5 bands at `1/age_dur`, widths [4,5,5,10,20] yr; birth = death + `delta` (population conserved). **Year = 365.2425 days**, camdl's + Gregorian-average constant (`caltime.rs` / `expander.ml`), so aging rate/day = + `1/(age_dur_yr·365.2425)`. +- **Forcing** `C_v(t)`: piecewise-constant, **left-continuous step** (value at + the largest knot ≤ t, clamped at ends), matching + `rust/crates/sim/src/propensity.rs::constant_value`. 129 shared knot dates, + per-village values, from `data/processed/mv_village_C.tsv`. +- **Observation**: beta-binomial over patent-prevalence cells. Patent + `p = (q1·I_infectious + q2·I_recovering + q3·I_immune)/N` (q1=q2=1, q3=0.7), + concentration `φ = 1/phi_inv`, `α = p·φ`, `β = (1−p)·φ`, `beta_binomial_lpmf`. +- **Fitted (4)**: g, r1, a2, phi_inv (matches `fits/ctl_bb.toml [estimate]` and + the recorded fit's `fit.meta.json`; the model file gives a1 a prior but the + fit **fixes** a1=0.002). Priors/bounds exactly as ctl_bb.toml: + g~lognormal(−2.50,0.8)∈[0.01,0.5]; r1~lognormal(−6.50,0.7)∈[5e-4,0.02]; + a2~lognormal(−8.50,1.0)∈[2e-5,2e-3]; phi_inv~half_normal(1.0)∈[1e-4,3.0]. +- **Solver**: hand-rolled **fixed-step RK4, dt = 1 day** (camdl + `[config] + dt=1.0`, `backend = ode`). Deliberately NOT `ode_rk45`: camdl + friction-log **F17** records that an adaptive solver restarts and takes tiny + steps at each of the ~130 forcing discontinuities, making it slower than + fixed-step RK4. Hand-rolled RK4 also gives a clean, predictable per-gradient + cost and the closest apples-to-apples with camdl. + +**Numeric faithfulness check.** Stan forward solve at camdl's recorded MAP +(g=0.0481, r1=0.001154, a2=0.000245, phi_inv=0.0900) reproduces the age +gradient: + +| age band | Stan pred p | observed p | +| -------- | ----------- | ---------- | +| a1_4 | 0.699 | 0.824 | +| a5_9 | 0.692 | 0.812 | +| a10_14 | 0.598 | 0.627 | +| a15_24 | 0.473 | 0.369 | +| a25_44 | 0.307 | 0.238 | +| overall | 0.554 | 0.574 | + +The declining-with-age patent prevalence (the model's signature mechanism) is +reproduced, and the _young-under / old-over_ residual matches the pattern +camdl's own report documents. This is a **structural-faithfulness** check, not +golden bit-parity (different sampler/RNG; standard Stan Jacobian on the bounded +transforms). + +**Simplifications**: none in the dynamics — full 60.9-yr burn-in, dt=1, all 5 +bands, both villages, no equilibrium shortcut. + +## 2. Environment / versions + +- **cmdstan 2.39.0** (built from source via `cmdstanpy.install_cmdstan`), + **cmdstanpy 1.3.0**, Python 3.12.11, Apple clang 17.0.0, macOS (M4 Max). +- **camdl 0.1.0+61551324** on PATH; the recorded `ctl_bb` fit used **camdl + 0.1.0+5c527817** (2026-07-10) — cited with that provenance below. + +## 3. Burn-in handling and its cost + +Origin **1910-01-01**; first observation **1970-11-17** = day 22235 = **60.88 +yr** burn-in; last obs 1973-12-11 = day 23355 = 63.94 yr total. At dt=1 the +faithful solve is **~23,355 RK4 steps** (×4 deriv evals/step ×70 states). The +Stan model integrates from 1910 with the seeded IC (S_naive=900, +I_infectious=100 in the youngest band per village), exactly as camdl does — no +equilibrium shortcut. The system reaches an annual limit cycle within a few +years (pre-1971 forcing is a repeating seasonal climatology), so the long +burn-in tail is wasted work in both frameworks; it is precisely the cost that +camdl's gh#396 "coarse burn-in dt" / equilibrium-warm-start work targets. + +## 4. Timing (measured, then extrapolated) + +All on M4 Max, cmdstan 2.39.0, dt=1 fixed-step RK4, single thread. + +### Measured anchors (real numbers) + +| quantity | value | how measured | +| ------------------------------------------- | ------------------------------------------ | ------------------------------------- | +| **forward solve** (double, no grad) | **0.055 s** | 50 `fixed_param` GQ evals (2.77 s) | +| **single reverse-mode gradient** | **1.0–1.36 s** | Stan "Gradient evaluation took …", ×2 | +| gradient / forward ratio | **~18–25×** | derived (memory-bound tape) | +| 1 warmup iter, max_depth=2 (≤4 leapfrog) | **2.2 s** | tiny run (real lower bound) | +| sampling iter, max_depth=5 (~10 leapfrog) | **9.1 s/iter** | completed 6-iter run | +| diagnostic (8 warmup + 6 sampling, depth≤5) | **163.3 s** (108.6 warmup + 54.7 sampling) | full short run | +| adaptation run (40 warmup, depth≤6) | killed at **222 s**, barely past warmup | ~60–90 s/iter with adapted step | + +The sampler moved to the right region even in the tiny run (g=0.047, r1=0.00133, +a2=0.000259, phi_inv=0.107 — all near camdl's MAP), so this is a cost problem, +not a wiring problem. + +### Extrapolation to a full fit (order-of-magnitude) + +Per-iteration cost ≈ (avg leapfrog steps) × (per-gradient cost). Gradient ≈ 1.2 +s. + +- **Optimistic floor — 10 leapfrog/iter**: 2000 × 10 × 1.2 s = **~6.7 + hr/chain**. (Stan's own printed estimate: ~5.6–7.6 hr for 2000.) +- **Realistic — ~128 leapfrog/iter** (depth ~7; the g–r1 ridge forces a small + adapted step ⇒ deep trees; camdl's own R̂≈2.2 documents the ridge): 2000 × 128 + × 1.2 s = **~3.6 days/chain**. + +Chains are independent, so with enough cores wall ≈ per-chain: **~7 hr (floor) +to ~3–4 days (realistic).** The killed 40-warmup run (≈60–90 s/iter) sits in the +realistic regime. + +## 5. Comparison to camdl (garki's recorded ctl_bb fit) + +Source: `results/fits/ctl_bb-b35f911a/.../run.json` (`wall_time_seconds`) + +`pmmh_summary.json`, produced by **camdl 0.1.0+5c527817**, 2026-07-10. + +| | camdl `ctl_bb` (mh-ode) | Stan `ctl_bb` (NUTS), extrapolated | +| ------------------- | -------------------------------------------------- | ------------------------------------- | +| sampler | adaptive MH, 1 forward solve/iter | NUTS, ~10–1000 gradient solves/iter | +| config | 6 chains × 2000 iters (burn_in 1000) | ~4 chains × ~2000 iters | +| **wall (full fit)** | **1682.8 s ≈ 28.0 min** | **~7 hr (floor) – ~3–4 days (ridge)** | +| slowdown vs camdl | 1× | **~15× (floor) – ~150–200× (ridge)** | +| convergence | R̂ g=2.22, r1=2.38, a2=1.97; phi_inv=1.06 (ESS 360) | not reached | +| ESS/sec (phi_inv) | 360 / 1682.8 = **0.21 ESS/s** | — | + +**camdl's 28 min ≠ converged.** Even camdl does not converge the ridge params in +28 min (R̂ ≈ 2 for g, r1, a2; only phi_inv, R̂=1.06, ESS=360, is usable). The +`ctl_bb` posterior is a genuinely hard g–r1 identifiability ridge with one stuck +chain. So the target is hard for _both_ samplers; NUTS's theoretical edge +(better ESS/iter on hard geometry) is swamped by the per-gradient burn-in cost. + +**Sub-finding — the solver gap (standalone camdl lever).** Stan's hand-rolled +forward solve (0.055 s) is _faster_ than camdl's implied per-mh-iteration +(1682.8 s / 2000 iters ≈ 0.84 s wall/iter across 6 parallel chains — mostly one +forward solve + MH overhead). Same algorithm (RK4 dt=1), same span: a tight +**compiled** RK4 beats camdl's **resolved-expression evaluator** per solve. If +camdl's per-solve cost closed toward Stan's, camdl's own ODE fits would speed up +roughly an order of magnitude — a lever independent of the NUTS-vs-MH question. +Needs its own measurement + profile (interpreted AST eval vs native codegen; cf. +the known `eval_resolved`-dominated sim profile). + +## 6. Honest caveats (not apples-to-apples) + +1. **Sampler.** camdl = random-walk adaptive MH (one forward solve/iter); Stan = + NUTS (many gradient solves/iter). The dominant, structural difference. +2. **Autodiff cost is the burn-in tape, not the model.** The 1.2 s gradient is + ~20× the forward solve because reverse-mode over a 23,355-step tape is + memory-bound. Shrink the tape and this collapses. +3. **Solver.** Both fixed-step RK4 dt=1, but camdl integrates via its + resolved-expression evaluator and Stan via compiled RK4; trajectories are not + bit-identical (dynamics verified at the MAP, not golden parity). No + adaptive-tolerance solver was used (rejected per F17). +4. **Convergence.** No full Stan fit ran, so no Stan R̂/ESS/divergence numbers; + the extrapolation is order-of-magnitude. +5. **The real lever is the burn-in, for both.** A fair Stan version would need + the same fix camdl is building: solve the **endemic / periodic steady state** + as the IC (a stroboscopic fixed point) instead of integrating 60 yr through + every gradient, or a coarse burn-in dt on the unscored warm-up (gh#396). With + the burn-in tape cut ~14× the gradient drops toward ~0.1 s and the full-fit + extrapolation falls into the ~0.5–4 hr range — competitive, but still an + unproven equilibrium-IC engineering exercise, and still fighting the g–r1 + ridge. + +## 7. Bottom line + +A faithful Stan re-implementation of `ctl_bb` **works** (dynamics validated at +camdl's MAP), but **NUTS through the 60.9-year burn-in is not viable as-is**: +0.055 s forward solve, ~1.2 s per gradient, full fit extrapolating to **~7 hr to +a few days** — vs camdl's mh-ode **~28 min** (which itself doesn't converge the +g–r1 ridge). Stan buys nothing here until the burn-in is killed by an +equilibrium / periodic steady-state IC (or a coarse warm-up dt) — exactly the +lever camdl's gh#396 is already pulling. Stan's per-solve engine is actually +fast; the wall-clock loss is entirely NUTS × burn-in, reproducing garki's F21 +conclusion from an independent toolchain. + +## Next + +- **Chase the §5 solver gap** — measure camdl's own `ctl_bb` forward solve + directly, profile where the per-step cost goes (resolved-expression eval vs + native), and scope whether a ~10× is reachable for camdl's ODE fits. +- (Optional) build the equilibrium-IC Stan variant to test the ~0.5–4 hr + hypothesis from the Stan side — a cross-check on gh#396's payoff. + +The Stan model + drivers (`ctl_bb.stan`, `ctl_bb_functions.stan`, +`forward_only.stan`, `prep_data.py`, `run_forward.py`, `run_fit.py`) were built +in a scratch working dir; the data reference is external +(`camdl-garki/data/processed/mv_village_C.tsv`), so this note is the durable +record rather than a committed reproduction.