From 0bf4bdb9d0fe87f5c3a6f2f94a98f821f3d8b181 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Thu, 22 Oct 2026 07:58:20 +0200
Subject: [PATCH] still battery: score R-free against the deposited data
 through the same dimple protocol

The published R-free comes from the depositors' own refinement (mostly phenix.refine), while the
battery refines with dimple's 8 REFMAC cycles, so the two are not like-for-like: on LCLS 222 (8TDQ)
the deposited intensities reach 0.1918 under dimple against a published 0.2012.

Every set now carries "R-free deposition (dimple)": the deposition's -sf.cif (intensity_meas, else
the mean of the anomalous pair, else F_meas_au squared) refined with the same refine() as our merge -
same deposited model, the deposition's free set, cut at its resolution. It is computed once per entry
and resolution and kept in the PDB cache (<pdb>-dimple-<d_cut>/reference.json); `run` and `report`
fill in a missing one, `reference` computes or lists them up front. The report gains the difference
(ours minus the deposition's) and a quality verdict (win below -0.01, draw within 0.01, loss
otherwise); the published R-free stays as an informational column.

Reproduces the hand-computed references: 8TDQ 0.1918, 9V3D 0.2220, 31JT 0.1945.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
---
 tools/battery/README.md        |  12 +++-
 tools/battery/still_battery.py | 108 ++++++++++++++++++++++++++++++---
 2 files changed, 108 insertions(+), 12 deletions(-)

