[Montana]

Fixes #1309.

Stacked on #1291, at its review head 2259a57, so it builds on the forced-dense klu_vs_dense. Until #1291 merges, this PR also shows #1291's commits. Only the last commit, 3f62c03, is this PR's. It adds one row to CHECKS (stable), as #1315, #1316 and #1318 do, so whichever lands second keeps both. I'll resolve it.

What it adds

A DIFF says two paths disagree, not which one is wrong, and a small own drift does not make an answer uncontested. On #1304's BIOMD0000000338, KLU moved 1e-7 at a tenfold tighter tolerance while dense and libRoadRunner sat 3.5% away.

Attribution. Every DIFF whose comparison kept its two trajectories goes to third arms. Verdict now carries the two trajectories, but they are never written to a report. The arms run once per model, however many DIFFs it has:

arm	what it runs
tighter	the base's own route at 100× tighter rtol and atol
dense / klu	the base on the other linear solver (dense for a KLU base, KLU otherwise)
roadrunner	libRoadRunner, for SBML, at the job's tolerance (species read as [id], concentrations, as rr_parity selects)

Each arm sides with the comparison's reference, its path, both, or neither, judged against the same bound as the DIFF. The row's extra["attribution"] records each arm's side and distances. Its verdict is the one side every arm that took a side took, or split. An arm that cannot run (libRoadRunner giving up, KLU declined) is left out.

Branch stability. A new stable check reruns the base:

• at n_points − 1 and + 1, read at the last point, which every grid shares;
• at twice the points, read at the base's;
• at rtol × 0.5 and × 2.

A jump past the bound is a DIFF with extra["contested"].

The report lists contested models in meta["contested"]: a model whose stable row reads DIFF, or one with a DIFF whose arms split (D.contested). The runner prints them.

#1304's two models, as the issue asks

model	klu_vs_dense	attribution	arms
BIOMD0000000338	DIFF 3.47e-2	split	tighter (KLU) → path; dense → reference; roadrunner → reference
BIOMD0000000339	DIFF 1.79e-2	split	tighter (KLU) → path; dense → reference; libRoadRunner gives up (CV_TOO_MUCH_WORK) and is left out

Both read contested. stable itself passes on both: KLU's answer does not move under jitter, as #1304 found. It is the split that marks them.

Over the corpus

Every check except locale, 1,353 models, local run on this commit.

13 contested:

• BIOMD0000000338, 339, 620, 621, 755, 1043, 1060;
• MODEL1109160000, MODEL1911100005, MODEL1912090001, MODEL2310250001, MODEL7743528808, MODEL7908934508.

7 attributed DIFFs:

