diff --git a/changelog.d/961.fixed.3.md b/changelog.d/961.fixed.3.md new file mode 100644 index 00000000..25fb18f6 --- /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. A model left with nothing to solve for, every species fixed, left out by a mask, or given by a law of its own (X in `$S + X -> $S + X`), was refused too, every column NaN: each row is a seed or follows from a total, and is returned. 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..0ff7a99e --- /dev/null +++ b/python/tests/test_steady_state_sensitivity_fixed_species.py @@ -0,0 +1,173 @@ +"""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 +""" + + +# Nothing is solved for once the fixed species is out: S alone, and S beside an X +# that S + X -> S + X leaves where it is, a law of its own. +ONLY_FIXED = """begin parameters + 1 S0 3.0 + 2 k 1.0 +end parameters +begin species + 1 $S() S0 +end species +begin reactions + 1 1 1 k +end reactions +""" + +NO_UNKNOWN_LEFT = """begin parameters + 1 S0 3.0 + 2 X0 2.0 + 3 k 1.0 +end parameters +begin species + 1 $S() S0 + 2 X() X0 +end species +begin reactions + 1 1,2 1,2 k +end reactions +""" + + +def _model(tmp_path, text): + path = tmp_path / "m.net" + path.write_text(text) + return bngsim.Model.from_net(str(path)) + + +def _columns(tmp_path, text, params, set_to=None): + model = _model(tmp_path, text) + for name, value in (set_to or {}).items(): + model.set_concentration(name, value) + out = bngsim.Simulator(model, 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) + + +def test_a_fixed_species_set_by_hand_is_set_by_no_parameter(tmp_path): + """Moved off its initial condition, S is where it was put whatever S0 is + (issue #113), with no law and no mask too: at S = 5, dX = [S/k2, -k1·S/k2², 0].""" + k1, k2, s = 2.0, 0.5, 5.0 + want = [[0.0, 0.0, 0.0], [s / k2, -k1 * s / k2**2, 0.0]] + got = _columns(tmp_path, NO_LAW, ["k1", "k2", "S0"], set_to={"S()": s}) + np.testing.assert_allclose(got, want, rtol=1e-6, atol=1e-12) + + +@pytest.mark.parametrize( + ("text", "params", "want"), + [ + (ONLY_FIXED, ["S0", "k"], [[1.0, 0.0]]), + (NO_UNKNOWN_LEFT, ["S0", "X0", "k"], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]), + ], + ids=["only_fixed", "the_rest_in_a_law_of_its_own"], +) +def test_a_model_with_no_unknown_left_is_solved(tmp_path, text, params, want): + """Each row is a seed or a total, with nothing to solve for. Every column + was NaN and refused as singular.""" + np.testing.assert_allclose(_columns(tmp_path, text, params), want, atol=1e-12) + + +@pytest.mark.parametrize("method", ["integration", "newton"]) +@pytest.mark.parametrize( + ("text", "params", "state"), + [ + (NO_LAW, ["k1", "k2", "S0"], [3.0, 12.0]), + (BESIDE_A_LAW, ["k1", "kf", "kr", "S0"], [3.0, 0.25 / 7.25, 7.0 / 7.25]), + ], + ids=["no_law", "beside_a_law"], +) +def test_the_state_is_the_one_found_without_the_columns(tmp_path, method, text, params, state): + """Control. The columns are solved at the state, and leave it where it was + found without them: X* = 12, and B* = a/(a + kr) beside the law.""" + + def solve(**kw): + sim = bngsim.Simulator(_model(tmp_path, text), method="ode") + return np.asarray(sim.steady_state(method=method, **kw).concentrations) + + alone = solve() + np.testing.assert_allclose(alone, state, rtol=1e-7) + np.testing.assert_array_equal(solve(sensitivity_params=params), alone) diff --git a/src/model_builder.cpp b/src/model_builder.cpp index 39903aa1..7c57e81d 100644 --- a/src/model_builder.cpp +++ b/src/model_builder.cpp @@ -1711,7 +1711,9 @@ ConservationLaws detect_conservation_laws(const std::vector &reactions } } - // Build independent species list (all non-dependent, non-fixed) + // Build independent species list: every species no law is solved for, a + // fixed one included (its consumers drop it where they need to, see + // ss_unknown_species) std::vector is_dep(ns, false); for (int d : cl.dependent) is_dep[d] = true; diff --git a/src/steady_state.cpp b/src/steady_state.cpp index 47f1bc10..13d6c6ae 100644 --- a/src/steady_state.cpp +++ b/src/steady_state.cpp @@ -1898,8 +1898,7 @@ static int kinsol_jac(N_Vector y, N_Vector fval, SUNMatrix J, void *ud, N_Vector // Which species a restricted steady-state system treats as UNKNOWNS. // -// Three callers have to agree on this set or they are answering questions about -// different systems: the KINSOL polish (solve_by_newton), the dY_ss/dp linear +// Three callers ask it: the KINSOL polish (solve_by_newton), the dY_ss/dp linear // solve (compute_ss_sensitivity), and the stability certificate at the accepted // root (issue #78). // @@ -1912,8 +1911,11 @@ static int kinsol_jac(N_Vector y, N_Vector fval, SUNMatrix J, void *ud, N_Vector // with the model. // // Unmasked WITH laws returns cl.independent verbatim (fixed species included), -// exactly as before #74. A caller with neither laws nor a mask builds no -// restricted system at all and does not come here. +// exactly as before #74, and the polish and the certificate take it so. The +// dY_ss/dp solve drops the fixed species from it, and comes here with neither +// laws nor a mask where the model has one (issue #961): solved for, its zero row +// made the columns singular. The polish and the certificate build no restricted +// system without laws or a mask, and do not come here then. static std::vector ss_unknown_species(const NetworkModel &model, const ConservationLaws &cl, const ResidualSubspace &sub) { if (!cl.empty() && !sub.masked) { @@ -5530,29 +5532,39 @@ 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 of the columns (issue #961): its row of J + // is 0, and solved for as one it made the solve singular, in a model with no + // law and no mask (the full solve below) as beside a law, among whose + // independents ss_unknown_species keeps it. 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 - // conservation-law independents, narrowed to the species the mask kept - // (issue #74), minus the `fixed` ones. A mask with no conservation laws - // lands here too, with every law loop below a no-op — one code path for - // "some species are not unknowns", mirroring solve_by_newton. - // A fixed species' dY_ss/dp is genuinely 0 — a boundary condition does not - // move with a rate constant — so its row stays 0 rather than NaN: it is - // the mask, not fixedness, that means "no answer". + // `solve_idx` is the unknown set: ss_unknown_species' (the + // conservation-law independents, narrowed to the species the mask kept, + // issue #74) without the `fixed` ones. A model with no law lands here + // too, masked or with a fixed species, with every law loop below a + // no-op — one code path for "some species are not unknowns", mirroring + // solve_by_newton. A fixed species' row is the seed of a parameter that + // sets it (set_by_a_parameter, below) and 0 for any other, not NaN: it + // is the mask, not fixedness, that means "no answer". // // Issue #995: a law that holds a masked-out species is not an equation // of the species the mask kept. Such a species is taken out of the // laws by one of them, which is dropped, and the species that law was // solved for is an unknown, with its own equation. See // ss_kept_species_system. Without a mask, or where no law holds a - // masked species, `cl` and `solve_idx` are the model's laws and - // ss_unknown_species' set, as they were. + // masked species, `cl` is the model's laws and `solve_idx` the unknowns + // as they come here. ConservationLaws cl; std::vector solve_idx; - ss_kept_species_system(model, all_laws, sub, ss_unknown_species(model, all_laws, sub), cl, - solve_idx); + 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 @@ -5585,14 +5597,11 @@ compute_ss_sensitivity(NetworkModel &model, SteadyStateRhs &rhs, SteadyStateResu } } } - if (solve_idx.empty()) { - // Nothing left to differentiate (every included species is pinned by - // a law). Report NaN rather than a zero matrix a fitter would read as - // "this parameter does not matter". - std::fill(result.sensitivity.begin(), result.sensitivity.end(), - std::numeric_limits::quiet_NaN()); - return; - } + // No unknown may be left: every species fixed, left out by the mask, or + // given by a law of its own (X in `$S + X -> $S + X`). Then nothing the + // columns are solved for moves, and each row is a seed or follows from + // a total, with no system to solve (solve_columns). They were NaN, and + // refused as singular where dS/dS0 = 1 (issue #961). const int n_ind = static_cast(solve_idx.size()); // Issue #704: a conserved total is set by the state the solve starts @@ -5653,11 +5662,49 @@ compute_ss_sensitivity(NetworkModel &model, SteadyStateRhs &rhs, SteadyStateResu } } + // The rows of column p that are given once the unknowns' are in `sens`: + // a fixed species a parameter sets is its seed, and each law's dependent + // follows from its law. + const auto fill_given_rows = [&](std::vector &sens, int p) { + for (const int i : set_by_a_parameter) { + sens[static_cast(i) * np + p] = dx0[static_cast(i) * np + p]; + } + + // Reconstruct dependent species sensitivity from conservation: + // Σ L[k,i] * dy_i/dp = dT_k/dp → + // dy_dep/dp = (dT_k/dp - Σ_{i≠dep} L[k,i] * dy_i/dp) / L[k,dep] + // (issue #704: the total's own derivative used to be taken as 0). + for (int k = 0; k < cl.n_laws; ++k) { + int dep = cl.dependent[k]; + double cd = cl.coefficients[k][dep]; + if (std::abs(cd) < 1e-15) + continue; + double s = 0.0; + for (int i = 0; i < ns; ++i) + if (i != dep) + s += cl.coefficients[k][i] * sens[i * np + p]; + // Written as it was where the total does not move, so that such + // a column is the same to the last bit, the sign of a 0 included. + const double moved_total = dT[static_cast(k) * np + p]; + sens[dep * np + p] = moved_total != 0.0 ? (moved_total - s) / cd : -s / cd; + } + }; + // The columns from a Jacobian `Jm`, into `sens` (species-major, ns×np), // and the pivot ratio of its factorization into `rcond` where that is // not null. A function of the matrix so that issue #1002's check can // solve them again with another (ss_confirm_differenced_columns). const auto solve_columns = [&](const double *Jm, std::vector &sens, double *rcond) { + if (n_ind == 0) { + // No system: every row is given (see n_ind above). + if (rcond != nullptr) { + *rcond = 1.0; + } + for (int p = 0; p < np; ++p) { + fill_given_rows(sens, p); + } + return; + } // The reduced Jacobian (n_ind × n_ind), projected through the // conservation-law reconstruction — see ss_reduce_jacobian, which the // stability certificate shares. @@ -5737,28 +5784,7 @@ compute_ss_sensitivity(NetworkModel &model, SteadyStateRhs &rhs, SteadyStateResu // Fill the solved species' sensitivity for (int i = 0; i < n_ind; ++i) sens[solve_idx[i] * np + p] = x_data[i]; - for (const int i : set_by_a_parameter) { - sens[static_cast(i) * np + p] = dx0[static_cast(i) * np + p]; - } - - // Reconstruct dependent species sensitivity from conservation: - // Σ L[k,i] * dy_i/dp = dT_k/dp → - // dy_dep/dp = (dT_k/dp - Σ_{i≠dep} L[k,i] * dy_i/dp) / L[k,dep] - // (issue #704: the total's own derivative used to be taken as 0). - for (int k = 0; k < cl.n_laws; ++k) { - int dep = cl.dependent[k]; - double cd = cl.coefficients[k][dep]; - if (std::abs(cd) < 1e-15) - continue; - double s = 0.0; - for (int i = 0; i < ns; ++i) - if (i != dep) - s += cl.coefficients[k][i] * sens[i * np + p]; - // Written as it was where the total does not move, so that such - // a column is the same to the last bit, the sign of a 0 included. - const double moved_total = dT[static_cast(k) * np + p]; - sens[dep * np + p] = moved_total != 0.0 ? (moved_total - s) / cd : -s / cd; - } + fill_given_rows(sens, p); } }; solve_columns(J.data(), result.sensitivity, &result.sens_jacobian_rcond);