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()); +}