[Montana]

Fixes #1279.

What was wrong

#1240 (#1228) decided once, at the run's start, whether a rate law's equality over state (if(Bobs==0, kb, 0)) was held, and refused the run if so. That was wrong in both directions.

Silent zeros: the start check missed them. These are #1279's case 1, plus one found while fixing it. dY/dθ, rtol 1e-10:

case	main	central difference	here
SBML washout: event at t ≥ 1 sets B = 0, B made again once A < q	[0, 0, 2.218876], no warning	[−6.4378, −1.0000, 2.2189]	refused where B leaves, t = 3.2189
the same event at t = 0 (t0=false)	[0, 0, 3.218873]	[−6.4378, −1.0000, 3.2189]	refused at t = 3.2189
event at t ≥ te sets B = 0 when B is made at once: on the surface for an instant, at a moving time	dY/dte = −1	−5.8e-6 (0)	refused at t = 4

The last row is the event's jump reading Y's flow on the surface, where B sits for no time. It was wrong on main, and is wrong on #1411's first version too.

Over-refusals: #1279's case 2. Main returned the right columns for every one of these, and #1240 refused them:

case	here	central difference
B never changes (fB() = 0)	[0, 0, 8]	[0, 0, 8]
#1228's model, run ends before t*	[0, 0, 2]	[0, 0, 2]
B starts at 0, nothing makes B, an event sets B = 1 at t = 4	[0, 0, 4]	[0, 0, 4]
if(fB()==0, kb, 0), a rooted condition inside the equality	[−6.437752, −1, 3.218876] (closed form)	[−6.4383, −0.9999, 3.2204]
if(Bobs==0 && Aobs<1, kb, 0), never true as B leaves	[0, 0, 0]	[0, 0, 0]
the same plus kb*Aobs, with B's production a clock switch at t = 3	closed form of kb·∫A	[−36.3359, 0, 19.6340]

#1228's own model over (0, 8), Bobs!=0, and the run from t = 5 are still refused, now at the point where the state leaves the surface.

The change: decide at the exit

• Core: rooting the exit. Each plain rate-law equality over state gets a CVODE root, after the state-switch roots. It reads 0.5 while the residual is exactly 0 and −0.5 otherwise, so it fires where the state leaves the surface or reaches it. An event's restart reads it afresh, so a stretch that only an event ends never fires it, and that event's own jump is right. Only runs with sensitivity columns register these roots.
• Core: judging an exit. At a leaving, judge_held_equality_exit compares the flux of the reactions that read the equality on the surface and just after (compute_flux_split, a gap > 1e-6 of the gross flux).
  • The on-surface state is the last one the root function read exactly 0 at, which bracketing puts within a root tolerance of the exit. CVODE's interpolant doesn't give a step's start back to the bit.
  • No gap: the run continues; nothing restarts.
  • A gap, covered only if the state was put on the surface at the run's start, or at an event whose time no column moves, and held it for no time. "No time" means within the time its species take to move a few ulps: A decaying from 10 reads Aobs==10 until t = 1.5e-16.
  • Otherwise: SensitivityUnsupportedError, naming the equality and the time.
• A state switch written inside the equality (if(fB()==0, …) with fB gated by q-Aobs>0): its jump reads the equality's readers either side of its own surface. A leaving at the instant it is crossed, in the same root return or right after its restart, is covered. Python passes these switch conditions per equality (set_held_equality_inner).
  • Whether that switch's restart lands on the equality's surface, so the leaving is rooted and judged, or lands past it, is rounding. A rate law that branches on an equality the state sits on (if(B==0, ...) with B held at 0) runs on the analytic sensitivity path and returns 0 for the columns that move the switch #1228's parameters landed past it on macOS and on it on Linux, where CI's first run of this design refused the case. A sweep of 186 parameter sets on macOS took the judged path in 99.
  • The test now runs 4 parameter sets, 3 of which take that path on macOS. Removing the coverage, in the core or in Python, fails those 3.
