[Montana]

Fixes #1340.

What was wrong

MODEL1812040005's species P is a rate-rule target at an exact equilibrium of 0, read by a_hat = a*P^z with z fractional. Run in the issue's window, with the whole state compared by name:

order	T, I1, I2, V agree to	P
as written	—	exactly 0 throughout
reversed	1e-13 of each species' peak	−3.6e-27 at t ≈ 0.15, so P^z is NaN and the run fails with CV_REPTD_RHSFUNC_ERR

The right-hand side is the same function in both orders. Summed in another order, its derivative of P at that equilibrium leaves a rounding residue a hair below 0 instead of exactly 0. The original order succeeds by an accident of rounding.

When the RHS goes non-finite, it is retried on a nonnegative-clamped state (GH #135). That clamp leaves alone every slot that may be negative by design (#706): a rate-rule target, such as a membrane potential at −61, a fixed or promoted slot, or one declared negative. So P was left at −3.6e-27, and the retry failed too.

The issue's hunch was right: the clamp mask's indexing (product_indices, 1-based) is fine. The mask does what it should; the gap is that "negative by design" also caught a value that is numerically zero.

The fix (cvode_simulator.cpp, clamp_state_nonneg)

On that retry only, a slot the clamp leaves alone is read as 0 where it is negative by less than its own absolute tolerance. One negative by more is left alone, and still fails where its rate law is undefined. This is the rule the sensitivity RHS already applies (GH #392, clamp_state_numerically_zero), now on the value RHS too. The tolerance is the run's own per-species atol (rhs_zero_atol), filled on both run paths: create_cvode_core and the warm cache.

-    for (int i = 0; i < ns; ++i)
-        s[i] = (y_ptr[i] > 0.0 || !mask[i]) ? y_ptr[i] : 0.0;
+    for (int i = 0; i < ns; ++i) {
+        const double v = y_ptr[i];
+        if (v > 0.0) {
+            s[i] = v;
+        } else if (mask[i]) {
+            s[i] = 0.0;
+        } else {
+            s[i] = (have_tol && v > -tol[i]) ? 0.0 : v;   // numerically zero (#1340)
+        }
+    }

The retry engages only after the unclamped RHS is non-finite. A model whose RHS stays finite is byte-identical, and the original order of MODEL1812040005 is one.

Verified

MODEL1812040005 in both orders now runs. Each species agrees across the orders to about 1e-12 of its peak, and P to 2e-16 absolute, against an atol of 1e-9.

Corpus A/B. Every one of the 1,300 BioModels ODE jobs, run on main's core and on this one:

	main	this PR
ran	1,246	1,247
ran on both: final state bit-identical	—	1,246 of 1,246
outcome changed	—	1: MODEL9811206584, fails on main, runs here

MODEL9811206584 is the same mechanism through a different kind of excluded slot: it has no rate rules. On main it fails at t = 95.8 on log(s) with s = −1.2e-12, inside atol (phi_b_s, phi_r_s NaN). Here s is read as 0 on the retry, so phi_b_s goes to 0, its limit as s → 0 from above. Its other negative species, r = −2.6e-8, is past atol and is not snapped. Over [0, 90] the new run agrees with RoadRunner to 1.4e-5 (rtol 1e-6), and with itself at 100× tighter tolerances to 2.4e-6. RoadRunner itself fails on this model at t = 91 (CV_TOO_MUCH_WORK), so past that neither engine is a reference.

Tests (python/tests/test_rate_rule_numerically_zero.py)

A deterministic model of the mechanism: a rate-rule target P' = -c from P = 0, read by P^0.5 in A's decay.

test	on main	here
c = 1e-30 (P sinks to −1e-29, inside atol), both declaration orders: runs, A stays 1, P itself unmoved	fail ×2	pass
the same, both orders: bit-identical	fail	pass
a run with sensitivity columns, which sets up on the other path, create_cvode_core (P^1.5, flat at 0)	fail	pass
control: c = 1e-3 (P to −1e-2, past atol), both orders: still fails	pass ×2	pass
MODEL1812040005 in either order agrees (skips without the corpus, "test data not found")	fail	pass

Mutations, each rebuilt and caught:

• no tolerance limit (snap every excluded negative): the control fails;
• the cold path's tolerance unfilled: the sensitivity test fails. It first survived, because the plain runs take the warm path; the sensitivity test is there for it;
• the warm path's unfilled: 5 fail.

Suites: python/tests and parity_checks/tests give 17,329 passed, 172 skipped, 1 xfailed. clang-format 17.0.6 and ruff 0.14.0 are clean.

[elyfmiller]

A heads-up before we review this one, which is ours this week: it conflicts with main since #1354 merged (c0877fa). The conflict is in the loop of clamp_state_nonneg that builds the clamped state.

• Refuse a run carried to a pole instead of clamping it there (BIOMD0000000338) #1354 rewrote that loop to leave a slot alone at a pole the run landed on (spare).
• This PR rewrites the same loop to read as 0 a slot negative by less than its atol.

Could you rebase it on main? The two rules need combining in that loop. #1370, which is on current main, conflicts with this branch in the same place.

We'll review the three PRs that change this clamp (#1300, #1346 and #1370) together, once all three are on current main.

[elyfmiller]

Thanks. We reviewed this independently, Opus 5.5 and Fable 5.1, each on our own build of d787eac (a fast-forward on main e425546; main has since moved to edf8876, with no change to src/ or include/, and the merge is clean), macOS arm64, KLU, MIR and LAPACK on.

For the issue's case the fix does what it says:

• MODEL1812040005 reversed fails on main and runs here, in both orders;
• a rate-rule target past its atol still fails;
• the sensitivity path reads the band too.

But the band is applied whenever the retry runs, to every excluded slot inside it, and that has two effects the tests don't see. One moves a slot that didn't make f non-finite and silently changes the answer (§1). The other turns a fast failure into a run that creeps until the wall-clock timeout (§2). What we checked comes first, so it needn't be redone.

What we checked

check	result
Fail before: the new file on main's core	4 failed, 2 passed, 1 skipped, under cc and MIR: the two inside-atol cases, the order case and the sensitivity case fail with CV_REPTD_RHSFUNC_ERR; the two past-atol controls pass. As your table.
MODEL1812040005 itself (from a local BioModels copy, md5 93a57e68)	Reversed: fails on main (CV_REPTD_RHSFUNC_ERR at t = 0.16, P = -1.4e-25), runs here. test_model1812040005_runs_in_either_order with the file at its path: fails on main, passes here.
After	6 passed, 1 skipped (no corpus) under cc, BNGSIM_CODEGEN_JIT=mir, BNGSIM_CODEGEN_THRESHOLD=1, and both; 7 passed with the corpus file.
Full python/tests, BNGSIM_SKIP_AUDIT=strict	16,420 passed, 167 skipped, 0 failed
parity_checks/tests	933 passed, 33 skipped
Corpus A/B: every local BioModels SBML (1,641), run(t_span=(0,1), n_points=101, rtol=1e-6, atol=1e-9), main's core against this PR	1,489 ran on both. 1,488 bit-identical; 1 differs (low 6). The one status change, MODEL1001200000 (timeout → SimulationError), was our 60 s alarm under load: given 400 s, main fails the same way.
Mutations, rebuilt, against the new file and the clamp-rule files	Caught: no atol limit (8 fail: your controls and the existing The GH #135 non-finite-RHS retry sets every negative state component to 0, not just the offending concentration. Legitimately negative rate-rule parameters (e.g. a membrane potential) and negative boundary species then follow the wrong... #706 …keeps_its_value and …follows_its_own_ode cases); the cold path's atol unfilled (the sensitivity case); the warm path's unfilled (3); a spared masked slot read at 0 (the three test_338_carried_to_its_pole… cases). Not caught: the band read from clamp_atol (the same numbers; low 4); tol[0] for tol[i], a band ten times wider, or a closed band >= (low 5).
Read	clamp_state_nonneg's callers: the RHS and codegen RHS retries, the four closed-form Jacobian fills, spare_at_pole, followed_below_zero, so the band reaches the Jacobian retries too; spare marks only masked slots, so the new branch never meets one; the warm cache is keyed on rtol, atol and atol_vec, so a run with a changed atol rebuilds and the band can't go stale; GH The zero-base exponent guard tests base == 0, and the solver hands it a base 24 orders below atol on the wrong side #392's clamp_state_numerically_zero, the same rule on the sensitivity RHS.
Gates	ruff check and format, clang-format 17.0.6 and ci/changelog.py check all clean
Overlap	Trial merges of the other 37 open PRs, on edf8876, onto main and onto main + this PR: it adds no conflict. Two touch the clamp: #1300 calls clamp_state_nonneg from the sensitivity RHS, so it inherits the band; #1370's run-test clamp is a copy that doesn't (§3).
CI	23 checks, all green. No leg that runs python/tests has the rr_parity corpus (low 7).

1. The band moves a slot that didn't make f non-finite

clamp_state_nonneg applies the band in the same pass as the masked clamp. So wherever any slot makes f non-finite, every excluded slot inside its atol is read as 0, even where the masked clamp alone already makes f finite. That covers a rate-rule target, a fixed or promoted slot, and a species declared negative. Main reads such a slot at its value, which is what the run has.

A model built for it, in Antimony:

species S = 1, A = 1, B = 0;
V = -5e-10; V' = 0;
c = 1; k = 0.1; g = 1e9;
J1: S -> ; c;
J2: A -> ; k*A*S^0.5;
J3: -> B; -g*V;

S is drained past 0 and read by S^0.5. It's masked, and its clamp alone clears the NaN. V is held at -5e-10 by design (a rate-rule target, declared negative), inside atol 1e-9. B grows at -g·V = 0.5, so exactly B(t) = 0.5·t. With run(t_span=(0, 10), n_points=11, rtol=1e-8, atol=1e-9), under auto and fd:

build	B at t = 0, 1, 2, 5, 10
main	0, 0.5, 1.0, 2.5, 5.0
this PR	0, 0.5, 0.5, 0.5, 0.5

V itself stays at -5e-10; every retry after S < 0 reads it as 0. Each of us built the same fix: a two-stage retry, evaluating at the masked-only clamp first and applying the band only if f is still non-finite. On the interpreted RHS, it passes the nine clamp-rule files (134 passed, 1 skipped, as the baseline) and gives B(10) = 5.0. The same staging would carry to the codegen RHS and the Jacobian retries. It doesn't address §2, where the slot in the band is the cause.

2. A slot that drains past -atol makes the run creep instead of fail

Inside the band the retry reads the slot as 0 and the step succeeds. An ulp outside, f is NaN, CVODE cuts h, and the cut step's prediction lands back inside the band and succeeds. h settles near 1e-16 to 1e-19. Each batch still moves t by a few ulps, which retry_while_advancing's stall test (!(t_after > window_from)) counts as progress. Nothing ends the run but the wall-clock timeout.

The test file's own shape, P' = -c from 0, read by k*A*P^0.5, with run(t_span=(0, t_end), rtol=1e-8):

case	P crosses -atol at	this PR	main
c = 1e-10, atol 1e-12	t = 0.01	still running at 150 s with max_steps 500 or 2000 (t_end 1) and 20000 (t_end 0.02); with a 20 s timeout (t_end 1), SimulationTimeout	fails in under 1 s
c = 1e-10, atol 1e-9, t_end 20	t = 10	still running at 150 s; SimulationTimeout at 20 s	fails in 0.0 s
c = 1e-6, atol 1e-9, t_end 10	t = 1e-3	"CVODE made no progress" in 0.0 s	fails in 0.0 s

So the outcome depends on where the crossing lands. The band is meant for rounding residue, which never leaves it. A slot genuinely draining past -atol is the "negative by more" case, which should still fail. Here it fails only by timeout, and only if one is set. Two ways out:

• a stall test that also counts h below CVODE's own "t + h = t" threshold (uround·|t|) as no progress, ending with the non-finite witness named;
• or a retry that refuses a slot it snapped on the previous step once that slot has gone past the band.

3. The steady-state run test

#1370's run test clamps with its own copy, clamp_nonneg_into, without this band. We built main + #1370 + this PR. On this PR's model shape, run() reaches 1e6 while steady_state(method="integration") is refused (residual 5e-31, run measure inf). It's in our #1370 comment, §2. Whichever lands second should carry the band to the run test, or say why it differs.

Low

4. rhs_zero_atol is clamp_atol.
   • record_clamp_tolerances, on the line after each new assignment, already stores the same per-slot atol (atol_v or the scalar), and the band read from it passes everything.
   • sens_zero_atol holds the same numbers too, and the band repeats clamp_state_numerically_zero's rule.
   • One per-slot atol and one "numerically zero" helper would keep the value and sensitivity RHS on the same band by construction.
5. The band's width and per-slot reading are unpinned.
   • Every case uses a scalar atol, with the inside cases at -1e-29 and the controls at -1e-2 against atol 1e-9. So any band between about 1e-28 and 1e-2 passes, as does reading slot 0's atol for every slot.
   • A vector-atol case pins the reading: with P at -5e-10, atol=[1e-12, 1e-6] runs. [1e-6, 1e-12] doesn't, but today that side is §2's creep rather than a failure.
   • A control a few atols out would pin the width.
6. The corpus claim is 1,488 of 1,489, not all. MODEL1210260004 ends 1.4e-10 apart on both our builds (5 steps / 10 RHS evaluations on main, 7 / 22 here; two slots end at -1e-42 and -2e-34). It's within tolerance and deterministic, and since the change acts only on the retry, the run met it. The body says none changed.
7. Two cases no CI run sees.
   • test_model1812040005_runs_in_either_order needs the rr_parity corpus, which no leg running python/tests has, so it skips on every PR. A copy of the model (5 species) under tests/data would run it, and keep it out of test_sens_default_columns' corpus sweep, which refuses a directory of one model.
   • Every new case puts a rate-rule target in the band. None has a fixed, promoted or declared-negative slot there, though the one corpus outcome you report, MODEL9811206584, has no rate rules.
8. The field comment and the changelog list "a rate-rule target, a fixed or promoted slot". The band also reaches a species declared negative, which is §1's V.

Not checked:

• Linux and Windows beyond CI.
• The two-stage retry on the codegen RHS and the Jacobian retries.
• sens: evaluate the sensitivity RHS at the value RHS's clamped state (#1298) #1300 with this PR on one build (sens: evaluate the sensitivity RHS at the value RHS's clamped state (#1298) #1300 is on its old base).
• Your RoadRunner comparison on MODEL9811206584.