Steady-state FD ∂f/∂p uses an absolute step floor, so a small rate constant is probed at 1500% (15.9x error on before_bunching)

[wshlavacek]

The finite-difference ∂f/∂p fallback in compute_ss_sensitivity
(src/steady_state.cpp:1077)
sizes its step with an absolute floor of 1:

double h = eps * std::max(std::abs(pval), 1.0);   // eps = sqrt(machine eps) = 1.49e-8

That floor is there to survive pval == 0, but it also overrides the relative
step for every parameter smaller than 1. For a rate constant of 1e-9 the probe
is 1.49e-8 — 1500% of the parameter — so the difference quotient is not a
derivative at that point at all, it is a secant across a decade and a half of the
rate law's curvature.

This only bites the issue #55 fallback path (sens_dfdp_source == "finite-difference", i.e. models whose rate laws are not all Elementary so the
codegen sensitivity RHS was not emitted). With an analytical ∂f/∂p the step is
irrelevant.

Reproducer

before_bunching has KD = 1e-9. Driving find_steady_state with no codegen
artifact, which is how the repo's own tests reach the FD fallback
(python/tests/test_steady_state_codegen.py::_core_sensitivity):

import numpy as np, bngsim
from bngsim._bngsim_core import SteadyStateOptions, find_steady_state, NetworkModel

P = "benchmarks/suites/ode_fullnet/nets/original__bngl_models__my_models__ode__before_bunching.bngl.net"

core = NetworkModel.from_net(P)
o = SteadyStateOptions()
o.tol, o.method, o.jacobian = 1e-12, "integration", "auto"
o.sensitivity_params = ["KD"]
r = find_steady_state(core, o)
assert r.sens_dfdp_source == "finite-difference"
s = np.array(r.sensitivity_data).reshape(len(r.species_names), -1)
print("formula:", np.abs(s).max())

# truth: central difference of the steady state itself
def root(v):
    m = bngsim.Model.from_net(P); m.set_param("KD", v); m.reset()
    return np.asarray(bngsim.Simulator(m, method="ode").steady_state(tol=1e-12).concentrations)

p0 = bngsim.Model.from_net(P).get_param("KD")
for frac in (1e-3, 1e-4, 1e-5):
    h = frac * abs(p0)
    print(f"truth (h/p={frac:.0e}):", np.abs((root(p0 + h) - root(p0 - h)) / (2 * h)).max())

	max abs dY_ss/dKD
formula, FD ∂f/∂p with the shipped step	6.754783e+10
truth, h/p = 1e-3	1.074089e+12
truth, h/p = 1e-4	1.074089e+12
truth, h/p = 1e-5	1.074091e+12

15.9x low. The reference is stable to 6 significant figures across three step
sizes, so this is the formula being wrong, not the reference being noisy. The
sens_jacobian_source here is analytical and ss.sens_jacobian_rcond is
1.18e-1 — the Jacobian factor is fine and the conditioning diagnostic is
healthy, so nothing on the result flags it. Recomputing the same ∂f/∂p with a
purely relative step eps*|p| brings the formula to within 4e-8 relative of the
reference.

Scope

262 of the 585 ode_fullnet models carry at least one parameter small enough
that the floor forces a probe larger than 0.1% of its value. The tail is extreme —
ode/AVdyn6 has epsilon = 1.2e-38, where the probe is ~1e32 times the
parameter, and ode/temp has kB = 1.4e-23. How many of those also land on the
FD ∂f/∂p path I have not measured; that intersection is the real blast radius
and is worth quantifying before picking a fix.

Suggested fix

Scale relatively when the parameter is nonzero, and keep an absolute floor only
for the genuinely-zero case:

double h = eps * (pval != 0.0 ? std::abs(pval) : 1.0);

A typical value per parameter would be better still, but bngsim has no such
notion and |p| itself is the standard stand-in. Note the same
eps * std::max(std::abs(...), 1.0) shape is used for the state perturbation
in the FD Jacobian a few lines above (h = eps * std::max(std::abs(y_ss[j]), 1.0)); there the floor is defensible, since a species at zero concentration is
common and a relative step would be exactly zero. It is worth a look regardless —
a model in units where all concentrations are ~1e-9 has the same problem.

Provenance

Found while investigating the issue #63 follow-up (whether
ss.sens_jacobian_rcond can gate a refusal). It first showed up as a false
positive in that analysis: before_bunching looked like a model returning a wrong
gradient despite excellent conditioning, and the cause turned out to be the step
size rather than the linear algebra.