steady_state, run: solve and stop early beside a function that reads the time and only reports (#961)

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. 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.

docs/user-guide/steady-state.md

@@ -26,15 +26,26 @@ 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
 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

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<int> 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<int> species;
     };

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.

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
@@ -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,52 @@ 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 _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).
+
+        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 +9320,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 +11386,7 @@ class SteadyStateResult:
         "_expression_values",
         "_rule_target_species",
         "_rule_target_expressions",
+        "_expressions_reading_clock",
     )
 
     def __init__(self, core) -> None:
@@ -11501,6 +11539,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 +11735,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):