Merge host-mem into rc175: the merge's working arrays handed back before the space-group search and after the all-observation merge, post-refinement frees its gather before the fit; output unchanged

Peak live heap: 4nwv (GPU) 17.1 -> 16.2 GiB, 8qaw 21.1 -> 20.0 GiB.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4

image_analysis/geom_refinement/PostRefine.cpp

@@ -473,10 +473,9 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
 
             // The observed spot positions, one entry per partial that produced a centroid. Count first,
             // then fill, exactly as the partial gather above does and for the same reason: this walks the
-            // same tens of millions of partials, and a pointer vector grown by push_back copies itself
-            // every time it doubles. new[] rather than a sized vector so the array is not zeroed on one
-            // thread before the parallel fill overwrites it. The fill lands in the order the serial loop
-            // produced, so the selection below sees the same sequence it always did.
+            // same tens of millions of partials, and a vector grown by push_back copies itself every
+            // time it doubles. The fill lands in the order the serial loop produced, so the selection
+            // below sees the same sequence it always did.
             const int n_obs_chunks = static_cast<int>(std::clamp<size_t>(nthreads, 1,
                                                                          std::max<size_t>(1, n_pts)));
             const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks;
@@ -491,28 +490,39 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
             });
             for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c];
             size_t n_obs = obs_offset[n_obs_chunks];
-            std::unique_ptr<const Partial *[]> obs(new const Partial *[n_obs]);
-            ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
-                size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
-                for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i];
-            });
+            // The fit reads at most MAX_OBS of them, the strongest by I/sigma. Those are copied out, and
+            // the gathered partials - gigabytes on a large cell, and read by nothing below - go back now
+            // rather than when the fit returns.
             constexpr size_t MAX_OBS = 20000;
+            std::vector<Partial> kept;
             if (n_obs > MAX_OBS) {
                 // By I/sigma, with the ratio worked out once per observation beside it rather than at
-                // every comparison through the pointer: the same comparisons, so the same selection in
-                // the same order.
+                // every comparison: the same comparisons, so the same selection in the same order.
                 struct Ranked { double i_over_sigma; const Partial *p; };
                 std::vector<Ranked> ranked(n_obs);
-                ParallelChunks(static_cast<int>(n_obs), nthreads, [&](int lo, int hi) {
+                ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
+                    size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
                     for (int i = lo; i < hi; ++i)
-                        ranked[i] = {obs[i]->I / std::max(1e-9, static_cast<double>(obs[i]->sigma)), obs[i]};
+                        if (keep_obs(i))
+                            ranked[at++] = {pts[i].I / std::max(1e-9, static_cast<double>(pts[i].sigma)), &pts[i]};
                 });
                 std::nth_element(ranked.begin(), ranked.begin() + MAX_OBS, ranked.end(),
                                  [](const Ranked &a, const Ranked &b) { return a.i_over_sigma > b.i_over_sigma; });
+                kept.resize(MAX_OBS);
                 for (size_t i = 0; i < MAX_OBS; ++i)
-                    obs[i] = ranked[i].p;
+                    kept[i] = *ranked[i].p;
                 n_obs = MAX_OBS;
+            } else {
+                kept.resize(n_obs);
+                ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
+                    size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
+                    for (int i = lo; i < hi; ++i) if (keep_obs(i)) kept[at++] = pts[i];
+                });
             }
+            observations.pts.reset();
+            std::unique_ptr<const Partial *[]> obs(new const Partial *[n_obs]);
+            for (size_t i = 0; i < n_obs; ++i)
+                obs[i] = &kept[i];
             result.obs_used = static_cast<int>(n_obs);
 
             // ===== The joint fit =====

image_analysis/scale_merge/RotationScaleMerge.h

@@ -183,6 +183,52 @@ public:
     // beside others with it off, at the point where it would have run in series.
     void WriteBackPerFrameScale() { FinalizePerFrameScale(last_cc, last_cc_n, last_frame_scaled); }
 
