diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 9e5d055d..f77c3a7a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1675,67 +1675,16 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool std::vector 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 samples; - std::vector 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(fulls.size()); - std::vector gs2(nf), gI2(nf), gdev2(nf); - std::vector 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 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(o.sigma) * o.corr; - const double w = 1.0 / (sigma_corr * sigma_corr); - sw[o.group] += w; swI[o.group] += w * (static_cast(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(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(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> 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 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 &pool) { + // Takes the samples by value (it sorts them in place). auto fit_ab = [&](std::vector smp) { if (smp.size() < static_cast(8 * n_bins)) return; @@ -1778,22 +1727,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 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(8 * n_bins) && kept.size() < samples.size()) + if (kept.size() >= static_cast(8 * n_bins) && kept.size() < pool.size()) fit_ab(std::move(kept)); } + }; + { + std::vector 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(fulls.size()); + std::vector gs2(nf), gI2(nf), gdev2(nf); + std::vector 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 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(o.sigma) * o.corr; + const double w = 1.0 / (sigma_corr * sigma_corr); + sw[o.group] += w; swI[o.group] += w * (static_cast(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(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(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> 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 @@ -1803,8 +1810,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // fractional scatter of its symmetry mates estimates the systematic term, and the robust median over // 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 = error_model_b; - if (error_model_active) { + 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 gs(n_groups); for (const auto &o : fulls) { @@ -1863,17 +1872,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // 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); - } - // 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 - 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); + }; + estimate_asymptote(); auto corrected_sigma = [&](float I_corr, float sigma_corr, int g) -> float { if (!error_model_active) return sigma_corr; @@ -1882,94 +1882,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(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 acc(n_groups); - size_t reject_count = 0; - std::vector 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 aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups), - aswh0(n_groups), aswh1(n_groups), ad(n_groups); - std::vector 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(anh0[g]); a.nh[1] = static_cast(anh1[g]); - a.d = static_cast(ad[g]); - reject_count += static_cast(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(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 acc; Result result; - result.isa = isa_reported; - std::vector merged_I(n_groups, NAN); - float d_min = std::numeric_limits::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(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(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(a.swIh[i] / a.swh[i]); - mr.sigma_half[i] = static_cast(1.0 / std::sqrt(a.swh[i])); + std::vector merged_I; + size_t reject_count = 0; + std::vector 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 aswI(n_groups), asw(n_groups), aswIh0(n_groups), aswIh1(n_groups), + aswh0(n_groups), aswh1(n_groups), ad(n_groups); + std::vector 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(anh0[g]); a.nh[1] = static_cast(anh1[g]); + a.d = static_cast(ad[g]); + reject_count += static_cast(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(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(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(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(a.swIh[i] / a.swh[i]); + mr.sigma_half[i] = static_cast(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 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 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.