Verify the assembled gradient against finite differences on every real gntr problem, not just fixtures

[wshlavacek]

Summary

The gradient path's FD oracles all run on purpose-built fixtures — small SBML/BNGL models written
to exercise one seam each. That suite is green and has been for the whole #385 epic.

The first time an FD check was pointed at a real model from the benchmark corpus, it found a bug
in under an hour: #534, where a free parameter bound by id that reached the trajectory only by seeding
another entity assembled to exactly 0 against a central difference of -10.4, stable at every step
size. No fixture had that shape, because the shape comes from how COPASI exports SBML, not from
anything a hand-written test would think to write.

Eleven slugs in the subset-I corpus run job_type = gntr. One of them has ever had its gradient
compared against finite differences. This issue is to do the other ten.

Why this is worth real time

• A wrong gradient column does not fail loudly. The objective stays correct, score.py still reports
  a valid OG, and the fit simply searches with less information. Laske solved at OG = 0.104 with a
  dead column and at -1e-06 without one, using an identical budget — the only visible symptom was
  a fit that stopped a bit short.
• The corpus is a far better generator of shapes than we are. It carries COPASI aliasing, mixed
  natural-log and linear noise, prediction-dependent σ, per-measurement σ sidecars, replicate
  observableParameters, t=0-only experiments, steady-state measurements, and pre-equilibration.
• Each check is cheap: ~2·k simulations plus one, i.e. a couple of minutes for most slugs.

The work

Run tools/fd_check.py (in the corpus, wshlavacek/BNGL-Models) at each slug's PEtab nominal point
and record the worst relative error per slug. The gntr slugs not yet checked:

• Armistead_CellDeathDis2024 (k=14) — prediction-dependent σ
• Blasi_CellSystems2016 (k=9) — steady-state (time = inf) measurements, lnnormal
• Boehm_JProteomeRes2014 (k=9) — ratio-formula observables
• Bruno_JExpBot2016 (k=13) — all 13 parameters routed through condition parameter references
• Crauste_CellSystems2017 (k=12) — fixed per-point σ
• Fiedler_BMCSystBiol2016 (k=22) — replicate observableParameters
• Rahman_MBS2016 (k=9)
• Raia_CancerResearch2011 (k=39) — prediction-dependent σ, largest k
• SalazarCavazos_MBoC2020 (k=6)
• Sneyd_PNAS2002 (k=15)
• Zhao_QuantBiol2020 (k=28) — per-measurement σ sidecars (ADR-0083)

