From 0a342fcb134876347dff1d59780e12b4fe788ae0 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 19 Sep 2026 10:34:48 +0200 Subject: [PATCH] rugnux: weight every CC1/2 by the information a reflection carries The resolution cutoff, the shell table and the overall CC1/2 counted every unique reflection equally. The merge itself is inverse-variance weighted, so an observation from a frame the crystal barely diffracted on enters it at 1/G^2 of a good one - honestly, with its sigma - but a reflection measured only on such frames is scaled-up noise that then counts as much as a well-measured pair in every Pearson CC1/2 read off the merge. On a sweep where half the frames are weak the curve collapses at every resolution: the cut lands at 3.7 A on a P1 crystal whose good frames reach 1.5 A, and the UNUSABLE verdict fires (CC1/2 0.40 beside I/sigma 10). Each merged reflection now carries cc_weight: the precision its half-sets would have had with every observation at the run's typical frame scale, over the precision they have. G_ref = sum G^3 / sum G^2 over the usable observations is the precision-weighted typical scale, which the dead frames cannot drag down however many there are; the factor per observation is max(1, (G_ref/G)^2), with G the frame's total scale (partial scale, flux and the fulls' own G) taken before the correction surfaces and before collapsed frames are dropped, so nothing intensity- or resolution-dependent enters it. On a sweep without a weak stretch every weight is 1 and the CC1/2 is the plain Pearson it was. The cutoff fit, the shell table and the overall CC1/2 (and so the UNUSABLE verdict and the report's shell checks) all read the same weighted statistic. The two extra per-group sums are accumulated on both merge paths, the host loop and MergeAccumKernel, from one per-frame factor array; a host recompute of the device sums agrees to 1e-15 relative on every merge after the corrected corr is uploaded, and a host-merge run gives the same cut, space group, CC1/2 and ISa on five sets. Weighting by the half-set error variance alone (1/(v0+v1)) is not this: v grows with the intensity, so it weights the weak end of the intensity distribution and biases homogeneous data coarser. Measured (written resolution, together with the weighted outlier median): a P1 sweep with a long weak stretch 3.73 -> 1.40 A against a 1.63 A XDS reference, UNUSABLE withdrawn (overall CC1/2 0.40 -> 0.97); a second 3.61 -> 3.00 A; one with most of the sweep out of beam 7.35 -> 5.20 A. Lysozyme, thaumatin and two insulin sets unchanged to 0.01 A (one lysozyme sweep with a weak wedge 1.13 -> 1.16 A, from the median), same space groups throughout. Co-Authored-By: Claude Opus 5 (1M context) --- common/CorrelationCoefficient.cpp | 17 +++--- common/CorrelationCoefficient.h | 5 +- common/Reflection.h | 5 ++ .../scale_merge/ResolutionCutoff.cpp | 22 +++++++- image_analysis/scale_merge/ResolutionCutoff.h | 8 +++ .../scale_merge/RotationScaleMerge.cpp | 39 ++++++++++++- .../scale_merge/RotationScaleMerge.h | 10 +++- .../scale_merge/RotationScaleMergeGPU.cu | 31 +++++++---- .../scale_merge/RotationScaleMergeGPU.h | 8 ++- tests/MergeScaleTest.cpp | 55 ++++++++++++++++++- 10 files changed, 168 insertions(+), 32 deletions(-) diff --git a/common/CorrelationCoefficient.cpp b/common/CorrelationCoefficient.cpp index 24925571c..7c21dbb26 100644 --- a/common/CorrelationCoefficient.cpp +++ b/common/CorrelationCoefficient.cpp @@ -6,19 +6,20 @@ #include "CorrelationCoefficient.h" -void CorrelationCoefficient::Add(double x, double y) { - sum_x += x; - sum_y += y; - sum_x2 += x * x; - sum_y2 += y * y; - sum_xy += x * y; +void CorrelationCoefficient::Add(double x, double y, double w) { + sum_x += w * x; + sum_y += w * y; + sum_x2 += w * x * x; + sum_y2 += w * y * y; + sum_xy += w * x * y; + sum_w += w; n_cc++; } double CorrelationCoefficient::GetCC() const { - if (n_cc < 2) + if (n_cc < 2 || !(sum_w > 0.0)) return NAN; - const double n = static_cast(n_cc); + const double n = sum_w; const double mean_x = sum_x / n; const double mean_y = sum_y / n; diff --git a/common/CorrelationCoefficient.h b/common/CorrelationCoefficient.h index 711284226..943a3936a 100644 --- a/common/CorrelationCoefficient.h +++ b/common/CorrelationCoefficient.h @@ -9,9 +9,12 @@ class CorrelationCoefficient { double sum_x2 = 0.0; double sum_y2 = 0.0; double sum_xy = 0.0; + double sum_w = 0.0; int n_cc = 0; public: - void Add(double x, double y); + void Add(double x, double y) { Add(x, y, 1.0); } + // Weighted pair: a plain Pearson when every weight is the same. + void Add(double x, double y, double w); [[nodiscard]] double GetCC() const; }; diff --git a/common/Reflection.h b/common/Reflection.h index 5e4f86321..2e304cd79 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -64,6 +64,11 @@ struct MergedReflection { float sigma = NAN; float I_half[2] = {NAN, NAN}; float sigma_half[2] = {NAN, NAN}; + // Weight of this reflection in a CC1/2: the precision its half-sets would have had with every + // observation at the run's typical frame scale, over the precision they have. 1 on a sweep + // without a weak stretch; small for a reflection measured only where the crystal barely + // diffracted, whose scaled-up noise would otherwise count as much as a well-measured pair. + 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 // written - the ring contaminates it, it does not make it absent - and this only marks it so diff --git a/image_analysis/scale_merge/ResolutionCutoff.cpp b/image_analysis/scale_merge/ResolutionCutoff.cpp index 071a32b56..236a0323c 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.cpp +++ b/image_analysis/scale_merge/ResolutionCutoff.cpp @@ -33,6 +33,26 @@ namespace { } } +std::vector CCHalfFrameFactors(const std::vector &frame_scale, + const std::vector &n_obs) { + auto g_of = [&](size_t f) { + return std::isfinite(frame_scale[f]) && frame_scale[f] > 0.0 ? frame_scale[f] : 1.0; + }; + double s3 = 0.0, s2 = 0.0; + for (size_t f = 0; f < frame_scale.size(); ++f) { + const double g = g_of(f); + s3 += n_obs[f] * g * g * g; + s2 += n_obs[f] * g * g; + } + const double g_ref = s2 > 0.0 ? s3 / s2 : 1.0; + std::vector factor(frame_scale.size()); + for (size_t f = 0; f < frame_scale.size(); ++f) { + const double r = g_ref / g_of(f); + factor[f] = std::max(1.0, r * r); + } + return factor; +} + ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, double cc_target, Logger &logger) { ResolutionCutoffResult result; @@ -82,7 +102,7 @@ ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector(m.d) * m.d); int b = static_cast((s - s_lo) / bin_w); b = std::clamp(b, 0, N_FIT_BINS - 1); - bin_cc[b].Add(m.I_half[0], m.I_half[1]); + bin_cc[b].Add(m.I_half[0], m.I_half[1], m.cc_weight); ++bin_n[b]; } diff --git a/image_analysis/scale_merge/ResolutionCutoff.h b/image_analysis/scale_merge/ResolutionCutoff.h index b6913aea3..730013d97 100644 --- a/image_analysis/scale_merge/ResolutionCutoff.h +++ b/image_analysis/scale_merge/ResolutionCutoff.h @@ -34,6 +34,14 @@ struct ResolutionCutoffResult { std::string note; // human-readable description of the decision (for logging) }; +// 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 = sum G^3 / sum G^2 over the observations (n_obs[f] of them on frame f). Precision- +// weighted, so frames the crystal barely diffracted on cannot drag the typical scale down however many +// of them there are. A frame without a finite positive scale counts as G = 1. +std::vector CCHalfFrameFactors(const std::vector &frame_scale, + const std::vector &n_obs); + ResolutionCutoffResult ComputeCCHalfLogisticCutoff(const std::vector &merged, double cc_target, Logger &logger); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7bb72e504..7457c5c05 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -3819,6 +3819,22 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // above, on all threads, so a second serial walk over the fat records only rebuilds it. AssignHalvesByRank(obs_key, mf.group, n_groups, merge_half, nthreads); + // ---- CC1/2 weight (MergedReflection::cc_weight): the precision each half-set WOULD have had with + // every observation at the run's typical frame scale, over the precision it has. An observation + // from a frame the crystal barely diffracted on was scaled up by G_ref/G, so it carries + // (G_ref/G)^2 more counting variance than the same observation at the typical scale; that + // factor, and nothing intensity-dependent, is what the pair is discounted by. G_ref = + // sum G^3 / sum G^2 over the usable observations: the precision-weighted typical scale, which + // the dead frames cannot drag down however many of them there are. A frame brighter than + // 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. ---- + std::vector n_obs_frame(frame_scale_fulls.size(), 0); + for (int i = 0; i < n_full; ++i) + if (mf.group[i] >= 0) ++n_obs_frame[fulls[i].frame]; + const std::vector frame_cc_factor = CCHalfFrameFactors(frame_scale_fulls, n_obs_frame); + // ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets, then // export. In a lambda because the automatic resolution cutoff below re-runs it once the error // model has been refitted on the reflections that survive the cut. ---- @@ -3852,8 +3868,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool ms.Resize(n_groups); gpu_->MergeAccum(error_model_a, error_model_b, error_model_active, reject_outliers, reject_nsigma, reject_median.data(), merge_half.data(), + frame_cc_factor.data(), ms.swI.data(), ms.sw.data(), ms.swIh0.data(), ms.swIh1.data(), - ms.swh0.data(), ms.swh1.data(), ms.nh0.data(), ms.nh1.data(), + ms.swh0.data(), ms.swh1.data(), ms.swh_typ0.data(), ms.swh_typ1.data(), + ms.nh0.data(), ms.nh1.data(), ms.d.data(), ms.rej.data(), ms.on_ice.data(), rejected_obs.data()); // Per-group and independent, so the unpack splits; the rejected count is an integer sum and // does not care in which order the chunks add theirs in. @@ -3864,6 +3882,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool Accum &a = acc[g]; a.swI = ms.swI[g]; a.sw = ms.sw[g]; a.swIh[0] = ms.swIh0[g]; a.swIh[1] = ms.swIh1[g]; a.swh[0] = ms.swh0[g]; a.swh[1] = ms.swh1[g]; + a.swh_typ[0] = ms.swh_typ0[g]; a.swh_typ[1] = ms.swh_typ1[g]; a.nh[0] = static_cast(ms.nh0[g]); a.nh[1] = static_cast(ms.nh1[g]); a.d = static_cast(ms.d[g]); a.on_ice = ms.on_ice[g] != 0; @@ -3893,6 +3912,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const int half = merge_half[&o - fulls.data()]; a.swI += wI; a.sw += w; a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++; + a.swh_typ[half] += w * frame_cc_factor[o.frame]; // Any contributing observation on an ice ring marks the merged reflection (as // MergeOnTheFly does on the stills path): the intensity is still merged and written, // the flag only lets a consumer that must not read ring signal as crystal signal - @@ -3922,6 +3942,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool mr.I_half[i] = static_cast(a.swIh[i] / a.swh[i]); mr.sigma_half[i] = static_cast(1.0 / std::sqrt(a.swh[i])); } + mr.cc_weight = static_cast((1.0 / a.swh_typ[0] + 1.0 / a.swh_typ[1]) + / (1.0 / a.swh[0] + 1.0 / a.swh[1])); } if (!std::isfinite(a.d) || a.d <= 0.0f) continue; merged_I[g] = mr.I; @@ -4142,8 +4164,8 @@ 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]); - cc_half_overall.Add(m.I_half[0], 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); } } } @@ -4792,6 +4814,17 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } logger.Info("Scaled fulls (XDS order, Unity model)"); } + // The total scale each frame's fulls were amplified by - partial scale, incident flux and the fulls' + // own G (corr = 1/G here) - taken BEFORE the correction surfaces fold per-observation terms into + // corr, and before a collapsed frame's fulls are zeroed: the resident fulls only get that zero with + // the corrected corr, so a merge on the device before then still sees them, and must weigh them by + // their own scale. Read by MergeAndStats for the CC1/2 weight. + frame_scale_fulls.assign(n_frames, NAN); + for (const auto &o : fulls) + if (std::isfinite(o.corr) && o.corr > 0.0f && !std::isfinite(frame_scale_fulls[o.frame])) { + const double gp = g_partial[o.frame] * frame_flux[o.frame]; + frame_scale_fulls[o.frame] = (std::isfinite(gp) && gp > 0.0 ? gp : 1.0) / o.corr; + } const bool rejected_full_scales = scale_fulls && DropCollapsedFullScales(combined_on_gpu && scaled_fulls_on_gpu); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 24a1fed00..3855f1a0f 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -273,6 +273,7 @@ private: std::vector fulls; // combined fulls (rebuilt each Run), sorted by frame std::vector fulls_frame_start, fulls_frame_count; // CSR ranges of `fulls` per frame std::vector g_full; // per-frame scale on the fulls + std::vector frame_scale_fulls; // per-frame total scale of the fulls before the surfaces (cc_weight) // The five fields the host-side walks over the fulls actually read, pulled out of the 80-byte record // once per merge (see MergeAndStats). `group` carries the usability decision: -1 means the full is @@ -306,15 +307,18 @@ private: // The merge accumulators (see MergeAndStats' run_merge): one entry per ASU group, plus the arrays // the device kernel fills that the host unpacks into them. Members for the same reason as // FullsStaging - a merge runs several times per run and this is a hundred megabytes between them. - struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}; size_t nh[2] = {0, 0}; float d = NAN; - bool on_ice = false; }; + // swh_typ: each half-set's precision had every observation been at the run's typical frame scale + // (see MergeAndStats, cc_weight). + struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}, swh_typ[2] = {0, 0}; + size_t nh[2] = {0, 0}; float d = NAN; bool on_ice = false; }; std::vector merge_acc; struct MergeAccumStaging { - std::vector swI, sw, swIh0, swIh1, swh0, swh1, d; + std::vector swI, sw, swIh0, swIh1, swh0, swh1, swh_typ0, swh_typ1, d; std::vector nh0, nh1, rej; std::vector on_ice; void Resize(int n) { swI.resize(n); sw.resize(n); swIh0.resize(n); swIh1.resize(n); swh0.resize(n); swh1.resize(n); + swh_typ0.resize(n); swh_typ1.resize(n); d.resize(n); nh0.resize(n); nh1.resize(n); rej.resize(n); on_ice.resize(n); } }; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 40e146dc5..2f51def78 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -464,13 +464,13 @@ namespace { const int32_t *group, *frame; const uint8_t *on_ice, *frame_cell_ok, *half; const int32_t *gperm, *gstart, *gcount; - const double *em_mean, *merged_I; + const double *em_mean, *merged_I, *cc_factor; // outputs double *sw, *swI, *em_mean_out; int32_t *cnt; double *s2, *I2, *dev2; uint8_t *valid; - double *a_swI, *a_sw, *a_swIh0, *a_swIh1, *a_swh0, *a_swh1, *a_d; + 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; @@ -544,7 +544,7 @@ namespace { __global__ void MergeAccumKernel(MergeParams p) { for (int g = blockIdx.x * blockDim.x + threadIdx.x; g < p.n_groups; g += gridDim.x * blockDim.x) { const int lo = p.gstart[g], hi = lo + p.gcount[g]; - double swI = 0, sw = 0, swIh0 = 0, swIh1 = 0, swh0 = 0, swh1 = 0, dd = NAN; + double swI = 0, sw = 0, swIh0 = 0, swIh1 = 0, swh0 = 0, swh1 = 0, swht0 = 0, swht1 = 0, dd = NAN; int nh0 = 0, nh1 = 0, rejected = 0, onice = 0; const float rmed = p.reject_median[g]; for (int q = lo; q < hi; ++q) { @@ -574,13 +574,16 @@ namespace { // a rank within the reflection cannot be read off a running count, and computing it // in one place is what makes this kernel and the host merge loop agree exactly. const int half = p.half[i]; + const double wt = w * p.cc_factor[p.frame[i]]; // see the host merge loop (cc_weight) swI += wI; sw += w; - if (half) { swIh1 += wI; swh1 += w; ++nh1; } else { swIh0 += wI; swh0 += w; ++nh0; } + if (half) { swIh1 += wI; swh1 += w; swht1 += wt; ++nh1; } + else { swIh0 += wI; swh0 += w; swht0 += wt; ++nh0; } if (p.on_ice[i]) onice = 1; // see the host merge loop if (!isfinite(dd) && isfinite(p.d[i]) && p.d[i] > 0.0f) dd = p.d[i]; } p.a_swI[g] = swI; p.a_sw[g] = sw; p.a_swIh0[g] = swIh0; p.a_swIh1[g] = swIh1; - p.a_swh0[g] = swh0; p.a_swh1[g] = swh1; p.a_nh0[g] = nh0; p.a_nh1[g] = nh1; p.a_d[g] = dd; + p.a_swh0[g] = swh0; p.a_swh1[g] = swh1; p.a_swht0[g] = swht0; p.a_swht1[g] = swht1; + p.a_nh0[g] = nh0; p.a_nh1[g] = nh1; p.a_d[g] = dd; p.a_rejected[g] = rejected; p.a_on_ice[g] = uint8_t(onice); } } @@ -657,7 +660,8 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr m_valid; CudaDevicePtr m_rejected; // per-full outlier-rejected flag (MergeAccum -> MergeRmeas) CudaDevicePtr m_half; // per-full CC1/2 half-set, assigned on the host - CudaDevicePtr a_swI, a_sw, a_swIh0, a_swIh1, a_swh0, a_swh1, a_d; // merge accum per group + CudaDevicePtr a_swI, a_sw, a_swIh0, a_swIh1, a_swh0, a_swh1, a_swht0, a_swht1, a_d; // merge accum per group + CudaDevicePtr cc_factor; // per-frame CC1/2 weight factor (length n_frames) CudaDevicePtr a_nh0, a_nh1, a_rejected; CudaDevicePtr a_on_ice; // per group: any contributing full on an ice ring CudaDevicePtr reject_median; @@ -902,9 +906,10 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_b, bool error_model_active, bool reject_outliers, double reject_nsigma, const float *reject_median, - const uint8_t *half, + const uint8_t *half, const double *frame_cc_factor, double *swI, double *sw, double *swIh0, double *swIh1, - double *swh0, double *swh1, int32_t *nh0, int32_t *nh1, double *d_out, + double *swh0, double *swh1, double *swh_typ0, double *swh_typ1, + int32_t *nh0, int32_t *nh1, double *d_out, int32_t *rejected, uint8_t *on_ice_out, uint8_t *rejected_obs) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; @@ -912,6 +917,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ d.a_swI = CudaDevicePtr(std::max(1, ng)); d.a_sw = CudaDevicePtr(std::max(1, ng)); d.a_swIh0 = CudaDevicePtr(std::max(1, ng)); d.a_swIh1 = CudaDevicePtr(std::max(1, ng)); d.a_swh0 = CudaDevicePtr(std::max(1, ng)); d.a_swh1 = CudaDevicePtr(std::max(1, ng)); + d.a_swht0 = CudaDevicePtr(std::max(1, ng)); d.a_swht1 = CudaDevicePtr(std::max(1, ng)); d.a_nh0 = CudaDevicePtr(std::max(1, ng)); d.a_nh1 = CudaDevicePtr(std::max(1, ng)); d.a_d = CudaDevicePtr(std::max(1, ng)); d.a_rejected = CudaDevicePtr(std::max(1, ng)); d.a_on_ice = CudaDevicePtr(std::max(1, ng)); @@ -920,6 +926,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ "zero m_rejected"); Upload(d.reject_median, reject_median, ng); Upload(d.m_half, half, d.n_fulls); + Upload(d.cc_factor, frame_cc_factor, d.n_frames); MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; @@ -934,9 +941,10 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.half = d.m_half.get(); p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); - p.em_mean = d.m_em_mean.get(); + p.em_mean = d.m_em_mean.get(); p.cc_factor = d.cc_factor.get(); p.a_swI = d.a_swI.get(); p.a_sw = d.a_sw.get(); p.a_swIh0 = d.a_swIh0.get(); p.a_swIh1 = d.a_swIh1.get(); - p.a_swh0 = d.a_swh0.get(); p.a_swh1 = d.a_swh1.get(); p.a_nh0 = d.a_nh0.get(); p.a_nh1 = d.a_nh1.get(); + p.a_swh0 = d.a_swh0.get(); p.a_swh1 = d.a_swh1.get(); + p.a_swht0 = d.a_swht0.get(); p.a_swht1 = d.a_swht1.get(); p.a_nh0 = d.a_nh0.get(); p.a_nh1 = d.a_nh1.get(); p.a_d = d.a_d.get(); p.a_rejected = d.a_rejected.get(); p.a_on_ice = d.a_on_ice.get(); p.rejected_obs = d.m_rejected.get(); @@ -947,7 +955,8 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ auto dl = [&](void *h, const CudaDevicePtr &s) { CudaCheck(cudaMemcpy(h, s.get(), size_t(ng) * sizeof(double), cudaMemcpyDeviceToHost), "dl accum"); }; dl(swI, d.a_swI); dl(sw, d.a_sw); dl(swIh0, d.a_swIh0); dl(swIh1, d.a_swIh1); - dl(swh0, d.a_swh0); dl(swh1, d.a_swh1); dl(d_out, d.a_d); + dl(swh0, d.a_swh0); dl(swh1, d.a_swh1); dl(swh_typ0, d.a_swht0); dl(swh_typ1, d.a_swht1); + dl(d_out, d.a_d); CudaCheck(cudaMemcpy(nh0, d.a_nh0.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl nh0"); CudaCheck(cudaMemcpy(nh1, d.a_nh1.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl nh1"); CudaCheck(cudaMemcpy(rejected, d.a_rejected.get(), size_t(ng) * sizeof(int32_t), cudaMemcpyDeviceToHost), "dl rej"); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index cd4ac9464..fb853b016 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -93,13 +93,15 @@ public: // kernel and the host merge loop cannot disagree about it. // Outputs length n_groups; rejected[g] counts outliers dropped (reject_median uploaded, NAN // where none). rejected_obs is the per-full flag (length n_fulls): the host needs it for the reductions - // it still does itself, above all the anomalous I(+)/I(-) split. + // it still does itself, above all the anomalous I(+)/I(-) split. frame_cc_factor (length n_frames) is + // each frame's (G_ref/G)^2 floored at 1; swh_typ0/1 are the half-set weights multiplied by it. // Requires MergeEmSamples first (em_mean resident). void MergeAccum(double error_model_a, double error_model_b, bool error_model_active, bool reject_outliers, double reject_nsigma, const float *reject_median, - const uint8_t *half, + const uint8_t *half, const double *frame_cc_factor, double *swI, double *sw, double *swIh0, double *swIh1, - double *swh0, double *swh1, int32_t *nh0, int32_t *nh1, double *d_out, + double *swh0, double *swh1, double *swh_typ0, double *swh_typ1, + int32_t *nh0, int32_t *nh1, double *d_out, int32_t *rejected, uint8_t *on_ice_out, uint8_t *rejected_obs); // Per-group R_meas accumulators (sum|I_corr-merged_I|, sum_I, n, and the count this looser walk diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 31dd7676a..4794b6a7e 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -306,14 +306,17 @@ namespace { // other (signal) or drawn independently (noise, CC1/2 ~ 0), so a whole CC1/2 curve can be built // band by band. void AddBand(std::vector &v, std::mt19937 &rng, - double s_from, double s_to, int n, bool correlated) { - std::normal_distribution g(0.0, 1.0); + double s_from, double s_to, int n, bool correlated, double sigma = 1.0) { + std::normal_distribution g(0.0, sigma); for (int j = 0; j < n; ++j) { MergedReflection m; m.d = static_cast(1.0 / std::sqrt(s_from + (j + 0.5) * (s_to - s_from) / n)); const double a = g(rng), b = g(rng); m.I_half[0] = static_cast(a); m.I_half[1] = static_cast(correlated ? a : b); + // The pair's weight in a CC1/2: the precision it would have had at the typical frame scale + // over the precision it has - 1/sigma^2 for a band scaled up by sigma from dead frames. + m.cc_weight = static_cast(1.0 / (sigma * sigma)); v.push_back(m); } } @@ -336,6 +339,54 @@ TEST_CASE("ResolutionCutoff_CleanFallOff") { CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15)); } +// A sweep with a long stretch where the crystal barely diffracted: the reflections measured only there +// are noise scaled up by 1/G, sigmas with them, and there are many of them at every resolution. Counted +// by the information they carry they must not hide where the well-measured reflections stop. +TEST_CASE("ResolutionCutoff_ScaledUpNoiseDoesNotHideTheFallOff") { + Logger logger("test"); + std::mt19937 rng(12345); + std::vector merged; + AddBand(merged, rng, 0.01, 0.25, 480, true); + AddBand(merged, rng, 0.25, 0.51, 520, false); + AddBand(merged, rng, 0.01, 0.51, 300, false, 50.0); + + const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger); + REQUIRE(rc.d_fit); + 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. + 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)); +} + +// The frame factor of the CC1/2 weight is 1 on a sweep whose frames all sit at one scale, however that +// scale is spread over the frames, and 1 for any frame brighter than typical. +TEST_CASE("CCHalfFrameFactors_UniformScaleIsOne") { + const auto f = CCHalfFrameFactors({1.7, 1.7, 1.7, 1.7}, {100, 3, 250, 0}); + for (double x : f) CHECK(x == 1.0); + const auto g = CCHalfFrameFactors({1.0, 2.0}, {100, 100}); + CHECK(g[1] == 1.0); + CHECK(g[0] == Catch::Approx(9.0 / 5.0 * 9.0 / 5.0)); // G_ref = (1 + 8) / (1 + 4) +} + +// Most of the sweep at 2% of the good frames' scale: the typical scale is still the good frames', not +// the run median, and an observation from a dead frame carries 2500x the variance of a good one. +TEST_CASE("CCHalfFrameFactors_DeadStretchDoesNotSetTheTypicalScale") { + std::vector scale; + std::vector n; + for (int f = 0; f < 100; ++f) { + scale.push_back(f < 60 ? 0.02 : 1.0); + n.push_back(50); + } + scale.push_back(NAN); // a frame without a scale counts as G = 1 + n.push_back(50); + const auto factor = CCHalfFrameFactors(scale, n); + CHECK(factor[0] == Catch::Approx(2500.0).epsilon(0.01)); + CHECK(factor[99] == Catch::Approx(1.0).epsilon(1e-3)); + CHECK(factor[100] == Catch::Approx(1.0).epsilon(1e-3)); +} + // 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