Add a stills validation harness, the analogue of rugnux_vs_xds.py
rugnux_vs_xds.py gates every rotation change in this project. Nothing gates a stills change, so until now one could not be validated at all - which is why the background-ring width below had to be measured before it could be argued about. The rotation harness does not transfer. A stills dataset is one file of 100k-500k images with no per-crystal reference, the space group is known in advance and never in question, and what a change moves is the merge. So this asks a different question - did this make stills data better or worse - by running rugnux over the same fixed image subset under two or more argument sets and tabulating the merge per resolution shell. Three things it does on purpose. It pins the resolution range per dataset so both arms share shell edges: without that the automatic cut-off moves between arms and the columns are not the same shells, which happened on two of seven datasets. It fixes and prints the image subset, contiguous rather than strided because a stride turns one sequential read into ten thousand seeks. And --repeat measures the control floor, which on six of seven datasets is bit-identical in every column of every shell. It reads the stable interfaces - the report's KEY= value lines and the mmCIF _reflns block and _reflns_shell. loop. The one console-log number, how many predicted reflections lost their background ring, is labelled as such; it is there because it is the only place the price of a wider ring is counted. A CrystFEL stream can be scored alongside as an optional reference, keyed on the global image serial number - Event: //N restarts at zero in every file and silently multiplies the count. No dataset list is shipped. A dataset directory name identifies a sample, and so does a run label or a pump-probe parameter in a filename, so the table committed here is a template carrying the reasoning and no paths. A real one belongs in a config kept outside the repository, which is why --config is required and has no default. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
Claude Opus 5
Executable
+848
@@ -0,0 +1,848 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
rugnux_stills_ab.py -- the stills analogue of rugnux_vs_xds.py.
|
||||
|
||||
rugnux_vs_xds.py scores a ROTATION change against an XDS reference: one crystal per
|
||||
directory, one CORRECT.LP to compare with, space-group agreement as the headline. None of
|
||||
that transfers to serial stills. A stills dataset is one 100k-500k image file with no
|
||||
per-crystal reference, the space group is known in advance and never in question, and the
|
||||
thing a change moves is the merge: how many frames indexed, how many observations survived
|
||||
integration, and what the merged statistics look like PER SHELL.
|
||||
|
||||
So this script answers a different question: "did this change make stills data better or
|
||||
worse". It runs rugnux over the SAME fixed image subset of the same datasets under two (or
|
||||
more) sets of arguments and tabulates, per dataset and per resolution shell:
|
||||
|
||||
indexing rate | integrated observations | merged reflections | completeness |
|
||||
multiplicity | <I/sigma> | R_meas | CC1/2
|
||||
|
||||
and the arm-to-arm delta on each. Pooled numbers are printed too, but the per-shell table
|
||||
is the one that decides: a change that improves the pooled R_meas by moving the resolution
|
||||
cut-off has not improved anything, and pooled scoring has reversed the verdict on real
|
||||
changes before.
|
||||
|
||||
THREE THINGS IT DOES ON PURPOSE
|
||||
-------------------------------
|
||||
1. **The resolution range is pinned per dataset** (`dmin`/`dmax` in the config). Both arms
|
||||
are handed the same --scaling-high-resolution / --scaling-low-resolution, so the shell
|
||||
edges - uniform in 1/d^2 between those two limits - are identical between arms and the
|
||||
per-shell columns line up. Without this the automatic cc-logistic cut-off moves between
|
||||
arms and the shells are not the same shells.
|
||||
2. **The image subset is fixed and recorded** (`start`/`end`/`stride`). These datasets are
|
||||
95k-500k images; processing all of them to score a one-float change is not affordable.
|
||||
Every run prints the exact -s/-e/-t it used, so a later run reproduces it.
|
||||
3. **A control arm can be repeated** (--repeat N). Two runs of the SAME arm establish the
|
||||
noise floor. A delta smaller than that floor is not a result. rugnux's stills path uses
|
||||
a per-image min-pix search and a threaded merge, so the floor is not automatically zero.
|
||||
|
||||
CRYSTFEL REFERENCE (optional, --stream)
|
||||
---------------------------------------
|
||||
Where a CrystFEL .stream exists for the same file, `--stream` restricts it to exactly the
|
||||
same events (the stream records `Event: //N`) and reports CrystFEL's indexing rate and
|
||||
reflection count over that same subset. With --crystfel-merge it also filters the stream to
|
||||
those events, runs process_hkl (--even-only / --odd-only) and compare_hkl, and adds
|
||||
CrystFEL's per-shell R_split and CC1/2. That needs CrystFEL on PATH; it is a bonus, and an
|
||||
internal before/after is enough to gate a change without it.
|
||||
|
||||
CONFIG FILE
|
||||
-----------
|
||||
Whitespace-separated, '#' comments. One dataset per line:
|
||||
|
||||
alias master.h5 start end stride dmin dmax [extra rugnux args...]
|
||||
|
||||
`alias` is what appears in the table. Use a neutral one: dataset directory names identify
|
||||
samples, and nothing here should carry a sample identity.
|
||||
|
||||
USAGE
|
||||
-----
|
||||
./rugnux_stills_ab.py --config stills_datasets.tsv \
|
||||
--workdir /data/tmp/stills1/run1 --rugnux /path/to/rugnux \
|
||||
--arm base= --arm r3_14='--integration-radius 6,8,14' --repeat 2
|
||||
|
||||
Self-bootstrapping: on first run it makes a private venv, pip-installs gemmi, re-execs.
|
||||
"""
|
||||
|
||||
import os
|
||||
import sys
|
||||
import subprocess
|
||||
import pathlib
|
||||
import venv
|
||||
|
||||
_SELF = os.path.abspath(__file__)
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# venv bootstrap (gemmi is the only third-party dependency)
|
||||
# --------------------------------------------------------------------------- #
|
||||
def _venv_dir():
|
||||
if os.environ.get("RUGNUX_STILLS_VENV"):
|
||||
return pathlib.Path(os.environ["RUGNUX_STILLS_VENV"])
|
||||
cache = os.environ.get("XDG_CACHE_HOME") or os.path.join(os.path.expanduser("~"), ".cache")
|
||||
return pathlib.Path(cache) / "rugnux_stills_ab" / "venv"
|
||||
|
||||
|
||||
def _bootstrap():
|
||||
vdir = _venv_dir()
|
||||
py = vdir / "bin" / "python"
|
||||
if not py.exists():
|
||||
venv.create(vdir, with_pip=True)
|
||||
ok = subprocess.run([str(py), "-c", "import gemmi"],
|
||||
stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL).returncode == 0
|
||||
if not ok:
|
||||
subprocess.run([str(py), "-m", "pip", "install", "-q", "--disable-pip-version-check", "gemmi"],
|
||||
check=True)
|
||||
os.environ["RUGNUX_STILLS_BOOTSTRAPPED"] = "1"
|
||||
os.execv(str(py), [str(py), _SELF, *sys.argv[1:]])
|
||||
|
||||
|
||||
if __name__ == "__main__" and not os.environ.get("RUGNUX_STILLS_BOOTSTRAPPED"):
|
||||
_bootstrap()
|
||||
|
||||
# ---- past this point we are running inside the venv -----------------------
|
||||
import argparse
|
||||
import re
|
||||
import shlex
|
||||
import shutil
|
||||
import time
|
||||
import json
|
||||
import datetime
|
||||
import gemmi
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# config
|
||||
# --------------------------------------------------------------------------- #
|
||||
class Dataset:
|
||||
__slots__ = ("alias", "master", "start", "end", "stride", "dmin", "dmax", "extra")
|
||||
|
||||
def __init__(self, alias, master, start, end, stride, dmin, dmax, extra):
|
||||
self.alias, self.master = alias, pathlib.Path(master)
|
||||
self.start, self.end, self.stride = int(start), int(end), int(stride)
|
||||
self.dmin, self.dmax = float(dmin), float(dmax)
|
||||
self.extra = extra
|
||||
|
||||
@property
|
||||
def n_images(self):
|
||||
return max(0, (self.end - self.start + self.stride - 1) // self.stride)
|
||||
|
||||
def events(self):
|
||||
"""The image indices this subset covers -- rugnux -s/-e/-t semantics, 0-based,
|
||||
end exclusive. Used to restrict a CrystFEL stream to the same frames."""
|
||||
return range(self.start, self.end, self.stride)
|
||||
|
||||
|
||||
def read_config(path):
|
||||
out = []
|
||||
for lineno, line in enumerate(pathlib.Path(path).read_text().splitlines(), 1):
|
||||
line = line.split("#", 1)[0].strip()
|
||||
if not line:
|
||||
continue
|
||||
f = line.split()
|
||||
if len(f) < 7:
|
||||
sys.exit(f"{path}:{lineno}: need at least "
|
||||
f"'alias master start end stride dmin dmax'")
|
||||
out.append(Dataset(f[0], f[1], f[2], f[3], f[4], f[5], f[6], " ".join(f[7:])))
|
||||
return out
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# rugnux result: <prefix>_report.txt (KEY= value) + <prefix>.cif (per shell)
|
||||
# --------------------------------------------------------------------------- #
|
||||
def parse_report(path):
|
||||
"""The report's 'KEY= value' lines are a documented stable interface (REPORT_VERSION
|
||||
says when it last changed) -- unlike the console log, which is meant for humans."""
|
||||
r = {}
|
||||
if not path.exists():
|
||||
return r
|
||||
for line in path.read_text(errors="replace").splitlines():
|
||||
m = re.match(r"^([A-Z][A-Z0-9_]*)=\s*(.*)$", line)
|
||||
if m:
|
||||
r.setdefault(m.group(1), m.group(2).strip())
|
||||
return r
|
||||
|
||||
|
||||
def _f(v, scale=1.0):
|
||||
try:
|
||||
return float(v) * scale
|
||||
except (TypeError, ValueError):
|
||||
return None
|
||||
|
||||
|
||||
def _i(v):
|
||||
try:
|
||||
return int(float(v))
|
||||
except (TypeError, ValueError):
|
||||
return None
|
||||
|
||||
|
||||
SHELL_COLS = ["d_res_high", "d_res_low", "number_measured_obs", "number_unique_obs",
|
||||
"pdbx_redundancy", "percent_possible_obs", "meanI_over_sigI_obs",
|
||||
"pdbx_Rrim_I_all", "pdbx_CC_half"]
|
||||
|
||||
|
||||
def parse_cif(path):
|
||||
"""Overall + per-shell merging statistics out of rugnux's mmCIF (its stable output).
|
||||
Fractions are converted to percent so every ratio in the table reads the same way."""
|
||||
block = gemmi.cif.read(str(path)).sole_block()
|
||||
|
||||
def val(tag):
|
||||
v = block.find_value(tag)
|
||||
return gemmi.cif.as_string(v) if v is not None else None
|
||||
|
||||
r = {
|
||||
"sg": _i(val("_symmetry.Int_Tables_number")),
|
||||
"uniq": _i(val("_reflns.number_obs")),
|
||||
"obs": _i(val("_reflns.pdbx_number_measured_all")),
|
||||
"mult": _f(val("_reflns.pdbx_redundancy")),
|
||||
"compl": _f(val("_reflns.percent_possible_obs")),
|
||||
"isig": _f(val("_reflns.pdbx_netI_over_sigmaI")),
|
||||
"rmeas": _f(val("_reflns.pdbx_Rrim_I_all"), 100.0),
|
||||
"cc": _f(val("_reflns.pdbx_CC_half"), 100.0),
|
||||
"isa": _f(val("_reflns.jfjoch_diffrn_ISa")),
|
||||
"dmin": _f(val("_reflns.d_resolution_high")),
|
||||
"dmax": _f(val("_reflns.d_resolution_low")),
|
||||
}
|
||||
shells = []
|
||||
for row in block.find("_reflns_shell.", SHELL_COLS):
|
||||
d_hi, d_lo = _f(row[0]), _f(row[1])
|
||||
if d_hi is None:
|
||||
continue
|
||||
shells.append({
|
||||
"d_hi": d_hi, "d_lo": d_lo,
|
||||
"obs": _i(row[2]), "uniq": _i(row[3]),
|
||||
"mult": _f(row[4]), "compl": _f(row[5]), "isig": _f(row[6]),
|
||||
"rmeas": _f(row[7], 100.0), "cc": _f(row[8], 100.0),
|
||||
})
|
||||
shells.sort(key=lambda s: -s["d_hi"]) # low resolution first
|
||||
r["shells"] = shells
|
||||
return r
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# running rugnux
|
||||
# --------------------------------------------------------------------------- #
|
||||
def find_rugnux(explicit):
|
||||
for cand in (explicit, os.environ.get("RUGNUX")):
|
||||
if cand and (shutil.which(cand) or os.path.exists(cand)):
|
||||
return shutil.which(cand) or cand
|
||||
return shutil.which("rugnux")
|
||||
|
||||
|
||||
def run_arm(ds, arm_name, arm_args, workdir, rugnux_bin, threads, timeout, reuse, common_args):
|
||||
"""One rugnux run. Returns (result-dict-or-None, error, seconds)."""
|
||||
wd = workdir / ds.alias / arm_name
|
||||
wd.mkdir(parents=True, exist_ok=True)
|
||||
prefix = "out"
|
||||
cif = wd / f"{prefix}.cif"
|
||||
report = wd / f"{prefix}_report.txt"
|
||||
|
||||
cmd = [rugnux_bin, "-o", prefix,
|
||||
"-s", str(ds.start), "-e", str(ds.end), "-t", str(ds.stride),
|
||||
"--scaling-high-resolution", f"{ds.dmin:.3f}",
|
||||
"--scaling-low-resolution", f"{ds.dmax:.3f}"]
|
||||
if threads:
|
||||
cmd += ["-N", str(threads)]
|
||||
if common_args:
|
||||
cmd += shlex.split(common_args)
|
||||
if ds.extra:
|
||||
cmd += shlex.split(ds.extra)
|
||||
if arm_args:
|
||||
cmd += shlex.split(arm_args)
|
||||
cmd.append(str(ds.master))
|
||||
|
||||
if reuse and cif.exists() and report.exists():
|
||||
prev = json.loads((wd / "timing.json").read_text()) if (wd / "timing.json").exists() else {}
|
||||
return collect(cif, report, wd / "rugnux.log"), None, prev.get("seconds")
|
||||
|
||||
(wd / "cmd.txt").write_text(" ".join(shlex.quote(c) for c in cmd) + "\n")
|
||||
log = wd / "rugnux.log"
|
||||
t0 = time.monotonic()
|
||||
with open(log, "w") as lf:
|
||||
lf.write("$ " + " ".join(shlex.quote(c) for c in cmd) + "\n\n")
|
||||
lf.flush()
|
||||
try:
|
||||
p = subprocess.run(cmd, cwd=str(wd), stdout=lf, stderr=subprocess.STDOUT,
|
||||
timeout=timeout)
|
||||
except subprocess.TimeoutExpired:
|
||||
return None, f"timeout after {timeout}s", float(timeout)
|
||||
dt = time.monotonic() - t0
|
||||
(wd / "timing.json").write_text(json.dumps({"seconds": dt}))
|
||||
if p.returncode != 0:
|
||||
return None, f"exit {p.returncode} (see {log})", dt
|
||||
if not cif.exists():
|
||||
return None, f"no .cif produced (see {log})", dt
|
||||
return collect(cif, report, log), None, dt
|
||||
|
||||
|
||||
RING_RE = re.compile(
|
||||
r"Integration at r1=([\d.]+) r2=([\d.]+) r3=([\d.]+): of (\d+) predicted reflections, "
|
||||
r"([\d.]+)% lost their background ring - ([\d.]+)% to neighbouring reflections")
|
||||
|
||||
|
||||
def parse_ring_loss(log):
|
||||
"""The integrator's own ring-starvation line. This one comes from the CONSOLE log, not the
|
||||
report, so it is not a stable interface -- but it is the only place the reflections a wider
|
||||
ring COSTS are counted, split into the two causes that matter here: crowding (a neighbour
|
||||
inside the annulus) and the detector (edge, gap, mask). A reflection that loses its ring is
|
||||
discarded, and on stills there is no second view of it, so this is the price side of r3."""
|
||||
if not log.exists():
|
||||
return {}
|
||||
for line in reversed(log.read_text(errors="replace").splitlines()):
|
||||
m = RING_RE.search(line)
|
||||
if m:
|
||||
return {"radii": f"{m.group(1)}/{m.group(2)}/{m.group(3)}",
|
||||
"predicted": int(m.group(4)),
|
||||
"ring_lost_pct": float(m.group(5)),
|
||||
"ring_lost_neigh_pct": float(m.group(6))}
|
||||
return {}
|
||||
|
||||
|
||||
def parse_image_dat(path):
|
||||
"""Per-image scale health from <prefix>_image.dat.
|
||||
|
||||
This is here because it is what tells a real result from a spurious one. A stills merge is
|
||||
ScaleOnTheFly's per-image scale G, and images whose G collapses more than 50x below the run
|
||||
median are dropped from the merge entirely. On a dataset whose scaling has converged, the
|
||||
median G is order 1 and that filter takes almost nothing. On one where it has NOT, the median
|
||||
G is orders of magnitude off, hundreds of images sit against the threshold, and a change of a
|
||||
couple of per cent anywhere in the pipeline moves which ones fall through - which shows up as
|
||||
a large change in observation count and <I/sigma> that has nothing to do with the change
|
||||
being tested. Print both arms' numbers and the difference is obvious."""
|
||||
if not path.exists():
|
||||
return {}
|
||||
g, n = [], 0
|
||||
for line in path.read_text(errors="replace").splitlines():
|
||||
if line.startswith("#"):
|
||||
continue
|
||||
f = line.split()
|
||||
if len(f) < 2:
|
||||
continue
|
||||
n += 1
|
||||
try:
|
||||
v = float(f[1])
|
||||
except ValueError:
|
||||
continue
|
||||
if v == v and v > 0:
|
||||
g.append(v)
|
||||
if not g:
|
||||
return {"images_seen": n, "images_scaled": 0}
|
||||
g.sort()
|
||||
return {"images_seen": n, "images_scaled": len(g),
|
||||
"scale_median": g[len(g) // 2],
|
||||
"scale_p5": g[max(0, int(0.05 * len(g)) - 1)]}
|
||||
|
||||
|
||||
def collect(cif, report, log=None):
|
||||
r = parse_cif(cif)
|
||||
r.update(parse_image_dat(pathlib.Path(str(cif)[:-4] + "_image.dat")))
|
||||
rep = parse_report(report)
|
||||
r["index_rate"] = _f(rep.get("INDEXING_RATE"), 100.0)
|
||||
r["images"] = _i(rep.get("IMAGES_PROCESSED"))
|
||||
if r["index_rate"] is not None and r["images"]:
|
||||
r["indexed"] = int(round(r["index_rate"] / 100.0 * r["images"]))
|
||||
r["cell"] = rep.get("UNIT_CELL_CONSTANTS")
|
||||
r["sg_report"] = _i(rep.get("SPACE_GROUP_NUMBER"))
|
||||
if log is not None:
|
||||
r.update(parse_ring_loss(log))
|
||||
return r
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# CrystFEL reference
|
||||
# --------------------------------------------------------------------------- #
|
||||
# The four line prefixes a subset scan needs, plus the chunk delimiter that resets the state.
|
||||
# "----- Begin chunk" MUST be in here: a chunk whose event line is missing or unparseable would
|
||||
# otherwise inherit the previous chunk's membership and be counted against the wrong subset.
|
||||
_STREAM_PAT = (r"^(----- Begin chunk|Image serial number: |Event: //|indexed_by = "
|
||||
r"|--- Begin crystal|num_reflections = )")
|
||||
|
||||
|
||||
def _stream_lines(sp):
|
||||
"""These streams are 0.8-17 GB; a Python loop over every line of one is minutes of pure
|
||||
interpreter overhead for four line types. grep drops 99%+ of them first."""
|
||||
if shutil.which("grep"):
|
||||
p = subprocess.Popen(["grep", "-a", "-E", _STREAM_PAT, str(sp)],
|
||||
stdout=subprocess.PIPE, text=True, errors="replace")
|
||||
yield from p.stdout
|
||||
p.stdout.close()
|
||||
p.wait()
|
||||
else:
|
||||
with open(sp, "r", errors="replace") as f:
|
||||
for line in f:
|
||||
if re.match(_STREAM_PAT, line):
|
||||
yield line
|
||||
|
||||
|
||||
def stream_index(line, prev):
|
||||
"""Global 0-based image index of a chunk, or `prev` if this line does not carry one.
|
||||
|
||||
`Image serial number` is the image's position in indexamajig's input list and is the key to
|
||||
use: it is global. `Event: //N` is NOT - where the list names the per-file data h5s rather
|
||||
than the master, N restarts at 0 in every file, and keying on it silently merges frame 5 of
|
||||
every one of 270 files into a single "event 5". (That is how a 10 000-frame subset first came
|
||||
back with 231 448 indexed frames.) Serial numbers are 1-based; Event is 0-based, and is kept
|
||||
only as a fallback for a stream that lacks the serial line.
|
||||
|
||||
This makes one assumption, and it is worth checking once per stream: that indexamajig's input
|
||||
list held the dataset's frames in order and in full, so that position-in-list is the frame
|
||||
number. Verify it by finding a chunk from the SECOND data file and checking that its serial
|
||||
equals its event plus the number of frames in the first file.
|
||||
"""
|
||||
if line.startswith("Image serial number: "):
|
||||
try:
|
||||
return int(line[21:].strip()) - 1
|
||||
except ValueError:
|
||||
return prev
|
||||
if line.startswith("Event: //") and prev is None:
|
||||
try:
|
||||
return int(line[9:].strip())
|
||||
except ValueError:
|
||||
return prev
|
||||
return prev
|
||||
|
||||
|
||||
def scan_stream(stream_paths, events):
|
||||
"""Indexing rate + reflection count over exactly the frames in `events`.
|
||||
|
||||
Chunks outside the subset are ignored, which is the whole point: comparing a 10 000-image
|
||||
rugnux run against a 269 307-chunk stream compares nothing.
|
||||
"""
|
||||
want = set(events)
|
||||
seen, tot = set(), [0, 0, 0] # indexed frames, crystals, reflections
|
||||
# One chunk is accumulated into `cur` and only added to the totals once the chunk ends, so
|
||||
# the counters never depend on whether the index line came before or after them.
|
||||
idx, cur = None, [0, 0, 0]
|
||||
|
||||
def flush():
|
||||
if idx is not None and idx in want:
|
||||
seen.add(idx)
|
||||
for i in range(3):
|
||||
tot[i] += cur[i]
|
||||
|
||||
for sp in stream_paths:
|
||||
for line in _stream_lines(sp):
|
||||
if line.startswith("----- Begin chunk"):
|
||||
flush()
|
||||
idx, cur = None, [0, 0, 0]
|
||||
elif line.startswith(("Image serial number: ", "Event: //")):
|
||||
idx = stream_index(line, idx)
|
||||
elif line.startswith("indexed_by = "):
|
||||
cur[0] += (line[13:].strip() != "none")
|
||||
elif line.startswith("--- Begin crystal"):
|
||||
cur[1] += 1
|
||||
elif line.startswith("num_reflections = "):
|
||||
cur[2] += int(line[18:].strip() or 0)
|
||||
flush()
|
||||
idx, cur = None, [0, 0, 0]
|
||||
return {"frames": len(seen), "indexed": tot[0], "crystals": tot[1], "obs": tot[2],
|
||||
"index_rate": (100.0 * tot[0] / len(seen)) if seen else None}
|
||||
|
||||
|
||||
def filter_stream(stream_paths, events, out_path):
|
||||
"""Write a stream holding only the chunks whose event is in `events` (header kept), so
|
||||
process_hkl merges exactly the frames rugnux processed."""
|
||||
want = set(events)
|
||||
with open(out_path, "w") as out:
|
||||
for i, sp in enumerate(stream_paths):
|
||||
with open(sp, "r", errors="replace") as f:
|
||||
buf, in_chunk, idx, seen_chunk = [], False, None, False
|
||||
for line in f:
|
||||
if line.startswith("----- Begin chunk"):
|
||||
in_chunk, buf, idx, seen_chunk = True, [line], None, True
|
||||
continue
|
||||
if not in_chunk:
|
||||
if i == 0 and not seen_chunk:
|
||||
out.write(line) # the header, once
|
||||
continue
|
||||
buf.append(line)
|
||||
if line.startswith(("Image serial number: ", "Event: //")):
|
||||
idx = stream_index(line, idx)
|
||||
elif line.startswith("----- End chunk"):
|
||||
if idx in want:
|
||||
out.writelines(buf)
|
||||
in_chunk, buf = False, []
|
||||
return out_path
|
||||
|
||||
|
||||
def crystfel_merge(sub_stream, wd, sym, cell, dmin, dmax, nshells):
|
||||
"""process_hkl (even / odd half-sets) + compare_hkl -> per-shell R_split and CC1/2.
|
||||
|
||||
compare_hkl needs the unit cell (`-p`) to put reflections in resolution shells at all, so a
|
||||
.cell file is required for this path; the indexamajig command line at the top of the stream
|
||||
names the one that produced it."""
|
||||
def hkl(tag, *extra):
|
||||
p = wd / f"cf_{tag}.hkl"
|
||||
cmd = ["process_hkl", "-i", str(sub_stream), "-o", str(p), "-y", sym, "--scale", *extra]
|
||||
subprocess.run(cmd, cwd=str(wd), stdout=subprocess.DEVNULL,
|
||||
stderr=open(wd / f"cf_{tag}.log", "w"), check=True)
|
||||
return p
|
||||
a, b = hkl("even", "--even-only"), hkl("odd", "--odd-only")
|
||||
out = {}
|
||||
for fom in ("Rsplit", "CC"):
|
||||
shellf = wd / f"cf_{fom}.dat"
|
||||
cmd = ["compare_hkl", str(a), str(b), "-y", sym, "-p", str(cell), f"--fom={fom}",
|
||||
"--shell-file=" + str(shellf), f"--nshells={nshells}",
|
||||
f"--highres={dmin}", f"--lowres={dmax}"]
|
||||
try:
|
||||
subprocess.run(cmd, cwd=str(wd), stdout=subprocess.DEVNULL,
|
||||
stderr=open(wd / f"cf_cmp_{fom}.log", "w"), check=True)
|
||||
except (subprocess.CalledProcessError, FileNotFoundError):
|
||||
continue
|
||||
# compare_hkl's shell file has a '#'-commented header; column 0 is the shell centre in
|
||||
# 1/d (nm^-1) and column 1 is the figure of merit. Converted to d (Angstrom) here so the
|
||||
# rows line up with every other resolution column in this report.
|
||||
rows = []
|
||||
for line in shellf.read_text().splitlines():
|
||||
if line.lstrip().startswith("#"):
|
||||
continue
|
||||
f = line.split()
|
||||
if len(f) >= 2 and re.match(r"^[\d.]+$", f[0]) and float(f[0]) > 0:
|
||||
rows.append((10.0 / float(f[0]), float(f[1])))
|
||||
out[fom] = rows
|
||||
return out
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# rendering
|
||||
# --------------------------------------------------------------------------- #
|
||||
def n(v, w, dec=1):
|
||||
return "-".rjust(w) if v is None else f"{v:.{dec}f}".rjust(w)
|
||||
|
||||
|
||||
def ni(v, w):
|
||||
return "-".rjust(w) if v is None else f"{v:,}".rjust(w)
|
||||
|
||||
|
||||
def delta(a, b, w, dec=2):
|
||||
if a is None or b is None:
|
||||
return "-".rjust(w)
|
||||
return f"{b - a:+.{dec}f}".rjust(w)
|
||||
|
||||
|
||||
SUMMARY_COLS = [
|
||||
("index%", "index_rate", 7, 2),
|
||||
("Obs", "obs", 11, 0),
|
||||
("Refl", "uniq", 8, 0),
|
||||
("Mult", "mult", 6, 2),
|
||||
("Compl%", "compl", 7, 1),
|
||||
("<I/sig>", "isig", 8, 2),
|
||||
("Rmeas%", "rmeas", 7, 2),
|
||||
("CC1/2%", "cc", 7, 2),
|
||||
("ISa", "isa", 6, 2),
|
||||
]
|
||||
|
||||
|
||||
def summary_header():
|
||||
h = f" {'arm':<12}"
|
||||
for name, _, w, _ in SUMMARY_COLS:
|
||||
h += " " + name.rjust(w)
|
||||
return h + " " + "time".rjust(8)
|
||||
|
||||
|
||||
def summary_row(label, r, secs):
|
||||
line = f" {label:<12}"
|
||||
for _, key, w, dec in SUMMARY_COLS:
|
||||
v = r.get(key) if r else None
|
||||
line += " " + (ni(v, w) if key in ("obs", "uniq") else n(v, w, dec))
|
||||
line += " " + ((f"{secs:.0f}s").rjust(8) if secs is not None else "-".rjust(8))
|
||||
return line
|
||||
|
||||
|
||||
SHELL_KEYS = [("obs", 9, 0), ("uniq", 7, 0), ("mult", 6, 2), ("compl", 7, 1),
|
||||
("isig", 8, 2), ("rmeas", 8, 2), ("cc", 7, 2)]
|
||||
SHELL_NAMES = ["Obs", "Refl", "Mult", "Compl%", "<I/sig>", "Rmeas%", "CC1/2%"]
|
||||
|
||||
|
||||
def print_shell_table(base_name, base, arm_name, arm, floor=None):
|
||||
"""Per shell: base | arm | delta. The delta line flags whether it clears the control
|
||||
floor measured from repeated runs of the base arm."""
|
||||
hdr = f" {'d_hi':>6} {'d_lo':>6} {'arm':>7}"
|
||||
for nm, (_, w, _) in zip(SHELL_NAMES, SHELL_KEYS):
|
||||
hdr += " " + nm.rjust(w)
|
||||
print(hdr)
|
||||
print(" " + "-" * (len(hdr) - 4))
|
||||
for i, sb in enumerate(base["shells"]):
|
||||
sa = arm["shells"][i] if i < len(arm["shells"]) else {}
|
||||
if abs(sb["d_hi"] - sa.get("d_hi", -1)) > 1e-3:
|
||||
print(f" SHELL EDGES DIFFER at {sb['d_hi']:.2f} -- resolution range not pinned?")
|
||||
return
|
||||
sig = " *" if (sb.get("isig") or 0) >= 1.0 else ""
|
||||
for tag, s in ((base_name, sb), (arm_name, sa)):
|
||||
line = f" {s['d_hi']:6.2f} {s['d_lo']:6.2f} {tag[:7]:>7}"
|
||||
for key, w, dec in SHELL_KEYS:
|
||||
v = s.get(key)
|
||||
line += " " + (ni(v, w) if key in ("obs", "uniq") else n(v, w, dec))
|
||||
print(line + (sig if tag == base_name else ""))
|
||||
line = f" {'':6} {'':6} {'delta':>7}"
|
||||
for key, w, dec in SHELL_KEYS:
|
||||
a, b = sb.get(key), sa.get(key)
|
||||
if key in ("obs", "uniq"):
|
||||
line += " " + ("-".rjust(w) if a is None or b is None
|
||||
else f"{b - a:+,}".rjust(w))
|
||||
else:
|
||||
line += " " + delta(a, b, w, dec)
|
||||
# <I/sigma> as a percentage, because that is the quantity a variance change moves and a
|
||||
# fixed number of counts means something different in the first shell and the last.
|
||||
a, b = sb.get("isig"), sa.get("isig")
|
||||
line += " " + (" -" if not a or b is None else f"{100.0 * (b / a - 1):+6.2f}%")
|
||||
if floor:
|
||||
flags = []
|
||||
for key, _, _ in SHELL_KEYS:
|
||||
fl = floor.get(i, {}).get(key)
|
||||
a, b = sb.get(key), sa.get(key)
|
||||
if fl is not None and a is not None and b is not None and abs(b - a) > fl:
|
||||
flags.append(key)
|
||||
line += " " + (">floor: " + ",".join(flags) if flags else "(within floor)")
|
||||
print(line)
|
||||
print()
|
||||
|
||||
|
||||
def control_floor(reps):
|
||||
"""Per-shell spread of repeated runs of the SAME arm -- max|x_i - x_j| per column.
|
||||
A delta at or below this is not a result. Returns {shell_index: {key: floor}}."""
|
||||
if len(reps) < 2:
|
||||
return {}
|
||||
out = {}
|
||||
nsh = min(len(r["shells"]) for r in reps)
|
||||
for i in range(nsh):
|
||||
d = {}
|
||||
for key, _, _ in SHELL_KEYS:
|
||||
vals = [r["shells"][i].get(key) for r in reps]
|
||||
vals = [v for v in vals if v is not None]
|
||||
d[key] = (max(vals) - min(vals)) if len(vals) > 1 else None
|
||||
out[i] = d
|
||||
return out
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# main
|
||||
# --------------------------------------------------------------------------- #
|
||||
def main():
|
||||
ap = argparse.ArgumentParser(
|
||||
description="Score a rugnux change on serial-stills data, per resolution shell.")
|
||||
ap.add_argument("--config", required=True, help="dataset table (see the docstring)")
|
||||
ap.add_argument("--workdir", required=True, help="where rugnux outputs go")
|
||||
ap.add_argument("--rugnux", help="rugnux binary (else $RUGNUX / PATH)")
|
||||
ap.add_argument("--arm", action="append", default=[], metavar="NAME=ARGS",
|
||||
help="an arm: a label and the extra rugnux arguments. Repeat. The FIRST "
|
||||
"arm is the baseline everything else is compared with. "
|
||||
"e.g. --arm base= --arm r3_14='--integration-radius 6,8,14'")
|
||||
ap.add_argument("--common-args", default="",
|
||||
help="arguments appended to EVERY arm (e.g. a fixed -S)")
|
||||
ap.add_argument("--repeat", type=int, default=1,
|
||||
help="run the BASELINE arm this many times to measure the control floor "
|
||||
"(default 1 = no floor). 2 is enough to say whether a delta is real.")
|
||||
ap.add_argument("--threads", type=int, help="rugnux -N")
|
||||
ap.add_argument("--timeout", type=int, default=14400, help="per run, seconds")
|
||||
ap.add_argument("--only", help="comma-separated aliases to run")
|
||||
ap.add_argument("--reuse", action="store_true", help="skip a run whose .cif already exists")
|
||||
ap.add_argument("--stream", action="append", default=[], metavar="ALIAS=GLOB",
|
||||
help="CrystFEL stream(s) for a dataset, for an indexing-rate reference "
|
||||
"over the same frames. Repeat per dataset.")
|
||||
ap.add_argument("--crystfel-merge", action="store_true",
|
||||
help="also merge the subset stream with process_hkl and compare_hkl "
|
||||
"(needs CrystFEL on PATH); adds R_split / CC1/2 per shell")
|
||||
ap.add_argument("--crystfel-sym", default="mmm",
|
||||
help="CrystFEL point group for process_hkl -y (default mmm)")
|
||||
ap.add_argument("--crystfel-cell", action="append", default=[], metavar="ALIAS=FILE.cell",
|
||||
help="unit-cell file for compare_hkl -p, per dataset. Required by "
|
||||
"--crystfel-merge; the indexamajig line at the top of the stream names "
|
||||
"the file that produced it")
|
||||
ap.add_argument("--shells", type=int, default=10, help="resolution shells (default 10)")
|
||||
ap.add_argument("--progress", action="store_true")
|
||||
args = ap.parse_args()
|
||||
|
||||
if not args.arm:
|
||||
sys.exit("give at least one --arm NAME=ARGS (the first is the baseline)")
|
||||
arms = []
|
||||
for a in args.arm:
|
||||
name, _, rest = a.partition("=")
|
||||
arms.append((name.strip(), rest.strip()))
|
||||
|
||||
rugnux_bin = find_rugnux(args.rugnux)
|
||||
if not rugnux_bin:
|
||||
sys.exit("rugnux binary not found -- pass --rugnux PATH or set $RUGNUX")
|
||||
|
||||
streams = {}
|
||||
for s in args.stream:
|
||||
alias, _, g = s.partition("=")
|
||||
g = g.strip()
|
||||
streams[alias.strip()] = sorted(pathlib.Path(g).parent.glob(pathlib.Path(g).name))
|
||||
cells = {}
|
||||
for s in args.crystfel_cell:
|
||||
alias, _, p = s.partition("=")
|
||||
cells[alias.strip()] = pathlib.Path(p.strip())
|
||||
|
||||
datasets = read_config(args.config)
|
||||
if args.only:
|
||||
want = {x.strip() for x in args.only.split(",")}
|
||||
datasets = [d for d in datasets if d.alias in want]
|
||||
workdir = pathlib.Path(args.workdir)
|
||||
common = (args.common_args + f" --resolution-shells {args.shells}").strip()
|
||||
|
||||
print()
|
||||
print(f" rugnux stills A/B . {datetime.date.today().isoformat()}")
|
||||
print(f" rugnux : {rugnux_bin}")
|
||||
print(f" arms : " + " | ".join(f"{nm} [{ar or 'defaults'}]" for nm, ar in arms))
|
||||
print(f" common : {common or '(none)'}")
|
||||
print(f" repeat : baseline x{args.repeat} (control floor)")
|
||||
print()
|
||||
|
||||
all_results = []
|
||||
for ds in datasets:
|
||||
if args.progress:
|
||||
print(f"[{ds.alias}] {ds.n_images} images "
|
||||
f"(-s {ds.start} -e {ds.end} -t {ds.stride}), "
|
||||
f"{ds.dmax}-{ds.dmin} A", file=sys.stderr, flush=True)
|
||||
per_arm = {}
|
||||
for ai, (name, arm_args) in enumerate(arms):
|
||||
nrep = args.repeat if ai == 0 else 1
|
||||
reps = []
|
||||
for k in range(nrep):
|
||||
label = name if k == 0 else f"{name}#{k + 1}"
|
||||
if args.progress:
|
||||
print(f" arm {label} ...", file=sys.stderr, flush=True)
|
||||
r, err, secs = run_arm(ds, label, arm_args, workdir, rugnux_bin,
|
||||
args.threads, args.timeout, args.reuse, common)
|
||||
if args.progress:
|
||||
print(f" arm {label} done ({err or 'ok'})", file=sys.stderr, flush=True)
|
||||
reps.append((label, r, err, secs))
|
||||
per_arm[name] = reps
|
||||
all_results.append((ds, per_arm))
|
||||
|
||||
# -------------------------------- report --------------------------------
|
||||
for ds, per_arm in all_results:
|
||||
print("=" * 120)
|
||||
print(f"* {ds.alias} images {ds.start}..{ds.end} step {ds.stride} "
|
||||
f"({ds.n_images} frames) {ds.dmax:.1f}-{ds.dmin:.2f} A")
|
||||
if ds.extra:
|
||||
print(f" dataset args: {ds.extra}")
|
||||
print()
|
||||
print(summary_header())
|
||||
print(" " + "-" * (len(summary_header()) - 2))
|
||||
for name, reps in per_arm.items():
|
||||
for label, r, err, secs in reps:
|
||||
if err:
|
||||
print(f" {label:<12} FAILED: {err}")
|
||||
else:
|
||||
print(summary_row(label, r, secs))
|
||||
print()
|
||||
print(f" {'arm':<12} {'radii':>10} {'predicted':>12} {'ring lost':>10} "
|
||||
f"{'to neighbours':>14} (the price side of r3: a reflection whose ring is gone "
|
||||
f"is discarded)")
|
||||
for name, reps in per_arm.items():
|
||||
for label, r, err, _ in reps:
|
||||
if err or not r or "ring_lost_pct" not in r:
|
||||
continue
|
||||
print(f" {label:<12} {r['radii']:>10} {r['predicted']:>12,} "
|
||||
f"{r['ring_lost_pct']:9.3f}% {r['ring_lost_neigh_pct']:13.3f}%")
|
||||
print()
|
||||
print(f" {'arm':<12} {'images scaled':>14} {'median G':>10} {'5th pct G':>11} "
|
||||
f"(a median G far from 1, or an arm-to-arm change in 'images scaled', means the")
|
||||
print(f" {'':<12} {'':>14} {'':>10} {'':>11} "
|
||||
f"per-image scaling - not this change - is moving the merge; see parse_image_dat)")
|
||||
warn_scale = False
|
||||
for name, reps in per_arm.items():
|
||||
for label, r, err, _ in reps:
|
||||
if err or not r or "scale_median" not in r:
|
||||
continue
|
||||
print(f" {label:<12} {r['images_scaled']:>14,} {r['scale_median']:>10.4f} "
|
||||
f"{r['scale_p5']:>11.5f}")
|
||||
if not (0.05 < r["scale_median"] < 20.0):
|
||||
warn_scale = True
|
||||
if warn_scale:
|
||||
print(" !! SCALING SUSPECT on this dataset: the median per-image scale is orders of "
|
||||
"magnitude from 1.")
|
||||
print(" !! Its merge is set by which images survive the 50x-collapse filter, so the "
|
||||
"deltas below")
|
||||
print(" !! measure that filter, not the change under test. Do not score this dataset.")
|
||||
sp = streams.get(ds.alias)
|
||||
if sp:
|
||||
cf = scan_stream(sp, ds.events())
|
||||
print(f" {'CrystFEL':<12} " + n(cf["index_rate"], 7, 2) +
|
||||
" " + ni(cf["obs"], 11) + f" ({cf['crystals']:,} crystals over "
|
||||
f"{cf['frames']:,} of the same frames; --int-radius=2,3,5)")
|
||||
print(" CrystFEL's reflection count is every predicted reflection it wrote, before "
|
||||
"any merge cut; rugnux's")
|
||||
print(" Obs is what survived into the merge. The two indexing rates ARE comparable "
|
||||
"- same frames, same index.")
|
||||
if args.crystfel_merge:
|
||||
wd = workdir / ds.alias / "crystfel"
|
||||
wd.mkdir(parents=True, exist_ok=True)
|
||||
sub = wd / "subset.stream"
|
||||
if not sub.exists():
|
||||
filter_stream(sp, ds.events(), sub)
|
||||
cell = cells.get(ds.alias)
|
||||
if not cell:
|
||||
print(" CrystFEL merge needs --crystfel-cell ALIAS=file.cell (compare_hkl "
|
||||
"cannot make shells without it)")
|
||||
cfm = {}
|
||||
else:
|
||||
try:
|
||||
cfm = crystfel_merge(sub, wd, args.crystfel_sym, cell,
|
||||
ds.dmin, ds.dmax, args.shells)
|
||||
except (subprocess.CalledProcessError, FileNotFoundError) as e:
|
||||
cfm = {}
|
||||
print(f" CrystFEL merge unavailable: {e}")
|
||||
if cfm:
|
||||
print()
|
||||
print(f" CrystFEL half-set FOMs over the same frames "
|
||||
f"(process_hkl --even-only/--odd-only, -y {args.crystfel_sym}):")
|
||||
print(f" {'d (A)':>8} {'Rsplit%':>9} {'CC1/2':>8}")
|
||||
rs = dict(cfm.get("Rsplit", []))
|
||||
cc = dict(cfm.get("CC", []))
|
||||
for d in sorted(set(rs) | set(cc), reverse=True):
|
||||
print(f" {d:8.2f} {rs.get(d, float('nan')):9.2f} "
|
||||
f"{cc.get(d, float('nan')):8.4f}")
|
||||
print()
|
||||
|
||||
base_name = arms[0][0]
|
||||
base_reps = [r for _, r, e, _ in per_arm[base_name] if r]
|
||||
if not base_reps:
|
||||
print(" baseline failed; nothing to compare")
|
||||
continue
|
||||
floor = control_floor(base_reps)
|
||||
all_zero = floor and all(not v for d in floor.values() for v in d.values() if v is not None)
|
||||
if floor and args.repeat > 1 and all_zero:
|
||||
print(f" CONTROL FLOOR: exactly zero in every column of every shell over "
|
||||
f"{len(base_reps)} runs of '{base_name}'.")
|
||||
print(f" This arm is bit-reproducible, so every delta below is a real difference in "
|
||||
f"the code path -")
|
||||
print(f" '>floor' then means 'reproducibly different', NOT 'large'. Read the sizes, "
|
||||
f"not the flags.")
|
||||
print()
|
||||
elif floor and args.repeat > 1:
|
||||
print(f" CONTROL FLOOR (spread over {len(base_reps)} runs of "
|
||||
f"'{base_name}', per shell, max-min):")
|
||||
hdr = f" {'d_hi':>6}"
|
||||
for nm, (_, w, _) in zip(SHELL_NAMES, SHELL_KEYS):
|
||||
hdr += " " + nm.rjust(w)
|
||||
print(hdr)
|
||||
for i in sorted(floor):
|
||||
line = f" {base_reps[0]['shells'][i]['d_hi']:6.2f}"
|
||||
for key, w, dec in SHELL_KEYS:
|
||||
v = floor[i].get(key)
|
||||
line += " " + (ni(None if v is None else int(v), w)
|
||||
if key in ("obs", "uniq") else n(v, w, dec))
|
||||
print(line)
|
||||
print()
|
||||
|
||||
base = base_reps[0]
|
||||
for name, reps in per_arm.items():
|
||||
if name == base_name:
|
||||
continue
|
||||
arm = reps[0][1]
|
||||
if not arm:
|
||||
continue
|
||||
print(f" PER SHELL {base_name} vs {name} "
|
||||
f"(* = signal-bearing, baseline <I/sigma> >= 1)")
|
||||
print()
|
||||
print_shell_table(base_name, base, name, arm, floor)
|
||||
print()
|
||||
|
||||
print("=" * 120)
|
||||
print()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,30 @@
|
||||
# Dataset table for rugnux_stills_ab.py -- TEMPLATE.
|
||||
#
|
||||
# Columns:
|
||||
# alias master.h5 start end stride dmin dmax [extra rugnux args]
|
||||
#
|
||||
# This file ships with NO real dataset paths, and none should be added to it. A dataset
|
||||
# directory name identifies a sample, and so does a run label or a pump-probe parameter in a
|
||||
# filename. Keep your own table outside the repository and point the harness at it:
|
||||
#
|
||||
# ./rugnux_stills_ab.py --config ~/my_stills_datasets.tsv --workdir /somewhere ...
|
||||
#
|
||||
# --config is a required argument with no default, so nothing depends on this file existing.
|
||||
#
|
||||
# Filling a row:
|
||||
# alias a neutral label; it is what appears in the results table. Not the sample.
|
||||
# dmin pin it, do not leave it to the automatic cut-off. Run a calibration pass on the
|
||||
# exact subset first and take the FINER of the two arms' cc-logistic cut-offs, so the
|
||||
# pin clips neither arm and the shell edges - uniform in 1/d^2 between dmax and dmin -
|
||||
# are identical in both. Without this the cut moves between arms and the columns being
|
||||
# compared are not the same shells.
|
||||
# dmax 50, the value XDS configurations use and rugnux's own --scaling-low-resolution default.
|
||||
# subset contiguous (stride 1), not strided: on a spinning disk a stride turns one sequential
|
||||
# read into ten thousand seeks.
|
||||
# -C/-S pin the cell and space group so neither arm can differ by having indexed a different
|
||||
# lattice. For plain lysozyme, the field's standard test specimen, the reference cell
|
||||
# 79.1,79.1,37.9,90,90,90 in P4(3)2(1)2 (-S 96) is the usual choice.
|
||||
#
|
||||
# Example row, with the path left as a placeholder to be filled in locally:
|
||||
#
|
||||
# still_a /path/to/your/dataset_master.h5 0 10000 1 1.66 50 -C 79.1,79.1,37.9,90,90,90 -S 96
|
||||
|
Reference in New Issue
Block a user