[Montana]

Fixes #1308.

Stacked on #1321. This PR's first nightly failed Path-agreement on TIMEOUTs only: models at the edge of the 900 s cap (MODEL1006230049, 0053, 0077) lost every row once a check was added to every model. #1321 bounds each check by the model's cap and keeps a killed worker's finished rows. It is merged in here, and this PR's new checks run after every check main runs. Until #1321 merges, this PR also shows its commits.

Stacked on #1291. #1308 is blocked by #1291, which adds parity_checks/differential/, and this branch is built on it. Until #1291 merges, this PR also shows #1291's commits. Only the last commit, 156a116, is this PR's. Its diff is below.

What it adds

A jacobian_values check. At the three states _states already collects, it compares Model.jacobian(y, t) (where source == "analytical"; otherwise SKIP) with central differences of Model.rhs. Columns are scaled as sparsity_covers_jacobian scales them, and each entry is judged on its column's own peak.

How an entry is judged. It is differenced over a ladder of steps, from 1e-2 to 1e-10 of its column's scale. Each step's bound is the larger of its gaps to the neighbouring steps. An entry is a finding only where all of these hold:

• every resolvable step contradicts it. A step is resolvable where its bound is under 1% of the entry. To contradict, it must miss the analytic entry by more than 10 bounds + 1e-4 of the column's peak + 1e-8 of the matrix's;
• no step bears it out. A step "heard" at a 10% bound, which agrees with the analytic entry within 2 bounds, vetoes the finding;
• it is not on a kink: the one-sided differences agree at the finest resolvable step.

A non-finite analytic entry is a finding wherever some step resolves the entry.

Why each rule is there

The first drafts each produced false alarms on the corpus, and each rule answers one:

false alarm	what it was	the rule that answers it
blbr	f ≈ 1e11: at small steps the differences rounded alike to 2.0345 where the entry is 2	gaps to the steps a decade either side bound the rounding
mm_tqssa_stiff	a column of 2e-19 against entries of order 1 reads "100% off" on its own peak	the 1e-8-of-the-matrix floor
12 models at t = 0	species at 0 on max(0, x) / if() kinks: central differences give the mean of the two slopes	the kink test
BIOMD0000000498, BIOMD0000000480	saturations far narrower than the coarse steps: the coarse steps settle on the secant, the fine ones on the slope (= analytic)	every resolvable step must contradict, not the best one
BIOMD0000000393	#533's plateau: in sqrt(1+u)−1 with u ≈ 1e-6, a step under ulp·N/4 rounds the term's slope away, so the fine steps agree on 0.661 where the entry is −1.96e-4; the coarse steps (2–3% noisy) read −1.96e-4	the veto by a step heard at 10%

Results over the whole corpus

The run covered 1,353 rows: the 53 .net jobs and the BioModels differential_run.py reads, run locally.

outcome	rows
PASS	1,279
SKIP (no analytic Jacobian: Model.jacobian differences the RHS)	72
DIFF	2
EXCEPTION / TIMEOUT	0

The two DIFFs, triaged:

