Skip to content

fix: steady-state expression output sensitivities drop the derived-parameter chain rule - #73

Merged
wshlavacek merged 1 commit into
mainfrom
claude/intelligent-cartwright-8fa4a5
Jul 27, 2026
Merged

fix: steady-state expression output sensitivities drop the derived-parameter chain rule#73
wshlavacek merged 1 commit into
mainfrom
claude/intelligent-cartwright-8fa4a5

Conversation

@wshlavacek

Copy link
Copy Markdown
Collaborator

The bug

compute_ss_output_sensitivity's explicit-parameter term ∂func/∂p is a finite difference taken by writing the perturbed value straight into the model's Parameter vector:

const_cast<std::vector<Parameter> &>(params)[pi].value = pval + h;
model.update_observables(y_ss);
model.evaluate_functions(0.0);

Neither update_observables nor evaluate_functions re-derives ConstantExpression parameters — only NetworkModel::set_param does. So when BNG2.pl encodes a rate law or function coefficient as a derived parameter (_rateLaw1 = chi*kon), perturbing the primary kon leaves _rateLaw1 at its nominal value and ∂func/∂_rateLaw1 · ∂_rateLaw1/∂kon silently vanishes from ss.output_sensitivities(["expression:..."]).

Same defect class as #2 / #41, and the same hole compute_ss_sensitivity's ∂f/∂p had until #63 — repaired there by routing parameter writes through SteadyStateRhs::sync_params(held). This site was left behind because it computes a different derivative (∂func/∂p, not ∂f/∂p).

Reproducer

tests/data/ss_expr_sens_derived.netderived_rate_const.net's A ⇌ B with _rateLaw1 = chi*kon, plus flux() = _rateLaw1*A_tot. Reversible, so unlike expr_sens_derived.net it has a non-degenerate steady state. With a = chi·kon and b = koff, conservation gives A+B = 1 and every quantity has a closed form:

A_ss = b/(a+b)          flux = a·A_ss = a·b/(a+b)
d flux/d kon = chi·b²/(a+b)²

which splits into a state-chain term a·dA_ss/dkon = -chi·a·b/(a+b)² and an explicit term A_ss·chi. The explicit term reaches kon only through _rateLaw1, so dropping the chain returns the bare state-chain term — wrong sign, 21× too large:

d flux/d kon d flux/d chi d flux/d koff
before −0.4535 −0.04535 0.9070
after / closed form 0.02268 0.002268 0.9070

Two details worth noting: koff was always right (it reaches the function only through the state), and so was a probe of _rateLaw1 itself (a direct write needs no chain). Only the primaries feeding a derived parameter were wrong — which is why this survived the existing coverage.

The fix

find_steady_state already constructs the SteadyStateRhs and hands it to compute_ss_sensitivity; it now goes to compute_ss_output_sensitivity too. sync_params(pi) runs after the perturb and after the restore, re-deriving every expression parameter except the one being probed — matching set_param's detach-then-refresh rule, so a probe of a derived parameter is not immediately undone.

The restore call matters independently: without it the derived parameters keep the values the last probe gave them and corrupt every later column and the caller's model. test_model_is_left_with_nominal_parameter_values guards that specifically.

The baseline f0 also takes a bare sync_params() first, matching the ordering compute_ss_sensitivity uses for its own f0.

Tests

