From 7762d8bd30f6fb75ca427406459fd12b07951e1a Mon Sep 17 00:00:00 2001 From: jungfrau Date: Sun, 23 Aug 2026 14:04:21 -0400 Subject: [PATCH] Build an observation only if the resolution range keeps it, and let the pre-pass skip what it discards Two changes to what the scaling stage does at all, rather than to how fast it does it. Ingest converted every integrated partial into the eighty-byte record the scaling works on, and then threw away whatever fell outside the requested resolution range. On a large cell that is sixty-three million records built and fifty-seven million discarded - five gigabytes written, most of it to be skipped by every consumer afterwards. It now emits a twenty-four byte key per observation, sorts and buckets those, decides from the runs which raw hkl the range keeps, and builds the full record only for the survivors. The flux meter still sums a frame's whole background in that frame's own order on one thread, because it is the number every merged intensity is divided by; the first usable d is still taken over the run rather than over the survivors; and the sort order was already total, so how the keys are filled cannot change it. The other is the pre-pass. It exists to choose a space group and post-refine the geometry, and its merged intensities are discarded - the second pass makes them again at the refined geometry. It was nonetheless fitting the decay, absorption and modulation surfaces, measuring radiation damage and sweep quality, assigning R-free flags, converting to amplitudes, walking the observations again for R_meas, splitting the anomalous pairs and analysing twinning, all for a result nobody reads. A flag threaded from the call site turns that off on the pre-pass, following the convention the anomalous split already used. What the pre-pass keeps is what is read later: the merge itself, the error model, the resolution cutoff, and the whole per-shell statistics block - because the second pass is judged against the first, and that guard needs the pre-pass's completeness and CC1/2. The flag that says the statistics exist is untouched and still set unconditionally; moving it is what disabled the guard entirely in an earlier attempt at this, and with it the completeness bound, the CC1/2 bound, the lattice-conflict test and the fall back to the header geometry. One consequence to be aware of: the pre-pass's CC1/2 is now measured without the correction surfaces while the second pass's is measured with them, so the two are no longer compared on quite the same footing. It moves the guard in the direction of firing less readily, never more, so it cannot roll back a good pass - but it is a small loss of sensitivity and the next commit removes it. Merged output byte-identical on a large-cell set, a high-multiplicity one and a small one; battery 6m28s against 7m55s, space group unchanged on all 24 crystals. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011n8riB6X59oRjkrSHzNPAU --- image_analysis/geom_refinement/PostRefine.cpp | 9 +- .../scale_merge/RotationScaleMerge.cpp | 441 +++++++++--------- .../scale_merge/RotationScaleMerge.h | 42 +- rugnux/Rugnux.cpp | 19 +- rugnux/rugnux_cli.cpp | 2 +- 5 files changed, 266 insertions(+), 247 deletions(-) diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 2af4a9d6..cedc68a8 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -637,11 +637,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorh, pp->k, pp->l, sys)), new ceres::CauchyLoss(0.02), beam, dist, const_cast(det_rot), rot_vec, p0, p1, p2); - p.SetParameterBlockConstant(const_cast(det_rot)); - p.SetParameterBlockConstant(rot_vec); - p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); } if (p.NumResidualBlocks() == 0) { beam_out[0] = beam_x0; beam_out[1] = beam_y0; dist_out = dist0; return false; } + // Everything but the beam and the distance is held at its step-A value. Once per problem: + // the block only has to exist, and repeating it per observation was up to MAX_OBS times + // five calls for the same five blocks. + p.SetParameterBlockConstant(const_cast(det_rot)); + p.SetParameterBlockConstant(rot_vec); + p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05); for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0); p.SetParameterUpperBound(beam, j, beam[j] + 15.0); } diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 42454be5..100822ab 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -231,23 +231,18 @@ void RotationScaleMerge::Ingest() { size_t total = 0; for (const auto &o : partials_out) total += o.reflections.size(); - partials.clear(); - partials.resize(total); - frame_start.assign(n_frames, 0); - frame_count.assign(n_frames, 0); frame_cell_ok.assign(n_frames, 1); g_partial.assign(n_frames, 1.0); - // The per-frame CSR up front. Every reflection of a frame becomes one observation, so where a - // frame's slice of `partials` starts is known before anything is copied - which is what lets the - // conversion below hand a frame to each thread instead of appending, frame after frame, on one. - // (Hence resize rather than reserve above: the slices are written, not pushed.) + // Where each frame's slice of the key array starts. Every reflection of a frame becomes one key, so + // this is known before anything is read - which is what lets the sweep below hand a frame to each + // thread instead of appending, frame after frame, on one. + std::vector src_start(n_frames, 0); { int32_t at = 0; for (int o = 0; o < n_frames; ++o) { - frame_start[o] = at; - frame_count[o] = static_cast(partials_out[o].reflections.size()); - at += frame_count[o]; + src_start[o] = at; + at += static_cast(partials_out[o].reflections.size()); } } @@ -259,6 +254,15 @@ void RotationScaleMerge::Ingest() { // which is what keeps the mean bit-exact. std::vector mean_bkg(n_frames, NAN); + // The h range comes out of the same sweep - the sort below buckets by h and needs to know how many + // buckets that is - and so does the narrow key of every observation. The fat Obs is NOT built here: + // the resolution range test below throws most of it away on a crystal that integrates far past what + // it merges, and eighty bytes written for an observation nothing can ever merge is eighty bytes of + // pure waste. Only the survivors are built, from the source reflections, once the runs are known. + constexpr int INT_LO = std::numeric_limits::min(); + constexpr int INT_HI = std::numeric_limits::max(); + std::vector unsorted(total); + std::vector h_lo(n_frames, INT_HI), h_hi(n_frames, INT_LO); ParallelFor(n_frames, nthreads, [&](int o) { if (reference_cell) { const auto cell = partials_out[o].latt.GetUnitCell(); @@ -266,68 +270,33 @@ void RotationScaleMerge::Ingest() { } double bkg_sum = 0.0; int bkg_n = 0; - int32_t at = frame_start[o]; + int lmin = INT_HI, lmax = INT_LO; + int32_t at = src_start[o]; for (const auto &r : partials_out[o].reflections) { if (std::isfinite(r.bkg) && r.bkg > 0.0f) { bkg_sum += r.bkg; ++bkg_n; } - Obs &obs = partials[at++]; - obs.h = r.h; obs.k = r.k; obs.l = r.l; - obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; obs.rlp = r.rlp; - obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; - obs.px = r.predicted_x; obs.py = r.predicted_y; - obs.image_number = r.image_number; - obs.frame = o; - obs.on_ice = r.on_ice_ring ? 1 : 0; - obs.corr = r.image_scale_corr; - obs.group = -1; + unsorted[at] = SortKey{r.h, r.k, r.l, r.image_number, at, r.d}; + lmin = std::min(lmin, r.h); + lmax = std::max(lmax, r.h); + ++at; } if (bkg_n > 0) mean_bkg[o] = bkg_sum / bkg_n; + h_lo[o] = lmin; h_hi[o] = lmax; }); MeasureIncidentFlux(mean_bkg); - // Per-obs AcceptReflection finiteness (immutable) - lets ComputeAsuGroups stamp the ASU-group id per - // obs from a flat 1-byte array instead of re-reading the fat Obs struct for every space group. - finite_ok.resize(partials.size()); - // The h range comes out of the same sweep - the sort below buckets by h and needs to know how - // many buckets that is, and this pass already reads every observation. int h_min = 0, h_max = 0; { - const int n_obs = static_cast(partials.size()); - const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n_obs))); - constexpr int INT_LO = std::numeric_limits::min(); - constexpr int INT_HI = std::numeric_limits::max(); - std::vector lo_of(nt, INT_HI), hi_of(nt, INT_LO); - const int chunk = (n_obs + nt - 1) / nt; - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int lmin = INT_HI, lmax = INT_LO; - const int lo = t * chunk, hi = std::min(n_obs, lo + chunk); - for (int i = lo; i < hi; ++i) { - auto &o = partials[i]; - // Divide the incident flux out of rlp here: this pass already has the observation - // in hand, and rlp is the numerator of corr, so the fits below see only what is - // left. Before the finiteness test, which is where it has always been. - o.rlp = static_cast(o.rlp / frame_flux[o.frame]); - finite_ok[i] = (std::isfinite(o.I) && std::isfinite(o.rlp) && o.rlp != 0.0f - && std::isfinite(o.sigma) && o.sigma > 0.0f) ? 1 : 0; - lmin = std::min(lmin, o.h); - lmax = std::max(lmax, o.h); - } - lo_of[t] = lmin; hi_of[t] = lmax; - } - }); int gmin = INT_HI, gmax = INT_LO; - for (int t = 0; t < nt; ++t) { - gmin = std::min(gmin, lo_of[t]); - gmax = std::max(gmax, hi_of[t]); + for (int o = 0; o < n_frames; ++o) { + gmin = std::min(gmin, h_lo[o]); + gmax = std::max(gmax, h_hi[o]); } if (gmin <= gmax) { h_min = gmin; h_max = gmax; } // otherwise there are no observations } // Sort ONCE by (raw h,k,l, image_number) and split into raw-hkl runs. This is the one expensive sort; // both the 3D combine (event split) and the per-space-group ASU grouping reuse this order. - // Sorting an index array whose comparator dereferences the 80-byte Obs is a cache miss per - // comparison over a multi-GB array, so sort a packed copy of the key fields instead. // // The observation's own index is the last key, which makes the order TOTAL. Two observations can // genuinely share (h,k,l) and image_number: the predictor emits both intersections of a @@ -337,7 +306,6 @@ void RotationScaleMerge::Ingest() { // member of an event, and its per-event sums are floating point. Ordering them by arrival makes // the result independent of the sort algorithm, which is what lets the order be reproduced by a // faster one. - struct SortKey { int32_t h, k, l; float image_number; int32_t idx; }; const auto key_less = [](const SortKey &a, const SortKey &b) { if (a.h != b.h) return a.h < b.h; if (a.k != b.k) return a.k < b.k; @@ -345,16 +313,16 @@ void RotationScaleMerge::Ingest() { if (a.image_number != b.image_number) return a.image_number < b.image_number; return a.idx < b.idx; }; - perm.resize(partials.size()); + // Two key arrays for a moment: the sweep above wrote them in index order, and where a key belongs + // in the bucketed order is not known until the whole histogram is. Both together are less than half + // of what the observation array they replace used to occupy on its own. + std::vector keys(total); { // Bucket by h first, then sort each bucket. h is the comparator's leading key, so the sorted // array is exactly the buckets laid end to end - and because the order is total (see above) // the permutation is the one a single std::sort produces, whichever algorithm gets there. // That is what makes this checkable: the output is byte-identical, not merely equivalent. - // - // The keys are built straight into their bucket slot, so this replaces the build pass rather - // than adding to it and the 20-bytes-per-observation array is never duplicated. - const int n = static_cast(partials.size()); + const int n = static_cast(total); const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n))); const int chunk = (n + nt - 1) / nt; const int H = h_max - h_min + 1; @@ -363,7 +331,7 @@ void RotationScaleMerge::Ingest() { ParallelChunks(nt, nthreads, [&](int tlo, int thi) { for (int t = tlo; t < thi; ++t) { const int lo = t * chunk, hi = std::min(n, lo + chunk); - for (int i = lo; i < hi; ++i) hist[t][partials[i].h - h_min]++; + for (int i = lo; i < hi; ++i) hist[t][unsorted[i].h - h_min]++; } }); @@ -381,18 +349,15 @@ void RotationScaleMerge::Ingest() { } bstart[H] = acc; - std::vector keys(partials.size()); ParallelChunks(nt, nthreads, [&](int tlo, int thi) { for (int t = tlo; t < thi; ++t) { std::vector fill = hist[t]; const int lo = t * chunk, hi = std::min(n, lo + chunk); - for (int i = lo; i < hi; ++i) { - const auto &o = partials[i]; - keys[fill[o.h - h_min]++] = - SortKey{o.h, o.k, o.l, o.image_number, static_cast(i)}; - } + for (int i = lo; i < hi; ++i) keys[fill[unsorted[i].h - h_min]++] = unsorted[i]; } }); + unsorted.clear(); + unsorted.shrink_to_fit(); // Largest bucket first: the tail of this is one bucket, so it should be the big one. std::vector order(H); @@ -405,10 +370,6 @@ void RotationScaleMerge::Ingest() { std::sort(keys.begin() + bstart[b], keys.begin() + bstart[b + 1], key_less); }); - ParallelChunks(n, nthreads, [&](int lo, int hi) { - for (int i = lo; i < hi; ++i) perm[i] = keys[i].idx; - }); - // The runs never cross an h boundary, so each bucket owns its own and they can be counted, // scanned and written without touching each other. Sizing the arrays from the count also // removes the repeated growth the push_backs paid for. @@ -440,7 +401,7 @@ void RotationScaleMerge::Ingest() { while (j < bstart[b + 1]) { if (keys[j].h != k0.h || keys[j].k != k0.k || keys[j].l != k0.l) break; if (!std::isfinite(d)) { - const float dj = partials[keys[j].idx].d; // resolution: only until one is usable + const float dj = keys[j].d; // resolution: only until one is usable if (std::isfinite(dj) && dj > 0.0f) d = dj; } ++j; @@ -454,7 +415,7 @@ void RotationScaleMerge::Ingest() { } }); } - DropOutOfRangeObservations(); + BuildInRangeObservations(keys); rawrun_group.assign(rawrun_start.size(), -1); logger.Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)", partials.size(), n_frames, rawrun_start.size()); @@ -497,84 +458,95 @@ void RotationScaleMerge::Ingest() { #endif } -// Compact `partials` (and the perm / raw-hkl tables that index it) down to the observations whose -// resolution is inside the requested range. See the header for why this cannot change any result: a raw -// hkl outside the range is group -1 in every pass, and dropping WHOLE runs on the same per-hkl resolution -// ComputeAsuGroups tests leaves rawrun_d, the surviving observations and their order untouched. +// Build `partials` (and the per-frame CSR, the finiteness mask and `perm`) from the source reflections, +// skipping the observations whose resolution is outside the requested range. See the header for why this +// cannot change any result: a raw hkl outside the range is group -1 in every pass, and skipping WHOLE +// runs on the same per-hkl resolution ComputeAsuGroups tests leaves rawrun_d, the surviving observations +// and their order untouched. // // It has to come after the incident flux has been measured (the flux meter is the frame's mean background // over ALL its reflections) and after the sort (which is what defines a raw-hkl run in the first place). -void RotationScaleMerge::DropOutOfRangeObservations() { - if (!d_min_limit && !d_max_limit) - return; - +void RotationScaleMerge::BuildInRangeObservations(const std::vector &keys) { const int n_run = static_cast(rawrun_start.size()); - const int n_obs = static_cast(partials.size()); + const int n_obs = static_cast(keys.size()); + const bool limited = d_min_limit.has_value() || d_max_limit.has_value(); - std::vector keep_run(n_run, 0); - std::vector keep(n_obs, 0); - ParallelChunks(n_run, nthreads, [&](int lo, int hi) { - for (int r = lo; r < hi; ++r) { - const float d = rawrun_d[r]; - if (!std::isfinite(d) || d <= 0.0f) continue; - if (d_min_limit && d < *d_min_limit) continue; - if (d_max_limit && d > *d_max_limit) continue; - keep_run[r] = 1; - const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r]; - for (int p = p0; p < p1; ++p) keep[perm[p]] = 1; // distinct runs own disjoint observations - } - }); - - // New index of every kept observation. Compacting in index order keeps each frame's observations - // contiguous and in the order they had, so the per-frame ranges still describe them and every - // per-frame sum keeps its sequence of roundings. - const int nt = static_cast(std::clamp(nthreads, 1, std::max(1, n_obs))); - const int chunk = (n_obs + nt - 1) / nt; - std::vector base(nt + 1, 0); - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int32_t c = 0; - for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) c += keep[i]; - base[t + 1] = c; - } - }); - for (int t = 0; t < nt; ++t) base[t + 1] += base[t]; - const int32_t n_keep = base[nt]; - if (n_keep == n_obs) - return; - - std::vector new_idx(n_obs); - std::vector kept_obs(n_keep); - std::vector kept_finite(n_keep); - ParallelChunks(nt, nthreads, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - int32_t at = base[t]; - for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) { - if (!keep[i]) continue; - new_idx[i] = at; - kept_obs[at] = partials[i]; - kept_finite[at] = finite_ok[i]; - ++at; + std::vector keep_run(n_run, 1); + std::vector keep(n_obs, 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]; + const bool in_range = std::isfinite(d) && d > 0.0f + && !(d_min_limit && d < *d_min_limit) && !(d_max_limit && d > *d_max_limit); + keep_run[r] = in_range ? 1 : 0; + if (!in_range) continue; + const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r]; + for (int p = p0; p < p1; ++p) keep[keys[p].idx] = 1; // distinct runs own disjoint observations } - } - }); - partials.swap(kept_obs); - finite_ok.swap(kept_finite); - kept_obs.clear(); kept_obs.shrink_to_fit(); - - frame_start.assign(n_frames, 0); - frame_count.assign(n_frames, 0); - for (int i = 0; i < n_keep; ) { - const int f = partials[i].frame; - int j = i; - while (j < n_keep && partials[j].frame == f) ++j; - frame_start[f] = i; - frame_count[f] = j - i; - i = j; + }); } - // The surviving runs, in the same order, each owning the same perm entries remapped to the new - // observation indices. + // Where each frame's reflections sit in the key numbering, and how many of them survive. Building in + // index order keeps each frame's observations contiguous and in the order they had, so the per-frame + // ranges still describe them and every per-frame sum keeps its sequence of roundings. + std::vector src_start(n_frames, 0); + { + int32_t at = 0; + for (int o = 0; o < n_frames; ++o) { + src_start[o] = at; + at += static_cast(partials_out[o].reflections.size()); + } + } + frame_count.assign(n_frames, 0); + ParallelFor(n_frames, nthreads, [&](int o) { + int32_t c = 0; + const int32_t n = static_cast(partials_out[o].reflections.size()); + for (int32_t j = 0; j < n; ++j) c += keep[src_start[o] + j]; + frame_count[o] = c; + }); + frame_start.assign(n_frames, 0); + int32_t n_keep = 0; + for (int o = 0; o < n_frames; ++o) { + frame_start[o] = n_keep; + n_keep += frame_count[o]; + } + + // The observations themselves, one frame per thread. finite_ok is the per-obs AcceptReflection + // finiteness (immutable) - it lets ComputeAsuGroups stamp the ASU-group id per obs from a flat + // 1-byte array instead of re-reading the fat Obs struct for every space group. + partials.clear(); + partials.resize(n_keep); + finite_ok.assign(n_keep, 0); + std::vector new_idx(n_obs); // only written where keep[i] + ParallelFor(n_frames, nthreads, [&](int o) { + int32_t src = src_start[o]; + int32_t at = frame_start[o]; + for (const auto &r : partials_out[o].reflections) { + if (!keep[src]) { ++src; continue; } + new_idx[src++] = at; + Obs &obs = partials[at]; + obs.h = r.h; obs.k = r.k; obs.l = r.l; + obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; + // Divide the incident flux out of rlp here: rlp is the numerator of corr, so the fits + // below see only what is left. + obs.rlp = static_cast(r.rlp / frame_flux[o]); + obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; + obs.px = r.predicted_x; obs.py = r.predicted_y; + obs.image_number = r.image_number; + obs.frame = o; + obs.on_ice = r.on_ice_ring ? 1 : 0; + obs.corr = r.image_scale_corr; + obs.group = -1; + finite_ok[at] = (std::isfinite(obs.I) && std::isfinite(obs.rlp) && obs.rlp != 0.0f + && std::isfinite(obs.sigma) && obs.sigma > 0.0f) ? 1 : 0; + ++at; + } + }); + + // The surviving runs, in the same order, each owning the same sorted keys remapped to the + // observation indices they were just built at. int n_keep_run = 0; std::vector dst_run(n_run, -1), dst_start(n_run, 0); { @@ -582,8 +554,8 @@ void RotationScaleMerge::DropOutOfRangeObservations() { for (int r = 0; r < n_run; ++r) if (keep_run[r]) { dst_run[r] = n_keep_run++; dst_start[r] = at; at += rawrun_count[r]; } } - std::vector new_perm(n_keep), ns(n_keep_run), nc(n_keep_run), - nh(n_keep_run), nk(n_keep_run), nl(n_keep_run); + perm.resize(n_keep); + std::vector ns(n_keep_run), nc(n_keep_run), nh(n_keep_run), nk(n_keep_run), nl(n_keep_run); std::vector nd(n_keep_run); ParallelChunks(n_run, nthreads, [&](int lo, int hi) { for (int r = lo; r < hi; ++r) { @@ -593,16 +565,16 @@ void RotationScaleMerge::DropOutOfRangeObservations() { ns[rr] = at; nc[rr] = rawrun_count[r]; nh[rr] = rawrun_h[r]; nk[rr] = rawrun_k[r]; nl[rr] = rawrun_l[r]; nd[rr] = rawrun_d[r]; for (int32_t p = 0; p < rawrun_count[r]; ++p) - new_perm[at + p] = new_idx[perm[rawrun_start[r] + p]]; + perm[at + p] = new_idx[keys[rawrun_start[r] + p].idx]; } }); - perm.swap(new_perm); rawrun_start.swap(ns); rawrun_count.swap(nc); rawrun_h.swap(nh); rawrun_k.swap(nk); rawrun_l.swap(nl); rawrun_d.swap(nd); - logger.Info("RotationScaleMerge: dropped {} of {} observations ({} of {} distinct hkl) outside the " - "requested resolution range - nothing downstream could have merged them", - n_obs - n_keep, n_obs, n_run - n_keep_run, n_run); + if (n_keep < n_obs) + logger.Info("RotationScaleMerge: dropped {} of {} observations ({} of {} distinct hkl) outside the " + "requested resolution range - nothing downstream could have merged them", + n_obs - n_keep, n_obs, n_run - n_keep_run, n_run); } namespace { @@ -2544,7 +2516,7 @@ namespace { } RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for_search, - bool fulls_resident) { + bool fulls_resident, bool full_stats) { // A full is usable for the merge / error model if it passes AddImage's filters (with the current // ice context). group >= 0 already encodes "not absent and passes AcceptReflection". auto usable_merge = [&](const Obs &o) { @@ -3020,9 +2992,11 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool em.a, em.b, em.isa, result.isa_asymptotic, 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. + // R-free flags and French-Wilson amplitudes are fields of the WRITTEN reflection, so a merge that + // is not written needs neither. French-Wilson is deferred until after the anomalous accumulator + // below has attached I(+)/I(-), so the two hands get their amplitudes in one pass. + if (full_stats) + AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction); if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); @@ -3069,71 +3043,77 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } } - // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate - // |I_i - | per reflection. - struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; - std::vector rmeas(n_groups); - bool did_gpu_rmeas = false; -#ifdef JFJOCH_USE_CUDA - if (use_gpu_merge) { - // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). - std::vector rabsdev(n_groups), rsumI(n_groups); - std::vector rn(n_groups), rnusable(n_groups); - gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); - for (int g = 0; g < n_groups; ++g) { - if (rnusable[g] == 0) continue; - const auto shell = shells.GetShell(acc[g].d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately - // ignores the ice flag on a search pass, so its count includes observations that never - // entered `unique` - which inflates the reported multiplicity of whatever shell they - // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a - // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above - // already declines it; this is the same statement made where it counts.) - sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); - if (std::isfinite(merged_I[g]) && rn[g] > 0) { - auto &r = rmeas[g]; - r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; - } - } - did_gpu_rmeas = true; - } -#endif - if (!did_gpu_rmeas) - for (const auto &o : fulls) { - if (o.group < 0) continue; - if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas - if (!frame_cell_ok[o.frame]) continue; - if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; - if (o.partiality < min_partiality) continue; - const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; - if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; - const auto shell = shells.GetShell(o.d); - if (!shell || *shell < 0 || *shell >= n_shells) continue; - // Only what the merge kept counts towards multiplicity - see the GPU branch above. The - // R_meas accumulation below keeps its own, wider filter. - if (usable_merge(o)) - sa[*shell].total_obs++; - if (std::isfinite(merged_I[o.group])) { - auto &r = rmeas[o.group]; - r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); - r.sum_I += I_corr; r.n++; r.shell = *shell; - } - } std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); double rmeas_num_all = 0.0, rmeas_den_all = 0.0; - for (const auto &r : rmeas) { - if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; - const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); - rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; - rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; + // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last + // full pass over the observations there is - so, like the correction surfaces, it is not made on a + // merge whose numbers nobody reads. + if (full_stats) { + // R_meas: re-walk the fulls (Mask = cell only; no ice / error-model), accumulate + // |I_i - | per reflection. + struct RmeasObs { double sum_abs_dev = 0, sum_I = 0; int n = 0, shell = -1; }; + std::vector rmeas(n_groups); + bool did_gpu_rmeas = false; +#ifdef JFJOCH_USE_CUDA + if (use_gpu_merge) { + // Per-group R_meas + usable count on the GPU; the shell is assigned per group (its fulls share d). + std::vector rabsdev(n_groups), rsumI(n_groups); + std::vector rn(n_groups), rnusable(n_groups); + gpu_->MergeRmeas(merged_I.data(), rabsdev.data(), rsumI.data(), rn.data(), rnusable.data()); + for (int g = 0; g < n_groups; ++g) { + if (rnusable[g] == 0) continue; + const auto shell = shells.GetShell(acc[g].d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Count the MERGED population, not the R_meas one. The R_meas re-walk deliberately + // ignores the ice flag on a search pass, so its count includes observations that never + // entered `unique` - which inflates the reported multiplicity of whatever shell they + // land in. acc[g].nh is what actually went into this group's mean, and it is zero for a + // group the merge dropped entirely. (acc[g].d is NaN for such a group, so GetShell above + // already declines it; this is the same statement made where it counts.) + sa[*shell].total_obs += static_cast(acc[g].nh[0] + acc[g].nh[1]); + if (std::isfinite(merged_I[g]) && rn[g] > 0) { + auto &r = rmeas[g]; + r.sum_abs_dev = rabsdev[g]; r.sum_I = rsumI[g]; r.n = rn[g]; r.shell = *shell; + } + } + did_gpu_rmeas = true; + } +#endif + if (!did_gpu_rmeas) + for (const auto &o : fulls) { + if (o.group < 0) continue; + if (rejected_obs[&o - fulls.data()]) continue; // outlier-rejected in the merge -> also out of R_meas + if (!frame_cell_ok[o.frame]) continue; + if (!(o.corr > 0.0f) || !std::isfinite(o.corr)) continue; + if (o.partiality < min_partiality) continue; + const float I_corr = o.I * o.corr, sigma_corr = o.sigma * o.corr; + if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) continue; + const auto shell = shells.GetShell(o.d); + if (!shell || *shell < 0 || *shell >= n_shells) continue; + // Only what the merge kept counts towards multiplicity - see the GPU branch above. The + // R_meas accumulation below keeps its own, wider filter. + if (usable_merge(o)) + sa[*shell].total_obs++; + if (std::isfinite(merged_I[o.group])) { + auto &r = rmeas[o.group]; + r.sum_abs_dev += std::fabs(static_cast(I_corr) - merged_I[o.group]); + r.sum_I += I_corr; r.n++; r.shell = *shell; + } + } + for (const auto &r : rmeas) { + if (r.n < 2 || r.shell < 0 || r.shell >= n_shells) continue; + const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); + rmeas_num[r.shell] += factor * r.sum_abs_dev; rmeas_den[r.shell] += r.sum_I; + rmeas_num_all += factor * r.sum_abs_dev; rmeas_den_all += r.sum_I; + } } // ---- Anomalous split (always, even when the merge is Friedel-averaged): for each acentric // reflection, the inverse-variance I(+)/I(-) from the SAME scaled fulls (so it never touches the // Friedel-merged IMEAN, scaling or error model above). Lets I(+)/I(-) and F(+)/F(-) be written by // default without scaling anomalously; a reflection with only one mate, or a centric, is left - // without the split. Skipped for the P1 search pass, which has no use for it. ---- + // without the split. Skipped for the P1 search pass, which has no use for it, and for a merge that + // is not written. ---- struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; }; std::unordered_map anom_export; // SigAno = <|I(+)-I(-)|>/ accumulated on the standard report shells (num=sum|dI|, @@ -3141,7 +3121,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool std::vector sig_num(n_shells, 0.0), sig_den(n_shells, 0.0); double sig_num_all = 0.0, sig_den_all = 0.0; size_t sig_n = 0; - if (!for_search) { + if (!for_search && full_stats) { const int sg_num = x.GetSpaceGroupNumber().value_or(1); const HKLKeyGenerator anom_keygen(/*merge_friedel=*/false, sg_num); const gemmi::GroupOps gops = gemmi::find_spacegroup_by_number(sg_num)->operations(); @@ -3233,15 +3213,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } // French-Wilson amplitudes for IMEAN and (now that they are attached) each Bijvoet hand. - FrenchWilsonOptions fw_opts; - fw_opts.num_threads = static_cast(nthreads); - ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1), fw_opts); + if (full_stats) { + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = static_cast(nthreads); + ApplyFrenchWilson(result.merged, x.GetSpaceGroupNumber().value_or(1), fw_opts); + } logger.Info("Merge complete ({} unique reflections)", result.merged.size()); return result; } -RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { +RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_stats) { const int sg_number = x.GetSpaceGroupNumber().value_or(1); HKLKeyGenerator keygen(merge_friedel, sg_number); @@ -3545,26 +3527,29 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { // diffracted-beam direction), each an alternating multiplicative fit of the fulls' corr against // the merged reference. Cheap host loops over the downloaded fulls; applied only on the final // in-symmetry merge, never the P1 space-group search pass (corrections there add risk and can - // perturb the symmetry determination). --- + // perturb the symmetry determination) and never on a merge that is not an output: a surface + // is a correction to the merged INTENSITIES, and a pass whose intensities are thrown away has + // nowhere to put one - the fitted gain lives on the fulls, which the next pass rebuilds. --- + const bool refine_surfaces = !for_search && full_stats; [[maybe_unused]] const bool corrections = - !for_search && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); + refine_surfaces && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); // Radiation-damage monitor: measure the per-batch relative-B on the scaled fulls BEFORE any correction // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass. - if (!for_search) + if (refine_surfaces) MeasureRadiationDamageB(n_groups); // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is // already out of it) and the per-frame CC computed above. Report-only; drops nothing. - if (!for_search) + if (refine_surfaces) MeasureSweepQuality(partial_scaled, cc, cc_n); - if (!for_search && refine_decay_b) + if (refine_surfaces && refine_decay_b) RefineDecay(n_groups); - if (!for_search && relative_b_deg > 0.0) + if (refine_surfaces && relative_b_deg > 0.0) RefineRelativeB(n_groups); // per-batch relative-B on top of the single decay slope - if (!for_search && absorption_iter > 0) + if (refine_surfaces && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); - if (!for_search && modulation_iter > 0) + if (refine_surfaces && modulation_iter > 0) RefineModulation(modulation_iter, n_groups); - if (!for_search && absorption_iter > 0) + if (refine_surfaces && absorption_iter > 0) RefineAbsorptionTime(absorption_iter, n_groups); // last: the static surfaces get first claim #ifdef JFJOCH_USE_CUDA // The corrections mutate the host fulls' corr; when the merge runs on the resident (GPU) fulls, push @@ -3580,6 +3565,6 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) { #endif // --- 5. Error model + merge + statistics. --- - auto r = MergeAndStats(n_groups, for_search, combined_on_gpu && scaled_fulls_on_gpu); + auto r = MergeAndStats(n_groups, for_search, combined_on_gpu && scaled_fulls_on_gpu, full_stats); return r; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 5255af16..35fadbbb 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -71,7 +71,14 @@ public: // for the space group currently set on the experiment, reusing the ingested buffers. // for_search: the de-novo P1 pass whose merged intensities feed the space-group search - ice-ring // reflections are dropped from the merge and the error model (kept otherwise, for completeness). - Result Run(bool for_search); + // full_stats: these merged intensities are an OUTPUT. False on the rotation two-pass geometry + // pre-pass, whose merge exists only to choose the space group and post-refine the geometry and + // whose reflections are never written: the correction surfaces, the report-only diagnostics, the + // R_meas re-walk, the anomalous split, the R-free flags and the French-Wilson amplitudes are then + // all skipped, because computing them fills in fields nothing reads. What the pre-pass IS read for + // - the merged intensities themselves, the error model, and the completeness / CC1/2 the second + // pass is judged against - is computed either way. + Result Run(bool for_search, bool full_stats); // Override the high-resolution cut for the next Run() - used to gate the de-novo P1 search pass at // >= 1 without cutting the final in-symmetry merge. Reset to the manual limit afterwards. @@ -99,6 +106,17 @@ private: int32_t group; // dense ASU-group id for the current space group; <0 = never mergeable }; + // 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 + // total - see the .cpp), and the resolution the range test reads. Twenty-four bytes against the + // Obs's eighty, and it is all the ingest needs before it knows which observations survive. + struct SortKey { + int32_t h, k, l; + float image_number; + int32_t idx; + float d; + }; + const DiffractionExperiment &x; std::vector &partials_out; // written back at the end of scaling std::optional reference_cell; @@ -224,14 +242,17 @@ private: // background is a usable flux meter and what it costs when it is not. void MeasureIncidentFlux(const std::vector &mean_bkg); - // Drop the observations whose resolution can never be in range: --scaling-high/low-resolution are the - // coarsest limits any Run() uses (the space-group search only ever RAISES d_min), and an out-of-range - // raw hkl gets group -1 in every pass, which keeps it out of the scaling reference, the per-frame fit, - // the combine, the merge and the error model alike. On a crystal that integrates to the detector - // corner and merges well short of it that is most of the array. Whole raw-hkl RUNS are dropped, on the - // same per-hkl resolution ComputeAsuGroups tests, so what survives - and the order of every sum formed - // over it - is exactly what it was. No-op when no manual limit was given. - void DropOutOfRangeObservations(); + // Build the flat `partials` array (and the per-frame CSR, the finiteness mask and `perm`) from the + // source reflections, skipping the observations whose resolution can never be in range: + // --scaling-high/low-resolution are the coarsest limits any Run() uses (the space-group search only + // ever RAISES d_min), and an out-of-range raw hkl gets group -1 in every pass, which keeps it out of + // the scaling reference, the per-frame fit, the combine, the merge and the error model alike. On a + // crystal that integrates to the detector corner and merges well short of it that is most of the + // array, and an eighty-byte record built for it is eighty bytes written and then thrown away. Whole + // raw-hkl RUNS are skipped, on the same per-hkl resolution ComputeAsuGroups tests, so what survives - + // and the order of every sum formed over it - is exactly what it would have been had the whole array + // been built and then filtered. Everything is built when no manual limit was given. + void BuildInRangeObservations(const std::vector &keys); // Compute the dense ASU-group id for the current space group by grouping the (pre-sorted) raw-hkl // runs by their ASU key - one gemmi ASU reduction per distinct raw hkl, not per observation. Fills @@ -347,5 +368,6 @@ private: // Error model + merge + statistics over the fulls (the last stage). n_groups is the fulls group count. // fulls_resident: the (scaled) fulls + their group CSR are still on the GPU, so the em-stats / samples // / merge-accumulate / R_meas reductions run there (only per-group + samples come back). - Result MergeAndStats(int n_groups, bool for_search, bool fulls_resident); + // full_stats: see Run(). + Result MergeAndStats(int n_groups, bool for_search, bool fulls_resident, bool full_stats); }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6255c59a..a8465297 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2090,7 +2090,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult { if (rsm) { phase("Scale/combine/merge (" + label + ")"); - auto r = rsm->Run(for_search); + // The geometry pre-pass merges only to choose the space group and to give the quality + // guard something to judge the second pass against; its reflections are never written, + // so the parts of the merge that only fill in an output file are skipped there. + auto r = rsm->Run(for_search, /*full_stats=*/!geometry_prepass); result.error_model_isa = r.isa; result.error_model_isa_asymptotic = r.isa_asymptotic; result.error_model_a = r.error_model_a; @@ -2564,11 +2567,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto twin_sg_number = experiment_.GetSpaceGroupNumber(); const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; - result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + // Not on the geometry pre-pass: the analysis goes into that pass's statistics text and its + // written reflections, and neither survives the run. The promotion flag below is a different + // thing - it is what the SEARCH did, the second pass reads it, and it is set either way. + if (!geometry_prepass) { + result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + stats_text << TwinningAnalysisToText(result.twinning) << "\n"; + } // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the // search itself made, so the text cannot claim "no twin law exists" on its own say-so. result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; - stats_text << TwinningAnalysisToText(result.twinning) << "\n"; // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the // crystal was mounted, which the user can change), warned about when the two nearly coincide. @@ -2639,8 +2647,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not - // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log. - { + // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log, none + // of which the geometry pre-pass produces. + if (!geometry_prepass) { const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); sm.statistics.wilson_b = wilson.b; sm.statistics.wilson_b_correlation = wilson.correlation; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index fedc5e09..fae3e0a6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1424,7 +1424,7 @@ static int RunRugnux(int argc, char **argv) { RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), scaling_iter, nthreads, logger); rsm.Ingest(); - auto r = rsm.Run(false); + auto r = rsm.Run(false, /*full_stats=*/true); merged_reflections = std::move(r.merged); merged_statistics = std::move(r.statistics); error_model_isa = r.isa;