1. BIOMD0000000608, d x₅/d x₅ at t = 0 (analytic 0.1, differences −0.11). At t = 0, x₃ = x₅ = 0, and the rate laws divide by x₃ throughout (x_5/x_3, /x_3/x_3). The state is a 0/0 point, where the rates are not differentiable, and a single-step probe there reads 0.1004. Not a bngsim defect, and not something a difference can judge. It stays in the baseline as a known row.
2. MODEL2310250001, d GL_e_0/d _6PGL_0 at t = 100 (analytic 1.78e16, differences exactly 0 at every step). If Model.rhs had that slope, a step of 4e-5 would move f by about 7e11, far over its rounding, and it does not move. So Model.jacobian and Model.rhs disagree. The two analytic entries in that row stand in the ratio −S₁/S₀ (1.78e16 / −1.18e16 = −1.503). That is the shape of v8·S₁/Stot², but with Stot acting like 6.6e-9 rather than 1.0, and it is the same whichever of rhs or jacobian is called first. The state lies past the blow-up both engines share (the job's known_artifact: "BOTH engines blow up to ~−5e16"), so it is a low-priority lead rather than a filed issue.

Power, on the 53 .net rows. The largest entry of the analytic Jacobian was corrupted, and the check run on the result:

corruption	caught	missed	no effect (entry already 0)
×1.1	46	0	7
×1.01	46	0	7
zeroed	46	0	7
a spurious nonzero entry	53	0	0

The seven "no effect" rows are the pulse and constant models, whose Jacobian is exactly zero (their rates read only the clock).

Tests (python/tests/test_differential_shapes.py)

• The shape tests pick up the new check automatically, as for every check in CHECKS, and all 15 shapes pass it.
• Negative controls, which must read DIFF:
  • a 10% wrong entry (named, and its size right);
  • a 1% wrong entry;
  • a spurious entry where the RHS reads nothing;
  • a non-finite entry;
  • a real shape (derived-log10) with its largest entry 10% off. The test asserts the corruption changes the matrix, since the triangle pulses' Jacobians are all 0.
• False-alarm controls, which must not read DIFF:
  • a differenced Jacobian, which returns None;
  • rounding in a large f;
  • a near-zero column at a steady state, where only the matrix-wide floor helps;
  • a kink, −|b−2| at b = 2;
  • a feature narrower than the coarse steps, 1e-9·tanh(b/1e-9), with a constant that makes the coarse plateau the tightest;
  • a slope rounded away at the fine steps (The FD self-check rejects BIOMD0000000393's correct analytical Jacobian: its rate laws' sqrt(1+ε)−1 term is accurate to 8e-10 in doubles, and the gate's two steps sit on the same wrong plateau #533's shape).

Mutations. Each guard was mutated, and every mutation fails its control:

• no matrix-wide floor;
• tolerance 10% instead of 1e-4;
• non-finite entries not flagged;
• a differenced Jacobian judged;
• no kink filter;
• no veto;
• a veto at ten bounds instead of two.

The one survivor is judging only the best step instead of every resolvable one. With the veto in place, the two differ only when a resolvable step is between 2 and 10 bounds from the analytic entry, so I kept the stricter rule. I also removed a rounding term (4·eps·|f|/h) that the ladder's gaps made redundant. Without it, every control and every hard corpus case reads the same.

Suites: test_differential_shapes.py and parity_checks/tests/test_differential_run.py: 194 passed, 1 xfailed (#1291's #1254 xfail). ruff 0.14.0 is clean. The branch merges cleanly with the current lanl main.

Known limits

• A feature narrower than any step that rounding lets the ladder resolve can still read as a finding. One example is a tanh of width 1e-10 under a constant of 10.
• So can rounding lost inside a law deeper than even the coarse steps can speak to. One example is The FD self-check rejects BIOMD0000000393's correct analytical Jacobian: its rate laws' sqrt(1+ε)−1 term is accurate to 8e-10 in doubles, and the gate's two steps sit on the same wrong plateau #533's shape at u = 1e-11 with nothing else in the law.
• Both are documented in the controls' docstrings.

Baseline

Recorded from this PR's own nightly run, 40448597362, in 8920baa, as #1308 asks. It holds the 1,353 jacobian_values rows that run added:

• 1,278 PASS;
• 68 SKIP (no analytic Jacobian);
• 5 TIMEOUT (the corpus's largest models);
• the two DIFFs triaged above.

CI matched the local run. Against the committed baseline the run had no regressions: no row dropped, and no existing row's outcome changed.

Changes (this PR's commit only)

diff --git a/parity_checks/differential/_differential.py b/parity_checks/differential/_differential.py
index 5bf10534..5d9fe08a 100644
--- a/parity_checks/differential/_differential.py
+++ b/parity_checks/differential/_differential.py
@@ -100,6 +100,7 @@ CHECKS = (
     "locale",
     "grid",
     "sparsity_covers_jacobian",
+    "jacobian_values",
     "backend_as_asked",
 )
 
@@ -521,6 +522,178 @@ def check_sparsity(base: Base) -> Verdict:
     return _guard(go)
 
 
+# jacobian_values (issue #1308). Each entry is differenced over a ladder of
+# steps, JAC_STEPS of its column's scale, and each step's difference carries
+# its own error bound: the larger of its gaps to the neighbouring steps', which
+# moves with truncation at the coarse end and with rounding at the fine end,
+# the steps being a decade apart. A step is resolvable where
+# that bound is under JAC_AGREE of the entry, and has a say where it is under
+# JAC_HEARD. The entry is off only where the analytic one contradicts the
+# difference at every resolvable step, by more than ten times that step's bound
+# plus JAC_TOL of its column's peak and JAC_FLOOR of the whole matrix's, and no
+# step with a say bears it out (within twice its bound).
+#
+# A say for a step a little noisier than a resolvable one, because rounding
+# can make a fine step look resolved and be wrong: in sqrt(1 + u) - 1 with
+# u ≈ 1e-6 (BIOMD0000000393's quadratic roots, #533), a step under
+# ulp(1)·N/4 leaves 1 + u as it was, so the term's whole slope rounds away and
+# the fine steps agree, from 1e-8 to 1e-10, on 0.661 where the entry is
+# -1.96e-4. The coarse steps, 2-3% apart from one another, read -1.96e-4. No
+# bound on rounding from |f| could see that loss, which happens inside the law.
+#
+# Every step, and not the best one: a feature narrower than the coarser steps
+# (a saturation k·s/(Km + s) with a Km far under the species' scale) holds two
+# plateaus, the secant across it at coarse steps and the slope at the point at
+# fine ones, and the analytic entry is the second. BIOMD0000000498's dS26/dS26
+# is -13.74 to 1e-6 and -13.11 at 1e-8, the analytic -13.11; BIOMD0000000480's
+# s1 sits at -1e-14 under a saturation with Km 1e-3, its analytic entry 2,000
+# off the coarse secant. A wrong entry is wrong at every scale.
+#
+# Rounding in a large f is what the gaps bound: blbr's f is 1e11, and at a
+# step of 1e-6 its differences came out 2.0345 at two steps a factor of two
+# apart, alike, where the entry is 2; a decade either side they do not. The
+# floor is what a near-zero column costs: mm_tqssa_stiff's carries
+# 2e-19 against entries of order 1, which on its own peak reads 100% off.
+JAC_STEPS = tuple(10.0**-k for k in range(2, 11))
+JAC_AGREE = 1e-2
+JAC_HEARD = 1e-1
+JAC_SMOOTH = 1e-3
+JAC_TOL = 1e-4
+JAC_FLOOR = 1e-8
+
+
+def jacobian_misfits(model, states, t=0.0, steps=JAC_STEPS) -> list[dict] | None:
+    """The entries of the analytic Jacobian that central differences of the
+    right-hand side contradict, at any of ``states``; ``None`` when
+    ``Model.jacobian`` does not give the analytic one there (nothing of the
+    closed form to judge).
+
+    Each column is scaled by the state's own size, as :func:`sparsity_misses`
+    scales it, and judged on its own peak (the larger of the analytic one and
+    the difference's), so an entry in units of one species is never held to
+    another's. An entry the differences cannot resolve at any step, or one on a
+    kink, is not judged (see the constants above). A non-finite analytic entry
+    is off wherever some step resolves the entry.
+    """
+    out: dict[tuple[int, int], dict] = {}
+    for y in states:
+        y = np.asarray(y, dtype=float)
+        n = y.size
+        if n == 0:
+            continue
+        J = model.jacobian(y, t)
+        if getattr(J, "source", None) != "analytical":
+            return None
+        big = float(np.abs(y).max())
+        scale = np.maximum(np.abs(y), max(big, 1.0) * 1e-3)
+        f0 = np.asarray(model.rhs(y, t))
+        cen, fwd, bwd = [], [], []
+        for h_rel in steps:
+            c_, f_, b_ = [], [], []
+            for j in range(n):
+                h = h_rel * scale[j]
+                up, dn = y.copy(), y.copy()
+                up[j] += h
+                dn[j] -= h
+                f_up, f_dn = np.asarray(model.rhs(up, t)), np.asarray(model.rhs(dn, t))
+                c_.append((f_up - f_dn) / (2 * h))
+                f_.append((f_up - f0) / h)
+                b_.append((f0 - f_dn) / h)
+            sc = scale[np.newaxis, :]
+            cen.append(np.stack(c_, axis=1) * sc)
+            fwd.append(np.stack(f_, axis=1) * sc)
+            bwd.append(np.stack(b_, axis=1) * sc)
+        cen, fwd, bwd = map(np.stack, (cen, fwd, bwd))  # (n_steps, n, n)
+        finite = np.isfinite(cen)
+        k_steps = len(steps)
+        # The truncation bound: the larger gap to a neighbouring step.
+        trunc = np.full_like(cen, np.inf)
+        for k in range(k_steps):
+            gaps = [np.abs(cen[k] - cen[m]) for m in (k - 1, k + 1) if 0 <= m < k_steps]
+            trunc[k] = np.max(np.stack(gaps), axis=0)
+        err = np.where(finite, trunc, np.inf)
+        a = np.asarray(J, dtype=float) * scale[np.newaxis, :]
+        mag = np.where(finite, np.abs(cen), 0.0)
+        peak = np.maximum(
+            mag.max(axis=(0, 1)), np.where(np.isfinite(a), np.abs(a), 0.0).max(axis=0)
+        )
+        top = float(peak.max()) if peak.size else 0.0
+        size = np.maximum(mag, JAC_FLOOR * peak[np.newaxis, np.newaxis, :])
+        resolvable = finite & (err <= JAC_AGREE * size)
+        any_res = resolvable.any(axis=0)
+        # A kink in the state (max(0, x) at x = 0, an if() on a species on its
+        # threshold) is the #368 trap: the central differences converge, at
+        # every step, to the mean of the two one-sided slopes, which no analytic
+        # entry is. The one-sided differences tell it at the finest resolvable
+        # step: they differ there by the slope's jump, and by O(h) where the law
+        # is smooth. MODEL2105110001's NEU sits at -4e-216 on a kink whose sides
+        # are 0.3% apart, and its analytic entry is the side it is on.
+        finest = np.where(any_res, k_steps - 1 - np.argmax(resolvable[::-1], axis=0), 0)
+        at = finest[np.newaxis]
+        fwd_f, bwd_f, size_f, err_f = (
+            np.take_along_axis(x, at, axis=0)[0] for x in (fwd, bwd, size, err)
+        )
+        kink = np.abs(fwd_f - bwd_f) > JAC_SMOOTH * size_f + 4.0 * err_f
+        slack = JAC_TOL * peak[np.newaxis, :] + JAC_FLOOR * top
+        contradicts = np.abs(a[np.newaxis] - cen) > 10.0 * err + slack[np.newaxis]
+        everywhere = np.all(~resolvable | contradicts, axis=0)
+        # A noisier step's say is a veto, and only where it bears the analytic
+        # entry out, within twice its own bound: one that merely fails to
+        # contradict it at ten times that bound would excuse any error.
+        heard = finite & (err <= JAC_HEARD * size)
+        vetoed = np.any(
+            heard & (np.abs(a[np.newaxis] - cen) <= 2.0 * err + slack[np.newaxis]), axis=0
+        )
+        off = any_res & (~np.isfinite(a) | (~kink & everywhere & ~vetoed))
+        best = np.argmin(np.where(resolvable, err, np.inf), axis=0)
+        d_best = np.take_along_axis(cen, best[np.newaxis], axis=0)[0]
+        for i, j in zip(*np.nonzero(off), strict=True):
+            finite_a = np.isfinite(a[i, j]) and peak[j] > 0
+            e = float(abs(a[i, j] - d_best[i, j]) / peak[j]) if finite_a else float("inf")
+            key = (int(i), int(j))
+            if key not in out or e > out[key]["off"]:
+                out[key] = {
+                    "entry": key,
+                    "t": float(t),
+                    "analytic": float(J[i, j]),
+                    "difference": float(d_best[i, j] / scale[j]),
+                    "off": e,
+                }
+    return sorted(out.values(), key=lambda m: -m["off"])
+
+
+def check_jacobian_values(base: Base) -> Verdict:
+    def go():
+        model = load(base.job.path)
+        states = _states(base.job)
+        worst: dict[tuple[int, int], dict] = {}
+        for t, y in states:
+            found = jacobian_misfits(model, [y], t=t)
+            if found is None:
+                return Verdict(
+                    "skip", note="no analytic Jacobian: Model.jacobian differences the RHS"
+                )
+            for m in found:
+                if m["entry"] not in worst or m["off"] > worst[m["entry"]]["off"]:
+                    worst[m["entry"]] = m
+        if worst:
+            names = list(model.species_names)
+            ranked = sorted(worst.values(), key=lambda m: -m["off"])
+            shown = ", ".join(
+                f"d{names[m['entry'][0]]}/d{names[m['entry'][1]]} {m['off']:.1e}"
+                for m in ranked[:5]
+            )
+            return Verdict(
+                "diff",
+                ranked[0]["off"],
+                f"{len(ranked)} analytic entries off the differences (of their column's peak): {shown}",
+                {"off": [{**m, "entry": list(m["entry"])} for m in ranked[:50]]},
+            )
+        return Verdict("pass", 0.0, f"every judged entry agrees at {len(states)} states")
+
+    return _guard(go)
+
+
 def comma_locale() -> str | None:
     """A comma-decimal LC_NUMERIC this machine has, or None. Asked of a child
     process, so this one's locale is never set."""
@@ -593,6 +766,7 @@ def check_all(job: Job, locale_name: str | None, checks=CHECKS) -> dict[str, Ver
         "cc_vs_mir": check_cc_vs_mir,
         "grid": check_grid,
         "sparsity_covers_jacobian": check_sparsity,
+        "jacobian_values": check_jacobian_values,
         "backend_as_asked": check_backend,
     }
     for name, fn in one.items():
diff --git a/parity_checks/differential/differential_run.py b/parity_checks/differential/differential_run.py
index 02cc1645..acc5b8d8 100644
--- a/parity_checks/differential/differential_run.py
+++ b/parity_checks/differential/differential_run.py
@@ -5,8 +5,10 @@ report (issue #1257).
 The checks are ``_differential.py``'s: a sensitivity run's state against a plain
 run's, the FD Jacobian against the analytic one, KLU against the dense solver,
 compiled against interpreted, cc against MIR, a comma-decimal locale against C,
-and three rules for bugs whose answers are right (a sensitivity run refused, a
-sparsity pattern that misses entries, a backend other than the one asked for).
+and four rules for bugs whose answers are right (a sensitivity run refused, a
+sparsity pattern that misses entries, analytic Jacobian entries that differences
+of the right-hand side contradict (#1308), a backend other than the one asked
+for).
 One row per (model, check), its ``method`` the check:
 
     the path agrees with its reference .............. PASS

Plus 198 lines of tests in python/tests/test_differential_shapes.py. No changelog fragment: this changes parity_checks/ and tests only, which ci/changelog.py does not ask an entry for.