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 (<I>/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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB

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 <cmath>

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*<I>^2 from symmetry-equivalent scatter. ----
     std::vector<double> em_mean(n_groups, NAN);
     std::vector<float> reject_median(n_groups, NAN);
+    double reject_log_sd = NAN;   // the measured multiplicative spread the outlier band uses (OutlierBand.h)
     std::vector<float> 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 (<I>/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<std::pair<double, double>> 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<double>(mf.I[i]) * mf.corr[i];
+                const double sc = static_cast<double>(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;

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();

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.