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

Merge branch 'fix23-work' into integration-variance-fixes

Browse Source
This commit is contained in:
leonarski_f
2026-08-09 20:58:41 +02:00
parent f4e281b2f5 97dbbc50b4
commit 1283e04ada
8 changed files with 268 additions and 235 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+233 -181
View File
@@ -1734,67 +1734,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
std::vector<float> reject_median(n_groups, NAN);
double error_model_a = 1.0, error_model_b = 0.0, error_model_chi2 = 0.0;
bool error_model_active = false;
{
struct Sample { double s2, I2, dev2; };
std::vector<Sample> samples;
std::vector<int32_t> cnt(n_groups, 0); // per-group usable count (both paths; feeds reject-median)
bool did_gpu = false;
#ifdef JFJOCH_USE_CUDA
if (use_gpu_merge) {
const int nf = static_cast<int>(fulls.size());
std::vector<double> gs2(nf), gI2(nf), gdev2(nf);
std::vector<uint8_t> gvalid(nf);
gpu_->MergeEmSamples(for_search, min_partiality, em_mean.data(), cnt.data(),
gs2.data(), gI2.data(), gdev2.data(), gvalid.data());
samples.reserve(nf);
for (int i = 0; i < nf; ++i)
if (gvalid[i]) samples.push_back({gs2[i], gI2[i], gdev2[i]});
did_gpu = true;
}
#endif
if (!did_gpu) {
// Per-group inverse-variance mean over usable fulls (>=2 obs), and the leverage-corrected samples.
std::vector<double> sw(n_groups, 0.0), swI(n_groups, 0.0);
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
const double sigma_corr = static_cast<double>(o.sigma) * o.corr;
const double w = 1.0 / (sigma_corr * sigma_corr);
sw[o.group] += w; swI[o.group] += w * (static_cast<double>(o.I) * o.corr); cnt[o.group]++;
}
for (int g = 0; g < n_groups; ++g)
if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g];
samples.reserve(fulls.size());
for (const auto &o : fulls) {
if (!usable_merge(o) || cnt[o.group] < 2) continue;
const double mean = em_mean[o.group];
if (!std::isfinite(mean)) continue;
const double sigma_corr = static_cast<double>(o.sigma) * o.corr;
const double s2 = sigma_corr * sigma_corr;
const double w = 1.0 / s2;
const double factor = 1.0 - w / sw[o.group];
if (factor < 0.05) continue;
const double resid = static_cast<double>(o.I) * o.corr - mean;
samples.push_back({s2, mean * mean, resid * resid / factor});
}
}
// Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on
// the GPU; cheap here, cnt >= 2 filter from the em pass). Fed to the merge accumulate.
if (reject_outliers) {
std::vector<std::vector<float>> iv(n_groups);
for (const auto &o : fulls)
if (usable_merge(o) && cnt[o.group] >= 2)
iv[o.group].push_back(o.I * o.corr);
for (int g = 0; g < n_groups; ++g)
if (!iv[g].empty()) {
std::nth_element(iv[g].begin(), iv[g].begin() + iv[g].size() / 2, iv[g].end());
reject_median[g] = iv[g][iv[g].size() / 2];
}
}
constexpr int n_bins = 16;
// 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).
// One leverage-corrected sample per usable full: its raw variance, its group's mean intensity, its
// squared deviation from that mean - and the resolution it sits at, because the fit is re-run below
// over the samples that survive the automatic resolution cutoff.
struct Sample { double s2, I2, dev2; float d; };
std::vector<Sample> samples;
constexpr int n_bins = 16;
// Fit (a, b) from the intensity-binned median deviations of a pool of samples, then refit on that
// pool's misfit-free subset. A lambda because the pool changes once the cutoff below is known.
auto fit_error_model = [&](const std::vector<Sample> &pool) {
// 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;
@@ -1837,22 +1786,80 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN;
}
};
fit_ab(samples);
fit_ab(pool);
// 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.
// 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) {
kept.reserve(pool.size());
for (const auto &s : pool) {
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())
if (kept.size() >= static_cast<size_t>(8 * n_bins) && kept.size() < pool.size())
fit_ab(std::move(kept));
}
};
{
std::vector<int32_t> cnt(n_groups, 0); // per-group usable count (both paths; feeds reject-median)
bool did_gpu = false;
#ifdef JFJOCH_USE_CUDA
if (use_gpu_merge) {
const int nf = static_cast<int>(fulls.size());
std::vector<double> gs2(nf), gI2(nf), gdev2(nf);
std::vector<uint8_t> gvalid(nf);
gpu_->MergeEmSamples(for_search, min_partiality, em_mean.data(), cnt.data(),
gs2.data(), gI2.data(), gdev2.data(), gvalid.data());
samples.reserve(nf);
for (int i = 0; i < nf; ++i)
if (gvalid[i]) samples.push_back({gs2[i], gI2[i], gdev2[i], fulls[i].d});
did_gpu = true;
}
#endif
if (!did_gpu) {
// Per-group inverse-variance mean over usable fulls (>=2 obs), and the leverage-corrected samples.
std::vector<double> sw(n_groups, 0.0), swI(n_groups, 0.0);
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
const double sigma_corr = static_cast<double>(o.sigma) * o.corr;
const double w = 1.0 / (sigma_corr * sigma_corr);
sw[o.group] += w; swI[o.group] += w * (static_cast<double>(o.I) * o.corr); cnt[o.group]++;
}
for (int g = 0; g < n_groups; ++g)
if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g];
samples.reserve(fulls.size());
for (const auto &o : fulls) {
if (!usable_merge(o) || cnt[o.group] < 2) continue;
const double mean = em_mean[o.group];
if (!std::isfinite(mean)) continue;
const double sigma_corr = static_cast<double>(o.sigma) * o.corr;
const double s2 = sigma_corr * sigma_corr;
const double w = 1.0 / s2;
const double factor = 1.0 - w / sw[o.group];
if (factor < 0.05) continue;
const double resid = static_cast<double>(o.I) * o.corr - mean;
samples.push_back({s2, mean * mean, resid * resid / factor, o.d});
}
}
// Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on
// the GPU; cheap here, cnt >= 2 filter from the em pass). Fed to the merge accumulate.
if (reject_outliers) {
std::vector<std::vector<float>> iv(n_groups);
for (const auto &o : fulls)
if (usable_merge(o) && cnt[o.group] >= 2)
iv[o.group].push_back(o.I * o.corr);
for (int g = 0; g < n_groups; ++g)
if (!iv[g].empty()) {
std::nth_element(iv[g].begin(), iv[g].begin() + iv[g].size() / 2, iv[g].end());
reject_median[g] = iv[g][iv[g].size() / 2];
}
}
fit_error_model(samples);
}
// Asymptotic I/sigma. ISa is by definition the I -> infinity limit of the signal-to-noise, i.e. the
// reproducibility of the strongest reflections (Diederichs, Acta Cryst. D66 (2010), 733-740). The
@@ -1860,12 +1867,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
// intensity raises its systematic term `b`, so 1/b understates that limit. Read it instead directly
// from the strong equivalents: for each well-measured reflection group the counting-subtracted
// fractional scatter of its symmetry mates estimates the systematic term, and the robust median over
// strong groups is the asymptote. The I/sigma threshold that selects "strong" is relaxed on weak /
// radiation-damaged data that has few strong reflections (rather than fall back to the higher
// whole-range b); it falls back only when even the relaxed set is too small. Host-side over the
// merged fulls, so CPU and GPU agree.
double error_model_b_asymptotic = error_model_b;
if (error_model_active) {
// strong groups is the asymptote. Report-only: it is the reported ISa, nothing downstream uses it.
// Host-side over the merged fulls, so CPU and GPU agree.
double error_model_b_asymptotic = 0.0;
auto estimate_asymptote = [&]() {
error_model_b_asymptotic = error_model_b;
if (!error_model_active) return;
struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; };
std::vector<GroupScatter> gs(n_groups);
for (const auto &o : fulls) {
@@ -1875,53 +1882,57 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
auto &g = gs[o.group];
g.sum += I_corr; g.sum_sq += I_corr * I_corr; g.sum_var += sigma_corr * sigma_corr; ++g.n;
}
// The counting term to subtract. The estimator below has always used the RAW reported
// sigma^2, but the (a, b) fit immediately above concluded the counting variance is
// a*sigma^2. Where a < 1 the raw subtraction overshoots: the per-group systematic is driven
// to zero, the median lands on the boundary, and 1/b_asy reports an impossible I/sigma
// (measured: one dataset reporting 64.6 where every other statistic supports ~16, and
// flipping between 10.9 and 64.6 on consecutive merges of the same data). Opt-in for now
// (env JFJOCH_ISA_ASY_A) so the default path stays bit-identical while the two are
// batteried against each other -- an earlier attempt at this subtraction was rejected for
// over-claiming on a > 1 data, which this change would reintroduce.
const double asy_counting_scale =
std::getenv("JFJOCH_ISA_ASY_A") != nullptr ? error_model_a : 1.0;
// Per-group counting-subtracted fractional systematic variance, paired with the group's I/sigma.
// Two conventions decide whether this measures anything:
// * the counting term is the error model's OWN, a*sigma^2. Subtracting the raw sigma^2 while the
// fit has just concluded the counting variance is a*sigma^2 leaves a bias of (1-a)*sigma^2,
// which at the I/sigma admitted here is the same size as the systematic being measured - so at
// a < 1 the per-group value collapses onto zero and 1/b_asy reports an impossible I/sigma.
// The I/sigma that selects "strong" uses the same counting term, so the gate and the
// subtraction agree about what a strong reflection is.
// * the sample variance of n observations is chi^2_(n-1)-distributed, and its MEDIAN lies below
// its mean (16% at n = 5, 5% at n = 13). Taking a median across groups of a variance minus an
// unbiased counting term therefore subtracts more than it should, again by an amount
// comparable to the systematic. Rescale each group's variance to be median-unbiased first
// (Wilson-Hilferty median(chi^2_k) = k*(1-2/9k)^3, exact to 0.4% for k >= 4) - the same
// median-to-mean conversion the (a, b) fit does with CHI2_1_MEDIAN.
// The clamp at zero is gone with them: median(max(x,0)) = max(median(x),0), so it never moved a
// positive median, and a negative median is the informative answer "not measurable here".
std::vector<std::pair<double, double>> group_scatter; // (systematic b^2, I/sigma)
for (const auto &g : gs) {
if (g.n < 5) continue;
const double mean = g.sum / g.n;
const double counting = asy_counting_scale * (g.sum_var / g.n);
const double counting = error_model_a * (g.sum_var / g.n);
if (mean <= 0.0 || counting <= 0.0) continue;
const double variance = (g.sum_sq - g.sum * g.sum / g.n) / (g.n - 1);
group_scatter.push_back({std::max(variance - counting, 0.0) / (mean * mean),
const double k = g.n - 1;
const double median_of_chi2 = std::pow(1.0 - 2.0 / (9.0 * k), 3);
const double variance = (g.sum_sq - g.sum * g.sum / g.n) / k / median_of_chi2;
group_scatter.push_back({(variance - counting) / (mean * mean),
mean / std::sqrt(counting)});
}
auto asymptote_above = [&](double snr_min, size_t min_groups) -> double {
// The threshold is relaxed on weak / radiation-damaged data that has too few strong reflections
// for the tight one. A tier that HAS enough groups gives its answer and is not retried lower
// because that answer came out small - the retry is what made the report flip between 1/b and an
// absurd value on consecutive merges of statistically identical data.
auto asymptote_above = [&](double snr_min, size_t min_groups) -> std::optional<double> {
std::vector<double> b2;
for (const auto &[b2_value, snr] : group_scatter)
if (snr >= snr_min) b2.push_back(b2_value);
return b2.size() >= min_groups ? std::sqrt(median_of(b2)) : 0.0;
return b2.size() >= min_groups ? std::optional(median_of(b2)) : std::nullopt;
};
double b_asy = asymptote_above(20.0, 100); // tight threshold on data that supports it
if (b_asy <= 0.0) b_asy = asymptote_above(10.0, 50); // relaxed for weak / damaged data
if (b_asy > 0.0) error_model_b_asymptotic = b_asy;
}
// Guard a degenerate low-multiplicity fit: with too few symmetry equivalents both the (a, b) fit and
// the per-group scatter collapse toward zero, and 1/b then reports an impossibly high asymptotic
// I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap report the asymptote
// as unmeasured rather than emit a spurious extreme. Only the REPORT is dropped: the sigma floor
// below keeps using the fitted value, because a degenerate fit is precisely where the raw
// 1/sqrt(sum_w) is least trustworthy and leaving the merged sigma uncapped would be the opposite of
// what this guard is for.
constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100
const double isa_reported = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B
? 1.0 / error_model_b_asymptotic : 0.0;
if (error_model_active)
logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b,
isa_reported, error_model_chi2);
auto b2_asy = asymptote_above(20.0, 100); // tight threshold on data that supports it
if (!b2_asy) b2_asy = asymptote_above(10.0, 50); // relaxed for weak / damaged data
// A non-positive median means the strong equivalents reproduce each other to within counting
// statistics: the asymptote is below what this data can resolve, so report the whole-range b
// rather than an extreme extrapolated from noise. The asymptote can also only ever REFINE 1/b
// upwards - that is the whole reason it is measured - so a strong-group estimate that comes out
// WORSE than the fit's own b has not measured an asymptote at all: it means "strong" was
// selected on a sigma scale that the fit itself rejects, which is what happens on data too weak
// to have strong reflections. Inert on healthy data (b_asy sits 0-40% below b there).
if (b2_asy && *b2_asy > 0.0)
error_model_b_asymptotic = std::min(std::sqrt(*b2_asy), error_model_b);
};
estimate_asymptote();
auto corrected_sigma = [&](float I_corr, float sigma_corr, int g) -> float {
if (!error_model_active) return sigma_corr;
@@ -1930,91 +1941,132 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
+ (error_model_b * I_for_b) * (error_model_b * I_for_b);
return v > 0.0 ? static_cast<float>(std::sqrt(v)) : sigma_corr;
};
// ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets. ----
// ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets, then
// export. In a lambda because the automatic resolution cutoff below re-runs it once the error
// model has been refitted on the reflections that survive the cut. ----
struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}; size_t nh[2] = {0, 0}; float d = NAN; };
std::vector<Accum> acc(n_groups);
size_t reject_count = 0;
std::vector<uint8_t> rejected_obs(fulls.size(), 0); // per-full outlier-rejected flag (both paths)
bool did_gpu_acc = false;
#ifdef JFJOCH_USE_CUDA
if (use_gpu_merge) {
std::vector<double> aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups),
aswh0(n_groups), aswh1(n_groups), ad(n_groups);
std::vector<int32_t> anh0(n_groups), anh1(n_groups), arej(n_groups);
gpu_->MergeAccum(error_model_a, error_model_b, error_model_active,
reject_outliers, reject_nsigma, reject_median.data(),
aswI.data(), asw.data(), aswIh0.data(), aswIh1.data(),
aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data(),
rejected_obs.data());
for (int g = 0; g < n_groups; ++g) {
Accum &a = acc[g];
a.swI = aswI[g]; a.sw = asw[g]; a.swIh[0] = aswIh0[g]; a.swIh[1] = aswIh1[g];
a.swh[0] = aswh0[g]; a.swh[1] = aswh1[g];
a.nh[0] = static_cast<size_t>(anh0[g]); a.nh[1] = static_cast<size_t>(anh1[g]);
a.d = static_cast<float>(ad[g]);
reject_count += static_cast<size_t>(arej[g]);
}
did_gpu_acc = true;
}
#endif
if (!did_gpu_acc)
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
const float I_corr = o.I * o.corr;
float sigma_corr = o.sigma * o.corr;
sigma_corr = corrected_sigma(I_corr, sigma_corr, o.group);
if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group])
&& std::fabs(I_corr - reject_median[o.group]) > reject_nsigma * sigma_corr) {
++reject_count;
rejected_obs[&o - fulls.data()] = 1;
continue;
}
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
const double wI = w * I_corr;
const int half = HalfForImage(o.frame);
auto &a = acc[o.group];
a.swI += wI; a.sw += w;
a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++;
if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d;
}
// ---- Export merged reflections (+ resolution-shell R-free flags). ----
std::vector<Accum> acc;
Result result;
result.isa = isa_reported;
std::vector<double> merged_I(n_groups, NAN);
float d_min = std::numeric_limits<float>::max(), d_max = 0.0f;
for (int g = 0; g < n_groups; ++g) {
const auto &a = acc[g];
if (a.sw <= 0.0) continue;
MergedReflection mr{};
mr.h = group_h[g]; mr.k = group_k[g]; mr.l = group_l[g];
mr.I = static_cast<float>(a.swI / a.sw);
// The systematic floor caps a high-multiplicity merged sigma at the reproducibility that the
// strongest reflections actually reach, i.e. the asymptotic term, so merged I/sigma approaches ISa.
mr.sigma = SigmaWithSystematicFloor(1.0 / std::sqrt(a.sw), mr.I, error_model_b_asymptotic);
mr.I_half[0] = mr.I_half[1] = NAN;
mr.sigma_half[0] = mr.sigma_half[1] = NAN;
mr.d = a.d;
if (a.nh[0] + a.nh[1] > 0 && a.swh[0] > 0.0 && a.swh[1] > 0.0) {
for (int i = 0; i < 2; ++i) {
mr.I_half[i] = static_cast<float>(a.swIh[i] / a.swh[i]);
mr.sigma_half[i] = SigmaWithSystematicFloor(1.0 / std::sqrt(a.swh[i]), mr.I_half[i], error_model_b_asymptotic);
std::vector<double> merged_I;
size_t reject_count = 0;
std::vector<uint8_t> rejected_obs; // per-full outlier-rejected flag (both paths)
auto run_merge = [&]() {
acc.assign(n_groups, Accum{});
result.merged.clear();
merged_I.assign(n_groups, NAN);
reject_count = 0;
rejected_obs.assign(fulls.size(), 0);
bool did_gpu_acc = false;
#ifdef JFJOCH_USE_CUDA
if (use_gpu_merge) {
std::vector<double> aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups),
aswh0(n_groups), aswh1(n_groups), ad(n_groups);
std::vector<int32_t> anh0(n_groups), anh1(n_groups), arej(n_groups);
gpu_->MergeAccum(error_model_a, error_model_b, error_model_active,
reject_outliers, reject_nsigma, reject_median.data(),
aswI.data(), asw.data(), aswIh0.data(), aswIh1.data(),
aswh0.data(), aswh1.data(), anh0.data(), anh1.data(), ad.data(), arej.data(),
rejected_obs.data());
for (int g = 0; g < n_groups; ++g) {
Accum &a = acc[g];
a.swI = aswI[g]; a.sw = asw[g]; a.swIh[0] = aswIh0[g]; a.swIh[1] = aswIh1[g];
a.swh[0] = aswh0[g]; a.swh[1] = aswh1[g];
a.nh[0] = static_cast<size_t>(anh0[g]); a.nh[1] = static_cast<size_t>(anh1[g]);
a.d = static_cast<float>(ad[g]);
reject_count += static_cast<size_t>(arej[g]);
}
did_gpu_acc = true;
}
if (!std::isfinite(a.d) || a.d <= 0.0f) continue;
d_min = std::min(d_min, a.d);
d_max = std::max(d_max, a.d);
merged_I[g] = mr.I;
result.merged.push_back(mr);
}
#endif
if (!did_gpu_acc)
for (const auto &o : fulls) {
if (!usable_merge(o)) continue;
const float I_corr = o.I * o.corr;
float sigma_corr = o.sigma * o.corr;
sigma_corr = corrected_sigma(I_corr, sigma_corr, o.group);
if (reject_outliers && error_model_active && std::isfinite(reject_median[o.group])
&& std::fabs(I_corr - reject_median[o.group]) > reject_nsigma * sigma_corr) {
++reject_count;
rejected_obs[&o - fulls.data()] = 1;
continue;
}
const double w = 1.0 / (static_cast<double>(sigma_corr) * sigma_corr);
const double wI = w * I_corr;
const int half = HalfForImage(o.frame);
auto &a = acc[o.group];
a.swI += wI; a.sw += w;
a.swIh[half] += wI; a.swh[half] += w; a.nh[half]++;
if (!std::isfinite(a.d) && std::isfinite(o.d) && o.d > 0.0f) a.d = o.d;
}
// ---- Export merged reflections. ----
for (int g = 0; g < n_groups; ++g) {
const auto &a = acc[g];
if (a.sw <= 0.0) continue;
MergedReflection mr{};
mr.h = group_h[g]; mr.k = group_k[g]; mr.l = group_l[g];
mr.I = static_cast<float>(a.swI / a.sw);
// Plain inverse-variance merged sigma. The error model's systematic term (b*I)^2 is measured
// from the scatter BETWEEN a reflection's symmetry equivalents, i.e. from exactly the part of
// the error that is NOT common to them, so it averages down over the multiplicity like the
// counting part and the merge must not hold it back. XDS behaves the same way: its merged
// I/sigma runs far above its own reported ISa in the strong low-resolution shells.
mr.sigma = static_cast<float>(1.0 / std::sqrt(a.sw));
mr.I_half[0] = mr.I_half[1] = NAN;
mr.sigma_half[0] = mr.sigma_half[1] = NAN;
mr.d = a.d;
if (a.nh[0] + a.nh[1] > 0 && a.swh[0] > 0.0 && a.swh[1] > 0.0) {
for (int i = 0; i < 2; ++i) {
mr.I_half[i] = static_cast<float>(a.swIh[i] / a.swh[i]);
mr.sigma_half[i] = static_cast<float>(1.0 / std::sqrt(a.swh[i]));
}
}
if (!std::isfinite(a.d) || a.d <= 0.0f) continue;
merged_I[g] = mr.I;
result.merged.push_back(mr);
}
};
run_merge();
// Automatic high-resolution cutoff (post-merge): trim the written reflections + reported shells to
// the CC1/2 fall-off. The scaling, combine and error model above already ran over the full range,
// and the per-image _process.h5 is written elsewhere from the partials, so no data is lost. A manual
// the CC1/2 fall-off. The scaling and combine above ran over the full range, and the per-image
// _process.h5 is written elsewhere from the partials, so no data is lost. A manual
// --scaling-high-resolution (d_min_limit) wins; the P1 search merge (for_search) is never cut, so
// the space-group search still sees the full range.
const std::optional<double> effective_d_min = ApplyResolutionCutoff(
result.merged, d_min_limit, resolution_cutoff_method, resolution_cc_target, for_search, logger);
// The error model has to be calibrated on the reflections that are kept, not on the ones that are
// thrown away: on a default run the cut can remove the majority of the measured range, and a fit
// spanning it is dominated by reflections that are not written (measured on one run: a = 0.28,
// b = 0.159 over the full range against a = 0.42, b = 0.119 over the kept one, on 40% of the
// observations the whole-range fit saw). A manual limit needs nothing here - it already
// restricted the observations at ingest - so this is the automatic cut catching up with it.
// The circularity is resolved by direction: the cutoff is read from the provisional merge, and
// CC1/2 is a correlation of the two half-set MEANS, which the sigma scale barely moves; the sigmas
// are then calibrated on the population the cutoff chose. One refinement, not an iteration.
if (effective_d_min && effective_d_min != d_min_limit) {
std::vector<Sample> in_range;
in_range.reserve(samples.size());
for (const auto &s : samples)
if (s.d >= *effective_d_min) in_range.push_back(s);
fit_error_model(in_range);
estimate_asymptote();
run_merge();
std::erase_if(result.merged, [&](const MergedReflection &m) {
return std::isfinite(m.d) && m.d < *effective_d_min;
});
}
// Guard a degenerate low-multiplicity fit: with too few symmetry equivalents both the (a, b) fit and
// the per-group scatter collapse toward zero, and 1/b then reports an impossibly high asymptotic
// I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap report the asymptote
// as unmeasured rather than emit a spurious extreme.
constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100
result.isa = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B ? 1.0 / error_model_b_asymptotic : 0.0;
if (error_model_active)
logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b,
result.isa, error_model_chi2);
AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction);
// French-Wilson (F, and F(+)/F(-) from the anomalous split) is deferred until after the anomalous
// accumulator below has attached I(+)/I(-), so the two hands get their amplitudes in one pass.