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) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed
1 parent
d397717ff7
commit
07b0f7e315
2 files changed
+194
No files matched your search
@@ -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.
|
||||
|
||||
@@ -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()
|
||||
Reference in new issue
Block a user