[Montana]

Fixes #1254.

Root cause

On the KLU route without an analytic Jacobian, CVODE has no difference quotient of its own for a sparse matrix, so bngsim fills the Jacobian by colored finite differences (colored_fd_jacobian, include/bngsim/sparse_jacobian.hpp). Each column's step was:

double h = sqrt_uround * std::max(std::abs(y[j]), 1.0);

That's a floor of 1 in absolute terms. A species at 1e-9 was stepped by 1.5e-8, fifteen times itself, so the matrix KLU factored was the Jacobian at a different state. Newton then converged slowly, failed to converge (CV_CONV_FAILURE), or was pushed to a state where the RHS was not finite (CV_REPTD_RHSFUNC_ERR).

The dense route doesn't have this problem, because CVODE's own difference quotient (cvLsDenseDQJac) takes each step on the error test's scale:

inc_j = max(√ε·|y_j|, minInc / ewt_j),   minInc = 1000·|h|·ε·N·‖f‖_wrms   (1 where ‖f‖ = 0)

Change

• sparse_jacobian.hpp: colored_fd_jacobian takes optional error weights and min_inc. When given them, it uses CVODE's rule above. Without them it keeps the old step, so the steady-state march's call is unchanged.
• cvode_simulator.cpp:
  • CvodeUserData carries the session's cvode_mem, set where it is registered.
  • The colored Jacobian callback reads the error weights into tmp1 (CVODE's scratch vector for this callback), along with the current step, and computes min_inc exactly as cvLsDenseDQJac does.

Why it is right

The colored method is exact column by column when the pattern is complete: no two columns in a color share a row. So with the same step rule, it gives the dense difference quotient's entries for the pattern. That is what this PR measures: on #1254's model and on the nanomolar chain below, KLU with the FD Jacobian now matches dense with the FD Jacobian.

Tests

python/tests/test_colored_fd_jacobian_steps.py (KLU-gated; added to the windows-klu.yml leg, which test_windows_klu_ci_coverage.py requires):

test	main	this PR
#1254's equality: KLU+FD's Y(10) is as close to 2 ln 5 as dense+FD's, at rtol 1e-8 and 1e-10, at 101 and 1001 samples (4 cases)	fails: 1.0e-3 to 1.2e-2 off, where dense is 1.4e-7 to 5.7e-6	passes: equal to dense, or closer (9.7e-10 at rtol 1e-10, 101 samples)
a dimerization chain at 1 nM: KLU+FD matches dense+FD and takes at most 1.2× its steps and Jacobians	fails: 781 steps and 100 Jacobians against 299 and 40	passes: 299 and 40
control: KLU with the analytic Jacobian is unchanged	passes	passes

Mutations (each rebuilt and run against the tests):

mutation	result
no min_inc term	5 failed
min_inc always 1, ignoring the step and norm	4 failed
weights never read, i.e. the old step	5 failed

What #1254 does and does not close

The KLU part is fixed: KLU with the FD Jacobian is now as right as dense with it.

Not fixed here: the equality itself remains fragile on every route. An equality on a state has no root to stop the solver at, so when its switch flips depends on where the steps fall. Even dense with the analytic Jacobian is 1.1e-3 off at rtol 1e-10 with 1001 samples:

rtol, samples	dense + analytic	dense + FD	KLU + FD, main	KLU + FD, this PR
1e-8, 101	1.1e-8	1.4e-7	1.0e-3	1.4e-7
1e-8, 1001	4.8e-8	2.6e-6	1.2e-2	2.6e-6
1e-10, 101	1.4e-9	5.7e-6	4.7e-3	9.7e-10
1e-10, 1001	1.1e-3	1.7e-7	3.7e-3	1.7e-7

That is a separate problem, an equality rate-law switch with no root registered, and I'll open it as its own issue if wanted.

The corpus: every model this can move

Models #1257's suite flagged (PR #1291), which run on KLU+FD only when forced:

model	main	this PR
BIOMD0000000827	CV_CONV_FAILURE	2.6e-6 of the peak from dense+FD
MODEL1006230034	CV_REPTD_RHSFUNC_ERR	3.2e-6
MODEL2002030001	CV_ERR_FAILURE	6.7e-6
MODEL1912090001	3.1e-2	1.6e-4
MODEL5952308332, MODEL2001200002	fail	still fail. Not this defect: both are fragile on every route. MODEL5952308332's dense analytic run goes from 4,104 to 598,459 steps at one tolerance tighter, and MODEL2001200002's dense FD run fails at rtol 1e-7.

Every model the automatic routing sends to KLU (163: 122 BioModels, 6 .net, 35 SBML suite and event models). Each was run on main and on this PR, against dense:

route	identical, or within 2×	closer to dense by >2×	farther by >2×	total steps
default (jacobian="auto")	157	2	0	69,253 → 69,282
jacobian="fd"	44	100	9	67,867 → 68,755 (+1.3%)

Some of the larger improvements, against dense+FD: BIOMD0000000152 9.6e-4 → 5.6e-11, metapop_sir_100 3.3e-5 → 1.9e-8, BIOMD0000000574 9.7e-5 → 1.7e-6, SHP2_base_model 3.5e-7 → 1.3e-12.

Of the 9 that drift farther from dense, 7 stay below their own rtol. The other two, measured against a reference at 10⁻⁴ × their tolerance:

• BIOMD0000000571: 6.1e-7 → 2.8e-6 at rtol 1e-6.
• MODEL1601140000: 8.2e-6 → 3.3e-5 at rtol 1e-6, where dense+FD is 6.4e-6. This is the one real, if modest, regression. My guess is that colored differences leave out the entries the pattern lacks, and dense differences don't.

BIOMD0000000153 moves from 3.1e-4 to 1.8e-1 from that reference. But every other route is already 1.8e-1 from it at rtol 1e-6: dense+FD, and KLU with the analytic Jacobian. The model isn't converged at its tolerance, and main's KLU+FD was nearer by chance. The fix makes it agree with every other route.

Suite

16,096 passed, 148 skipped, on macOS against main (023ea50), with the 3 failures below.

• test_a_relative_interpreter_stays_in_its_venv is my scratch venv's /tmp reading as /private/tmp; it fails the same way on main.
• test_skip_audit and test_windows_klu_ci_coverage flagged the first version of the new test file, which was right of them. Both are fixed: the test no longer needs Antimony, and the file is in windows-klu.yml's run list and paths filter. KLU_PASSED moves from 20 to 26, floor 24, slack 5 (below the smallest file, 6), per the floors file's own rule. All 115 tests of those meta-checks and the new file pass.

clang-format (pinned 17.0.6), ruff and ci/changelog.py check are clean.

Interaction with #1291

#1291's shape test marks equality-dwell / klu_vs_dense as a strict xfail on #1254. Whichever of the two merges second will see it unexpectedly pass, and must drop that xfail. I'll do it in whichever PR is left.

Remaining risks

• One model is measurably worse: MODEL1601140000, at 33× its rtol against a tight reference.
• The steady-state march is unchanged. Its colored Jacobian still takes the old step. The differenced state Jacobian steps every species by 1.5e-8 of the largest concentration: dY_ss/dp with jacobian="fd" is 7% off beside an unrelated species at 1e7 and 44% at 1e8, with no error #1002 is that path's own issue, and this PR doesn't touch it.
• Platform: checked on macOS arm64 only. windows-klu.yml will run the tests on Windows.

Maintainer review: the three strongest objections

1. "Why CVODE's rule rather than a smaller floor?" Any fixed floor is wrong for some model's scale. CVODE's rule is the one the dense route already uses, so the two routes now agree by construction rather than by tuning. It's also the rule the error test itself measures accuracy by.
2. "The step now depends on CVODE's internal state." It reads the public accessors CVodeGetErrWeights and CVodeGetCurrentStep, which give exactly what cvLsDenseDQJac reads internally. Without a session, the old step stays.
3. "One model got worse." Yes: MODEL1601140000, from 8e-6 to 3.3e-5 at rtol 1e-6, and the PR says so. Against that, 100 models come substantially closer to dense, three crashing models run, and An exact equality on a state, if(Bobs==0, kb, 0), is silently wrong on KLU with the FD Jacobian: Y(10) off by 1e-3 at rtol 1e-8 and 4.7e-3 at 1e-10, where dense and the analytic Jacobian are right #1254 is fixed. A follow-up could look at the pattern's completeness for that model.

Diff (source; the test, the workflow, the floors and the changelog fragment are not shown)

diff --git a/include/bngsim/sparse_jacobian.hpp b/include/bngsim/sparse_jacobian.hpp
index f1b9f635..943bdb9c 100644
--- a/include/bngsim/sparse_jacobian.hpp
+++ b/include/bngsim/sparse_jacobian.hpp
@@ -207,10 +207,19 @@ inline void install_csc_structure(IndexT *col_ptrs, IndexT *row_indices,
 // scratch buffers are ns long — the caller owns them so a per-step Jacobian
 // allocates nothing. `eval_rhs(t, y_pert, fy_pert)` evaluates the SAME right-
 // hand side the session is integrating.
+//
+// The step. With the session's error weights `ewt` (and `min_inc` from them), a
+// column's step is CVODE's own dense difference quotient's (cvLsDenseDQJac):
+// max(sqrt(eps)*|y_j|, min_inc/ewt_j), on the scale the error test judges y_j by.
+// Without them, sqrt(eps)*max(|y_j|, 1): a floor of 1 in absolute terms, which
+// moved a species of 1e-9 by fifteen times itself, so the matrix KLU factored
+// was the Jacobian somewhere else, and Newton failed to converge or reached a
+// state with no finite rate where dense differences, on the weights, ran (#1254).
 template <typename EvalRhs>
 inline void colored_fd_jacobian(const JacobianSparsity &sp, int ns, double t, const double *y,
                                 const double *fy, double *jac_data, double *y_pert, double *fy_pert,
-                                double *h_vals, EvalRhs &&eval_rhs) {
+                                double *h_vals, EvalRhs &&eval_rhs, const double *ewt = nullptr,
+                                double min_inc = 1.0) {
     // Finite difference perturbation scale
     const double sqrt_uround = 1.4901161193847656e-8; // sqrt(machine epsilon)
 
@@ -222,7 +231,8 @@ inline void colored_fd_jacobian(const JacobianSparsity &sp, int ns, double t, co
         std::memcpy(y_pert, y, static_cast<std::size_t>(ns) * sizeof(double));
 
         for (int j : group) {
-            double h = sqrt_uround * std::max(std::abs(y[j]), 1.0);
+            double h = ewt != nullptr ? std::max(sqrt_uround * std::abs(y[j]), min_inc / ewt[j])
+                                      : sqrt_uround * std::max(std::abs(y[j]), 1.0);
             h_vals[j] = h;
             y_pert[j] += h;
         }
diff --git a/src/cvode_simulator.cpp b/src/cvode_simulator.cpp
index 8840c0cc..ac8693c4 100644
--- a/src/cvode_simulator.cpp
+++ b/src/cvode_simulator.cpp
@@ -394,6 +394,10 @@ struct StateSwitchHold {
 
 struct CvodeUserData {
     NetworkModel *model;
+    // The session these callbacks belong to, set where this is registered with
+    // CVodeSetUserData: the colored FD Jacobian reads its error weights and
+    // step (#1254). Null where nothing has registered it.
+    void *cvode_mem = nullptr;
     // The largest |x| of each species at any state the right-hand side is
     // asked at, where a run keeps it (issue #1025): an error made where a
     // species was large is carried to where it is small, between outputs as
@@ -1942,7 +1946,7 @@ static int cvode_codegen_sparse_jac(sunrealtype t, N_Vector y, N_Vector fy, SUNM
 // (issue #128); what stays here is the CVODE plumbing around it.
 
 static int cvode_colored_jac_body(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
-                                  void *user_data, N_Vector /*tmp1*/, N_Vector /*tmp2*/,
+                                  void *user_data, N_Vector tmp1, N_Vector /*tmp2*/,
                                   N_Vector /*tmp3*/) {
 
     auto *data = static_cast<CvodeUserData *>(user_data);
@@ -1977,10 +1981,24 @@ static int cvode_colored_jac_body(sunrealtype t, N_Vector y, N_Vector fy, SUNMat
     double *fy_pert = data->colored_jac_fy_pert.data();
     double *h_vals = data->colored_jac_h_vals.data();
 
-    colored_fd_jacobian(sp, ns, static_cast<double>(t), y_data, fy_data, jac_data, y_pert, fy_pert,
-                        h_vals, [model](double tt, const double *yy, double *ydot) {
-                            model->compute_derivs(tt, yy, ydot);
-                        });
+    // Each column's step on the scale the error test judges it by, as CVODE's
+    // dense difference quotient (cvLsDenseDQJac) takes it: the session's error
+    // weights, in tmp1 (CVODE's scratch for this callback), and its step size.
+    const double *ewt = nullptr;
+    double min_inc = 1.0;
+    if (data->cvode_mem != nullptr && CVodeGetErrWeights(data->cvode_mem, tmp1) == CV_SUCCESS) {
+        ewt = N_VGetArrayPointer(tmp1);
+        sunrealtype h = 0.0;
+        CVodeGetCurrentStep(data->cvode_mem, &h);
+        const double fnorm = N_VWrmsNorm(fy, tmp1);
+        constexpr double uround = 2.220446049250313e-16;
+        min_inc = fnorm != 0.0 ? 1000.0 * std::abs(h) * uround * ns * fnorm : 1.0;
+    }
+
+    colored_fd_jacobian(
+        sp, ns, static_cast<double>(t), y_data, fy_data, jac_data, y_pert, fy_pert, h_vals,
+        [model](double tt, const double *yy, double *ydot) { model->compute_derivs(tt, yy, ydot); },
+        ewt, min_inc);
 
     return 0; // success
 }
@@ -3953,6 +3971,7 @@ Result CvodeSimulator::Impl::run_warm(const TimeSpec &times, const SolverOptions
         w.user_data.atol_tracking = make_atol_tracking(rtol, atol_v, opts.atol_track_decades);
         apply_cvode_tolerances(w.cvode_mem, w.ctx, rtol, atol, atol_v, ns,
                                w.user_data.atol_tracking.active() ? cvode_tracking_ewt : nullptr);
+        w.user_data.cvode_mem = w.cvode_mem;
         CVodeSetUserData(w.cvode_mem, &w.user_data);
         CVodeSetMaxNumSteps(w.cvode_mem, max_steps);
         if (opts.max_step_size > 0) {
@@ -4343,6 +4362,7 @@ void CvodeSimulator::Impl::create_cvode_core(const TimeSpec &times, const Solver
     apply_cvode_tolerances(cvode_mem, ctx, rtol, atol, atol_v, ns,
                            user_data.atol_tracking.active() ? cvode_tracking_ewt : nullptr);
 
+    user_data.cvode_mem = cvode_mem;
     flag = CVodeSetUserData(cvode_mem, &user_data);
     flag = CVodeSetMaxNumSteps(cvode_mem, max_steps);