[Montana]

Closes #1257, scoped as the maintainer review asked (#1257 (comment)).

What it does

A recurring bngsim bug is two code paths that should give the same answer and don't (#1182, #1227, #1228, #1235, #1254). This runs those probes for every model, every time:

check	path	reference
sens_vs_plain	the state of a sensitivity run	a plain run
ajac_vs_fdjac	jacobian="fd"	jacobian="auto"
klu_vs_dense	force_sparse_linear_solver=True	the dense solver, on the same Jacobian
cc_vs_interp	codegen=True	codegen=False
cc_vs_mir	BNGSIM_CODEGEN_JIT=mir	cc (on a MIR build)
locale	the run under de_DE.UTF-8, in a child process	under C

It also applies three rules for bugs whose answers are right, which no comparison of answers sees:

• sens_runs: the sensitivity run doesn't complete where the plain run does (The issue #952 hold check refuses a sensitivity run at a bend: a continuous narrow pulse (max/abs or if()) is refused at rtol <= 1e-10, blamed on a jump it does not have #1229). Reported as UNSUPPORTED; the baseline turns a new one into a finding.
• sparsity_covers_jacobian: every entry that central differences of the RHS find nonzero has a slot in jacobian_sparsity (Jacobian sparsity has no edges for a rate_of__ token, so a written-out rate declines the analytic Jacobian (#1124) #1239).
• backend_as_asked: codegen=False runs the interpreter and codegen=True runs compiled code (On a MIR build, Simulator(method='ssa', codegen=False) still JITs the propensities: the C++ selector reads only BNGSIM_SSA_NO_CODEGEN #1213).

Where it lives, per the review

• One finite-difference oracle. fd_sens gains the tolerance sweep: an analytic column that still moves between rtol 1e-10 and 1e-12 where the trajectory settles is a DIFF. It also gains eight pulse models over the clock (max/abs, clamp, if(); narrow w = 0.01 and wide w = 0.5). Of the proposal's five rules, fd_sens already had three (the relative step, step agreement, noise floors). The sweep supplies rules 3 and 5, and it judges only the columns fd_sens's floors judge.
• Shared checks: parity_checks/differential/_differential.py.
• Shapes as pytest: python/tests/test_differential_shapes.py, 158 tests in about 7 s, on every PR. It includes the shapes the proposal promised to add: an event at a sample time, a moving compartment, rateOf across compartments, and an event-set pulse train. Known open bugs are strict xfails naming their issue: 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 (KLU with FD on an equality) and log10(p) or log2(p) in a derived rate constant declines the analytic sensitivity RHS for the whole model, which falls back to the difference quotient #1214 (log10 declines the analytic sensitivity RHS).
• Corpus as nightly:
  • differential_run.py writes a _core report, one row per (model, check), compared by compare.sh against baselines/differential.json and the previous night.
  • The corpus is fd_sens's 53 models and the BioModels under 1 MB, from the corpus job's cache.
  • The job runs locale-gen de_DE.UTF-8. The runner checks for a comma-decimal locale once and exits 2 without one, so the missing report alerts.

How the oracle stays careful

Every false positive seen so far came from the checking, not from bngsim. A disagreement counts only past max(1e3·rtol, 10 × the base's own drift at a tenfold tighter tolerance). Each cell is judged on the trajectory's own scale, past an absolute floor of the run's atol.

Running over the larger corpus first, as asked, changed four things:

1. The sweep judges only columns above the solver's resolution. It first flagged mwc (its phi column) and the narrow pulses' position and width. Those columns are truly 0 and hold only the integration's error.
2. A base not converged to 0.1% at a tighter tolerance is no reference, and its comparisons are skipped, saying why. MODEL2310250001 drives a species through a pole: its plain run moves 2.6% at a tenfold tighter tolerance, and libRoadRunner gives −1.4 to −4.6 at three tolerances.
3. KLU is compared with dense on the same Jacobian. Against dense with the analytic Jacobian, an FD-Jacobian weakness read as a KLU one.
4. A model that doesn't load, or whose plain run doesn't complete, has no reference. That belongs to rr_ode and sbml_semantic, so it's SKIP here; a PASS turning into SKIP still alerts.

The pytest file holds the gate to account both ways. Lorenz over t = 25 has an FD-Jacobian path 3.4e-4 from the base, past the 1e-5 that 1e3·rtol allows, and both are right. With the gate it isn't a finding; without it, it is.

Tests

• Shapes: 134 tests in python/tests/test_differential_shapes.py:
  • every shape through every check;
  • the sensitivity RHS taken;
  • negative controls: a path off by more than its bound, a sparsity pattern missing a slot (on a fake model), a lying backend, a refusal, a trajectory smaller than atol, a shape or non-finite mismatch, and the conditioning gate.
• Runner: 11 tests in parity_checks/tests/test_differential_run.py: rows and outcomes, a model with no rows, the manifests, the BioModels size cap, the locale stop, an unknown check.
• Sweep: 7 tests in parity_checks/tests/test_fd_sens.py: a moving column, an unsettled trajectory, the rounding floor, judged columns only, a converged column, and how the sweep's verdict applies.

Mutations. The library: 9 of 9 caught.

mutation	caught by
compare always passes	the controls
no conditioning gate	the Lorenz test
no atol floor	the small-trajectory test
shape mismatch not caught	the mismatch test
KLU never forced	the #1254 strict xfail, which then passes
sparsity ignores the pattern	the fake-model test
backend never checked	the backend test
sensitivity RHS always reported analytic	the #1214 strict xfail
a refusal read as a crash	the refusal tests

The sweep: 4 of 4 caught (no settled gate, judged columns ignored, no rounding floor, sweep never applied).

Negative control on #1227. Against the commit just before #1242, fd_sens with the sweep flags pulse_triangle_clamp_narrow (k_clock moves 107%) and pulse_ramp_clamp_narrow (1.15e-3). Without the sweep, all 16 pulse rows pass. On main, all 16 pass with it.

Suites. parity_checks/tests: 921 passed, 18 skipped. The CI-config tests, the skip audit and the shapes: 295 passed, 2 xfailed. ruff is clean on everything touched.

False positives on the larger corpus

3,437 models, 30,933 rows: every .net in benchmarks/models/net, the 214 SBML event models, the 1,823 semantic cases and the 1,323 BioModels. 23,361 PASS, 7,288 SKIP, 217 UNSUPPORTED, 36 TIMEOUT (4 models at the 300 s cap), 12 EXCEPTION, 19 DIFF.

One model is a false positive: MODEL7743528808, 2 rows. Its paths are 2.15e-3 apart against a bound of 2.1e-3, and the difference is gone at a tenfold tighter tolerance. I left the gate as it is rather than fit it to one model. The case would sit in the baseline and never alert.

What it found on main

finding	models	tracked as
comma-locale misreads	13	#1235, fixed by #1238 (merged)
jacobian="fd" alone refused: "force_sparse_linear_solver: … Drop the flag", a flag never set. The automatic routing picks KLU, and the pattern has no nonzeros to color.	semantic 01466, 01530–01533, 01564	#1296, fix in #1292
KLU's colored FD Jacobian steps each column by at least 1.5e-8 in absolute terms, so it fails or reads 3% off where dense FD and KLU analytic both run	4 BioModels, and #1254's equality	#1254, fix in #1295
the answer changes with the output grid: 8× the points gives 5.3e7 of the peak off the converged run, with no error (the dense FD run, which first read 2.3e17, now fails loudly)	BIOMD0000000810	#1297
a sensitivity run's state departs about 1% from a converged plain run, with no events	MODEL2001200002	#1298
codegen copies an ExprTk implicit product into C	BIOMD0000001042	#1212

Since opening

• A grid check. 8× the output points, read at the base's own sample times, must give the base. On the commit before stops: place a schedule's stops stretch by stretch between the events that set it (#1244) #1263 it finds A pulse schedule whose width an event writes gets no stops: X(240) = 1 for 17, at every tolerance #1244's pulse train 7.5× the trajectory's peak off; on main it passes. Shapes now include the width set inside a pulse.
• cc against MIR runs. The shapes run on mir.yml's legs, the one place a MIR build exists, which also exercises On a MIR build, Simulator(method='ssa', codegen=False) still JITs the propensities: the C++ selector reads only BNGSIM_SSA_NO_CODEGEN #1213's check (codegen=False on a MIR build). A missing C compiler is a skip there, not a crash.
• 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's xfail holds only where KLU is built. The MIR legs build without it.
• The tolerance sweep runs on every PR, on the four narrow pulses (parity_checks/tests/test_fd_sens.py), not only nightly.
• A 900 s cap per model in the nightly (300 s let nine models flip with runner speed, compiling cold); the remaining 5 timeouts are the same every night.
• Baselines from this PR's own run 40390175020: 13,530 differential cases and 106 fd_sens cases, each comparing clean against it.
• log10(p) or log2(p) in a derived rate constant declines the analytic sensitivity RHS for the whole model, which falls back to the difference quotient #1214's xfail is removed, since struct: differentiate a derived parameter over log10, log2 and ExprTk's other calls (#1214) #1224 fixed it. The test-count floors follow the added tests.

Compatibility

• Nothing in python/bngsim changes, so no changelog fragment is required (ci/changelog.py required).
• ci/test_floors.env: the floors move with the added tests, as the file asks. Python is 16,137 on main, 16,271 here; floor 16,261. Parity is 921 on main, 939 here; floor 929.
• Baselines: baselines/differential.json comes from this PR's own nightly run (ALLOW_MISSING_BASELINE is set on pull requests), and fd_sens.json is refreshed from it for the eight new models. Until then they are listed as new, not alerted.

Performance

• Pytest: 134 shape tests, about 7 s.
• Nightly: about 3,400 s of worker time on an M3 Max for the job's corpus. On a 4-core runner that should be roughly half an hour, inside the job's 180-minute limit.
• fd_sens: two more sensitivity runs per job for the sweep.

Remaining risks

• The gate's constants (1e3·rtol, 10× drift, 1e-3 unconverged) were measured on these 3,437 models. A corpus of another kind could show another false positive. If so, it would land in the baseline, not alert.
• The refusal text tests. _sparse_declined and codegen's "declined this model" are recognized by their text: bngsim raises no type of its own for them. A reworded message would turn a SKIP or UNSUPPORTED into an EXCEPTION, which alerts.
• MIR: cc_vs_mir runs only on a MIR build; the nightly's runner has none, so it reads SKIP there.
• Platform: checked on macOS arm64. The first CI night on Linux is the first run there.

Maintainer review: the three strongest objections

1. "A path comparison with a tolerance gate will either miss bugs or flag noise." Every constant was set from a measured false positive, and each check has a test that it fails on the defect it is there to see (9 of 9 library mutations caught). On 3,437 models it reported one false positive. It also found three bugs no one had filed and two that were already open.
2. "This duplicates fd_sens." It doesn't touch sensitivities against differences: those stay in fd_sens, which gains the sweep. The new suite compares bngsim's paths with each other, which no existing suite does.
3. "Classifying refusals by message text is fragile." Agreed. The classifier is used only where bngsim has no exception type for the refusal, and a misread message fails loudly (EXCEPTION alerts) rather than silently. A SparseJacobianUnavailable type and a codegen-declined type would make it exact. That's a small follow-up I can open.

Diff (fd_sens_run.py, the one existing runner this changes; the new files are listed above)

diff --git a/parity_checks/fd_sens/fd_sens_run.py b/parity_checks/fd_sens/fd_sens_run.py
index d358cac6..06b8a989 100644
--- a/parity_checks/fd_sens/fd_sens_run.py
+++ b/parity_checks/fd_sens/fd_sens_run.py
@@ -106,6 +106,26 @@ REFERENCE_ENGINE = "finite_difference"
 # starts at the child's spawn, before the worker's.
 CAP_SHARE = 0.9
 
+# The tolerance sweep (issue #1257): the analytic run again at two tolerances
+# far tighter than the suite's. An analytic column converges with the
+# tolerance; a column on CVODES' difference quotient that does not converge
+# still moves there, which is how the narrow clamp pulse of #1227 showed (5x
+# off, and a difference of plain runs at one tolerance did not see it). Judged
+# only where the trajectory itself settles between the two: a model that
+# amplifies every difference moves its columns with the tolerance too, and no
+# oracle can judge those (LV).
+SWEEP_RTOLS = (1e-10, 1e-12)
+# atol at each, relative to its rtol
+SWEEP_ATOL_RATIO = 1e-2
+# The trajectory settled: the two runs' states this close on their peak.
+SWEEP_SETTLED = 1e-6
+# A column that moves more than this on its own peak between the two has not
+# converged.
+SWEEP_TOL = 1e-3
+# A column under this of the state's peak is zero to the arithmetic, and is
+# not asked to converge: erk_activation's 3.9e-23 beside a state of 3e6.
+SWEEP_ROUNDING = 1e-14
+
 
 def load_model(path: Path):
     import bngsim
@@ -359,6 +379,79 @@ def column_verdicts(analytic, reference, gap, x_max: float, values, atol: float,
     return out
 
 
+def tolerance_sweep(runs, params, judged=None) -> dict:
+    """The sweep's finding on ``runs``, ``[(species, sensitivities)]`` at
+    :data:`SWEEP_RTOLS` (loosest first): ``settled`` (the trajectory converged
+    between them), ``state_gap``, and per parameter column of ``judged`` (all of
+    ``params`` by default) above the state's rounding, how far it moved on its
+    own peak (``column_gaps``). ``moving`` names the columns past
+    :data:`SWEEP_TOL`, only where the trajectory settled."""
+    judged = set(params if judged is None else judged)
+    (x1, s1), (x2, s2) = (
+        (np.asarray(x, dtype=float), np.asarray(sx, dtype=float)) for x, sx in runs
+    )
+    if x1.shape != x2.shape or s1.shape != s2.shape:
+        return {"settled": False, "state_gap": float("inf"), "column_gaps": {}, "moving": []}
+    peak = float(np.nanmax(np.abs(x2))) if x2.size else 0.0
+    state_gap = float(np.nanmax(np.abs(x1 - x2)) / peak) if peak > 0.0 else 0.0
+    settled = bool(np.isfinite(state_gap) and state_gap <= SWEEP_SETTLED)
+    gaps = {}
+    for j, p in enumerate(params):
+        if p not in judged:
+            continue
+        col = s2[..., j]
+        size = float(np.nanmax(np.abs(col))) if col.size else 0.0
+        if not size > SWEEP_ROUNDING * max(peak, 1.0):
+            continue
+        gaps[p] = float(np.nanmax(np.abs(s1[..., j] - col)) / size)
+    moving = [p for p, g in gaps.items() if settled and not g <= SWEEP_TOL]
+    return {"settled": settled, "state_gap": state_gap, "column_gaps": gaps, "moving": moving}
+
+
+def _sweep_runs(path: Path, params: list[str], spec: dict):
+    import bngsim
+
+    out = []
+    for rtol in SWEEP_RTOLS:
+        model = load_model(path)
+        run = bngsim.Simulator(
+            model,
+            method="ode",
+            sensitivity_params=list(params),
+            sensitivity_method=spec["sens_method"],
+        ).run(
+            t_span=(spec["t_start"], spec["t_end"]),
+            n_points=spec["n_points"],
+            rtol=rtol,
+            atol=rtol * SWEEP_ATOL_RATIO,
+            max_steps=spec.get("max_steps", 100_000),
+        )
+        out.append((run.species, run.sensitivities))
+    return out
+
+
+def _apply_sweep(res: dict, sweep: dict) -> None:
+    """The sweep's verdict over the difference's, into ``res``: a column that
+    does not converge with the tolerance where the trajectory does is a DIFF,
+    whatever the difference said (it may pass, or have been unable to judge)."""
+    res["sweep_settled"] = sweep["settled"]
+    res["sweep_state_gap"] = sweep["state_gap"]
+    res["sweep_worst"] = max(sweep["column_gaps"].values(), default=0.0)
+    if sweep["moving"] and res.get("status") in ("pass", "bad_test", "diff"):
+        worst = max(sweep["column_gaps"][p] for p in sweep["moving"])
+        res["columns_diff"] = sorted(set(res.get("columns_diff") or []) | set(sweep["moving"]))
+        res.update(
+            status="diff",
+            value=max(float(res.get("value") or 0.0), worst),
+            why=(
+                f"the analytic columns {', '.join(sweep['moving'])} still move by "
+                f"{worst:.2e} of their peak between rtol {SWEEP_RTOLS[0]:g} and "
+                f"{SWEEP_RTOLS[1]:g}, where the trajectory settles to "
+                f"{sweep['state_gap']:.1e}: they do not converge with the tolerance"
+            ),
+        )
+
+
 def _judge(res: dict, analytic, reference, gap, species, values, spec: dict) -> None:
     """The verdict on two same-shaped tensors, into ``res``: which columns are
     judged, whether the reference converged where the verdict looks, then
@@ -512,6 +605,20 @@ def _worker(spec: dict, q) -> None:
         else:
             species = np.asarray(run.species, dtype=float)
             _judge(res, analytic, reference, gap, species, values, spec)
+        # The sweep asks of the columns the verdict judged, those above the
+        # solver's resolution: a column that is 0 to the integration (the
+        # position and width of a pulse whose area holds) reads its noise, and
+        # its noise moves with the tolerance.
+        judged = [p for p in params if p not in (res.get("columns_not_judged") or ())]
+        if (
+            spec.get("sweep", True)
+            and judged
+            and res.get("status") in ("pass", "bad_test", "diff")
+        ):
+            try:
+                _apply_sweep(res, tolerance_sweep(_sweep_runs(path, params, spec), params, judged))
+            except Exception as e:  # noqa: BLE001 - a sweep that does not run judges nothing
+                res["sweep_error"] = f"{type(e).__name__}: {str(e)[:200]}"
     except Exception as e:  # noqa: BLE001 - one model's failure is that row's EXCEPTION
         captured = capture_exception("bngsim", e)
         res.update(status="exception", exception=captured.text, exception_class=captured.cls)
@@ -549,6 +656,7 @@ def make_specs(jobs: list[dict], args) -> list[dict]:
                     # atol is absolute, and the atol/|p| floor with it.
                     "atol": float(job.get("atol", args.atol)),
                     "fd_tol": args.fd_tol,
+                    "sweep": not getattr(args, "no_sweep", False),
                     "param_cap": int(job.get("param_cap", args.param_cap)),
                     "sens_method": method,
                     "timeout": float(job.get("timeout", args.timeout)),
@@ -584,6 +692,10 @@ def to_result(spec: dict, r: dict, stamp: dict) -> JobResult:
         "resolvable_columns": r.get("resolvable_columns"),
         "columns_diff": r.get("columns_diff"),
         "sens_rhs_decline_reason": r.get("sens_rhs_decline_reason"),
+        "sweep_settled": r.get("sweep_settled"),
+        "sweep_state_gap": r.get("sweep_state_gap"),
+        "sweep_worst": r.get("sweep_worst"),
+        "sweep_error": r.get("sweep_error"),
         **(r.get("verdict") or {}),
     }
     return JobResult(
@@ -623,6 +735,7 @@ def main(argv=None) -> int:
         "is declared unconverged (BAD_TEST)",
     )
     ap.add_argument("--param-cap", type=int, default=DEFAULT_PARAM_CAP)
+    ap.add_argument("--no-sweep", action="store_true", help="leave the tolerance sweep out")
     ap.add_argument("--models", default="", help="comma-separated model_id substrings")
     args = ap.parse_args(argv)
 
@@ -659,6 +772,7 @@ def main(argv=None) -> int:
             "fd_atol": FD_ATOL,
             "fd_rel_step": FD_REL_STEP,
             "fd_tol": args.fd_tol,
+            "sweep_rtols": None if args.no_sweep else list(SWEEP_RTOLS),
             "param_cap": args.param_cap,
         },
         "hardware": asens.hardware_info(),

[wilton-williams]

This rocks dude!