diff --git a/tools/battery/README.md b/tools/battery/README.md
index 4fe3d9b79..37ce2a877 100644
--- a/tools/battery/README.md
+++ b/tools/battery/README.md
@@ -949,8 +949,10 @@ once and its row says what the merge is worth and what it cost:
 | facility, beamline, detector, dataset, PDB | from the manifest |
 | images | frames in the input files (rugnux's `IMAGES_PROCESSED`); images (deposition) the deposition's own count, where it is known |
 | indexed %, d_min, R_split | rugnux's report (`INDEXING_RATE`, `INCLUDE_RESOLUTION_RANGE`, `R_SPLIT`: whole frames go to one half, as in CrystFEL) |
-| R_free rugnux | `dimple` of the deposited model against our merge, reindexed into the deposition's space group, with the deposition's own free set, at the deposition's resolution |
-| R_free deposited | the entry's published R-free (dark and ground-state entries only) |
+| R_free rugnux (dimple) | `dimple` of the deposited model against our merge, reindexed into the deposition's space group, with the deposition's own free set, at the deposition's resolution |
+| R_free deposition (dimple) | the same protocol on the deposited data (the entry's `-sf.cif`: `intensity_meas`, else the mean of the anomalous pair, else `F_meas_au` squared): the reference |
+| difference, quality | ours minus the deposition's; win where ours is lower by more than 0.01, draw within 0.01, loss otherwise |
+| R_free published | the entry's own R-free (dark and ground-state entries only), from its own refinement program: for information, not like-for-like |
 | time, ms/image | wall clock of the rugnux command, which ends with the model validation |
 | I/O floor | the input's bytes over `--io-gbps` (default 5 GB/s): the time reading alone would take from an SSD |
 | page cache | share of the input already in memory when the run started (`fincore`): cold, warm or partial |
@@ -966,9 +968,15 @@ still_battery.py run --rugnux build/rugnux/rugnux [--only pal90,sacla236 | --tie
 still_battery.py compare RUN_A RUN_B
 still_battery.py report RUN [--out DIR]
 still_battery.py list
+still_battery.py reference [--only pal90,sacla236]
 still_battery.py keep RUN --as workstation|server
 ```
 
+The deposited-data reference is computed once per entry and resolution and kept in the PDB cache
+(`<pdb_cache>/<pdb>-dimple-<d_cut>/reference.json`, with the dimple run beside it); `run` and
+`report` compute a missing one on the way (so `report` of an older run gains the column), and
+`reference` computes or lists them all up front.
+
 A run writes `<runs_root>/<YYYYMMDD-HHMM>_<rugnux git>_still-<label>/` with `manifest.json`,
 `results.json`, `report.md` and `work/<set>/` (the rugnux outputs and the dimple run), made
 read-only when done; `keep` copies the three files into `results/<role>/<run>/`, like a rotation run.
diff --git a/tools/battery/still_battery.py b/tools/battery/still_battery.py
index 57bb93a0b..625df9dd9 100644
--- a/tools/battery/still_battery.py
+++ b/tools/battery/still_battery.py
@@ -5,8 +5,9 @@ merge, report, compare.
 The rotation battery (battery.py) scores a lattice and a space group. For stills both are given
 (-C/-S, the way a serial user runs), so this battery reports what the merge is worth and what it
 costs instead: indexing rate, R_split, the R-free of the deposited model refined against our merge
-(dimple, with the deposition's own free set) beside the deposition's published R-free, and the wall
-time of the rugnux command with the input's I/O floor next to it.
+(dimple, with the deposition's own free set) beside the R-free of the same protocol on the deposited
+data (the reference the verdict is taken against; the published R-free is shown for information), and
+the wall time of the rugnux command with the input's I/O floor next to it.
 
 Usage (README.md in this directory, section "The still battery", is the guide):
   still_battery.py run --rugnux BIN [--label L] [--only a,b | --tier smoke] [--threads N]
@@ -14,6 +15,7 @@ Usage (README.md in this directory, section "The still battery", is the guide):
   still_battery.py compare RUN_A RUN_B
   still_battery.py report RUN [--out DIR]          re-render a run's report.md
   still_battery.py list                            still runs under the runs root
+  still_battery.py reference [--only a,b]          compute (or show) the deposited-data dimple references
   still_battery.py keep RUN --as workstation|server  copy report.md, results.json, manifest.json
                                                    into results/<role>/<run>/ (committed); refuses a
                                                    private run
@@ -48,11 +50,13 @@ import numpy as np  # noqa: E402
 MANIFEST = os.path.join(HERE, "still_manifest.json")
 KEPT_FILES = ("report.md", "results.json", "manifest.json")
 IO_GBPS = 5.0          # sequential read speed the I/O floor is quoted at (a PCIe 4 NVMe SSD)
+DRAW = 0.01            # R-free difference to the deposition's (same dimple protocol) still called a draw
 
 COLUMNS = [("facility", "facility"), ("beamline", "beamline"), ("detector", "detector"),
            ("dataset", "set"), ("PDB", "pdb"), ("images", "images"), ("images (deposition)", "images_deposited"),
            ("indexed %", "indexed_pct"), ("d_min", "d_min"), ("R_split", "r_split"),
-           ("R_free rugnux", "r_free"), ("R_free deposited", "r_free_deposited"),
+           ("R_free rugnux (dimple)", "r_free"), ("R_free deposition (dimple)", "r_free_dimple_deposited"),
+           ("difference", "r_free_diff"), ("quality", "quality"), ("R_free published", "r_free_deposited"),
            ("time s", "wall_s"), ("ms/image", "ms_per_image"), ("I/O floor s", "io_floor_s"),
            ("page cache", "cache")]
 
@@ -233,6 +237,62 @@ def refine(wd, e, cache):
     return float(m[-1][0]), float(m[-1][1]), f"free set {free}"
 
 
+def deposited_mtz(sf_cif, out_mtz):
+    """IMEAN/SIGIMEAN of the deposition from its -sf.cif: intensity_meas where deposited, else the mean
+    of the anomalous pair pdbx_I_plus/minus, else F_meas_au squared."""
+    rb = gemmi.as_refln_blocks(gemmi.cif.read(sf_cif))[0]
+    cols = rb.column_labels()
+    if "intensity_meas" in cols:
+        i, s = rb.make_float_array("intensity_meas"), rb.make_float_array("intensity_sigma")
+    elif "pdbx_I_plus" in cols:
+        ip, sp, im, sm = (rb.make_float_array(c) for c in ("pdbx_I_plus", "pdbx_I_plus_sigma",
+                                                           "pdbx_I_minus", "pdbx_I_minus_sigma"))
+        n = np.isfinite(ip).astype(int) + np.isfinite(im)
+        with np.errstate(invalid="ignore"):
+            i = (np.nan_to_num(ip) + np.nan_to_num(im)) / n
+            s = np.sqrt(np.nan_to_num(sp) ** 2 + np.nan_to_num(sm) ** 2) / n
+    else:
+        f, sf = rb.make_float_array("F_meas_au"), rb.make_float_array("F_meas_sigma_au")
+        i, s = f ** 2, 2 * f * sf
+    ok = np.isfinite(i)
+    mtz = gemmi.Mtz(with_base=True)
+    mtz.spacegroup = rb.spacegroup
+    mtz.cell = rb.cell
+    mtz.add_dataset("deposited")
+    mtz.add_column("IMEAN", "J")
+    mtz.add_column("SIGIMEAN", "Q")
+    mtz.set_data(np.column_stack([rb.make_miller_array()[ok], i[ok], s[ok]]).astype(np.float32))
+    mtz.write_to_file(out_mtz)
+
+
+def deposited_reference(e, cache):
+    """R-free of the deposited data through the same dimple protocol as our merge (refine()): what the
+    deposition reaches under 8 REFMAC cycles, so the like-for-like reference rather than the published
+    R-free. Computed once per entry and resolution and kept in the PDB cache. {r_work, r_free, note}."""
+    pdb = e["pdb"].lower()
+    wd = os.path.join(cache, f"{pdb}-dimple-{e['d_cut']}")
+    done = os.path.join(wd, "reference.json")
+    if os.path.exists(done):
+        return json.load(open(done))
+    sf = model_check.fetch(model_check.RCSB + pdb + "-sf.cif", os.path.join(cache, pdb + "-sf.cif"))
+    shutil.rmtree(wd, ignore_errors=True)
+    os.makedirs(wd)
+    deposited_mtz(sf, os.path.join(wd, "p.mtz"))
+    r_work, r_free, note = refine(wd, e, cache)
+    ref = {"r_work": r_work, "r_free": r_free, "note": note}
+    if r_free is not None:
+        battery.write_json(done, ref)
+    return ref
+
+
+def quality(r_free, reference):
+    """win / draw / loss of our dimple R-free against the deposition's under the same protocol."""
+    if r_free is None or reference is None:
+        return None
+    d = r_free - reference
+    return "win" if d < -DRAW else "draw" if d <= DRAW else "loss"
+
+
 def run_one(run_dir, binary, e, st, opts):
     wd = os.path.join(run_dir, "work", e["id"])
     os.makedirs(wd)
@@ -270,6 +330,7 @@ def run_one(run_dir, binary, e, st, opts):
         return dict(row, verdict="crash" if rc else "fail", reason=battery.failure_note(os.path.join(wd, "run.log")))
     if opts["refine"] and e.get("pdb"):
         row["r_work"], row["r_free"], row["refine_note"] = refine(wd, e, st["pdb_cache"])
+        row["r_free_dimple_deposited"] = deposited_reference(e, st["pdb_cache"])["r_free"]
     return dict(row, verdict="done", reason="")
 
 
@@ -280,8 +341,10 @@ def fmt(v, key):
         return "-"
     if key == "r_split":
         return f"{100 * v:.1f}%"
-    if key in ("r_free", "r_free_deposited"):
+    if key in ("r_free", "r_free_deposited", "r_free_dimple_deposited"):
         return f"{v:.3f}"
+    if key == "r_free_diff":
+        return f"{v:+.3f}"
     if key == "images_deposited" and isinstance(v, (int, float)):
         return f"{int(v):,}"
     if key == "images":
@@ -289,9 +352,18 @@ def fmt(v, key):
     return str(v)
 
 
-def render(run_dir):
+def render(run_dir, cache):
+    """report.md of a run. A run from before the deposited-data reference existed gets it from the
+    PDB cache (computed there if missing)."""
     man = json.load(open(os.path.join(run_dir, "manifest.json")))
     res = json.load(open(os.path.join(run_dir, "results.json")))
+    sets = {e["id"]: e for e in man["sets"]}
+    for r in res:
+        if r.get("r_free") is not None and r.get("r_free_dimple_deposited") is None:
+            r["r_free_dimple_deposited"] = deposited_reference(sets[r["set"]], cache)["r_free"]
+        if r.get("r_free") is not None and r.get("r_free_dimple_deposited") is not None:
+            r["r_free_diff"] = r["r_free"] - r["r_free_dimple_deposited"]
+        r["quality"] = quality(r.get("r_free"), r.get("r_free_dimple_deposited"))
     b = man["binary"]
     hw = man.get("hardware") or {}
     lines = [f"# Rugnux still battery - {os.path.basename(os.path.normpath(run_dir))}", "",
@@ -301,8 +373,11 @@ def render(run_dir):
              f"{', '.join(g['name'] for g in hw.get('gpus', [])) or 'none'}",
              f"- start {man['start']}, end {man.get('end') or 'unfinished'}; subset: {man.get('subset') or 'all sets'}"
              + ("; ABORTED" if man.get("aborted") else ""),
-             f"- R_free rugnux: dimple of the deposited model against our merge, with the deposition's free set, "
-             f"at its resolution; R_free deposited: the entry's published value. Time: wall clock of the rugnux "
+             f"- R_free rugnux (dimple): dimple of the deposited model against our merge, with the deposition's "
+             f"free set, at its resolution; R_free deposition (dimple): the same protocol on the deposited data; "
+             f"difference: ours minus the deposition's; quality: win where ours is lower by more than {DRAW}, "
+             f"draw within {DRAW}, loss otherwise; R_free published: the entry's own value (phenix etc., "
+             f"not like-for-like, for information). Time: wall clock of the rugnux "
              f"command (to model validation); I/O floor: input bytes / {man['options']['io_gbps']} GB/s; page "
              f"cache: share of the input already in memory when the run started.", "",
              "| " + " | ".join(c for c, _ in COLUMNS) + " |", "|" + "---|" * len(COLUMNS)]
@@ -356,7 +431,7 @@ def cmd_run(a, st):
         sys.exit(f"aborted; {run_dir} is marked as such")
     man.update(end=datetime.datetime.now().isoformat(timespec="seconds"), complete=True)
     battery.write_json(os.path.join(run_dir, "manifest.json"), man)
-    print(f"report {render(run_dir)}", flush=True)
+    print(f"report {render(run_dir, st['pdb_cache'])}", flush=True)
     battery.make_read_only(run_dir)
     if any(r["verdict"] in ("crash", "fail") for r in results):
         sys.exit("some sets failed: " + ", ".join(r["set"] for r in results if r["verdict"] in ("crash", "fail")))
@@ -393,7 +468,16 @@ def cmd_report(a, st):
         for f in ("manifest.json", "results.json"):
             shutil.copyfile(os.path.join(src, f), os.path.join(tmp, f))
         src = tmp
-    print(render(src))
+    print(render(src, st["pdb_cache"]))
+
+
+def cmd_reference(a, st):
+    print("| set | PDB | d_cut | R_work | R_free deposition (dimple) | R_free published |")
+    print("|---|---|---|---|---|---|")
+    for e in select(st, a.only, None, a.private):
+        ref = deposited_reference(e, st["pdb_cache"])
+        print(f"| {e['id']} | {e['pdb']} | {e['d_cut']} | {ref['r_work']} | {ref['r_free']} | "
+              f"{e.get('r_free_deposited')} |", flush=True)
 
 
 def cmd_keep(a, st):
@@ -431,12 +515,16 @@ def main():
     p.add_argument("run")
     p.add_argument("--out")
     sub.add_parser("list")
+    f = sub.add_parser("reference")
+    f.add_argument("--only")
+    f.add_argument("--private", action="store_true")
     k = sub.add_parser("keep")
     k.add_argument("run")
     k.add_argument("--as", dest="role", required=True, choices=["workstation", "server"])
     a = ap.parse_args()
     st = load_site(a.site)
-    {"run": cmd_run, "compare": cmd_compare, "list": cmd_list, "report": cmd_report, "keep": cmd_keep}[a.cmd](a, st)
+    {"run": cmd_run, "compare": cmd_compare, "list": cmd_list, "report": cmd_report, "reference": cmd_reference,
+     "keep": cmd_keep}[a.cmd](a, st)
 
 
 if __name__ == "__main__":