+    // Hand back the merge's working arrays - MergeAndStats' accumulators, R_meas sums, ISa scatter,
+    // error-model pools and per-full fields. Each is resized and refilled before it is read, so the
+    // next merge starts from them as an engine's first merge does, and gives the same answer. On a
+    // large cell in P1 they are gigabytes, kept between merges only to save faulting them in again;
+    // this is for a gap between merges where something else sits on the run's memory peak.
+    void ReleaseMergeScratch() {
+        merge_fields = {};
+        merge_acc = {};
+        merge_accum = {};
+        merge_merged_I = {};
+        merge_rejected = {};
+        merge_half_set = {};
+        merge_obs_key = {};
+        merge_wobs = {};
+        merge_rmeas = {};
+        merge_rmeas_staging = {};
+        asymptote_scatter = {};
+        em_samples = {};
+        em_cut_pool = {};
+        em_fit_pool = {};
+        em_chi2 = {};
+    }
+
+    // For an engine that will not merge again: everything on the host but the per-frame state
+    // WriteBackPerFrameScale reads (partials_out, g_partial, frame_flux, mos_smooth, last_cc,
+    // last_cc_n, last_frame_scaled). Its card memory goes with the engine.
+    void ReleaseAllButPerFrameScale() {
+        ReleaseMergeScratch();
+        partials = {};
+        fulls = {};
+        fulls_staging = {};
+        fulls_group_start = {};
+        fulls_group_perm = {};
+        perm = {};
+        rawrun_start = {};
+        rawrun_count = {};
+        rawrun_h = {};
+        rawrun_k = {};
+        rawrun_l = {};
+        rawrun_d = {};
+        rawrun_group = {};
+        finite_ok = {};
+        corr_ingested = {};
+        obs_position.reset();
+    }
+
     // A second engine on this one's ingest, on `device` and logging to `log`, for a merge to run beside
     // this engine's: the host half of the ingest is copied, and the copy's IngestDevice(this) is still
     // to be called. Call after this engine's IngestDevice(). No merge depends on the merges run

rugnux/RugnuxScaleMerge.cpp

@@ -1183,8 +1183,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
             all_arm_ahead = std::async(std::launch::async, [&rsm_all, &all_arm_log] {
                 pin_gpu(rsm_all->GetGpuDevice());   // the card and its NUMA node, as the guard's merge has them
                 TimingMark mark(all_arm_log, "merge P1, all observations");
-                return rsm_all->Run(/*for_search=*/true, /*full_stats=*/true,
-                                    /*measure_cc_before_corrections=*/false);
+                auto r = rsm_all->Run(/*for_search=*/true, /*full_stats=*/true,
+                                      /*measure_cc_before_corrections=*/false);
+                // This engine merges once: what is left for it to do is the write-back of its
+                // per-frame scales, so its observations and its merge's working arrays go now,
+                // rather than through the search the others are about to run.
+                rsm_all->ReleaseAllButPerFrameScale();
+                return r;
             });
         }
     }
@@ -1593,6 +1598,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
         postrefine_gathered_future.wait();   // the search re-indexes the reflections the gather reads
     TimingMark search_mark(logger, "space-group search");
     if (search_space_group) {
+        // The engine's next merge is in the group the search picks; until then its last merge's
+        // working arrays only sit beside the searches' own (gigabytes on a large cell in P1).
+        if (rsm)
+            rsm->ReleaseMergeScratch();
         SearchSpaceGroupOptions sg_opts;
         sg_opts.nthreads = static_cast<size_t>(std::max(1, config_.nthreads));
         sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel();
@@ -1690,6 +1699,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
                 rsm->SetSearchMinZeta(0.0);
                 sm_all = scale_and_merge("P1, all observations", true);
                 rsm->SetSearchMinZeta(zeta);
+                rsm->ReleaseMergeScratch();   // as before the first arm's search
                 sg_opts.merge_isa = result.error_model_isa;   // this arm's own error model
                 TimingMark all_arm_mark(logger, "SearchSpaceGroup (all-observation arm)");
                 alt = SearchSpaceGroup(sm_all.merged, sg_opts);
@@ -3300,6 +3310,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
     }
     adoption_mark.Stop();
     TimingMark analyses_mark(logger, "analyses");
+    // The answer's merge is made and in `sm`; a merge below (twin law, ledger) refills the arrays.
+    if (rsm)
+        rsm->ReleaseMergeScratch();
 
     // --finalist-ledger (see where it is read, below) is folded from a merge of its own: the same
     // observations in P1. It needs every node of the lattice integrated, which a fixed centred group