From be750018339dc810b54cf2bda1b1f803e51bc04c Mon Sep 17 00:00:00 2001 From: jungfrau Date: Mon, 24 Aug 2026 05:32:12 -0400 Subject: [PATCH] Build the merge's per-frame quantities once per frame, and count what a shell can hold only when it is reported Three passes over the ingested observations were doing more than they needed. The smoothed-geometry pass rebuilt a CrystalLattice and its three reciprocal vectors for every observation, each one a cross product and a cell volume, for a value that depends only on which frame the observation came from. On a large sweep that is tens of millions of constructions against a couple of thousand distinct answers. The completeness column counts how many unique reflections a shell could hold. It is read off a merge that gets written out, never off the ones the space-group search runs on the way there - and those are the expensive ones to count, because the search merges in P1, where the list is the whole hemisphere rather than an asymmetric unit of it. The keep flags were written over the whole observation array as 1 and then immediately over it again as 0 whenever a resolution limit is set, which the default low-resolution limit always does. They are filled once now. The merged reflections also get their capacity up front rather than doubling their way to it several times per pass. Merged intensities are unchanged. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/RotationScaleMerge.cpp | 32 +++++++++++++++---- 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index d2448f8a..a14de6e9 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -472,9 +472,10 @@ void RotationScaleMerge::BuildInRangeObservations(const std::vector &ke const bool limited = d_min_limit.has_value() || d_max_limit.has_value(); std::vector keep_run(n_run, 1); - std::vector keep(n_obs, 1); + // Filled once, with the value the branch below actually wants: a resolution-limited run used to + // write the whole array as 1 and then immediately write it again as 0. + std::vector keep(n_obs, limited ? 0 : 1); if (limited) { - keep.assign(n_obs, 0); ParallelChunks(n_run, nthreads, [&](int lo, int hi) { for (int r = lo; r < hi; ++r) { const float d = rawrun_d[r]; @@ -733,6 +734,16 @@ void RotationScaleMerge::SmoothGeometry() { const Coord m3 = (m1 % m2).Normalize(); const float m2_S0 = m2 * S0, m3_S0 = m3 * S0; + // The reciprocal basis follows the frame's smoothed lattice, not the observation, so build it once + // per frame. Each of Astar/Bstar/Cstar is a cross product and a cell volume, and there are tens of + // millions of observations against a couple of thousand frames. + std::vector > frame_recip(n_frames); + for (int f = 0; f < n_frames; ++f) { + if (!ok[f]) continue; + const CrystalLattice sl(smoothed[f][0], smoothed[f][1], smoothed[f][2]); + frame_recip[f] = {sl.Astar(), sl.Bstar(), sl.Cstar()}; + } + std::atomic changed{0}; ParallelChunks(static_cast(partials.size()), nthreads, [&](int lo, int hi) { int64_t local = 0; @@ -740,9 +751,9 @@ void RotationScaleMerge::SmoothGeometry() { auto &o = partials[i]; if (o.frame < 0 || o.frame >= n_frames || !ok[o.frame]) continue; - const CrystalLattice sl(smoothed[o.frame][0], smoothed[o.frame][1], smoothed[o.frame][2]); - const Coord p0 = sl.Astar() * static_cast(o.h) + sl.Bstar() * static_cast(o.k) - + sl.Cstar() * static_cast(o.l); + const std::array &r = frame_recip[o.frame]; + const Coord p0 = r[0] * static_cast(o.h) + r[1] * static_cast(o.k) + + r[2] * static_cast(o.l); const float phi_near = o.delta_phi * static_cast(PI) / 180.0f; const auto phi = SolvePhi(p0, S0, m1, m2, m3, m2_S0, m3_S0, phi_near); if (!phi || !std::isfinite(*phi)) @@ -2857,6 +2868,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool auto run_merge = [&]() { 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 + // already holds each time. + result.merged.reserve(n_groups); merged_I.assign(n_groups, NAN); reject_count = 0; rejected_obs.assign(fulls.size(), 0); @@ -3023,7 +3038,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool }; std::vector sa(n_shells); std::vector possible(n_shells, 0); - if (reference_cell) + // Completeness is a reported quantity: it is read off the merge a pass writes out, never off the + // ones the space-group search runs on the way there. Those are also the expensive ones to count - + // the search merges in P1, where the list of possible reflections is the whole hemisphere rather + // than an asymmetric unit of it - so enumerating it there is the largest single piece of work in + // 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); for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s];