model	check	verdict	arms
BIOMD0000000338	klu_vs_dense	split	tighter→path, dense→ref, roadrunner→ref
BIOMD0000000339	klu_vs_dense	split	tighter→path, dense→ref
MODEL7743528808	ajac_vs_fdjac	split	tighter→ref, klu→ref, roadrunner→path
MODEL7743528808	klu_vs_dense	split	tighter→path, klu→path, roadrunner→ref
MODEL7743528808	grid	reference	tighter→ref, klu→ref, roadrunner→neither
MODEL1912090001	klu_vs_dense	reference	tighter→ref, klu→ref
MODEL2001200002	sens_vs_plain	reference	tighter→ref, klu→ref (#1298's model; #1300 fixes it)

stable: 1,284 PASS, 10 DIFF, 54 SKIP, 4 EXCEPTION, 1 TIMEOUT. The DIFFs include:

• BIOMD0000000620 and 621, under n_points ± 1, by the whole value: at t = 10, 101 points read Anakinrasc 14,222, while 100 or 102 points read 107 as the largest species. The model has one event, and whether it fires depends on a sample landing on it: the event-at-a-sample class (A sensitivity run records the row before an event that a plain run records after it (#1133's read-ahead is off with sensitivities) #1182, A pulse schedule whose width an event writes gets no stops: X(240) = 1 for 17, at every tolerance #1244). No other check here sees it.
• Tolerance jumps: MODEL1911100005 (1.6e234 at rtol ×2), MODEL7908934508 (4.0e3), BIOMD0000000755 (4.3e2), MODEL1109160000 (3.1e2), and smaller ones.

Tests (python/tests/test_differential_shapes.py, appended)

• Attribution logic:
  • arms on the reference's side, the path's and neither give split, and contested;
  • arms that agree attribute the DIFF to one side.
• End to end: a planted cc_vs_interp DIFF (the path moved 50%) is attributed by real arms to the reference.
• Stable: an answer that jumps at rtol ×2, and one that jumps at n_points + 1, read DIFF and contested. All 15 shapes pass stable.
• BIOMD0000000338 and 339: KLU and the dense solver take different branches past a step-size collapse (Fg 311 on KLU, 0.73 on dense and libRoadRunner) #1304: BIOMD0000000338 reads split with tighter → path and dense → reference. It skips if the corpus is absent, as the repo's other corpus tests do.

Mutations, each caught:

• attribution always "reference";
• check_all not attributing;
• stable without its tolerance jitter;
• stable without its grid jitter;
• contested ignoring a split;
• no other-solver arm.

Suites: parity_checks/tests and the shapes give 1,105 passed, 28 skipped, 1 xfailed. ruff 0.14.0 is clean. The branch merges cleanly with lanl main.

Cost

• stable: 5 more runs per model.
• Attribution: 2–3 runs, on DIFF models only (a handful a night).

Baseline

parity_checks/nightly/baselines/differential.json will be re-recorded from this PR's own nightly run, in a follow-up commit here. It adds the stable rows; the contested list rides in the report's meta.

Changes (this PR's commit, less the tests)

diff --git a/parity_checks/differential/_differential.py b/parity_checks/differential/_differential.py
index 735a1e8f..3b28e7be 100644
--- a/parity_checks/differential/_differential.py
+++ b/parity_checks/differential/_differential.py
@@ -104,6 +104,8 @@ CHECKS = (
     "grid",
     "sparsity_covers_jacobian",
     "backend_as_asked",
+    # #1309: the base's answer under a jittered grid and tolerance.
+    "stable",
 )
 
 
@@ -113,6 +115,10 @@ class Verdict:
     value: float | None = None
     note: str = ""
     extra: dict = field(default_factory=dict)
+    # The two trajectories a comparison judged, kept for attribution (#1309)
+    # and never written to a report.
+    ref_x: np.ndarray | None = field(default=None, repr=False, compare=False)
+    path_x: np.ndarray | None = field(default=None, repr=False, compare=False)
 
 
 @dataclass
@@ -212,6 +218,12 @@ class Base:
         result, sim = run(job)
         self.job = job
         self.x = np.asarray(result.species, dtype=float)
+        # The columns of x, by name: a Result's species, which leave out a
+        # parameter an event promotes into the state (Model.species_names has it).
+        self.species_names = list(result.species_names)
+        # The linear solver the base ran on (0 dense, 1 KLU, 2 BLAS dense), so
+        # attribution can run it on the other one (#1309).
+        self.linear_solver = int(result.solver_stats.get("linear_solver", 0))
         self.backend = sim.codegen_backend
         # A tenfold tighter run that does not complete leaves the base no
         # reference, not the model unrun: the checks that compare no answers
@@ -257,7 +269,10 @@ class Base:
                 f"{label}nothing could count: the trajectory's peak {peak:.1e} is under the "
                 f"absolute floor {floor:.1e} ({PATH_RTOL_FACTOR:g}·atol); {note}",
             )
-        return Verdict("pass" if g <= self.limit else "diff", g, f"{label}{note}")
+        status = "pass" if g <= self.limit else "diff"
+        return Verdict(
+            status, g, f"{label}{note}", ref_x=self.x, path_x=np.asarray(x, dtype=float)
+        )
 
 
 def _first_line(e: BaseException) -> str:
@@ -609,6 +624,147 @@ def check_sparsity(base: Base) -> Verdict:
     return _guard(go)
 
 
+# Attribution and branch stability (issue #1309). A DIFF says two paths
+# disagree, not which is wrong, and a small own drift does not make an answer
+# uncontested: on #1304's BIOMD0000000338 the KLU base moved 1e-7 at a tenfold
+# tighter tolerance while dense and libRoadRunner sat 3.5% away, and dense once
+# landed on KLU's branch.
+
+ATTRIBUTION_TIGHTEN = 100.0
+JITTER_RTOL = (0.5, 2.0)
+
+
+def _arm_runs(base: Base) -> dict[str, np.ndarray]:
+    """The third arms, each a trajectory in the base's columns: the base's own
+    route at ``ATTRIBUTION_TIGHTEN`` x tighter tolerances, the base on the
+    other linear solver, and libRoadRunner for SBML where it is installed. An
+    arm that does not run, or does not apply, is left out."""
+    import bngsim
+
+    job, out = base.job, {}
+
+    def take(name, fn):
+        # An arm that cannot run (libRoadRunner gives up, KLU is declined) says
+        # nothing, and is left out of the attribution.
+        with contextlib.suppress(Exception):
+            out[name] = fn()
+
+    f = ATTRIBUTION_TIGHTEN
+    tight = Job(**{**job.__dict__, "rtol": job.rtol / f, "atol": job.atol / f})
+    take("tighter", lambda: _species(run(tight)[0], base))
+    if base.linear_solver == 1:
+        take("dense", lambda: _species(run(job, force_dense_linear_solver=True)[0], base))
+    elif bngsim._bngsim_core.HAS_KLU:
+        take("klu", lambda: _species(run(job, force_sparse_linear_solver=True)[0], base))
+    if Path(job.path).suffix != ".net":
+        take("roadrunner", lambda: _roadrunner(base))
+    return out
+
+
+def _roadrunner(base: Base) -> np.ndarray:
+    import roadrunner
+
+    job = base.job
+    r = roadrunner.RoadRunner(str(job.path))
+    r.integrator.relative_tolerance = job.rtol
+    r.integrator.absolute_tolerance = job.atol
+    # A species' concentration, as bngsim reports it, is [id] in libRoadRunner
+    # (the bare id is its amount); a rate rule's parameter is its id. As
+    # rr_parity selects them.
+    species = set(r.model.getFloatingSpeciesIds()) | set(r.model.getBoundarySpeciesIds())
+    sel = [f"[{n}]" if n in species else n for n in base.species_names]
+    out = r.simulate(job.t_start, job.t_end, job.n_points, sel)
+    return np.asarray(out, dtype=float)
+
+
+def _species(result, base: Base) -> np.ndarray:
+    """A Result's species in the base's column order."""
+    names = list(result.species_names)
+    x = np.asarray(result.species, dtype=float)
+    return x[:, [names.index(n) for n in base.species_names]]
+
+
+def contested(rows: dict[str, dict]) -> bool:
+    """Whether a model's rows mark it contested: its ``stable`` row reads DIFF
+    (the answer jumps under a jittered grid or tolerance), or a DIFF's third
+    arms split between its two sides (#1304's BIOMD0000000338: KLU at a
+    hundredfold tighter tolerance sides with KLU, the dense solver and
+    libRoadRunner with dense). Its DIFFs are real but unattributed, its PASSes
+    luck. ``rows`` maps a check to ``{"status", "extra"}``, as a report holds
+    them."""
+    for check, row in rows.items():
+        extra = row.get("extra") or {}
+        if check == "stable" and row.get("status") == "diff":
+            return True
+        if (extra.get("attribution") or {}).get("verdict") == "split":
+            return True
+    return False
+
+
+def attribute(base: Base, v: Verdict, arms: dict[str, np.ndarray] | None = None) -> dict:
+    """Which side of a DIFF each third arm agrees with: the comparison's
+    ``reference`` (the base, or the path it was judged against), its
+    ``path``, ``both`` (the two are within the bound of the arm), or
+    ``neither``. Its ``verdict`` is the one side every arm that took one
+    took, ``split`` where they differ (the model is contested), or
+    ``unattributed``."""
+    arms = _arm_runs(base) if arms is None else arms
+    atol = PATH_RTOL_FACTOR * base.job.atol
+    sides = {}
+    for name, x in arms.items():
+        to_ref, to_path = gap(x, v.ref_x, atol), gap(x, v.path_x, atol)
+        near_ref, near_path = to_ref <= base.limit, to_path <= base.limit
+        side = {(True, False): "reference", (False, True): "path", (True, True): "both"}.get(
+            (near_ref, near_path), "neither"
+        )
+        sides[name] = {"side": side, "to_reference": to_ref, "to_path": to_path}
+    taken = {a["side"] for a in sides.values() if a["side"] in ("reference", "path")}
+    verdict = taken.pop() if len(taken) == 1 else ("split" if taken else "unattributed")
+    return {"verdict": verdict, "arms": sides}
+
+
+def check_stable(base: Base) -> Verdict:
+    """The base's answer under a jittered grid (one point fewer, one more,
+    twice as many) and tolerance (rtol x0.5, x2): a jump past the bound marks
+    the model contested, its DIFFs real but unattributed and its PASSes luck.
+    A grid of another length is read at its last point, which every grid
+    shares; twice as many at every other point, which are the base's."""
+
+    def go():
+        job, jumps = base.job, {}
+        atol = PATH_RTOL_FACTOR * job.atol
+        n = job.n_points
+        for label, m in (
+            ("n_points-1", n - 1),
+            ("n_points+1", n + 1),
+            ("2x the points", 2 * n - 1),
+        ):
+            if m < 2:
+                continue
+            x = _species(run(Job(**{**job.__dict__, "n_points": m}))[0], base)
+            g = gap(x[::2], base.x, atol) if m == 2 * n - 1 else gap(x[-1:], base.x[-1:], atol)
+            if g > base.limit:
+                jumps[label] = g
+        for f in JITTER_RTOL:
+            x = _species(
+                run(Job(**{**job.__dict__, "rtol": job.rtol * f, "atol": job.atol * f}))[0], base
+            )
+            g = gap(x, base.x, atol)
+            if g > max(base.limit, f * base.limit):
+                jumps[f"rtol x{f:g}"] = g
+        if jumps:
+            shown = ", ".join(f"{k} {g:.1e}" for k, g in jumps.items())
+            return Verdict(
+                "diff",
+                max(jumps.values()),
+                f"contested: the answer jumps past its bound {base.limit:.1e} under {shown}",
+                {"contested": True, "jumps": jumps},
+            )
+        return Verdict("pass", 0.0, "the answer holds under a jittered grid and tolerance")
+
+    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."""
@@ -682,10 +838,22 @@ def check_all(job: Job, locale_name: str | None, checks=CHECKS) -> dict[str, Ver
         "grid": check_grid,
         "sparsity_covers_jacobian": check_sparsity,
         "backend_as_asked": check_backend,
+        "stable": check_stable,
     }
     for name, fn in one.items():
         if name in checks:
             out[name] = fn(base)
     if "locale" in checks:
         out["locale"] = check_locale(base, locale_name)
+    # Every DIFF with two trajectories says which side the third arms take
+    # (#1309), the arms run once per model however many DIFFs it has.
+    diffs = [v for v in out.values() if v.status == "diff" and v.path_x is not None]
+    if diffs:
+        try:
+            arms = _arm_runs(base)
+            for v in diffs:
+                v.extra["attribution"] = attribute(base, v, arms)
+        except Exception as e:  # noqa: BLE001 - attribution is advice, never a verdict
+            for v in diffs:
+                v.extra["attribution"] = {"verdict": "unattributed", "error": str(e)[:300]}
     return {c: out[c] for c in checks if c in out}
diff --git a/parity_checks/differential/differential_run.py b/parity_checks/differential/differential_run.py
index 0ac3d2e4..d5bf54b3 100644
--- a/parity_checks/differential/differential_run.py
+++ b/parity_checks/differential/differential_run.py
@@ -274,6 +274,13 @@ def main(argv=None) -> int:
         specs, _worker, workers=args.workers, timeout_of=lambda s: s["timeout"], on_done=on_done
     )
     stamp = versions.stamp(REFERENCE_ENGINE)
+    # #1309: a model whose answer is contested, by its stable row or by a
+    # DIFF whose third arms split, listed by name.
+    contested = sorted(
+        spec["job"]["model_id"]
+        for spec, r in zip(specs, raw, strict=False)
+        if D.contested(r.get("rows") or {})
+    )
     results = [x for spec, r in zip(specs, raw, strict=False) for x in to_results(spec, r, stamp)]
     meta = {
         "suite": "differential",
@@ -282,6 +289,7 @@ def main(argv=None) -> int:
         "git_rev": git_rev(str(HERE)),
         "versions": stamp,
         "tally": tally(r.outcome for r in results),
+        "contested": contested,
         "n_jobs": len(results),
         "n_models": len(specs),
         "elapsed_sec": round(time.perf_counter() - t0, 2),
@@ -299,6 +307,8 @@ def main(argv=None) -> int:
     out.parent.mkdir(parents=True, exist_ok=True)
     write_report(out, results, meta=meta)
     print(f"report: {out}  {meta['tally']}")
+    if contested:
+        print(f"contested ({len(contested)}): {', '.join(contested)}")
     return 0
 
 

No changelog fragment: parity_checks/ and tests only.

The nightly's cap

Stacked on #1321, which fixes the cause on main for every PR here. Until it merges, this PR also shows its commits.

The first nightly on this PR failed Path-agreement check with 33 regressions, every one a TIMEOUT ("cap of 900 s"). The pool kills a model's worker at its cap, and every row of that model goes with it:

• BIOMD0000000810: its checks take 8 s, but one stable rerun ground on past the cap (CVODE error-test failures, at up to 100,000 steps per output point);
• MODEL1006230049, 0053 and 0077: these already took 716–751 s, and stable's five reruns tipped them over.

#1315, #1318 and #1320 failed the same way. #1321 bounds every run in the harness by the model's cap, and keeps a killed worker's finished rows. Here (merged into this branch, replacing 6357188's narrower fix) the additions run after every check main runs:

1. the attribution's arms, each run bounded; an arm that does not finish is left out;
2. then stable.

Each starts only with time left, and a stable the cap stops is TIMEOUT, its own row only.

BIOMD0000000810 end to end through the runner, at a 120 s cap: every row returns in 109 s. Its three baseline DIFFs are attributed (ajac_vs_fdjac → reference, klu_vs_dense → path, grid → reference), and stable is TIMEOUT, "not finished within the model's cap".

Tests:

• with the cap spent after main's checks, stable is TIMEOUT and every other row matches a run without a cap;
• a planted DIFF's attribution in a spent cap asks no arm.

Mutations:

• run unbounded: caught;
• stable starting without its own check of the time left: an equivalent mutant, since its first run refuses to start past the deadline anyway.

Suites: the shapes and the skip audit give 233 passed, 1 xfailed.