• Removed:
  • equalities_held_at_start and the start-time refusal;
  • sensitivity: root where the state leaves a rate-law equality, and decide there (#1279) #1411's first-version check at events;
  • the corpus filter that version needed (a surface only events move doesn't root).
• Kept from the first version:
  • the typed bngsim::SensitivityUnsupported, mapped to SensitivityUnsupportedError;
  • compute_all_sensitivities passing that refusal through instead of wrapping it in SimulationError.
• Python: _apply_held_equalities passes the equalities (rate_law_equalities), their reactions, and the switches written inside them (set_held_equality_conditions, set_held_equality_reactions, set_held_equality_inner).
• Docs and changelog: docs/user-guide/sensitivities.md describes the refusal as it now stands, and changelog.d/1279.fixed.md is added.

The diff (src/cvode_simulator.cpp, abridged)

@@ src/cvode_simulator.cpp — cvode_event_root_fn: a root where the residual leaves exactly 0
+    if (data->held_equality_roots != nullptr) {
+        [...]
+        for (int j = 0; j < static_cast<int>(held.size()); ++j) {
+            const bool on = eval.evaluate(held[ju]) == 0.0;
+            gout[base + j] = on ? 0.5 : -0.5;
+            if (on && data->held_on_t != nullptr) {  // the last state read on the surface
+                (*data->held_on_t)[ju] = static_cast<double>(t);
+                std::copy(y_ptr, y_ptr + nsp, [...]);
+            }
+        }
+    }
@@ judge_held_equality_exit — at a leaving
+        bool changed = !on_surface;
+        if (on_surface && !held.reactions.empty()) {
+            model.compute_flux_split(t_b, x_b, &held.reactions, net_b.data(), gross_b.data());
+            model.compute_flux_split(t_x, x_after, &held.reactions, net_a.data(), gross_a.data());
+            for (std::size_t i = 0; i < net_b.size() && !changed; ++i) {
+                changed = std::fabs(net_a[i] - net_b[i]) > 1e-6 * (gross_a[i] + gross_b[i]);
+            }
+        }
+        [...]  // at_once: the time the residual's species take to move a few ulps
+        bool covered =
+            !changed || held_since(t_start) ||
+            (std::isfinite(last_event_t) && !last_event_moved && held_since(last_event_t));
+        [...]  // or a state switch written inside the equality, crossed at this instant
+        if (!covered) {
+            throw SensitivityUnsupported("[...] leaves it at t=[...], where the rate laws "
+                                         "that read it change. [...] (issue #1228; issue #1279)");
+        }
@@ the root return: no crossing, restarts nothing
+            if (flag == CV_ROOT_RETURN && n_held_equality > 0) {
+                [...]
+                    equality_left.push_back(eq_base + k);
+                    impl_->judge_held_equality_exit(ns, k, t_ret, y_data, [...]);
+            }
-            bool only_held_left = !held_left.empty();
+            bool only_held_left = !held_left.empty() || !equality_left.empty();
@@ python/bngsim/_simulator.py — _apply_state_switch_sens
-            held = equalities_held_at_start(core, ctx, t0=t_start)
-        [...]
-        if held:
-            raise SensitivityUnsupportedError("[...] the state starts on that equality's surface [...]")
+        self._apply_held_equalities(opts, core, ctx, conditions)   # conditions, reactions, inner switches

Tests

test_equality_held_over_state.py (#1228's file, updated):

• the 5 lifted .net stretches above, at rtol 1e-6 against the closed form;
• the inner-switch stretch at 4 parameter sets, against the closed form [−ln(10/q)/a², −1/(a·q), ln(10/q)/a];
• every former refusal (Bobs==0, Bobs!=0, the run from t = 5, run_batch, compute_all_sensitivities) is still refused, now at the exit;
• the test of the removed start check is gone.

test_equality_landed_on_by_event.py:

• refused: the washout (mid-run and at t = 0), and the moving-time event;
• run: a stretch an event starts that outlasts the run, and one only an event ends;
• controls held to closed-form columns: an event off the surface; one on it for an instant at a fixed time; one that writes what the equality doesn't read; a flag only events move, as a parameter and as a boundary species;
• every entry point raises the typed refusal.

Mutation tests (the core rebuilt for each C++ mutant): every one is killed except the last.

mutant	killed by
judge never refuses	the refusal tests
no flux-gap test	the conjunction rows
no at-once at the start	Aobs==10 / Aobs!=10
no at-once at a fixed-time event	"on it for an instant"
event time moving ignored	the moving-time test
no ulp window (time tolerance only)	Aobs==10
on-surface state not recorded	the controls
event batch not recorded	"on it for an instant"
equalities not passed (Python)	11 tests
reactions not passed (Python; all reactions read instead)	the clock-switch row, added for this
chunk wraps the refusal	compute_all_sensitivities
no inner-switch coverage (core), or inner conditions not passed (Python)	3 of the 4 inner-switch parameter sets
a fixed species counted as moved, from the first version	equivalent; that code is gone with the filter

No change where the state never sits on a surface

• Synthetic: 4 models whose equality is rooted but never held (B only crosses 0, Bobs!=0 the same, Aobs==1000, and an equality only an output reads). Species and sensitivities are bit-identical to main, with the same step counts.
• Corpus: 3 of 1,323 BioModels have a rate-law equality.
  • BIOMD0000000301 and MODEL2003190004: bit-identical to main, with 8 parameters, at t_end = 20 and 200.
  • BIOMD0000000446: refused for an unrelated reason (no analytic sensitivity right-hand side), on main and here alike.
• Plain runs: nothing registers.

Checks (macOS arm64)

• python/tests: 16,871 passed, 0 failed.
• parity_checks/tests: 1,015 passed. The 2 failures (test_bngsim_provenance and the reference-archive replay) also fail on main here: their subprocess can't import bngsim without my local import shim.
• ci/changelog.py required and check pass: the first push lacked the fragment, which is what turned "fragment staged" red.

pre-commit run --all-files and sphinx -W are clean.

Not verified

• Linux: the at-once window and the exact-zero readings are rounding-level. CI's Linux legs are the check.
• Still over-refused: a leaving at once after any other restart. One example is B made by a literal clock switch, fB() = if(time()>5, kp, 0), while Y reads Bobs==0. It is refused, but central differences give the exact [0, 0, 5]. Main refused it too, at the start.
  • The judge accepts "held for no time" only after the run's start or an event, since only those restarts record whether their time moves with a column.
  • Lifting it means each crossing-stop restart site recording the same flag. That is a follow-up.