[Montana]

Fixes #1296. Found by the differential suite's ajac_vs_fdjac check (#1257, PR #1291).

Root cause

route_to_sparse_linear_solver (include/bngsim/sparse_jacobian.hpp) sends a model to KLU when it has 50 or more species and a pattern density below 10%. A structurally empty Jacobian (nnz == 0: no reaction moves any species, and an event assigns them) has a density of 0, so it passes, and the automatic routing picks KLU.

KLU then has nothing to factor. Under jacobian="fd" the colored differences have nothing to color, and the run reaches the backstop refusal in cvode_simulator.cpp. That refusal was written for force_sparse_linear_solver and says so:

SimulationError: force_sparse_linear_solver: this model cannot supply a sparse Jacobian — its analytical
Jacobian is unavailable (suppressed by jacobian="fd") and its Jacobian sparsity pattern has no structural
nonzeros to color for finite differences, … Drop the flag to use the dense solver, …

The caller never set the flag. Codegen's own reading of the same rule already required a nonzero: its fallback for a core that doesn't report routes_sparse is ns >= 50 and nnz > 0 and density < 0.10. The C++ rule had drifted from it.

Change

The size-and-density branch of the rule also requires sp.nnz > 0. That's one condition, in the one header the time-course solver, the steady-state march and the codegen emitter all read (routes_sparse), so the three stay in agreement. force_sparse_linear_solver is untouched: it still asks for KLU, and is still refused by name where nothing can fill the matrix, which is right for a caller who forced it.

Why it is right

A structurally empty Jacobian is J = 0, so the Newton matrix is I − γJ = I. The dense solver differences the full matrix with CVODE's difference quotient. That needs no pattern and is exact for J = 0. KLU had no advantage to lose: it factors the identity at best, and here it couldn't even be filled.

Tests

python/tests/test_empty_jacobian_routes_dense.py, on a 60-species model that no reaction moves, with an event at t = 5 assigning S_i = S_i + 1:

test	main	this PR
the pattern is empty (nnz == 0)	pass	pass
routes_sparse is False for it	fail	pass
a 60-species diagonal model still routes to KLU (control)	pass	pass
jacobian="auto" gives the closed form, S_i = i then i + 1	pass	pass
jacobian="fd" gives the closed form	fail (refused)	pass
force_sparse_linear_solver=True is still refused by name	pass	pass

Blast radius: every model whose route changes

The engine's own routes_sparse was compared on main and on this branch for every model on disk:

corpus	models	routed to KLU on main	route changed
BioModels	1,242	122	0
SBML Test Suite (semantic)	1,614	—	6: 01466, 01530–01533, 01564
SBML event models	212	—	0
.net corpus	68	—	0

The six are exactly the models with an empty pattern and 50 or more species. Against the SBML Test Suite's expected results:

case	main, auto	main, fd	this PR, auto	this PR, fd
01466	0	refused	0	0
01530	1.3e-8	refused	1.3e-8	1.3e-8
01531	5.5e-16	refused	5.5e-16	5.5e-16
01532	3.2e-16	refused	3.2e-16	3.2e-16
01533	0	refused	0	0
01564	1.8e-8	refused	1.8e-8	1.8e-8

Suite

16,071 passed, 148 skipped, on macOS against main (0062aeb). The 1 failure, test_a_relative_interpreter_stays_in_its_venv, is my scratch venv's /tmp reading as /private/tmp, and fails the same way on main. All 28 test files that use KLU or the routing pass (454 tests). clang-format (pinned 17.0.6), ruff and ci/changelog.py check are clean.

Remaining risks

• No C++ unit test. The native C++ suite builds with KLU off, where this rule always returns false, so a C++ test of the nnz condition would never run in CI. The Python tests exercise it where KLU exists: through the engine's own routes_sparse and end to end.
• Platform: checked on macOS arm64 only.

Maintainer review: the three strongest objections

1. "Better to fix the refusal's message than the routing." The message is right for the caller it was written for, someone who forced KLU, and stays as it is. The bug is that the automatic routing put a model there at all. Fixing only the message would turn a misleading refusal into an accurate one, for a run that has no reason to be refused.
2. "nnz == 0 is a corner case." It's 6 SBML Test Suite cases, and every model whose species only events move. More to the point, codegen already assumed the rule excluded it. This makes the engine agree with the code that reads it.
3. "Could this move a real model off KLU?" Only a model with no structural nonzero, which by construction has nothing for KLU to do. The table above shows the route changes for no BioModel, no .net model and no event model.

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