(Laske_PLOSComputBiol2019 is done: worst 2.5e-03 across 13 columns, after #534.)

Anything that turns up gets its own issue and fix here; the corpus side is tracked at
wshlavacek/BNGL-Models#38.

Doing it correctly

tools/README.md in the corpus records two failure modes that both produce a plausible wrong
answer rather than an error, and I hit both:

• Score the Result before assembling. Otherwise the objective's _pset_values is unbound and its
  measurement layer unmaterialized, and every free parameter living in an observable formula or a
  noise scale reads back as exactly 0.0 — indistinguishable from a real defect.
• Sweep the step size. On a stiff model h = 1e-7 is pure roundoff. A real defect does not move
  with h; FD noise does.

Longer term

If this sweep finds more than one further bug, a standing FD gate over a handful of real models —
rather than only fixtures — is probably worth adding to the recovery tier.

[wshlavacek]

Done — with one correction to the issue's premise and one product bug.

The corpus has 18 gntr slugs, not 11

The six the checklist misses are Bertozzi_PNAS2020, Brannmark_JBC2010, Giordano_Nature2020,
Oliveira_NatCommun2021, Schwen_PONE2014, Weber_BMC2015. All eighteen are checked below.

A gradient check is only a test where the gradient is large

Running fd_check.py at each slug's PEtab nominal point, as this issue specifies, flags six
slugs red. Every one of those flags is an artefact, and the reason is structural: twelve of the
eighteen slugs have OG_nominal ≈ 0 — their nominal point is the optimum, where the true gradient
vanishes and "assembled vs central difference" compares two quantities that are both noise. The
correlation is exact: the slugs that came back clean at the nominal point were precisely the ones
far from the optimum (Blasi at OG 448, Zhao at 276, Laske at 40, Raia at 0.78).

Fiedler_BMCSystBiol2016 sits at OG_nominal = -0.0022, and its real defect was
indistinguishable from those eleven artefacts at the nominal point. Re-evaluating every slug at a
point displaced a few percent of each box width separated them in a single pass.

Two further harness traps, both of which manufacture a plausible-looking dead column:

• Box bounds. FreeParameter.set_value clamps an out-of-box value, so a parameter sitting on
  a bound has f(u+h) == f(u): its central difference is a half-step, or exactly zero when both
  sides clamp. Four of the five remaining red rows were this and nothing else.
• Sweep h upward, not just downward. On these models the small-h end is the roundoff end. A
  column drifting monotonically away from the assembled value as h shrinks is converging at the
  other end — Fiedler's k2 reads 216.07 / 212.04 / 210.02 at h = 1e-4 / 3e-4 / 1e-3 against an
  assembled 210.95, and lands on it at 1e-2.

All three are now written into tools/README.md beside the two it already records, and fd_check.py
grew a --disp flag and keeps the evaluation point 8h clear of both bounds.

The bug: Fiedler_BMCSystBiol2016, seven of twenty-two columns

Fiedler initialises all six species from initialAssignments that are the closed-form steady state
of the network — long rational expressions over K_1, k2 … k11. Those same symbols are the rate
constants in every kinetic law. ADR-0096 had promoted "an id that only seeds species initial
conditions gets no parameter axis of its own" to the general rule, on the premise that such an id is
absent from the ODE right-hand side. Here it is not, so every one of those parameters assembled from
its six IC seed terms alone and never entered sensitivity_params at all:

param     assembled   central diff        after
K_1         +95.022       -339.83       -341.04
k2          -95.027       +212.04       +210.95
k3         +162.37        -339.74       -341.13
k4         -162.37        +270.79       +269.52
k5         +168.62        +597.44       +599.36
k6         -168.62        -657.72       -658.66
k10         +95.027       -174.85       -174.05

Several are sign-reversed — a gradient fit was being steered uphill on them. This is the common
case in the PEtab benchmark collection, not an exotic shape.

Fixed in ADR-0097: the axis is dropped only when the model states the right-hand side never reads
the id (ode_rhs_symbols(), on both backends), and that answer is a veto on dropping, never a
new reason to drop — so the change can only ever add a column that was missing, never remove one
that is there. (Replacing the rule outright was tried first and rejected: it makes every parameter's
axis contingent on that set being complete, and an initialAssignment on a compartment size is
exactly the symbol that would go missing.)

Results: all 18 slugs, at a displaced, bounds-clear point, h ∈ {1e-2, 3e-3, 1e-3}

worst |Δ|/‖∇‖ is the largest disagreement between the assembled column and the finest-step central
difference, relative to the gradient's infinity norm — i.e. how much the descent direction is
distorted.

slug	k	OG_nominal	‖∇‖∞	worst |Δ|/‖∇‖
Armistead_CellDeathDis2024	14	5.8e-06	4.81e+03	4.2e-05
Bertozzi_PNAS2020	8	5.09e-06	2.68e+03	7.5e-05
Blasi_CellSystems2016	9	448	9.48e+03	1.1e-05
Boehm_JProteomeRes2014	9	3.16e-05	352	3.7e-07
Brannmark_JBC2010	22	1.53e+03	1.90e+03	2.4e-04
Bruno_JExpBot2016	13	3.23e-06	3.16e+03	3.2e-06
Crauste_CellSystems2017	12	0.509	173	5.8e-05
Fiedler_BMCSystBiol2016	22	-0.00223	1.20e+03	4.4e-03 †
Giordano_Nature2020	50	3.78e+03	106	1.3e-04
Laske_PLOSComputBiol2019	13	39.9	271	7.0e-06
Oliveira_NatCommun2021	12	9.62e+06	6.74e+08	2.2e-04
Rahman_MBS2016	9	3.9e-06	8.60e+03	1.2e-05
Raia_CancerResearch2011	39	0.78	2.27e+04	7.2e-06
SalazarCavazos_MBoC2020	6	0.326	1.28e+03	1.6e-05
Schwen_PONE2014	30	-8.42	6.08e+06	1.2e-04 ‡
Sneyd_PNAS2002	15	0.000582	546	5.5e-05
Weber_BMC2015	36	1.37e+04	4.95e+03	2.0e-06
Zhao_QuantBiol2020	28	276	1.69e+18	2.4e-04

† Fiedler's residual is tau2, whose own central difference spans −335 to −348 across h with no
convergence; the assembled −338 sits inside that band.
‡ Schwen's IR_obs_std disagrees only when the bound-clearing margin pushes it toward its box edge,
where the σ-NLL is violently nonlinear; at any interior point it agrees to 7e-09.

One thing that did not reproduce

Raia_CancerResearch2011's init_Rec_i read exactly 2× its central difference on a single run
(−45475 against −22738, FD stable to five digits). It has not recurred in five subsequent runs at
that same point, including two with byte-identical routing, and Raia is clean above. A
non-reproducible factor-of-two on a sensitivity column is worth naming rather than burying —
lanl/bngsim#117 (the printer-thread race behind #525) is the obvious suspect. I have not filed
anything, since I could not reproduce it.

Longer term

The issue asks whether a standing FD gate over real models is worth adding to the recovery tier if
this sweep finds more than one bug. It found one product bug and three ways for the check itself to
lie. On that evidence the higher-value gate is not "run FD on more models" but "run FD on a model
with a steady-state initialAssignment" — which is now a fixture-level FD oracle in
test_gradient_sbml.py, running in ordinary CI, and it fails loudly (−5.20 against +1343.88) under
the old rule.

[wshlavacek]

The unreproduced Raia 2x reading noted above is now filed as #537, with the ruled-out causes and the five failed reproduction attempts.