mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases 67 Wiki Activity

Remove dead per-crystal deltaCChalf rejection

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-25 15:05:06 +02:00
co-authored by Claude Opus 4.8
parent ecf79af018
commit f72b4484e2
3 changed files with 0 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/Merge.cpp
-101
View File
@@ -478,107 +478,6 @@ void MergeOnTheFly::RefineErrorModel(const std::vector<IntegrationOutcome> &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<char> MergeOnTheFly::DeltaCChalfReject(const std::vector<IntegrationOutcome> &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<uint64_t, Acc> acc;
std::vector<int> 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<double>(r.I) * r.image_scale_corr;
const double s = static_cast<double>(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<int64_t>(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<double> 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<uint64_t, std::pair<double, double>> mine; // key -> (sum wI, sum w), count via .first
std::unordered_map<uint64_t, size_t> 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<size_t>(1, delta.size());
for (double d : delta) dv += (d - dm) * (d - dm);
const double dstd = std::sqrt(dv / std::max<size_t>(1, delta.size()));
const double cut = dm - nsigma * dstd;
std::vector<char> 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();
image_analysis/scale_merge/Merge.h
-6
View File
@@ -178,12 +178,6 @@ public:
// scaling and before RefineErrorModel/AddImage. Returns the applied held-out gain fraction (0 = no-op).
double RefineModulation(std::vector<IntegrationOutcome> &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<char> DeltaCChalfReject(const std::vector<IntegrationOutcome> &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.
rugnux/rugnux_cli.cpp
-7
View File
@@ -113,7 +113,6 @@ void print_usage() {
std::cout << " --min-captured-fraction <num> 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 <num> Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
std::cout << " --reject-outliers <num> 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 <num> 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 <num> Per-image CC limit in percent (default: no limit)" << std::endl;
std::cout << " --scaling-iterations <num> Number of scaling iterations with no reference data (default: 3)" << std::endl;
std::cout << " -z, --reference-mtz <file> 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<IntegratorMode> integrator_mode; // --integrator boxsum|gaussian|empirical
bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement
std::optional<double> outlier_reject_nsigma; // merge per-observation outlier rejection
std::optional<double> 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;