diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d3dbd56e..1c2cd348 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,7 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`) instead of a fixed 100, derived per crystal offline from the refined cell (`--max-hkl` overrides); the broker keeps a fixed bootstrap so a live acquisition has a predictable per-image cost. * rugnux: De-novo **space-group search** substantially more robust - centering ranked by net absences and judged on absent-class strength, merohedral-twin over-promotion vetoed, and genuine high-symmetry groups recovered on weak data. * rugnux: The space-group search takes systematic absences from the merge of all observations, needs at least three control reflections on an axial row to claim a **screw axis**, and no longer alters the production merge. -* rugnux: The per-image mosaicity is fitted from the strongest 250 spots, so the indexing spot budget no longer sets it. +* rugnux: The per-image mosaicity and profile radius are fitted from the strongest 250 spots, so the indexing spot budget no longer sets them. * rugnux: Stills **partiality post-refinement** added, on by default (`--simple-stills` disables); several non-helping stills scaling/detection knobs removed. * rugnux: **Scaling** hardened against a collapsed per-frame scale on the stills path as well as rotation, and the merged-sigma systematic floor is kept when ISa is too degenerate to report. * rugnux: `--min-image-cc` now works for rotation data (opt-in); new `--search-min-zeta` (rotation default 0.85); reports how close a symmetry axis lies to the spindle. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 10b48ed1..d74ebd40 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -731,7 +731,7 @@ A reference dataset (`--reference-mtz`) supplies known intensities for the same ### 11.1 Profile radius (intrinsic excitation-error width) -The “profile radius” is the intrinsic angular width of a reflection — crystal mosaicity plus beam divergence — estimated from the spread of $\Delta_\mathrm{Ewald}$ over indexed spots, +The “profile radius” is the intrinsic angular width of a reflection — crystal mosaicity plus beam divergence — estimated from the spread of $\Delta_\mathrm{Ewald}$ over the **strongest 250 indexed spots** (weaker spots sit further off the Ewald sphere, so an uncapped RMS grows with the indexing spot budget — measured ~60% over a 150-to-unlimited budget — and describes the spot list rather than the crystal), $ R \approx \sqrt{\tfrac{1}{N}\sum_i \Delta_{\mathrm{Ewald},i}^2}. $ diff --git a/image_analysis/indexing/FitProfileRadius.cpp b/image_analysis/indexing/FitProfileRadius.cpp index 43ef5f2f..322a9890 100644 --- a/image_analysis/indexing/FitProfileRadius.cpp +++ b/image_analysis/indexing/FitProfileRadius.cpp @@ -11,7 +11,20 @@ std::optional FitProfileRadius(const std::vector& spots, double sum_bw_var = 0.0; // energy-bandwidth contribution to subtract out int count = 0; - for (const auto &s: spots) { + // Fitted from the strongest spots only. This is an RMS over whatever spots were kept, and weaker + // spots sit further off the Ewald sphere, so it grows with the depth of the list: measured over a + // 150 -> unlimited spot budget it rises ~60% on both a hard and a clean rotation dataset. That + // made it a function of the INDEXING budget (--max-spots) rather than of the crystal. Its one + // consumer treats it as a membership gate (ewald_dist_cutoff = 2x it), where reflections at the + // cutoff carry ~zero partiality, so the integrated data barely moved (0.2% of partials across + // that range) - but it is also reported per image as a diagnostic, where a number that slides + // with an unrelated setting is simply misleading. FilterSpotsByCount leaves the list + // strongest-first, so taking the head selects the spots a smaller --max-spots would. + constexpr size_t PROFILE_RADIUS_FIT_SPOTS = 250; + const size_t n_fit = std::min(spots.size(), PROFILE_RADIUS_FIT_SPOTS); + + for (size_t si = 0; si < n_fit; ++si) { + const auto &s = spots[si]; if (!s.indexed) continue; sum_squares += static_cast(s.dist_ewald_sphere) * s.dist_ewald_sphere; diff --git a/rugnux_vs_xds.py b/rugnux_vs_xds.py index d56e2dd1..3423bf16 100755 --- a/rugnux_vs_xds.py +++ b/rugnux_vs_xds.py @@ -111,12 +111,60 @@ def parse_xds(correct_lp): # (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 + + # Mosaicity for comparison only. XDS refines this per run and the last value is the one it + # settled on; it is a useful sanity check on rugnux's own estimate but NOT ground truth -- + # XDS produces outliers here too, so read a big disagreement as "look at this crystal", + # not as "rugnux is wrong". + mos = re.findall(r"REFLECTING_RANGE_E\.S\.D\.=\s*([0-9.]+)", txt) + r["mosaicity"] = float(mos[-1]) if mos else None return r +def xds_integrate_mosaicity(integrate_lp): + """Mean per-batch SIGMAR from INTEGRATE.LP - XDS's INTEGRATION-STAGE estimate. + + This, not CORRECT.LP's REFLECTING_RANGE_E.S.D., is the apples-to-apples comparison for + rugnux's per-image mosaicity: both are fitted from the images during integration. The + CORRECT.LP value is XDS's POST-REFINED mosaicity, and the two can differ by more than 2x. + """ + if not integrate_lp.exists(): + return None + vals = [] + for line in integrate_lp.read_text(errors="replace").splitlines(): + if line.strip().startswith("SIGMAR (degree)"): + for f in line.split()[2:]: + try: + vals.append(float(f)) + except ValueError: + pass + return sum(vals) / len(vals) if vals else None + + # --------------------------------------------------------------------------- # # rugnux result (mmCIF) # --------------------------------------------------------------------------- # +def rugnux_mosaicity(workdir, name): + """Median per-image mosaicity (deg) from rugnux's _image.dat, or None.""" + import statistics + hits = sorted(workdir.glob(f"{name}*_image.dat")) + if not hits: + return None + vals = [] + for line in hits[-1].read_text(errors="replace").splitlines(): + if line.startswith("#"): + continue + f = line.split() + if len(f) >= 3: + try: + v = float(f[2]) + except ValueError: + continue + if v == v and v > 0: # skip nan / unfitted frames + vals.append(v) + return statistics.median(vals) if vals else None + + 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() @@ -242,7 +290,7 @@ def sym_key(n): 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) +W = dict(src=6, refl=8, obs=9, sg=14, rmeas=7, cc=6, rlo=8, chi=7, isa=6, mos=15, time=8) def fmt_dur(s): @@ -273,6 +321,13 @@ def _sg(n, w): return (f"{sg_name(n)} ({n})" if n else "-").ljust(w) +def _mos(v, w, v2=None): + a = "-" if v is None else f"{v:.4f}" + if v2 is not None: + a += f"|{v2:.4f}" + return a.rjust(w) + + def _dur(v, w): return (fmt_dur(v) or "-").rjust(w) @@ -288,6 +343,7 @@ def row_line(src, d, elapsed=None): f"{_pct(d.get('rmeas_lo'), W['rlo'])} " f"{_pct(d.get('cc_hi'), W['chi'])} " f"{_isa(d.get('isa'), W['isa'])} " + f"{_mos(d.get('mosaicity'), W['mos'], d.get('mosaicity_mle'))} " f"{_dur(elapsed, W['time'])}") @@ -295,7 +351,7 @@ 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']}}") + f"{'Mosaic':>{W['mos']}} {'Time':>{W['time']}}") # --------------------------------------------------------------------------- # @@ -344,6 +400,7 @@ def main(): results = [] for i, (name, master, correct) in enumerate(crystals, 1): xds = parse_xds(correct) + xds["mosaicity_mle"] = xds_integrate_mosaicity(correct.parent / "INTEGRATE.LP") rug, err, elapsed = ({}, None, None) if not args.xds_only: if args.progress: @@ -356,6 +413,7 @@ def main(): if cif: try: rug = parse_rugnux_cif(cif) + rug["mosaicity"] = rugnux_mosaicity(workdir / name, name) 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)) @@ -370,6 +428,10 @@ def main(): 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(f" Mosaic (deg): XDS cell = post-refined (CORRECT.LP) | integration-stage MLE (INTEGRATE.LP);") + print(f" rugnux cell = median per-image mosaicity. Compare rugnux against the INTEGRATION-STAGE") + print(f" number - the post-refined one measures something else and can differ by >2x.") + print(f" Reference only: XDS produces outliers here too, badly so on weak / low-resolution data.") print() print(hdr) print(rule)