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];