tools: add rugnux_vs_xds.py — rugnux vs XDS regression table

Self-contained script (self-bootstrapping venv + gemmi) that runs rugnux
de-novo over a directory of test crystals and tabulates its merging statistics
against the XDS CORRECT.LP reference: number of reflections, observations,
space group (H-M + number), R_meas / CC1/2 overall, R_meas low shell,
CC1/2 high shell, ISa, and rugnux wall-clock time — with a per-crystal
point-group-level space-group match (OK/DIFF).

rugnux stats are read from its mmCIF output (the stable interface, not the
console log); resolution limit and Friedel/anomalous setting are matched to
XDS for a like-for-like comparison. Default output dir rugnux_cmp/ is ignored.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

.gitignore

@@ -3,3 +3,6 @@ cmake-build-release/
 build*/
 python-client/
 openapi-generator-cli.jar
+
+# rugnux_vs_xds.py default output dir
+rugnux_cmp/

rugnux_vs_xds.py

@@ -0,0 +1,399 @@
+#!/usr/bin/env python3
+"""
+rugnux_vs_xds.py -- run rugnux over a directory of test crystals and tabulate its
+merging statistics against the XDS reference (CORRECT.LP) next to each dataset.
+
+For every <data-root>/<crystal>/ that has both a *_master.h5 and an XDS CORRECT.LP it:
+  * reads the XDS reference (space group, resolution range, Friedel setting, merge
+    stats, ISa);
+  * runs rugnux de-novo with the SAME high-resolution limit and Friedel/anomalous
+    setting as XDS, so the two are comparable -- the space group is NOT forced, it is
+    the thing we compare;
+  * parses rugnux's own mmCIF output (the stable interface -- NOT the console log,
+    which is meant for humans and may be trimmed/reformatted);
+  * prints one compact table: reflections, observations, space group (H-M + number),
+    R_meas / CC1/2 overall, R_meas low shell, CC1/2 high shell, ISa -- XDS vs rugnux.
+
+Self-contained: on first run it builds a private venv, pip-installs gemmi and re-execs.
+All rugnux console output goes to a per-crystal log; only the table is printed.
+"""
+
+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_CMP_VENV"):
+        return pathlib.Path(os.environ["RUGNUX_CMP_VENV"])
+    cache = os.environ.get("XDG_CACHE_HOME") or os.path.join(os.path.expanduser("~"), ".cache")
+    return pathlib.Path(cache) / "rugnux_vs_xds" / "venv"
+
+
+def _bootstrap():
+    vdir = _venv_dir()
+    py = vdir / "bin" / "python"
+    if not py.exists():
+        venv.create(vdir, with_pip=True)
+    have_gemmi = subprocess.run([str(py), "-c", "import gemmi"],
+                                stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL).returncode == 0
+    if not have_gemmi:
+        subprocess.run([str(py), "-m", "pip", "install", "-q", "--disable-pip-version-check", "gemmi"],
+                       check=True)
+    os.environ["RUGNUX_CMP_BOOTSTRAPPED"] = "1"
+    os.execv(str(py), [str(py), _SELF, *sys.argv[1:]])
+
+
+if __name__ == "__main__" and not os.environ.get("RUGNUX_CMP_BOOTSTRAPPED"):
+    _bootstrap()
+
+# ---- past this point we are running inside the venv -----------------------
+import argparse
+import re
+import shutil
+import time
+import datetime
+import gemmi
+
+
+# --------------------------------------------------------------------------- #
+#  XDS reference (CORRECT.LP)
+# --------------------------------------------------------------------------- #
+def parse_xds(correct_lp):
+    """Pull the final merge statistics out of an XDS CORRECT.LP."""
+    txt = pathlib.Path(correct_lp).read_text(errors="replace")
+    r = {}
+
+    # The first SPACE_GROUP_NUMBER is only the INTEGRATE-step group (often P1); the LAST one is
+    # the symmetry CORRECT actually merged/scaled in and wrote to XDS_ASCII.HKL -- that is the
+    # XDS space group we compare (and the one the stats below are computed in).
+    sgs = re.findall(r"SPACE_GROUP_NUMBER=\s*(\d+)", txt)
+    r["sg"] = int(sgs[-1]) if sgs else None
+
+    m = re.search(r"FRIEDEL'S_LAW=\s*(TRUE|FALSE)", txt)
+    r["anomalous"] = (m.group(1) == "FALSE") if m else False
+
+    m = re.search(r"INCLUDE_RESOLUTION_RANGE=\s*[\d.]+\s+([\d.]+)", txt)
+    include_high = float(m.group(1)) if m else 0.0
+
+    m = re.search(r"^\s*a\s+b\s+ISa\s*\n\s*[\d.Ee+-]+\s+[\d.Ee+-]+\s+([\d.]+)", txt, re.M)
+    r["isa"] = float(m.group(1)) if m else None
+
+    # Final "... AS FUNCTION OF RESOLUTION" table (last occurrence); rows run from the
+    # lowest-resolution shell down to the highest, then a "total" line.
+    seg = txt[txt.rfind("AS FUNCTION OF RESOLUTION"):]
+    num = lambda x: float(x.rstrip("*%"))
+    shells, total = [], None
+    for line in seg.splitlines():
+        f = line.split()
+        if not f:
+            continue
+        if f[0] == "total" and len(f) >= 11:
+            total = f
+            break
+        if re.match(r"^\d+\.\d+$", f[0]) and len(f) >= 11:
+            shells.append(f)
+    if total:
+        r["obs"], r["uniq"] = int(total[1]), int(total[2])
+        r["rmeas"], r["cc"] = num(total[9]), num(total[10])
+    if shells:
+        r["rmeas_lo"] = num(shells[0][9])    # first row  = lowest resolution
+        r["cc_hi"] = num(shells[-1][10])     # last row   = highest resolution
+    # Effective XDS high-resolution limit to hand to rugnux: prefer the explicit
+    # INCLUDE_RESOLUTION_RANGE, but fall back to the finest shell actually tabulated
+    # (XDS.INP often leaves the high limit at 0.0 = "use the full detector range").
+    table_high = float(shells[-1][0]) if shells else 0.0
+    r["dmin"] = (include_high if include_high > 0 else table_high) or None
+    return r
+
+
+# --------------------------------------------------------------------------- #
+#  rugnux result (mmCIF)
+# --------------------------------------------------------------------------- #
+def parse_rugnux_cif(cif_path):
+    """Pull the same statistics out of rugnux's mmCIF (fractions -> percent)."""
+    block = gemmi.cif.read(str(cif_path)).sole_block()
+
+    def val(tag):
+        v = block.find_value(tag)
+        return gemmi.cif.as_string(v) if v is not None else None
+
+    def fnum(tag, scale=1.0):
+        v = val(tag)
+        try:
+            return float(v) * scale
+        except (TypeError, ValueError):
+            return None
+
+    def inum(tag):
+        v = val(tag)
+        try:
+            return int(float(v))
+        except (TypeError, ValueError):
+            return None
+
+    r = {
+        "sg": inum("_symmetry.Int_Tables_number"),
+        "uniq": inum("_reflns.number_obs"),
+        "obs": inum("_reflns.pdbx_number_measured_all"),
+        "rmeas": fnum("_reflns.pdbx_Rrim_I_all", 100.0),
+        "cc": fnum("_reflns.pdbx_CC_half", 100.0),
+        "isa": fnum("_reflns.jfjoch_diffrn_ISa"),
+    }
+
+    shells = []
+    for row in block.find("_reflns_shell.", ["d_res_high", "pdbx_Rrim_I_all", "pdbx_CC_half"]):
+        try:
+            shells.append((float(row[0]), float(row[1]), float(row[2])))
+        except ValueError:
+            pass
+    if shells:
+        shells.sort(key=lambda s: s[0])
+        r["cc_hi"] = shells[0][2] * 100.0      # smallest d_res_high = highest resolution
+        r["rmeas_lo"] = shells[-1][1] * 100.0  # largest  d_res_high = lowest resolution
+    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
+    if shutil.which("rugnux"):
+        return shutil.which("rugnux")
+    repo = pathlib.Path(_SELF).parent
+    for rel in ("build/tools/rugnux", "build_gpu/tools/rugnux",
+                "cmake-build-release-nogpu/tools/rugnux"):
+        if (repo / rel).exists():
+            return str(repo / rel)
+    return None
+
+
+def run_rugnux(master, workdir, name, xds, rugnux_bin, threads, timeout, reuse):
+    """Run rugnux; return (cif_path, error). Resolution + anomalous matched to XDS."""
+    workdir.mkdir(parents=True, exist_ok=True)
+    cif = workdir / f"{name}.cif"
+    if reuse and cif.exists():
+        return cif, None, None   # no fresh timing when we skip the run
+
+    cmd = [rugnux_bin, "-o", name, "--scaling-output", "cif"]
+    if xds.get("anomalous"):
+        cmd.append("-A")
+    if xds.get("dmin"):
+        cmd += ["--scaling-high-resolution", f"{xds['dmin']:.3f}"]
+    if threads:
+        cmd += ["-N", str(threads)]
+    cmd.append(str(master))
+
+    log = workdir / "rugnux.log"
+    t0 = time.monotonic()
+    with open(log, "w") as lf:
+        lf.write("$ " + " ".join(cmd) + "\n\n")
+        lf.flush()
+        try:
+            p = subprocess.run(cmd, cwd=str(workdir), stdout=lf,
+                               stderr=subprocess.STDOUT, timeout=timeout)
+        except subprocess.TimeoutExpired:
+            return None, f"timeout after {timeout}s", float(timeout)
+    dt = time.monotonic() - t0
+    if p.returncode != 0:
+        return None, f"exit {p.returncode} (see {log})", dt
+    if cif.exists():
+        return cif, None, dt
+    hits = sorted(workdir.glob("*.cif"))
+    return (hits[-1], None, dt) if hits else (None, f"no .cif produced (see {log})", dt)
+
+
+# --------------------------------------------------------------------------- #
+#  rendering
+# --------------------------------------------------------------------------- #
+def sg_name(n):
+    sg = gemmi.find_spacegroup_by_number(n) if n else None
+    if not sg:
+        return "?"
+    try:
+        return sg.short_name()
+    except Exception:
+        return sg.hm.replace(" ", "")
+
+
+def sym_key(n):
+    """Merge-relevant symmetry key = point group + lattice centring. Treats screw-axis /
+    enantiomorph variants as equal (XDS reports e.g. P422=89 where rugnux finds P4_3 2_1 2=96;
+    I23=197 vs I2_1 3=199) while keeping genuinely different lattices (P2 vs C2) distinct."""
+    sg = gemmi.find_spacegroup_by_number(n) if n else None
+    return f"{sg.point_group_hm()}{sg.hm[0]}" if sg else None
+
+
+W = dict(src=6, refl=8, obs=9, sg=14, rmeas=7, cc=6, rlo=8, chi=7, isa=6, time=8)
+
+
+def fmt_dur(s):
+    """Compact wall-clock: 42s / 2m05s / 1h03m."""
+    if s is None:
+        return None
+    s = int(round(s))
+    if s < 60:
+        return f"{s}s"
+    if s < 3600:
+        return f"{s // 60}m{s % 60:02d}s"
+    return f"{s // 3600}h{(s % 3600) // 60:02d}m"
+
+
+def _int(v, w):
+    return "-".rjust(w) if v is None else f"{v:>{w}d}"
+
+
+def _pct(v, w):
+    return "-".rjust(w) if v is None else f"{v:.1f}%".rjust(w)
+
+
+def _isa(v, w):
+    return "-".rjust(w) if v is None else f"{v:.2f}".rjust(w)
+
+
+def _sg(n, w):
+    return (f"{sg_name(n)} ({n})" if n else "-").ljust(w)
+
+
+def _dur(v, w):
+    return (fmt_dur(v) or "-").rjust(w)
+
+
+def row_line(src, d, elapsed=None):
+    d = d or {}
+    return (f"  {src:<{W['src']}} "
+            f"{_int(d.get('uniq'), W['refl'])} "
+            f"{_int(d.get('obs'), W['obs'])}  "
+            f"{_sg(d.get('sg'), W['sg'])} "
+            f"{_pct(d.get('rmeas'), W['rmeas'])} "
+            f"{_pct(d.get('cc'), W['cc'])} "
+            f"{_pct(d.get('rmeas_lo'), W['rlo'])} "
+            f"{_pct(d.get('cc_hi'), W['chi'])} "
+            f"{_isa(d.get('isa'), W['isa'])} "
+            f"{_dur(elapsed, W['time'])}")
+
+
+def header_line():
+    return (f"  {'Source':<{W['src']}} {'Refl':>{W['refl']}} {'Obs':>{W['obs']}}  "
+            f"{'Space group':<{W['sg']}} {'Rmeas':>{W['rmeas']}} {'CC1/2':>{W['cc']}} "
+            f"{'Rmeas_lo':>{W['rlo']}} {'CC1/2_hi':>{W['chi']}} {'ISa':>{W['isa']}} "
+            f"{'Time':>{W['time']}}")
+
+
+# --------------------------------------------------------------------------- #
+#  main
+# --------------------------------------------------------------------------- #
+def main():
+    ap = argparse.ArgumentParser(description="Compare rugnux merge stats against XDS CORRECT.LP.")
+    ap.add_argument("data_root", nargs="?", default="/data/rotation_test",
+                    help="directory of <crystal>/ subdirs (default: /data/rotation_test)")
+    ap.add_argument("--rugnux", help="path to the rugnux binary (else $RUGNUX / PATH / build*/tools)")
+    ap.add_argument("--workdir", default="rugnux_cmp", help="where rugnux outputs go (default: ./rugnux_cmp)")
+    ap.add_argument("--only", help="comma-separated substrings; only matching crystals are run")
+    ap.add_argument("--threads", type=int, help="rugnux -N (default: all hardware threads)")
+    ap.add_argument("--timeout", type=int, default=3600, help="per-crystal timeout, seconds (default: 3600)")
+    ap.add_argument("--reuse", action="store_true", help="skip rugnux where its .cif already exists")
+    ap.add_argument("--xds-only", action="store_true", help="only parse+print the XDS side (no rugnux run)")
+    ap.add_argument("--progress", action="store_true", help="print per-crystal progress to stderr")
+    args = ap.parse_args()
+
+    root = pathlib.Path(args.data_root)
+    if not root.is_dir():
+        sys.exit(f"data root not found: {root}")
+
+    rugnux_bin = None if args.xds_only else find_rugnux(args.rugnux)
+    if not args.xds_only and not rugnux_bin:
+        sys.exit("rugnux binary not found -- pass --rugnux PATH or set $RUGNUX")
+
+    filters = [s for s in (args.only or "").split(",") if s]
+    workdir = pathlib.Path(args.workdir)
+
+    crystals, skipped = [], []
+    for name in sorted(os.listdir(root)):
+        d = root / name
+        if not d.is_dir():
+            continue
+        if filters and not any(f in name for f in filters):
+            continue
+        correct = d / "CORRECT.LP"
+        master = sorted(d.glob("*_master.h5"))
+        if not correct.exists() or not master:
+            skipped.append(name)
+            continue
+        crystals.append((name, master[0], correct))
+
+    results = []
+    for i, (name, master, correct) in enumerate(crystals, 1):
+        xds = parse_xds(correct)
+        rug, err, elapsed = ({}, None, None)
+        if not args.xds_only:
+            if args.progress:
+                print(f"[{i}/{len(crystals)}] rugnux {name} "
+                      f"(dmin={xds.get('dmin')}, anom={xds.get('anomalous')}) ...",
+                      file=sys.stderr, flush=True)
+            cif, err, elapsed = run_rugnux(master, workdir / name, name, xds,
+                                           rugnux_bin, args.threads, args.timeout, args.reuse)
+            if cif:
+                try:
+                    rug = parse_rugnux_cif(cif)
+                except Exception as e:  # keep the table alive if one cif is unreadable
+                    err = f"cif parse failed: {e}"
+        results.append((name, xds, rug, err, elapsed))
+
+    # -------- print the table --------
+    hdr = header_line()
+    rule = "  " + "-" * (len(hdr) - 2)
+    today = datetime.date.today().isoformat()
+    print()
+    print(f"  rugnux vs XDS   ·   {root}   ·   {today}")
+    if not args.xds_only:
+        print(f"  rugnux: {rugnux_bin}")
+    print(f"  rugnux run de-novo, resolution + Friedel matched to XDS.  XDS SG = symmetry")
+    print(f"  CORRECT merged in; match is point-group level (screw/enantiomorph ignored).  [ - = missing ]")
+    print()
+    print(hdr)
+    print(rule)
+
+    n_match = 0
+    for name, xds, rug, err, elapsed in results:
+        if args.xds_only:
+            mark = ""
+        elif err:
+            mark = "  FAIL"
+        else:
+            xk, rk = sym_key(xds.get("sg")), sym_key(rug.get("sg"))
+            same = xk is not None and xk == rk
+            n_match += same
+            mark = "  SG OK  ✅" if same else "  SG DIFF ❌"
+        print(f"● {name}{mark}")
+        print(row_line("XDS", xds))
+        if args.xds_only:
+            continue
+        if err:
+            print(f"  {'rugnux':<{W['src']}} (failed: {err})  [{fmt_dur(elapsed) or '-'}]")
+        else:
+            print(row_line("rugnux", rug, elapsed))
+    print(rule)
+
+    if not args.xds_only:
+        n_fail = sum(1 for *_, e, _ in results if e)
+        total = sum(e for *_, e in results if e)
+        print(f"  {len(results)} crystals   ·   space group OK {n_match}/{len(results)}"
+              f"   ·   rugnux failed {n_fail}   ·   total rugnux time {fmt_dur(total) or '-'}")
+    if skipped:
+        print(f"  skipped (no CORRECT.LP or no *_master.h5): {', '.join(skipped)}")
+    print()
+
+
+if __name__ == "__main__":
+    main()