diff --git a/rugnux_stills_ab.py b/rugnux_stills_ab.py new file mode 100755 index 00000000..aaafb277 --- /dev/null +++ b/rugnux_stills_ab.py @@ -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 | | 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: _report.txt (KEY= value) + .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 _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 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), + ("", "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%", "", "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) + # 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 >= 1)") + print() + print_shell_table(base_name, base, name, arm, floor) + print() + + print("=" * 120) + print() + + +if __name__ == "__main__": + main() diff --git a/stills_datasets.tsv b/stills_datasets.tsv new file mode 100644 index 00000000..529f5ecb --- /dev/null +++ b/stills_datasets.tsv @@ -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