[Montana]

Closes #1357.

What was wrong

steady_state(method="integration") on glioblastoma_egfrviii_seed.net, on main 9a26f5b:

codegen	jacobian	converged	residual	Jacobian	continued_run_distance
False	auto	True	2.7e-12	analytical	0.091
True	auto	False	2.1e-39	codegen	inf
True	fd	True	6.7e-12	finite-difference

What it was not

It was not a wrong compiled Jacobian. I called bngsim_codegen_jac through ctypes and compared it with the interpreter's Model.jacobian, at the steady state and at three states perturbed by 10%. The two are identical to the bit (max difference 0), and both match a central difference to 1e-16 of the matrix's scale. A time course on the compiled path from the steady state also runs to t = 1e6, with the closed form and with differences.

What it was

The compiled right-hand side rounds a few ulps differently from the interpreter. That led the march and polish to a slightly different state, still at residual 2.1e-39. The scaled criterion's run test (ss_continued_run, #1001) then takes the state on for up to 1e6 at its tolerances. From that state, on the closed-form Jacobian, CVODE's corrector failed at t = 3.9e5, so ss_continued_run_once returned +inf, and the state was refused.

The solve's own retry on difference quotients (issue #127, !result.converged && integrator_failed && jacobian == "auto") runs before the run test, and the run test had none of its own. Started from that same state in a fresh solve, the compiled path converges. So the state is a steady state, and only this one run on the closed form fails.

The fix (src/steady_state.cpp)

The run test now gets the retry that the march (#127) and the sensitivity hold run already have. Under jacobian="auto" with a closed form installed, a run test that comes back inf is taken once more on difference quotients:

             result.continued_run_distance =
                 ss_continued_run(model, rhs, *effective, sub, result.concentrations);
+            if (!std::isfinite(result.continued_run_distance) && effective->jacobian == "auto" &&
+                ss_install_solver_jacobian(rhs, *effective)) {
+                SteadyStateOptions by_differences = *effective;
+                by_differences.jacobian = "fd";
+                result.continued_run_distance =
+                    ss_continued_run(model, rhs, by_differences, sub, result.concentrations);
+            }
             result.converged = result.continued_run_distance <= 1.0;

• solver_jacobian_retried is not set. Its documented meaning is that the answer is a retry's, and solver_jacobian_source then reads "finite-difference". Here the returned state is still the closed-form solve's; only the run that judges it is differenced. That matches the hold run, which retries the same way and does not set the field either.
• jacobian="analytical" is still not second-guessed, as steady_state: neither tier installs the Jacobian it already has (CVodeSetJacFn / KINSetJacFn) #127 leaves it.

Tests (test_steady_state_scaled_criterion.py)

• test_a_run_test_the_closed_form_cannot_take_is_taken_on_differences: glioblastoma with codegen=True converges, continued_run_distance <= 1, the source is codegen, solver_jacobian_retried is False, and the concentrations equal the interpreted solve's at the march's own tolerances (rtol and atol 1e-8). It fails on main.
• test_an_explicit_closed_form_is_not_called_off_in_the_run_test: control. The same model with jacobian="analytical" is still refused, with a run measure of inf.

The fix is the only difference between main's core and this one, so main's core is the mutation check.

With #1350

#1350 marks two tests with strict xfails citing #1357, under a forced threshold:

• test_a_march_whose_integrator_fails_returns_the_state_it_passed[integration] (glioblastoma) now passes under BNGSIM_CODEGEN_THRESHOLD=1. Whichever PR merges second drops that mark.
• test_a_march_whose_integrator_fails_is_still_retried_on_difference_quotients (MODEL1712050008) is not this bug. On every path, main and here give identical results. The compiled path converges, with residual 1.05e-15, continued_run_distance 7.3e-7, and a state within 1.6e-6 of the interpreter's, without a retry. The test asserts the arm64 interpreter's route (solver_jacobian_retried), so under a forced threshold it is a path assumption. I'll move it to tests: declare what a forced-JIT run cannot pass, so the run is a signal (#1217) #1350's threshold_forced skips.

Measured

• Every .net model in the repo (177: tests/data, benchmarks, parity_checks), steady_state() interpreted and compiled, on main's core and this one: 339 identical results and 14 identical errors. The one change is glioblastoma compiled, refused (crd=inf) → converged at the same residual 2.15e-39.
• Full suite on this core: 16,398 passed. The run loaded this file before the comparison above was widened to the march's own tolerances (rtol and atol 1e-8; one species reads 0 here and 9.1e-10 interpreted), so the new test failed in it. Rerun after that change, test_steady_state_scaled_criterion.py passes in full.