From 07b0f7e31569156bc016b59057e72f1ef9bf3418 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 25 Sep 2026 15:53:47 +0200 Subject: [PATCH] battery: paper_plots.py draws the open-arm R-free and resolution figures from a run The figures, their value pairs, the excluded sets with reasons, the gnuplot script and the runs' provenance are written together, so a published figure can be audited and redrawn from the battery run it came from. Facility and beamline counts come from the entries' _diffrn_source. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- tools/battery/README.md | 13 +++ tools/battery/paper_plots.py | 181 +++++++++++++++++++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 tools/battery/paper_plots.py diff --git a/tools/battery/README.md b/tools/battery/README.md index d236c01e6..f793fd403 100644 --- a/tools/battery/README.md +++ b/tools/battery/README.md @@ -442,3 +442,16 @@ whole range; the manifest keeps the rule used as `dmin_rule`. `{"old id": {"id": "...", "input": "..."}}`. Each old id is recorded under `aliases` in the manifest, so runs made before the rename still compare set by set. Without `--write`, the command only shows the result and flags inputs that do not exist. + +## Figures for publication + +`paper_plots.py RUN [RUN-REDO ...] --out DIR` draws the two open-arm figures from a run's +`results.json`: R-free of the deposited model after the REFMAC check (`--model-check`) on the +depositor's data against Rugnux's data (same model, protocol, resolution and the depositor's free +set), and the deposited high-resolution limit against Rugnux's cut. Later runs fill in sets an +earlier one did not run (a `-redo`). Beside the PDF/PNG figures it writes the value pairs +(`rfree.dat`, `dmin.dat`, with each set's beamline), every excluded set with its reason +(`excluded.txt`), the gnuplot script (`figures.gp`) and `summary.txt` with the runs' binary and +runner provenance and the structure, facility and beamline counts. Facility and beamline are read +from each entry's cached mmCIF (`_diffrn_source`). Needs gnuplot and gemmi. + diff --git a/tools/battery/paper_plots.py b/tools/battery/paper_plots.py new file mode 100644 index 000000000..40eb58937 --- /dev/null +++ b/tools/battery/paper_plots.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python3 +"""paper_plots.py -- the two open-arm figures: R-free and high-resolution limit, deposition vs Rugnux. + + tools/battery/paper_plots.py RUN [RUN ...] --out DIR + +RUN is a battery run directory (or its name under the site's runs_root); later runs fill in sets +an earlier one lacks or that failed to run there (a `-redo` of CUDA casualties). Only the open arm +is read. Everything the figures are drawn from is written beside them, so a figure can be audited +and redrawn: rfree.dat, dmin.dat (value pairs with the set and its beamline), excluded.txt (every +set left out, with the reason), summary.txt (the runs, their binary and the counts), figures.gp +(the gnuplot script) and the rendered figures (PDF for print, PNG for screens). gnuplot must be +on PATH; gemmi reads the facility and beamline from each entry's cached mmCIF (_diffrn_source). + +Left figure: R-free of the deposited model after the battery's REFMAC check (model_check.py: the +same model, protocol, resolution and the depositor's free set) on the depositor's data (x) and on +Rugnux's data (y). Right figure: the deposited high-resolution limit (x) against Rugnux's own cut +(y), log axes so the low-resolution entries show. Small molecules (no deposited model) and sets +whose verdict is `fail` (a wrong lattice: nothing to compare) are left out. +""" +import argparse +import json +import math +import os +import re +import shutil +import subprocess +import sys + +import gemmi + +RUNS_ROOT = "/data/battery/runs" +RFREE_BAND = 0.05 # |R-free Rugnux - R-free depositor| counted as "within" + +GNUPLOT = r""" +set encoding utf8 +set size ratio 1 +set grid lc rgb "#cccccc" lw 1 +set border lw 2 +set tics scale 1.5 +set key off +set style line 1 lc rgb "black" lw 3 dt 2 +set style line 2 pt 7 ps 1.1 lc rgb "#1f77b4" +set style line 3 pt 7 ps 1.1 lc rgb "#d62728" +label_box = "Structure number: %(n)d\nFacilities number: %(nf)d\nBeamline number: %(nb)d" + +# --- R-free +set output "%(out)s/rfree.%(ext)s" +set xlabel "R_{free}, depositor data (REFMAC)" +set ylabel "R_{free}, Rugnux data (REFMAC)" +set xrange [%(rlo).2f:%(rhi).2f]; set yrange [%(rlo).2f:%(rhi).2f] +set xtics 0.05; set ytics 0.05 +set format x "%%.2f"; set format y "%%.2f" +set label 1 sprintf(label_box) at graph 0.04, graph 0.96 left front +set label 2 "Within ±%(band).2f: %(within)d of %(n)d" at graph 0.96, graph 0.06 right front +plot x ls 1, "%(out)s/rfree.dat" using 1:2 ls 2 +unset label 1; unset label 2 + +# --- high-resolution limit +set output "%(out)s/dmin.%(ext)s" +set logscale xy +set xlabel "d_{min}, deposition (Å)" +set ylabel "d_{min}, Rugnux (Å)" +set xrange [%(dlo).2f:%(dhi).2f]; set yrange [%(dlo).2f:%(dhi).2f] +set xtics (0.5, 0.7, 1, 1.5, 2, 3, 4, 6, 8) +set ytics (0.5, 0.7, 1, 1.5, 2, 3, 4, 6, 8) +set format x "%%g"; set format y "%%g" +set label 1 sprintf("Structure number: %(dn)d\nFacilities number: %(dnf)d\nBeamline number: %(dnb)d") at graph 0.04, graph 0.96 left front +set label 3 "Rugnux reaches further (%(deeper)d)" at graph 0.96, graph 0.06 right front +set label 4 "Deposition reaches further (%(shallower)d)" at graph 0.04, graph 0.70 left front +plot x ls 1, "%(out)s/dmin.dat" using 1:2 ls 3 +""" + + +def load_open(runs): + """Rows of the open arm, a later run's row replacing one that did not run.""" + rows, manifests = {}, [] + for run in runs: + path = run if os.path.isdir(run) else os.path.join(RUNS_ROOT, run) + manifests.append((path, json.load(open(os.path.join(path, "manifest.json"))))) + for r in json.load(open(os.path.join(path, "results.json"))): + if r["arm"] != "open": + continue + if r["set"] not in rows or rows[r["set"]].get("exit_code"): + rows[r["set"]] = r + return rows, manifests + + +def source(model_path): + """(facility, beamline) from the entry's _diffrn_source, facility without its operator suffix + ("PETRA III, EMBL c/o DESY" -> "PETRA III"); None where the entry does not record one.""" + block = gemmi.cif.read(model_path).sole_block() + + def value(tag): + v = block.find_value(tag) + return None if v in (None, "?", ".") else gemmi.cif.as_string(v).strip() + + site = value("_diffrn_source.pdbx_synchrotron_site") + if not site: + return None + beamline = value("_diffrn_source.pdbx_synchrotron_beamline") or "?" + facility = site.split(",")[0].strip().upper() + return facility, facility + " " + re.sub(r"\s+", "", beamline.upper()) + + +def counts(sets, where): + known = [where[s] for s in sets if where[s]] + return len(sets), len({f for f, _ in known}), len({b for _, b in known}) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) + ap.add_argument("runs", nargs="+", help="battery run directories or names (later ones fill in)") + ap.add_argument("--out", required=True, help="output directory") + a = ap.parse_args() + if not shutil.which("gnuplot"): + sys.exit("gnuplot not found on PATH") + os.makedirs(a.out, exist_ok=True) + out = os.path.abspath(a.out) + + rows, manifests = load_open(a.runs) + excluded, rfree, dmin, where = [], [], [], {} + for s, r in sorted(rows.items()): + if not r.get("model"): + excluded.append((s, "both", "no deposited model (small molecule or unpublished)")) + continue + if r.get("verdict") == "fail": + excluded.append((s, "both", "verdict fail: " + str(r.get("reason")))) + continue + where[s] = source(r["model"]) + dep, rug = r.get("refmac_rfree_depdata"), r.get("refmac_rfree_depflags") + if dep and rug: + rfree.append((dep, rug, s)) + else: + excluded.append((s, "rfree", "no REFMAC pair: " + str(r.get("refmac_reason") + or r.get("dep_reason") or "depositor data not scored"))) + if r.get("d_min_ref") and r.get("d_min"): + dmin.append((r["d_min_ref"], r["d_min"], s)) + else: + excluded.append((s, "dmin", "no deposited or Rugnux d_min")) + + bl = lambda s: (where[s][1] if where[s] else "not-recorded").replace(" ", "_") + with open(os.path.join(out, "rfree.dat"), "w") as fh: + fh.write("# rfree_depositor_data rfree_rugnux_data set beamline\n") + fh.writelines(f"{x:.4f} {y:.4f} {s} {bl(s)}\n" for x, y, s in rfree) + with open(os.path.join(out, "dmin.dat"), "w") as fh: + fh.write("# dmin_deposition dmin_rugnux set beamline\n") + fh.writelines(f"{x:.3f} {y:.3f} {s} {bl(s)}\n" for x, y, s in dmin) + with open(os.path.join(out, "excluded.txt"), "w") as fh: + fh.writelines(f"{s}\t{which}\t{why}\n" for s, which, why in excluded) + + n, nf, nb = counts([s for _, _, s in rfree], where) + dn, dnf, dnb = counts([s for _, _, s in dmin], where) + within = sum(abs(y - x) <= RFREE_BAND for x, y, _ in rfree) + deeper = sum(y < x for x, y, _ in dmin) + values = [v for x, y, _ in rfree for v in (x, y)] + params = dict(out=out, n=n, nf=nf, nb=nb, dn=dn, dnf=dnf, dnb=dnb, within=within, band=RFREE_BAND, + deeper=deeper, shallower=sum(y > x for x, y, _ in dmin), + rlo=0.05 * math.floor((min(values) - 0.01) / 0.05), rhi=0.05 * math.ceil((max(values) + 0.01) / 0.05), + dlo=0.5, dhi=8.0) + script = "" + for term, ext in (('pdfcairo size 5in,5in font "Helvetica,14"', "pdf"), + ('pngcairo size 1500,1500 font "Helvetica,30"', "png")): + script += f"set terminal {term}\n" + GNUPLOT % dict(params, ext=ext) + "\nreset\n" + gp = os.path.join(out, "figures.gp") + open(gp, "w").write(script) + subprocess.run(["gnuplot", gp], check=True) + + with open(os.path.join(out, "summary.txt"), "w") as fh: + for path, m in manifests: + b = m.get("binary", {}) + fh.write(f"run {path}\n label {m.get('label')} binary git {b.get('source_head')} " + f"dirty {b.get('source_dirty')} sha256 {b.get('sha256')}\n flags {b.get('flags')} " + f"runner git {m.get('runner_git')} dirty {m.get('runner_dirty')}\n") + fh.write(f"R-free: {n} structures, {nf} facilities, {nb} beamlines; within +-{RFREE_BAND}: {within}\n" + f"d_min: {dn} structures, {dnf} facilities, {dnb} beamlines; Rugnux reaches further on " + f"{deeper}\nexcluded (set, figure, reason): excluded.txt\n") + print(open(os.path.join(out, "summary.txt")).read(), end="") + + +if __name__ == "__main__": + main()