[Montana]

Fixes #1298.

The bug

A sensitivity run on MODEL2001200002 returned a different trajectory from its plain run, and columns that were wrong. The only sign was a warning.

	plain run	sensitivity run, main	this PR
state, vs a plain run at 1e-3× the tolerance	1.3e-5	1.2e-2	1.3e-5 (identical to the plain run)
∂/∂sigma_nk, vs a central difference of plain runs		0.28 of its peak	2.0e-7
∂/∂g_reg		0.96	8.9e-7
∂/∂p_IC		1.0	4.4e-5

The state row is the corpus job (rtol 1e-6, atol 1e-9). The column rows are the test's fixture at rtol 1e-9, atol 1e-12.

Why

L decays to about 1e-16 by t = 25 and sits a hair below zero from there. There D = d_im*(L/T)^1.36/(s + (L/T)^1.36)*T is NaN.

1. The value RHS retries at the clamped state (Extract the two forward-sensitivity jump handlers out of CvodeSimulator::run (follow-up to #109) #135), so the run integrates f(clamp(y)). The analytical Jacobian does the same, so the Newton matrix is J(clamp(y)).
2. The sensitivity RHS had only 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 narrower retry, which snaps a component inside (-atol_i, 0). That retry is what lets the run step a hair past zero: a state with L in the band is accepted.
3. Once L is past zero, the value RHS's clamp holds it at the value it reached rather than restoring it (the clamped RHS has no -λ·Treg·L term). It drifts past -atol_i (L = -1.39e-12 against atol 1e-12).
4. From then on the sensitivity RHS refused every state. CVODES cut h until it gave up, and the Four Linux-only test failures, found by #169's first whole-suite run: the GH #176 FD fallback dies at a second crossing, and a steady-state fixture is exactly singular under reference LAPACK #176 retry returned a run with the finite-difference Jacobian instead.

BNGSIM_SENS_CLAMP_NUMERIC_ZERO=0, which turns #392's retry off, gives the right answer on main too, because every NaN then rejects the step and keeps L ≥ 0. So this PR doesn't revert #392, which fixed BIOMD0000000833. It gives the run that #392 lets past zero a consistent way to continue.

The fix

After #392's retry declines, the sensitivity RHS retries at clamp_state_nonneg's state, the one the value RHS integrated. It does this only where f(y) itself is non-finite.

That condition is the point. Where f(y) is finite at a negative base, the run integrates f(y) there, and a column read at clamp(y) would belong to another state. That is #392's pow(x, 4.0) beside log(x). Such a NaN still goes to the #384/#386 refusal. The new control test covers this case, and a mutation that drops the condition fails it.

On #392's reasoning

#392 kept its band narrow for a stated reason: a base more than atol negative, such as V_membrane = -61 in BIOMD0000000374, has no ∂/∂n base^n, so its NaN should be refused rather than answered with a confident 0.

That reason still holds, and this PR keeps it, in two ways:

• The clamp mask. The zero-base exponent guard tests base == 0, and the solver hands it a base 24 orders below atol on the wrong side #392 landed on 2026-08-17. 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 (2026-09-25) later restricted clamp_state_nonneg to slots that are concentrations by construction. A rate-rule target, a promoted parameter, a fixed species and a species declared negative are never touched. So V_membrane keeps its NaN.
• The f(y) condition. The retry fires only where the run is integrating the clamped flow anyway.

I've added a note to #392's comment saying so.

Changes

diff --git a/src/cvode_simulator.cpp b/src/cvode_simulator.cpp
index 8840c0cc..9f5a2dbb 100644
--- a/src/cvode_simulator.cpp
+++ b/src/cvode_simulator.cpp
@@ -772,7 +772,10 @@ static inline double *clamp_state_nonneg(CvodeUserData *data, const double *y_pt
 // claim entirely — `BIOMD0000000374` carries a species `V_membrane = -61` — and
 // `∂/∂n base^n` does not exist there at all, so that component is left alone and
 // its NaN survives to the GH #384/#386 refusal rather than being answered with a
-// confident `0`.
+// confident `0`. (That case is now also the clamp mask's to keep (#706): a slot
+// that is negative by design is one clamp_state_nonneg never touches, so the
+// sensitivity RHS's second retry, at the value RHS's own clamped state (issue
+// #1298), leaves it to the same refusal.)
 //
 // This is not the clamp GH #93 removed. That one was unconditional and lived on
 // one side of a pair: the interpreted Michaelis-Menten RHS clamped `sFree` where
@@ -1563,6 +1566,22 @@ static void capture_sens_witness(CvodeUserData *data, sunrealtype t, double *y_p
     }
 }
 
+// Whether the value RHS is non-finite at y itself, before any clamp: the test
+// that decides which state its clamp-retry integrated (issue #1298). Written
+// into `scratch`, through the same function the value RHS callback calls.
+static bool value_rhs_nonfinite_at(CvodeUserData *data, sunrealtype t, double *y_ptr,
+                                   N_Vector scratch) {
+    double *f = N_VGetArrayPointer(scratch);
+    if (data->codegen_fn != nullptr) {
+        if (data->codegen_fn(static_cast<double>(t), y_ptr, f, &data->codegen_so_data) != 0) {
+            return true;
+        }
+    } else {
+        data->model->compute_derivs(static_cast<double>(t), y_ptr, f);
+    }
+    return rhs_has_nonfinite(data, f);
+}
+
 static int cvode_codegen_sens_rhs_body(int Ns, sunrealtype t, N_Vector y, N_Vector ydot, int iS,
                                        N_Vector yS, N_Vector ySdot, void *user_data, N_Vector tmp1,
                                        N_Vector tmp2) {
@@ -1661,6 +1680,37 @@ static int cvode_codegen_sens_rhs_body(int Ns, sunrealtype t, N_Vector y, N_Vect
                 return 0;
             }
         }
+        // Issue #1298: then at the state the value RHS and the Jacobian were
+        // evaluated at, where that is a different state. Where f(y) is
+        // non-finite, cvode_rhs integrates f(clamp(y)) and the Newton matrix
+        // is J(clamp(y)), so the columns of that same flow are
+        // J(clamp(y))·s + f_p(clamp(y)). The retry above lets a run step a
+        // hair past zero; once there, the value RHS's clamp holds the species
+        // at the value it reached rather than restoring it, it drifts past
+        // -atol_i, and this RHS refused every state from then on: CVODES cut
+        // h until it gave up, and the GH #176 retry returned another run's
+        // trajectory (MODEL2001200002: 2.6e-3 off in the state, every column
+        // wrong).
+        //
+        // Only where f(y) itself is non-finite: where it is finite the run
+        // integrates f(y), at the negative base, and a column read at
+        // clamp(y) would be another state's (#392's `pow(x, 4.0)` beside
+        // `log(x)`), so that NaN still goes to the refusal. So does a slot
+        // nonneg_clamp_mask leaves alone, negative by design (#706). f(y) is
+        // written to tmp1, the scratch CVODES gives this call, before the
+        // sensitivity function reuses it. Off with
+        // BNGSIM_SENS_CLAMP_NUMERIC_ZERO=0, like the retry above.
+        if (!data->sens_zero_atol.empty() && value_rhs_nonfinite_at(data, t, y_ptr, tmp1)) {
+            rc = sens_fn(Ns, static_cast<double>(t), clamp_state_nonneg(data, y_ptr),
+                         N_VGetArrayPointer(ydot), iS, N_VGetArrayPointer(yS), ySdot_ptr, &so_data,
+                         N_VGetArrayPointer(tmp1), N_VGetArrayPointer(tmp2));
+            if (rc != 0) {
+                return rc;
+            }
+            if (!rhs_has_nonfinite(data, ySdot_ptr)) {
+                return 0;
+            }
+        }
         capture_sens_witness(data, t, y_ptr, Ns, iS, ySdot_ptr, so_data);
         return 1; // recoverable -> CVODES cuts h and retries
     }

Plus python/tests/test_sens_rhs_at_the_clamped_state.py and changelog.d/1298.fixed.md.

Tests (test_sens_rhs_at_the_clamped_state.py)

The fixture is MODEL2001200002 (Mpekris 2017, BioModels, CC0) as a 7-species Antimony string with its unused parameters dropped, so the test needs no corpus file.

test	main	this PR
sensitivity run's state equals the plain run's (staggered)	fails: 3.4e-3 of the peak	passes
same, simultaneous	passes (already right here)	passes
columns vs _fd_sens.central_difference (sigma_nk)	fails: 122 cells	passes
columns vs _fd_sens.central_difference (g_reg)	fails: 183 cells	passes
control: an exponent column at a negative base where f(y) is finite is still refused	passes	passes
control: a column finite there still runs (dA/dv = -4 at t = 4, exact)	passes	passes

The columns are judged by the repo's oracle (python/tests/_fd_sens.py). Its steps are chosen from the measured noise floor, and it refuses kinks rather than comparing them. It judges the species that stay above atol (C, N, Treg). T, Csc and I fall to 1e-34 and L sits below zero, where a plain run's own error exceeds the value and the oracle (correctly) refuses those cells as kinks.

Mutation: removing the f(y) condition makes only the control fail.

Corpus A/B

Every BioModel under 1 MB was run as a sensitivity run on three nonzero primary parameters (compartment sizes excluded), at its job's tolerance, with main's behaviour and with this PR's. That's about 1,190 runs; the rest are refused identically on both sides for unrelated reasons.

	models
bit-identical (state and every column)	1,187
MODEL2001200002	state 1.2e-2 → 1.3e-5 off the reference, now identical to the plain run
BIOMD0000001065	columns move by ≤ 1e-8 of their peak; state unchanged against the reference (2.65e-6)
MODEL5952308332	columns move by ≤ 2.1e-6 of their peak; state 1.0e-5 → 3.4e-5 off the reference, both better than its own plain run's 3.9e-4

The final code gives exactly the prototype's results on all three changed models.

Full suite

16,131 passed, 148 skipped. One failure, test_perf_ab.py::test_a_relative_interpreter_stays_in_its_venv, fails identically on main locally: macOS spells /tmp as /private/tmp. clang-format 17.0.6 and ruff are clean.

CI: Python suite · ubuntu-latest · 1/2 failed one timing check. test_codegen_gen_scaling.py::test_large_model_generation_is_not_quadratic took 22.7 s against its 20 s wall-clock budget. This PR doesn't touch code generation, the test it times. Locally the test takes 4.2 s on this branch and 3.9 s on main, and the same shard passed on ubuntu-latest-lapack and macos-14. It reads as a slow runner. A rerun of that job should clear it; I don't have the rights to trigger one.

Compatibility and performance

• Only runs whose sensitivity RHS was already non-finite after 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 retry are affected. Every other run takes exactly the path it took before.
• The retry costs one value-RHS evaluation, plus one sensitivity-RHS evaluation, on that path only.
• BNGSIM_SENS_CLAMP_NUMERIC_ZERO=0 still turns both retries off.

Not fixed here: out-of-step reads after the refusal

On main, after the sensitivity RHS's NaN, this model's run reads its state at times outside the current step. Examples:

• CVodeGetDky at t = 49 when the step spans [49.0087, 49.0261];
• CVodeGetSensDky1 at t = 2.6e-312, a garbage time.

There are 56 such errors in one run, before it dies with CV_ILL_INPUT. This reproduces on an untouched main build. With this PR the path is no longer reached on this model, but the defect is still there for any run that reaches it. Filed as #1301.

Risks and likely objections

1. "A NaN should be refused, not retried around." The run is already integrating f(clamp(y)) at these states, so the retry doesn't invent a state. Refusing there doesn't stop the run either: it ends in the Four Linux-only test failures, found by #169's first whole-suite run: the GH #176 FD fallback dies at a second crossing, and a steady-state fixture is exactly singular under reference LAPACK #176 retry returning another trajectory, with a warning.
2. "This undoes 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 narrowness." It doesn't widen 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 band. Where f(y) is finite, nothing changes. A slot negative by design is the mask's to keep (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).
3. "The test model is one corpus model." The corpus A/B is the broader evidence: about 1,190 BioModels, three changed, none for the worse.

Remaining risk: a model whose value RHS is non-finite at a negative concentration and whose sensitivity was relying on the refusal to stay non-negative. With the fix, such a model follows its plain run instead. None appears in the corpus. The trajectory it would follow is the plain run's, so the sensitivity run would now match it, and #1297 is where that plain-run behaviour belongs.

[elyfmiller]

A heads-up before we review this one, which is ours this week: since #1354 merged (main c0877fa), it no longer builds against main.

• The type error. clamp_state_nonneg now takes (CvodeUserData *, N_Vector y, const char *spare = nullptr), and nothing else accepts a state array. The new retry in the sensitivity RHS calls clamp_state_nonneg(data, y_ptr) with a double *.
• Why GitHub shows it clean: the two changes don't overlap in text, and this PR's CI ran on 10-27, before Refuse a run carried to a pole instead of clamping it there (BIOMD0000000338) #1354.

Could you rebase it on main? #1354 also changed what the clamp does:

• it leaves a slot alone at a pole the run landed on (spare);
• it judges its band where the step started;
• it refuses a run carried far below 0.

So the state this retry evaluates the sensitivity RHS at may want the same spare as the value RHS.

We'll review it with #1346 and #1370, which change the same clamp, once all three are on current main.