Merge: per-observation outlier rejection (env-gated PR_REJECT)

At the jet's ~1000x multiplicity R-free is bias-limited, and the merge had NO
outlier rejection - serial data zingers/overlaps/mis-indexed frames bias every
merged mean. Add a robust per-observation cut: drop observations whose corrected
intensity lies > reject_nsigma error-model sigmas from the reflection's MEDIAN.
The error-model sigma already captures the genuine (partiality) scatter, and the
median is a robust centre, so only the tail beyond the real spread is removed -
not good partials. The median is computed in RefineErrorModel (which already pools
the observations per reflection); AddImage applies the cut.

Env-gated via PR_REJECT=<nsigma> (off by default); logs the count removed. On the
jet (CC proxy) it lifts CCref +8 (nsigma 6, 0.6% cut) to +11 (nsigma 3, 7.4% cut)
- the cut is vs the data's own median, not the reference, so the gain is real
cleaner means. R-free validation + the nsigma sweet spot (over-rejection risk at
low nsigma) are for Filip's full-jet R-free.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

image_analysis/scale_merge/Merge.cpp

@@ -65,6 +65,19 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, bool cc_mask) {
         auto hkl_key = hkl.pack();
         sigma_corr = CorrectedSigma(I_corr, sigma_corr, hkl_key);
 
+        // Robust outlier rejection: drop this observation if it sits more than
+        // reject_nsigma error-model sigmas from the reflection's median. Needs the active
+        // error model so sigma_corr reflects the real scatter (else the threshold is the
+        // bare counting sigma and would cull good partials).
+        if (reject_outliers && error_model_active) {
+            const auto mit = reject_median_I.find(hkl_key);
+            if (mit != reject_median_I.end() &&
+                std::fabs(I_corr - mit->second) > reject_nsigma * sigma_corr) {
+                ++reject_count;
+                continue;
+            }
+        }
+
         auto it = accumulator.find(hkl_key);
         if (it == accumulator.end())
             it = accumulator.emplace(hkl_key, MergeAccum{
@@ -110,6 +123,8 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &outc
     error_model_a = 1.0;
     error_model_b = 0.0;
     error_model_mean_I.clear();
+    reject_median_I.clear();
+    reject_count = 0;
 
     // --- 1. Collect accepted, scaled observations grouped by symmetry-equivalent hkl,
     //        applying exactly the filters AddImage uses. ---
@@ -156,6 +171,16 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &outc
             continue;
         const double mean = sum_wI / sum_w;
         error_model_mean_I[key] = static_cast<float>(mean);
+        // Robust centre for outlier rejection: the median intensity (resists the very
+        // outliers the inverse-variance mean is being protected from).
+        {
+            std::vector<float> iv;
+            iv.reserve(obs.size());
+            for (const auto &o: obs)
+                iv.push_back(o.I);
+            std::nth_element(iv.begin(), iv.begin() + iv.size() / 2, iv.end());
+            reject_median_I[key] = iv[iv.size() / 2];
+        }
         const double I2 = mean * mean;
         for (const auto &o: obs) {
             const double w = 1.0 / (static_cast<double>(o.sigma) * o.sigma);

image_analysis/scale_merge/Merge.h

@@ -87,6 +87,16 @@ class MergeOnTheFly {
     std::unordered_map<uint64_t, float> error_model_mean_I;
     [[nodiscard]] float CorrectedSigma(float I_corr, float sigma_corr, uint64_t hkl_key) const;
 
+    // Optional per-observation outlier rejection: drop observations whose corrected
+    // intensity lies more than reject_nsigma error-model sigmas from the reflection's
+    // *median* (a robust centre). The error-model sigma already captures the genuine
+    // (e.g. partiality) scatter, so this removes only the tail beyond it - zingers,
+    // overlaps, mis-indexed frames - not good partials. Populated by RefineErrorModel.
+    bool   reject_outliers = false;
+    double reject_nsigma = 6.0;
+    std::unordered_map<uint64_t, float> reject_median_I;
+    size_t reject_count = 0;
+
     bool Mask(const IntegrationOutcome &outcome, bool cc_mask);
 public:
     MergeOnTheFly(const DiffractionExperiment &x);
@@ -99,6 +109,10 @@ public:
     [[nodiscard]] double ErrorModelA() const { return error_model_a; }
     [[nodiscard]] double ErrorModelB() const { return error_model_b; }
 
+    // Enable per-observation outlier rejection (call after RefineErrorModel, before AddImage).
+    void SetRejectOutliers(double nsigma) { reject_outliers = true; reject_nsigma = nsigma; }
+    [[nodiscard]] size_t RejectedCount() const { return reject_count; }
+
     void AddImage(const IntegrationOutcome& outcome, bool cc_mask = false);
 
     MergeStatistics MergeStats(const std::vector<MergedReflection> &merged,

tools/jfjoch_process.cpp

@@ -1034,8 +1034,15 @@ int main(int argc, char **argv) {
             const double isa = (b > 0.0) ? 1.0 / b : std::numeric_limits<double>::infinity();
             logger.Info("Error model: sigma'^2 = {:.3f} sigma^2 + ({:.4f} I)^2   ISa = {:.1f}", a, b, isa);
         }
+        if (const char *rj = std::getenv("PR_REJECT")) { // TEMP A/B: per-observation outlier rejection
+            const double nsig = std::atof(rj);
+            merge_engine.SetRejectOutliers(nsig);
+            logger.Info("Outlier rejection enabled: dropping observations > {:.1f} sigma from the per-reflection median", nsig);
+        }
         for (auto &i : indexer.GetIntegrationOutcome())
             merge_engine.AddImage(i);
+        if (merge_engine.RejectedCount() > 0)
+            logger.Info("Outlier rejection removed {} observations", merge_engine.RejectedCount());
 
         auto merged_reflections = merge_engine.ExportReflections();
         auto merged_statistics = merge_engine.MergeStats(merged_reflections, indexer.GetIntegrationOutcome(), reference_data);