diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 985ce4692..b54a9fa7d 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -28,7 +28,9 @@ ADD_LIBRARY(JFJochScaleMerge ScalingResult.h ScalingResult.cpp WilsonOutliers.cpp - WilsonOutliers.h) + WilsonOutliers.h + ErrorModel.cpp + ErrorModel.h) TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon fftw3f) IF (JFJOCH_CUDA_AVAILABLE) diff --git a/image_analysis/scale_merge/ErrorModel.cpp b/image_analysis/scale_merge/ErrorModel.cpp new file mode 100644 index 000000000..4e68cdac6 --- /dev/null +++ b/image_analysis/scale_merge/ErrorModel.cpp @@ -0,0 +1,140 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ErrorModel.h" + +#include +#include +#include +#include + +#include "../../common/ParallelFor.h" + +ErrorModelFit FitErrorModel(const std::vector &pool, std::vector &smp, + size_t nthreads) { + constexpr int n_bins = 16; + ErrorModelFit fit; + if (pool.size() < static_cast(8 * n_bins)) + return fit; + + // The rank key is computed once, here, and carried with the sample: a key recomputed inside the + // comparator can round differently at different inlined call sites, and nth_element then sees an + // inconsistent order. + smp.resize(pool.size()); + ParallelChunks(static_cast(pool.size()), ThreadsForWork(pool.size(), nthreads), [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) + smp[i] = {pool[i].I2 / pool[i].s2, pool[i].s2, pool[i].I2, pool[i].dev2, 0.0}; + }); + + // Equal-count bins, each boundary put in place with nth_element; the split runs down the middle so + // each level halves the range it works on. Total order, so which of two equal keys lands in which bin + // is a property of the samples, not of the order they arrived in. The two halves of a split are + // disjoint ranges, so they can run side by side while the range is still worth a thread. + const size_t per = smp.size() / n_bins; + const auto by_snr = [](const ErrorModelBinned &a, const ErrorModelBinned &b) { + if (a.snr2 != b.snr2) return a.snr2 < b.snr2; + if (a.s2 != b.s2) return a.s2 < b.s2; + if (a.I2 != b.I2) return a.I2 < b.I2; + return a.dev2 < b.dev2; + }; + constexpr size_t SPLIT_MIN_PARALLEL = 1u << 17; + const std::function split = + [&](size_t lo, size_t hi, int b0, int b1, int depth) { + if (b0 >= b1) return; + const int mid = (b0 + b1) / 2; + const size_t k = static_cast(mid) * per; + std::nth_element(smp.begin() + lo, smp.begin() + k, smp.begin() + hi, by_snr); + if (depth > 0 && hi - lo > SPLIT_MIN_PARALLEL) { + auto left = std::async(std::launch::async, + [&, lo, k, b0, mid, depth] { split(lo, k, b0, mid, depth - 1); }); + split(k, hi, mid + 1, b1, depth - 1); + left.get(); + } else { + split(lo, k, b0, mid, 0); + split(k, hi, mid + 1, b1, 0); + } + }; + split(0, smp.size(), 1, n_bins, nthreads > 1 ? 3 : 0); + auto bin_lo = [&](int bin) { return static_cast(bin) * per; }; + auto bin_hi = [&](int bin) { return bin == n_bins - 1 ? smp.size() : static_cast(bin + 1) * per; }; + + // b is identified only if the strongest bin reaches where b* can outgrow the counting term: when + // fewer reflections are strong than one bin holds, that bin sits at I/sigma ~ 1 and a two-term fit + // hands b the bins' own noise instead (measured: b = 5.6, ISa 0.5, on a crystal XDS gives 6.1). Then + // fit a alone and leave b at zero. Over the rotation battery the two crystals this fires on sit at + // (I/sigma)^2 0.22 and 0.84 while the next is 31.7, so the threshold is not delicate. + constexpr double B_LEVER_MIN = 4.0; + { + const auto top_mid = smp.begin() + (bin_lo(n_bins - 1) + bin_hi(n_bins - 1)) / 2; + std::nth_element(smp.begin() + bin_lo(n_bins - 1), top_mid, smp.end(), by_snr); + fit.b_measured = top_mid->snr2 > B_LEVER_MIN; + } + + // The median of a chi^2 with one degree of freedom: a normalised squared deviation that is right + // on average has half of its values below this. + constexpr double CHI2_1_MEDIAN = 0.454936; + std::vector bs2(n_bins), bI2(n_bins), bd2(n_bins); + double a = 1.0, b2 = 0.0; + for (int iter = 0; iter < 30; ++iter) { + // Per bin: the terms of var = a*s2 + b^2*I2, each divided by the current model's variance, and + // the median normalised squared deviation, expressed as the sum it stands for. Reordering inside + // a bin does not change which samples it holds. One bin per task, so the result does not depend + // on threads. + ParallelChunks(n_bins, ThreadsForWork(smp.size(), nthreads, 8 * 32768), [&](int blo, int bhi) { + for (int bin = blo; bin < bhi; ++bin) { + double s = 0.0, I = 0.0; + for (size_t i = bin_lo(bin); i < bin_hi(bin); ++i) { + auto &x = smp[i]; + const double v = a * x.s2 + b2 * x.I2; + s += x.s2 / v; I += x.I2 / v; + x.z2 = x.dev2 / v; + } + const size_t n = bin_hi(bin) - bin_lo(bin); + const auto mid = smp.begin() + bin_lo(bin) + n / 2; + std::nth_element(smp.begin() + bin_lo(bin), mid, smp.begin() + bin_hi(bin), + [](const ErrorModelBinned &p, const ErrorModelBinned &q) { return p.z2 < q.z2; }); + bs2[bin] = s; bI2[bin] = I; bd2[bin] = mid->z2 / CHI2_1_MEDIAN * static_cast(n); + } + }); + + // Every bin weighs the same in RELATIVE terms, so the weak bins still constrain a. + double Ass = 0, AsI = 0, AII = 0, Bs = 0, BI = 0; + for (int bin = 0; bin < n_bins; ++bin) { + if (!(bd2[bin] > 0.0)) continue; + const double w = 1.0 / (bd2[bin] * bd2[bin]); + Ass += w * bs2[bin] * bs2[bin]; AsI += w * bs2[bin] * bI2[bin]; AII += w * bI2[bin] * bI2[bin]; + Bs += w * bs2[bin] * bd2[bin]; BI += w * bI2[bin] * bd2[bin]; + } + const double det = Ass * AII - AsI * AsI; + double a_new, b2_new; + if (!fit.b_measured) { + if (!(Ass > 0.0)) return fit; + a_new = std::clamp(Bs / Ass, 0.25, 100.0); + b2_new = 0.0; + } else { + if (!(det > 1e-10 * Ass * AII)) return fit; + const double a_fit = (Bs * AII - BI * AsI) / det; + const double b2_fit = (Ass * BI - AsI * Bs) / det; + a_new = std::clamp(a_fit, 0.25, 100.0); + b2_new = std::max(b2_fit, 0.0); + // Leverage is not enough: b^2 can have it and still come out at zero within its own error, + // and 1/b then reads the noise in b as an I/sigma. The standard error of b^2 comes from the + // bins' own scatter about the fitted line; this only decides what is printed. + double chi = 0.0; + for (int bin = 0; bin < n_bins; ++bin) { + if (!(bd2[bin] > 0.0)) continue; + const double r = (bd2[bin] - a_fit * bs2[bin] - b2_fit * bI2[bin]) / bd2[bin]; + chi += r * r; + } + fit.b_resolved = b2_fit >= 2.0 * std::sqrt(chi / (n_bins - 2) * Ass / det); + } + const bool settled = std::fabs(a_new - a) <= 1e-4 * a && std::fabs(b2_new - b2) <= 1e-4 * b2; + a = a_new; + b2 = b2_new; + fit.active = true; + fit.a = a; + fit.b2 = b2; + if (settled) break; + } + return fit; +} diff --git a/image_analysis/scale_merge/ErrorModel.h b/image_analysis/scale_merge/ErrorModel.h new file mode 100644 index 000000000..89716bf9f --- /dev/null +++ b/image_analysis/scale_merge/ErrorModel.h @@ -0,0 +1,48 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// The error model of the rotation merge: var = a * s2 + b^2 * ^2, fitted so that the scatter of +// symmetry equivalents matches it - in each bin of counting I/sigma, the median of dev2 / var is the +// median of a chi^2 with one degree of freedom. +// +// One sample per observation: s2 is its counting variance as the merge rebuilds it at the reflection's +// mean intensity, I2 that mean squared, dev2 its leverage-corrected squared deviation from the mean. +// +// Two choices, each against a measured failure of the fit this replaces, which took the medians of s2, +// I2 and dev2 separately over bins of I2: +// * Bins are equal counts in COUNTING I/sigma (I2 / s2), not in I2. b is identified by how far the +// bins reach into the regime where b* dominates the counting term, and that is what this ranks +// on; ranking on I2 mixes weak high-resolution reflections into the strong bins. +// * The median is taken of the NORMALISED deviation, each sample divided by the variance the previous +// iteration gave it. A median of dev2 over observations whose variances differ is not 0.455 times +// their mean variance, and the ratio of three medians read a too low and b too high: on most of 28 +// sets the core of the normalised deviations then scattered at up to 1.8 times its stated variance +// in the weak and middle bins and at 0.1-0.7 times it in the strongest. +// A median, not a mean: on heavy-tailed data (split and powder-contaminated crystals) a mean calibrates +// the sigmas on the tails, the merge's outlier test widens with them, and CC1/2 fell 0.80 -> 0.71 on one +// such set where this fit raises it to 0.83. +struct ErrorModelSample { + double s2, I2, dev2; + float d; +}; + +struct ErrorModelFit { + bool active = false; // the pool was large enough to fit + double a = 1.0, b2 = 0.0; + bool b_measured = true; // false: no bin reached where b could be seen; b is held at 0 + bool b_resolved = true; // b^2 stood clear of two of its standard errors +}; + +// Scratch: the partitioned copy of the pool (reused between calls to save the allocation), with the rank +// key and the current normalised deviation. +struct ErrorModelBinned { + double snr2, s2, I2, dev2, z2; +}; + +ErrorModelFit FitErrorModel(const std::vector &pool, std::vector &scratch, + size_t nthreads); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index a558e74ee..b57d21126 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4047,8 +4047,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool double error_model_a = 1.0, error_model_b = 0.0; bool error_model_active = false; // What the reported chi2 is taken from, below: the pool the LAST successful fit was handed, with - // that fit's own a and b^2. fit_ab works on a partitioned copy of its pool and a median does not - // depend on the order, so the copy is not the thing to keep; b^2 is the fit's own because + // that fit's own a and b^2. FitErrorModel works on a partitioned copy of its pool and a median does + // not depend on the order, so the copy is not the thing to keep; b^2 is the fit's own because // error_model_b is its square root and squaring it back up is not the same number. const std::vector *chi2_pool = nullptr; double chi2_a = 0.0, chi2_b2 = 0.0; @@ -4056,159 +4056,29 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool bool error_model_b_resolved = true; // b^2 stood clear of its own standard error std::vector &samples = em_samples; samples.clear(); - 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. - std::vector &fit_scratch = em_fit_pool; + // A full's counting variance at intensity I_for_b. Its own sigma carries its own Poisson fluctuation: + // a full that came out low has a smaller sigma, so 1/sigma^2 weights it up and the merged mean drifts + // below . Rebuild the variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*, + // the linear model the combine measured - so the weight no longer knows this full's own fluctuation. + // Mirrors MergeOnTheFly::CorrectedSigma on the stills path. The error model is fitted on this same + // variance, so that its a scales what the merge applies it to. + auto counting_variance = [](const Obs &o, double I_for_b, double own_s2) { + const double base = static_cast(o.corr) * o.corr * o.var_bkg + + static_cast(o.corr) * o.var_per_I * std::max(0.0, I_for_b); + return base > 0.0 ? base : own_s2; + }; + // Fit (a, b) to a pool of samples (ErrorModel.h). A lambda because the pool changes once the cutoff + // below is known. + std::vector &fit_scratch = em_fit_pool; auto fit_error_model = [&](const std::vector &pool) { - // Works on a scratch copy, and the caller's pool keeps its order. That matters: the refit - // below builds its own pool by walking this one, and the partition is not a total order, so - // handing the refit a reordered pool moves which of two samples of equal I2 lands in which - // bin. One reused buffer rather than a fresh vector per call - the copy is a memcpy and costs - // a few milliseconds against the hundreds the sort it replaces used to. - auto fit_ab = [&](const std::vector &in) { - fit_scratch.assign(in.begin(), in.end()); - std::vector &smp = fit_scratch; - if (smp.size() < static_cast(8 * n_bins)) - return; - const size_t per = smp.size() / n_bins; - // The bins are equal COUNTS in I2 rank, and all that comes out of one is three medians - so - // the bin has to be the right SET, not a sorted one. Put each boundary in place with - // nth_element instead, splitting the boundaries down the middle so each level halves the - // range it works on: order n per level and four levels, against n log n for the sort. The - // medians are then taken off the bin's own span, by the field wanted, with no copy at all - - // the three vectors this used to build were the whole pool again, three times over, and - // median_of takes the lower median, which is what nth_element leaves at that index. - // Total, so which of two samples of equal I2 falls in which bin is a property of the - // samples rather than of the algorithm that partitioned them. The full sort this replaces - // was unstable, so it resolved those ties by whatever order it happened to receive - and - // the refit below is handed a different order from the first fit. Ordering on the other - // two fields costs nothing (they are already in the cache line) and settles it for good. - const auto by_I2 = [](const Sample &a, const Sample &b) { - if (a.I2 != b.I2) return a.I2 < b.I2; - if (a.s2 != b.s2) return a.s2 < b.s2; - if (a.dev2 != b.dev2) return a.dev2 < b.dev2; - return a.d < b.d; - }; - // The two halves a split leaves behind are disjoint ranges whose contents the nth_element - // above them has already fixed, so a child sees the same sequence whether it runs before, - // after or beside its sibling and the bins come out the same either way. Run them side by - // side while the range is still worth a thread; below that they stay on this one. - constexpr size_t SPLIT_MIN_PARALLEL = 1u << 17; - const std::function split = - [&](size_t lo, size_t hi, int b0, int b1, int depth) { - if (b0 >= b1) return; - const int mid = (b0 + b1) / 2; - const size_t k = static_cast(mid) * per; - std::nth_element(smp.begin() + lo, smp.begin() + k, smp.begin() + hi, by_I2); - if (depth > 0 && hi - lo > SPLIT_MIN_PARALLEL) { - auto left = std::async(std::launch::async, - [&, lo, k, b0, mid, depth] { split(lo, k, b0, mid, depth - 1); }); - split(k, hi, mid + 1, b1, depth - 1); - left.get(); - } else { - split(lo, k, b0, mid, 0); - split(k, hi, mid + 1, b1, 0); - } - }; - split(0, smp.size(), 1, n_bins, nthreads > 1 ? 3 : 0); - std::vector bs2(n_bins), bI2(n_bins), bd2(n_bins); - ParallelChunks(n_bins, ThreadsForWork(smp.size(), nthreads, 8 * 32768), [&](int blo, int bhi) { - for (int bin = blo; bin < bhi; ++bin) { - const size_t lo = static_cast(bin) * per; - const size_t hi = (bin == n_bins - 1) ? smp.size() : lo + per; - const auto mid = smp.begin() + lo + (hi - lo) / 2; - std::nth_element(smp.begin() + lo, mid, smp.begin() + hi, - [](const Sample &a, const Sample &b) { return a.s2 < b.s2; }); - bs2[bin] = mid->s2; - std::nth_element(smp.begin() + lo, mid, smp.begin() + hi, by_I2); - bI2[bin] = mid->I2; - std::nth_element(smp.begin() + lo, mid, smp.begin() + hi, - [](const Sample &a, const Sample &b) { return a.dev2 < b.dev2; }); - bd2[bin] = mid->dev2 / CHI2_1_MEDIAN; - } - }); - // `b` is identified ONLY by the spread of I^2/sigma^2 across the bins, and the bins hold equal - // COUNTS - so when fewer reflections are strong than one bin holds (1/16 of the pool), the top - // bin's median sits at an intensity where b cannot be measured, and the fit assigns it the - // bins' own noise-selection slope instead. Measured: a crystal with 2.2% of its fulls above - // I/sigma 2 reports b = 5.6, i.e. sigma -> 2*I at the strong end and ISa 0.50, while XDS gets - // 6.13 from the same images; synthetic data with b_true = 0 and 1.8% strong reproduces that to - // two digits, and recovers the truth as soon as the strong fraction passes 1/16. sigma^2 is - // then a*sigma^2 + (b*mean)^2 with a per-group constant, which caps merged |I/sigma| at - // sqrt(n)/b - reporting an impossible systematic error on data whose CC1/2 is 99%. - // - // So refuse to report what was not measured: if the strongest bin's own (I/sigma)^2 is small, - // fit `a` alone and leave b at zero. The threshold is not delicate - over the rotation battery - // the two crystals this fires on sit at 0.22 and 0.84 while the next is 31.7 and a healthy one - // is 342, so anything from 4 to 25 selects the same two and nothing else changes. - constexpr double B_LEVER_MIN = 4.0; - const bool b_measurable = bI2.back() > B_LEVER_MIN * bs2.back(); - std::vector bd2_sorted = bd2; - const double dev2_floor = std::max(1e-30, 1e-3 * median_of(bd2_sorted)); - double Ass = 0, AsI = 0, AII = 0, Bs = 0, BI = 0; - for (int bin = 0; bin < n_bins; ++bin) { - const double s2 = bs2[bin], I2 = bI2[bin], d2 = bd2[bin]; - const double d2w = std::max(d2, dev2_floor); - const double wgt = 1.0 / (d2w * d2w); - Ass += wgt * s2 * s2; AsI += wgt * s2 * I2; AII += wgt * I2 * I2; - Bs += wgt * s2 * d2; BI += wgt * I2 * d2; - } - const double det = Ass * AII - AsI * AsI; - if (!b_measurable) { - // One parameter, same weights: dev2 = a*sigma^2 alone. - if (Ass > 0.0) { - error_model_a = std::clamp(Bs / Ass, 0.25, 100.0); - error_model_b = 0.0; - error_model_b_unmeasured = true; - error_model_active = true; - chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = 0.0; - } - } else if (det > 1e-10 * Ass * AII) { - const double a_fit = (Bs * AII - BI * AsI) / det; - const double b2_fit = (Ass * BI - AsI * Bs) / det; - error_model_a = std::clamp(a_fit, 0.25, 100.0); - const double b2 = std::max(b2_fit, 0.0); - error_model_b = std::sqrt(b2); - error_model_b_unmeasured = false; - // Leverage is not enough: b^2 can have it and still come out at zero within its own error - - // the strong bins scatter no more than counting says, or there are too few samples for a - // bin median to mean anything. 1/b then reads the noise in b as an I/sigma, any value from - // ~30 to infinity between merges of the same data. The standard error of b^2 is taken from - // the bins' own scatter about the fitted line; the fit and the sigmas are left as they are, - // and so is the fitted ISa every decision reads - this only decides what is printed. - // Measured: z = 13.6 and 22 on two crystals whose b is plainly visible, 0.1-0.5 on a - // 3000-sample low-multiplicity set and -1.2 to -6 on one whose strongest bins scatter - // LESS than counting says (both printed ISa 96-130). - double chi = 0.0; - for (int bin = 0; bin < n_bins; ++bin) { - const double d2w = std::max(bd2[bin], dev2_floor); - const double r = bd2[bin] - a_fit * bs2[bin] - b2_fit * bI2[bin]; - chi += r * r / (d2w * d2w); - } - const double b2_se = std::sqrt(chi / (n_bins - 2) * Ass / det); - error_model_b_resolved = b2_fit >= 2.0 * b2_se; - error_model_active = true; - chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = b2; - } - }; - 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. - if (reject_outliers && error_model_active) { - const double ns2 = reject_nsigma * reject_nsigma, b2 = error_model_b * error_model_b; - std::vector &kept = em_refit_pool; - kept.clear(); - 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() < pool.size()) - fit_ab(kept); - } + const ErrorModelFit fit = FitErrorModel(pool, fit_scratch, nthreads); + if (!fit.active) return; + error_model_a = fit.a; + error_model_b = std::sqrt(fit.b2); + error_model_b_unmeasured = !fit.b_measured; + error_model_b_resolved = fit.b_resolved; + error_model_active = true; + chi2_pool = &pool; chi2_a = fit.a; chi2_b2 = fit.b2; }; // Friedel pair of each group when the hands are merged separately (-A); filled with the medians below. std::vector pair_of_group; @@ -4227,7 +4097,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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], mf.d[i]}); + if (gvalid[i]) + samples.push_back({counting_variance(fulls[i], em_mean[fulls[i].group], gs2[i]), + gI2[i], gdev2[i], mf.d[i]}); did_gpu = true; } #endif @@ -4267,7 +4139,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const double factor = 1.0 - w / (on_hand ? swh[s] : 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}); + samples.push_back({counting_variance(o, em_mean[o.group], s2), mean * mean, + resid * resid / factor, o.d}); } } @@ -4373,16 +4246,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // The full's sigma under the error model, with the variance evaluated at intensity I_for_b. auto model_sigma = [&](const Obs &o, double I_for_b, float sigma_corr) -> float { if (!error_model_active) return sigma_corr; - // A full's own sigma carries its own Poisson fluctuation: a full that came out low has a - // smaller sigma, so 1/sigma^2 weights it up and the merged mean drifts below . Rebuild the - // variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*, the linear model - // the combine measured - so the weight no longer knows this full's own fluctuation. Mirrors - // MergeOnTheFly::CorrectedSigma on the stills path. - double a_var = static_cast(sigma_corr) * sigma_corr; - const double base = static_cast(o.corr) * o.corr * o.var_bkg - + static_cast(o.corr) * o.var_per_I * std::max(0.0, I_for_b); - if (base > 0.0) - a_var = base; + const double a_var = counting_variance(o, I_for_b, static_cast(sigma_corr) * sigma_corr); const double v = error_model_a * a_var + (error_model_b * I_for_b) * (error_model_b * I_for_b); return v > 0.0 ? static_cast(std::sqrt(v)) : sigma_corr; @@ -4715,7 +4579,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const double I_corr = static_cast(mf.I[i]) * mf.corr[i]; const double sigma_corr = static_cast(mf.sigma[i]) * mf.corr[i]; auto &g = gs[gi]; - g.sum += I_corr; g.sum_sq += I_corr * I_corr; g.sum_var += sigma_corr * sigma_corr; ++g.n; + g.sum += I_corr; g.sum_sq += I_corr * I_corr; + g.sum_var += counting_variance(fulls[i], em_mean[gi], sigma_corr * sigma_corr); + ++g.n; } // Per-group counting-subtracted fractional systematic variance, paired with the group's I/sigma. // Two conventions decide whether this measures anything: diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 5f5e7a6e0..116378579 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -16,6 +16,7 @@ #include "../../common/UnitCell.h" #include "../IntegrationOutcome.h" +#include "ErrorModel.h" #include "Merge.h" // MergedReflection, MergeStatistics #ifdef JFJOCH_USE_CUDA #include @@ -144,10 +145,11 @@ private: int32_t group; // dense ASU-group id for the current space group; <0 = never mergeable }; - // One leverage-corrected error-model 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 over the samples that survive the automatic resolution cutoff. See MergeAndStats. - struct Sample { double s2, I2, dev2; float d; }; + // One leverage-corrected error-model sample per usable full: its counting variance at the group's + // mean, that mean squared, its squared deviation from the mean - and the resolution it sits at, + // because the fit is re-run over the samples that survive the automatic resolution cutoff. See + // MergeAndStats and ErrorModel.h. + using Sample = ErrorModelSample; // The narrow per-observation record the ingest sort orders: the raw hkl the runs are cut on, the // frame position that breaks a tie inside one, the observation's own index (which makes the order @@ -355,11 +357,12 @@ private: std::vector asymptote_scatter; // The error model's working pools (see MergeAndStats): the samples themselves, the scratch copy each - // fit partitions, the misfit-free subset the refit uses, the subset inside the resolution cutoff, and + // fit partitions, the subset inside the resolution cutoff, and // the per-sample chi2 the reported number is the median of. Members for the same reason as FullsStaging - one is 32 bytes per full and // there are two dozen fits per run, so as locals this is gigabytes of pages faulted in and handed // straight back. Every one of them is cleared and refilled before it is read. - std::vector em_samples, em_fit_pool, em_refit_pool, em_cut_pool; + std::vector em_samples, em_cut_pool; + std::vector em_fit_pool; std::vector em_chi2; // Set by FitPerFrameG: which frames were fitted this call (so corr/G is updated only there). diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 4794b6a7e..7be51a64e 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -4,6 +4,7 @@ #include #include +#include "../image_analysis/scale_merge/ErrorModel.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" @@ -406,3 +407,72 @@ TEST_CASE("ResolutionCutoff_RaggedFallOffIsReadOffTheBins") { REQUIRE(rc.d_cut); CHECK(*rc.d_cut > 2.30); // coarser than the band that correlates again } + +namespace { + // Samples of var = a*s2 + b^2*I2 over four decades of counting I/sigma, with counting variances + // spread over a factor of 100 at every intensity, and a fraction of gross outliers. + std::vector SyntheticErrorModelSamples(double a, double b, double max_snr, + double outlier_fraction) { + std::mt19937 rng(7); + std::uniform_real_distribution u(0.0, 1.0); + std::normal_distribution z(0.0, 1.0); + std::vector out; + for (int i = 0; i < 200000; ++i) { + const double s2 = std::pow(10.0, 2.0 * u(rng)); + const double snr = max_snr * std::pow(10.0, -4.0 * u(rng)); + const double I2 = snr * snr * s2; + const double e = z(rng); + double dev2 = (a * s2 + b * b * I2) * e * e; + if (u(rng) < outlier_fraction) dev2 *= 400.0; + out.push_back({s2, I2, dev2, 2.0f}); + } + return out; + } +} + +// The fit recovers a and b from observations whose counting variances differ widely inside every +// intensity bin (where a ratio of bin medians does not), with 0.3% gross outliers in the pool. +TEST_CASE("ErrorModel_RecoversAAndBThroughOutliers") { + const auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003); + std::vector scratch; + const auto fit = FitErrorModel(pool, scratch, 4); + REQUIRE(fit.active); + CHECK(fit.b_measured); + CHECK(fit.b_resolved); + CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03)); + CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05)); +} + +// A few observations whose counting variance is hugely overstated - their scatter is ordinary - do not +// set the fit: each sample counts by its deviation relative to its own variance. +TEST_CASE("ErrorModel_OverstatedCountingVarianceDoesNotSetTheFit") { + auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.0); + for (size_t i = 0; i < pool.size(); i += 200) + pool[i].s2 *= 1e4; + std::vector scratch; + const auto fit = FitErrorModel(pool, scratch, 4); + REQUIRE(fit.active); + CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03)); + CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05)); +} + +// The same pool on one thread and on many gives the same numbers. +TEST_CASE("ErrorModel_SameOnAnyNumberOfThreads") { + const auto pool = SyntheticErrorModelSamples(1.1, 0.05, 100.0, 0.001); + std::vector scratch; + const auto one = FitErrorModel(pool, scratch, 1); + const auto many = FitErrorModel(pool, scratch, 16); + CHECK(one.a == many.a); + CHECK(one.b2 == many.b2); +} + +// Data that never reach where b could be seen: a is fitted alone and b held at 0. +TEST_CASE("ErrorModel_BNotMeasuredOnWeakData") { + const auto pool = SyntheticErrorModelSamples(0.9, 0.05, 1.5, 0.0); + std::vector scratch; + const auto fit = FitErrorModel(pool, scratch, 4); + REQUIRE(fit.active); + CHECK(!fit.b_measured); + CHECK(fit.b2 == 0.0); + CHECK(fit.a == Catch::Approx(0.9).epsilon(0.03)); +}