From 4ec197c68c86f8d932cc476cd11e5ea2d1541dee Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 4 Oct 2026 20:51:01 +0200 Subject: [PATCH] Rotation merge: outlier band's multiplicative spread measured, not taken from b The multiplicative outlier band (c8de5f8d6) took its width from the error model's b. b is a Gaussian width in I and grows with whatever tail the fulls have: on a low-ISa hexagonal small-molecule set with a population of fulls lost to near zero, b = 0.54 against a core ln-spread of 0.25, and exp(6 b) = 25 let the high outliers of the weak reflections through (SHELXL R1 0.191 -> 0.333 there, 0.111 -> 0.124 on its 25 keV sweep). The spread is now measured on the merge's own fulls: the weighted median of |ln(I / median)|, weighted by (/sigma_counting)^2 so the fulls whose ratio counting noise does not blur carry it; b remains the fallback where nothing can be measured. Battery (29 sets: 13 small-molecule, 16 protein), SHELXL R1 against the median fix alone: the strongly absorbing cubic set 0.116 -> 0.0998 (b-band 0.103); the low-ISa hexagonal set 0.191 -> 0.174 at 20 keV, 0.111 -> 0.101 at 25 keV; organics within +-0.0003 or better (cytidine 0.0617 -> 0.0613, lalanine 0.0541 -> 0.0533). Proteins unchanged except two low-ISa sets whose CC1/2 cutoff moves: 6yqf 3.32 -> 2.89 A (placement R-free at a fixed 3.32 A 0.4329 -> 0.4326, at 2.89 A 0.4642 -> 0.4630), and the split myoglobin set 1.77 -> 1.94 A (low-resolution R_meas 33.0% -> 27.3%). Private subset: two sets' R_meas lower, the rest unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- image_analysis/scale_merge/OutlierBand.h | 7 ++-- .../scale_merge/RotationScaleMerge.cpp | 42 +++++++++++++++++-- .../scale_merge/RotationScaleMergeGPU.cu | 7 ++-- .../scale_merge/RotationScaleMergeGPU.h | 3 +- 4 files changed, 48 insertions(+), 11 deletions(-) diff --git a/image_analysis/scale_merge/OutlierBand.h b/image_analysis/scale_merge/OutlierBand.h index c8cfb0610..8eed99152 100644 --- a/image_analysis/scale_merge/OutlierBand.h +++ b/image_analysis/scale_merge/OutlierBand.h @@ -13,9 +13,10 @@ // linear band reaches far below the median and only a little above it. On a strongly absorbing crystal // (b = 0.34) every full the test removed read 2-7x its median - the least absorbed fulls, the ones // nearest the true intensity - and none read low, so the merged intensities were cut from above only. -// Here the multiplicative part of the band is exp(+-n*b) about the expected intensity, which is the -// linear band to first order in n*b, and the counting part stays linear. -// up = exp(n*b) - 1 and down = 1 - exp(-n*b) are per-merge constants, taken on the host and handed to +// Here the multiplicative part of the band is exp(+-n*s) about the expected intensity, s the spread of +// ln(I / median) measured on the merge's own fulls (MergeAndStats), and the counting part stays linear. +// With s = b it is the linear band to first order in n*b. +// up = exp(n*s) - 1 and down = 1 - exp(-n*s) are per-merge constants, taken on the host and handed to // the device, so that both sides compute them from the same library call. #include diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 12613a3f2..ff01d4037 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4525,6 +4525,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // ---- Error model: fit dev2 = a*sigma^2 + b^2*^2 from symmetry-equivalent scatter. ---- std::vector em_mean(n_groups, NAN); std::vector reject_median(n_groups, NAN); + double reject_log_sd = NAN; // the measured multiplicative spread the outlier band uses (OutlierBand.h) std::vector reject_var_add; // per group: the shell's measured share of dI^2/4 double error_model_a = 1.0, error_model_b = 0.0; bool error_model_active = false; @@ -4751,6 +4752,37 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (int g = 0; g < n_groups; ++g) reject_median[g] = cnt[g] >= 3 ? set_median[g] : !pair_needed.empty() ? set_median[n_groups + pair_of_group[g]] : NAN; + + // The spread of ln(I / median) the band's multiplicative part is drawn from: its weighted median + // absolute value, the weight (/sigma_counting)^2 so the fulls whose ratio counting noise does + // not blur carry it. Measured here rather than taken from the error model's b, which is a + // Gaussian width in I and grows with whatever tail the fulls have: a run with a population of + // fulls lost to near zero fits b = 0.54 against a core spread of 0.25, and exp(6 b) then let + // every high outlier of the weak reflections through. + std::vector> lv; + for (int i = 0; i < n_full; ++i) { + const int g = mf.group[i]; + if (g < 0 || !std::isfinite(reject_median[g]) || !(reject_median[g] > 0.0f) + || !std::isfinite(em_mean[g]) || !(em_mean[g] > 0.0)) + continue; + const double I_corr = static_cast(mf.I[i]) * mf.corr[i]; + const double sc = static_cast(mf.sigma[i]) * mf.corr[i]; + const double var = counting_variance(fulls[i], em_mean[g], sc * sc); + if (!(I_corr > 0.0) || !(var > 0.0)) continue; + lv.emplace_back(std::fabs(std::log(I_corr / reject_median[g])), em_mean[g] * em_mean[g] / var); + } + if (!lv.empty()) { + std::sort(lv.begin(), lv.end()); + double total = 0.0; + for (const auto &[v, w] : lv) total += w; + double run = 0.0; + size_t j = 0; + for (; j + 1 < lv.size(); ++j) { + run += lv[j].second; + if (run >= 0.5 * total) break; + } + reject_log_sd = 1.4826 * lv[j].first; + } } fit_error_model(samples); } @@ -4940,11 +4972,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (size_t j = 0; j < wobs.size(); ++j) if (wilson.rejected[j]) rejected_obs[wilson_full[j]] = 1; } + // The multiplicative part of the outlier band (OutlierBand.h): the measured spread, or the error + // model's b where nothing could be measured. + const double mult_sd = std::isfinite(reject_log_sd) ? reject_log_sd : error_model_b; + const double reject_up = std::expm1(reject_nsigma * mult_sd); + const double reject_down = -std::expm1(-reject_nsigma * mult_sd); bool did_gpu_acc = false; #ifdef JFJOCH_USE_CUDA if (use_gpu_merge) { gpu_->MergeAccum(error_model_a, error_model_b, error_model_active, - reject_outliers, reject_nsigma, reject_median.data(), + reject_outliers, reject_nsigma, mult_sd, reject_median.data(), reject_var_add.empty() ? nullptr : reject_var_add.data(), merge_half.data(), frame_cc_factor.data(), rejected_obs.data()); // Downloaded a slice of groups at a time and unpacked straight into acc, so the host never @@ -4982,9 +5019,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool did_gpu_acc = true; } #endif - // The multiplicative part of the outlier band (OutlierBand.h), from the error model as it stands. - const double reject_up = std::expm1(reject_nsigma * error_model_b); - const double reject_down = -std::expm1(-reject_nsigma * error_model_b); if (!did_gpu_acc) for (const auto &o : fulls) { if (!usable_merge(o)) continue; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 4681bbd0f..cf19a59ee 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -1197,7 +1197,8 @@ 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, + bool reject_outliers, double reject_nsigma, double reject_mult_sd, + const float *reject_median, const float *reject_var_add, const uint8_t *half, const double *frame_cc_factor, uint8_t *rejected_obs) { @@ -1225,8 +1226,8 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.error_model_a = error_model_a; p.error_model_b = error_model_b; p.error_model_active = error_model_active ? 1 : 0; p.reject_outliers = reject_outliers ? 1 : 0; p.reject_nsigma = reject_nsigma; - p.reject_up = std::expm1(reject_nsigma * error_model_b); - p.reject_down = -std::expm1(-reject_nsigma * error_model_b); + p.reject_up = std::expm1(reject_nsigma * reject_mult_sd); + p.reject_down = -std::expm1(-reject_nsigma * reject_mult_sd); p.reject_median = d.reject_median.get(); p.reject_var_add = d.reject_var_add.get(); p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 4a75a7874..0f319e651 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -105,8 +105,9 @@ public: // half-set weights multiplied by it. Requires MergeEmSamples first (em_mean resident). // reject_var_add (n_groups) widens the pooled cut by the shell's own measured Bijvoet // variance; null leaves the plain n-sigma test. + // reject_mult_sd is the multiplicative spread of the outlier band (OutlierBand.h). void MergeAccum(double error_model_a, double error_model_b, bool error_model_active, - bool reject_outliers, double reject_nsigma, const float *reject_median, + bool reject_outliers, double reject_nsigma, double reject_mult_sd, const float *reject_median, const float *reject_var_add, const uint8_t *half, const double *frame_cc_factor, uint8_t *rejected_obs); // Download groups [g0, g0 + n) of what MergeAccum left on the device; every output has length n.