From 6a1d859d499356a6c287ef46da98faafe5e4ddfd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 20 Oct 2026 14:52:07 +0200 Subject: [PATCH] Stills merge: R_split and a frame-split CC1/2, per shell and overall Serial crystallography quotes R_split and a CC1/2 whose half-sets keep each crystal whole (CrystFEL's partialator), and rugnux's stills report had neither, so its numbers could not be set beside CrystFEL's. MergeOnTheFly now also sums every frame's observations into the frame's own half (the existing per-image hash preference), and MergeStats reports from those two means R_split = (1/sqrt 2) sum|I1-I2| / (1/2 sum(I1+I2)) and their CC1/2. They are printed as a "Serial statistics" table under the shell table (also in the --report-resolution table) and as R_SPLIT / CC_HALF_SPLIT and their REFRES_ keys. The balanced per-reflection split behind CC_HALF is unchanged, and so is everything decided on it (resolution cut, space-group search, pass gates). Rotation output is unchanged: RotationScaleMerge does not fill the new fields and the table is not printed. On a 4000-frame jet subset over 30-1.5 A: CC1/2 0.830, R_split 31.8%, against 0.834 / 31.3% computed by hand from two separate even/odd runs and 0.842 / 30.5% from CrystFEL on the same frames. Co-Authored-By: Claude Opus 5.5 --- common/Reflection.h | 3 ++ docs/RUGNUX_REPORT.md | 9 +++++ image_analysis/scale_merge/Merge.cpp | 54 ++++++++++++++++++++++++++++ image_analysis/scale_merge/Merge.h | 16 +++++++++ rugnux/ResultReport.cpp | 9 +++++ tests/MergeScaleTest.cpp | 51 ++++++++++++++++++++++++++ 6 files changed, 142 insertions(+) diff --git a/common/Reflection.h b/common/Reflection.h index a6b5ad9c5..1b8b5caec 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -85,6 +85,9 @@ struct MergedReflection { float sigma = NAN; float I_half[2] = {NAN, NAN}; float sigma_half[2] = {NAN, NAN}; + // Half-set means with every frame's observations kept whole in one half (stills: R_split and the + // CrystFEL-style CC1/2). NaN where a half saw none, and on rotation. + float I_split[2] = {NAN, NAN}; // Weight of this reflection in a CC1/2: the precision its half-sets would have had with every // observation at the run's typical frame scale, over the precision they have. 1 on a sweep // without a weak stretch; small for a reflection measured only where the crystal barely diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 49736a946..bf18df563 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -67,6 +67,15 @@ the number with its noise; a large gap between the two says weak frames are kept absent where CC1/2 is not positive. All three are also columns of the shell tables (`R_merge`, `R_pim`, `CC*`) and are written per shell into the merged mmCIF. +**`R_SPLIT=`** and **`CC_HALF_SPLIT=`** (and their `REFRES_` forms) are written for stills only, with +a "Serial statistics" table beside the shell table. Every frame's observations go whole into one of two +half-sets - as CrystFEL's partialator deals out crystals - and R_split = (1/sqrt 2) sum|I_1 - I_2| / +(1/2 sum (I_1 + I_2)) (White et al., 2012) and CC1/2 are taken over the reflections both halves +measured, so both compare directly with CrystFEL's `compare_hkl`. `CC_HALF=` deals each reflection's +observations out evenly instead, as XDS and cctbx do, and is the one the processing decides on; the two +agree overall and can differ in a single shell. On stills `R_MEAS=` measures the scatter of partial +observations, not the error of the merged intensity, and is not the number to quote. + **`OBSERVATIONS_REJECTED=`** counts the observations the merge's outlier tests removed, and **`OBSERVATIONS_REJECTED_WILSON=`** the part of them the Wilson test took (each listed in the developer report). **`OBSERVATIONS_REJECTED_OVERLOAD=`** is separate from both: the rotation measurements left diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 039924825..e72daef80 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -150,6 +150,8 @@ void MergeOnTheFly::Accumulate(const std::vector &contributions, i it->second.sum_wI_half[half] += c.wI; it->second.sum_w_half[half] += c.w; it->second.n_half[half]++; + it->second.sum_wI_split[prefer] += c.wI; + it->second.sum_w_split[prefer] += c.w; if (c.on_ice_ring) it->second.on_ice_ring = true; @@ -452,6 +454,10 @@ std::vector MergeOnTheFly::ExportIntensities() { } } + for (int i = 0; i < 2; ++i) + if (accum.sum_w_split[i] > 0.0) + mr.I_split[i] = static_cast(accum.sum_wI_split[i] / accum.sum_w_split[i]); + if (!std::isfinite(accum.d) || accum.d <= 0.0f) continue; @@ -488,6 +494,11 @@ struct ShellAccum { CorrelationCoefficient cc_half; CorrelationCoefficient cc_ref; + + CorrelationCoefficient cc_split; + double split_abs_diff = 0.0; + double split_sum = 0.0; + int split_pairs = 0; }; std::pair ImageReferenceCC(std::span reflections, @@ -663,6 +674,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m CorrelationCoefficient cc_half_overall; CorrelationCoefficient cc_ref_overall; + ShellAccum split_overall; for (const auto &m: merged) { const auto shell = shells.GetShell(m.d); @@ -690,6 +702,15 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m cc_half_overall.Add(m.I_half[0], m.I_half[1]); } + if (std::isfinite(m.I_split[0]) && std::isfinite(m.I_split[1])) { + for (auto *a: {&acc[s], &split_overall}) { + a->cc_split.Add(m.I_split[0], m.I_split[1]); + a->split_abs_diff += std::fabs(m.I_split[0] - m.I_split[1]); + a->split_sum += 0.5 * (m.I_split[0] + m.I_split[1]); + a->split_pairs++; + } + } + } } } @@ -773,6 +794,16 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m rmeas_den_all += ra.sum_I; } + // R_split as White et al. (2012) J. Appl. Cryst. 45, 335-341 define it (CrystFEL's compare_hkl). + auto set_split = [](const ShellAccum &a, MergeStatisticsShell &ss) { + ss.n_split_pairs = a.split_pairs; + if (a.split_pairs == 0) + return; + ss.cc_half_split = a.cc_split.GetCC(); + if (a.split_sum > 0.0) + ss.r_split = a.split_abs_diff / a.split_sum / std::sqrt(2.0); + }; + MergeStatistics out; out.shells.resize(n_shells); @@ -795,6 +826,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; ss.r_merge = rmeas_den[s] > 0.0 ? rmerge_num[s] / rmeas_den[s] : NAN; ss.r_pim = rmeas_den[s] > 0.0 ? rpim_num[s] / rmeas_den[s] : NAN; + set_split(sa, ss); } auto &overall = out.overall; @@ -824,6 +856,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; overall.r_merge = rmeas_den_all > 0.0 ? rmerge_num_all / rmeas_den_all : NAN; overall.r_pim = rmeas_den_all > 0.0 ? rpim_num_all / rmeas_den_all : NAN; + set_split(split_overall, overall); return out; } @@ -923,6 +956,26 @@ std::string ShellTableRule() { "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""); } +std::string SerialShellTable(const MergeStatistics &in) { + if (in.overall.n_split_pairs == 0) + return ""; + std::string o = "\n Serial statistics: every frame whole in one half-set (as CrystFEL's partialator splits by crystal)\n"; + o += fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}\n", "d_min", "N_pairs", "CC1/2", "CC*", "R_split"); + o += fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}\n", "", "", "", "", ""); + auto row = [](const MergeStatisticsShell &sh) { + return fmt::format("{:8d} {:>8s} {:>8s} {:>8s}\n", sh.n_split_pairs, + Measured(sh.cc_half_split * 100.0, "%"), + Measured(CCStar(sh.cc_half_split) * 100.0, "%"), + Measured(sh.r_split * 100.0, "%")); + }; + for (const auto &sh: in.shells) + if (sh.n_split_pairs > 0) + o += fmt::format(" {:8.2f} ", sh.d_min) + row(sh); + o += fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}\n", "", "", "", "", ""); + o += fmt::format(" {:>8s} ", "Overall") + row(in.overall); + return o; +} + std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << std::endl; output << ShellTableHeader() << std::endl; @@ -939,6 +992,7 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << fmt::format(" {:>8s} ", "Overall"); output << in.overall; output << std::endl; + output << SerialShellTable(in); if (std::isfinite(in.wilson_b) && in.wilson_b > 0.0) output << fmt::format(" Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, fitted over {:.2f}-{:.2f} A)", in.wilson_b, in.wilson_b_correlation, in.wilson_b_d_low_A, in.wilson_b_d_high_A) diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index cfb5cd615..6d73040f0 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -58,6 +58,16 @@ struct MergeStatisticsShell { // it. NaN when no reflection has an anomalous split (the stills path), which is not the same // statement as zero - zero means measured and absent. double cc_anom = NAN; + + // The serial-crystallography pair (stills only, NaN on rotation): every frame's observations go + // whole to one half-set, as CrystFEL's partialator deals out crystals, and the two half-set means + // give CC1/2 and R_split = (1/sqrt 2) sum|I_1 - I_2| / (1/2 sum (I_1 + I_2)) (White et al. (2012) + // J. Appl. Cryst. 45, 335-341), over the reflections both halves measured (n_split_pairs). This + // CC1/2 is not cc_half above, which deals each reflection's observations out evenly and is the + // one the processing decides on. + double r_split = NAN; + double cc_half_split = NAN; + int n_split_pairs = 0; }; // Why a stretch of the sweep came out much weaker than the rest of the run. The codes are the field's @@ -227,6 +237,9 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in); // elsewhere (the report's reference-range table). std::string ShellTableHeader(); std::string ShellTableRule(); +// The serial table of R_split and the frame-split CC1/2, per shell and overall; empty where they were +// not computed (rotation). +std::string SerialShellTable(const MergeStatistics &in); struct MergeAccum { int32_t h = 0; @@ -239,6 +252,9 @@ struct MergeAccum { double sum_wI_half[2] = {0.0, 0.0}; double sum_w_half[2] = {0.0, 0.0}; size_t n_half[2] = {0, 0}; + // The same sums with each frame's observations kept whole in the frame's half (R_split). + double sum_wI_split[2] = {0.0, 0.0}; + double sum_w_split[2] = {0.0, 0.0}; }; // XDS's error-model convention (Diederichs, Acta Cryst. D66 (2010) 733) is sigma^2 = a*(sigma0^2 + diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index e1fc9a660..bc7bb3c64 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -330,6 +330,7 @@ namespace { } o << ShellTableRule() << "\n"; o << fmt::format(" {:>8s} ", "Overall") << rs.overall << "\n"; + o << SerialShellTable(rs); return o.str(); } @@ -1155,6 +1156,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("CC_HALF", o.cc_half, "{:.4f}")); if (std::isfinite(CCStar(o.cc_half))) Add(s, KeyReal("CC_STAR", CCStar(o.cc_half), "{:.4f}")); + if (std::isfinite(o.r_split)) + Add(s, KeyReal("R_SPLIT", o.r_split, "{:.4f}")); + if (std::isfinite(o.cc_half_split)) + Add(s, KeyReal("CC_HALF_SPLIT", o.cc_half_split, "{:.4f}")); if (std::isfinite(o.abs_diff_over_sigma_anomalous)) Add(s, KeyReal("SIGANO", o.abs_diff_over_sigma_anomalous, "{:.3f}")); if (std::isfinite(o.cc_anom)) @@ -1301,6 +1306,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("REFRES_R_MEAS_WEIGHTED", ro.r_meas_weighted, "{:.4f}")); if (std::isfinite(ro.cc_half)) Add(s, KeyReal("REFRES_CC_HALF", ro.cc_half, "{:.4f}")); + if (std::isfinite(ro.r_split)) + Add(s, KeyReal("REFRES_R_SPLIT", ro.r_split, "{:.4f}")); + if (std::isfinite(ro.cc_half_split)) + Add(s, KeyReal("REFRES_CC_HALF_SPLIT", ro.cc_half_split, "{:.4f}")); if (std::isfinite(ro.abs_diff_over_sigma_anomalous)) Add(s, KeyReal("REFRES_SIGANO", ro.abs_diff_over_sigma_anomalous, "{:.3f}")); if (std::isfinite(ro.cc_anom)) diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 6bc6e0ff2..70f08b06b 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -863,3 +863,54 @@ TEST_CASE("SublatticeRelabel_SignalPicksTheAxisThatWasDoubled") { CHECK((ScatteringVector(outcomes[0].latt, r) - ScatteringVector(before.latt, b)).Length() < 1e-5f); } } + +// R_split's half-sets keep every frame whole: each frame measures the same reflections at its own +// intensity, so every reflection's two half-set means are the same two numbers. R_split is then +// (1/sqrt 2) |I_1 - I_2| / ((I_1 + I_2) / 2) in every shell, and the frame-split CC1/2 is undefined +// (the half-set means do not vary) - the per-reflection balanced split, which the processing decides on, +// is left as it was. +TEST_CASE("MergeOnTheFly_RsplitKeepsEveryFrameWhole") { + DiffractionExperiment x; + x.SetSpaceGroup(TestSpaceGroup()); + + std::vector outcomes(64); + for (size_t f = 0; f < outcomes.size(); ++f) { + for (int h = 1; h <= 4; ++h) + for (int k = 1; k <= 4; ++k) { + Reflection r{}; + r.h = h; + r.k = k; + r.l = 1; + r.d = 20.0f / std::sqrt(static_cast(h * h + k * k + 1)); + r.I = 100.0f + static_cast(f); + r.sigma = 10.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + r.image_scale_corr = 1.0f; + outcomes[f].reflections.push_back(r); + } + } + + MergeOnTheFly merge(x); + merge.AddImages(outcomes); + const auto merged = merge.ExportReflections(); + REQUIRE(merged.size() == 16); + + const float a = merged[0].I_split[0], b = merged[0].I_split[1]; + REQUIRE(std::isfinite(a)); + REQUIRE(std::isfinite(b)); + CHECK(a != b); + for (const auto &m: merged) { + CHECK(m.I_split[0] == a); + CHECK(m.I_split[1] == b); + } + + const auto stats = merge.MergeStats(merged, outcomes); + const double expected = std::fabs(a - b) / (0.5 * (a + b)) / std::sqrt(2.0); + CHECK(stats.overall.n_split_pairs == 16); + CHECK(stats.overall.r_split == Catch::Approx(expected).epsilon(1e-5)); + for (const auto &sh: stats.shells) + if (sh.n_split_pairs > 0) + CHECK(sh.r_split == Catch::Approx(expected).epsilon(1e-5)); + CHECK_FALSE(SerialShellTable(stats).empty()); +}