diff --git a/include/bngsim/sparse_jacobian.hpp b/include/bngsim/sparse_jacobian.hpp
index f1b9f635..9655df03 100644
--- a/include/bngsim/sparse_jacobian.hpp
+++ b/include/bngsim/sparse_jacobian.hpp
@@ -142,6 +142,15 @@ inline bool fill_routes_to_dense(const JacobianSparsity &sp, int ns, double lu_f
 //
 // `lu_fill` is routing_lu_fill(model, ns, jacobian_strategy): a model the size
 // and density test sends to KLU goes dense when fill_routes_to_dense says so.
+//
+// A structurally empty Jacobian (nnz == 0: species no reaction moves, an event's
+// targets) has a density of 0, which the density test reads as sparse, but it
+// gives KLU nothing to factor and, under jacobian="fd", nothing to color, and the
+// run was refused for force_sparse_linear_solver, a flag the caller never set
+// (SBML Test Suite 01530). So the size and density test sends it to the dense
+// solver, as codegen's own reading of this rule already did (_codegen.py's
+// fallback for a core that does not report routes_sparse). force_sparse still
+// asks for KLU, and is refused by name where nothing can fill the matrix.
 inline bool route_to_sparse_linear_solver(const JacobianSparsity &sp, int ns,
                                           const std::string &jacobian_strategy, bool force_dense,
                                           bool force_sparse, double lu_fill) {
@@ -149,7 +158,7 @@ inline bool route_to_sparse_linear_solver(const JacobianSparsity &sp, int ns,
     const bool sparse_ok = (jacobian_strategy != "jax") && !sp.empty();
     return sparse_ok && !force_dense &&
            (force_sparse ||
-            ((ns >= SPARSE_THRESHOLD) && (sp.density < SPARSE_DENSITY_MAX) &&
+            ((ns >= SPARSE_THRESHOLD) && (sp.nnz > 0) && (sp.density < SPARSE_DENSITY_MAX) &&
              !fill_routes_to_dense(sp, ns, lu_fill, jacobian_strategy, force_dense, force_sparse)));
 #else
     (void)sp;

[elyfmiller]

Thanks. We reviewed this independently, Opus 5.5 and Fable 5.1, each on our own build of main c0877fa with f66d1e4 merged (clean), on macOS arm64 with KLU on. Neither of us found a correctness defect.

The rule change is right where it sits. The one conjunct sp.nnz > 0 is inside the size-and-density branch that the two force flags straddle, so it changes only the automatic route. All three readers go through route_to_sparse_linear_solver:

• the time course (choose_use_sparse);
• the steady-state march (ss_use_sparse_linsol);
• codegen (codegen_jacobian_plan().routes_sparse).

_codegen.py:7831's fallback already had nnz > 0. What we checked is recorded below so it needn't be redone, followed by five low items, all docs, comments and tests.

What we checked

check	result
The new test file on main, which is this PR with the conjunct removed	2 failed, 4 passed: test_an_empty_jacobian_routes_dense and test_the_run_is_the_closed_form[fd], the rows the body names. So removing the conjunct is caught. With this PR: 6 passed.
Mutation the other way: nnz > 0 also gating force_sparse_linear_solver	test_forcing_klu_is_still_refused_by_name fails, so over-tightening onto the forced path is caught too
Full python/tests on the merged build	16,393 passed, 166 skipped, 0 failed, with BNGSIM_SKIP_AUDIT=strict (Fable). 15,944 passed on the second build; its one failure is test_run_tests_sh's walk-up probe, which fails the same way on main there because that venv lacks jax, diffrax and roadrunner.
The nine test files that read the routing	135 passed, 2 skipped (KLU-off only)
SBML Test Suite 01466, 01530–01533 and 01564 (pinned 473e119d), auto and fd, graded by the suite's own |x − e| ≤ atol + rtol·|e|	all 12 runs: 0 cells out of tolerance. Worst err/tol is 2.6e-4 (01564). Each case has n ≥ 51, nnz = 0, routes_sparse False, linear_solver 0. Forcing sparse on 01530 is still refused by name.
Route changes across the models on our machines: 1,641 BioModels SBML and 177 .net	No model has nnz = 0 with n ≥ 50, so none changes route. routes_sparse for the 177 .net models is identical on main and this PR.
The codegen path on the test's 60-species model	codegen=True with jacobian="fd" is refused on main. Here it runs on the dense solver, with the exact answer.
The steady-state march: 60 species, no reactions, no events	auto, fd and forced sparse all converge on "dense", unchanged from main
routing_lu_fill on an empty pattern	safe: the estimator pads the diagonal, and lu_fill_matters needs n ≥ 256
Gates	changelog fragment 884 characters and well-formed; ruff and clang-format 17.0.6 clean; "KLU not compiled" is a declared skip

Low items

1. docs/user-guide/solvers.md:341. The KLU row of the linear-solver table doesn't have the new condition. steady-state.md:266 already says "and a structural nonzero to factor"; the same phrase there would keep the two pages and the code in step.
2. python/bngsim/_simulator.py:645–672. The force_sparse_linear_solver docstring now misstates the rule.
   • It says "Only those two gates are bypassed"; with nnz > 0 there are three.
   • It justifies steady state's dense exception with "the steady-state auto rule can route such a model to KLU on its own (density 0 < 10%)", which this PR makes untrue.
3. src/steady_state.cpp:756–762. The comment above ss_use_sparse_linsol makes the same claim. The march's own sp.nnz > 0 check now matters only under force_sparse_linear_solver, where the march goes dense silently and run() refuses by name. That asymmetry predates this PR and is fine to leave, but the comment should describe the case that remains. Likewise, "the three stay in agreement" in the body holds for the automatic rule, not the forced one.
4. test_empty_jacobian_routes_dense.py:70–77, against the project's test-writing rules:
   • The route. The run tests pin the closed form but not the route the run took, which is what the fix changes. assert result.solver_stats["linear_solver"] != 1 on the fd row would pin it end to end; the plan's routes_sparse reads only the default options.
   • Codegen. A codegen=True case would cover the codegen path, which this PR also fixes.
   • Tolerances. pytest.approx(i) has no stated reason, while the answer is exact (f ≡ 0, and the event adds exactly 1). For S0 it rests on approx's default abs=1e-12. A one-line reason, or exact equality with the exactness stated, would do.
   • Optional. The module builds its SBML, and calls antimony.clearPreviousLoads(), at import time. A module-scoped fixture would keep that out of collection.
5. A note on the blast-radius section, not a change request: the dense route stores the Newton matrix as n × n. On a reaction-free model with an event at t = 5, the same SBML file on both builds, jacobian="auto", every answer exact:

species	main (KLU)	this PR (dense)
3,000	2.97 s, 147 MB peak	0.52 s, 309 MB
6,000	5.75 s, 185 MB	1.45 s, 720 MB

So it's faster at these sizes and grows as n² in memory. Nothing on disk is that shape, so it's worth a line in the body rather than a change.

Not checked: a KLU-off build (the #else branch is unchanged), the full 1,242-model and 1,614-case sweeps end to end, and Linux or Windows on today's main.

[elyfmiller]

Following up on our review above: we've pushed its low items 1–4 as one commit on top, 2c5f165, and are merging now that CI is green. Your change itself is untouched: sparse_jacobian.hpp and the changelog fragment are as you wrote them.

What the commit changes

• docs/user-guide/solvers.md. The KLU row of the linear-solver table now names the structural nonzero, as steady-state.md already did.
• The force_sparse_linear_solver docstring.
  • It now lists three gates; the third is an empty pattern. It also says the LU-fill test is bypassed, as the next paragraph already did.
  • The steady-state note no longer says the auto rule can send an empty pattern to KLU. Instead it says what remains under the flag: run() takes KLU with the analytical Jacobian and refuses with jacobian="fd", while the march stays dense.
• The comment above ss_use_sparse_linsol. The same correction, comment only.
• test_empty_jacobian_routes_dense.py.
  • The run rows pin the route (solver_stats["linear_solver"] != 1) as well as the answer.
  • The answer is compared with ==, since it's exact: f ≡ 0, so the solver never moves the state, and the event adds 1 to an integer.
  • A codegen=True case runs under "auto" and "fd".
  • The SBML is built in module-scoped fixtures instead of at import.

What we checked

check	result
The new test file on main 05ee803's core	5 failed, 3 passed, under cc and MIR. The 5: the route test, both run rows ([auto] on the route pin, linear_solver 1; [fd] refused) and both codegen rows ([auto] on the pin; [fd] refused). Your version failed 2 there. The 3 that pass are the controls.
The same file on main + this PR + the commit	8 passed, under cc and MIR
The docstring's claims, probed on both builds	forced, analytical Jacobian: KLU, exact; forced, fd: refused by name; the march on 60 species with no reactions: "dense" and converged, forced or not.
Full python/tests, BNGSIM_SKIP_AUDIT=strict, on main + this PR + the commit, before its last amend (one assertion and one docstring sentence)	16,788 passed, 140 skipped, 0 failed on one build. 16,351 passed and 430 skipped on the other; its one failure is test_run_tests_sh's walk-up probe, which needs jax and diffrax, and that venv has neither.
Gates	ruff 0.14.0 check and format, clang-format 17.0.6, ci/changelog.py check: clean
Overlap	Trial merges of the other 39 open PRs onto main + this PR, with and without the commit: no conflict either way. Eleven share a file with it; the nearest hunk is about 50 lines from the ones it edits (#1372 and #1385 in _simulator.py; #1370's in steady_state.cpp is about 60).
CI	22 checks, all green, on 2c5f165

Item 5 of the review, the dense route's memory, stays a note rather than a change. The dense route stores the Newton matrix as n²: at 6,000 species with no reactions, 720 MB against KLU's 185 MB, though faster at that size. Nothing on disk is that shape.

This lands ahead of #1336, as our review there suggested. Without it, an unflagged empty-pattern model run with fd raises #1336's new type.