diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 34f3e437..6679fd6c 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -16,7 +16,7 @@ #include #include -#include +#include #include #include "HKLKey.h" @@ -2287,16 +2287,28 @@ void RotationScaleMerge::SmoothG(std::vector &obs, std::vector &g, // partials guard passed. Detection and remedy are both DropCollapsedScales': the frame's fulls leave the // merge. Merging them unscaled - the old behaviour - asserts G = 1 for a frame whose scale is // demonstrably not 1, which is worse than leaving it out altogether. -bool RotationScaleMerge::DropCollapsedFullScales() { +bool RotationScaleMerge::DropCollapsedFullScales(bool from_staging) { // The Unity model leaves corr = 1/G, and every full of a frame carries that frame's G. A full whose // frame was NOT fitted still carries the 1.0 the combine gave it, on both the host and the device // path - those are not measurements of anything and must stay out of the median, or a run with many // unfitted frames drags the median toward 1 and the floor with it. + // Frame and corr are the only two fields of the eighty this needs, and after a device scale-fulls + // they are still sitting in the staging arrays the download filled, in the same order - so read the + // eight bytes from there rather than striding the whole fulls array. The host scale-fulls path fits + // corr onto the fulls themselves and has no staging to read. std::vector g_frame(n_frames, NAN); - for (const auto &o : fulls) - if (o.frame >= 0 && o.frame < n_frames && std::isfinite(o.corr) && o.corr > 0.0f - && o.corr != 1.0f) - g_frame[o.frame] = 1.0 / static_cast(o.corr); + if (from_staging) + for (size_t i = 0; i < fulls.size(); ++i) { + const int f = fulls_staging.frame[i]; + const float c = fulls_staging.corr[i]; + if (f >= 0 && f < n_frames && std::isfinite(c) && c > 0.0f && c != 1.0f) + g_frame[f] = 1.0 / static_cast(c); + } + else + for (const auto &o : fulls) + if (o.frame >= 0 && o.frame < n_frames && std::isfinite(o.corr) && o.corr > 0.0f + && o.corr != 1.0f) + g_frame[o.frame] = 1.0 / static_cast(o.corr); std::vector fitted; fitted.reserve(g_frame.size()); @@ -2586,25 +2598,49 @@ void RotationScaleMerge::FinalizePerFrameScale(const std::vector &cc, co namespace { // Possible unique reflections per shell for the completeness column - mirrors CalcPossibleReflections. + // gemmi::for_all_reflections is inlined here rather than called: it visits every point of the + // (2h+1)(2k+1)(2l+1) box, which is tens of millions on a large cell at high resolution, and building + // a Miller vector of the survivors first is that many more megabytes for a column of ten integers. + // Inlined, the h planes are independent, so one thread takes each and they are summed back in h + // order - integer addition, so the counts are the serial ones whatever the split. void PossiblePerShell(int space_group_number, const UnitCell &cell, double d_min, double d_max, - const ResolutionShells &shells, bool merge_friedel, std::vector &possible) { - gemmi::UnitCell gemmi_cell = cell; + const ResolutionShells &shells, bool merge_friedel, std::vector &possible, + size_t nthreads) { + const gemmi::UnitCell gemmi_cell = cell; const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); if (sg == nullptr) return; - const std::vector hkls = gemmi::make_miller_vector(gemmi_cell, sg, d_min, d_max, true); const gemmi::GroupOps gops = sg->operations(); + const gemmi::ReciprocalAsu asu(sg); + const gemmi::Miller lim = gemmi_cell.get_hkl_limits(d_min); + const double inv_dmin2 = 1.0 / (d_min * d_min); + const double inv_dmax2 = d_max > 0 ? 1.0 / (d_max * d_max) : 0.0; CrystalLattice lattice(cell); const auto astar = lattice.Astar(), bstar = lattice.Bstar(), cstar = lattice.Cstar(); - for (const auto &hkl : hkls) { - const auto q = hkl[0] * astar + hkl[1] * bstar + hkl[2] * cstar; - const auto qlen = q.Length(); - if (qlen < 1e-6) continue; - const auto shell = shells.GetShell(1.0 / qlen); - if (!shell.has_value()) continue; - const int s = *shell; - if (s >= 0 && s < static_cast(possible.size())) - possible[s] += (merge_friedel || gops.is_reflection_centric(hkl)) ? 1 : 2; - } + const int n_shells = static_cast(possible.size()); + const int n_h = 2 * lim[0] + 1; + std::vector> per_h(n_h, std::vector(n_shells, 0)); + ParallelFor(n_h, std::min(nthreads, n_h), [&](int h_index) { + std::vector &p = per_h[h_index]; + gemmi::Miller hkl; + hkl[0] = -lim[0] + h_index; + for (hkl[1] = -lim[1]; hkl[1] <= lim[1]; ++hkl[1]) + for (hkl[2] = -lim[2]; hkl[2] <= lim[2]; ++hkl[2]) { + if (!asu.is_in(hkl)) continue; + const double inv_d2 = gemmi_cell.calculate_1_d2(hkl); + if (!(inv_d2 <= inv_dmin2 && inv_d2 > inv_dmax2)) continue; + if (gops.is_systematically_absent(hkl)) continue; + const auto q = hkl[0] * astar + hkl[1] * bstar + hkl[2] * cstar; + const auto qlen = q.Length(); + if (qlen < 1e-6) continue; + const auto shell = shells.GetShell(1.0 / qlen); + if (!shell.has_value()) continue; + const int sh = *shell; + if (sh >= 0 && sh < n_shells) + p[sh] += (merge_friedel || gops.is_reflection_centric(hkl)) ? 1 : 2; + } + }); + for (const auto &p : per_h) + for (int sh = 0; sh < n_shells; ++sh) possible[sh] += p[sh]; } } @@ -2649,8 +2685,14 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // ---- Error model: fit dev2 = a*sigma^2 + b^2*^2 from symmetry-equivalent scatter. ---- std::vector em_mean(n_groups, NAN); std::vector reject_median(n_groups, NAN); - double error_model_a = 1.0, error_model_b = 0.0, error_model_chi2 = 0.0; + 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 + // 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; bool error_model_b_unmeasured = false; // b had no leverage; ISa is not reported std::vector &samples = em_samples; samples.clear(); @@ -2760,14 +2802,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool error_model_b = 0.0; error_model_b_unmeasured = true; error_model_active = true; - std::vector &chi2 = em_chi2; - chi2.clear(); - chi2.reserve(smp.size()); - for (const auto &s : smp) { - const double v = error_model_a * s.s2; - if (v > 0.0) chi2.push_back(s.dev2 / v); - } - error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN; + chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = 0.0; } } else if (det > 1e-10 * Ass * AII) { error_model_a = std::clamp((Bs * AII - BI * AsI) / det, 0.25, 100.0); @@ -2775,14 +2810,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool error_model_b = std::sqrt(b2); error_model_b_unmeasured = false; error_model_active = true; - std::vector &chi2 = em_chi2; - chi2.clear(); - 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; + chi2_pool = ∈ chi2_a = error_model_a; chi2_b2 = b2; } }; fit_ab(pool); @@ -2879,79 +2907,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } 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 - // (a, b) fit above spans the whole intensity range, and a mild excess of scatter at intermediate - // 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. 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 gs(n_groups); - for (int i = 0; i < n_full; ++i) { - const int gi = mf.group[i]; - if (gi < 0) continue; - 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; - } - // 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> 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 = error_model_a * (g.sum_var / g.n); - if (mean <= 0.0 || counting <= 0.0) continue; - 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)}); - } - // 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 { - std::vector b2; - for (const auto &[b2_value, snr] : group_scatter) - if (snr >= snr_min) b2.push_back(b2_value); - return b2.size() >= min_groups ? std::optional(median_of(b2)) : std::nullopt; - }; - 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 = [&](const Obs &o, float I_corr, float sigma_corr) -> float { if (!error_model_active) return sigma_corr; @@ -2973,14 +2928,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // ---- 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; + std::vector &acc = merge_acc; Result result; 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{}); + // The device path's unpack below writes every field of every group - MergeAccumKernel fills the + // empty ones too - so there this only has to be the right size. The host path accumulates into + // it with += and does need the zeroes. + if (use_gpu_merge) acc.resize(n_groups); + else acc.assign(n_groups, Accum{}); result.merged.clear(); // One reflection per group at most, and a merge runs several times per pass. Without this the // vector doubles its way up to a few hundred megabytes on a large cell, copying everything it @@ -2992,22 +2950,33 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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); + // Resized, not sized vectors: the kernel writes every one of these for every group, so + // value-initialising them first is a hundred megabytes of zeroes faulted in on one thread + // and overwritten on the next line - the same trap Ingest calls out. Members for the same + // reason as FullsStaging. + MergeAccumStaging &ms = merge_accum; + ms.Resize(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]); - } + ms.swI.data(), ms.sw.data(), ms.swIh0.data(), ms.swIh1.data(), + ms.swh0.data(), ms.swh1.data(), ms.nh0.data(), ms.nh1.data(), + ms.d.data(), ms.rej.data(), rejected_obs.data()); + // Per-group and independent, so the unpack splits; the rejected count is an integer sum and + // does not care in which order the chunks add theirs in. + std::atomic rejected_total{0}; + ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int glo, int ghi) { + size_t rej = 0; + for (int g = glo; g < ghi; ++g) { + Accum &a = acc[g]; + a.swI = ms.swI[g]; a.sw = ms.sw[g]; a.swIh[0] = ms.swIh0[g]; a.swIh[1] = ms.swIh1[g]; + a.swh[0] = ms.swh0[g]; a.swh[1] = ms.swh1[g]; + a.nh[0] = static_cast(ms.nh0[g]); a.nh[1] = static_cast(ms.nh1[g]); + a.d = static_cast(ms.d[g]); + rej += static_cast(ms.rej[g]); + } + rejected_total += rej; + }); + reject_count = rejected_total; did_gpu_acc = true; } #endif @@ -3084,13 +3053,86 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool 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; }); } + // 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 + // (a, b) fit above spans the whole intensity range, and a mild excess of scatter at intermediate + // 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. Report-only: it is the reported ISa, nothing downstream uses it. + // Host-side over the merged fulls, so CPU and GPU agree. Taken once, on the error model as it + // finally stands: nothing between the fits reads it, so an estimate made before the resolution + // cutoff's refit is only overwritten, and it costs a pass over every full scattered into a + // per-group array as long as the group count. + double error_model_b_asymptotic = error_model_b; + if (error_model_active) { + std::vector &gs = asymptote_scatter; + gs.assign(n_groups, GroupScatter{}); + for (int i = 0; i < n_full; ++i) { + const int gi = mf.group[i]; + if (gi < 0) continue; + 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; + } + // 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> 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 = error_model_a * (g.sum_var / g.n); + if (mean <= 0.0 || counting <= 0.0) continue; + 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)}); + } + // 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 { + std::vector b2; + for (const auto &[b2_value, snr] : group_scatter) + if (snr >= snr_min) b2.push_back(b2_value); + return b2.size() >= min_groups ? std::optional(median_of(b2)) : std::nullopt; + }; + 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); + } + // 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 @@ -3103,6 +3145,20 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool result.error_model_a = em.a; result.error_model_b = em.b; } + // Report-only, and only the LAST successful fit's chi2 is ever printed - so its median is taken + // here, once, instead of in every fit: the pool is one sample per full and the median is a pass over + // it plus an nth_element, ten times a run. + double error_model_chi2 = 0.0; + if (chi2_pool) { + std::vector &chi2 = em_chi2; + chi2.clear(); + chi2.reserve(chi2_pool->size()); + for (const auto &s : *chi2_pool) { + const double v = chi2_a * s.s2 + chi2_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; + } if (error_model_active) { // Reported in XDS's convention so the numbers can be read against a CORRECT.LP directly. // Two ISa are printed and they are different quantities: the whole-range 1/sqrt(a*b), which @@ -3160,7 +3216,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // the merge that nothing goes on to read. if (reference_cell && !for_search) PossiblePerShell(x.GetSpaceGroupNumber().value_or(1), *reference_cell, d_min_pad, d_max_pad, - shells, merge_friedel, possible); + shells, merge_friedel, possible, nthreads); for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; CorrelationCoefficient cc_half_overall; @@ -3695,7 +3751,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st } logger.Info("Scaled fulls (XDS order, Unity model)"); } - const bool rejected_full_scales = scale_fulls && DropCollapsedFullScales(); + const bool rejected_full_scales = + scale_fulls && DropCollapsedFullScales(combined_on_gpu && scaled_fulls_on_gpu); // --- 4b. Optional correction surfaces (decay = resolution x time; absorption = goniometer-frame // diffracted-beam direction), each an alternating multiplicative fit of the fulls' corr against diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 2974556b..90e46cf8 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include #include @@ -220,6 +221,43 @@ private: }; MergeFields merge_fields; + // One host array per field for the fulls download: the device hands back an array per field and the + // host gathers them into `fulls`. Members rather than locals in Run() because the whole + // scale->combine->merge chain runs several times per run and these are a few hundred megabytes + // between them, so as locals every chain allocates, faults in and zeroes the lot again. + struct FullsStaging { + std::vector h, k, l, frame, group; + std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I; + std::vector on_ice; + void Resize(int n) { + h.resize(n); k.resize(n); l.resize(n); frame.resize(n); group.resize(n); + I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); + px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); + on_ice.resize(n); + } + }; + FullsStaging fulls_staging; + + // The merge accumulators (see MergeAndStats' run_merge): one entry per ASU group, plus the arrays + // the device kernel fills that the host unpacks into them. Members for the same reason as + // FullsStaging - a merge runs several times per run and this is a hundred megabytes between them. + 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 merge_acc; + struct MergeAccumStaging { + std::vector swI, sw, swIh0, swIh1, swh0, swh1, d; + std::vector nh0, nh1, rej; + void Resize(int n) { + swI.resize(n); sw.resize(n); swIh0.resize(n); swIh1.resize(n); swh0.resize(n); swh1.resize(n); + d.resize(n); nh0.resize(n); nh1.resize(n); rej.resize(n); + } + }; + MergeAccumStaging merge_accum; + + // Per-group scatter for the strong-reflection ISa asymptote (see MergeAndStats). A member for the + // same reason: 32 bytes a group, once per merge. + struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; }; + 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 // 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 @@ -254,23 +292,6 @@ private: // the CPU loops as the bit-parity fallback. Built in Ingest. std::unique_ptr gpu_; bool gpu_active_ = false; - - // One host array per field for the fulls download: the device hands back an array per field and the - // host gathers them into `fulls`. Members rather than locals in Run() because the whole - // scale->combine->merge chain runs several times per run and these are a few hundred megabytes - // between them, so as locals every chain allocates, faults in and zeroes the lot again. - struct FullsStaging { - std::vector h, k, l, frame, group; - std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I; - std::vector on_ice; - void Resize(int n) { - h.resize(n); k.resize(n); l.resize(n); frame.resize(n); group.resize(n); - I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); - px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); - on_ice.resize(n); - } - }; - FullsStaging fulls_staging; #endif // --- helpers (each a flat pass; see the .cpp) --- @@ -334,9 +355,10 @@ private: // Drop the fulls of any frame whose scale collapsed toward zero. The fulls are scaled with the Unity // model, so their corr IS 1/G and a collapsed G multiplies every intensity on that frame without - // bound. Reads the host fulls, so it covers the CPU and GPU scaling paths alike. Returns true if + // bound. Covers the CPU and GPU scaling paths alike; `from_staging` says the fulls' frame and corr + // are still in fulls_staging, which is where the scan reads them from when they are. Returns true if // anything was dropped (the caller then has to push the corrected corr back to the device). - bool DropCollapsedFullScales(); + bool DropCollapsedFullScales(bool from_staging); // Post-scale-fulls correction surfaces, each an alternating multiplicative fit of the host fulls' corr // against the merged reference (cheap host loops; the corrected corr is re-uploaded to the resident diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 21b3fe61..9e6bc3e9 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -562,8 +562,8 @@ SearchSpaceGroupResult SearchSpaceGroup( std::string refused_pg_hm, refused_why; std::vector pg_cands; double chi2_ref = std::numeric_limits::infinity(); - // Which point groups their own operators confirm. Serial: this fills the operator cache, and it - // is only a handful of cache lookups per group either way. + // Which point groups their own operators confirm. Serial: the operators were all scored in + // parallel above, so this is a handful of cache lookups per group. std::vector confirmed; std::vector confirmed_cc; for (const auto& pg : point_groups) {