diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 43a0921fb..6c2bce937 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/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(std::clamp(nthreads, 1, std::max(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 obs(new const Partial *[n_obs]); - ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { - size_t at = obs_offset[static_cast(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 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(n_obs); - ParallelChunks(static_cast(n_obs), nthreads, [&](int lo, int hi) { + ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { + size_t at = obs_offset[static_cast(lo) / obs_chunk]; for (int i = lo; i < hi; ++i) - ranked[i] = {obs[i]->I / std::max(1e-9, static_cast(obs[i]->sigma)), obs[i]}; + if (keep_obs(i)) + ranked[at++] = {pts[i].I / std::max(1e-9, static_cast(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(n_pts), nthreads, [&](int lo, int hi) { + size_t at = obs_offset[static_cast(lo) / obs_chunk]; + for (int i = lo; i < hi; ++i) if (keep_obs(i)) kept[at++] = pts[i]; + }); } + observations.pts.reset(); + std::unique_ptr obs(new const Partial *[n_obs]); + for (size_t i = 0; i < n_obs; ++i) + obs[i] = &kept[i]; result.obs_used = static_cast(n_obs); // ===== The joint fit ===== diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 5c313d8d2..b1aa6d658 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/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 diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index f072be960..658fadb27 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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(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