From 72d4eee4aa369e31611a751358805e9adced6f33 Mon Sep 17 00:00:00 2001 From: Montana Date: Fri, 30 Oct 2026 18:11:48 -0700 Subject: [PATCH] steady_state: solve the sensitivity of a model with a fixed species, not refuse it as singular (#961) A $-fixed species has no equation, and its row of J is 0. The columns were solved with it among the unknowns where a model had no law and no mask, and beside a law, whose independents it was among: J was singular and dY_ss/dp was refused as not existing. It is taken out of the columns' unknowns; the solve for the state is unchanged. 44 corpus models solve where they were refused, each agreeing with central differences, and the 566 that solved before are bit-identical. --- changelog.d/961.fixed.3.md | 1 + ..._steady_state_sensitivity_fixed_species.py | 105 ++++++++++++++++++ src/steady_state.cpp | 17 ++- 3 files changed, 120 insertions(+), 3 deletions(-) create mode 100644 changelog.d/961.fixed.3.md create mode 100644 python/tests/test_steady_state_sensitivity_fixed_species.py diff --git a/changelog.d/961.fixed.3.md b/changelog.d/961.fixed.3.md new file mode 100644 index 00000000..b2afeede --- /dev/null +++ b/changelog.d/961.fixed.3.md @@ -0,0 +1 @@ +- **`steady_state(sensitivity_params=...)` refused a model with a `$`-fixed species as singular** ([#961](https://github.com/lanl/bngsim/issues/961)). A fixed species has no equation, and its row of J is 0. With no conservation law and no mask, the columns were solved on every species; beside a law, the fixed species was among the law's independents. Either way the solve was singular, and the columns were refused as "dY_ss/dp does not exist". The fixed species is taken out of the columns' unknowns now, and its own row is what a parameter sets it to. On the BioModels corpus and the repository's `.net` files, 44 models solve where they were refused, each agreeing with central differences of its steady state, and the 566 that solved before are bit-identical. The solve for the state itself is unchanged. diff --git a/python/tests/test_steady_state_sensitivity_fixed_species.py b/python/tests/test_steady_state_sensitivity_fixed_species.py new file mode 100644 index 00000000..d6bbaa44 --- /dev/null +++ b/python/tests/test_steady_state_sensitivity_fixed_species.py @@ -0,0 +1,105 @@ +"""The steady-state sensitivity of a model with a ``$``-fixed species is +solved, not refused as singular (issue #961). + +A fixed species has no equation: its row of J is 0. The columns are solved on +the unknowns, and a fixed species was kept among them in two cases: a model with +no conservation law and no mask, which took the full solve, and one with a law, +whose independent species it was among. Either way J was singular there, and +``steady_state(sensitivity_params=...)`` raised "dY_ss/dp does not exist" for +columns that exist. The fixed species is taken out of the unknowns now; its own +row is what a parameter sets it to. +""" + +from __future__ import annotations + +import bngsim +import numpy as np +import pytest + +# X is made at k1 per S and lost at k2: X* = k1·S/k2. +NO_LAW = """begin parameters + 1 k1 2.0 + 2 k2 0.5 + 3 S0 3.0 +end parameters +begin species + 1 $S() S0 + 2 X() 0 +end species +begin reactions + 1 1 1,2 k1 + 2 2 0 k2 +end reactions +""" + +# A <-> B at kf, kr, and S turns A into B at k1: with a = kf + k1·S and +# A + B = 1, B* = a/(a + kr). +BESIDE_A_LAW = """begin parameters + 1 k1 2.0 + 2 kf 1.0 + 3 kr 0.25 + 4 S0 3.0 +end parameters +begin species + 1 $S() S0 + 2 A() 1 + 3 B() 0 +end species +begin reactions + 1 2 3 kf + 2 3 2 kr + 3 1,2 1,3 k1 +end reactions +""" + + +def _columns(tmp_path, text, params): + path = tmp_path / "m.net" + path.write_text(text) + out = bngsim.Simulator(bngsim.Model.from_net(str(path)), method="ode").steady_state( + sensitivity_params=params + ) + assert out.converged + return np.asarray(out.sensitivity) + + +def test_a_fixed_species_with_no_law_is_solved(tmp_path): + """Refused as singular. dS/d[k1, k2, S0] = [0, 0, 1]; dX = [S/k2, -k1·S/k2², k1/k2].""" + k1, k2, s = 2.0, 0.5, 3.0 + want = [[0.0, 0.0, 1.0], [s / k2, -k1 * s / k2**2, k1 / k2]] + np.testing.assert_allclose( + _columns(tmp_path, NO_LAW, ["k1", "k2", "S0"]), want, rtol=1e-6, atol=1e-12 + ) + + +def test_a_fixed_species_beside_a_conservation_law_is_solved(tmp_path): + """Refused as singular, the fixed species among the law's independents.""" + k1, kf, kr, s = 2.0, 1.0, 0.25, 3.0 + a = kf + k1 * s + d = (a + kr) ** 2 + db = [s * kr / d, kr / d, -a / d, k1 * kr / d] # d/d[k1, kf, kr, S0] + want = [[0.0, 0.0, 0.0, 1.0], [-x for x in db], db] + np.testing.assert_allclose( + _columns(tmp_path, BESIDE_A_LAW, ["k1", "kf", "kr", "S0"]), want, rtol=1e-6, atol=1e-12 + ) + + +def test_a_model_with_no_fixed_species_is_solved_as_before(tmp_path): + """Control. 0 -> X at k1, X -> 0 at k2: dX*/d[k1, k2] = [1/k2, -k1/k2²].""" + text = ( + "begin parameters\n 1 k1 2.0\n 2 k2 0.5\nend parameters\n" + "begin species\n 1 X() 0\nend species\n" + "begin reactions\n 1 0 1 k1\n 2 1 0 k2\nend reactions\n" + ) + np.testing.assert_allclose(_columns(tmp_path, text, ["k1", "k2"]), [[2.0, -8.0]], rtol=1e-6) + + +@pytest.mark.parametrize("method", ["integration", "newton"]) +def test_the_state_itself_is_what_it_was(tmp_path, method): + """The unknowns of the solve for the state are not changed: X* = 12.""" + path = tmp_path / "m.net" + path.write_text(NO_LAW) + out = bngsim.Simulator(bngsim.Model.from_net(str(path)), method="ode").steady_state( + method=method + ) + np.testing.assert_allclose(out.concentrations, [3.0, 12.0], rtol=1e-7) diff --git a/src/steady_state.cpp b/src/steady_state.cpp index e12f959f..7469be00 100644 --- a/src/steady_state.cpp +++ b/src/steady_state.cpp @@ -5157,7 +5157,12 @@ compute_ss_sensitivity(NetworkModel &model, SteadyStateRhs &rhs, SteadyStateResu column_system.param_values.push_back(params[static_cast(pidx[p])].value); } - if (!all_laws.empty() || sub.masked) { + // A `$`-fixed species is no unknown either (issue #961): solved for as one, + // its zero row made J singular in a model with no law and no mask, and the + // columns were refused as "dY_ss/dp does not exist". + const bool any_fixed = std::any_of(model.species().begin(), model.species().end(), + [](const Species &s) { return s.fixed; }); + if (!all_laws.empty() || sub.masked || any_fixed) { // Reduced-space sensitivity solve. // // `solve_idx` is the unknown set (ss_unknown_species): the @@ -5178,8 +5183,14 @@ compute_ss_sensitivity(NetworkModel &model, SteadyStateRhs &rhs, SteadyStateResu // ss_unknown_species' set, as they were. ConservationLaws cl; std::vector solve_idx; - ss_kept_species_system(model, all_laws, sub, ss_unknown_species(model, all_laws, sub), cl, - solve_idx); + // A fixed species is no unknown of the columns (issue #961): its row of + // J is 0, and kept, it made the solve singular, beside a law as well as + // with none. Its own row of dY_ss/dp is what a parameter sets it to. + std::vector unknowns = ss_unknown_species(model, all_laws, sub); + unknowns.erase(std::remove_if(unknowns.begin(), unknowns.end(), + [&](int i) { return model.species()[i].fixed; }), + unknowns.end()); + ss_kept_species_system(model, all_laws, sub, unknowns, cl, solve_idx); // Issue #995: a species the mask left out, that no law gives, is held // where the solve left it. That is its part in the derivative only // where no equation the columns are solved on reads it: `A <-> M` with