[Montana]

Closes #1357.

What was wrong

steady_state(method="integration") on glioblastoma_egfrviii_seed.net, compiled (codegen=True), on main:

platform	converged	residual	continued_run_distance	where the run test stopped
arm64 macOS	False	2.1e-39	inf	t = 3.9e5
x86-64 Linux	False	1.6e-40	inf	t = 0, its first step

The interpreted solve converges on both platforms.

What it was

The model's v_cancer_drive raises Total_pAKT, which is AKT(s~P), to n_hill = 3.5. A fractional power of any negative number is NaN.

On x86-64 Linux, the compiled march ends with AKT(sP) at exactly 0. On macOS it ends within 1e-52 of it. The scaled criterion's run test (#1001) then takes that state on for 1e6. The corrector's first iterate puts AKT(sP) an ulp below 0, f is NaN, and CVODE fails on any step size. I measured this on main's core:

• Model.rhs is finite at AKT(s~P) = 0 and 1e-60. It is NaN at -1e-60 and at -1e-300.
• From the Linux state, the march and the run test fail at t = 0 on the closed form and on difference quotients.
• Simulator.run from the same state, at rtol and atol of 1e-8 and of 1e-13, runs to t = 1e6 on every Jacobian and moves nothing, by 8e-38 at most. The time course retries a non-finite f at the nonnegative-clamped state (GH Extract the two forward-sensitivity jump handlers out of CvodeSimulator::run (follow-up to #109) #135). The steady-state march has never had that retry.

So the run test refused a steady state over a NaN that only its own integrator met.

The first version of this PR was wrong

It retried the run test on difference quotients. That passed on macOS, where the FD run happened to step around the NaN. It failed on Linux CI, where from the exact-0 state the FD run also fails at t = 0. I reproduced the Linux failure in an x86-64 container, and on macOS by starting the solve from the Linux state's bits. That retry is removed.

The fix

src/steady_state.cpp. The run test now integrates as the time course does. Where f is NaN, with no infinity, it is taken again at the nonnegative-clamped state, and only then. A model whose f is finite runs as before, to the bit.

An infinite f is not retried. Since this branch was opened, main's time course spares a slot it reads at a pole (spare_at_pole). There, BIOMD0000000338's XIIa and K were clamped into a wrong answer, -1. A NaN is still clamped at any depth. So the run test takes only the NaN case, and a run test that meets an infinity is refused, as on main.

+static void ss_eval_rhs(SteadyStateUserData *data, double t, const double *y, double *ydot) {
+    data->rhs->eval(t, y, ydot);
+    const auto un = static_cast<std::size_t>(data->model->n_species());
+    if (data->clamp_nonfinite && ss_nan_and_no_inf(ydot, un))
+        data->rhs->eval(t, ss_clamped_state(data, y), ydot);
+}
 ...
         SteadyStateMarcher run(model, rhs, run_opts, sub);
         run.allow_error_test_failures(kHoldErrorTestFailures);
+        if (end == nullptr)
+            run.clamp_nonfinite_as_the_time_course();

cvode_ss_rhs and the colored difference quotient (KLU) both go through ss_eval_rhs, so the quotient is formed from the same f as fy. The dense quotient calls cvode_ss_rhs, so it gets the retry with no extra change.

The closed-form Jacobians, dense and sparse, get the same NaN-only retry in the run test. This is the Jacobian half of GH #135. CVODE evaluates the Jacobian at its prediction. At the Linux state, AKT(s~P) is exactly 0 with f = -1.1e-39, so the first prediction puts it near -1e-43, and d/dx x^3.5 = 3.5·x^2.5 is NaN there.

include/bngsim/nonneg_clamp.hpp (new). This holds the rule for which slots may be clamped (#706: not fixed, reported, not declared negative, not a rate-rule target). It moved verbatim out of cvode_simulator.cpp, whose nonneg_clamp_mask now calls it. Main's clamp_state_nonneg, with its spare slots, is untouched. The time course is unchanged.

What stays unclamped, and why:

• The march. I tried clamping the whole march first (mutant M1 below). It accepts states with AKT(s~P) a hair below 0, where the unclamped residual is NaN. Four tests fail, including two existing march tests.

• The take-on run (ss_continued_run_once(..., &end)), whose end state the solve adopts. That run belongs to the march, so it stays as the march is. Clamping it too (M5) changes no test and no corpus result. The exclusion is conservative, not forced.

• Note on the Jacobian retry. The previous push left it out, because no macOS test needed it (M3 survived). The Linux shards then failed test_a_compiled_march_to_a_fractional_power_at_zero_is_not_refused: from the Linux state, the closed-form run test still failed at t = 0. The three pinned-state tests passed there. I could not reach that run test from Python on macOS: from the same state, the march's own first step fails and the polish moves the state. So I used a temporary build hook that runs the run test from the solve's initial state, and measured it from the Linux bits:

  Jacobian	without the Jacobian retry	with it
  closed form, interpreted	inf	1.4e-16
  closed form, compiled	inf	1.4e-16
  fd, interpreted	1.3e-8	1.3e-8
  fd, compiled	2.0e-10	2.0e-10

  The hook is not in the PR. The compiled test is the regression test for this half, and it fails on x86-64 Linux without it.

Tests (test_steady_state_scaled_criterion.py)

• test_a_compiled_march_to_a_fractional_power_at_zero_is_not_refused: glioblastoma with codegen=True converges, with run measure ≤ 1. The source is codegen, solver_jacobian_retried is False, and the concentrations match the interpreted solve's at the march's tolerances, 1e-8.
• test_a_run_test_from_a_fractional_power_at_zero_is_taken_to_its_time[auto | fd | fd colored]: a solve started from the Linux state's exact bits. It converges with run measure ≤ 1, at that state to 1e-8. fd colored forces KLU and asserts it. This test needs no Linux: on main it fails on macOS too.

On main's core: three of the four fail, each with converged=False. [auto] passes there; it is a check, not a regression test.

The previous control, test_an_explicit_closed_form_is_not_called_off_in_the_run_test, asserted that jacobian="analytical" is still refused. The cause was never the Jacobian, so that solve now converges too, and the test is removed.

Mutations (C++ rebuilt for each; all of test_steady_state_scaled_criterion.py)

mutant	result
M1: clamp the whole march, not just the run test	4 fail: both new compiled/fd tests, plus ..._returns_the_state_it_passed[integration] and ..._still_retried_on_difference_quotients
M2: no retry in ss_eval_rhs	3 fail: the compiled test, [fd] and [fd colored]
M4: the run test never switches it on	3 fail: the same
M6: colored quotient evaluates f unclamped	1 fails: [fd colored]
M3: no retry in the dense closed-form Jacobian	passes on macOS, which never reaches it (above). On x86-64 Linux, the compiled test fails: that was this branch's previous push
M5: also clamp the take-on run	227 pass. Excluded on purpose (above)
M7: retry an infinite f too, after the merge of main	238 pass, with test_clamp_far_below_zero.py. No test or corpus model meets an infinity in the run test; NaN-only follows main's time course (above)

Measured

Corpus. Every .net model in the repo (177, from tests/data, benchmarks and parity_checks), steady_state() interpreted and compiled. Before the merge of main, against main's core:

Jacobian	identical to main	changed
default	353 of 354	glioblastoma compiled: refused (crd inf) → converged, at the same residual 2.15e-39
fd	353 of 354	glioblastoma interpreted: converged before and after, same state and residual; the run measure goes 0.022 → 0.57, since the run now continues through the NaN on clamped f

M5's sweep is identical to this PR's on all 354. On the final code (main merged, Jacobian retry in), against the same tree with the run test's switch off (M4, main's steady-state behaviour): default Jacobian, 353 identical, and glioblastoma compiled goes refused → converged at 2.15e-39. fd, 353 identical, and glioblastoma interpreted has the same state and residual, with run measure 0.022 → 0.57, as above.

With #1350. Under BNGSIM_CODEGEN_THRESHOLD=1:

• test_a_march_whose_integrator_fails_returns_the_state_it_passed[integration], which tests: declare what a forced-JIT run cannot pass, so the run is a signal (#1217) #1350 marks as a strict xfail citing steady_state(method='integration') on the compiled RHS and Jacobian refuses a state at residual 2e-39 and does not retry on difference quotients #1357, now passes. Whichever PR merges second drops that mark.
• test_a_march_whose_integrator_fails_is_still_retried_on_difference_quotients (MODEL1712050008) fails identically on main and here: converged=True, residual 1.05e-15, solver_jacobian_retried=False. It pins the arm64 interpreter's route, not this bug, and tests: declare what a forced-JIT run cannot pass, so the run is a signal (#1217) #1350 moves it to its threshold_forced skips.

Full suite on the final code (main c0877fa merged): 16,419 passed, 0 failed, 138 skipped.

Linux. I reproduced the original failure in an x86-64 Linux container. Docker would not start for the later versions, so this PR's Linux CI shards are the x86-64 check. On the previous push they passed the three pinned-state cases and failed only the compiled one, for the Jacobian reason above.

[elyfmiller]

Thanks. We reviewed this independently, Opus 5.5 and Fable 5.1, each on our own build of 7bb8fb1 (on main c0877fa; the merge is a fast-forward, tree 3584813b), macOS arm64, KLU and MIR on. Main has since moved to e425546 (#1315, #1325), which touches none of this PR's files; the merge is clean. The fix does what it says:

• with the clamp on for the run test alone, glioblastoma's compiled steady state is accepted at its residual of 2.15e-39;
• a solve from the Linux state's bits is accepted under fd and under fd with KLU;
• a run test whose f stays finite runs as before, which the corpus shows bit for bit.

We found no defect in the diff. There is one pre-existing gap the issue asks about, one merge-order point with #1346, one test that fails on a KLU-off build, and lows. What we checked comes first, so it needn't be redone.

What we checked

check	result
Fail before: the new tests on main's core	3 failed, 1 passed, under cc and MIR. The compiled test fails at residual 2.15e-39; [fd] and [fd colored] fail at 1.56e-40, CVODE failing at t = 0. [auto] passes on main, as you say.
After	4 passed under cc, BNGSIM_CODEGEN_JIT=mir, BNGSIM_CODEGEN_THRESHOLD=1, and both. The whole file: 228 passed, 2 skipped (BIOMD0000000358 not present).
Full python/tests, BNGSIM_SKIP_AUDIT=strict	16,391 passed, 166 skipped, 0 failed (your 16,419 / 138; this environment skips 28 more, BNG2.pl-gated among them)
parity_checks/tests	926 passed, 33 skipped
Corpus: steady_state(method="integration") on every .net in tests/data, benchmarks and parity_checks, interpreted and compiled, auto and fd, main's core against this PR	708 solves each side. 706 identical, bit for bit in every concentration, the residual, the run measure and the Jacobian source. The 2 that differ are your two glioblastoma rows: compiled auto goes from refused to converged at the same residual, and interpreted fd keeps its state and residual with the run measure 0.022 → 0.568.
Mutations, rebuilt, the whole file run	Caught: no retry in ss_eval_rhs (3 fail); the run test never switches it on (3); the colored quotient takes f unclamped ([fd colored]); every f in the run test clamped, not only a NaN one (4 existing tests fail, so "only then" is pinned); the shared mask rule with a rate-rule target clamped (4 fail in test_nonneg_clamp_scope.py). Not caught: no retry in the dense closed-form Jacobian (your M3, which macOS doesn't reach); no retry in the sparse closed-form Jacobian (low 4); the run test's own clamp ignoring the mask (low 4); an infinity retried too, and the take-on run clamped (your M7 and M5, by design).
KLU off (the configuration mir.yml builds)	1 failed, 227 passed, 2 skipped: [fd colored] (§3).
With #1350	Under BNGSIM_CODEGEN_THRESHOLD=1, its strict xfail citing #1357 on test_a_march_whose_integrator_fails_returns_the_state_it_passed[integration] is XPASS(strict), so whichever of the two lands second drops _COMPILED_MARCH_REFUSES, as both bodies say. At the default threshold the file passes.
Read	ss_continued_run_once and its three call sites: the two of the run test, which clamp, and the take-on run with &end, which doesn't; every Jacobian callback the marcher installs, each reaching the retry or ss_eval_rhs; KINSOL's callbacks, which share the user data but never switch it on; nonneg_clamp_mask_of against main's nonneg_clamp_mask, line for line; the time course's retry (cvode_rhs_body, read_values, spare_at_pole, clamp_state_nonneg) against the run test's.
Gates	ruff check and format, clang-format 17.0.6 and ci/changelog.py check all clean
Overlap	Trial merges of the other 41 open PRs onto main e425546 + this PR: it adds no conflict, and none of the 41 conflicts with main. Nine share cvode_simulator.cpp, #1346 among them, and #1350 shares the test file. None touches steady_state.cpp or the new header.
CI	23 checks, all green. No leg runs the new tests under MIR, on Windows, or without KLU.

1. The issue's second ask still stands

#1357 also says the closed form's failure in the run test "does not retry on difference quotients … Whatever fails here is not recorded as integrator_failed." That is still so, on main and here. The solve's quotient retry (steady_state.cpp:6123) is decided before the run test (6295), and keys on integrator_failed, which the run test never sets. This PR removes the NaN-undershoot cause. A run test the closed form can't take for any other reason still refuses, under jacobian="auto", a state fd accepts. Two two-species models, each started at S = 0 exactly (residual 0), with a sink read through a law whose closed-form derivative has no value at 0:

model	auto	analytical	fd	codegen=True	run() from the same state
k*exp(-1/S)	refused, run measure inf, retried=False	refused	converged, run measure 0	refused	to 1e6, S = 0
kp*S^0.5	refused, run measure inf, retried=False	refused	converged, run measure 0	refused	to 1e6, S = 0

The clamp can't help here, because the state is already at 0. The time course survives by falling back to the quotient (GH #176). One more run test on fd when the closed form's returns inf under auto, or setting the flag the 6123 retry reads, would close it. Otherwise, it would be worth not closing #1357 or opening a new issue for it, so it stays tracked.

2. Land with #1346 in mind: the run test's clamp is a copy

The new header moves the mask, not the clamp. clamp_nonneg_into repeats clamp_state_nonneg's loop (minus spare), and ss_nan_and_no_inf repeats read_values(...) == ClampReading::nan_only. #1346 changes the time course's copy: on the non-finite retry, it reads a slot the mask leaves alone, negative by less than its atol, as 0 (MODEL1812040005's rate-rule target P under a*P^z). The run test's copy won't follow, so a run test from such a state would still meet the NaN that run() no longer does, which is #1357's asymmetry again.

With #1346 rebased this morning (d787eac), we built main e425546 + this PR + #1346. Both PRs' test files pass together: 234 passed, 3 skipped (MODEL1812040005 not materialized here), under cc and MIR. And the asymmetry shows up on #1346's own model shape: P' = -1e-30 from P = 0, with A decaying at k*A*P^0.5.

build	run() to 1e6	steady_state(method="integration"), auto and fd
main + this PR	fails at t = 5e5 (CV_REPTD_RHSFUNC_ERR)	refused: residual 5e-31, run measure inf
main + this PR + #1346	A = 1, P = -1e-24	still refused: residual 5e-31, run measure inf

It's the same in either species order. P moves 1e-24 over the run, far inside its tolerance.

One clamp and one reading in nonneg_clamp.hpp, with spare (and #1346's band) as inputs, would keep "one rule, two integrators" true by construction. Failing that, whichever of the two lands second should carry the change to the run test, or say why it differs.

3. [fd colored] fails on a KLU-off build

force_sparse_linear_solver is a no-op without KLU (_simulator.py:648), so assert (ss.linear_solver == "klu") == sparse is False == True. KLU-off is a supported configuration (conftest.py:288). We built this PR with KLU off and MIR on, as mir.yml does: the file gives 1 failed, 227 passed. The case as

pytest.param("fd", True, marks=pytest.mark.skipif(not bngsim.HAS_KLU, reason="KLU not compiled"), id="fd colored"),

skips cleanly, and the strict audit takes the reason.

Low

4. Two retries and the run test's mask are reached by no test here.
   • The sparse closed-form Jacobian retry (cvode_ss_sparse_jac_body) is reached by no test on any platform: [fd colored] runs the colored quotient, and the compiled test's march is routed dense. A KLU leg of the compiled test would reach it, at least where the dense one is reached.
   • The dense one is pinned on x86-64 Linux only, as you say.
   • The run test's own clamp ignoring the 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 mask passes the file. The rule is pinned through the time course, but not that the run test applies it.
5. The compiled test's skip isn't the compiler bngsim uses. _CC_ON_PATH looks for cc, gcc or clang on PATH. _find_c_compiler takes $CC first, and on Windows cl then gcc, and under MIR no compiler is needed. So the test skips on an MSVC-only host, with CC off PATH, or on a MIR-only machine, though codegen would work there, and "no C compiler…" is declared _LOCAL_ONLY. We read this; we didn't run it on Windows. shutil.which(cg._find_c_compiler()[0]) in a try asks bngsim itself. The _has_cc() pattern elsewhere is no better on Unix, where _find_c_compiler returns ["cc"] unchecked.
6. A still-NaN value goes to CVODE as a success (pre-existing). cvode_ss_rhs_body and the two closed-form Jacobian bodies return 0 when the retry leaves a NaN. The time course returns 1 there (recoverable; cvodes: an analytical Jacobian still non-finite after the GH #135 clamp is used anyway, and silently freezes that species' Newton correction #867, "not a silent success"), and CVODE cuts the step.
7. Names.
   • clamp_nonfinite and clamp_nonfinite_as_the_time_course retry a NaN only, where the time course also retries an infinity (sparing poles). The body says so, but the names and the title don't.
   • Whether the run clamps is decided by end == nullptr, whether the caller wants the end state. It's right for the three call sites today. An explicit argument would read as the comment does.

Not checked:

• x86-64 Linux beyond CI. The dense closed-form half is pinned only there, and CI's Linux legs passed.
• Windows (no leg runs the file).
• Your build-hook measurements.

[wshlavacek]

A note for whoever merges this: #1350 (#1217) lands first, and this PR has to undo one of its marks.

#1350 strict-xfails test_a_march_whose_integrator_fails_returns_the_state_it_passed[integration] in test_steady_state_scaled_criterion.py under a lowered BNGSIM_CODEGEN_THRESHOLD. It cites #1357 and lives in _COMPILED_MARCH_REFUSES. This PR fixes #1357, so once it has main merged in, that case passes under BNGSIM_CODEGEN_THRESHOLD=1 and the strict xfail fails the forced run. No workflow lowers the threshold, so this PR's CI won't show it. #1358's xfail went stale the same way.

Once #1350 is in main, this PR should:

• delete _COMPILED_MARCH_REFUSES and restore @pytest.mark.parametrize("method", ["integration", "newton"]);
• delete xfail_when_threshold_forced from python/tests/_codegen_backend.py and from that file's import, since nothing else uses it.

Check: BNGSIM_CODEGEN_THRESHOLD=1 python -m pytest python/tests/test_steady_state_scaled_criterion.py -k returns_the_state_it_passed passes with no xfail.