feat(codegen): analytic df/dp for Functional rate laws, emitted and gated off (#66)

[wshlavacek]

Closes #66. Stage 2 of #55, on top of #65 (obs[]/func[] plumbing). Blocks #67.

What this is

generate_sens_from_model returns None on the first non-Elementary reaction,
so one Functional rate law out of 689 puts a whole model on CVODES' internal
difference quotient. This derives the missing half.

Nothing changes at runtime. The != "elementary" → return None gate stays
shut: the new path is reachable only through a keyword-only functional= that
defaults to False and that no production caller sets. #67 flips it.

The derivative needed no emitter surgery

A Functional rate is stat_r · func_r(obs, p, t) · [af_r · ∏R_r], the species
factor is parameter-free, and an observable is a fixed linear combination of
species (∂obs/∂p = 0 at fixed y), so

df_i/dp = Σ_r  stat_r · netstoich_ir · (dfunc_r/dp) · af_r · prod(R)_r

is the same (derivative expression) × geometry product _emit_sens_rhs_body
already emits for a derived rate constant (#15/#41). The terms ride that existing
derived_terms channel with ∂func_r/∂p in place of ∂p_d/∂primary. What is
new is only the differentiation — ~90 lines, because _inline_functions,
_exprtk_to_sympy, sympy_to_c and _derived_param_jacobian_checked all apply
as-is, exactly as #55's reuse table predicted.

A derived parameter is deliberately not inlined into the rate law: it keeps
its own column as the direct partial (matching how the Elementary path emits
case iP for a _rateLaw{N}), and the chain rule to its primaries is added as
extra terms — which is also what keeps a nested derived parameter reachable.

The SIR exemplar, betaI() = beta*I with beta = 1/S0 # ConstantExpression:

case 0:   /* S0   — reachable ONLY through beta = 1/S0 */
    v = ((obs[1]) * (-1/pow(p[0], 2))) * y[0];
    dfdp_out[0] -= v;  dfdp_out[1] += v;
case 2:   /* beta — the direct partial */
    v = (obs[1]) * y[0];
    dfdp_out[0] -= v;  dfdp_out[1] += v;

Decline loudly

CVodeSensInit1 takes one callback for every column, so a decline is necessarily
model-wide — there is no per-reaction fallback to mix in. Every refusal names
what blocked it, and an empty term list stays a success (a genuinely zero
column) rather than being confused with a failure; that ambiguity is what #56 was
filed about.

The construct gate reuses _jacobian.unsupported_expr_construct rather than
re-spelling it. Two spellings drift, which is exactly how the #53/#56 hole
survived a release. The conditional class it rejects is #68's to lift, and only
behind the shared switch-time guard.

One decline the per-reaction loop cannot see coming — every law differentiates,
but the obs[]/func[] values they are written in need data->tfun_eval, which
CodegenSensUserData does not have (#65) — gets its own warning rather than a
quiet None.

Oracle: solver-free, corpus-wide

Central finite difference of the emitted bngsim_codegen_rhs against the
emitted ∂f/∂p, both through ctypes, yS = 0 so J·yS vanishes exactly. A
primary-parameter perturbation goes through set_param, which re-derives every
ConstantExpression parameter — so the #15/#41 chain rule is under test rather
than assumed.

Three well-separated steps, because cancellation noise falls as 1/h while
truncation grows as h². This is not belt-and-braces: a first pass at one step
reported its own noise as codegen bugs. Worst offender, an RHS with
‖f‖∞ = 3.4e7 where the analytic answer is exactly 1:

rel. step	1e-8	1e-7	1e-6	1e-5	1e-4	1e-3	1e-2
FD error	1.0e0	1.0e0	4.6e-2	1.4e-3	1.0e-3	7.2e-5	9.5e-7

Monotone 1/h decay — the textbook noise signature, not a wrong derivative.

Measured: 585 .net files under benchmarks/suites/ode_fullnet/nets

	models
all-Elementary (untouched)	400
Functional / MM	185
├ emitted — 5 315 ∂func/∂p terms	137
├ declined, with a reason	48
└ errors	0

14 127 (parameter, state) columns differenced. Worst disagreement 0.22 of the
per-component tolerance — every emitted model agrees to FD precision, with a
4.4× margin at the worst point.

Declines: 38 if() · 6 max() · 2 min() · 1 derived parameter whose own
flattened expression carries a min() · 1 Michaelis–Menten (out of scope; worth
exactly one extra model). The 137/48 split matches the 140 CLEAN / 45 blocked
partition #55's staging measured independently.

Derivation cost: 2.6 s for all 185 models, 0.31 s worst case, memoized on the
rate-law text since a rule-generated network reuses one expression across many
reactions.

The no-behavior-change claim is verified, not argued

generate_sens_rhs_c, generate_sens_from_model and the full combined source
are byte-identical to main on all 585 models, and the 189-model fallback set
is unchanged (396 models get a sens RHS before and after). test_codegen_functional_dfdp.py
pins the same claim per-model, including that an Elementary model's #15 derived
chain is untouched with the keyword set.

Not covered, and declined rather than approximated

Cross-compartment (per-species volume scaling) Functional reactions: their
∂f/∂p scatter would need the per-row (possibly live, #171) volume divide the
emitted switch has no form for. No .net in the corpus is one.

The .net emitter generate_sens_rhs_c is deliberately untouched, for #65's
reason: the .net parse carries neither the observable-entry nor the
table-function shapes _emit_observable_lines / _emit_function_lines consume,
which is also why generate_combined_c already sources its Jacobian from the
model. Michaelis–Menten stays out of scope — 11 reactions across 2 models, one
of which is Functional-blocked anyway, so it unblocks exactly one extra model.

For #67

bngsim_jac_vec is still Elementary-only, so a Functional model's
bngsim_codegen_sens_rhs is incomplete — the emitted C says so in its own
header, not only in the Python docstring:

/* GH #66: bngsim_dfdp below also covers Functional rate laws, but
   bngsim_jac_vec is still Elementary-only (GH #67) — so this
   bngsim_codegen_sens_rhs is INCOMPLETE and must not be installed as
   a model's sensitivity RHS. ... */

Testing

• 30 new tests in python/tests/test_codegen_functional_dfdp.py (gate-shut,
  byte-identity, emitted terms, seven decline classes, four FD-oracle models).
• Full suite: 2 694 passed, 89 skipped; pre-push hook green.
• ruff, ruff-format, mypy clean.