From f97ac7bf3ea839c9c733279ee463a14acf293b2f Mon Sep 17 00:00:00 2001 From: Montana Date: Fri, 30 Oct 2026 11:18:11 -0700 Subject: [PATCH 1/3] steady_state: solve a model whose only read of the time is a function that reports (#961) #710 refused steady_state() wherever any function read the time, which refused minutes := time/60 in BIOMD0000000739 and 951 though nothing in the dynamics reads it. The refusal now asks what the time reaches: a rate law (rate_laws_read_clock, which does not count a rateOf), or a rule a species or compartment size is reported through. Where it reaches neither, the state is solved for and each function that reads the clock (functions_reading_clock) is reported as NaN, with a NaN sensitivity row output_sensitivities refuses by that reason. --- changelog.d/961.changed.md | 1 + docs/user-guide/steady-state.md | 14 ++- include/bngsim/model.hpp | 13 +++ python/bngsim/_bngsim_core.pyi | 10 ++ python/bngsim/_simulator.py | 54 ++++++++-- ...steady_state_of_a_model_that_reads_time.py | 98 +++++++++++++++++-- src/_bngsim_core.cpp | 17 ++++ src/model.cpp | 26 ++++- src/steady_state.cpp | 28 +++++- 9 files changed, 239 insertions(+), 22 deletions(-) create mode 100644 changelog.d/961.changed.md diff --git a/changelog.d/961.changed.md b/changelog.d/961.changed.md new file mode 100644 index 000000000..2ff432927 --- /dev/null +++ b/changelog.d/961.changed.md @@ -0,0 +1 @@ +- **`steady_state()` solves a model whose only read of the time is a function nothing in the dynamics reads** ([#961](https://github.com/lanl/bngsim/issues/961)). `minutes := time/60` in BIOMD0000000739 and 951, and `parameter_1`'s pulse in BIOMD0000000498, made #710 refuse the whole model, though the state does not depend on them. The refusal now asks what the time reaches: a rate law (`NetworkModel.rate_laws_read_clock`, which does not count a `rateOf`), or a rule a species or compartment size is reported through. Where it reaches neither, the state is solved for, and each function that reads the clock (`functions_reading_clock`) is reported as NaN with a NaN sensitivity row, which `output_sensitivities` refuses by name. BIOMD0000000739's steady state matches a run to t = 1e9 at rtol 1e-12 to 3e-7; 498 and 951 report `converged=False`, as their runs still drift. diff --git a/docs/user-guide/steady-state.md b/docs/user-guide/steady-state.md index 1eb8482ac..6c2214df7 100644 --- a/docs/user-guide/steady-state.md +++ b/docs/user-guide/steady-state.md @@ -26,9 +26,17 @@ a root of `f(y)` with everything read at `t = 0`, and fire no event: `A -> B` at state as converged, and a model whose event switches a production rate on came back at the state from before it. A reported quantity that reads the time has no steady value either: an assignment-rule species `S := B*(1 - exp(-time))` -came back 0 for 2/3. A `time()` call or a table function indexed by time, in a -rate law, a rule or a function that is only reported, is a read of the time; -`rateOf` is not. +came back 0 for 2/3, and so has a compartment sized by a rule on the time. A +`time()` call or a table function indexed by time, in a rate law or in a rule a +species or compartment size is reported through, directly or through what it +reads, is a read of the time; `rateOf` is not. + +A function that reads the time and that only reports, `minutes := time/60` as +BIOMD0000000739 and 951 have it, is not refused (issue #961): the state does +not depend on it and is solved for as in a model without it. The function has +no steady value, so its entry in the result's expression values, and its row of +the expression sensitivities, is NaN, and `output_sensitivities` refuses it by +that reason. Integrate such a model with `run()` over a span long enough for the trajectory to settle. `run(..., steady_state=True)` does not stop early for it: the diff --git a/include/bngsim/model.hpp b/include/bngsim/model.hpp index a37c2eca5..3f684a773 100644 --- a/include/bngsim/model.hpp +++ b/include/bngsim/model.hpp @@ -337,6 +337,18 @@ class NetworkModel { // A trigger that does not can change only when the state does. bool event_trigger_reads_time(int event_idx0) const; + // Does any reaction's rate read the clock (issue #961)? The same walk as + // reaction_rate_reads_time, without a rate accessor, which reads the running + // derivatives and not the time: at a steady state they are 0 whatever the + // time. A model whose rates read no clock has an f(y) = 0 to solve even + // where a function that only reports reads it. + bool rate_laws_read_clock() const; + + // The 0-based indices of the functions whose value moves with the clock, + // through `time`, a time-indexed table, or a function or derived parameter + // that does (issue #961). A steady state has no value for one. + std::vector functions_reading_clock() const; + // ─── Rate-law switch conditions that read model state (issue #150) ─────── // // The rate-law twin of the state-dependent event trigger issue #144 covers. @@ -941,6 +953,7 @@ class NetworkModel { // What a rate or a trigger reads through the model's definitions (#719). struct RateDeps { bool time = false; // the clock, directly or through a definition + bool clock = false; // the same, without a rate accessor (issue #961) bool unknown = false; // something whose reads cannot be named std::set species; }; diff --git a/python/bngsim/_bngsim_core.pyi b/python/bngsim/_bngsim_core.pyi index abc6034ea..5dd00bbd4 100644 --- a/python/bngsim/_bngsim_core.pyi +++ b/python/bngsim/_bngsim_core.pyi @@ -384,6 +384,11 @@ class NetworkModel: Names of the model's functions, in declaration order. Each also names a parameter slot holding its evaluated value, which ``param_is_internal`` flags and ``primary_param_names`` omits (issue #227). """ @property + def functions_reading_clock(self) -> list[str]: + """ + The functions whose value moves with the clock (issue #961), by name: a steady state reports each as NaN. + """ + @property def functions_use_time(self) -> bool: """ Whether any function's value can move with simulation time alone (issue #654). Gates the SSA's time-dependent sub-stepping. @@ -482,6 +487,11 @@ class NetworkModel: Parameter names labeling the columns of pending_sensitivity_seed() (GH #210). """ @property + def rate_laws_read_clock(self) -> bool: + """ + Whether a reaction's rate reads the clock, through time, a time-indexed table, or a function or derived parameter that does; a rate accessor is not the clock (issue #961). A model whose rates read none has a steady state to solve for. + """ + @property def rhs_eval_count(self) -> int: """ Number of RHS (compute_derivs_core) evaluations so far. diff --git a/python/bngsim/_simulator.py b/python/bngsim/_simulator.py index 2e0cc2ec9..c81bee403 100644 --- a/python/bngsim/_simulator.py +++ b/python/bngsim/_simulator.py @@ -7893,6 +7893,7 @@ def steady_state( raise SimulationError(f"Steady-state computation failed: {e}") from e result = SteadyStateResult(core_result) + result._expressions_reading_clock = frozenset(self._model._core.functions_reading_clock) self._note_ss_unsupported(result) self._apply_ss_ar_report_map(result) self._apply_ss_ar_sensitivity_map(result) @@ -8589,7 +8590,8 @@ def _raise_if_no_steady_state_to_solve_for(self, where: str) -> None: ``t = 0`` beside a state marked converged: an assignment-rule species ``S := B*(1 - exp(-time))`` as 0 for 2/3, a species in a compartment sized by a rule on the time at 0.5 for 0.25. A function that reads the - time and that nothing reads is refused with the rest. + time and that only reports is not refused: the state does not depend on + it, and its own value is NaN (issue #961). A ``rateOf`` accessor is not a read of the time: it reads the state's own derivatives, which are 0 at a steady state. @@ -8625,18 +8627,41 @@ def _raise_if_no_steady_state_to_solve_for(self, where: str) -> None: "jump (issue #710). Integrate the model with run() over a span long enough for " "the trajectory to settle, and read the state where it ends." ) - if core.functions_use_time: + if core.functions_use_time and self._ss_reads_the_time(core): raise SimulationError( f"{where} is not supported for this model: it reads the time, through time() " - "or a table function indexed by time, in a rate law, a rule or a reported " - "function. The steady-state solvers look for a root of f(y) with everything " - "read at t = 0, which is the model's steady state only where nothing changes " - "with the time after that: a rate that is 0 at t = 0 returned the initial " - "state as converged, and a reported function its value at t = 0 (issue #710). " - "Integrate the model with run() over a span long enough for the trajectory to " - "settle, and read the state where it ends." + "or a table function indexed by time, in a rate law, a rule a species or a " + "compartment size is reported through, or what one of those reads. The " + "steady-state solvers look for a root of f(y) with everything read at t = 0, " + "which is the model's steady state only where nothing changes with the time " + "after that: a rate that is 0 at t = 0 returned the initial state as " + "converged, and a species set by a rule on the time its value at t = 0 " + "(issue #710). Integrate the model with run() over a span long enough for the " + "trajectory to settle, and read the state where it ends." ) + def _ss_reads_the_time(self, core) -> bool: + """Whether the time reaches what a steady state solves for or reports as a + state (issues #710, #961). + + A rate that reads the clock has no ``f(y) = 0`` without a time being chosen. + A rule on the clock that a species is reported through (``S := B*(1 - + exp(-time))``) or a compartment is sized by has no steady value, and the + state reported through it is not one. A function that reads the time and + that only reports, ``minutes := time/60`` in BIOMD0000000739 and 951, is + neither: the state is solved for as in a model without it, and the + function's value is NaN (``find_steady_state``). + """ + if core.rate_laws_read_clock: + return True + clocked = set(core.functions_reading_clock) + if not clocked: + return False + amap = getattr(self._model, "_ar_report_map", None) or {} + if any(entry[1] in clocked for entry in amap.values()): + return True + return any(name in clocked for name in self._model.compartment_size_params) + def _ss_mask_excluding_ar_species(self, mask_selector: list[int] | None) -> list[int] | None: """Drop AssignmentRule-target species from the steady-state subspace (#247). @@ -9284,6 +9309,7 @@ def _run_one(i): except RuntimeError as e: raise SimulationError(f"Batch {i} failed: {e}") from e result = SteadyStateResult(core_result) + result._expressions_reading_clock = frozenset(clone._core.functions_reading_clock) # GH #247 — the same two passes steady_state() runs, against the # CLONE: a scan that moves a compartment volume moves the divisor with # it, so resolving vdiv off self._model would report every entry @@ -11349,6 +11375,7 @@ class SteadyStateResult: "_expression_values", "_rule_target_species", "_rule_target_expressions", + "_expressions_reading_clock", ) def __init__(self, core) -> None: @@ -11501,6 +11528,9 @@ def _named(names: Sequence[str], index: int) -> str | None: self._expression_values = np.asarray( getattr(core, "expression_values", []), dtype=np.float64 ) + # Issue #961: the functions that read the time, reported as NaN, which + # the solve fills in from the model. + self._expressions_reading_clock: frozenset[str] = frozenset() @property def concentrations(self): @@ -11694,6 +11724,12 @@ def _no_derivative_reason(self, meta: dict[str, Any], row: np.ndarray) -> str | return None kind = meta["kind"] values = self._expression_values if kind == "expression" else self._observable_values + if kind == "expression" and meta["name"] in self._expressions_reading_clock: + return ( + f"output_sensitivities: {meta['selector']!r} has no steady-state " + "derivative: it reads the time, so it has no steady value, and it is " + "reported as NaN (issue #961)." + ) if kind != "species" and meta["index"] < values.size: value = float(values[meta["index"]]) if not np.isfinite(value): diff --git a/python/tests/test_steady_state_of_a_model_that_reads_time.py b/python/tests/test_steady_state_of_a_model_that_reads_time.py index 34ea435c6..60e7b9867 100644 --- a/python/tests/test_steady_state_of_a_model_that_reads_time.py +++ b/python/tests/test_steady_state_of_a_model_that_reads_time.py @@ -12,8 +12,13 @@ A rate that reads the time has no ``f(y) = 0`` without a time being chosen, and an event needs a trajectory. Both are refused now, by name, as is a model -where a reported quantity reads the time and no rate does: that came back at -its value at ``t = 0`` beside a state marked converged. +whose species or compartment size is reported through a rule on the time: that +came back at its value at ``t = 0`` beside a state marked converged. + +A function that reads the time and that only reports, ``minutes := time/60``, +was refused with them, though the state does not depend on it (issue #961). +The state is solved for now, and the function reported as NaN, having no +steady value. The early stop of ``run(steady_state=True)`` tests ``‖f(t, y)‖`` at an output point, which for such a model says nothing of the next: a rate switched on at @@ -161,11 +166,92 @@ def test_the_run_the_message_points_at_reaches_the_steady_state(tmp_path): assert np.asarray(out.species)[-1, names.index("A")] == pytest.approx(1.0, abs=1e-6) -def test_a_reported_function_that_reads_the_time_is_refused(tmp_path): - """``seen() = 2*time()`` is read by no reaction, and is reported: it came - back 0, its value at t = 0, beside a state marked converged.""" +def _expression(result, name): + return float(np.asarray(result._expression_values)[list(result.expression_names).index(name)]) + + +@pytest.mark.parametrize("method", ["integration", "newton"]) +def test_a_function_that_reads_the_time_and_only_reports_is_nan(tmp_path, method): + """``seen() = 2*time()`` is read by no reaction: A -> B at k settles at + [0, 1] whatever the time. #710 refused it, and before that it came back 0, + its value at t = 0 (issue #961).""" sim = bngsim.Simulator(_ramped(tmp_path, law="k"), method="ode") - _refused(sim.steady_state, "it reads the time, through") + out = sim.steady_state(method=method) + assert out.converged + np.testing.assert_allclose(out.concentrations, [0.0, 1.0], atol=1e-8) + assert np.isnan(_expression(out, "seen")) + assert _expression(out, "kt") == pytest.approx(1.0) + + +def test_its_sensitivity_row_is_nan_and_the_state_s_is_not(tmp_path): + """dB/dk = 0 at a steady state of A -> B; seen's row has no value.""" + sim = bngsim.Simulator(_ramped(tmp_path, law="k"), method="ode") + out = sim.steady_state(sensitivity_params=["k"]) + np.testing.assert_allclose(np.asarray(out.sensitivity).ravel(), [0.0, 0.0], atol=1e-8) + names = list(out.expression_names) + rows = np.asarray(out.sensitivities_expressions) + assert np.isnan(rows[names.index("seen")]).all() + assert np.isfinite(rows[names.index("kt")]).all() + with pytest.raises(ValueError, match=r"'expression:seen'.*reads the time.*issue #961"): + out.output_sensitivities("expression:seen") + + +def test_a_batch_of_such_a_model_is_solved(tmp_path): + rows = bngsim.Simulator(_ramped(tmp_path, law="k"), method="ode").steady_state_batch( + [{"k": 1.0}, {"k": 2.0}] + ) + for row in rows: + np.testing.assert_allclose(row.concentrations, [0.0, 1.0], atol=1e-8) + assert np.isnan(_expression(row, "seen")) + + +def test_a_function_that_reads_the_time_through_another_is_nan(tmp_path): + path = tmp_path / "nested_report.net" + path.write_text( + RAMPED.format(law="k", rate="kt").replace( + " 2 seen() 2*time()\n", " 2 seen() 2*time()\n 3 later() 1 + seen()\n" + ) + ) + out = bngsim.Simulator(bngsim.Model.from_net(str(path)), method="ode").steady_state() + assert out.converged + assert np.isnan(_expression(out, "later")) + + +def test_a_table_on_the_time_that_only_reports_is_nan(tmp_path): + path = tmp_path / "table_report.net" + path.write_text( + RAMPED.format(law="k", rate="kt").replace( + " 2 seen() 2*time()\n", " 2 seen() tfun([0,1,2,50],[0,0.5,1,1],time)\n" + ) + ) + out = bngsim.Simulator(bngsim.Model.from_net(str(path)), method="ode").steady_state() + np.testing.assert_allclose(out.concentrations, [0.0, 1.0], atol=1e-8) + assert np.isnan(_expression(out, "seen")) + + +def test_an_sbml_rule_on_the_time_that_only_reports_is_nan(): + """BIOMD0000000739 and 951 report the time in minutes, ``minutes := + time/60``, and nothing reads it.""" + text = "species A, B; A = 1; B = 0; kf = 1; kr = 0.5; minutes := time/60\n" + text += "J1: A -> B; kf*A; J2: B -> A; kr*B\n" + sim = bngsim.Simulator(bngsim.Model.from_antimony_string(text), method="ode") + out = sim.steady_state() + names = list(out.species_names) + got = np.asarray(out.concentrations) + assert got[names.index("A")] == pytest.approx(1 / 3, rel=1e-6) + assert got[names.index("B")] == pytest.approx(2 / 3, rel=1e-6) + assert np.isnan(_expression(out, "minutes")) + + +def test_the_core_solver_reports_such_a_function_as_nan(tmp_path): + from bngsim._bngsim_core import SteadyStateOptions, find_steady_state + + core = _ramped(tmp_path, law="k")._core + assert core.functions_use_time and not core.rate_laws_read_clock + assert core.functions_reading_clock == ["seen"] + out = find_steady_state(core, SteadyStateOptions()) + values = dict(zip(out.expression_names, out.expression_values, strict=False)) + assert np.isnan(values["seen"]) and values["kt"] == pytest.approx(1.0) def test_the_time_written_bare_is_refused(tmp_path): diff --git a/src/_bngsim_core.cpp b/src/_bngsim_core.cpp index 78dc319cd..c1c027b4f 100644 --- a/src/_bngsim_core.cpp +++ b/src/_bngsim_core.cpp @@ -1795,6 +1795,23 @@ PYBIND11_MODULE(_bngsim_core, m) { .def_property_readonly("functions_use_time", &bngsim::NetworkModel::functions_use_time, "Whether any function's value can move with simulation time alone " "(issue #654). Gates the SSA's time-dependent sub-stepping.") + .def_property_readonly( + "rate_laws_read_clock", &bngsim::NetworkModel::rate_laws_read_clock, + "Whether a reaction's rate reads the clock, through time, a time-indexed table, or a " + "function or derived parameter that does; a rate accessor is not the clock (issue " + "#961). A model whose rates read none has a steady state to solve for.") + .def_property_readonly( + "functions_reading_clock", + [](const bngsim::NetworkModel &self) { + std::vector names; + const auto &funcs = self.functions(); + for (const int f : self.functions_reading_clock()) { + names.push_back(funcs[static_cast(f)].name); + } + return names; + }, + "The functions whose value moves with the clock (issue #961), by name: a steady " + "state reports each as NaN.") .def_property_readonly( "ssa_reads_clock", [](const bngsim::NetworkModel &self) { diff --git a/src/model.cpp b/src/model.cpp index 727864c65..1430a0d38 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -2187,14 +2187,17 @@ NetworkModel::rate_dependencies_(std::vector params, const std::vector(e[k]))) ++k; const bool call = k < e.size() && e[k] == '('; - if (clock && id == "time" && (call || !time_declared)) + if (clock && id == "time" && (call || !time_declared)) { deps.time = true; + deps.clock = true; + } auto tit = tables.find(id); if (tit != tables.end()) { const std::string &idx = tit->second->index_name(); const std::string key = strip_paren_suffix(idx); if (is_time_index(idx)) { deps.time = true; + deps.clock = true; } else if (auto pit = sd.param_name_to_idx.find(key); pit != sd.param_name_to_idx.end()) { params.push_back(pit->second); @@ -2311,6 +2314,27 @@ bool NetworkModel::reaction_rate_reads_time(int rxn_idx0, return reaction_rate_dependencies_(rxn_idx0, pc_functions).time; } +bool NetworkModel::rate_laws_read_clock() const { + for (int r = 0; r < static_cast(impl_->shared->reactions.size()); ++r) { + if (reaction_rate_dependencies_(r).clock) + return true; + } + return false; +} + +std::vector NetworkModel::functions_reading_clock() const { + std::vector out; + for (int f = 0; f < static_cast(impl_->functions.size()); ++f) { + const Function &fn = impl_->functions[f]; + std::vector params; + if (fn.evaluator_id >= 0) + expression_support(fn.evaluator_id, nullptr, ¶ms); + if (rate_dependencies_(std::move(params), {fn.expression, fn.eval_expression}).clock) + out.push_back(f); + } + return out; +} + bool NetworkModel::event_trigger_reads_time(int event_idx0) const { const auto &events = impl_->events; if (event_idx0 < 0 || event_idx0 >= static_cast(events.size())) diff --git a/src/steady_state.cpp b/src/steady_state.cpp index e12f959ff..57fb36091 100644 --- a/src/steady_state.cpp +++ b/src/steady_state.cpp @@ -5933,10 +5933,13 @@ SteadyStateResult find_steady_state(NetworkModel &model, const SteadyStateOption "steady-state solvers never make (issue #1054, issue #710). Integrate it with " "run()."); } - if (model.n_events() > 0 || model.functions_use_time()) { + // A function that reads the time and that no rate reads, `minutes := + // time/60`, leaves f(y) = 0 one to solve: it is reported as NaN below, + // having no steady value (issue #961). + if (model.n_events() > 0 || model.rate_laws_read_clock()) { throw std::runtime_error( - "Cannot find steady state: the model has an event or reads the time, and the " - "steady-state solvers read the right-hand side at t = 0 and fire no event " + "Cannot find steady state: the model has an event or a rate that reads the time, " + "and the steady-state solvers read the right-hand side at t = 0 and fire no event " "(issue #710). Integrate it with run()."); } @@ -6455,6 +6458,25 @@ SteadyStateResult find_steady_state(NetworkModel &model, const SteadyStateOption ss_newton_distance(model, rhs, *effective, sub, result.concentrations.data()); } + // A function that moves with the clock has no steady value: what was read + // is its value at t = 0 (issue #961). No rate reads one, or the model was + // refused above. + if (model.functions_use_time()) { + const double nan = std::numeric_limits::quiet_NaN(); + const std::size_t n_func = result.function_values.size(); + const std::size_t np = n_func == 0 ? 0 : result.function_sensitivity.size() / n_func; + for (const int f : model.functions_reading_clock()) { + const auto fu = static_cast(f); + if (fu < n_func) { + result.function_values[fu] = nan; + } + for (std::size_t p = 0; p < np && (fu + 1) * np <= result.function_sensitivity.size(); + ++p) { + result.function_sensitivity[fu * np + p] = nan; + } + } + } + return result; } From a0bbb10c4b721631ad036e0fd7f5636730d3ff21 Mon Sep 17 00:00:00 2001 From: Montana Date: Fri, 30 Oct 2026 11:37:15 -0700 Subject: [PATCH 2/3] ci: rerun after the macOS 14 runner brownout From d4d08c4aa9d695ab45fa016d822ccde55de9484f Mon Sep 17 00:00:00 2001 From: Montana Date: Fri, 30 Oct 2026 14:30:27 -0700 Subject: [PATCH 3/3] run: keep the early stop beside a function that reads the time and only reports (#961) #710 held off run(steady_state=True)'s early stop wherever any function read the time. Its test, ||f(y)|| at an output point, says as much in a model whose rates read no clock as in one with no such function, so the solver's two gates now test rate_laws_read_clock. The Python layer still holds the stop off, in run and in each run_batch row, where the time reaches a species an assignment rule reports or a compartment size. --- changelog.d/961.changed.md | 2 +- docs/user-guide/steady-state.md | 3 ++ python/bngsim/_simulator.py | 15 ++++++- ...steady_state_of_a_model_that_reads_time.py | 45 +++++++++++++++++++ src/cvode_simulator.cpp | 9 +++- 5 files changed, 69 insertions(+), 5 deletions(-) diff --git a/changelog.d/961.changed.md b/changelog.d/961.changed.md index 2ff432927..8404f86b7 100644 --- a/changelog.d/961.changed.md +++ b/changelog.d/961.changed.md @@ -1 +1 @@ -- **`steady_state()` solves a model whose only read of the time is a function nothing in the dynamics reads** ([#961](https://github.com/lanl/bngsim/issues/961)). `minutes := time/60` in BIOMD0000000739 and 951, and `parameter_1`'s pulse in BIOMD0000000498, made #710 refuse the whole model, though the state does not depend on them. The refusal now asks what the time reaches: a rate law (`NetworkModel.rate_laws_read_clock`, which does not count a `rateOf`), or a rule a species or compartment size is reported through. Where it reaches neither, the state is solved for, and each function that reads the clock (`functions_reading_clock`) is reported as NaN with a NaN sensitivity row, which `output_sensitivities` refuses by name. BIOMD0000000739's steady state matches a run to t = 1e9 at rtol 1e-12 to 3e-7; 498 and 951 report `converged=False`, as their runs still drift. +- **`steady_state()` solves a model whose only read of the time is a function nothing in the dynamics reads** ([#961](https://github.com/lanl/bngsim/issues/961)). `minutes := time/60` in BIOMD0000000739 and 951, and `parameter_1`'s pulse in BIOMD0000000498, made #710 refuse the whole model, though the state does not depend on them. The refusal now asks what the time reaches: a rate law (`NetworkModel.rate_laws_read_clock`, which does not count a `rateOf`), or a rule a species or compartment size is reported through. Where it reaches neither, the state is solved for, and each function that reads the clock (`functions_reading_clock`) is reported as NaN with a NaN sensitivity row, which `output_sensitivities` refuses by name. The early stop of `run(steady_state=True)` and `run_batch(steady_state=True)` is lifted for the same models, on the same test, and stops where the model without that function would. BIOMD0000000739's steady state matches a run to t = 1e9 at rtol 1e-12 to 3e-7; 498 and 951 report `converged=False`, as their runs still drift. diff --git a/docs/user-guide/steady-state.md b/docs/user-guide/steady-state.md index 6c2214df7..4006ad9fd 100644 --- a/docs/user-guide/steady-state.md +++ b/docs/user-guide/steady-state.md @@ -43,6 +43,9 @@ to settle. `run(..., steady_state=True)` does not stop early for it: the criterion is `‖f(t, y)‖` at an output point, and a rate that is switched on at `t = 5` is 0 at `t = 1`. The run goes to the end of its span, and `steady_state_reached` stays 0. +A function that reads the time and that only reports leaves the criterion as +good as in a model without it, and the run stops early as that model would +(issue #961); the rows it records read the time where they were taken. **`method="newton"`**: the two-tier integrate-first solver. Tier 1 is the *same* CVODE burst as `"integration"`, carrying the state into the physical diff --git a/python/bngsim/_simulator.py b/python/bngsim/_simulator.py index c81bee403..6094216a7 100644 --- a/python/bngsim/_simulator.py +++ b/python/bngsim/_simulator.py @@ -5150,7 +5150,7 @@ def run( opts.force_sparse_linear_solver = self._force_sparse_linear_solver opts.continues_trajectory = self._continues_trajectory opts.timeout_seconds = timeout_seconds - opts.steady_state = bool(steady_state) + opts.steady_state = bool(steady_state) and not self._early_stop_reads_time() opts.steady_state_tol = ss_tol_value opts.carry_sensitivities = bool(carry_sensitivities) # Seed for random equal-priority event tie-breaking (GH #242). @@ -6615,7 +6615,7 @@ def _run_single_batch( opts.max_steps = max_steps opts.jacobian = self._jacobian opts.timeout_seconds = timeout_seconds - opts.steady_state = steady_state + opts.steady_state = steady_state and not self._early_stop_reads_time() opts.steady_state_tol = steady_state_tol if max_step is not None and max_step > 0.0: opts.max_step_size = max_step @@ -8640,6 +8640,17 @@ def _raise_if_no_steady_state_to_solve_for(self, where: str) -> None: "trajectory to settle, and read the state where it ends." ) + def _early_stop_reads_time(self) -> bool: + """Whether a time course's early stop is held off for the time (issue #961). + + The core holds it off where a rate reads the clock. A species reported + through a rule on the clock, or a compartment it sizes, still moves when + ``‖f(y)‖`` has settled, so the stop is held off here for those too. A + function that only reports leaves it on. + """ + core = self._model._core + return bool(core.functions_use_time) and self._ss_reads_the_time(core) + def _ss_reads_the_time(self, core) -> bool: """Whether the time reaches what a steady state solves for or reports as a state (issues #710, #961). diff --git a/python/tests/test_steady_state_of_a_model_that_reads_time.py b/python/tests/test_steady_state_of_a_model_that_reads_time.py index 60e7b9867..44c86e076 100644 --- a/python/tests/test_steady_state_of_a_model_that_reads_time.py +++ b/python/tests/test_steady_state_of_a_model_that_reads_time.py @@ -394,6 +394,51 @@ def test_a_time_course_of_a_model_that_reads_no_time_stops_early(tmp_path): assert len(species) < 201 and out.solver_stats["steady_state_reached"] +@pytest.mark.parametrize("sensitivity", [False, True], ids=["plain", "with a sensitivity"]) +def test_a_time_course_stops_early_beside_a_function_that_only_reports(tmp_path, sensitivity): + """``seen() = 2*time()`` is read by no rate: A -> B settles as it would + without it, and the run stops there. #710 ran it to the end of its span + (issue #961). The rows it reports read the time where they were taken.""" + kw = {"sensitivity_params": ["k"]} if sensitivity else {} + sim = bngsim.Simulator(_ramped(tmp_path, law="k"), method="ode", **kw) + out, species = _early(sim) + assert len(species) < 201 and out.solver_stats["steady_state_reached"] + np.testing.assert_allclose(species[-1], [0.0, 1.0], atol=1e-6) + seen = np.asarray(out.expressions)[:, list(out.expression_names).index("seen")] + np.testing.assert_allclose(seen, 2 * np.asarray(out.time), rtol=1e-12) + + +def test_a_batch_row_stops_early_beside_a_function_that_only_reports(tmp_path): + rows = bngsim.Simulator(_ramped(tmp_path, law="k"), method="ode").run_batch( + params=[{"k": 1.0}, {"k": 2.0}], t_span=(0.0, 200.0), n_points=201, steady_state=True + ) + for row in rows: + assert len(np.asarray(row.species)) < 201 + + +@pytest.mark.parametrize( + "text", + [ + RULE_ON_TIME, + "compartment C; C := 1 + time/(1+time); substanceOnly species A in C, B in C;" + " A = 1; B = 0; k = 1\nJ: A -> B; k*A\nJ2: B -> A; k*B\n", + ], + ids=["a species set by a rule on the time", "a compartment sized by one"], +) +def test_a_time_course_does_not_stop_while_a_reported_state_reads_the_time(text): + """No rate reads the time, and what is reported as a state does: S, or the + concentrations in C, are still moving where ``‖f(y)‖`` has settled.""" + out, species = _early(bngsim.Simulator(bngsim.Model.from_antimony_string(text), method="ode")) + assert len(species) == 201 and not out.solver_stats["steady_state_reached"] + rows = bngsim.Simulator(bngsim.Model.from_antimony_string(text), method="ode").run_batch( + params=[{"kf" if "kf" in text else "k": 1.0}], + t_span=(0.0, 200.0), + n_points=201, + steady_state=True, + ) + assert all(len(np.asarray(row.species)) == 201 for row in rows) + + def test_a_rate_that_reads_its_own_derivative_is_not_refused(): """Control. ``rateOf`` reads the state's derivatives, which are 0 at a steady state whatever the time: B settles at 1 + rateOf(A) = 1.""" diff --git a/src/cvode_simulator.cpp b/src/cvode_simulator.cpp index ac9286fa1..7dbb330be 100644 --- a/src/cvode_simulator.cpp +++ b/src/cvode_simulator.cpp @@ -4625,7 +4625,10 @@ Result CvodeSimulator::Impl::run_warm(const TimeSpec ×, const SolverOptions // and the run stopped there with the initial state marked steady. Such a // run is integrated to the end of its span, and never marked. (A model // with an event is not on this path.) - const bool check_ss = opts.steady_state && !model.functions_use_time(); + // A function that reads the time and that no rate reads, `minutes := + // time/60`, leaves the test as good as in a model without it, and the stop + // is kept (issue #961). + const bool check_ss = opts.steady_state && !model.rate_laws_read_clock(); const double ss_tol = (opts.steady_state_tol > 0.0) ? opts.steady_state_tol : atol; std::vector ss_derivs; if (check_ss) { @@ -16270,8 +16273,10 @@ Result CvodeSimulator::run(const TimeSpec ×, const SolverOptions &opts) { // to fire: `if(time() > 5, k, 0)` is 0 at the first output point, and the // run stopped there with the initial state marked steady. Such a run is // integrated to the end of its span, and never marked. + // A function that reads the time and that no rate reads leaves the test as + // good as in a model without it, and the stop is kept (issue #961). const bool check_ss = - opts.steady_state && !(model.functions_use_time() || model.n_events() > 0); + opts.steady_state && !(model.rate_laws_read_clock() || model.n_events() > 0); const double ss_tol = (opts.steady_state_tol > 0.0) ? opts.steady_state_tol : atol; std::vector ss_derivs; if (check_ss) {