From fd5a7cee5d63b3fa993ebd87cbfaaf48219994b8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 1 Nov 2026 21:45:35 +0100 Subject: [PATCH] Report the standard unweighted CC1/2; keep the weighted form internal Every reported CC1/2 on rotation - CC_HALF, CC_STAR, REFRES_CC_HALF, the CC1/2 and CC* columns of the shell tables, the mmCIF and the viewer - is now the plain Pearson correlation of the 1/sigma^2-weighted half-set means, as XDS, Aimless and DIALS report it. The per-reflection cc_weight (frame scale against the typical scale) stays in what the processing decides on: the automatic resolution cut, the batch delta-CC1/2 ledger, the two-pass and beam-centre guards (MergeStatisticsShell::cc_half_internal), and the report's judgement of whether FITTED_RESOLUTION is quotable. Stills were already unweighted. Adds the inert diagnostic report keys FITTED_RESOLUTION_PLAIN and RESOLUTION_CUT_PLAIN (rotation): the same logistic fit on the plain CC1/2 and the limit it would have written. Nothing is cut at them. Merged and unmerged MTZ and .hkl are bit-identical to before on three rotation sets with weak stretches and one stills set; the written d_min is unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- common/Reflection.h | 2 ++ docs/CPU_DATA_ANALYSIS_DECISIONS.md | 2 ++ docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 2 +- docs/RUGNUX_REPORT.md | 11 +++++- image_analysis/scale_merge/Merge.cpp | 2 ++ image_analysis/scale_merge/Merge.h | 9 +++++ .../scale_merge/ResolutionCutoff.cpp | 6 ++-- image_analysis/scale_merge/ResolutionCutoff.h | 7 ++-- .../scale_merge/RotationScaleMerge.cpp | 35 ++++++++++++++----- .../scale_merge/RotationScaleMerge.h | 10 ++++-- rugnux/ResultReport.cpp | 29 +++++++++++---- rugnux/Rugnux.h | 6 ++++ rugnux/RugnuxScaleMerge.cpp | 10 ++++-- tests/MergeScaleTest.cpp | 4 +++ tests/ResultReportTest.cpp | 29 ++++++++++----- tests/RotationScaleMergeSegmentsTest.cpp | 24 +++++++++++++ 16 files changed, 156 insertions(+), 32 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index 9c864b8f0..9e8f9abe6 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -100,6 +100,8 @@ struct MergedReflection { // 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 // diffracted, whose scaled-up noise would otherwise count as much as a well-measured pair. + // Internal, not reported: it enters the CC1/2 the processing decides on (the automatic resolution + // cut, MergeStatisticsShell::cc_half_internal), never the reported CC1/2, which is the plain Pearson. float cc_weight = 1.0f; float d = 0.0; // Any observation of this reflection sat on an ice ring. The intensity is still merged and diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index baac5c2f9..abb73189b 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -82,6 +82,8 @@ The median needs three observations, so a reflection measured once or twice — By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off: a logistic is fitted to $\mathrm{CC}_{1/2}(s)$, and the limit is set **one reported-shell width past** the point where the fit crosses 0.30 — deliberately "one shell too far", so weak-but-real data below the crossing are kept rather than discarded. The extension is measured over the range that is actually kept, not the full measured range, so a detector reaching far past where the crystal diffracts cannot inflate it. `--scaling-high-resolution` overrides the limit and `--resolution-cutoff off` disables it. +**The fit runs on an internal, weighted $\mathrm{CC}_{1/2}$, not on the reported one.** An observation from a frame at scale $G$ below the run's typical scale $G_\mathrm{ref}$ was scaled up by $G_\mathrm{ref}/G$ and carries $(G_\mathrm{ref}/G)^2$ more counting variance than the same observation at the typical scale. On rotation each reflection's half-set pair therefore enters the fit's Pearson sums with a weight in $(0, 1]$: the variance its two half-set means would have had with every observation at the typical scale, over the variance they have. It is 1 on a sweep without a weak stretch, and small for a reflection measured only where the crystal barely diffracted, whose scaled-up noise would otherwise pull the curve down as much as a well-measured pair holds it up. The same per-frame factor sets the information each observation carries in the ΔCC1/2 frame ledger (§ frame disposition in [integration](CPU_DATA_ANALYSIS_INTEGRATION.md)), and the two-pass and beam-centre guards read this weighted form too. The **reported** $\mathrm{CC}_{1/2}$ — `CC_HALF`, `CC_STAR`, the shell tables and the mmCIF — is the standard one of XDS, Aimless and DIALS: the plain Pearson correlation of the $1/\sigma^2$-weighted half-set means, every reflection counted alike. Stills weight nothing in either. + **Ice sits out of this fit, and only this fit.** A powder ring is reproducible: past the crystal's own limit the ice is still there and still the same in both half-sets, so the two halves agree *about the ring*, and a Pearson $\mathrm{CC}_{1/2}$ cannot tell that agreement from diffraction. Shells have been measured at $\mathrm{CC}_{1/2}=0.765$ where $\langle I/\sigma\rangle$ is $-0.1$ and $R_\mathrm{meas}$ is 470 %, inside the written range because the logistic was fitted through them; against archived references, runs whose frames carry ice were written a mean twelve per cent finer than the reference where runs without ice were written four per cent finer, and every over-claim past twenty per cent has rings. The same ring drags the curve the other way where the data are good — one shell of forty falling to 0.51 because the two strongest ice lines cross it, with the shells either side at 0.99 — which is a real defect of those reflections but narrower than the shell it defames. Both signs are the same cause and both leave the fit: the reflections flagged in §3.3 are still merged, written and counted in the shell table wherever the per-ring decision of §10.10 keeps them, and what changes is only that the resolution **decision** now reads the same curve that scaling, the error model and the space-group search already read. **A cut the shell table refutes is not quoted.** Two things keep the number inside what the bins show. The fitted crossing has to lie **inside the fitted bins**: on a real fall-off it sits between the last bin at the target and the first bin below it, with the extension bins to spare, and a crossing past the last fitted bin is an extrapolation of a fall-off these data never showed — the crossing is then read off the bins themselves, between the two that straddle the target. And no single resolution is quoted at all where the curve is not a fall-off: where $\mathrm{CC}_{1/2}$ never reaches the target in any shell, where the fit reads finer than the finest shell whose $\mathrm{CC}_{1/2}$ still reaches it, or where $\mathrm{CC}_{1/2}$ **climbs back over the target** after falling below it — a noise shell that climbs back is not the edge of the data. The report then says which of the three happened and points at the shell table instead. diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 82aee398c..b5da3d702 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -432,7 +432,7 @@ Each surface is **cross-validated** on the half-set agreement it is meant to imp Each batch's $B$ is fitted on **resolution-shell means**, not on single observations: $\ln(I_\mathrm{ref}/I_\mathrm{obs})$ of one weak observation is unbounded and biased downwards — the observation appears in the response and in its own weight, and the logarithm needs $I_\mathrm{obs} > 0$, which keeps only the upward half of the noise — and on decayed data that bias grows with dose until it reverses the sign of the answer. The shells are laid inside the range the run actually diffracted to, and the fit carries an intercept as well as a slope, so a batch that is merely *dimmer* than the run (an attenuated beam, a mis-fitted frame scale) is not reported as damage. A batch whose shells are too weak to fit, or whose solved value reaches the bound the smoothing solve clamps to, is reported as **absent** rather than as a number. The first→last headline is reported only where a straight line describes the curve: radiation damage is progressive, so a curve that dips and recovers is a disturbance rather than dose, and is left to the sweep-quality report ([the results report](RUGNUX_REPORT.md)) to name. -**Frame disposition (rotation).** After the correction surfaces are fitted, and on the corrected fulls, Rugnux measures **ΔCC1/2** — the overall CC1/2 of the merged data with a group of images minus the CC1/2 without it, over the reflections that group touches — for each 10° batch of the sweep and for each stretch the sweep-quality diagnostic flagged. It is computed in the σ-τ form (no random half-dataset split, so the answer is the same every run) with each reflection's error variance taken from the **observed** scatter of its own observations rather than from the error model's σ's: a bad stretch claims the same σ's as a good one, so an error-model estimate would read a batch that adds noise as a batch that adds precision, inverting the sign of the measurement. The leave-one-out is a subtraction of the group's own $(n, \sum I, \sum I^2)$ from the per-reflection totals, so measuring a group costs one pass over its observations rather than a re-merge, and reflections the group holds the only observations of are excluded from both sides — the published statistic's own restriction. Its standard error is reported beside it as that of a single CC1/2 on the same reflection count, $(1-CC^2)/\sqrt{n_\mathrm{refl}-3}$, and a ΔCC1/2 smaller than that says nothing; a group whose ΔCC1/2 is positive or near zero is not evidence of harm. A batch is removed only where its ΔCC1/2 is **both** several standard errors below zero (Fisher-transformed) **and** well below the run's own per-batch distribution (median − 3 robust σ, measured once before anything is removed), **and** where the per-image CC to the merge — an independent channel, measured on the partials one frame at a time — also says the stretch agrees with the run worse than a typical frame does, because selecting frames by their disagreement with the merge and then reporting that the merge improved is circular; its edges are then slid frame by frame with the whole stretch re-measured at each position — so the range is located finely while the count of reflections the verdict rests on stays that of the stretch — the worst range is removed, the reference is re-formed and the scan repeats, never past a quarter of the sweep and never over a stretch narrower than one rocking event. The resulting per-frame `merged` / `downgraded` / `rejected` ledger is described in [the results report](RUGNUX_REPORT.md). +**Frame disposition (rotation).** After the correction surfaces are fitted, and on the corrected fulls, Rugnux measures **ΔCC1/2** — the overall CC1/2 of the merged data with a group of images minus the CC1/2 without it, over the reflections that group touches — for each 10° batch of the sweep and for each stretch the sweep-quality diagnostic flagged. It is computed in the σ-τ form (no random half-dataset split, so the answer is the same every run), with each observation carrying the information $\min(1, (G/G_\mathrm{ref})^2)$ of its frame's scale against the typical one — the internal weighting of the automatic resolution cut ([decisions](CPU_DATA_ANALYSIS_DECISIONS.md) §13.4), not part of the reported $\mathrm{CC}_{1/2}$ — and with each reflection's error variance taken from the **observed** scatter of its own observations rather than from the error model's σ's: a bad stretch claims the same σ's as a good one, so an error-model estimate would read a batch that adds noise as a batch that adds precision, inverting the sign of the measurement. The leave-one-out is a subtraction of the group's own $(n, \sum I, \sum I^2)$ from the per-reflection totals, so measuring a group costs one pass over its observations rather than a re-merge, and reflections the group holds the only observations of are excluded from both sides — the published statistic's own restriction. Its standard error is reported beside it as that of a single CC1/2 on the same reflection count, $(1-CC^2)/\sqrt{n_\mathrm{refl}-3}$, and a ΔCC1/2 smaller than that says nothing; a group whose ΔCC1/2 is positive or near zero is not evidence of harm. A batch is removed only where its ΔCC1/2 is **both** several standard errors below zero (Fisher-transformed) **and** well below the run's own per-batch distribution (median − 3 robust σ, measured once before anything is removed), **and** where the per-image CC to the merge — an independent channel, measured on the partials one frame at a time — also says the stretch agrees with the run worse than a typical frame does, because selecting frames by their disagreement with the merge and then reporting that the merge improved is circular; its edges are then slid frame by frame with the whole stretch re-measured at each position — so the range is located finely while the count of reflections the verdict rests on stays that of the stretch — the worst range is removed, the reference is re-formed and the scan repeats, never past a quarter of the sweep and never over a stretch narrower than one rocking event. The resulting per-frame `merged` / `downgraded` / `rejected` ledger is described in [the results report](RUGNUX_REPORT.md). ### 10.7 R-free test-set flags diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ead78ecef..312f781da 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -60,6 +60,14 @@ scaled-up counting error gives it, and is in `R_MEAS_WEIGHTED=` at that weight t the number with its noise; a large gap between the two says weak frames are kept at low weight. `MULTIPLICITY=` stays a count of observations. +**`CC_HALF=`** is the standard CC1/2, as XDS, Aimless and DIALS report it and as the mmCIF item +`_reflns.pdbx_CC_half` defines it: the observations of each reflection are dealt into two half-sets, +each half is averaged with 1/sigma^2 weights, and CC1/2 is the plain Pearson correlation of the two +half-set means, every reflection counted alike. The same holds for `REFRES_CC_HALF=`, the CC1/2 and CC\* +columns of every shell table and the mmCIF. The automatic resolution cut is fitted to a different, +internal form of CC1/2, in which a reflection measured only on frames the crystal barely diffracted on +counts less (see `FITTED_RESOLUTION`); on a sweep without a weak stretch the two agree. + **`R_MERGE=`** and **`R_PIM=`** are the other two merging R-factors over the same observations as `R_MEAS=`: R_merge without the multiplicity correction, and the precision-indicating R_pim (Weiss, 2001), which falls as the multiplicity rises. **`CC_STAR=`** is CC\* = sqrt(2 CC1/2 / (1 + CC1/2)) @@ -237,7 +245,8 @@ The keys most scripts read, by section; the rest of this page explains the ones | `MOSAICITY_DEG` | Rotation only ([below](#mosaicity)) | | `MERGE` | Section 3: `PERFORMED` or `NOT_PERFORMED` (`--no-merge`) | | `INCLUDE_RESOLUTION_RANGE` | `dmax dmin` of the written reflections (Å), as XDS's keyword | -| `FITTED_RESOLUTION` | Where the CC1/2 fall-off crosses the `--resolution-cc-target`; absent where the curve is not a fall-off (a `RESOLUTION_FIT` warning says why) | +| `FITTED_RESOLUTION` | Where the CC1/2 fall-off crosses the `--resolution-cc-target`; absent where the curve is not a fall-off (a `RESOLUTION_FIT` warning says why). It is read off the automatic cut's own fit, which weights each reflection by the scale of the frames it was measured on so that a weak stretch of the sweep does not drag it in; on such a sweep the unweighted `CC_HALF` of the shell table can cross the target earlier | +| `FITTED_RESOLUTION_PLAIN`, `RESOLUTION_CUT_PLAIN` | Diagnostic, rotation only: the same automatic fit made on the unweighted CC1/2 of the shell table - its crossing, and the limit the cut would have written from it. Nothing is cut at them; they show what the weighting of the fit is worth on this run | | `UNIQUE_REFLECTIONS`, `TOTAL_OBSERVATIONS`, `COMPLETENESS`, `MULTIPLICITY` | Over the written range; completeness in percent | | `I_OVER_SIGMA`, `R_MEAS`, `CC_HALF`, `WILSON_B` | The overall merging statistics; the shell table beside them is in the comment lines. `R_MEAS` (and `R_MERGE`, `R_PIM`) rotation only, `R_SPLIT` stills only | | `SIGANO`, `CC_ANOM`, `FRIEDELS_LAW` | The anomalous signal, and whether Friedel pairs were merged (`FALSE` with `-A`) | diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index ef01f58ca..124549ab7 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -945,6 +945,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m : 0.0; ss.cc_half = sa.cc_half.GetCC(); + ss.cc_half_internal = ss.cc_half; // stills weight nothing ss.cc_ref = sa.cc_ref.GetCC(); set_split(sa, ss); } @@ -972,6 +973,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m overall.unique_reflections = all_unique; overall.mean_i_over_sigma = n_i_over_sigma > 0 ? sum_i_over_sigma / n_i_over_sigma : 0.0; overall.cc_half = cc_half_overall.GetCC(); + overall.cc_half_internal = overall.cc_half; overall.cc_ref = cc_ref_overall.GetCC(); set_split(split_overall, overall); diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index b762577b8..84c7c82ea 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -29,7 +29,16 @@ struct MergeStatisticsShell { double mean_i_over_sigma = 0.0; + // The reported CC1/2: the plain Pearson correlation of the two half-set means (each mean itself + // 1/sigma^2-weighted over its observations), every reflection counted alike - the CC1/2 of XDS, + // Aimless and DIALS, and the one the mmCIF item pdbx_CC_half defines. double cc_half = 0.0f; + // Internal, not reported: the CC1/2 the processing DECIDES on. On rotation each reflection's pair + // enters with its MergedReflection::cc_weight, so a pair measured only on frames the crystal barely + // diffracted on does not count as much as a well-measured one; equal to cc_half where every weight + // is 1 (a sweep without a weak stretch), and on stills, which weight nothing. The automatic + // resolution cut is fitted to this same weighted form (ComputeCCHalfLogisticCutoff). + double cc_half_internal = 0.0f; double cc_ref = NAN; // Redundancy-independent merging R-factor (Diederichs & Karplus 1997), computed over the diff --git a/image_analysis/scale_merge/ResolutionCutoff.cpp b/image_analysis/scale_merge/ResolutionCutoff.cpp index 9f7b5aa62..6f4a7dd3f 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.cpp +++ b/image_analysis/scale_merge/ResolutionCutoff.cpp @@ -63,7 +63,7 @@ std::vector CCHalfFrameFactors(const std::vector &frame_scale, } ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, - double cc_target, Logger &logger) { + double cc_target, Logger &logger, bool weighted) { ResolutionCutoffResult result; if (!(cc_target > 0.0 && cc_target < 1.0)) { result.note = "invalid CC target"; @@ -130,7 +130,9 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector(m.h) * 73856093u ^ static_cast(m.k) * 19349663u ^ static_cast(m.l) * 83492791u; - group_cc[b][hash % bin_groups[b]].Add(m.I_half[0], m.I_half[1], m.cc_weight); + // Weighted by cc_weight: the cut is an internal decision and runs on the weighted CC1/2 + // (MergeStatisticsShell::cc_half_internal), not on the plain one the report prints. + group_cc[b][hash % bin_groups[b]].Add(m.I_half[0], m.I_half[1], weighted ? m.cc_weight : 1.0f); } // Usable bins (enough counts), in ascending-s order, with their bin-centre s. diff --git a/image_analysis/scale_merge/ResolutionCutoff.h b/image_analysis/scale_merge/ResolutionCutoff.h index 8252fbb1b..7c73aa658 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.h +++ b/image_analysis/scale_merge/ResolutionCutoff.h @@ -44,12 +44,15 @@ double TypicalFrameScale(const std::vector &frame_scale, const std::vect // The per-frame factor of the CC1/2 weight (MergedReflection::cc_weight): max(1, (G_ref/G)^2), the // variance an observation from a frame at scale G carries over the same observation at the typical -// scale G_ref (TypicalFrameScale). +// scale G_ref (TypicalFrameScale). Internal decisions only (the resolution cut, the batch delta-CC1/2 +// ledger); the reported CC1/2 is unweighted. std::vector CCHalfFrameFactors(const std::vector &frame_scale, const std::vector &n_obs); +// weighted: each reflection's pair enters with its cc_weight, the internal form the cut is decided on; +// false fits the plain (reported) CC1/2 instead - a diagnostic only (FITTED_RESOLUTION_PLAIN). ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, - double cc_target, Logger &logger); + double cc_target, Logger &logger, bool weighted = true); // Resolve the effective high-resolution limit and trim `merged` to it, in one place shared by the // stills merge, the rotation merge and the offline --scale path. A manual limit (manual_limit) wins; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 117f120a9..b0041a413 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3450,7 +3450,8 @@ void RotationScaleMerge::MeasureBatchDeltaCCHalf(int n_groups) { const std::vector frame_factor = CCHalfFrameFactors(frame_scale_fulls, n_obs_frame); ParallelChunks(n_full, ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) - obs_info[i] = 1.0 / frame_factor[obs_frame[i]]; // the information the observation carries, <= 1 + // the information the observation carries, <= 1 (the cc_weight factor; internal, not reported) + obs_info[i] = 1.0 / frame_factor[obs_frame[i]]; }); for (int f = 0; f < n_frames; ++f) fstart[f + 1] += fstart[f]; // The usable observations once more, in frame order and with the three fields a measurement reads @@ -5835,7 +5836,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // typical earns nothing (the factor is floored at 1), so on a sweep without a weak stretch // every weight is 1 and the CC1/2 is the plain Pearson it was. G is the frame's scale taken // before the correction surfaces (frame_scale_fulls), which are per-observation and - // resolution-dependent and are not a property of the frame. ---- + // resolution-dependent and are not a property of the frame. Internal, not reported: the weight + // enters the CC1/2 the processing decides on (cc_half_internal, the resolution cut); the + // reported CC1/2 is the plain Pearson of the half-set means. ---- // Counts, so the blocks may add theirs up in any split. std::vector n_obs_frame(frame_scale_fulls.size(), 0); { @@ -6125,6 +6128,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // _process.h5 is written elsewhere from the partials, so no data is lost. A manual // --scaling-high-resolution (d_min_limit) wins; the P1 search merge (for_search) is never cut, so // the space-group search still sees the full range. + // Diagnostic only: the same fit on the plain CC1/2, before the cut trims the reflections. It is + // reported beside the applied cut and decides nothing. + if (!d_min_limit && !for_search && resolution_cutoff_method == ResolutionCutoffMethod::CCHalfLogistic) { + const auto plain = ComputeCCHalfLogisticCutoff(result.merged, resolution_cc_target, *logger, + /*weighted=*/false); + result.resolution_fit_plain_A = plain.d_fit; + result.resolution_cut_plain_A = plain.d_cut; + logger->Info("Resolution fit on the plain CC1/2 (diagnostic, not applied): crossing {}, cut {}", + plain.d_fit ? fmt::format("{:.2f} A", *plain.d_fit) : "none", + plain.d_cut ? fmt::format("{:.2f} A", *plain.d_cut) : "none"); + } const std::optional effective_d_min = ApplyResolutionCutoff( result.merged, d_min_limit, resolution_cutoff_method, resolution_cc_target, for_search, *logger, &result.resolution_fit_A); @@ -6407,7 +6421,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; + CorrelationCoefficient cc_half; // plain Pearson: the reported CC1/2 + CorrelationCoefficient cc_half_internal; // cc_weight-weighted: internal decisions, not reported }; std::vector sa(n_shells); std::vector possible(n_shells, 0); @@ -6422,7 +6437,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; sub("st merged walk"); - CorrelationCoefficient cc_half_overall; + CorrelationCoefficient cc_half_overall, cc_half_internal_overall; for (const auto &m : merged) { const auto shell = shells.GetShell(m.d); if (!shell || *shell < 0 || *shell >= n_shells) continue; @@ -6431,8 +6446,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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], m.cc_weight); - cc_half_overall.Add(m.I_half[0], m.I_half[1], m.cc_weight); + s.cc_half.Add(m.I_half[0], m.I_half[1]); + cc_half_overall.Add(m.I_half[0], m.I_half[1]); + s.cc_half_internal.Add(m.I_half[0], m.I_half[1], m.cc_weight); + cc_half_internal_overall.Add(m.I_half[0], m.I_half[1], m.cc_weight); } } } @@ -6835,6 +6852,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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_half_internal = sa[s].cc_half_internal.GetCC(); ss.cc_ref = NAN; ss.r_meas = urmeas_den[s] > 0.0 ? urmeas_num[s] / urmeas_den[s] : NAN; ss.r_meas_weighted = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; @@ -6853,6 +6871,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } overall.mean_i_over_sigma = n_ios > 0 ? sum_ios / n_ios : 0.0; overall.cc_half = cc_half_overall.GetCC(); + overall.cc_half_internal = cc_half_internal_overall.GetCC(); overall.cc_ref = NAN; overall.r_meas = urmeas_den_all > 0.0 ? urmeas_num_all / urmeas_den_all : NAN; overall.r_meas_weighted = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; @@ -7734,7 +7753,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool if (measure_cc_before_corrections && corrections) { cc_half_before_corrections = MergeAndStats(n_groups, for_search, combined_on_gpu && scaled_fulls_on_gpu, - /*full_stats=*/false).statistics.overall.cc_half; + /*full_stats=*/false).statistics.overall.cc_half_internal; logger->Info("CC1/2 before the correction surfaces: {:.4f} (the merge just above; the pass this " "one is judged against fits no surfaces, so this is what the two are compared on)", *cc_half_before_corrections); @@ -7810,7 +7829,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // With the surfaces off (or none of them adopted) the merge that just ran IS the uncorrected one, so // its CC1/2 is the answer without a second merge. A caller that did not ask keeps the NaN. if (measure_cc_before_corrections) - r.cc_half_before_corrections = cc_half_before_corrections.value_or(r.statistics.overall.cc_half); + r.cc_half_before_corrections = cc_half_before_corrections.value_or(r.statistics.overall.cc_half_internal); // Each segment's scale: the median of its frames' scales, at the run's gauge. for (size_t k = 0; k < segments.size(); ++k) { std::vector g; diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 337c40728..fe9d668a2 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -119,17 +119,23 @@ public: double isa_high = std::numeric_limits::quiet_NaN(); double error_model_a = 0.0; // XDS convention: sigma^2 = a*(sigma0^2 + b*I^2) double error_model_b = 0.0; - // The overall CC1/2 as it stood BEFORE the correction surfaces were folded in. That is what the + // The overall CC1/2 (internal weighted form, cc_half_internal; not reported) as it stood BEFORE the + // correction surfaces were folded in. That is what the // two-pass quality guard compares one pass against the other with: a pass whose intensities are // discarded fits no surfaces (see full_stats), so judging the pass that does fit them by its // corrected CC1/2 would set two different measurements against each other. Equal to - // statistics.overall.cc_half whenever no surface was fitted, and NaN when the caller did not + // statistics.overall.cc_half_internal whenever no surface was fitted, and NaN when the caller did not // ask for it (see measure_cc_before_corrections) - it costs a merge, so it is not measured on // spec, and a caller that did not ask must not be handed a number that looks measured. double cc_half_before_corrections = std::numeric_limits::quiet_NaN(); // 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; + // Diagnostic, inert: the same automatic fit made on the plain (reported) CC1/2 instead of the + // weighted one the cut is decided on - its crossing and the limit it would have written. Nothing + // is cut at them. Empty where the automatic cut did not run or that fit found no crossing. + std::optional resolution_fit_plain_A; + std::optional resolution_cut_plain_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 diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 2889e0814..fd8ed5ec2 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -50,6 +50,10 @@ namespace { // 18: stills report their error model as ERROR_MODEL_A / _B / _C (C the partiality term) and no // ISA (nor ISa in the summary's Signal row); no R_MERGE / R_MEAS / R_PIM (nor REFRES_R_MEAS) // and no R_merge / R_meas / R_pim columns in their shell tables. + // 19: on rotation, CC_HALF, CC_STAR, REFRES_CC_HALF and the CC1/2 / CC* columns of the shell tables + // (and the mmCIF) are the standard unweighted Pearson of the half-set means, as XDS, Aimless and + // DIALS report it; before, each reflection was weighted by its frames' scale (cc_weight), which + // is now internal only (the automatic cut and FITTED_RESOLUTION). Stills were already unweighted. // 8: every line that is not `KEY= value` and not blank starts with `#`, warnings included // (`# WARNING:`); the header gains RUGNUX_DOWNLOAD and BUILD_CXX_FLAGS, and RUGNUX_GIT is // the full build-time hash. @@ -395,12 +399,16 @@ namespace { Add(s, KeyReal((prefix + "CC_HALF").c_str(), po.cc_half, "{:.4f}")); } + // The three helpers below judge the automatic cut's fit against the shells, so they read the curve the + // fit was made on - the internal weighted CC1/2 (cc_half_internal) - not the plain one the shell table + // prints. On a sweep without a weak stretch the two are the same number. + // The finest shell whose CC1/2 still reaches the cut's own target: where the table says the signal // is, for the completeness quoted beside the written range. double FinestSupportedShell(const MergeStatistics &ms, double target) { double finest = 0.0; for (const auto &sh : ms.shells) - if (std::isfinite(sh.cc_half) && sh.cc_half >= target && (finest == 0.0 || sh.d_min < finest)) + if (std::isfinite(sh.cc_half_internal) && sh.cc_half_internal >= target && (finest == 0.0 || sh.d_min < finest)) finest = sh.d_min; return finest; // 0.0 = no shell reaches the target at all } @@ -415,9 +423,9 @@ namespace { double CrossingShellLimit(const MergeStatistics &ms, double target) { double finest = 0.0; for (const auto &sh : ms.shells) { - if (!std::isfinite(sh.cc_half)) + if (!std::isfinite(sh.cc_half_internal)) continue; - if (sh.cc_half < target) + if (sh.cc_half_internal < target) return finest == 0.0 ? 0.0 : sh.d_min; finest = sh.d_min; } @@ -433,9 +441,9 @@ namespace { double CCHalfClimbsBackAboveTarget(const MergeStatistics &ms, double target) { bool below = false; for (const auto &sh : ms.shells) { - if (!std::isfinite(sh.cc_half)) + if (!std::isfinite(sh.cc_half_internal)) continue; - if (sh.cc_half < target) + if (sh.cc_half_internal < target) below = true; else if (below) return sh.d_min; @@ -1485,6 +1493,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (fitted_resolution) Add(s, KeyReal("FITTED_RESOLUTION", *fitted_resolution, "{:.2f}")); + // Diagnostic: the same automatic fit on the plain CC1/2 (rotation), its crossing and the limit + // it would have written. Not applied; printed as fitted, without the suppression above. + if (result.resolution_fit_plain_A) + Add(s, KeyReal("FITTED_RESOLUTION_PLAIN", *result.resolution_fit_plain_A, "{:.2f}")); + if (result.resolution_cut_plain_A) + Add(s, KeyReal("RESOLUTION_CUT_PLAIN", *result.resolution_cut_plain_A, "{:.2f}")); Add(s, KeyBool("FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel())); Add(s, KeyInt("UNIQUE_REFLECTIONS", o.unique_reflections)); @@ -1625,7 +1639,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " INCLUDE_RESOLUTION_RANGE is the range the reflections were WRITTEN to;" " FITTED_RESOLUTION\n is where the CC1/2 fall-off crosses {:.2f}, and is the number" " to quote. The data are kept\n one shell past it on purpose: an included shell can" - " still be downweighted by refinement,\n while a truncated one cannot be put back.", + " still be downweighted by refinement,\n while a truncated one cannot be put back." + " The fit weights each reflection by the scale of the\n frames it was measured on," + " so a weak stretch of the sweep does not drag it in; the CC1/2 in\n the shell table" + " is the standard unweighted one, and can read lower where there is such a stretch.", cc_target))); if (completeness_narrower && o.possible_unique_reflections > 0) Add(s, Prose(fmt::format( diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index ef411059c..42691107b 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -484,6 +484,7 @@ public: // The numbers of a search merge (ProcessResult::search_merge_*) the beam-centre arms are judged on // (MeasuredCentreWins), with the counts their noise is read from. +// The CC1/2 here are the internal weighted ones (MergeStatisticsShell::cc_half_internal), not the reported. struct SearchMergeSignal { int64_t strong = 0; // reflections at I/sigma >= 2 int64_t reflections = 0; @@ -581,6 +582,10 @@ struct ProcessResult { // Truncating cannot be undone and refinement can downweight, so the data are cut generously and this // is the number to quote. Empty when no automatic cut ran (a manual limit, or the method turned off). std::optional resolution_fit_A; + // Diagnostic (rotation): the same fit on the plain CC1/2 - crossing and the cut it would give; never + // applied (RotationScaleMerge::Result::resolution_fit_plain_A). + std::optional resolution_fit_plain_A; + std::optional resolution_cut_plain_A; // Detector tilt (PONI rot1/rot2, degrees) as rotation indexing MEASURED it, which is not what the // run integrated at - nothing writes a refined tilt back onto the experiment. Reported so the value // can be checked against a powder calibration or handed to another program; empty on stills and @@ -688,6 +693,7 @@ struct ProcessResult { // False when the merge had no cell to enumerate possible reflections against, so the guard is told // not to test a number it does not have and the log says so instead of printing a zero. bool search_merge_completeness_measured = false; + // Internal, not reported: the weighted CC1/2 the decisions read (MergeStatisticsShell::cc_half_internal). double search_merge_cc_half = NAN; // ... and in its lowest-resolution shell alone, where every reflection is strong and noise cannot // pull the number down: what reads low there is equivalents that disagree, i.e. a geometry that diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index e246a54f6..99964bf9e 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1819,6 +1819,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.error_model_isa_low = r.isa_low; result.error_model_isa_high = r.isa_high; if (r.resolution_fit_A) result.resolution_fit_A = r.resolution_fit_A; + if (r.resolution_fit_plain_A) result.resolution_fit_plain_A = r.resolution_fit_plain_A; + if (r.resolution_cut_plain_A) result.resolution_cut_plain_A = r.resolution_cut_plain_A; result.error_model_a = r.error_model_a; result.error_model_b = r.error_model_b; result.scaling_iterations_partials = r.scaling_iterations_partials; @@ -1884,7 +1886,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { out.merged = std::move(sm.merged); out.statistics = std::move(sm.statistics); out.reference_statistics = std::move(sm.reference_statistics); - out.cc_half_before_corrections = out.statistics.overall.cc_half; // no surfaces on this path + out.cc_half_before_corrections = out.statistics.overall.cc_half_internal; // no surfaces on this path logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label); Note(fmt::format("merge in {}: {} unique", label, out.merged.size())); return out; @@ -1985,7 +1987,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { ? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0; into.search_merge_cc_half = counted.cc_half_before_corrections; if (!counted.statistics.shells.empty()) { - into.search_merge_inner_cc_half = counted.statistics.shells.front().cc_half; + into.search_merge_inner_cc_half = counted.statistics.shells.front().cc_half_internal; into.search_merge_inner_reflections = counted.statistics.shells.front().unique_reflections; } // ... and the reflections it holds on the principal axial rows, which the guard needs @@ -3655,6 +3657,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const double em_a = result.error_model_a; const double em_b = result.error_model_b; const auto res_fit = result.resolution_fit_A; + const auto res_fit_plain = result.resolution_fit_plain_A; + const auto res_cut_plain = result.resolution_cut_plain_A; const auto manual_d_min = experiment_.GetScalingSettings().GetHighResolutionLimit_A(); phase("Merging under the refused point group " @@ -3726,6 +3730,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.error_model_a = em_a; result.error_model_b = em_b; result.resolution_fit_A = res_fit; + result.resolution_fit_plain_A = res_fit_plain; + result.resolution_cut_plain_A = res_cut_plain; } } } diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 91afc4693..cbaf1cf1b 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -396,9 +396,13 @@ TEST_CASE("ResolutionCutoff_ScaledUpNoiseDoesNotHideTheFallOff") { CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15)); // Counted as equals, the same reflections hide it: the curve reads noise from the first bin on. + // The plain fit (the FITTED_RESOLUTION_PLAIN diagnostic) is that, with the weights left in place. + const auto plain = ComputeCCHalfLogisticCutoff(merged, 0.30, logger, /*weighted=*/false); for (auto &m : merged) m.cc_weight = 1.0f; const auto unweighted = ComputeCCHalfLogisticCutoff(merged, 0.30, logger); CHECK((!unweighted.d_fit || *unweighted.d_fit > 3.0)); + CHECK(plain.d_fit == unweighted.d_fit); + CHECK(plain.d_cut == unweighted.d_cut); } // The frame factor of the CC1/2 weight is 1 on a sweep whose frames all sit at one scale, however that diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index dedae9003..344180b1c 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -734,7 +734,7 @@ TEST_CASE("ResultReport_UnusableMergeAndCompleteness", "[Diagnostics]") { MergeStatisticsShell shell; shell.d_max = 50.0f; shell.d_min = 2.0f; - shell.cc_half = 0.98; + shell.cc_half = shell.cc_half_internal = 0.98; shell.unique_reflections = 300; shell.possible_unique_reflections = 2000; thin.merge_statistics.shells.push_back(shell); @@ -763,13 +763,13 @@ TEST_CASE("ResultReport_FittedResolutionSuppressed", "[Diagnostics]") { MergeStatisticsShell inner; inner.d_max = 6.0f; inner.d_min = 0.90f; - inner.cc_half = 0.95; + inner.cc_half = inner.cc_half_internal = 0.95; inner.unique_reflections = 250; inner.possible_unique_reflections = 260; MergeStatisticsShell outer = inner; outer.d_max = 0.90f; outer.d_min = 0.80f; - outer.cc_half = 0.02; + outer.cc_half = outer.cc_half_internal = 0.02; result.merge_statistics.shells.push_back(inner); result.merge_statistics.shells.push_back(outer); // ...and a fit that claims 0.70 A, finer than any shell that still reaches the target. @@ -792,15 +792,28 @@ TEST_CASE("ResultReport_FittedResolutionSuppressed", "[Diagnostics]") { != std::string::npos); result.resolution_fit_A = 0.95; + // The fit is judged against the curve it was made on, the internal weighted CC1/2, not the plain one + // the table prints: a weak stretch can pull the plain shell CC1/2 under the target where the + // weighted one, and the fit, stay above it. + ProcessResult weak = result; + weak.merge_statistics.shells.front().cc_half = 0.25; + CHECK(RenderResultReport("p", "in.h5", x, weak).find("\nFITTED_RESOLUTION= 0.95\n") != std::string::npos); + // The fit on the plain CC1/2 is a diagnostic beside it, printed as fitted. + weak.resolution_fit_plain_A = 1.10; + weak.resolution_cut_plain_A = 1.02; + const auto weak_text = RenderResultReport("p", "in.h5", x, weak); + CHECK(weak_text.find("\nFITTED_RESOLUTION_PLAIN= 1.10\n") != std::string::npos); + CHECK(weak_text.find("\nRESOLUTION_CUT_PLAIN= 1.02\n") != std::string::npos); + // CC1/2 that drops through the target and climbs back above it is not a fall-off either, and the // finest shell reaching the target is then noise correlating with noise - so a fit the test above // accepts must still be suppressed, and the user told so in the default report. ProcessResult climbs = result; - climbs.merge_statistics.shells.back().cc_half = 0.05; + climbs.merge_statistics.shells.back().cc_half = climbs.merge_statistics.shells.back().cc_half_internal = 0.05; MergeStatisticsShell noise = climbs.merge_statistics.shells.back(); noise.d_max = 0.80f; noise.d_min = 0.70f; - noise.cc_half = 0.60; + noise.cc_half = noise.cc_half_internal = 0.60; climbs.merge_statistics.shells.push_back(noise); const auto climbs_text = RenderResultReport("p", "in.h5", x, climbs); CHECK(climbs_text.find("FITTED_RESOLUTION=") == std::string::npos); @@ -917,13 +930,13 @@ TEST_CASE("ResultReport_CompletenessTwoRanges", "[Diagnostics]") { MergeStatisticsShell inner; // signal here inner.d_max = 50.0f; inner.d_min = 1.50f; - inner.cc_half = 0.95; + inner.cc_half = inner.cc_half_internal = 0.95; inner.unique_reflections = 700; inner.possible_unique_reflections = 720; MergeStatisticsShell outer = inner; // written, but past the fall-off outer.d_max = 1.50f; outer.d_min = 1.20f; - outer.cc_half = 0.02; + outer.cc_half = outer.cc_half_internal = 0.02; outer.unique_reflections = 200; outer.possible_unique_reflections = 480; result.merge_statistics.shells.push_back(inner); @@ -937,7 +950,7 @@ TEST_CASE("ResultReport_CompletenessTwoRanges", "[Diagnostics]") { // Detector-limited: every shell still reaches the target, so the range where the signal is IS // the written range and the second number would repeat the first. - result.merge_statistics.shells[1].cc_half = 0.95; + result.merge_statistics.shells[1].cc_half = result.merge_statistics.shells[1].cc_half_internal = 0.95; result.merge_statistics.shells[1].possible_unique_reflections = 200; const auto same = RenderResultReport("p", "in.h5", x, result); CHECK(same.find(" % to ") == std::string::npos); diff --git a/tests/RotationScaleMergeSegmentsTest.cpp b/tests/RotationScaleMergeSegmentsTest.cpp index 92fa4325e..01f601b17 100644 --- a/tests/RotationScaleMergeSegmentsTest.cpp +++ b/tests/RotationScaleMergeSegmentsTest.cpp @@ -10,6 +10,7 @@ #include #include +#include "../common/CorrelationCoefficient.h" #include "../common/CUDAWrapper.h" #include "../common/Logger.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" @@ -244,3 +245,26 @@ TEST_CASE("RotationScaleMerge_NoFullsIsAnErrorNotACrash", "[RotationScale]") { x.ImportScalingSettings(s); CHECK_THROWS_AS(Merge(data, x, {}), JFJochException); } + +// The reported CC1/2 is the plain Pearson of the half-set means, every reflection counted alike, as +// XDS, Aimless and DIALS report it; the cc_weight-weighted form is internal (cc_half_internal). A sweep +// with a weak stretch makes the weights differ from 1, so the two are not the same number. +TEST_CASE("RotationScaleMerge_ReportedCCHalfIsThePlainPearson", "[RotationScale]") { + const std::vector sweeps{{60, 1.0, 0.0}, {60, 0.1, 0.0}}; + auto data = MakeSweeps(sweeps, 19); + const auto r = Merge(data, SyntheticExperiment(false), {}); + + CorrelationCoefficient plain, weighted; + int down_weighted = 0; + for (const auto &m: r.merged) + if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0f + && std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) { + plain.Add(m.I_half[0], m.I_half[1]); + weighted.Add(m.I_half[0], m.I_half[1], m.cc_weight); + if (m.cc_weight < 0.9f) ++down_weighted; + } + REQUIRE(down_weighted > 100); + CHECK(r.statistics.overall.cc_half == Catch::Approx(plain.GetCC()).epsilon(1e-9)); + CHECK(r.statistics.overall.cc_half_internal == Catch::Approx(weighted.GetCC()).epsilon(1e-9)); + CHECK(std::fabs(r.statistics.overall.cc_half - r.statistics.overall.cc_half_internal) > 1e-6); +}