fix(inference): IF2 Algorithm-1 conformance — perturb before the step (gh#365), per-particle x₀ (gh#364)

[vsbuffalo]

Two IF2 correctness fixes, one commit each, both in
rust/crates/sim/src/inference/if2.rs. Both issues were tagged SUSPECTED; both
reproduce. Together they bring the per-iteration loop into conformance with
Algorithm 1 of Ionides et al. (2015).

Closes #365
Closes #364

The reference, read rather than recalled

Ionides, E. L., Nguyen, D., Atchadé, Y., Stoev, S. & King, A. A. (2015).
Inference for dynamic and latent variable models via iterated, perturbed Bayes
maps. PNAS 112(3), 719–724.
doi:10.1073/pnas.1410597112.
Algorithm 1:

Θ^F_{0,j} ~ h_0(θ | Θ^{m-1}_j; σ_m)                       for j in 1:J
X^F_{0,j} ~ f_{X_0}(x_0; Θ^F_{0,j})                       for j in 1:J
for n in 1:N
  Θ^P_{n,j} ~ h_n(θ | Θ^F_{n-1,j}, σ_m)
  X^P_{n,j} ~ f_{X_n|X_{n-1}}(x_n | X^F_{n-1,j}; Θ^P_{n,j})
  w_{n,j}   = f_{Y_n|X_n}(y*_n | X^P_{n,j}; Θ^P_{n,j})
  resample (Θ, X) jointly

pomp 6.4.0.2 (kingaa/pomp@0eaf3c018f4d5857741fe29b88ef755f06b64220),
R/mif2.R lines 461–493, pomp:::mif2_pfilter:

for (nt in seq_len(ntimes)) {
  ## perturb parameters
  pmag <- cooling.fn(nt,mifiter)$alpha*rw.sd[,nt]
  params <- .Call(P_randwalk_perturbation,params,pmag)
  tparams <- partrans(object,params,dir="fromEst",.gnsi=gnsi)
  ## get initial states
  if (nt == 1L) { x <- rinit(object,params=tparams) }
  ## advance the state variables according to the process model
  X <- rprocess(object,x0=x,t0=times[nt],times=times[nt+1],params=tparams,...)
  ## determine the weights
  weights <- dmeasure(object,y=...,x=X,times=times[nt+1],params=tparams,log=TRUE,...)

tparams is an npars × Np matrix and src/mif2.c's
randwalk_perturbation jitters it column by column, so both facts below are
visible in one loop: perturb precedes rprocess, and rinit receives the
per-particle perturbed matrix.

gh#365 — perturb before the process step

Real. Within each observation window IF2 ran propagate → fold → perturb →
weight, so X_n was simulated at Θ^F_{n-1} while g(y_n | X_n; θ) was scored
at the freshly perturbed Θ^P_n. Harmless (a vanishing phase offset) for a
parameter living only in the process or only in the measurement; a genuine
coupling error for one living in both.

The fix moves the observation-time perturbation to the top of the loop.
global_step moves with it, so the cooling schedule is untouched — the SD
applied at observation k is still per_step_cooling^(k+1).

Red/green is exact, not statistical: a mock ProcessModel stamps the θ it was
stepped with into the particle state, and the observation model compares it
against the θ it is handed for weighting.

before: 240/240 particle-observations were weighted at a θ different from the
        θ that generated their state. First offender: obs 0 stepped
        θ=0.538370508862 but weighted θ=1.498434206341
after:  ok. 1 passed

docs/methods/particle-methods.md:143 already documented the correct order;
docs/inference.md:526 documented the implemented (wrong) one. Classified
doc-vs-code with the code as the loser — code fixed, docs/inference.md
updated.

gh#364 — every particle gets its own initial state

Real, and it silently disables a documented feature. initial_state was
evaluated once per iteration from the iteration-mean parameters and copied to
every particle, so the per-particle t=0 perturbation never reached the initial
compartment counts.

The issue asked to verify the precondition first: ivp parameters do not
re-enter step. Confirmed — rg -n '\.ivp\b' rust/crates/sim/src/inference/
returns a single engine site, if2.rs:522, which only skips the
observation-time perturbation. A pure-IC ivp parameter (the documented
S0/E0/I0 case, and simplex-group compositions) reaches the data through
x₀ and nothing else. With one shared x₀ the weights are independent of it,
resampling is a blind subsample, and its filter mean is never selected.

ic_free = true validates that such a parameter exists precisely to guarantee
the t=0 spread the engine was discarding, and fit/mod.rs:493 prints "initial
state spread from ivp params: [...]" to the user for spread that did not exist.

Measured on a chain-binomial SIR whose i0 sets I₀ and appears in no rate and
in no observation model:

	before	after
mean weighted_var_ratio for i0	0.9978	0.9365
i0 estimate (start 0.01000, truth 0.00600)	0.01370	0.00716

weighted_var_ratio ≈ 1 is the issue's claim measured directly: the weights are
blind to i0. After the fix i0 locks on within two iterations and holds for
the remaining 28.

One trap worth recording, because it made an earlier version of this test pass
against the broken code: an ivp parameter perturbed symmetrically on the logit
scale but averaged back on the natural scale picks up a data-free Jensen
drift toward the midpoint of its declared bounds. In the first draft that
drift walked i0 straight through the truth and the test went green on a filter
that had learned nothing (weighted_var_ratio was 0.9986 the whole time). The
fixture now starts i0 on the far side of the truth from the bounds midpoint,
so the drift pushes away and only selection can close the gap. (pomp has the
same natural-scale averaging — mif2.R:323,
start <- apply(.paramMatrix,1L,mean) — so this is inherent to IF2 as
published, not a camdl divergence.)

Seam

ProcessModel::initial_state was already the producer being called; the fix
routes each particle through it instead of the swarm mean, so no second
implementation appears. PGAS's per-particle Binomial(N₀, θ) draw
(pgas.rs::csmc_as) is deliberately not reused: it exists because PGAS
needs a tractable initial-state density p(x₀|θ) for the complete-data
likelihood, it is keyed off IVPMappings built by finite-differencing a
&CompiledModel that ProcessModel only optionally exposes, and it would
inject Monte-Carlo variance Algorithm 1 does not ask for. IF2 needs a draw, not
a density — which is also exactly what pomp's rinit is here.

What is NOT verified

The IF2-vs-mif2 MLE cross-check the issue asks for. It needs R + pomp
driving the same synthetic data, which this PR does not run. What is verified is
the step order and the rinit argument against pomp's source and against
Algorithm 1, plus the two exact structural tests. The numerical cross-check
remains open — worth its own issue if you want it.

Baselines

gh#365 reorders per-particle RNG consumption (each stream now yields its
perturbation draws before its process draws), so every IF2 log-likelihood and
MLE moves. Expected and legitimate. gh#364 is RNG-neutral — initial_state
draws nothing, and it moved to after a perturbation that already ran in the same
order — so results move only for models whose initial_conditions mention an
estimated parameter; Explicit initial conditions are bit-identical.

No pinned IF2 baseline exists to update: gate_inference_baseline.rs has no IF2
arm, and every IF2 test in the suite asserts structure or a tolerance band, not
bit values. Content-addressed fit artifacts re-key automatically — the fit CAS
digest folds engine_version (cli/src/version.rs::VERSION_SHORT, which
carries the git hash), so no cached pre-fix result can be served under a
post-fix run_id.

Tests

Targeted suites, all green on the final tree (make test not run locally — CI
is the gate):

gh364_if2_per_particle_initial_state   2 passed
gh365_if2_perturb_before_step          1 passed
if2                                   12 passed
if2_cooling                            3 passed
gh216_cursor_firing                   15 passed
pf_watchdog                            5 passed   (incl. θ̂ thread-count invariance)
gate_licm_ab                           3 passed

Adjacent things noticed, not touched

• camdl collapses the particle cloud to its mean between IF2 iterations
  (current_params = param_means, then every particle is re-seeded from it).
  Algorithm 1 sets Θ^m_j = Θ^F_{N,j} — particle j carries its own value
  forward — and pomp passes paramMatrix from iteration to iteration without
  collapsing it (mif2.R:389). camdl's behaviour is deliberate and documented
  (docs/methods/particle-methods.md), but it is a third divergence from the
  published algorithm and it is what makes the Jensen drift above accumulate.
  Worth a decision, out of scope here.
• The bootstrap PF has no t=0 spread mechanism at all (particle_filter.rs
  copies one initial_state(params) to every particle, and θ is global there),
  so the ic_free precondition documented in SMCConfig cannot be met on that
  path the way it now can on IF2's.
• docs/methods/particle-methods.md carries a stale cooling exponent
  (σ_n = σ_0 · c^{2n/N} — the exponent-2 form gh#363 fixed) and a stale
  if2.rs:191 line reference.

---

CI round 2 — two cli-crate failures, both consequences of gh#365

Both were verified against origin/main first (if2.rs reverted, everything
else identical): both pass on main and fail here, so neither is a
pre-existing flake. Neither assertion was weakened and no threshold was
relaxed.

synthetic_fit_grid::v2_if2_chains_diverge_at_iter_0_when_no_starts_from — class (a)

The recorded quantity's relationship to the chain start changed.
write_chain_starts' own docstring (runner.rs:2583) already says
parameter_traces.tsv iteration 0 is "post-first-filter (already perturbed)"
and cannot answer "did the chains get distinct starts?" — chain_starts.tsv is
the instrument for that, and assertion 1 against it still passes.

Not (b). Per-chain init spread is intact. Measured on the branch build:

chain	start_beta	iter0_beta
1	0.0366	1.9054
2	0.9311	0.9385
3	0.0214	2.2312 ← min start
4	0.8730	0.7535
5	0.0509	2.0929
6	0.0553	1.5353
7	0.0604	1.3580
8	1.7323	1.6384 ← max start

Chains starting near the shared fallback (2, 4, 8 — starts 0.87…1.73) stay in
0.75…1.64; chains starting near the lower bound (1, 3, 5, 6, 7 — starts
0.02…0.06) are carried to 1.36…2.23. A run where every chain started from the
shared start = 1.0 could not produce the second group. Swarm-wide iter-0
spread is 1.477.

The defect was in the test: it hard-coded "chain 1 vs chain 8", assuming
chain_starts came out ordered. It does not — the LHS draws land in arbitrary
chain order, so the pair can be any two of the eight, and chains 1 and 8
happened to land 0.267 apart. The pair is now selected by actual extremity of
start (chains 3 and 8 here): |2.2312 − 1.6384| = 0.593. The > 0.3
threshold is unchanged.

Added sim/tests/if2_honours_per_chain_initial.rs — the exact, threshold-free
version of the same property at the level it lives: same base_params, seed and
config, differing only in EstimatedParam::initial, must give a different first
θ and a different MLE; identical .initial must reproduce bit-for-bit. That is
the real guard for the 2026-04-18 incident and no change to filter dynamics can
retune it.

On the run's 2/2 parameters have R̂ > 1.1 (max 5.07): that is a 2-iteration,
50-particle smoke fixture whose 8 chains start anywhere in [0.01, 5.0]; R̂ is
large on main too. It is not evidence of (b) — the chains have not converged
because they were never asked to.

profile_multi_stream::profile_family_root_sums_all_expanded_streams

Root cause: the R0 = 15 cell returned -inf, which serialises as
"best_loglik": null (JSON has no ±inf), so the point dropped out of the test's
collect_logliks — "expected 2 grid points, got [-1443.65…]".

The cell is not genuinely ruled out; it is under-sampled. The sweep straddles
the truth (true_params has R0 = 20), and a 100-particle bootstrap filter on
this 5-patch spatial SEIR loses the whole swarm that far off the optimum (cell
wall time ~8 ms). Measured:

--particles	R0 = 15 loglik
100	null (−inf)
500	−1876.13
2000	−1610.70

Two changes, kept separate:

1. The silent drop is now loud (fix(profile): …). A non-finite cell loglik
   emitted nothing — the leaf was still written and the run still printed "2
   cells written", so the point just vanished from the curve. It now warns on
   stderr, names the grid point, and names --particles as the lever. Verified
   it fires once at 100 particles and zero times at 500. stderr only; the
   recorded inputs JSON is byte-unchanged, so no CAS leaf re-keys.
2. The fixture gets enough particles to answer (100 → 500, the smallest
   measured value that scores every cell). Every assertion is untouched: still
   assert_eq!(len, 2), still finite-and-negative at both points, still the ≥3×
   family-vs-single magnitude check. The sibling test that asserts an
   argument-validation error stays at 100 — it never runs a filter.

Verification

sim  if2_honours_per_chain_initial    2 passed
cli  synthetic_fit_grid               8 passed
cli  profile_multi_stream             2 passed

Plus every IF2-touching cli binary (rg -l 'if2|IF2' rust/crates/cli/tests/,
27 binaries) swept green — profile_diagnostics, profile_grid_identity,
profile_priors, profile_greedy_order_invariant, profile_indexed_long_obs,
profile_pmmh, prior_posterior_invariants, calendar_fit_summary,
gh226_degenerate_fit, fit_predict_e2e, acceptance_fit_cas, cas_integration (33),
fit_sparse_holes, fit_survey_denominator, fit_experiment_management,
backend_provenance, compare_auto_prequential, contrasts_e2e,
fit_table_quantities, pfilter_indexed_long_obs, polio_afp_es_multicadence,
ode_dt_check, state_grad_gating, survey_top_k_pgas, survey_top_k_pmmh.