[Montana]

Closes #1335.

#1335 lists three paths that #1302 (#1002) left without a test or a decision. No wrong answer was known on any of them, and this PR finds none. It pins each one and documents the third.

1 and 2: the two steady-state quotients with no jump mend

#1302 mends a difference-quotient row whose step crosses a jump for runs (dense and colored) and for the dense steady-state march. Two quotients were left as they were:

• KINSOL's own quotient: method="newton" with jacobian="fd". A declined law sends newton to the march instead, so that path is reached only with "fd".
• The steady-state march's colored sparse quotient: KLU.

New tests in test_difference_quotient_jump.py run each path, plus KINSOL on KLU, from:

• a start on a threshold, and 1e-9 above and below it, on X' = s − kx·k(X)·X with k = 0.04 above 5 and 0.02 at or below. With s = 0.08 the one steady state is X = 4, below the threshold; with s = 0.4 it is X = 10, above;
• beside a species held at 0, using the two 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 models that broke runs: the gate (Y = 0.5) and the inhibitor (Y = 20).

Each is checked against its closed form (< 1e-4 relative, the file's OK), and the held species against exactly 0.

Each test also asserts the path it took: method_used, linear_solver (dense/klu) and solver_jacobian_source == "finite-difference". Prototyping showed why that matters: with a declined law, method="newton" silently converges by the march instead, and would have "tested" KINSOL without touching it.

	cases	on main
threshold starts × 2 sources × 3 paths	18	all pass
held species × 2 models × 3 paths	6	all pass

The KLU paths skip on a KLU-off build, with the reason conftest.py already declares.

These pin the current behaviour rather than fix a defect. They are the "nothing pins that" #1335 asks for, and they agree with #1302's 72 and 24 variants that found no wrong answer.

3: SUNDIALS other than the vendored 7.2

include/bngsim/dq_jacobian.hpp reproduces cvLsDenseDQJac of the vendored SUNDIALS (v7.2.1) bit for bit. With -DBNGSIM_USE_SYSTEM_SUNDIALS=ON and another release, it is still a valid quotient, but not necessarily that release's own. #1335 asks for a CI leg or a documented refusal. This does the documented half:

• CMakeLists.txt: in the system-SUNDIALS branch, compare SUNDIALS_VERSION's major.minor with VENDOR.json's authoritative_release_tag, and warn on a mismatch:

  CMake Warning: System SUNDIALS 7.3.0 is not the vendored v7.2.1.  bngsim's dense
  difference quotient (jacobian="fd", or a law the closed form declines) reproduces
  CVODE 7.2's.  With this release it is still a valid quotient, but can differ from
  this CVODE's own in the last digits, and no CI leg tests the combination (issue #1335).
  

• docs/installation.md: the same note beside BNGSIM_USE_SYSTEM_SUNDIALS.

Verified by running the block in cmake -P with SUNDIALS_VERSION set to each of these:

version	warns
7.2.1	no
7.2.0	no
7.3.0	yes
6.7.0	yes
unset	yes

The default (vendored) configure does not enter the branch. Every CI leg uses the vendored path, so the full configure is exercised there.

A warning, not a refusal: a system SUNDIALS of another 7.x is supported (find_package(SUNDIALS 7.0 …)), and refusing it would drop HPC builds that link a site SUNDIALS.

Measured

test_difference_quotient_jump.py with the new cases, plus test_skip_audit.py, test_ci_run_list_coverage.py and test_docs_api_examples.py: 120 passed.

[elyfmiller]

Thanks. We reviewed this independently, Opus 5.5 and Fable 5.1, each on our own build of c4e8d10 (already on main c0877fa), on macOS arm64 with KLU on. We found no defect in what it changes:

• the configure warning works;
• the docs note is right;
• the new tests' oracles, closed forms and an exact 0 for the held species, are sound.

Three medium points are about what the new tests claim against what they run. Four low ones follow. What we checked is first, so it needn't be redone.

What we checked

check	result
The 24 new cases, and test_difference_quotient_jump.py	24 pass (and pass on main, as the body says). The file gives 58 passed. With the steady-state linear-solver tests and the skip audit: 101 passed.
Full python/tests on the merged build, BNGSIM_SKIP_AUDIT=strict	16,411 passed, 166 skipped, 0 failed
A real configure against a "system" SUNDIALS: the vendored 7.2.1 installed to a prefix, its SUNDIALSConfigVersion.cmake edited	reported 7.2.1: no warning; 7.3.0: the warning, word for word as in the body; 7.2.0: no warning. Configure succeeds in all three, and find_package sets SUNDIALS_VERSION itself.
The block run with cmake -P against the real VENDOR.json	the body's five cases as stated. 7.10.0 warns (the regex is anchored); 7.2 and 7.2.1-dev don't.
Upstream SUNDIALS v7.2.0 to v7.2.1	neither cvode_ls.c nor cvodes_ls.c changed, so comparing major.minor is the right granularity for that pair
Gates	changelog fragment 620 characters and well-formed; ruff clean; "build has no SuiteSparse/KLU" is a declared skip

1. The KINSOL threshold rows don't form KINSOL's quotient near the jump

method="newton" is two-tier by default (solve_by_newton_two_tier, steady_state.cpp:2563). A CVODE burst first carries the state into the basin of X* = 4 or 10, and KINSOL polishes from there in 2–3 steps. So the threshold rows test the burst plus a polish away from the threshold, not #1335's item 1.

The held-at-zero rows do reach KINSOL beside the held species.

Started on the threshold directly, KINSOL passes: steady_state_batch(states=[ic], method="newton", newton_burst=False), dense and forced-KLU, the six threshold starts plus gate and inhibitor. 16 of 16 converge to the exact root, in 8–13 steps from on or below the threshold and 2 from above. For s = 0.08, X₀ = 5, the default takes 3 steps and the direct start 8, both landing on X = 4.

We'd suggest the threshold rows use that call, so they pin KINSOL's own quotient.

2. "KINSOL, KLU" is the burst on KLU, then KINSOL dense

KINSOL's polish always builds a SUNDenseMatrix (steady_state.cpp:2299–2301); KLU exists only in the march. result.linear_solver is the march tier's answer (steady_state.cpp:6085), which types.hpp:1425 documents as describing the march.

So this row isn't "KINSOL on KLU", as the section comment and the body describe it. On all 8 models it gives the same result bits and the same counters as "KINSOL, dense", and differs only in linear_solver. That's because the burst's colored and dense quotients form one matrix (#1254). We'd suggest relabelling the row as the burst's route, or dropping it.

Underneath, and pre-existing: with newton_burst=False, the field still reports "klu" for a solve that built no sparse matrix.

3. The rows can't see a broken colored quotient

As pins of the march's colored quotient, they have little power:

• Negated (cvode_ss_colored_jac_body, rebuilt): 23 of 24 still pass; only gate-colored march fails.
• Minimum increment ×10⁶, so every column steps across the threshold: 24 of 24 pass.

The solver's counts do see the negated quotient, though. Against the restored build:

• RHS evaluations ×2.0–3.1 on the threshold rows;
• the gate goes from 183 RHS evaluations to 7,896, unconverged;
• the inhibitor ×45.

A per-row budget on n_steps or n_rhs_evals, say 3× the current counts, would catch what the final-state oracle can't. At least, the section comment should say the rows pin the root and the path rather than the quotient.

Low

4. The bound. The steady-state rows reuse OK = 1e-4, whose stated reason is about runs. All 24 roots are exact here, relative error 0 on every path (measured by both of us), so they need a reason of their own or a tighter bound.

5. Wording. The changelog and installation.md say a system SUNDIALS "other than the vendored release" warns, but the block compares major.minor, so 7.2.0 beside the vendored 7.2.1 is silent. That's right, given the upstream check above; the wording should match.

6. CMakeLists.txt:289. The system branch reads VENDOR.json with no EXISTS check, where the vendored branch checks.

   • A missing file is a raw file(READ) error, so an advisory warning becomes fatal for a packager that strips the vendored tree.
   • A tag-less file disables the check silently.
   • An unset version reads "System SUNDIALS is not…".

   if(EXISTS …) around the block would keep it best-effort.

7. Optional. A cross-reference from the fd Jacobian docs to the new installation note.

Not checked: a real system SUNDIALS of another release (we edited the version file instead), a mutation of KINSOL's own vendored difference quotient, and Linux or Windows beyond CI.

[elyfmiller]

Following up on our review above: we've pushed its items 1–7 as one commit on top, e699834, and are merging now that CI is green. The configure warning and the docs note stand as you wrote them. The commit changes what the tests run and assert, makes the check best-effort, and says major.minor in the docs.

What the commit changes

• test_difference_quotient_jump.py, the Tests: KINSOL's and the steady-state march's difference quotients on a jump, and SUNDIALS other than 7.2, are not covered #1335 section. Three paths, each still asserting the path it took:
  • "KINSOL" starts KINSOL on the given state: steady_state_batch(params=[{}], method="newton", newton_burst=False), dense. The default newton marched a burst into the root's basin first, so KINSOL's quotient was formed away from the jump (item 1).
  • "burst on KLU, then KINSOL" is the old "KINSOL, KLU" row, named for what it runs: a burst on the march's colored quotient, then KINSOL's polish, which always factors densely (item 2).
  • "colored march" is as before.
  • The march's rows also hold it to 350 RHS evaluations (item 3). On these rows it takes at most 236 with the right quotient, and at least 414 with it negated. The negated quotient still lands every threshold row on its root.
  • The roots are held to ROOT = 1e-12, with its reason (item 4): f is affine in X on each side of the threshold, so every root is exact.
• CMakeLists.txt (item 6). An if(EXISTS …) guard. Metadata it can't read, missing or without a tag, skips the check with a STATUS note instead of failing the configure or skipping silently. An unset system version is named "of an unreported release", and still warns.
• The changelog and installation.md say major.minor, and that a patch release of the vendored one doesn't warn (item 5).
• solvers.md links the "fd" paragraph to the installation note (item 7).

What we checked

check	result
Fail-before, with the march's colored quotient negated (cvode_ss_colored_jac_body, rebuilt)	The file with the commit: 8 failed, every colored-march row, on the budget (414–573 RHS evaluations on the threshold rows; the gate 7,896 and unconverged; the inhibitor 9,125). Your version: 1 failed, the gate.
KINSOL's own quotient negated (kinLsDenseDQJac, rebuilt)	The 8 direct KINSOL rows: none converges. The 8 burst rows: KINSOL fails, the solve falls back to integration, and the path assertion catches it.
The file on main 2eae53d + this PR + the commit	58 passed. The march's counts are unchanged under BNGSIM_CODEGEN_JIT=mir, under BNGSIM_CODEGEN_THRESHOLD=1, and under both.
The release check, run with cmake -P against the real VENDOR.json	As before: 7.2.1, 7.2.0, 7.2 and 7.2.1-dev silent; 7.3.0, 6.7.0, 7.10.0 and unset warn. A tag-less file: silent before, the note now. A missing file: a hard configure error before, the note and success now.
Full python/tests, BNGSIM_SKIP_AUDIT=strict, on main + this PR + the commit	16,834 passed, 143 skipped, 0 failed on one build. 16,397 passed and 433 skipped on the other build; its one failure is test_run_tests_sh's walk-up probe, which needs jax and diffrax, and its audit flags only pybind11, which that venv lacks.
Gates	ruff 0.14.0 check and format, ci/changelog.py check, git diff --check: clean
Overlap	Trial merges of the other 34 open PRs onto main + this PR, with and without the commit: the same either way. The one conflict is #1356 with main itself, since #1343 merged. #1379 shares solvers.md and merges clean.
CI	21 checks, all green, on e699834

Not changed here: with newton_burst=False and KLU forced, linear_solver still reports "klu" for a solve that built no sparse matrix (item 2's note). The new KINSOL rows run dense, so they don't lean on it. On a KLU-off build the march's rows, and their budget, skip.