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

rugnux: exclude merge-rejected outliers from the error-model fit

Browse Source 
The rotation error-model (a,b) fit built its sample pool from all usable
fulls, while the merge separately dropped symmetry outliers (> Nsigma from
the per-reflection median). So the fitted sigmas were calibrated against a
dirtier pool than the reflections that actually entered the merge.

Refactor the intensity-binned median fit into a `fit_ab` lambda and re-run
it once on the misfit-free pool (samples with dev2 <= reject_nsigma^2 *
model variance), matching the merge's own rejection. It operates on the
shared `samples` array, so the CPU and GPU error-model paths stay
bit-identical.

Neutral on lyso_ref (the median-binned fit is already robust to the few
outliers) - a correctness fix that matters more on dirtier data. Gated on
reject_outliers (default 6 sigma for rot3d).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-12 10:04:43 +02:00
co-authored by Claude Opus 4.8
parent a854523fee
commit dee679380b
2 changed files with 28 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CPU_DATA_ANALYSIS.md
+1 -1
View File
@@ -604,6 +604,6 @@ A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of
- **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted (within an excitation-error window) and scaled with unit partiality.
- **Space-group determination.** When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities.
- **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ are computed and written to the merged mmCIF as a twinning diagnostic.
- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default).
- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool.
- **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it.
- **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream.
image_analysis/scale_merge/RotationScaleMerge.cpp
+27 -8
View File
@@ -913,17 +913,20 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
}
}
constexpr int n_bins = 16;
if (samples.size() >= static_cast<size_t>(8 * n_bins)) {
std::sort(samples.begin(), samples.end(),
[](const Sample &a, const Sample &b) { return a.I2 < b.I2; });
// Fit (a, b) from the intensity-binned median deviations. Factored into a lambda so it can be
// re-run on a misfit-free pool below; takes the samples by value (it sorts them in place).
auto fit_ab = [&](std::vector<Sample> smp) {
if (smp.size() < static_cast<size_t>(8 * n_bins))
return;
std::sort(smp.begin(), smp.end(), [](const Sample &a, const Sample &b) { return a.I2 < b.I2; });
std::vector<double> bs2, bI2, bd2;
const size_t per = samples.size() / n_bins;
const size_t per = smp.size() / n_bins;
for (int bin = 0; bin < n_bins; ++bin) {
const size_t lo = bin * per;
const size_t hi = (bin == n_bins - 1) ? samples.size() : lo + per;
const size_t hi = (bin == n_bins - 1) ? smp.size() : lo + per;
std::vector<double> vs2, vI2, vd2;
for (size_t i = lo; i < hi; ++i) {
vs2.push_back(samples[i].s2); vI2.push_back(samples[i].I2); vd2.push_back(samples[i].dev2);
vs2.push_back(smp[i].s2); vI2.push_back(smp[i].I2); vd2.push_back(smp[i].dev2);
}
bs2.push_back(median_of(vs2));
bI2.push_back(median_of(vI2));
@@ -946,13 +949,29 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
error_model_b = std::sqrt(b2);
error_model_active = true;
std::vector<double> chi2;
chi2.reserve(samples.size());
for (const auto &s : samples) {
chi2.reserve(smp.size());
for (const auto &s : smp) {
const double v = error_model_a * s.s2 + b2 * s.I2;
if (v > 0.0) chi2.push_back(s.dev2 / v);
}
error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN;
}
};
fit_ab(samples);
// Refit on a misfit-free pool: the merge drops symmetry outliers (|I - median| > reject_nsigma *
// sigma) from the merged intensity, so drop the equivalent samples (dev2 > reject_nsigma^2 * model
// variance) from the error-model fit too, keeping the fitted sigmas consistent with the reflections
// that actually survive. Operates on the shared `samples`, so CPU and GPU stay bit-identical.
if (reject_outliers && error_model_active) {
const double ns2 = reject_nsigma * reject_nsigma, b2 = error_model_b * error_model_b;
std::vector<Sample> kept;
kept.reserve(samples.size());
for (const auto &s : samples) {
const double v = error_model_a * s.s2 + b2 * s.I2;
if (v > 0.0 && s.dev2 <= ns2 * v) kept.push_back(s);
}
if (kept.size() >= static_cast<size_t>(8 * n_bins) && kept.size() < samples.size())
fit_ab(std::move(kept));
}
}
if (error_model_active)