Four in TestDerivedParameterChain:

  • test_expression_sensitivity_matches_closed_form — the full 4-parameter row including _rateLaw1 itself, plus the state and species sensitivity it rides on, so a failure localizes.
  • test_state_chain_alone_is_not_accepted — pins the failure mode, not just the answer. The two candidates have opposite signs, so an rtol-only assertion on a near-cancelling sum could pass for the wrong reason.
  • test_matches_long_cvode_run — independent oracle: the compiled bngsim_codegen_output_sens chain rule (#198) at the last point of a converged forward-sensitivity run. Agrees with the closed form to 1e-4.
  • test_model_is_left_with_nominal_parameter_values — the restore path.

Three of the four fail without the C++ change (verified by git stash push src/steady_state.cpp + rebuild). Full suite after: 2581 passed, 130 skipped.

Not done here

Replacing the FD block with the compiled bngsim_codegen_output_sens evaluator, which already carries this chain rule analytically. It is emitted only when the model carries _want_output_sens, which Simulator.__init__ sets from its constructor sensitivity_params; steady_state() takes its own as a method argument. Measured consequence:

ctor sensitivity_params=None      → output_sens symbol present=False
ctor sensitivity_params=['kon']   → output_sens symbol present=True

So in the ordinary Simulator(m, method="ode").steady_state(sensitivity_params=[...]) call the symbol is not in the artifact at all. Wiring it up needs the #205 re-prep dance compute_all_sensitivities does (set the flag, drop the plain artifact, regenerate), the dY_ss/dp rows transposed into the per-column pointers the ABI wants, and NaN-sentinel handling for functions the codegen marks unsupported — which this block would still have to answer for. It becomes a preferred path with FD as its fallback, structured exactly as ∂f/∂p is above it, rather than a replacement. Rationale recorded in the comment block so it does not need re-deriving.

…rameter chain rule

`compute_ss_output_sensitivity`'s explicit-parameter term ∂func/∂p is a finite
difference taken by writing the perturbed value straight into the model's
Parameter vector. Neither `update_observables` nor `evaluate_functions`
re-derives ConstantExpression parameters — only `set_param` does. So when
BNG2.pl encodes a rate law or function coefficient as a derived parameter
(`_rateLaw1 = chi*kon`), perturbing the primary `kon` left `_rateLaw1` at its
nominal value and ∂func/∂_rateLaw1 · ∂_rateLaw1/∂kon silently vanished from
`ss.output_sensitivities(["expression:..."])`.

Same defect class as issues #2 / #41, and the same hole `compute_ss_sensitivity`'s
∂f/∂p had until #63 — repaired there by routing parameter writes through
`SteadyStateRhs::sync_params(held)`. This site was left behind because it computes
a different derivative.

The probe now takes the same route: `find_steady_state` passes the
`SteadyStateRhs` it already builds, and `sync_params(pi)` runs after the perturb
AND after the restore, re-deriving every expression parameter except the one being
probed (matching `set_param`'s detach-then-refresh rule, so a probe of a derived
parameter is not immediately undone). The restore call matters independently:
without it the derived parameters keep the values the last probe gave them and
corrupt every later column and the caller's model.

What came back before was not a small error. On `flux() = _rateLaw1*A_tot` over
A <-> B (a = chi*kon, b = koff, flux = a*b/(a+b)), the returned value was the bare
state-chain term a*dA_ss/dp — wrong sign and 21x too large:

| | d flux/d kon | d flux/d chi | d flux/d koff |
|---|---:|---:|---:|
| before | -0.4535 | -0.04535 | 0.9070 |
| after / closed form | 0.02268 | 0.002268 | 0.9070 |

`koff` was always right (it reaches the function only through the state), and so
was a probe of `_rateLaw1` itself (a direct write needs no chain) — only the
primaries feeding a derived parameter were wrong.

Regression coverage: `tests/data/ss_expr_sens_derived.net` is
`derived_rate_const.net`'s A <-> B with `flux() = _rateLaw1*A_tot`. Reversible, so
unlike `expr_sens_derived.net` it has a non-degenerate steady state. Four tests
pin the closed form, the failure mode explicitly (the two candidates have opposite
signs, so an rtol-only assertion could pass for the wrong reason), the CVODES
`bngsim_codegen_output_sens` chain rule as an independent oracle, and parameter
restoration. Three of the four fail without the fix.

Not done here: replacing the FD block with the compiled `bngsim_codegen_output_sens`
evaluator. It is emitted only when the model carries `_want_output_sens`, which
`Simulator.__init__` sets from its CONSTRUCTOR `sensitivity_params`;
`steady_state()` takes its own as a METHOD argument, so the artifact reaching the
solver in the ordinary call has no such symbol. Rationale recorded in the comment
block.
@wshlavacek
wshlavacek force-pushed the claude/intelligent-cartwright-8fa4a5 branch from e083da5 to 47e5111 Compare July 27, 2026 01:47
@wshlavacek
wshlavacek merged commit abe8cc4 into main Jul 27, 2026
2 checks passed
@wshlavacek
wshlavacek deleted the claude/intelligent-cartwright-8fa4a5 branch July 27, 2026 01:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant