feat(codegen): analytic sensitivity RHS for Functional rate laws (#67) - #84
Conversation
A single Functional reaction put a whole model on CVODES' internal difference quotient — `CVodeSensInit1` installs one callback for every sensitivity column, and `generate_sens_rhs_c` declined on the first non-Elementary rate law. On the 585-model corpus that was 187 models (32%), blocked by 1.1% of reactions. #66 supplied the analytic df/dp behind a keyword. This supplies the other half of `ySdot = J*yS + df/dp` and turns the keyword on. J*yS is not a second derivation: it is the same per-species chain rule and per-observable product rule `generate_jacobian_from_model` already emits. Both callers now share one reconstruction (`_functional_jacobian_groups`) and differ only in where a contribution lands. The sensitivity RHS passes a scatter that fuses the matvec — `Jv_out[i] += coeff*dj*v[j]` instead of `jac[j*n+i] += coeff*dj` — so no n*n scratch buffer is formed inside the callback (267 KiB on the widest Functional corpus model, memset per column per step), the work is O(nnz) not O(n^2), and CodegenSensUserData does not widen. The analytical Jacobian's emitted C is byte-identical on all 585 corpus models, as is every Elementary model's sensitivity RHS. `generate_combined_c` falls back to the model-based emitter when the .net one declines: that is the path a .net-loaded model actually takes, so without it #67 would reach only the SBML/Antimony entry points. Corpus: 535/585 models get the analytic sens RHS, up from 398 (+137), 0 lost, 137/137 compile and export bngsim_codegen_sens_rhs. The 50 that still decline are Michaelis-Menten and the Functional laws carrying a condition (#68) or a non-smooth builtin; every decline warns and names what blocked it. Sampled end to end on 46 gained models: 46/46 agree with the difference quotient, median relative disagreement 6.1e-08. Cost per analytic column at 16 parameters is 0.59-1.04 plain-solve equivalents against the difference quotient's 9.0-36.7; DQ needed a median 295x more steps, and on 2 of 6 benchmarked models never converged at all — it exhausted mxstep and returned a gradient of exactly zero. BNGSIM_NO_FUNCTIONAL_SENS_RHS=1 restores the previous behaviour; it is part of the .net cache key so a hatched run cannot collide with a .so compiled without it. steady_state.cpp reads the same symbol at yS = 0, so a Functional model's ss.sens_dfdp_source moves from "finite-difference" to "codegen", retiring the sqrt(eps) step floor #76 describes for those models.
…ackend Both files already *trigger* the MIR job (the paths filter is python/tests/test_codegen*.py) but neither was in the list the job runs, so #66's and #67's tests fired the matrix without ever executing in it. The pre-push hook was their only gate. Running them under BNGSIM_CODEGEN_JIT=mir also covers something no local run does: c2mir compiling the newly-emitted Functional df/dp and J*v, rather than only cc. The one test that reads the installed artifact now accepts either shape (.so symbol or JIT C source), since the JIT path produces no .so and asserting on it alone would make the JIT legs a false failure.
|
Added a second commit ( Both files already triggered that job — its paths filter is One test needed adjusting for that: |
…he MIR backend" Adding these two files to the MIR job turned all four legs red, and the cause is not #67: c2mir rejects `bngsim_codegen_output_sens` (GH #198), which this branch emits byte-identically to main. Verified on a probe branch off clean main carrying only a Functional model run with sensitivity_params under the JIT — it fails there too, with the same "syntax error on double (expected '<statement>')". Nothing in the MIR job's current list reaches that combination: the one Functional test there pre-sets `_codegen_so_path`, which makes the Simulator skip auto-codegen entirely, so the JIT never sees the module. So the mir.yml line belongs with the c2mir fix, not with #67. Filed separately; keeping the JIT-aware assertion in test_the_analytic_rhs_is_actually_the_one_installed, which is correct either way.
|
Reverted the mir.yml commit ( c2mir rejects Nothing in the MIR job's list reached that combination before: the one Functional Kept from that commit: the JIT-aware assertion in Worth flagging for review: #85 grows in scope once this lands. Today it bites |
…ing on (#87) test_codegen_functional_dfdp.py (#66), test_codegen_functional_sens_rhs.py (#67) and test_codegen_switch_condition_sens.py (#68) all match the python/tests/test_codegen*.py paths filter in mir.yml, so editing them fires the job — but none was in its pytest invocation, so the matrix went green having never run them. The pre-push hook was their only gate, and --no-verify removes even that. #66's file was reverted out of #84 and #68's left out of #86 for the same reason, which is #85; this lands all three together so the gap retires in one place. Measured under a local -DBNGSIM_ENABLE_MIR=ON -DBNGSIM_ENABLE_KLU=OFF build, which is the mir.yml configuration exactly (HAS_MIR True, HAS_KLU False): dfdp under mir 29 passed | under cc 29 passed sens_rhs under mir 7 FAILED | under cc 27 passed switch_condition under mir 3 FAILED | under cc 27 passed Nothing is cc-assertion specific — dfdp never builds a Simulator, so it is backend-agnostic. All 10 failures are #85, open and pre-existing: c2mir rejects the #198 bngsim_codegen_output_sens block that cc accepts, so a Functional model constructed with sensitivity_params cannot compile under the JIT. The whole of mir.yml's list is now 229 passed / 10 xfailed. blocked_on_mir_jit quarantines exactly those 10. In sens_rhs it is per-test: test_steady_state_dfdp_stops_differencing is deliberately NOT marked, because it reaches the sensitivity codegen without ever setting _want_output_sens and so passes under the JIT. In switch_condition it is on the class, where every test goes through _run_sens. xfail(strict=True), not skipif, so the quarantine retires itself: verified that with the marker forced on while the tests pass, all 7 report XPASS and the run goes red. A skipif would have sat there until someone remembered. raises= RuntimeError keeps it from absorbing an unrelated failure (MirJit raises RuntimeError; SimulationError subclasses it), so an assertion regression behind the quarantine stays red rather than reading as expected.
Closes #67. Stage 3 of #55; follows #65 (PR #80) and #66 (PR #83).
The problem
A single Functional reaction put a whole model on CVODES' internal difference
quotient.
CVodeSensInit1installs one callback for every sensitivitycolumn, and
generate_sens_rhs_creturnedNoneon the first non-Elementaryrate law — so 1.1% of reactions blocked 187 of 585 corpus models (32%).
#66 derived the analytic
∂f/∂pand left it gated off, because the other half ofySdot = J·yS + ∂f/∂pwas still Elementary-only. This supplies that half andopens the gate.
What it does
J·ySis not a second derivation. It is the same per-species chain rule andper-observable product rule
generate_jacobian_from_modelalready emits andvalidates. Both callers now share one reconstruction,
_functional_jacobian_groups, and differ only in where a contribution lands:Both indices are known at emit time, so the matvec fuses into the scatter. No
n×nscratch buffer inside the CVODES callback (267 KiB onatg_model_v2, andbngsim_codegen_jacmemsets it per column per step), O(nnz) instead of O(n²),and no
CodegenSensUserDatawidening — the alternative §1.3 of the plan costedout.
bngsim_jac_vectakest/obs/funconly when its body reads them, so anElementary model's signature and whole emitted source are unchanged.
generate_combined_cfalls back to the model-based emitter when the .net onedeclines. The .net emitter reads rate laws as text and has no expression to
differentiate; the built model does. This is the path a
.net-loaded modelactually takes, so without the hook #67 would have reached only the
SBML/Antimony entry points.
Validation
Byte-identity.
generate_jacobian_from_modelandgenerate_sens_from_modelemit byte-identical C on all 585 corpus models, before vs. after the extraction
(both emitters loaded in one process and string-compared).
Corpus.
The 50 still declining are Michaelis–Menten and the Functional laws carrying a
condition (#68 owns those) or a non-smooth builtin. Every decline warns and names
what blocked it.
Correctness. Four independent oracles, none circular:
bngsim_codegen_rhsin each species, againstJ·ySat unityS— solver-free, through ctypes, over both reconstructionbranches (per-observable and per-species).
bngsim_jac_vecagainst the densebngsim_codegen_jacfrom the same .so,times
vin Python. This is what tests the matvec fusion specifically; an FDoracle cannot separate a bad fusion from a bad derivative.
p ± h— the only non-circular end-to-end groundtruth, since the difference quotient is the thing under suspicion.
above 100 reactions, plus a random 24): 46/46 agree, median relative
disagreement 6.1e-08, worst 6.0e-05.
Performance, warm, cost per analytic column in plain-solve equivalents
(
(t_sens − t_plain) / (k · t_plain)):#55's claimed 0.75–0.95 band lands at realistic
k; atk = 4the per-columnfixed cost still shows on sub-millisecond solves. The DQ needed a median 295×
more steps across the 46 sampled models, up to 177 060×.
And on 2 of the 6 benchmarked models the difference quotient never converged
at all — it exhausted
mxstepand returned a gradient of exactly zero, whichreads downstream as a converged answer. The SIR fixture in the new test file is
one: 30 920 562 DQ steps returning zeros, versus 1 777 analytic steps landing on
the re-solved trajectory to 7 digits. Verified identical on
main, so this is apre-existing silent-wrong-answer that #67 fixes rather than anything introduced
here.
Downstream
steady_state.cppreads the samebngsim_codegen_sens_rhsatyS = 0for the∂f/∂pfactor ofdY_ss/dp = -J⁻¹·(∂f/∂p), so a Functional model'sss.sens_dfdp_sourcemoves"finite-difference"→"codegen"— retiring the√eps step floor #76 describes, for those models. Checked against a re-solved
steady state (9 significant digits) and pinned by a test.
Notes for review
BNGSIM_NO_FUNCTIONAL_SENS_RHS=1restores the previous behaviour for an A/B.It is part of the .net codegen cache key and the in-process memo key, so a
hatched run cannot collide with a
.socompiled without it.question) —
CVodeSensInit1is all-or-nothing, so there is nothing finer tokey on.
unsupported_expr_constructthrough #55b — analytic df/dp for Functional rate laws (emitted, gated off) #66's_functional_dfdp_terms, so #55d — extend the switch-condition guard so condition-bearing Functional laws can use the analytic sens RHS #68 hasone place to lift, behind the shared switch-time predicate.
src/steady_state.cppis a comment-only change (its "every Functional/MMmodel" claim is now wrong).