feat(codegen): analytic sensitivity RHS for clock-conditional Functional rate laws (#68) - #86
Merged
Merged
Conversation
…nal rate laws (#68) #67 declined every Functional rate law carrying an `if()`. Some of those declines were necessary and some were not, and telling them apart is all this stage does. if(t >= sigma, beta, 0) admitted. d(f)/d(sigma) is a clean 0 in both branch interiors, and that IS the answer there - #48 supplies the rest as the crossing jump s+ = s- + (f- - f+)*dt*/dp. if(I >= thresh, beta, 0) refused. The same clean 0 comes back, but the crossing moves with the trajectory and nothing supplies that term. The two are indistinguishable downstream: sympy differentiates the Piecewise to 0 with no Dirac delta, so `_is_emittable` never rejects it, the C compiles, and the gradient is silently zero. No oracle separates them either - the emitted derivative and a naive finite difference agree with each other and both miss the crossing. So the gate is a token pre-scan, and what is tested is the gate. One recognizer, two callers. `_clock_threshold_split` is lifted out of `compute_switch_time_sens` unchanged and reached through one `switch_condition_scope(core)` (one RHS clock probe, one parameter-table pass), so the gate cannot admit a condition the detector would not compensate. `unsupported_expr_construct` gains `allow_conditions=`, with the conditional class marked in the same table it already spells - the one-line relaxation #66 left room for, not a third predicate. abs/min/max/floor/ceil/round are kinks in the state with no crossing time and nothing to compensate them, so they stay rejected for every caller. The shipped rule is stricter than the issue's. #68 asks that every *requested* sensitivity parameter's occurrence in a condition be a clock threshold, which admits `if(Xtot > 0, k*A, 0)` - no parameter in the condition at all. But that crossing still moves with every parameter through the trajectory: it is the rate-law twin of the state-dependent event trigger `NetworkModel::event_sensitivity_unsupported_reason` refuses for #52, and two guards answering one question opposite ways is the drift #68 was filed to prevent. So an atom is admissible only if it is a recognized clock threshold, or names no symbol at all (`0>0`, which real corpus models carry). 8 models instead of 14. Measured before deciding: on a fixture with closed-form truth (dY/dk = -30.807) CVODES' internal difference quotient returns exactly 0 for a state-dependent crossing, the same as the analytic path would. So admitting that class was never a regression - it was extending the analytic path into a class we know is incomplete, for nothing. That limitation is pre-existing on both paths and unchanged here. Corpus: 535 -> 543 models with an analytic sensitivity RHS (+8, 0 lost), and the emitted C is byte-identical to main on the other 577 across both emitter paths. The 8 are exactly the counter-clock switch-time family #48 exists for - EID x2, m15, person_v1, CaM x2, m_vs_p_v1, diversity - and all 8 compile. Both halves carry the condition: d(f)/d(beta) emits `obs[1]*((p[4] <= obs[3]) ? (1.0) : (0.0))` and the fused J*v emits `d0 = ((p[4] <= obs[3]) ? (p[2]) : (0.0))`. sigma gets no d(f)/d(p) column at all, which is correct - the jump is #48's to apply, and an in-branch term would double-count it. End to end on a counter-clock switch model, analytic against the difference quotient: 2.4e-07 relative on every column including sigma, which is entirely the #48 jump on both paths. A trajectory finite difference cannot referee this - across a moving discontinuity it is only O(h) accurate and both paths miss it by the same 1.5e-03 - so the agreement number is the claim, and the smooth part is already covered by #66's and #67's oracles. Also fixes a latent ordering bug this stage exposed. `differentiate_rate_law` ordered its {observable: d(rate)/d(obs)} result by iterating a *set*, so a rate law reading two or more observables ordered every `double d0`/`d1` temporary and the `v[j]` column each is scattered into by PYTHONHASHSEED. Latent since #76: the saturable twin (`_saturable_jacobian.differentiate_rate_law_native`) already sorted for exactly this reason and the sympy branch never did, and the single-observable laws that dominate the corpus hid it. Four of the eight models admitted here read >= 2 observables, so their emitted C re-hashed on every process and the content-addressed .so could never hit. The new PYTHONHASHSEED test needs `sin()` in its fixture, or the already-sorted native path takes the law and it passes with the bug in place. `_switch_sensitivity` joins `_CODEGEN_SOURCE_MODULES`: the gate lives there now, so that module determines which models get a sensitivity RHS, and #51's digest did not cover it. _CODEGEN_VERSION 23 -> 24. `_switch_sensitivity.py` is added to mir.yml's paths filter for the same reason. The test file is deliberately NOT added to that job's pytest list: its end-to-end cases construct a Simulator with sensitivity_params on a Functional model, which is the combination #67 found c2mir rejects (GH #198, `bngsim_codegen_output_sens`), and #67 reverted its own mir.yml line over it. That line belongs with the c2mir fix. What #68 asks for and this does not do: fold `event_sensitivity_unsupported_reason` into the shared predicate. It answers the same question but on ExprTk *addresses*, and #52's fix was precisely that observable-vs-species name binding makes text analysis wrong - `v > 30` binds to the observable total, not the concentration - so re-expressing it against the Python text predicate would reintroduce the bug #52 closed. The two also carry deliberately different policies: the event path refuses `time >= p2` because event-time sensitivity is unimplemented, while the rate-law path supports it via #48. Cross references are in both directions in the source. Closes #68.
wshlavacek
added a commit
that referenced
this pull request
Jul 28, 2026
…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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #68. Last stage of #55, following #65 / #66 / #67.
#67 declined every Functional rate law carrying an
if(). Some of those declines were necessary and some were not, and telling them apart is all this stage does.if(t >= sigma, beta, 0)∂f/∂sigmais a clean0in both branch interiors, and that is the answer there — #48 supplies the rest as the crossing jumps⁺ = s⁻ + (f⁻−f⁺)·∂t*/∂sigmaif(I >= thresh, beta, 0)0comes back, but the crossing moves with the trajectory and nothing supplies that termThe two are indistinguishable downstream.
sympy.diffof thePiecewisereturns0with no Dirac delta, so_is_emittablenever rejects it, the C compiles, the solver converges, and the gradient is silently zero. There is no oracle that separates them either: the emitted derivative and a naive finite difference agree with each other and both miss the crossing. So the gate is a token pre-scan, and what is under test is the gate.One recognizer, two callers
_clock_threshold_splitis lifted out ofcompute_switch_time_sensunchanged and reached through oneswitch_condition_scope(core)— one RHS clock probe, one parameter-table pass — so the gate cannot admit a condition the detector would not compensate.unsupported_expr_constructgainsallow_conditions=, with the conditional class marked in the table it already spells: the one-line relaxation #66 left room for, not a third predicate.abs/min/max/floor/ceil/roundare kinks in the state, with no crossing time and nothing to compensate them, so they stay rejected for every caller.TestTheGateAndTheDetectorAgreeis what actually holds the two together. It asserts behavioural agreement — gate admits ⟺ detector emits a compensating record — rather than "they share a function", because a shared function can be refactored apart silently and behavioural agreement is the property #56 lost when #53 fixed one copy of a predicate and not the other.The shipped rule is stricter than the issue's
#68 asks that every requested sensitivity parameter's occurrence in a condition be a recognized clock threshold. That admits
if(Xtot > 0, k*A, 0)— no parameter in the condition at all — but that crossing still moves with every parameter through the trajectory. It is the rate-law twin of the state-dependent event triggerNetworkModel::event_sensitivity_unsupported_reasonrefuses for #52, and two guards answering one question opposite ways is exactly the drift #68 was filed to prevent.So an atom is admissible only if it is a recognized clock threshold, or names no symbol at all (
0>0, which real corpus models carry). That is 8 models instead of 14.Measured before deciding, not after. On a fixture with closed-form truth (
dY/dk = −30.807), CVODES' internal difference quotient returns exactly 0 for a state-dependent crossing — the same as the analytic path would. So admitting that class was never a regression; it was extending the analytic path into a class we know is incomplete, for nothing. That limitation is pre-existing on both paths and unchanged here.Corpus
535 → 543 models with an analytic sensitivity RHS (+8, 0 lost), and the emitted C is byte-identical to
mainon the other 577 across both emitter paths. All 8 compile. They are exactly the counter-clock switch-time family #48 exists for:EID_v2_current_code_params·EID_vs_RBM_18Feb2022·m15·person_v1·CaM_Ca_interaction_v1_1·CaM_KRas4B_PI3K_v5·m_vs_p_v1·diversityBoth halves carry the condition —
∂f/∂betaemitsobs[1]*((p[4] <= obs[3]) ? (1.0) : (0.0)), and the fusedJ·vemitsd0 = ((p[4] <= obs[3]) ? (p[2]) : (0.0)).sigmagets no∂f/∂pcolumn at all, which is correct: the jump is #48's to apply and an in-branch term would double-count it.Declines name the parameter and the condition, e.g.
End to end
Analytic against the difference quotient on a counter-clock switch model: 2.4e-07 relative on every column, including
sigma, which is entirely the #48 jump on both paths.A trajectory finite difference cannot referee this one — across a moving discontinuity it is only O(h) accurate, and both paths miss it by the same 1.5e-03, so it separates nothing. The agreement number is the claim; the smooth part is already covered by #66's and #67's oracles.
Also: a latent ordering bug this stage exposed
differentiate_rate_lawordered its{observable: ∂rate/∂obs}result by iterating a set, so a rate law reading two or more observables ordered everydouble d0/d1temporary — and thev[j]column each is scattered into — byPYTHONHASHSEED.Latent since #76. The saturable twin (
_saturable_jacobian.differentiate_rate_law_native) already sorted for exactly this reason and the sympy branch never did; the single-observable laws that dominate the corpus hid it. Four of the eight models admitted here read ≥ 2 observables, so their emitted C re-hashed on every process and the content-addressed.socould never hit.Fixed, with a
PYTHONHASHSEEDtest intest_codegen_determinism.pythat produces 4 distinct hashes across 4 seeds without it. The fixture needssin()— otherwise the already-sorted native path takes the law and the test passes with the bug in place.Cache key
_switch_sensitivityjoins_CODEGEN_SOURCE_MODULES: the gate lives there now, so that module determines which models get a sensitivity RHS, and #51's digest did not cover it._CODEGEN_VERSION23 → 24.Two deliberate omissions
The test file is not added to mir.yml's pytest list. Its end-to-end cases construct a
Simulatorwithsensitivity_paramson a Functional model — the combination #67 found c2mir rejects (#198,bngsim_codegen_output_sens), and #67 reverted its own mir.yml line over it. That line belongs with the c2mir fix. Only thepaths:filter gains_switch_sensitivity.py, which carries no such risk.event_sensitivity_unsupported_reasonis not folded into the shared predicate, which #68 asks for literally. It answers the same question but on ExprTk addresses, and #52's fix was precisely that observable-vs-species name binding makes text analysis wrong —v > 30binds to the observable total, not the concentration — so re-expressing it against the Python text predicate would reintroduce the bug #52 closed. The two also carry deliberately different policies: the event path refusestime >= p2because event-time sensitivity is unimplemented, while the rate-law path supports it via #48. Cross-references are in both directions in the source.Test plan
python/tests/test_codegen_switch_condition_sens.py— 27 new tests: the rule, derived thresholds, the gate/detector agreement matrix, and end to end. Includes #55d — extend the switch-condition guard so condition-bearing Functional laws can use the analytic sens RHS #68's definition-of-done regression: a fitted state-dependent threshold is refused rather than answered with thePiecewise-derivative zero.mainand on this branch, both emitter paths.PYTHONHASHSEED0/1/2 — 0 differences (4 before the fix).