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

[elyfmiller]

Following up on our review above: we've merged main in and pushed the review items on top, as c4239d2, 688da77 and 532985e, and we're merging now that CI is green on them. The band itself, read at each slot's own atol, both run paths, and the sensitivity case stand as you wrote them.

What c4239d2, 688da77 and 532985e change

• The merge (c4239d2). steady state: the run test retries a NaN f at the clamped state, as the time course does (#1357) #1370 moved the clamp's loop into clamp_nonneg_into (include/bngsim/nonneg_clamp.hpp), which the steady-state run test shares; that was the conflict. The merge keeps main's call and reads your band over its result, the same at every slot. Its diff against main is this PR's three files.
• §1: the band is a second reading. read_near_zero_into, in the shared header, puts at 0 each slot the mask leaves alone that is below 0 by less than its atol. At all six retries (the RHS and the dense and sparse Jacobians, interpreted and compiled), it's read only where the masked clamp's reading is still non-finite. On §1's model, B(10) is 5.0 again under auto and fd, as on main and in the exact answer.
• §2: a slot drained through the band is refused where it leaves it. drained_through_band, beside followed_below_zero, ends the run where three things hold: the step started at the band's edge (within a relative 2^-20 of -atol), the slot is past the band now, and the reading turns finite with it at 0 too. The message names the slot, where the step started it and where it is, and the time. The three creep cases are refused at t = atol/c, in under a second each. It's the edge, not anywhere in the band, because a step tried too long for the band is cut and the shorter one accepted. Refused from anywhere inside, the t_end = 20 case went at t = 3.8, on a step tried out to t = 38, though P leaves the band at t = 10.
• §3: the run test reads the band. Its f and both closed-form Jacobian fills take the same second reading, at the march's own atol per slot. On your model's shape, steady_state(method="integration") now converges in both orders, under auto and fd, where it was refused (run measure inf). The run test is taken over a budget of steps (take_on_for), so it needs no drain refusal of its own.
• The closed-form Jacobians (532985e). Every other case reads P through P^0.5, whose slope at 0 is infinite under any reading, so the Jacobian falls back to the difference quotient and none of the four fills' second readings was reached. P^1.5 under jacobian="analytical" reaches each one: dense and sparse, interpreted and compiled. The sparse cases carry 60 inert species, since codegen emits its sparse fill only past the sparse gate.
• Low 4. The band reads clamp_atol, which record_clamp_tolerances already fills on both run paths, and rhs_zero_atol is gone. The sensitivity RHS keeps its own sens_zero_atol, which its opt-out clears.
• Low 5. A vector-atol case reads P against its own atol: inside at 1e-6, past at 1e-12, with A's the other way round. A width case puts P at half an atol (inside), one atol and three (past).
• Low 7. MODEL1812040005 is copied to tests/data (16 KB, beside MODEL1812040006), so its order test runs on every leg. It's outside the rr_parity directory, so test_sens_default_columns' sweep doesn't see it. A species declared below 0 and a fixed one each get a case in the band.
• Lows 6 and 8. The changelog and the comments name a species declared negative. The fragment no longer says every corpus model runs bit-identical. Read in one pass, the band moved MODEL1210260004 by 1.4e-10. Read second, it's bit-identical again, as the corpus row below shows.

One change on main matters here. Since #1434 (#1420), a scalar atol is a ceiling per species, min(atol, rtol·max(|x_i(0)|, floor)). So a slot that starts at -5e-10 is held to about 5e-18, past its own band, and is left alone, as a slot negative by design should be. The new cases that hold a slot below 0 pass bngsim.ExactAtol(1e-9) or a vector.

What we checked

check	result
The merge	c4239d2 against main 3616b39 is this PR's three files, with your changes. Our two independent merges read the same at every slot. On it, the nine clamp-rule files and test_steady_state_scaled_criterion.py: 362 passed, 3 skipped.
Fail before: the new file on the merge's core	9 failed, 17 passed, under cc and MIR: §1 (auto, fd), the three drains (two run to the 20 s timeout; the third stops with "CVODE made no progress" at the crossing, without naming P), and the run test (4). The 17 that pass there are your six, MODEL1812040005, the six pins on width, per-slot atol and the declared-negative and fixed slots, and the four Jacobian cases, which your one-pass band reaches too. The mutations below check those.
After	26 passed under default, BNGSIM_CODEGEN_JIT=mir, BNGSIM_CODEGEN_THRESHOLD=1 and both. With the clamp files and the run-test file: 382 passed, 2 skipped.
Mutations, each rebuilt	13 of 14 caught: the band read in one pass with the clamp: the V case (2); no second reading on the interpreted RHS: 11 (under BNGSIM_CODEGEN_THRESHOLD=1 the compiled RHS takes them); on the compiled RHS: 12 under that threshold, and the sensitivity case by default; none at each of the four Jacobian fills: that fill's own case; a drain refused from anywhere in the band: the three drains, on their times; no drain refusal: the three, at their 20 s timeouts; the band read from slot 0's atol: the per-slot case; a closed band: the one-atol case; a band ten atols wide: the one- and three-atol cases and the three drains; the run test without its second reading: its four. Not caught: a drain refused without the reading at 0 that confirms the slot is the cause, since no case has a second cause of a non-finite reading beside the slot.
Full python/tests, BNGSIM_SKIP_AUDIT=strict	On main 3616b39 + 532985e: 17,130 passed, 435 skipped, 1 xfailed (#1334's RNG pin), 1 failed: test_run_tests_sh's walk-up probe, since that venv lacks jax, diffrax and h5py and their modules skip. On main 3616b39 + 688da77 in a complete environment: 17,564 passed, 144 skipped, 1 xfailed, 0 failed.
Corpus A/B: every local BioModels SBML (1,641), run(t_span=(0, 1), n_points=101, rtol=1e-6, atol=1e-9, timeout=120), main's core against 688da77's	1,489 ran on both, and all 1,489 are bit-identical, at the same steps and RHS evaluations. No status changes, and no model reaches the band or the drain refusal. MODEL1210260004, which the one-pass band moved (5 steps and 10 evaluations on main, 7 and 22 there), is back at main's.
Gates	clang-format 17.0.6, ruff 0.14.0 check and format, ci/changelog.py check: clean. git diff --check flags only the copied SBML's own trailing whitespace (9 lines), which the file has as BioModels ships it, as MODEL1812040006 beside it does.
Overlap	We trial-merged the other 12 open PRs onto main fcf07d1 + 532985e. None conflicts with this PR; #1382 conflicts with main on its own, the same with or without it. Four share a file and merge clean: #1300 and #1398 (cvode_simulator.cpp), #1424 (steady_state.cpp), and #1412 (both). #1300 calls clamp_state_nonneg from the sensitivity RHS, so it gets the masked clamp alone; we'll take that up on #1300. #1411, merged today, also changes cvode_simulator.cpp (the event root and the sensitivity jump, not the retries): on main 047c12a + 532985e, built, the new file, the clamp files and #1411's two files give 422 passed, 2 skipped under cc and MIR. #1421, merged since, changes none of this PR's files. 532985e merges clean onto main fcf07d1.
CI	23 checks passed on 532985e, none failed, the MIR codegen legs among them. The suite counts: 17,758 collected.