RotationScaleMerge: MergeAndStats' per-pass arrays as members, reset on all threads

The rejected flags and the observation keys and halves (per full), the merged
intensities (per group) and the Wilson observation records were locals: on a
large cell every one of the run's merges allocated, page-faulted and zeroed
them again on one thread - hundreds of megabytes per merge on 8qaw. They are
now members, only resized, and reset with ParallelChunks where a pass needs
them reset. The error-model samples are no longer cleared before the device
download overwrites them (the host path still clears before it appends).

p.mtz md5 unchanged on myob, cytc, thau, kdp, 8a1a and 8qaw (GPU build).

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -4696,7 +4696,6 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     bool error_model_b_unmeasured = false;   // b had no leverage; ISa is not reported
     bool error_model_b_resolved = true;      // b^2 stood clear of its own standard error
     std::vector<Sample> &samples = em_samples;
-    samples.clear();
     // 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 <I>. Rebuild the variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*<I>,
@@ -4780,6 +4779,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
             for (int g = 0; g < n_groups; ++g)
                 if (cnt[g] >= 2 && sw[g] > 0.0) em_mean[g] = swI[g] / sw[g];
 
+            samples.clear();
             samples.reserve(fulls.size());
             for (const auto &o : fulls) {
                 if (!usable_merge(o) || cnt[o.group] < 2) continue;
@@ -5079,8 +5079,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     // and whether it is usable at all are fixed here, and neither depends on the error model that the
     // merge below refits. Rank parity within the reflection - see AssignHalvesOfMembers - so the halves are
     // balanced and identical on the host and device paths.
-    std::vector<uint64_t> obs_key(fulls.size());
-    std::vector<uint8_t> merge_half(fulls.size(), 0);
+    std::vector<uint64_t> &obs_key = merge_obs_key;
+    std::vector<uint8_t> &merge_half = merge_half_set;
+    obs_key.resize(fulls.size());
+    merge_half.resize(fulls.size());
     bool did_gpu_halves = false;
 #ifdef JFJOCH_USE_CUDA
     // On the device where the merge is: integer work, the same halves.
@@ -5092,8 +5094,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     if (!did_gpu_halves) {
         ParallelChunks(static_cast<int>(fulls.size()), ThreadsForWork(fulls.size(), nthreads),
                        [&](int lo, int hi) {
-            for (int i = lo; i < hi; ++i)
+            for (int i = lo; i < hi; ++i) {
                 obs_key[i] = ObservationKey(fulls[i].h, fulls[i].k, fulls[i].l, fulls[i].frame);
+                merge_half[i] = 0;
+            }
         });
         // mf.group is already this array: it was filled with the same `usable_merge(o) ? o.group : -1`
         // above, on all threads, so a second serial walk over the fat records only rebuilds it.
@@ -5139,9 +5143,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     //      model has been refitted on the reflections that survive the cut. ----
     std::vector<Accum> &acc = merge_acc;
     Result result;
-    std::vector<double> merged_I;
+    std::vector<double> &merged_I = merge_merged_I;
     size_t reject_count = 0;
-    std::vector<uint8_t> rejected_obs;   // per-full outlier-rejected flag (both paths)
+    std::vector<uint8_t> &rejected_obs = merge_rejected;   // per-full outlier-rejected flag (both paths)
 
     // ---- The Wilson test (WilsonOutliers.h): every observation against the intensity distribution of
     //      its shell, beside the median test. It reaches what the median cannot - a reflection measured
@@ -5180,15 +5184,22 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
         // 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);
+        merged_I.resize(n_groups);
+        rejected_obs.resize(fulls.size());
+        ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int lo, int hi) {
+            std::fill(merged_I.begin() + lo, merged_I.begin() + hi, NAN);
+        });
+        ParallelChunks(n_full, ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) {
+            std::fill(rejected_obs.begin() + lo, rejected_obs.begin() + hi, uint8_t{0});
+        });
         reject_count = 0;
-        rejected_obs.assign(fulls.size(), 0);
         if (wilson_test) {
             // Each observation with its OWN sigma, the error model evaluated at its own intensity: the
             // reflection's mean, which the merge weights use, is exactly what an artefact inflates.
             // Which fulls take part, in order, and then each one's record on its own - on all threads.
             wilson_full = KeptIndices(n_full, nthreads, [&](int i) { return mf.group[i] >= 0; });
-            std::vector<WilsonObservation> wobs(wilson_full.size());
+            std::vector<WilsonObservation> &wobs = merge_wobs;
+            wobs.resize(wilson_full.size());
             ParallelChunks(static_cast<int>(wobs.size()), ThreadsForWork(wobs.size(), nthreads),
                            [&](int lo, int hi) {
                 for (int j = lo; j < hi; ++j) {

image_analysis/scale_merge/RotationScaleMerge.h

@@ -17,6 +17,7 @@
 #include "../IntegrationOutcome.h"
 
 #include "ErrorModel.h"
+#include "WilsonOutliers.h"
 #include "Merge.h"        // MergedReflection, MergeStatistics
 #include "SphericalHarmonicSurface.h"
 #ifdef JFJOCH_USE_CUDA
@@ -427,6 +428,15 @@ private:
     };
     MergeAccumStaging merge_accum;
 
+    // MergeAndStats' other per-full and per-group arrays, members for the same reason: five merges a
+    // run, each over up to tens of millions of fulls and millions of groups, and as locals each one
+    // allocated, faulted in and zeroed them again on one thread. As members they are only resized, and
+    // reset on all threads where a pass needs them reset.
+    std::vector<double> merge_merged_I;
+    std::vector<uint8_t> merge_rejected, merge_half_set;
+    std::vector<uint64_t> merge_obs_key;
+    std::vector<WilsonObservation> merge_wobs;
+
     // 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; };