From 94ed6aec045b1147064496a0a053c27cc2369b90 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 25 Sep 2026 04:44:50 +0200 Subject: [PATCH 1/2] Scaling: a frame below the credible floor does not hold the loop's step open ComputeSmoothGWindow leaves a frame below MIN_CREDIBLE_SCALE_RATIO of the median out of every window, and DropCollapsedScales drops it after the loop, but RunScalingLoop still counted its step. Such a frame has lost its vote in its own references and keeps moving, so the loop was stopped as "unsettled". On 8xtf (hq-pool battery) the first-pass partials loop stopped after 16 rounds at rms dlogG 1.8e-01 (base 9.5e-04); with the frame left out of the step it settles in 23 rounds (8.0e-04), and every other loop of that run that had stopped unsettled (pass 2 at 1.6e-03) now settles too. Every loop of that sweep that stopped unsettled was one that went on to drop a frame. The final 8xtf merge is bit-identical; the false SCALING_NOT_CONVERGED warning goes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- docs/CHANGELOG.md | 1 + image_analysis/scale_merge/RotationScaleMerge.cpp | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a2928f583..9504b1a74 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -7,6 +7,7 @@ * Rugnux rejects single observations that Wilson statistics make implausible (typically a hot pixel, zinger or ice spot), including pairs and single measurements that the equivalents test cannot judge; the count is reported as `OBSERVATIONS_REJECTED_WILSON`. * Rugnux reports the X-ray bandwidth it measures from spot shapes (for multilayer and pink beams); the value is reported only and does not change processing. * Rugnux accepts a detector-modulation or absorption correction surface when held-out data support it on Fisher's z, and fits them in the order modulation, time, goniometer frame; this recovers corrections that were wrongly refused. +* Rugnux no longer stops per-frame scaling as unsettled because of a single frame that is about to be dropped as blank. * `rugnux --polarization` is documented as the polarization degree (XDS `FRACTION_OF_POLARIZATION` = (1 + p)/2); the default 0.99 suits undulators. ### 1.0.0-rc.172 diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index b57d21126..e480d2aee 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1568,10 +1568,23 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunScalingLoop( } rescale(scaled_any, ratio); if (it > 0) { + // A frame already below the credible floor is not part of the fit's answer: the window + // leaves it out (ComputeSmoothGWindow) and the guard after the loop drops it. Left in the + // step, a single such frame held the rms at 0.18 and the loop was stopped as unsettled after + // 16 rounds; with it left out the same loop settles in 23. On that sweep every loop that + // stopped unsettled was one that went on to drop a frame. + std::vector fitted; + for (int f = 0; f < n_frames; ++f) + if (scaled_any[f] && g[f] > 0.0) fitted.push_back(g[f]); + double g_floor = 0.0; + if (!fitted.empty()) { + std::nth_element(fitted.begin(), fitted.begin() + fitted.size() / 2, fitted.end()); + g_floor = fitted[fitted.size() / 2] * MIN_CREDIBLE_SCALE_RATIO; + } double ss = 0.0, sw = 0.0; int n = 0; for (int f = 0; f < n_frames; ++f) - if (scaled_any[f] && g_prev[f] > 0.0 && g[f] > 0.0) { + if (scaled_any[f] && g_prev[f] > 0.0 && g[f] > 0.0 && g[f] >= g_floor) { const double d = std::log(g[f] / g_prev[f]); ss += frame_obs_count[f] * d * d; sw += frame_obs_count[f]; From 5c5eb00904468c9e912a53db6f3788e4253d3d5b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 25 Sep 2026 04:52:18 +0200 Subject: [PATCH 2/2] Merge statistics: R_meas weighted as the merge weights each observation The delta-CC1/2 disposition keeps a weak stretch in the merge (downgraded) wherever removing it would not raise CC1/2, and the merge then carries each of its observations at the small 1/sigma^2 its scaled-up counting error gives it. R_meas counted those observations at full weight, so the frames that add next to nothing to the intensities set the number. hq-pool battery: 8xtf kept 115 weak frames (scale 0.1-0.2 of the run's) with the same CC1/2, and CC_model per shell as rc173 with 84 frames rejected, and R_meas was 1.3-1.5x higher in every shell; 9w3y (33 deg rejected -> 0, every model metric better) and lyso_x10sa_strong read the same way. The disposition itself is not segmentation-dependent: conviction is on the batch grid, not the ledger ranges, and those sets' dispositions changed because the corrected data changed. R_meas now weights each observation by its merge weight v = 1/sigma^2 under the error model (corrected_sigma on the host, ModelSigma on the GPU, with the error model of the last MergeAccum), normalised per reflection to Kish's effective count, so equal sigmas give the ordinary formula (WeightedRmeasTerms). Where the proportional term dominates - strong reflections - frames are weighted alike, as in the merge: 8xtf's lowest shell reads 15.0%, the rc173 run with 84 frames rejected 15.2%. The per-hand table is weighted the same way. R_MEAS_UNWEIGHTED / REFRES_R_MEAS_UNWEIGHTED keep the XDS/AIMLESS convention and are what to set beside XDS; the battery scorer records them. MULTIPLICITY stays a count. Effect (weighted / unweighted): lyso_x06da_ref 0.0454 / 0.0479, 8xtf 0.225 / 0.907, insu_I_x06da_ref REFRES 0.072 / 0.232. GPU and CPU paths agree to the last printed digit on lyso_x06da_ref. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- docs/CHANGELOG.md | 1 + docs/RUGNUX_REPORT.md | 8 +++ image_analysis/scale_merge/Merge.cpp | 13 ++++ image_analysis/scale_merge/Merge.h | 15 +++- .../scale_merge/RotationScaleMerge.cpp | 68 ++++++++++++++----- .../scale_merge/RotationScaleMergeGPU.cu | 43 +++++++++--- .../scale_merge/RotationScaleMergeGPU.h | 11 +-- rugnux/ResultReport.cpp | 4 ++ tests/MergeScaleTest.cpp | 21 ++++++ tools/battery/score.py | 5 +- 10 files changed, 158 insertions(+), 31 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9504b1a74..0fd229f20 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -7,6 +7,7 @@ * Rugnux rejects single observations that Wilson statistics make implausible (typically a hot pixel, zinger or ice spot), including pairs and single measurements that the equivalents test cannot judge; the count is reported as `OBSERVATIONS_REJECTED_WILSON`. * Rugnux reports the X-ray bandwidth it measures from spot shapes (for multilayer and pink beams); the value is reported only and does not change processing. * Rugnux accepts a detector-modulation or absorption correction surface when held-out data support it on Fisher's z, and fits them in the order modulation, time, goniometer frame; this recovers corrections that were wrongly refused. +* Rugnux weights each observation in `R_MEAS` as the merge weights it, so weak frames kept in the merge no longer inflate it; the unweighted value is reported as `R_MEAS_UNWEIGHTED`. * Rugnux no longer stops per-frame scaling as unsettled because of a single frame that is about to be dropped as blank. * `rugnux --polarization` is documented as the polarization degree (XDS `FRACTION_OF_POLARIZATION` = (1 + p)/2); the default 0.99 suits undulators. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index f1b68b8c9..6f0fbd676 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -51,6 +51,14 @@ not `-A` was given, but where no Bijvoet pair could be split in both hands — a `-A`, or too few pairs — the quantity does not exist and the key is absent; `COMPLETENESS=`, `MULTIPLICITY=`, `I_OVER_SIGMA=`, `R_MEAS=`, `CC_HALF=` and `WILSON_B=` follow the same rule. +**`R_MEAS=` weights each observation as the merge weights it** - by 1/sigma^2 under the error model - +with each reflection's weights normalised to their effective number, so where a reflection's +observations share one sigma it is the ordinary R_meas. A stretch the crystal barely diffracted over is +in the merge at the small weight its scaled-up counting error gives it, and is in `R_MEAS=` at that +weight too, rather than setting the number with its noise. **`R_MEAS_UNWEIGHTED=`** +(and `REFRES_R_MEAS_UNWEIGHTED=`) is the same statistic with every observation counted once, as XDS +and AIMLESS count it - the one to set beside theirs. `MULTIPLICITY=` stays a count of observations. + The last blocks before `END OF REPORT` are the authorship and the acknowledgement: who wrote rugnux, its licence (GPLv3 — free to use for academic institutions and commercial companies alike) and where releases are published, then the credit to the X-ray research community whose methods diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index dea9f923f..4c63aa15c 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -756,6 +756,19 @@ const char *SweepQualityReasonText(SweepQualityReason reason) { return "unknown"; } +bool WeightedRmeasTerms(double sum_v_absdev, double sum_v_I, double sum_v, double sum_v2, + double &num, double &den) { + if (!(sum_v > 0.0) || !(sum_v2 > 0.0)) + return false; + const double n_eff = sum_v * sum_v / sum_v2; + if (!(n_eff > 1.0)) + return false; + const double scale = n_eff / sum_v; + num = std::sqrt(n_eff / (n_eff - 1.0)) * scale * sum_v_absdev; + den = scale * sum_v_I; + return true; +} + const char *FrameDispositionCode(FrameDisposition disposition) { switch (disposition) { case FrameDisposition::Merged: return "merged"; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index adb81c7ce..da17c9ce3 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -33,8 +33,13 @@ struct MergeStatisticsShell { double cc_ref = NAN; // Redundancy-independent merging R-factor (Diederichs & Karplus 1997), computed over the - // observations that enter the merge: R_meas = sum_hkl sqrt(n/(n-1)) sum_i|I_i-| / sum I_i. + // observations that enter the merge: R_meas = sum_hkl sqrt(n/(n-1)) sum_i|I_i-| / sum I_i. The + // rotation merge weights each observation as the merge weights it, 1/sigma^2 under the error model + // (WeightedRmeasTerms); with equal sigmas that is the formula as written. double r_meas = NAN; + // The same with every observation at weight one, as XDS and AIMLESS count it - for comparing with + // them. NaN where the two are the same by construction (the stills merge). + double r_meas_unweighted = NAN; // Anomalous signal-to-noise (XDS "SigAno" / mmCIF pdbx_absDiff_over_sigma_anomalous): // <|I(+)-I(-)|> / over acentric reflections measured in both hands. NaN when @@ -77,6 +82,14 @@ enum class FrameDisposition { const char *FrameDispositionCode(FrameDisposition disposition); // "merged", "downgraded", "rejected" +// One reflection's share of R_meas when its observations carry weights v_i: given sum v|I_i - |, +// sum v I_i, sum v and sum v^2, the numerator sqrt(n_eff/(n_eff-1)) sum v'|I_i - | and the +// denominator sum v' I_i, with the weights rescaled to v' = v n_eff / sum v so they add up to Kish's +// effective count n_eff = (sum v)^2 / sum v^2. Equal weights give the ordinary R_meas terms. False when +// the reflection has no second effective observation (n_eff <= 1), which then counts in neither sum. +bool WeightedRmeasTerms(double sum_v_absdev, double sum_v_I, double sum_v, double sum_v2, + double &num, double &den); + struct SweepQualityRange { // Inclusive, in processed-image ordinals - the numbering of _plot.txt and of every other // per-image array rugnux writes. With -s/--stride the source image is start + ordinal * stride. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e480d2aee..10a3580d3 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4768,21 +4768,37 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; + std::vector urmeas_num(n_shells, 0.0), urmeas_den(n_shells, 0.0); + double urmeas_num_all = 0.0, urmeas_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; }; + // R_meas: re-walk the fulls (Mask = cell only; no ice), accumulate |I_i - | per + // reflection - each observation weighted as the merge weights it, v = 1/sigma^2 under the + // error model (corrected_sigma). The merge carries an observation from a frame the crystal + // barely diffracted on at the small weight its scaled-up counting error gives it; counted at + // full weight in R_meas, that frame's noise set the number while adding next to nothing to + // the intensities it describes. Where the error model's proportional term dominates - the + // strong reflections - the merge weights every frame alike, and so does this. Measured: a + // sweep keeping 115 weak frames (scale 0.1-0.2 of the run's) read the same CC1/2, + // and CC to the model as the same sweep with 84 frames rejected, and R_meas 1.3-1.5x higher + // in every shell. Per reflection the weights are normalised to Kish's effective count + // n_eff = (sum v)^2 / sum v^2, which is n where the sigmas are equal, and there this is the + // ordinary R_meas = sum sqrt(n/(n-1)) sum|I_i - | / sum I_i. + // The unweighted sums are kept beside them for r_meas_unweighted. + struct RmeasObs { double sum_abs_dev = 0, sum_I = 0, w_abs_dev = 0, w_I = 0, sum_v = 0, sum_v2 = 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 rabsdev(n_groups), rsumI(n_groups), rwabsdev(n_groups), rwsumI(n_groups), + rsumv(n_groups), rsumv2(n_groups); std::vector rn(n_groups), rnusable(n_groups); - gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); + gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rwabsdev.data(), rwsumI.data(), + rsumv.data(), rsumv2.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); @@ -4796,7 +4812,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; + r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.w_abs_dev = rwabsdev[g]; r.w_I = rwsumI[g]; + r.sum_v = rsumv[g]; r.sum_v2 = rsumv2[g]; + r.n = rn[g]; r.shell = *shell; } } did_gpu_rmeas = true; @@ -4819,15 +4837,24 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; + const double sc = corrected_sigma(o, I_corr, sigma_corr); + const double v = 1.0 / (sc * sc); + const double dev = std::fabs(static_cast(I_corr) - merged_I[o.group]); + r.sum_abs_dev += dev; r.sum_I += I_corr; + r.w_abs_dev += v * dev; r.w_I += v * I_corr; r.sum_v += v; r.sum_v2 += v * v; + r.n++; r.shell = *shell; } } for (const auto &r : rmeas) { if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; + double num, den; + if (WeightedRmeasTerms(r.w_abs_dev, r.w_I, r.sum_v, r.sum_v2, num, den)) { + rmeas_num[r.shell] += num; rmeas_den[r.shell] += den; + rmeas_num_all += num; rmeas_den_all += den; + } 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; + urmeas_num[r.shell] += factor * r.sum_abs_dev; urmeas_den[r.shell] += r.sum_I; + urmeas_num_all += factor * r.sum_abs_dev; urmeas_den_all += r.sum_I; } } @@ -5000,7 +5027,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // observations the merge kept (anom_cell carries that filter), so an over-rejecting run // looks better here for the same reason it does in the table above. if (ha_measured) { - std::vector row_absdev(row_mean.size(), 0.0), row_sumI(row_mean.size(), 0.0); + // Weighted as the merge weights each observation, as the table above is. + std::vector row_absdev(row_mean.size(), 0.0), row_sumI(row_mean.size(), 0.0), + row_sumv(row_mean.size(), 0.0), row_sumv2(row_mean.size(), 0.0); std::vector row_n(row_mean.size(), 0); for (size_t i = 0; i < fulls.size(); ++i) { const int cell = anom_cell[i]; @@ -5008,14 +5037,18 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const int row = row_of_cell[cell]; if (!std::isfinite(row_mean[row])) continue; const float I_corr = fulls[i].I * fulls[i].corr; - row_absdev[row] += std::fabs(static_cast(I_corr) - row_mean[row]); - row_sumI[row] += I_corr; row_n[row]++; + const double sc = corrected_sigma(fulls[i], I_corr, fulls[i].sigma * fulls[i].corr); + const double v = 1.0 / (sc * sc); + row_absdev[row] += v * std::fabs(static_cast(I_corr) - row_mean[row]); + row_sumI[row] += v * I_corr; row_sumv[row] += v; row_sumv2[row] += v * v; row_n[row]++; } for (size_t row = 0; row < row_mean.size(); ++row) { if (row_n[row] < 2 || row_shell[row] < 0) continue; - const double factor = std::sqrt(static_cast(row_n[row]) / (row_n[row] - 1)); - ha[row_shell[row]].rmeas_num += factor * row_absdev[row]; - ha[row_shell[row]].rmeas_den += row_sumI[row]; + double num, den; + if (WeightedRmeasTerms(row_absdev[row], row_sumI[row], row_sumv[row], row_sumv2[row], num, den)) { + ha[row_shell[row]].rmeas_num += num; + ha[row_shell[row]].rmeas_den += den; + } } // The completeness denominator counted the same way: two per acentric, one per centric. if (reference_cell) { @@ -5041,6 +5074,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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.r_meas_unweighted = urmeas_den[s] > 0.0 ? urmeas_num[s] / urmeas_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(); } @@ -5056,6 +5090,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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.r_meas_unweighted = urmeas_den_all > 0.0 ? urmeas_num_all / urmeas_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(); @@ -5075,6 +5110,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool ? ha[sh].sum_i_over_sigma / ha[sh].n_i_over_sigma : 0.0; hs.cc_half = ha[sh].cc_half.GetCC(); hs.r_meas = ha[sh].rmeas_den > 0.0 ? ha[sh].rmeas_num / ha[sh].rmeas_den : NAN; + hs.r_meas_unweighted = NAN; out.per_hand.overall.total_observations += ha[sh].total_obs; out.per_hand.overall.unique_reflections += ha[sh].unique; out.per_hand.overall.possible_unique_reflections += ha[sh].possible; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index ccce57c9d..b4fcc4352 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -450,7 +450,7 @@ namespace { double *a_swI, *a_sw, *a_swIh0, *a_swIh1, *a_swh0, *a_swh1, *a_swht0, *a_swht1, *a_d; int32_t *a_nh0, *a_nh1, *a_rejected; uint8_t *a_on_ice; - double *r_absdev, *r_sumI; + double *r_absdev, *r_sumI, *r_wabsdev, *r_wsumI, *r_sumv, *r_sumv2; int32_t *r_n, *r_nusable; uint8_t *rejected_obs; // per-full: set by MergeAccum (outlier-rejected), read by MergeRmeas }; @@ -588,8 +588,10 @@ namespace { } } - // One thread per group: R_meas accumulators (sum|I_corr - merged_I|, sum_I, n) + the count of - // observations this looser walk accepted. Mirrors MergeAndStats' R_meas re-walk (cell-only filter). + // One thread per group: R_meas accumulators (sum|I_corr - merged_I|, sum I_corr, the same weighted by + // v = the observation's merge weight 1/sigma^2 under the error model (as MergeAccumKernel weights it), + // sum v, sum v^2, n) + the count of observations this looser walk accepted. + // Mirrors MergeAndStats' R_meas re-walk (cell-only filter). // That count is NOT the per-shell total_observations - it is wider than the merge, so it would // over-report multiplicity; the host uses it only to skip empty groups. __global__ void MergeRmeasKernel(MergeParams p) { @@ -597,7 +599,7 @@ namespace { const int lo = p.gstart[g], hi = lo + p.gcount[g]; const double mI = p.merged_I[g]; const bool have = isfinite(mI); - double absdev = 0, sumI = 0; + double absdev = 0, sumI = 0, wabsdev = 0, wsumI = 0, sumv = 0, sumv2 = 0; int n = 0, nusable = 0; for (int q = lo; q < hi; ++q) { const int i = p.gperm[q]; @@ -607,10 +609,19 @@ namespace { ++nusable; if (have) { const double I_corr = double(p.I[i]) * p.corr[i]; - absdev += fabs(I_corr - mI); sumI += I_corr; ++n; + const float sigma_raw = p.sigma[i] * p.corr[i]; + const double sc = p.error_model_active + ? ModelSigma(i, p, isfinite(p.em_mean[g]) ? p.em_mean[g] : I_corr, sigma_raw) + : sigma_raw; + const double v = 1.0 / (sc * sc); + const double dev = fabs(I_corr - mI); + absdev += dev; sumI += I_corr; + wabsdev += v * dev; wsumI += v * I_corr; sumv += v; sumv2 += v * v; ++n; } } - p.r_absdev[g] = absdev; p.r_sumI[g] = sumI; p.r_n[g] = n; p.r_nusable[g] = nusable; + p.r_absdev[g] = absdev; p.r_sumI[g] = sumI; p.r_wabsdev[g] = wabsdev; p.r_wsumI[g] = wsumI; + p.r_sumv[g] = sumv; p.r_sumv2[g] = sumv2; + p.r_n[g] = n; p.r_nusable[g] = nusable; } } @@ -670,10 +681,12 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr a_on_ice; // per group: any contributing full on an ice ring CudaDevicePtr reject_median; CudaDevicePtr reject_var_add; // the pooled test's widening - CudaDevicePtr merged_I, r_absdev, r_sumI; // R_meas per group (merged_I uploaded) + CudaDevicePtr merged_I, r_absdev, r_sumI, r_wabsdev, r_wsumI, r_sumv, r_sumv2; // R_meas per group (merged_I uploaded) CudaDevicePtr r_n, r_nusable; int merge_for_search = 0; // filter context for one MergeAndStats call double merge_min_part = 0.0; + double merge_em_a = 1.0, merge_em_b = 0.0; // the error model of the last MergeAccum, for MergeRmeas + int merge_em_active = 0; // combine: extra per-obs inputs + the one-time raw-hkl run layout CudaDevicePtr bkg, var_bkg, image_number, d_obs, px_obs, py_obs; @@ -954,6 +967,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ Upload(d.m_half, half, d.n_fulls); Upload(d.cc_factor, frame_cc_factor, d.n_frames); + d.merge_em_a = error_model_a; d.merge_em_b = error_model_b; d.merge_em_active = error_model_active ? 1 : 0; MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; @@ -998,13 +1012,16 @@ void RotationScaleMergeGPU::MergeAccumRange(int g0, int n, double *swI, double * dl(nh0, d.a_nh0); dl(nh1, d.a_nh1); dl(rejected, d.a_rejected); dl(on_ice_out, d.a_on_ice); } -void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, double *sumI, - int32_t *n, int32_t *nusable) { +void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, double *sumI, double *wabsdev, + double *wsumI, double *sumv, double *sumv2, int32_t *n, + int32_t *nusable) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int ng = d.n_groups; Upload(d.merged_I, merged_I, ng); d.r_absdev = CudaDevicePtr(std::max(1, ng)); d.r_sumI = CudaDevicePtr(std::max(1, ng)); + d.r_wabsdev = CudaDevicePtr(std::max(1, ng)); d.r_wsumI = CudaDevicePtr(std::max(1, ng)); + d.r_sumv = CudaDevicePtr(std::max(1, ng)); d.r_sumv2 = CudaDevicePtr(std::max(1, ng)); d.r_n = CudaDevicePtr(std::max(1, ng)); d.r_nusable = CudaDevicePtr(std::max(1, ng)); MergeParams p{}; @@ -1014,6 +1031,10 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.merged_I = d.merged_I.get(); p.r_absdev = d.r_absdev.get(); p.r_sumI = d.r_sumI.get(); p.r_n = d.r_n.get(); p.r_nusable = d.r_nusable.get(); + p.r_wabsdev = d.r_wabsdev.get(); p.r_wsumI = d.r_wsumI.get(); + p.r_sumv = d.r_sumv.get(); p.r_sumv2 = d.r_sumv2.get(); + p.error_model_a = d.merge_em_a; p.error_model_b = d.merge_em_b; p.error_model_active = d.merge_em_active; + p.em_mean = d.m_em_mean.get(); p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.rejected_obs = d.m_rejected.get(); const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK); @@ -1022,6 +1043,10 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d CudaCheck(cudaDeviceSynchronize(), "merge rmeas sync"); CudaCheck(cudaMemcpy(absdev, d.r_absdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl absdev"); CudaCheck(cudaMemcpy(sumI, d.r_sumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumI"); + CudaCheck(cudaMemcpy(wabsdev, d.r_wabsdev.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl wabsdev"); + CudaCheck(cudaMemcpy(wsumI, d.r_wsumI.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl wsumI"); + CudaCheck(cudaMemcpy(sumv, d.r_sumv.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumv"); + CudaCheck(cudaMemcpy(sumv2, d.r_sumv2.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl sumv2"); CudaCheck(cudaMemcpy(n, d.r_n.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rn"); CudaCheck(cudaMemcpy(nusable, d.r_nusable.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rnusable"); } diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 4aad3fce8..730a79add 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -114,10 +114,13 @@ public: int32_t *nh0, int32_t *nh1, double *d_out, int32_t *rejected, uint8_t *on_ice_out); - // Per-group R_meas accumulators (sum|I_corr-merged_I|, sum_I, n, and the count this looser walk - // accepted - which the host uses only to skip empty groups, NOT as the per-shell - // total_observations); merged_I is uploaded. All arrays length n_groups. - void MergeRmeas(const double *merged_I, double *absdev, double *sumI, int32_t *n, int32_t *nusable); + // Per-group R_meas accumulators: sum|I_corr-merged_I| and sum I, the same with each observation + // weighted by its merge weight v = 1/sigma^2 under the error model of the last MergeAccum, + // sum v, sum v^2, n, and the count this looser walk accepted - which the host uses only to skip + // empty groups, NOT as the per-shell total_observations; merged_I is uploaded. All arrays length + // n_groups. + void MergeRmeas(const double *merged_I, double *absdev, double *sumI, double *wabsdev, double *wsumI, + double *sumv, double *sumv2, int32_t *n, int32_t *nusable); // Post-smooth per-frame diagnostic CC: recompute the group means from the resident (smoothed) corr // and the Pearson CC of each frame's I*corr vs its group mean, downloading only the per-frame cc / diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 5c2b47175..246f889e8 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -803,6 +803,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("I_OVER_SIGMA", o.mean_i_over_sigma, "{:.2f}")); if (std::isfinite(o.r_meas)) Add(s, KeyReal("R_MEAS", o.r_meas, "{:.4f}")); + if (std::isfinite(o.r_meas_unweighted)) + Add(s, KeyReal("R_MEAS_UNWEIGHTED", o.r_meas_unweighted, "{:.4f}")); if (std::isfinite(o.cc_half)) Add(s, KeyReal("CC_HALF", o.cc_half, "{:.4f}")); if (std::isfinite(o.abs_diff_over_sigma_anomalous)) @@ -937,6 +939,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (std::isfinite(ro.r_meas)) Add(s, KeyReal("REFRES_R_MEAS", ro.r_meas, "{:.4f}")); + if (std::isfinite(ro.r_meas_unweighted)) + Add(s, KeyReal("REFRES_R_MEAS_UNWEIGHTED", ro.r_meas_unweighted, "{:.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)) diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 7be51a64e..408bce534 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -388,6 +388,27 @@ TEST_CASE("CCHalfFrameFactors_DeadStretchDoesNotSetTheTypicalScale") { CHECK(factor[100] == Catch::Approx(1.0).epsilon(1e-3)); } +// Equal weights give the ordinary R_meas terms, whatever their common value. +TEST_CASE("WeightedRmeas_EqualWeightsAreTheOrdinaryRmeas") { + // Three observations 9, 10, 11 of a reflection with = 10: sum|dev| = 2, sum I = 30. + double num, den; + REQUIRE(WeightedRmeasTerms(0.5 * 2.0, 0.5 * 30.0, 0.5 * 3, 0.25 * 3, num, den)); + CHECK(num == Catch::Approx(std::sqrt(3.0 / 2.0) * 2.0)); + CHECK(den == Catch::Approx(30.0)); +} + +// An observation that carries almost no information counts for almost nothing: the terms tend to those +// of the reflection without it, and a reflection left with one effective observation has none. +TEST_CASE("WeightedRmeas_NegligibleWeightDropsOut") { + // 9 and 11 at weight 1, = 10; a third at 40 with weight 1e-6. + const double v = 1e-6; + double num, den; + REQUIRE(WeightedRmeasTerms(2.0 + v * 30.0, 20.0 + v * 40.0, 2.0 + v, 2.0 + v * v, num, den)); + CHECK(num == Catch::Approx(std::sqrt(2.0) * 2.0).epsilon(1e-4)); + CHECK(den == Catch::Approx(20.0).epsilon(1e-4)); + CHECK_FALSE(WeightedRmeasTerms(0.0, 10.0, 1.0, 1.0, num, den)); +} + // A fall-off region a logistic cannot follow: CC1/2 drops through the target and comes straight back // up. The fitted crossing is then an extrapolation far past the bins it was made over, and reading // the cut off it writes the data deep into the noise; the crossing the bins themselves show is where diff --git a/tools/battery/score.py b/tools/battery/score.py index c44b3bb40..bc9d4d1fc 100644 --- a/tools/battery/score.py +++ b/tools/battery/score.py @@ -52,7 +52,7 @@ def read_report(path): # REFRES_* keys of the report read as numbers, stored lower-case under the same name REFRES_KEYS = ("refres_shells_past_limit", "refres_unique_reflections", "refres_completeness", "refres_multiplicity", "refres_i_over_sigma", "refres_r_meas", "refres_cc_half", - "refres_isa") + "refres_isa", "refres_r_meas_unweighted") def fnum(s): @@ -186,6 +186,9 @@ def judge(entry, rep, run_note): "d_min": None, "d_min_ref": ref.get("dmin"), "d_min_ref_rule": ref.get("dmin_rule"), "d_min_xds": ref.get("dmin_xds"), "res_gain_pct": None, "r_meas": fnum(rep.get("R_MEAS")), "cc_half": fnum(rep.get("CC_HALF")), + # R_MEAS weights each observation by its frame's information, as the merge does; the + # unweighted one is how XDS counts it (absent from reports before the weighting) + "r_meas_unweighted": fnum(rep.get("R_MEAS_UNWEIGHTED")), "isa": fnum(rep.get("ISA")), "completeness": fnum(rep.get("COMPLETENESS")), "multiplicity": fnum(rep.get("MULTIPLICITY")), "i_over_sigma": fnum(rep.get("I_OVER_SIGMA")), "indexing_rate": fnum(rep.get("INDEXING_RATE")), "images": fnum(rep.get("IMAGES_PROCESSED")),