fix(steady-state): size the finite-difference probes relative to what they perturb (#76)

[wshlavacek]

Closes #76.

dY_ss/dp = -J⁻¹·(∂f/∂p) finite-differences either factor when no closed form is
available, and both probes sized their step as eps·max(|x|, 1) with
eps = sqrt(machine eps) = 1.49e-8. The floor is there to survive x == 0, but it
also overrides the relative step for everything below 1, and then the probe is no
longer small compared with what it perturbs.

The reported case

ode/before_bunching carries KD = 1e-9 with the forward rate constant derived
as kon/KD, so the probe is 1500% of the parameter — it drags kf from 1.0
to 0.063, and the quotient is a secant across a decade and a half of the rate
law's curvature.

max abs dY_ss/dKD
formula, shipped step	6.754783e+10
formula, relative step	1.074089e+12
truth (re-solved central difference, h/p = 1e-3 / 1e-4 / 1e-5)	1.074089e+12 / 1.074089e+12 / 1.074091e+12

The fix

• parameters — eps·|p|. The old absolute step survives only for a
  parameter that is exactly zero (no scale of its own) or subnormal (no relative
  step survives the addition, which would divide by zero). Parameters have no
  common unit, so |p| is the only scale the model offers.
• species — eps·max(|y_j|, max|y|): relative to the species, floored at the
  state's magnitude rather than at 1.0. Unlike parameters, every species is a
  concentration in one unit, so the state has a typical scale to probe a zero
  species against. The absolute 1.0 was wrong in both directions — a nanomolar
  model probed at 1 molar, and a molecule-count model (~1e6) probed at 1e-14 of
  its own state, which is cancellation noise, not a derivative.
• the state scale skips species a mask= excluded (steady_state(): no way to exclude write-only accumulator species from the convergence test #74): a write-only
  accumulator holds whatever integration left it at (7.5e8 on Barua 2013 against
  1e-10 for the settled species) and would otherwise set the probe for everyone
  else.
• both quotients divide by the realized step (x + h) - x, not the requested
  one.

All four sites take it: the FD Jacobian and FD ∂f/∂p in compute_ss_sensitivity,
and the state-chain and explicit-parameter halves of
compute_ss_output_sensitivity's ∂func/∂p fallback.

Blast radius (the issue asks for this before picking a fix)

The FD ∂f/∂p is reachable from Simulator.steady_state(sensitivity_params=…)
only when codegen emits no bngsim_codegen_sens_rhs — since #67/#89 that is
Michaelis-Menten and Functional laws carrying a condition or a non-smooth
builtin: 41 of the 585 ode_fullnet models. 15 of those also carry a
parameter the floor probes at more than 0.1% of its value (up to 1.27e30 for
ode/AVdyn6's epsilon = 1.18e-38). Driving find_steady_state without a
codegen artifact — how the repo's own tests reach the fallback, and how
before_bunching was found — exposes any of the 267 corpus models carrying
such a parameter. The species probe is reachable wherever the analytical Jacobian
is incomplete, and through jacobian="fd" everywhere.

Those 15 do not, today, yield a case the re-solve oracle can score: 11 never
reach a steady state at all (tol=1e-11, max_time=1e7), and on the rest the
small parameters do not move the root (truth ≡ 0) or the re-solve is not stable
in h. So the quantitative evidence below comes from the corpus A/B and from
closed-form models, not from those 15.

Corpus A/B: FD ∂f/∂p against the analytical one, at the same root

Every model solved twice — once with the compiled ∂f/∂p, once forced onto the
fallback — over its four smallest parameters plus the largest as a control.
1,201 columns over 358 models that converge and have an analytic column to be
scored against. Errors are normalized by the model's sensitivity scale, not
the column's own: dimer_phos2's dY/dNA peaks at 1.9e-41 against 1.1e-05 for
the same model's dY/dV, and normalizing per column makes that dust count as
much as the real answer.

Of the 850 columns where the step rule changes (|p| < 1) and the sensitivity
system is not itself degenerate (sens_jacobian_rcond > 1e-8):

	before	after
columns wrong by > 1e-3	103	57
columns wrong by > 1e-1	70	42
crossing 1e-3	—	51 fixed, 5 broken
crossing 1e-1	—	39 fixed, 11 broken

The 278 |p| >= 1 control columns move by at most 2.2e-9 (the realized-step
rounding), which is what confirms the change is confined to where the step
differs. The largest class of fixes is compartment volumes: 37 columns over 35
models — Vcyt / V / Vecf at 1e-12, probed at 15,000x themselves — were
100% wrong and now land within 1e-5.

The regressions are real and I am not hiding them. Where a parameter's own
term is a tiny fraction of the RHS component it sits in, a probe relative to the
parameter moves f by less than its roundoff, where the old wide step did not.
Worst case FceRI_viz's kp1 = 1.7e-6, 4.2e-07 → 5.0e-01. Fixing that needs a
step chosen from the measured response rather than from |p| alone — filed as
#123 rather than folded in here, because it is a different algorithm
(probe, measure the RHS's own jitter, widen per component) and wants its own
validation.

Corpus A/B: the species probe against the analytical Jacobian

1,066 model-states (every corpus model with a complete analytical Jacobian, at
its initial condition and at a generic all-species-populated state), scored on
‖J_fd − J_an‖∞ / ‖J_an‖∞:

rule	median	> 1e-6	> 1e-3	> 1e-1	max
`eps·max(	y_j	, 1)` (shipped)	1.5e-08	157	28
`eps·max(	y_j	, max	y	)` (this PR)	1.3e-09
`eps·	y_j	` (no floor)	6.9e-09	29	6

Head to head at 10x, the shipped rule is beaten 511 to 54; the floor-free
relative step is beaten 291 to 13 — the floor is what keeps a species far below
the state's scale out of the cancellation noise.

Tests

python/tests/test_codegen_steady_state_fd_step.py, 7 tests, 6 fail on main
(the seventh is the control asserting the well-scaled parameters did not move).
Three new fixtures, each with its closed-form steady state and closed-form
gradient derived in the .net header, so the assertions owe the solver nothing:

• reciprocal_rate_const.net — A <-> B with kf = kon/KD, KD = 1e-9
  (before_bunching's shape, and the reason the RHS is non-linear in KD — a
  bare mass-action constant is linear, so any step is exact for it).
  dA*/dKD = 2.5e8; the shipped step returns 1.572e7, off by 93.7%. Its
  perKD() = A_tot/KD function carries the same defect into
  compute_ss_output_sensitivity: -1.572e16 against -2.5e17 exact.
• nanomolar_dimer.net — A + A <-> B at A* = 1e-9 under jacobian="fd".
  Homodimerization is quadratic in A, so the species step does not cancel out of
  the Jacobian column; the shipped floor leaves a spurious -kdim·h worth 7.5x
  the real term and the gradient comes back 78.8% low.
• nanomolar_dimer_sink.net — the same model plus a write-only accumulator at
  1e3, masked out. Pins that a species with no steady value does not set the
  probe scale for the species that have one.

The file is added to mir.yml's run list (it matches the existing
test_codegen*.py path filter, but that list is explicit — see the note in that
workflow).

Full python/tests suite green locally against the patched core: 2982 passed, 114 skipped.