From f72b4484e2a253e02de1be6dde9a854e4a0ce272 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 25 Jul 2026 15:05:06 +0200 Subject: [PATCH] Remove dead per-crystal deltaCChalf rejection Drops the --reject-delta-cchalf flag and MergeOnTheFly::DeltaCChalfReject. The CLI value was parsed but never consumed (the method had no call site), so the flag was already a no-op. Wiring it up and testing against an external reference structure showed it is confirmation bias: on a spurious-crystal flood it raised internal CC1/2 while CCref (correlation to the true structure) fell, and it never improved R_meas. The merge weights are already correct; per-crystal merge-side rejection has no genuine lever here. Co-Authored-By: Claude Opus 4.8 --- image_analysis/scale_merge/Merge.cpp | 101 --------------------------- image_analysis/scale_merge/Merge.h | 6 -- rugnux/rugnux_cli.cpp | 7 -- 3 files changed, 114 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 0b7ec2cb..82f91431 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -478,107 +478,6 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_chi2 = chi2.empty() ? 0.0 : median(chi2) / CHI2_1_MEDIAN; } -// Per-crystal CC1/2-delta rejection (CrystFEL deltaCChalf style). Each image is assigned -// to one CC1/2 half, so removing an image only perturbs that half's per-reflection means. -// deltaCChalf_i = CC1/2(all) - CC1/2(without image i): a NEGATIVE value means removing the -// image RAISES CC1/2, i.e. it is inconsistent with the consensus. We flag images whose -// deltaCChalf is a low-side statistical outlier (< mean - nsigma*stddev). Reference-free. -// Two passes over the (retained) outcomes; per-image contributions are re-derived, not -// stored, so memory stays O(unique reflections + images) for full 200k-frame datasets. -std::vector MergeOnTheFly::DeltaCChalfReject(const std::vector &outcomes, - double nsigma) const { - struct Acc { double swI[2] = {0, 0}; double sw[2] = {0, 0}; size_t n[2] = {0, 0}; }; - std::unordered_map acc; - std::vector img_half(outcomes.size(), 0); - - // ---- pass 1: accumulate half-set sums, record each image's half ---- - auto contribution = [&](const Reflection &r, uint64_t &key, double &wI, double &w) -> bool { - if (generator.IsSystematicallyAbsent(r)) return false; - if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) return false; - if (!AcceptReflection(r, high_resolution_limit)) return false; - if (r.partiality < min_partiality) return false; - const double I = static_cast(r.I) * r.image_scale_corr; - const double s = static_cast(r.sigma) * r.image_scale_corr; - if (!std::isfinite(I) || !std::isfinite(s) || s <= 0.0) return false; - w = 1.0 / (s * s); - wI = w * I; - key = generator(r).pack(); - return true; - }; - - for (size_t i = 0; i < outcomes.size(); ++i) { - // Same deterministic half-set as the merge (HalfForImage), so deltaCChalf is measured on the - // exact CC1/2 split the statistics report - not an independent, order-dependent RNG draw. - const int half = HalfForImage(static_cast(i)); - img_half[i] = half; - for (const auto &r : outcomes[i].reflections) { - uint64_t key; double wI, w; - if (!contribution(r, key, wI, w)) continue; - auto &a = acc[key]; - a.swI[half] += wI; a.sw[half] += w; a.n[half]++; - } - } - - // ---- baseline Pearson over reflections present in BOTH halves (x=half0, y=half1) ---- - auto pearson = [](double N, double Sx, double Sy, double Sxx, double Syy, double Sxy) -> double { - const double cov = N * Sxy - Sx * Sy; - const double vx = N * Sxx - Sx * Sx, vy = N * Syy - Sy * Sy; - const double den = std::sqrt(vx * vy); - return den > 0.0 ? cov / den : 0.0; - }; - double N = 0, Sx = 0, Sy = 0, Sxx = 0, Syy = 0, Sxy = 0; - for (const auto &kv : acc) { - const auto &a = kv.second; - if (a.n[0] == 0 || a.n[1] == 0) continue; - const double x = a.swI[0] / a.sw[0], y = a.swI[1] / a.sw[1]; - N += 1; Sx += x; Sy += y; Sxx += x * x; Syy += y * y; Sxy += x * y; - } - const double cc_base = pearson(N, Sx, Sy, Sxx, Syy, Sxy); - - // ---- pass 2: leave-one-out deltaCChalf per image ---- - std::vector delta(outcomes.size(), 0.0); - for (size_t i = 0; i < outcomes.size(); ++i) { - const int h = img_half[i]; - // aggregate this image's contributions per reflection key (an image may, rarely, - // touch the same ASU reflection twice) - std::unordered_map> mine; // key -> (sum wI, sum w), count via .first - std::unordered_map mine_n; - for (const auto &r : outcomes[i].reflections) { - uint64_t key; double wI, w; - if (!contribution(r, key, wI, w)) continue; - auto &p = mine[key]; p.first += wI; p.second += w; mine_n[key]++; - } - double n = N, sx = Sx, sy = Sy, sxx = Sxx, syy = Syy, sxy = Sxy; - for (const auto &m : mine) { - const auto &a = acc.at(m.first); - if (a.n[0] == 0 || a.n[1] == 0) continue; // reflection not in CC1/2 - const double x0 = a.swI[0] / a.sw[0], y0 = a.swI[1] / a.sw[1]; - n -= 1; sx -= x0; sy -= y0; sxx -= x0 * x0; syy -= y0 * y0; sxy -= x0 * y0; - const double swI_h = a.swI[h] - m.second.first; - const double sw_h = a.sw[h] - m.second.second; - if (a.n[h] - mine_n[m.first] == 0 || sw_h <= 0.0) continue; // reflection drops half h - const double mean_h = swI_h / sw_h; - const double xnew = (h == 0) ? mean_h : x0; - const double ynew = (h == 1) ? mean_h : y0; - n += 1; sx += xnew; sy += ynew; sxx += xnew * xnew; syy += ynew * ynew; sxy += xnew * ynew; - } - delta[i] = cc_base - pearson(n, sx, sy, sxx, syy, sxy); - } - - // ---- reject low-side outliers: delta < mean - nsigma*stddev ---- - double dm = 0, dv = 0; - for (double d : delta) dm += d; - dm /= std::max(1, delta.size()); - for (double d : delta) dv += (d - dm) * (d - dm); - const double dstd = std::sqrt(dv / std::max(1, delta.size())); - const double cut = dm - nsigma * dstd; - - std::vector reject(outcomes.size(), 0); - for (size_t i = 0; i < outcomes.size(); ++i) - reject[i] = (outcomes[i].reflections.empty() ? 0 : (delta[i] < cut ? 1 : 0)); - return reject; -} - bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome, bool cc_mask) { if (reference_cell) { auto cell = outcome.latt.GetUnitCell(); diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index fdcd42c4..5b9dc093 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -178,12 +178,6 @@ public: // scaling and before RefineErrorModel/AddImage. Returns the applied held-out gain fraction (0 = no-op). double RefineModulation(std::vector &outcomes); - // Per-crystal CC1/2-delta rejection (CrystFEL deltaCChalf): returns a per-image flag - // marking images whose removal would raise CC1/2 by a low-side outlier amount - // (deltaCChalf < mean - nsigma*stddev). Skip the flagged images when merging. - [[nodiscard]] std::vector DeltaCChalfReject(const std::vector &outcomes, - double nsigma) const; - // d_min_override, when set, is the effective high-resolution limit for the shell table (used for // the automatic resolution cutoff computed by the caller); otherwise the manual // ScalingSettings high-resolution limit stands. The number of shells is ScalingSettings::ReportShellCount. diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 182d41b2..a490bf93 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -113,7 +113,6 @@ void print_usage() { std::cout << " --min-captured-fraction rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl; std::cout << " --mosaicity Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl; std::cout << " --reject-outliers Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl; - std::cout << " --reject-delta-cchalf Per-crystal CC1/2-delta rejection: drop images with deltaCChalf below mean - N*stddev (default: off; e.g. 2.5)" << std::endl; std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; @@ -173,7 +172,6 @@ enum { OPT_INTEGRATION_RADIUS, OPT_BACKGROUND_TRIM, OPT_REJECT_OUTLIERS, - OPT_REJECT_DELTA_CCHALF, OPT_REFERENCE_COLUMN, OPT_MODEL, OPT_DUMP_OBSERVATIONS, @@ -290,7 +288,6 @@ static option long_options[] = { {"simple-stills", no_argument, nullptr, OPT_SIMPLE_STILLS}, {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, - {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, {nullptr, 0, nullptr, 0} }; @@ -570,7 +567,6 @@ int main(int argc, char **argv) { std::optional integrator_mode; // --integrator boxsum|gaussian|empirical bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement std::optional outlier_reject_nsigma; // merge per-observation outlier rejection - std::optional delta_cchalf_nsigma; // per-crystal CC1/2-delta rejection if (argc == 1) { print_usage(); @@ -879,9 +875,6 @@ int main(int argc, char **argv) { case OPT_REJECT_OUTLIERS: outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger); break; - case OPT_REJECT_DELTA_CCHALF: - delta_cchalf_nsigma = parse_double_arg(optarg, "--reject-delta-cchalf", logger); - break; case OPT_MIN_IMAGE_CC: min_image_cc = parse_double_arg(optarg, "--min-image-cc", logger); break;