From 8a0da57941696f5c2ef3c28c453c06a0be8f409d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 19 Sep 2026 10:56:26 +0200 Subject: [PATCH] rugnux: --report-resolution, a second statistics table at a reference range Comparing a run with another program's table has meant running rugnux AT that program's resolution range (--scaling-high-resolution), which is a different run: the range moves the cut, the space-group decision and everything after them, so the comparison buys itself a different answer. --report-resolution [,] instead leaves the run alone and adds a second table to section 3 of the report - the REFRES_* keys and a shell table - binned from the same merged reflections over the range given, with the completeness denominator enumerated over that range and the shells in equal steps of 1/d^2 so they read row for row against a CORRECT.LP at the same range. Report-only: the merged files and every decision are byte-identical with and without it. The table holds only what the run kept. Where the reference range is finer than the run's own limit, the shells past it are printed as not merged (with their possible count) rather than as zeros, REFRES_SHELLS_PAST_LIMIT counts them so a consumer can tell "not merged" from a measured zero, REFRES_ COMPLETENESS counts their reflections as missing, and the other overall numbers are over the shells the run reached; nothing is read from the observations the run judged to carry no signal. REFRES_ISA is the error model refitted on the reflections of the table alone, in XDS's convention (rotation only; the stills model is fitted over the whole range already). On the rotation path the statistics block of MergeAndStats becomes a lambda over a shell grid, called once for the run's own grid and once for the reference one; the reference call floors every observation-level count at the cut by group d, the rule the erase applied. The stills MergeStats takes a declared range, whose bounds are the grid's whether or not any reflection reaches them. Both --mode mx and --mode scale report it, the viewer's command line echoes it, and the docs describe the keys. Co-Authored-By: Claude Fable 5.1 --- common/ScalingSettings.cpp | 12 + common/ScalingSettings.h | 13 + docs/RUGNUX_ADVANCED.md | 1 + docs/RUGNUX_REPORT.md | 31 + image_analysis/scale_merge/Merge.cpp | 40 +- image_analysis/scale_merge/Merge.h | 10 +- .../scale_merge/RotationScaleMerge.cpp | 555 ++++++++++-------- .../scale_merge/RotationScaleMerge.h | 9 + rugnux/ResultReport.cpp | 88 +++ rugnux/Rugnux.cpp | 13 +- rugnux/Rugnux.h | 9 + rugnux/RugnuxCommandLine.cpp | 5 + rugnux/rugnux_cli.cpp | 46 ++ tests/MergeScaleTest.cpp | 41 ++ tests/ResultReportTest.cpp | 55 ++ 15 files changed, 661 insertions(+), 267 deletions(-) diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 773ed9fa6..2016afba5 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -303,3 +303,15 @@ ScalingSettings &ScalingSettings::ReportShellCount(int input) { int ScalingSettings::GetReportShellCount() const { return report_shell_count; } + +ScalingSettings& ScalingSettings::ReportResolutionRange_A(std::optional input) { + if (input && !(input->d_min > 0.0 && input->d_max > input->d_min)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Report resolution range must be positive and d_max coarser than d_min"); + report_resolution_range = input; + return *this; +} + +std::optional ScalingSettings::GetReportResolutionRange_A() const { + return report_resolution_range; +} diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index ac25888c9..1d2e5a9bb 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -11,6 +11,12 @@ // CCHalfLogistic = fit the CC1/2 fall-off and cut one shell past cc_target (DIALS-style, generous). enum class ResolutionCutoffMethod { Off, CCHalfLogistic }; +// A resolution range in A, d_max coarse and d_min fine; d_max may be infinite (no low-resolution bound). +struct ReportResolutionRange { + double d_min = 0.0; + double d_max = 0.0; +}; + class ScalingSettings { bool refine_wedge = false; @@ -109,6 +115,11 @@ class ScalingSettings { double resolution_cc_target = 0.30; // CC1/2 value defining the fall-off limit before the +1 shell int report_shell_count = 9; // resolution shells in the reported statistics table (XDS's count, // over the same equal-1/d^2 bins, so the two tables read row for row) + // A second statistics table over this REFERENCE range (--report-resolution), binned from the same + // scaled and merged observations as the run's own table and reported beside it, so a run can be + // read against another program's numbers at that program's range. Report-only: the cut, the + // scaling, the error model and every decision are the run's own whether or not it is set. + std::optional report_resolution_range; public: ScalingSettings& RefineRotationWedge(bool input); @@ -139,6 +150,7 @@ public: ScalingSettings& ResolutionCutoff(ResolutionCutoffMethod input); ScalingSettings& ResolutionCCTarget(double input); ScalingSettings& ReportShellCount(int input); + ScalingSettings& ReportResolutionRange_A(std::optional input); [[nodiscard]] bool GetRefineWedge() const; @@ -177,4 +189,5 @@ public: [[nodiscard]] ResolutionCutoffMethod GetResolutionCutoff() const; [[nodiscard]] double GetResolutionCCTarget() const; [[nodiscard]] int GetReportShellCount() const; + [[nodiscard]] std::optional GetReportResolutionRange_A() const; }; diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 6a8a885ed..7984cfd5e 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -323,6 +323,7 @@ Scaling and merging: | `--resolution-cutoff ` | Automatic high-resolution cutoff for the written reflections and reported shells: `cc-logistic` \| `off` (default: `cc-logistic`; ignored when `--scaling-high-resolution` is set) | | `--resolution-cc-target ` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) | | `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 9). The bins are equal steps in 1/d² between the lowest- and highest-resolution reflection merged, which is XDS's rule, and 9 is XDS's count — so at the same resolution limits the two tables have the same shells and can be read row for row | +| `--report-resolution [,]` | Also report the merging statistics over this resolution range, Å (either order; `dmax` defaults to the run's own low-resolution limit), as a second table beside the run's own — the `REFRES_*` keys of the [results report](RUGNUX_REPORT.md#the-reference-range-table) — binned from the same merged reflections, so a run can be read against another program's table at that program's range. Report-only: the cut, the scaling, the error model and every decision stay the run's own, and shells finer than the run's own limit are reported as not merged rather than filled from data the run did not keep | | `--min-partiality ` | Minimum partiality to accept a reflection (default: 0.02) | | `--ice-min-score ` | Ice-presence gate: the measured per-run ice score (1 = no ice) a dataset must reach before **any** ice handling is applied — the flagging and the exclusion from scaling (default: 1.5; 0 = no gate). The fixed hexagonal bands cover 16–26 % of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal only costs completeness | | `--ice-min-spot-ratio ` | The second ice-presence channel: found **spots** on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so `--ice-min-score` alone is blind to it (default: 2.0; 0 disables this channel) | diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index b7ed36671..264bf1cb9 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -148,6 +148,37 @@ and a run that produced a textbook data set read identically for their first thr Then the numbered sections: **1. DATA SET AND GEOMETRY**, **2. CRYSTAL**, **3. MERGED DATA**, **4. DIAGNOSTICS**, and **5. MODEL VALIDATION** only with `--model`. +## The reference-range table + +`--report-resolution [,]` adds a second block of merging statistics to section 3, the +`REFRES_*` keys and a second shell table, over the resolution range it is given rather than the range +the run chose. It exists for comparison: another program's table is at that program's range, and +running rugnux *at* that range (`--scaling-high-resolution`) is not the same run — the cut moves, and +with it the symmetry decision and everything downstream of it. The reference table is instead the +same merged reflections binned again, so nothing is processed differently whether or not it is +asked for; the merged files and every decision are byte-for-byte the run's own. + +- **`REFRES_RANGE=`** is the requested range (`dmax dmin`, as `INCLUDE_RESOLUTION_RANGE`; `INF` when + the run has no low-resolution limit and none was given). **`REFRES_MEASURED_RANGE=`** is the + coarsest and finest merged reflection that actually landed in it. +- **`REFRES_COMPLETENESS=`**, **`REFRES_MULTIPLICITY=`**, **`REFRES_I_OVER_SIGMA=`**, + **`REFRES_R_MEAS=`**, **`REFRES_CC_HALF=`**, **`REFRES_SIGANO=`**, **`REFRES_CC_ANOM=`**, + **`REFRES_UNIQUE_REFLECTIONS=`** and **`REFRES_TOTAL_OBSERVATIONS=`** are the overall numbers over + that range, under the same absent-when-unmeasured rule as their section-3 namesakes. + **`REFRES_ISA=`** is the error model refitted on the reflections of this table alone, in XDS's + convention, so it reads against an `ISa` produced at that range; the table itself is merged under + the run's own model. +- **The table holds only what the run kept.** Where the reference range is finer than the run's own + limit, the shells past that limit are empty by the run's decision — it judged them to carry no + signal and did not merge them — and are printed as `past the run's own limit of X A: not merged + (N possible)` rather than as zeros; the shell the limit falls inside is marked. **`REFRES_SHELLS_PAST_LIMIT=`** + counts those shells (0 when the range lies within the run's own), so a consumer can tell *not + merged* from a measured zero. `REFRES_COMPLETENESS` counts their reflections as missing; every + other overall number is over the shells the run reached. + +The shells are equal steps in 1/d² between the two bounds, as XDS's are, so at XDS's +`INCLUDE_RESOLUTION_RANGE` the two tables read row for row. + ## The developer report `--developer` renders the same report in full. The default report carries what a person deciding diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 2aeebdb20..ce7b87f2e 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -512,13 +512,19 @@ void CalcPossibleReflections(const gemmi::SpaceGroup &sg, MergeStatistics MergeOnTheFly::MergeStats(const std::vector &merged, const std::vector &integration_outcome, const std::vector &reference, - std::optional d_min_override) { + std::optional d_min_override, + std::optional declared_range) { const int n_shells = scaling_settings.GetReportShellCount(); auto d_min_limit_A = d_min_override.has_value() ? d_min_override : scaling_settings.GetHighResolutionLimit_A(); - const auto d_max_limit_A = scaling_settings.GetLowResolutionLimit_A(); + auto d_max_limit_A = scaling_settings.GetLowResolutionLimit_A(); + if (declared_range) { + d_min_limit_A = declared_range->d_min; + d_max_limit_A = std::isfinite(declared_range->d_max) + ? std::optional(declared_range->d_max) : std::nullopt; + } std::unordered_map reference_intensities; if (!reference.empty()) { @@ -547,7 +553,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m d_max = std::max(d_max, m.d); } - if (!(d_min < d_max && d_min > 0.0f)) + // A declared range is a grid whether or not anything landed in it. + if (!(d_min < d_max && d_min > 0.0f) && !declared_range) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "MergeStats: Error in resolution calculation"); @@ -560,7 +567,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m // infinite d_max is how that is spelled. The high end stays the finest d reached anywhere - // that is the resolution the run reports, and the full sphere down to it already shows // anisotropic loss. - const float d_min_pad = d_min * 0.999f; // GetShell rejects d <= d_min; the finest reflection needs the room + const float d_min_pad = declared_range ? static_cast(declared_range->d_min) + : d_min * 0.999f; // GetShell rejects d <= d_min; the finest reflection needs the room const float d_max_declared = d_max_limit_A ? static_cast(*d_max_limit_A) : std::numeric_limits::infinity(); // The innermost shell is labelled by the bound it reaches, except when that bound is infinite @@ -696,7 +704,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m } auto &overall = out.overall; - overall.d_min = d_min; + overall.d_min = d_min < d_max ? d_min : 0.0f; // 0: nothing landed in a declared range overall.d_max = d_max; int all_possible = 0; @@ -792,14 +800,21 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { return output; } +std::string ShellTableHeader() { + return fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref", + "CCanom", "SigAno"); +} + +std::string ShellTableRule() { + return fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", + "", "", "", "", "", "", "", "", "", "", "", ""); +} + std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << std::endl; - output << fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", - "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref", - "CCanom", "SigAno") - << std::endl; - output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "", "", "") << std::endl; + output << ShellTableHeader() << std::endl; + output << ShellTableRule() << std::endl; for (const auto &sh: in.shells) { if (sh.unique_reflections == 0) continue; @@ -807,8 +822,7 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << sh; output << std::endl; } - output << fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "", "", "") << std::endl; + output << ShellTableRule() << std::endl; output << fmt::format(" {:>8s} ", "Overall"); output << in.overall; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 32144fad0..025dfe7cb 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -165,6 +165,10 @@ struct MergeStatistics { std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in); std::ostream &operator<<(std::ostream &output, const MergeStatistics &in); +// The header and the rule of the table operator<< prints, for a table with the same columns rendered +// elsewhere (the report's reference-range table). +std::string ShellTableHeader(); +std::string ShellTableRule(); struct MergeAccum { int32_t h = 0; @@ -283,10 +287,14 @@ public: // d_min_override, when set, is the effective high-resolution limit for the shell table (used for // the automatic resolution cutoff computed by the caller); otherwise the manual // ScalingSettings high-resolution limit stands. The number of shells is ScalingSettings::ReportShellCount. + // declared_range: bin over exactly this range instead - its bounds are the shell grid's and the + // completeness denominator's whether or not any reflection reaches them, and the limits above + // are ignored. The reference-range table of --report-resolution. MergeStatistics MergeStats(const std::vector &merged, const std::vector &reflections, const std::vector &reference = {}, - std::optional d_min_override = std::nullopt); + std::optional d_min_override = std::nullopt, + std::optional declared_range = std::nullopt); std::vector ExportReflections(); }; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 8d4f2da49..b95fac06b 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -250,6 +250,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, resolution_cutoff_method = s.GetResolutionCutoff(); resolution_cc_target = s.GetResolutionCCTarget(); report_shell_count = s.GetReportShellCount(); + report_range = s.GetReportResolutionRange_A(); if (const auto forced = s.GetForcedMosaicity(); forced.has_value()) mosaicity_deg = *forced; else @@ -4080,8 +4081,274 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); - // ---- Statistics (report_shell_count shells): completeness, multiplicity, , R_meas, CC1/2. ---- + // ---- Statistics (report_shell_count shells): completeness, multiplicity, , R_meas, CC1/2, + // SigAno and CCanom, over a shell grid. A lambda because the same merged reflections are + // binned a second time, over the reference range, when --report-resolution asked for that + // table (see the end of this function). The grid runs from grid_d_min (exclusive) to + // grid_d_max, and the completeness denominator is enumerated over exactly that. The fulls + // hold every group that was ingested, cut or not, so a grid reaching past the cut needs + // d_floor - the cut, as the erase above applied it, by group d - to keep the observations + // of the groups the run did not keep out of every count; the run's own grid ends at the + // finest kept reflection and needs none. ---- const int n_shells = report_shell_count; + struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; }; + std::unordered_map anom_export; // filled by the run's own table only + auto statistics_over = [&](const std::vector &merged, float grid_d_min, + float grid_d_max, float d_max_label, float d_floor, + bool own_table) -> MergeStatistics { + ResolutionShells shells(grid_d_min, grid_d_max, n_shells); + const auto shell_mean_1_d2 = shells.GetShellMeanOneOverResSq(); + const auto shell_min_res = shells.GetShellMinRes(); + + struct ShellAcc { + int unique = 0, total_obs = 0, possible = 0; + double sum_i_over_sigma = 0.0; int n_i_over_sigma = 0; + CorrelationCoefficient cc_half; + }; + std::vector sa(n_shells); + std::vector possible(n_shells, 0); + // Completeness is a reported quantity: it is read off the merge a pass writes out, never off the + // ones the space-group search runs on the way there. Those are also the expensive ones to count - + // the search merges in P1, where the list of possible reflections is the whole hemisphere rather + // than an asymmetric unit of it - so enumerating it there is the largest single piece of work in + // the merge that nothing goes on to read. + if (reference_cell && !for_search) + PossiblePerShell(x.GetSpaceGroupOrP1(), *reference_cell, grid_d_min, grid_d_max, + shells, merge_friedel, possible, nthreads); + for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; + + CorrelationCoefficient cc_half_overall; + for (const auto &m : merged) { + const auto shell = shells.GetShell(m.d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0) { + auto &s = sa[*shell]; + s.unique++; + s.sum_i_over_sigma += m.I / m.sigma; s.n_i_over_sigma++; + if (std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) { + s.cc_half.Add(m.I_half[0], m.I_half[1]); + cc_half_overall.Add(m.I_half[0], m.I_half[1]); + } + } + } + + std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); + double rmeas_num_all = 0.0, rmeas_den_all = 0.0; + // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last + // full pass over the observations there is - so, like the correction surfaces, it is not made on a + // merge whose numbers nobody reads. + if (full_stats) { + // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate + // |I_i - | per reflection. + struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; + std::vector rmeas(n_groups); + bool did_gpu_rmeas = false; +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) { + // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). + std::vector rabsdev(n_groups), rsumI(n_groups); + std::vector rn(n_groups), rnusable(n_groups); + gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); + for (int g = 0; g < n_groups; ++g) { + if (rnusable[g] == 0 || acc[g].d < d_floor) continue; + const auto shell = shells.GetShell(acc[g].d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately + // ignores the ice flag on a search pass, so its count includes observations that never + // entered `unique` - which inflates the reported multiplicity of whatever shell they + // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a + // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above + // already declines it; this is the same statement made where it counts.) + sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); + if (std::isfinite(merged_I[g]) && rn[g] > 0) { + auto &r = rmeas[g]; + r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; + } + } + did_gpu_rmeas = true; + } +#endif + if (!did_gpu_rmeas) + for (const auto &o : fulls) { + if (o.group < 0 || acc[o.group].d < d_floor) continue; + if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas + if (!frame_cell_ok[o.frame]) continue; + if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; + if (o.partiality < min_partiality) continue; + const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; + if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; + const auto shell = shells.GetShell(o.d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Only what the merge kept counts towards multiplicity - see the GPU branch above. The + // R_meas accumulation below keeps its own, wider filter. + if (usable_merge(o)) + sa[*shell].total_obs++; + if (std::isfinite(merged_I[o.group])) { + auto &r = rmeas[o.group]; + r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); + r.sum_I += I_corr; r.n++; r.shell = *shell; + } + } + for (const auto &r : rmeas) { + if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; + const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); + rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; + rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; + } + } + + // ---- Anomalous split (always, even when the merge is Friedel-averaged): for each acentric + // reflection, the inverse-variance I(+)/I(-) from the SAME scaled fulls (so it never touches the + // Friedel-merged IMEAN, scaling or error model above). Lets I(+)/I(-) and F(+)/F(-) be written by + // default without scaling anomalously; a reflection with only one mate, or a centric, is left + // without the split. Skipped for the P1 search pass, which has no use for it, and for a merge that + // is not written. ---- + // SigAno = <|I(+)-I(-)|>/ accumulated on the standard report shells (num=sum|dI|, + // den=sum sigma(dI)); mmCIF pdbx_absDiff_over_sigma_anomalous. NaN'd below when there is no split. + std::vector sig_num(n_shells, 0.0), sig_den(n_shells, 0.0); + double sig_num_all = 0.0, sig_den_all = 0.0; + size_t sig_n = 0; + // CCanom on the same shells, and over the whole range as one correlation (not a mean of shells). + std::vector cc_anom_shell(n_shells); + CorrelationCoefficient cc_anom_overall; + size_t cc_anom_n = 0; + if (!for_search && full_stats) { + const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupOrP1()); + const gemmi::GroupOps gops = x.GetSpaceGroupOrP1().operations(); + // [hand][half]: hand 0 = I(+), 1 = I(-). The half sets are what CCanom needs - the anomalous + // difference formed twice, once from each half of the observations - and their sums add back + // to the whole-mate I(+)/I(-) that SigAno and the export read, so there is one accumulation. + struct AnomAcc { + double swI[2][2] = {}; double sw[2][2] = {}; int32_t nh[2][2] = {}; + int32_t h = 0, k = 0, l = 0; float d = NAN; + }; + // The Friedel-merged key an observation lands on, and which hand it is, are properties of its + // ASU GROUP, not of the observation - group_h/k/l is that group's SIGNED representative, so the + // same reduction applied to it returns both. Reduce once per group (a few hundred thousand) + // instead of once per observation (millions), and index a dense accumulator from there, which + // is the rule ComputeAsuGroups states for itself. The hand only follows the group when the + // merge distinguishes the hands; a Friedel-merged run holds both in one group and must ask. + std::vector pair_of_group(n_groups, -1); + std::vector hand_of_group(merge_friedel ? 0 : n_groups, 0); + std::vector pair_hkl; + { + std::unordered_map pair_id; + pair_id.reserve(static_cast(n_groups) + 1); + for (int g = 0; g < n_groups; ++g) { + const HKLKey gk = anom_keygen(group_h[g], group_k[g], group_l[g]); + if (!merge_friedel) hand_of_group[g] = gk.plus ? 0 : 1; + const uint64_t key = HKLKey{gk.h, gk.k, gk.l, true}.pack(); + const auto [it, fresh] = pair_id.emplace(key, static_cast(pair_hkl.size())); + if (fresh) pair_hkl.push_back(HKLKey{gk.h, gk.k, gk.l, true}); + pair_of_group[g] = it->second; + } + } + std::vector anom(pair_hkl.size()); + // The half-sets CCanom needs are balanced within the MATE, not merely within the reflection: + // it correlates one half's I(+)-I(-) against the other's, so each hand separately needs both + // of its halves filled, and a mate left entirely in one half loses the whole pair rather than + // half of it. So the cell here is (pair, hand), where the merge above used the reflection. + std::vector anom_cell(fulls.size(), -1); + for (size_t i = 0; i < fulls.size(); ++i) { + const Obs &o = fulls[i]; + if (!usable_merge(o) || rejected_obs[i]) continue; // as the merge saw it + if (acc[o.group].d < d_floor) continue; + const int hand = merge_friedel ? (anom_keygen(o.h, o.k, o.l).plus ? 0 : 1) + : hand_of_group[o.group]; + anom_cell[i] = pair_of_group[o.group] * 2 + hand; + } + std::vector anom_half; + AssignHalvesByRank(obs_key, anom_cell, static_cast(pair_hkl.size()) * 2, anom_half, nthreads); + for (size_t i = 0; i < fulls.size(); ++i) { + const int cell = anom_cell[i]; + if (cell < 0) continue; + const Obs &o = fulls[i]; + const float I_corr = o.I * o.corr; + const float sigma_corr = corrected_sigma(o, I_corr, o.sigma * o.corr); + if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue; + const double w = 1.0 / (static_cast(sigma_corr) * sigma_corr); + AnomAcc &a = anom[cell >> 1]; + a.d = o.d; + a.swI[cell & 1][anom_half[i]] += w * static_cast(I_corr); + a.sw[cell & 1][anom_half[i]] += w; + a.nh[cell & 1][anom_half[i]]++; + } + for (size_t p = 0; p < anom.size(); ++p) { + AnomAcc a = anom[p]; + a.h = pair_hkl[p].h; a.k = pair_hkl[p].k; a.l = pair_hkl[p].l; + const uint64_t fkey = pair_hkl[p].pack(); + // Centrics have I(+)=I(-) by symmetry; leave them without an anomalous split. + if (gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l})) continue; + // Whole-mate sums: the two halves put back together. I(+)/I(-) and SigAno are read off + // these, so splitting the accumulation above changed neither. + const double swp = a.sw[0][0] + a.sw[0][1], swm = a.sw[1][0] + a.sw[1][1]; + const double swIp = a.swI[0][0] + a.swI[0][1], swIm = a.swI[1][0] + a.swI[1][1]; + AnomExport ex; + if (swp > 0.0) { ex.Ip = static_cast(swIp / swp); ex.sIp = static_cast(1.0 / std::sqrt(swp)); } + if (swm > 0.0) { ex.Im = static_cast(swIm / swm); ex.sIm = static_cast(1.0 / std::sqrt(swm)); } + if (own_table && (std::isfinite(ex.Ip) || std::isfinite(ex.Im))) anom_export[fkey] = ex; + // A pair outside the grid is outside the table, its overall line included: the fulls + // hold every pair that was ingested, and the ones past the written range are not + // the data these numbers describe. + const auto shell = shells.GetShell(a.d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + if (swp > 0.0 && swm > 0.0) { + const double abs_di = std::fabs(swIp / swp - swIm / swm); + const double sig_di = std::sqrt(1.0 / swp + 1.0 / swm); + sig_num[*shell] += abs_di; sig_den[*shell] += sig_di; + sig_num_all += abs_di; sig_den_all += sig_di; ++sig_n; + } + // CCanom: the anomalous difference measured twice, once per half, correlated over the + // pairs where all four cells were filled - which needs both mates measured at least + // twice, so it asks more of the data than SigAno does. + if (a.sw[0][0] > 0.0 && a.sw[0][1] > 0.0 && a.sw[1][0] > 0.0 && a.sw[1][1] > 0.0) { + const double d_first = a.swI[0][0] / a.sw[0][0] - a.swI[1][0] / a.sw[1][0]; + const double d_second = a.swI[0][1] / a.sw[0][1] - a.swI[1][1] / a.sw[1][1]; + cc_anom_shell[*shell].Add(d_first, d_second); + cc_anom_overall.Add(d_first, d_second); + ++cc_anom_n; + } + } + } + + MergeStatistics out; + out.shells.resize(n_shells); + for (int s = 0; s < n_shells; ++s) { + auto &ss = out.shells[s]; + ss.mean_one_over_d2 = shell_mean_1_d2[s]; + ss.d_min = shell_min_res[s]; + ss.d_max = s == 0 ? d_max_label : shell_min_res[s - 1]; + ss.total_observations = sa[s].total_obs; + ss.unique_reflections = sa[s].unique; + ss.possible_unique_reflections = sa[s].possible; + ss.mean_i_over_sigma = sa[s].n_i_over_sigma > 0 ? sa[s].sum_i_over_sigma / sa[s].n_i_over_sigma : 0.0; + ss.cc_half = sa[s].cc_half.GetCC(); + ss.cc_ref = NAN; + ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; + ss.abs_diff_over_sigma_anomalous = sig_den[s] > 0.0 ? sig_num[s] / sig_den[s] : NAN; + ss.cc_anom = cc_anom_shell[s].GetCC(); + } + auto &overall = out.overall; + double sum_ios = 0.0; int n_ios = 0; + for (int s = 0; s < n_shells; ++s) { + overall.total_observations += sa[s].total_obs; + overall.unique_reflections += sa[s].unique; + overall.possible_unique_reflections += sa[s].possible; + sum_ios += sa[s].sum_i_over_sigma; n_ios += sa[s].n_i_over_sigma; + } + overall.mean_i_over_sigma = n_ios > 0 ? sum_ios / n_ios : 0.0; + overall.cc_half = cc_half_overall.GetCC(); + overall.cc_ref = NAN; + overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; + overall.abs_diff_over_sigma_anomalous = sig_den_all > 0.0 ? sig_num_all / sig_den_all : NAN; + overall.cc_anom = cc_anom_overall.GetCC(); + if (own_table && std::isfinite(overall.abs_diff_over_sigma_anomalous)) + logger.Info("Anomalous signal SigAno = {:.2f} ({} acentric pairs), CCanom = {:.3f} ({} pairs " + "split in both hands)", + overall.abs_diff_over_sigma_anomalous, sig_n, overall.cc_anom, cc_anom_n); + return out; + }; + float sd_min = std::numeric_limits::max(), sd_max = 0.0f; for (const auto &m : result.merged) { if (!std::isfinite(m.d) || m.d <= 0.0f) continue; @@ -4102,258 +4369,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // The innermost shell is labelled by the bound it reaches, except when that bound is infinite // and the coarsest reflection measured is the only finite thing to say. const float d_max_label = std::isfinite(d_max_declared) ? d_max_declared : sd_max; - ResolutionShells shells(d_min_pad, d_max_declared, n_shells); - const auto shell_mean_1_d2 = shells.GetShellMeanOneOverResSq(); - const auto shell_min_res = shells.GetShellMinRes(); - - struct ShellAcc { - int unique = 0, total_obs = 0, possible = 0; - double sum_i_over_sigma = 0.0; int n_i_over_sigma = 0; - CorrelationCoefficient cc_half; - }; - std::vector sa(n_shells); - std::vector possible(n_shells, 0); - // Completeness is a reported quantity: it is read off the merge a pass writes out, never off the - // ones the space-group search runs on the way there. Those are also the expensive ones to count - - // the search merges in P1, where the list of possible reflections is the whole hemisphere rather - // than an asymmetric unit of it - so enumerating it there is the largest single piece of work in - // the merge that nothing goes on to read. - if (reference_cell && !for_search) - PossiblePerShell(x.GetSpaceGroupOrP1(), *reference_cell, d_min_pad, d_max_declared, - shells, merge_friedel, possible, nthreads); - for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; - - CorrelationCoefficient cc_half_overall; - for (const auto &m : result.merged) { - const auto shell = shells.GetShell(m.d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0) { - auto &s = sa[*shell]; - s.unique++; - s.sum_i_over_sigma += m.I / m.sigma; s.n_i_over_sigma++; - if (std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) { - s.cc_half.Add(m.I_half[0], m.I_half[1]); - cc_half_overall.Add(m.I_half[0], m.I_half[1]); - } - } - } - - std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); - double rmeas_num_all = 0.0, rmeas_den_all = 0.0; - // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last - // full pass over the observations there is - so, like the correction surfaces, it is not made on a - // merge whose numbers nobody reads. - if (full_stats) { - // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate - // |I_i - | per reflection. - struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; - std::vector rmeas(n_groups); - bool did_gpu_rmeas = false; -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) { - // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). - std::vector rabsdev(n_groups), rsumI(n_groups); - std::vector rn(n_groups), rnusable(n_groups); - gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); - for (int g = 0; g < n_groups; ++g) { - if (rnusable[g] == 0) continue; - const auto shell = shells.GetShell(acc[g].d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately - // ignores the ice flag on a search pass, so its count includes observations that never - // entered `unique` - which inflates the reported multiplicity of whatever shell they - // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a - // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above - // already declines it; this is the same statement made where it counts.) - sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); - if (std::isfinite(merged_I[g]) && rn[g] > 0) { - auto &r = rmeas[g]; - r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; - } - } - did_gpu_rmeas = true; - } -#endif - if (!did_gpu_rmeas) - for (const auto &o : fulls) { - if (o.group < 0) continue; - if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas - if (!frame_cell_ok[o.frame]) continue; - if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; - if (o.partiality < min_partiality) continue; - const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; - if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; - const auto shell = shells.GetShell(o.d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Only what the merge kept counts towards multiplicity - see the GPU branch above. The - // R_meas accumulation below keeps its own, wider filter. - if (usable_merge(o)) - sa[*shell].total_obs++; - if (std::isfinite(merged_I[o.group])) { - auto &r = rmeas[o.group]; - r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); - r.sum_I += I_corr; r.n++; r.shell = *shell; - } - } - for (const auto &r : rmeas) { - if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; - const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); - rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; - rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; - } - } - - // ---- Anomalous split (always, even when the merge is Friedel-averaged): for each acentric - // reflection, the inverse-variance I(+)/I(-) from the SAME scaled fulls (so it never touches the - // Friedel-merged IMEAN, scaling or error model above). Lets I(+)/I(-) and F(+)/F(-) be written by - // default without scaling anomalously; a reflection with only one mate, or a centric, is left - // without the split. Skipped for the P1 search pass, which has no use for it, and for a merge that - // is not written. ---- - struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; }; - std::unordered_map anom_export; - // SigAno = <|I(+)-I(-)|>/ accumulated on the standard report shells (num=sum|dI|, - // den=sum sigma(dI)); mmCIF pdbx_absDiff_over_sigma_anomalous. NaN'd below when there is no split. - std::vector sig_num(n_shells, 0.0), sig_den(n_shells, 0.0); - double sig_num_all = 0.0, sig_den_all = 0.0; - size_t sig_n = 0; - // CCanom on the same shells, and over the whole range as one correlation (not a mean of shells). - std::vector cc_anom_shell(n_shells); - CorrelationCoefficient cc_anom_overall; - size_t cc_anom_n = 0; - if (!for_search && full_stats) { - const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, x.GetSpaceGroupOrP1()); - const gemmi::GroupOps gops = x.GetSpaceGroupOrP1().operations(); - // [hand][half]: hand 0 = I(+), 1 = I(-). The half sets are what CCanom needs - the anomalous - // difference formed twice, once from each half of the observations - and their sums add back - // to the whole-mate I(+)/I(-) that SigAno and the export read, so there is one accumulation. - struct AnomAcc { - double swI[2][2] = {}; double sw[2][2] = {}; int32_t nh[2][2] = {}; - int32_t h = 0, k = 0, l = 0; float d = NAN; - }; - // The Friedel-merged key an observation lands on, and which hand it is, are properties of its - // ASU GROUP, not of the observation - group_h/k/l is that group's SIGNED representative, so the - // same reduction applied to it returns both. Reduce once per group (a few hundred thousand) - // instead of once per observation (millions), and index a dense accumulator from there, which - // is the rule ComputeAsuGroups states for itself. The hand only follows the group when the - // merge distinguishes the hands; a Friedel-merged run holds both in one group and must ask. - std::vector pair_of_group(n_groups, -1); - std::vector hand_of_group(merge_friedel ? 0 : n_groups, 0); - std::vector pair_hkl; - { - std::unordered_map pair_id; - pair_id.reserve(static_cast(n_groups) + 1); - for (int g = 0; g < n_groups; ++g) { - const HKLKey gk = anom_keygen(group_h[g], group_k[g], group_l[g]); - if (!merge_friedel) hand_of_group[g] = gk.plus ? 0 : 1; - const uint64_t key = HKLKey{gk.h, gk.k, gk.l, true}.pack(); - const auto [it, fresh] = pair_id.emplace(key, static_cast(pair_hkl.size())); - if (fresh) pair_hkl.push_back(HKLKey{gk.h, gk.k, gk.l, true}); - pair_of_group[g] = it->second; - } - } - std::vector anom(pair_hkl.size()); - // The half-sets CCanom needs are balanced within the MATE, not merely within the reflection: - // it correlates one half's I(+)-I(-) against the other's, so each hand separately needs both - // of its halves filled, and a mate left entirely in one half loses the whole pair rather than - // half of it. So the cell here is (pair, hand), where the merge above used the reflection. - std::vector anom_cell(fulls.size(), -1); - for (size_t i = 0; i < fulls.size(); ++i) { - const Obs &o = fulls[i]; - if (!usable_merge(o) || rejected_obs[i]) continue; // as the merge saw it - const int hand = merge_friedel ? (anom_keygen(o.h, o.k, o.l).plus ? 0 : 1) - : hand_of_group[o.group]; - anom_cell[i] = pair_of_group[o.group] * 2 + hand; - } - std::vector anom_half; - AssignHalvesByRank(obs_key, anom_cell, static_cast(pair_hkl.size()) * 2, anom_half, nthreads); - for (size_t i = 0; i < fulls.size(); ++i) { - const int cell = anom_cell[i]; - if (cell < 0) continue; - const Obs &o = fulls[i]; - const float I_corr = o.I * o.corr; - const float sigma_corr = corrected_sigma(o, I_corr, o.sigma * o.corr); - if (!(sigma_corr > 0.0f) || !std::isfinite(sigma_corr)) continue; - const double w = 1.0 / (static_cast(sigma_corr) * sigma_corr); - AnomAcc &a = anom[cell >> 1]; - a.d = o.d; - a.swI[cell & 1][anom_half[i]] += w * static_cast(I_corr); - a.sw[cell & 1][anom_half[i]] += w; - a.nh[cell & 1][anom_half[i]]++; - } - for (size_t p = 0; p < anom.size(); ++p) { - AnomAcc a = anom[p]; - a.h = pair_hkl[p].h; a.k = pair_hkl[p].k; a.l = pair_hkl[p].l; - const uint64_t fkey = pair_hkl[p].pack(); - // Centrics have I(+)=I(-) by symmetry; leave them without an anomalous split. - if (gops.is_reflection_centric(gemmi::Op::Miller{a.h, a.k, a.l})) continue; - // Whole-mate sums: the two halves put back together. I(+)/I(-) and SigAno are read off - // these, so splitting the accumulation above changed neither. - const double swp = a.sw[0][0] + a.sw[0][1], swm = a.sw[1][0] + a.sw[1][1]; - const double swIp = a.swI[0][0] + a.swI[0][1], swIm = a.swI[1][0] + a.swI[1][1]; - AnomExport ex; - if (swp > 0.0) { ex.Ip = static_cast(swIp / swp); ex.sIp = static_cast(1.0 / std::sqrt(swp)); } - if (swm > 0.0) { ex.Im = static_cast(swIm / swm); ex.sIm = static_cast(1.0 / std::sqrt(swm)); } - if (std::isfinite(ex.Ip) || std::isfinite(ex.Im)) anom_export[fkey] = ex; - // A pair outside the grid is outside the table, its overall line included: the fulls - // hold every pair that was ingested, and the ones past the written range are not the - // data these numbers describe. - const auto shell = shells.GetShell(a.d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - if (swp > 0.0 && swm > 0.0) { - const double abs_di = std::fabs(swIp / swp - swIm / swm); - const double sig_di = std::sqrt(1.0 / swp + 1.0 / swm); - sig_num[*shell] += abs_di; sig_den[*shell] += sig_di; - sig_num_all += abs_di; sig_den_all += sig_di; ++sig_n; - } - // CCanom: the anomalous difference measured twice, once per half, correlated over the - // pairs where all four cells were filled - which needs both mates measured at least - // twice, so it asks more of the data than SigAno does. - if (a.sw[0][0] > 0.0 && a.sw[0][1] > 0.0 && a.sw[1][0] > 0.0 && a.sw[1][1] > 0.0) { - const double d_first = a.swI[0][0] / a.sw[0][0] - a.swI[1][0] / a.sw[1][0]; - const double d_second = a.swI[0][1] / a.sw[0][1] - a.swI[1][1] / a.sw[1][1]; - cc_anom_shell[*shell].Add(d_first, d_second); - cc_anom_overall.Add(d_first, d_second); - ++cc_anom_n; - } - } - } - + result.statistics = statistics_over(result.merged, d_min_pad, d_max_declared, d_max_label, 0.0f, true); MergeStatistics &out = result.statistics; - out.shells.resize(n_shells); - for (int s = 0; s < n_shells; ++s) { - auto &ss = out.shells[s]; - ss.mean_one_over_d2 = shell_mean_1_d2[s]; - ss.d_min = shell_min_res[s]; - ss.d_max = s == 0 ? d_max_label : shell_min_res[s - 1]; - ss.total_observations = sa[s].total_obs; - ss.unique_reflections = sa[s].unique; - ss.possible_unique_reflections = sa[s].possible; - ss.mean_i_over_sigma = sa[s].n_i_over_sigma > 0 ? sa[s].sum_i_over_sigma / sa[s].n_i_over_sigma : 0.0; - ss.cc_half = sa[s].cc_half.GetCC(); - ss.cc_ref = NAN; - ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; - ss.abs_diff_over_sigma_anomalous = sig_den[s] > 0.0 ? sig_num[s] / sig_den[s] : NAN; - ss.cc_anom = cc_anom_shell[s].GetCC(); - } - auto &overall = out.overall; - overall.d_min = sd_min; overall.d_max = sd_max; - double sum_ios = 0.0; int n_ios = 0; - for (int s = 0; s < n_shells; ++s) { - overall.total_observations += sa[s].total_obs; - overall.unique_reflections += sa[s].unique; - overall.possible_unique_reflections += sa[s].possible; - sum_ios += sa[s].sum_i_over_sigma; n_ios += sa[s].n_i_over_sigma; - } - overall.mean_i_over_sigma = n_ios > 0 ? sum_ios / n_ios : 0.0; - overall.cc_half = cc_half_overall.GetCC(); - overall.cc_ref = NAN; - overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; - overall.abs_diff_over_sigma_anomalous = sig_den_all > 0.0 ? sig_num_all / sig_den_all : NAN; - overall.cc_anom = cc_anom_overall.GetCC(); - if (std::isfinite(overall.abs_diff_over_sigma_anomalous)) - logger.Info("Anomalous signal SigAno = {:.2f} ({} acentric pairs), CCanom = {:.3f} ({} pairs " - "split in both hands)", - overall.abs_diff_over_sigma_anomalous, sig_n, overall.cc_anom, cc_anom_n); + out.overall.d_min = sd_min; out.overall.d_max = sd_max; // Radiation-damage monitor (measured before any correction by MeasureRadiationDamageB): carry the // first->last relative-B change and the per-batch curve into the reported statistics / mmCIF. @@ -4363,6 +4381,39 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool out.radiation_damage_batch_deg = rad_damage_batch_deg; out.sweep_quality = sweep_quality; + // The reference-range table (--report-resolution): the reflections the run KEPT, binned over the + // reference range instead of the run's own. Nothing past the run's own cut is merged into it - + // the run decided those reflections are not signal, and a table that read them anyway would + // describe data the run's output does not contain - so where the reference range is finer than + // the cut its outer shells are empty and the overall numbers are over what was kept. The ISa is + // the error model refitted on the samples of exactly that table, in XDS's convention, so it can + // be read against a program that was run at that range; report-only, and fitted last because + // nothing after it reads the model. + if (report_range && full_stats && !for_search) { + const float d_lo = static_cast(report_range->d_min); + const float d_hi = static_cast(report_range->d_max); + float in_min = std::numeric_limits::max(), in_max = 0.0f; + for (const auto &m : result.merged) { + if (!std::isfinite(m.d) || m.d <= d_lo || m.d > d_hi) continue; + in_min = std::min(in_min, m.d); in_max = std::max(in_max, m.d); + } + MergeStatistics ref = statistics_over(result.merged, d_lo, d_hi, std::isfinite(d_hi) ? d_hi : in_max, + static_cast(effective_d_min.value_or(0.0)), false); + ref.overall.d_min = in_min < in_max ? in_min : 0.0f; + ref.overall.d_max = in_max; + result.reference_statistics = std::move(ref); + + std::vector pool; + pool.reserve(samples.size()); + for (const auto &s : samples) + if (s.d > d_lo && s.d <= d_hi && !(effective_d_min && s.d < *effective_d_min)) + pool.push_back(s); + error_model_active = false; // a pool too small to fit leaves no model, not the run's own + fit_error_model(pool); + if (error_model_active) + result.reference_isa = ToXdsErrorModel(error_model_a, error_model_b).isa; + } + // Attach the per-reflection anomalous split so the writer can emit I(+)/I(-) by default (each merged // reflection maps to its Friedel-ASU key; in an anomalous merge both mates map to the same key). if (!anom_export.empty()) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index f7cdccf40..24a1fed00 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -62,6 +62,14 @@ public: // Where the automatic CC1/2 fit says the data reach, without the deliberate one-shell // extension the reflections are written to. Empty when no automatic cut ran. std::optional resolution_fit_A; + // The merged reflections the run kept, binned over the reference range of + // ScalingSettings::ReportResolutionRange_A (--report-resolution) instead of the run's own, + // and the error model refitted on the observations of that table alone, for its ISa. Nothing + // past the run's own cut is in it: those shells are empty, and the overall numbers are over + // what the run kept. Absent when no range was given, on a search merge, and on a merge whose + // numbers are not written. + std::optional reference_statistics; + double reference_isa = 0.0; }; // experiment: read live (its space group is changed by the caller between Run() calls). @@ -205,6 +213,7 @@ private: ResolutionCutoffMethod resolution_cutoff_method = ResolutionCutoffMethod::Off; double resolution_cc_target = 0.30; int report_shell_count = 9; + std::optional report_range; // the second, reference-range table // Flat buffers, allocated once by Ingest() and reused across Run() calls. Exactly one of the // two observation arrays is built: the narrow resident record when the device pipeline is diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index af0fddec6..f3d6b25e4 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -181,6 +181,29 @@ namespace { return out; } + // The reference-range table of --report-resolution, with the columns of the run's own. A shell + // finer than the run's own limit holds nothing BY THE RUN'S DECISION, and says so in words rather + // than as a row of zeros, so a reader cannot take it for a measured total failure. + std::string ReferenceShellTable(const MergeStatistics &rs, double own_d_min) { + std::ostringstream o; + o << "\n" << ShellTableHeader() << "\n" << ShellTableRule() << "\n"; + for (const auto &sh : rs.shells) { + if (sh.unique_reflections == 0) { + if (sh.d_max <= own_d_min) + o << fmt::format(" {:8.2f} past the run's own limit of {:.2f} A: not merged ({} possible)\n", + sh.d_min, own_d_min, sh.possible_unique_reflections); + continue; + } + o << fmt::format(" {:8.2f} ", sh.d_min) << sh; + if (sh.d_min < own_d_min) + o << fmt::format(" <- the run's own limit, {:.2f} A, is inside this shell", own_d_min); + o << "\n"; + } + o << ShellTableRule() << "\n"; + o << fmt::format(" {:>8s} ", "Overall") << rs.overall << "\n"; + return o.str(); + } + // The finest shell whose CC1/2 still reaches the cut's own target. A fitted resolution finer than // this was read off a curve that never behaved like a fall-off, and asserting it states a // resolution the shell table on the same page refutes. @@ -705,6 +728,71 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Prose(shells.str())); } + // The reference-range table (--report-resolution): the same kept reflections binned over + // a range given from outside, so the run can be read against another program's table at + // that program's range without having been RUN at it - which moves the cut, the symmetry + // decision and everything after them. Nothing past the run's own limit is in it: those + // shells are empty by the run's decision, REFRES_COMPLETENESS counts their reflections as + // missing, and REFRES_SHELLS_PAST_LIMIT says how many there are, so a consumer can tell + // "not merged" from a measured zero. + if (result.reference_merge_statistics) { + const auto &rs = *result.reference_merge_statistics; + const auto &ro = rs.overall; + const auto range = experiment.GetScalingSettings().GetReportResolutionRange_A(); + const double d_max_req = range ? range->d_max : rs.shells.front().d_max; + const double d_min_req = range ? range->d_min : rs.shells.back().d_min; + int past = 0; + for (const auto &sh : rs.shells) + if (sh.d_max <= o.d_min) + ++past; + Add(s, Blank()); + Add(s, KeyText("REFRES_RANGE", fmt::format( + "{} {:.3f}", std::isfinite(d_max_req) ? fmt::format("{:.3f}", d_max_req) : "INF", + d_min_req))); + if (ro.d_min > 0.0f) + Add(s, KeyText("REFRES_MEASURED_RANGE", fmt::format("{:.3f} {:.3f}", ro.d_max, ro.d_min))); + Add(s, KeyInt("REFRES_SHELLS_PAST_LIMIT", past)); + Add(s, KeyInt("REFRES_UNIQUE_REFLECTIONS", ro.unique_reflections)); + Add(s, KeyInt("REFRES_TOTAL_OBSERVATIONS", ro.total_observations)); + if (ro.possible_unique_reflections > 0) + Add(s, KeyReal("REFRES_COMPLETENESS", + 100.0 * ro.unique_reflections / ro.possible_unique_reflections, "{:.1f}")); + if (ro.unique_reflections > 0) { + Add(s, KeyReal("REFRES_MULTIPLICITY", + static_cast(ro.total_observations) / ro.unique_reflections, "{:.2f}")); + Add(s, KeyReal("REFRES_I_OVER_SIGMA", ro.mean_i_over_sigma, "{:.2f}")); + } + if (std::isfinite(ro.r_meas)) + Add(s, KeyReal("REFRES_R_MEAS", ro.r_meas, "{:.4f}")); + if (std::isfinite(ro.cc_half)) + Add(s, KeyReal("REFRES_CC_HALF", ro.cc_half, "{:.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)) + Add(s, KeyReal("REFRES_CC_ANOM", ro.cc_anom, "{:.4f}")); + if (result.reference_isa > 0.0) + Add(s, KeyReal("REFRES_ISA", result.reference_isa, "{:.2f}")); + Add(s, Blank()); + Add(s, Prose(fmt::format( + " REFRES_* is the same merged data binned over REFRES_RANGE, the reference range given\n" + " with --report-resolution, so it can be read against another program's table at that\n" + " program's range; nothing was processed differently for it. It holds only what the run\n" + " kept: {}{}", + past > 0 + ? fmt::format("the run's own limit is {:.2f} A (INCLUDE_RESOLUTION_RANGE), the {} shell{}\n" + " finer than it hold{} nothing because the run judged them to carry no signal and did\n" + " not merge them, REFRES_COMPLETENESS counts their reflections as missing, and every\n" + " other REFRES_ number is over the shells the run reached.", + o.d_min, past, past == 1 ? "" : "s", past == 1 ? "s" : "") + : std::string("the range lies within the run's own, so this is a subset of the table\n" + " above."), + result.reference_isa > 0.0 + ? "\n REFRES_ISA is the error model refitted on the reflections of this table alone, in\n" + " XDS's convention." + : ""))); + Add(s, Prose(ReferenceShellTable(rs, o.d_min))); + } + // A merge with no signal in it. Neither of these is caught anywhere else: CC1/2 near zero // beside a healthy I/sigma is the signature of observations of DIFFERENT reflections being // averaged together (mixed indexing hands, a wrong point group), and a rejection count diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index bf5518985..011d69512 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4705,6 +4705,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // quality guard in Run() compares the passes on. Equal to the reported CC1/2 wherever no // surface was fitted, which is every merge on the stills path. double cc_half_before_corrections = NAN; + // The reference-range table and its ISa (--report-resolution); see ProcessResult. + std::optional reference_statistics; + double reference_isa = 0.0; }; // The reference path computes each image's G once (per-image scaling against the // reference); the scaling loop below is skipped, so G is stable across the two passes. @@ -4917,7 +4920,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.error_model_a = r.error_model_a; result.error_model_b = r.error_model_b; return ScaleMergeResult{std::move(r.merged), std::move(r.statistics), - r.cc_half_before_corrections}; + r.cc_half_before_corrections, + std::move(r.reference_statistics), r.reference_isa}; } // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the // running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even @@ -4993,6 +4997,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b phase("Computing statistics"); out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min); + // The same kept reflections binned over the reference range (--report-resolution); the + // stills error model is fitted over the whole range already, so no ISa of its own. + if (const auto range = cut_ss.GetReportResolutionRange_A(); range && !for_search) + out.reference_statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, + std::nullopt, range); { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); result.error_model_isa = em.isa; @@ -6412,6 +6421,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.merge_statistics_text = stats_text.str(); result.has_merge_statistics = true; result.merge_statistics = sm.statistics; + result.reference_merge_statistics = sm.reference_statistics; + result.reference_isa = sm.reference_isa; result.cc_half_before_corrections = sm.cc_half_before_corrections; // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8365b84e1..e992ffa3b 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -243,6 +243,15 @@ struct ProcessResult { double error_model_b = 0.0; bool has_reference = false; + // The reflections the run kept, binned over the reference range of --report-resolution + // (ScalingSettings::ReportResolutionRange_A) instead of the run's own range: a second table, not + // a second processing, so a run can be read against another program's numbers at that program's + // range. Shells finer than the run's own limit are empty. reference_isa is the error model + // refitted on that table's reflections alone (rotation; 0 where not measured). Absent when no + // range was given. + std::optional reference_merge_statistics; + double reference_isa = 0.0; + // The overall CC1/2 as it stood before the scaling correction surfaces were folded in (NaN when no // merge ran). The rotation two-pass quality guard judges the refined pass on this rather than on the // reported CC1/2: the geometry pre-pass fits no surfaces, so only the uncorrected number is the same diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 2f83c90a1..3896fbd6c 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -4,6 +4,7 @@ #include "RugnuxCommandLine.h" #include "../common/DiffractionExperiment.h" +#include #include #include @@ -221,6 +222,10 @@ std::string RugnuxCommandLine(const ProcessConfig &config, const auto d_max_default = ScalingSettings().GetLowResolutionLimit_A(); if (d_max != d_max_default) args.push_back("--scaling-low-resolution=" + num(d_max.value_or(0.0))); + // An infinite d_max is the default of a run with no low limit, so d_min alone restates it. + if (const auto r = sc.GetReportResolutionRange_A()) + args.push_back("--report-resolution=" + num(r->d_min) + + (std::isfinite(r->d_max) ? "," + num(r->d_max) : "")); } else { args.emplace_back("--no-merge"); } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 20ee52170..ecd6e07cc 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -164,6 +164,7 @@ void print_usage() { std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl; std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl; std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 9, as XDS)" << std::endl; + std::cout << " --report-resolution [,] Also report the merging statistics over this resolution range (A, either order; dmax defaults to the run's own low-resolution limit) as a second table beside the run's own - the REFRES_* keys of the report - binned from the same merged reflections, so a run can be read against another program's table at that program's range. Report-only: the cut, the scaling, the error model and every decision stay the run's own, and shells finer than the run's own limit are empty rather than filled from data the run did not keep" << std::endl; std::cout << " --ice-min-score Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging and the exclusion from scaling (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl; std::cout << " --ice-min-spot-ratio Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl; std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; @@ -239,6 +240,7 @@ enum { OPT_RESOLUTION_CUTOFF, OPT_RESOLUTION_CC_TARGET, OPT_RESOLUTION_SHELLS, + OPT_REPORT_RESOLUTION, OPT_SINGLE_PASS_ROTATION, OPT_FORCE_ROTATION_LATTICE, OPT_MIN_INDEXED_SPOTS, @@ -407,6 +409,7 @@ static option long_options[] = { {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET}, {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS}, + {"report-resolution", required_argument, nullptr, OPT_REPORT_RESOLUTION}, {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"adaptive-integration-radius", optional_argument, nullptr, OPT_ADAPTIVE_INTEGRATION_RADIUS}, @@ -494,6 +497,30 @@ float parse_float_arg(const char *arg, const char *option_name, Logger &logger) return parse_number_arg(arg, option_name, logger); } +// --report-resolution [,], in A and in either order; d_max defaults to the run's own +// low-resolution limit. Report-only - see ScalingSettings::ReportResolutionRange_A. +ReportResolutionRange ParseReportResolution(const std::string &arg, const ScalingSettings &scaling, + Logger &logger) { + std::vector v; + std::stringstream ss(arg); + std::string tok; + while (std::getline(ss, tok, ',')) { + trim_in_place(tok); + if (!tok.empty()) + v.push_back(parse_number_arg(tok.c_str(), "--report-resolution", logger, 0.1, 100000.0)); + } + if (v.empty() || v.size() > 2 || (v.size() == 2 && v[0] == v[1])) { + logger.Error("--report-resolution expects d_min or d_min,d_max in A: '{}'", arg); + print_usage(); + exit(EXIT_FAILURE); + } + ReportResolutionRange r; + r.d_min = v.size() == 2 ? std::min(v[0], v[1]) : v[0]; + r.d_max = v.size() == 2 ? std::max(v[0], v[1]) + : scaling.GetLowResolutionLimit_A().value_or(std::numeric_limits::infinity()); + return r; +} + bool parse_on_off(const char *arg, bool &out) { std::string s = arg ? arg : ""; std::transform(s.begin(), s.end(), s.begin(), @@ -764,6 +791,7 @@ static int RunRugnux(int argc, char **argv) { std::optional resolution_cutoff_method; // --resolution-cutoff cc-logistic|off std::optional resolution_cc_target; // --resolution-cc-target std::optional report_shell_count; // --resolution-shells + std::optional report_resolution_arg; // --report-resolution std::optional integration_radius_arg; // --adaptive-integration-radius: unset means the per-workflow default below (on for rotation). std::optional adaptive_integration_radius_arg; @@ -1255,6 +1283,9 @@ static int RunRugnux(int argc, char **argv) { case OPT_RESOLUTION_CC_TARGET: resolution_cc_target = parse_double_arg(optarg, "--resolution-cc-target", logger); break; + case OPT_REPORT_RESOLUTION: + report_resolution_arg = optarg; + break; case OPT_RESOLUTION_SHELLS: report_shell_count = atoi(optarg); if (report_shell_count.value() < 1) { @@ -1543,6 +1574,9 @@ static int RunRugnux(int argc, char **argv) { if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); + if (report_resolution_arg) + scaling_settings.ReportResolutionRange_A( + ParseReportResolution(*report_resolution_arg, scaling_settings, logger)); scaling_settings.MergeFriedel(!anomalous_mode); scaling_settings.MinPartiality(min_partiality); scaling_settings.MinCapturedFraction( @@ -1696,6 +1730,8 @@ static int RunRugnux(int argc, char **argv) { std::optional resolution_fit_A; // CC1/2 crossing, without the one-shell extension double error_model_a = 0.0; double error_model_b = 0.0; + std::optional reference_statistics; // --report-resolution + double reference_isa = 0.0; // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine -> // merge, including the default-on decay + absorption correction surfaces. It does not support @@ -1722,6 +1758,8 @@ static int RunRugnux(int argc, char **argv) { error_model_a = r.error_model_a; error_model_b = r.error_model_b; resolution_fit_A = r.resolution_fit_A; + reference_statistics = std::move(r.reference_statistics); + reference_isa = r.reference_isa; } else { // Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports // CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image @@ -1768,6 +1806,9 @@ static int RunRugnux(int argc, char **argv) { merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, effective_d_min); + if (const auto range = cut_ss.GetReportResolutionRange_A()) + reference_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data, + std::nullopt, range); { // stills: no asymptote tier, and its ISa already IS the whole-range 1/b const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB()); error_model_isa = em.isa; @@ -1952,6 +1993,8 @@ static int RunRugnux(int argc, char **argv) { scale_result.error_model_a = error_model_a; scale_result.error_model_b = error_model_b; scale_result.resolution_fit_A = resolution_fit_A; + scale_result.reference_merge_statistics = reference_statistics; + scale_result.reference_isa = reference_isa; scale_result.has_reference = !reference_data.empty(); scale_result.twinning = twinning; scale_result.model_validation = model_validation; @@ -2439,6 +2482,9 @@ static int RunRugnux(int argc, char **argv) { if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); + if (report_resolution_arg) + scaling_settings.ReportResolutionRange_A( + ParseReportResolution(*report_resolution_arg, scaling_settings, logger)); scaling_settings.MergeFriedel(!anomalous_mode); scaling_settings.MinPartiality(min_partiality); // Drop edge-of-sweep truncated fulls (rocking curve captured < this fraction) from the rot3d combine. diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 5ca2294cb..31dd7676a 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -259,6 +259,47 @@ TEST_CASE("MergeStats_NoReferenceCellLeavesCompletenessUnmeasured") { } } +// The reference-range table (--report-resolution) is binned over the range it is GIVEN, whether or +// not the data reach it: its last shell ends at the declared d_min, and past the run's own limit +// nothing is read from anywhere - the reflections the run did not keep count as missing. +TEST_CASE("MergeStats_ReferenceRangeIsBinnedOverTheDeclaredRange") { + const auto &sg = TestSpaceGroup(); + const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true); // what the run kept + DiffractionExperiment x; + x.SetSpaceGroup(sg); + MergeOnTheFly merge(x); + merge.ReferenceCell(TEST_CELL); + const auto own = merge.MergeStats(merged, {}); + + // A range coarser than the run's own: a complete subset, on a grid ending at the declared bound. + const auto within = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{3.0, 50.0}); + CHECK(within.shells.back().d_min == Catch::Approx(3.0)); + CHECK(within.shells.front().d_max == Catch::Approx(50.0)); + int coarser_than_3 = 0; + for (const auto &m: merged) + if (m.d > 3.0f) ++coarser_than_3; + CHECK(within.overall.unique_reflections == coarser_than_3); + CHECK(within.overall.unique_reflections < own.overall.unique_reflections); + CHECK(Completeness(within.overall) > 99.0); + + // A range finer than the run kept: the grid still ends at 1.5 A, the outer shell is empty and its + // possible reflections are in the denominator, and the measured range says where the data stop. + const auto beyond = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.5, 50.0}); + CHECK(beyond.shells.back().d_min == Catch::Approx(1.5)); + CHECK(beyond.overall.unique_reflections == own.overall.unique_reflections); + CHECK(beyond.overall.possible_unique_reflections > own.overall.possible_unique_reflections); + CHECK(Completeness(beyond.overall) < Completeness(own.overall)); + CHECK(beyond.shells.back().unique_reflections == 0); + CHECK(beyond.shells.back().possible_unique_reflections > 0); + CHECK(beyond.overall.d_min == Catch::Approx(own.overall.d_min)); + + // A range the data never reach at all is still a grid, not an error. + const auto empty = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.0, 1.5}); + CHECK(empty.overall.unique_reflections == 0); + CHECK(empty.overall.possible_unique_reflections > 0); + CHECK(empty.overall.d_min == 0.0f); +} + // ---------------------------------------------------------------- the automatic resolution cutoff namespace { // Half-set pairs spread uniformly in s = 1/d^2 over [s_from, s_to), either correlated with each diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 9dcec514a..5bafc0525 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -157,6 +157,61 @@ TEST_CASE("ResultReport_Render", "[Diagnostics]") { != std::string::npos); } +// The reference-range table: REFRES_* keys beside the run's own, the shells past the run's own limit +// reported as not merged rather than as zeros, and no REFRES_ key at all when no range was given. +TEST_CASE("ResultReport_ReferenceRangeTable", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + ProcessResult result; + result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; + result.space_group = *gemmi::find_spacegroup_by_number(96); + result.has_merge_statistics = true; + auto &own = result.merge_statistics.overall; + own.d_max = 45.0f; own.d_min = 1.6f; // the run's own limit + own.unique_reflections = 1000; own.possible_unique_reflections = 1000; own.total_observations = 5000; + own.mean_i_over_sigma = 10.0; own.r_meas = 0.1; own.cc_half = 0.99; + + CHECK(RenderResultReport("prefix", "in.h5", x, result).find("REFRES_") == std::string::npos); + + ScalingSettings s = x.GetScalingSettings(); + s.ReportResolutionRange_A(ReportResolutionRange{1.2, 50.0}); + x.ImportScalingSettings(s); + MergeStatistics rs; + auto shell = [](float d_max, float d_min, int unique, int possible, int obs) { + MergeStatisticsShell sh; + sh.d_max = d_max; sh.d_min = d_min; + sh.unique_reflections = unique; sh.possible_unique_reflections = possible; sh.total_observations = obs; + sh.mean_i_over_sigma = unique > 0 ? 8.0 : 0.0; + sh.cc_half = unique > 0 ? 0.9 : NAN; + sh.r_meas = unique > 0 ? 0.2 : NAN; + return sh; + }; + rs.shells = {shell(50.0f, 2.0f, 600, 600, 3000), // within the run's own range + shell(2.0f, 1.5f, 300, 500, 1200), // the run's limit of 1.6 A falls inside + shell(1.5f, 1.2f, 0, 700, 0)}; // past it: not merged + rs.overall = shell(45.0f, 1.6f, 900, 1800, 4200); + result.reference_merge_statistics = rs; + result.reference_isa = 12.5; + + const auto text = RenderResultReport("prefix", "in.h5", x, result); + CHECK(text.find("\nREFRES_RANGE= 50.000 1.200\n") != std::string::npos); + CHECK(text.find("\nREFRES_MEASURED_RANGE= 45.000 1.600\n") != std::string::npos); + CHECK(text.find("\nREFRES_SHELLS_PAST_LIMIT= 1\n") != std::string::npos); + CHECK(text.find("\nREFRES_UNIQUE_REFLECTIONS= 900\n") != std::string::npos); + CHECK(text.find("\nREFRES_COMPLETENESS= 50.0\n") != std::string::npos); + CHECK(text.find("\nREFRES_MULTIPLICITY= 4.67\n") != std::string::npos); + CHECK(text.find("\nREFRES_R_MEAS= 0.2000\n") != std::string::npos); + CHECK(text.find("\nREFRES_CC_HALF= 0.9000\n") != std::string::npos); + CHECK(text.find("\nREFRES_ISA= 12.50\n") != std::string::npos); + // The run's own keys are untouched by the second table. + CHECK(text.find("\nCOMPLETENESS= 100.0\n") != std::string::npos); + CHECK(text.find("\nINCLUDE_RESOLUTION_RANGE= 45.000 1.600\n") != std::string::npos); + // The shell past the limit is words, not zeros; the shell holding the limit is marked. + CHECK(text.find("1.20 past the run's own limit of 1.60 A: not merged (700 possible)") != std::string::npos); + CHECK(text.find("<- the run's own limit, 1.60 A, is inside this shell") != std::string::npos); + CHECK(text.find("the 1 shell\n") != std::string::npos); +} + TEST_CASE("ResultReport_RenderEmpty", "[Diagnostics]") { // A clean run and a run that never looked must be distinguishable: both have a count of 0, and // only the STATUS key separates them. This is the property a consumer relies on.