From bbb341a566b8ad6fbbe7fa80d7815d960993e96b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 18 Oct 2026 10:34:41 +0200 Subject: [PATCH] RotationScaleMerge: TIMING marks inside the surfaces and the merge, decay terms built once, merge scratch not zeroed on one thread Prototype for the tail campaign (output unchanged, checked on 4nwv --mode scale --no-gpu: same .hkl, same .mtz data): - rsm/s and rsm/m TIMING sub-marks for each correction surface, the disposition and each stage of MergeAndStats (the pieces the tail is spent in had no marks of their own). - RefineDecay takes the terms MeasureRadiationDamageB just built (nothing between them changes corr) instead of building them again, and evaluates its decay factor once per term where it used it twice. - The merge's per-group and per-full scratch (MergeFields, FullsStaging, the accumulators, the R_meas sums and their staging, the ISa scatter) default-initialised: each is written in full before it is read, and std::vector zeroed - and faulted in - gigabytes of it on one thread per P1 merge. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4 --- .../scale_merge/RotationScaleMerge.cpp | 69 +++++++++++++++---- .../scale_merge/RotationScaleMerge.h | 45 +++++++----- 2 files changed, 83 insertions(+), 31 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index a3caf75a3..fa02a8d48 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -2125,8 +2125,14 @@ void RotationScaleMerge::RefineDecay(int n_groups) { // fat records once, with the group-ordered copy that lets the per-group references run on all // threads. See BuildDecayTerms: the compacted array is walked in fulls order and a group's terms are // still summed in fulls order, so every sum here is formed exactly as the serial walks formed it. + // The monitor's terms where it has just built them (decay_terms_kept). DecayTerms t; - BuildDecayTerms(n_groups, t); + if (decay_terms_kept_valid) + t = std::move(decay_terms_kept); + else + BuildDecayTerms(n_groups, t); + decay_terms_kept = {}; + decay_terms_kept_valid = false; const int nt = static_cast(ThreadsForWork(t.term.size(), nthreads)); std::vector sw(n_groups), swI(n_groups); @@ -2177,8 +2183,9 @@ void RotationScaleMerge::RefineDecay(int n_groups) { for (int k = t.g_start[g]; k < t.g_start[g + 1]; ++k) { const DecayTerm &o = t.gterm[k]; if ((o.frame & 1) != parity) continue; - const double Is = static_cast(o.I) * o.corr * decay_factor(o, slope); - const double sc = static_cast(o.sigma) * o.corr * decay_factor(o, slope); + const double f = decay_factor(o, slope); + const double Is = static_cast(o.I) * o.corr * f; + const double sc = static_cast(o.sigma) * o.corr * f; const double w = 1.0 / (sc * sc); s_w += w; s_wI += w * Is; } @@ -2354,9 +2361,10 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { auto batch_of = [&](const DecayTerm &o) { return std::min(n_batch - 1, o.frame / frames_per_batch); }; // The fields the three walks below read, pulled out of the fat records once, with the group-ordered - // copy the low-dose reference runs over (BuildDecayTerms). - DecayTerms t; + // copy the low-dose reference runs over (BuildDecayTerms), kept for RefineDecay. + DecayTerms &t = decay_terms_kept; BuildDecayTerms(n_groups, t); + decay_terms_kept_valid = true; const int nt = static_cast(ThreadsForWork(t.term.size(), nthreads)); // Low-dose reference: inverse-variance group means over the first ~10% of frames (>= 1 batch). @@ -4819,6 +4827,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // the device's error-model samples - read four or five fields of an eighty-byte record, so each of // them dragged the whole fulls array through the cache to touch twenty bytes an observation. Pull // those fields out once, on all threads, and let the walks run over flat arrays instead. + std::unique_ptr sub_mark; + auto sub = [&](const char *n) { sub_mark.reset(); sub_mark = std::make_unique(*logger, std::string("rsm/m ") + n); }; + sub("fields+hands"); const int n_full = static_cast(fulls.size()); MergeFields &mf = merge_fields; mf.Resize(n_full); @@ -4916,6 +4927,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool }; // Friedel pair of each group when the hands are merged separately (-A); filled with the medians below. std::vector pair_of_group; + sub("em samples"); { std::vector cnt(n_groups, 0); // per-group usable count (both paths; feeds reject-median) bool did_gpu = false; @@ -4991,6 +5003,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Per-group outlier-rejection median of I*corr (host both paths - a per-group median is awkward on // the GPU; cheap here, cnt >= 3 filter from the em pass). Fed to the merge accumulate. + sub("reject medians"); if (reject_outliers) { // One flat array with a per-group span, not a vector per group: n_groups is over a million // on a P1 pass, so a vector each is a million allocations, their growth, and a million frees @@ -5155,7 +5168,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool reject_log_sd = 1.4826 * lv[j].first; } } + sub("em fit"); fit_error_model(samples); + sub("bijvoet excess"); } // The full's sigma under the error model, with the variance evaluated at intensity I_for_b. @@ -5270,6 +5285,7 @@ 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. + sub("halves+factors"); std::vector &obs_key = merge_obs_key; std::vector &merge_half = merge_half_set; obs_key.resize(fulls.size()); @@ -5332,7 +5348,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // ---- Merge: per-group inverse-variance sums with corrected sigma + deterministic half sets, then // export. In a lambda because the automatic resolution cutoff below re-runs it once the error // model has been refitted on the reflections that survive the cut. ---- - std::vector &acc = merge_acc; + FillVector &acc = merge_acc; Result result; std::vector &merged_I = merge_merged_I; size_t reject_count = 0; @@ -5369,7 +5385,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // empty ones too - so there this only has to be the right size. The host path accumulates into // it with += and does need the zeroes. if (use_gpu_merge) acc.resize(n_groups); - else acc.assign(n_groups, Accum{}); + else acc.assign(n_groups, ACCUM_EMPTY); merged_I.resize(n_groups); rejected_obs.resize(fulls.size()); ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int lo, int hi) { @@ -5414,6 +5430,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const double mult_sd = std::isfinite(reject_log_sd) ? reject_log_sd : error_model_b; const double reject_up = std::expm1(reject_nsigma * mult_sd); const double reject_down = -std::expm1(-reject_nsigma * mult_sd); + sub("mg accum"); bool did_gpu_acc = false; #ifdef JFJOCH_USE_CUDA if (use_gpu_merge) { @@ -5504,6 +5521,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool }, 512); for (size_t r : blk_rejected) reject_count += r; } + sub("mg d+export"); // A merged reflection's d from the cell the completeness denominator is enumerated with, not // from whichever frame first observed it: a per-frame cell puts a reflection near a shell edge // on the other side of it from its possible twin, and a shell then reads above 100 %. @@ -5555,7 +5573,9 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } }); }; + sub("merge 1"); run_merge(); + sub("cutoff+refit+merge 2"); // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported shells to // the CC1/2 fall-off. The scaling and combine above ran over the full range, and the per-image @@ -5622,16 +5642,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Not measured on a merge whose numbers nothing reads (the geometry pre-pass, the discarded // pre-correction CC1/2 merge): it is a report-only diagnostic, and the group-scatter array alone is // a few hundred megabytes value-initialised on one thread at a million groups. + sub("asymptote+chi2+rfree"); double error_model_b_asymptotic = error_model_b; if (full_stats && error_model_active) { - std::vector &gs = asymptote_scatter; - gs.assign(n_groups, GroupScatter{}); + FillVector &gs = asymptote_scatter; + gs.resize(n_groups); bool did_gpu_scatter = false; #ifdef JFJOCH_USE_CUDA // On the device where the merge is, to the same bits. if (use_gpu_merge) { - std::vector sum(n_groups), sum_sq(n_groups), sum_var(n_groups); - std::vector n(n_groups); + FillVector sum(n_groups), sum_sq(n_groups), sum_var(n_groups); + FillVector n(n_groups); gpu_->MergeGroupScatter(sum.data(), sum_sq.data(), sum_var.data(), n.data()); ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int lo, int hi) { for (int g = lo; g < hi; ++g) gs[g] = {sum[g], sum_sq[g], sum_var[g], n[g]}; @@ -5643,6 +5664,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool if (!did_gpu_scatter) { ParallelBlocks(n_groups, nthreads, [&](int, int g0, int g1) { for (int gi = g0; gi < g1; ++gi) { + gs[gi] = {0.0, 0.0, 0.0, 0}; for (int p = fulls_group_start[gi]; p < fulls_group_start[gi + 1]; ++p) { const int i = fulls_group_perm[p]; if (mf.group[i] != gi) continue; @@ -5672,6 +5694,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // median-to-mean conversion the (a, b) fit does with CHI2_1_MEDIAN. // The clamp at zero is gone with them: median(max(x,0)) = max(median(x),0), so it never moved a // positive median, and a negative median is the informative answer "not measurable here". + sub("asy groups"); std::vector> group_scatter; // (systematic b^2, I/sigma) for (const auto &g : gs) { if (g.n < 5) continue; @@ -5724,6 +5747,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // here, once, instead of in every fit: the pool is one sample per full and the median is a pass over // it plus an nth_element, ten times a run. double error_model_chi2 = 0.0; + sub("chi2+rfree"); if (full_stats && chi2_pool) { std::vector &chi2 = em_chi2; chi2.clear(); @@ -5787,6 +5811,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // it bottoms 0.1% below the finest kept reflection, and the groups in that band were erased // from the merge, so without the floor they would count in N_obs and R_meas but not in // N_uniq. ---- + sub("statistics"); const int n_shells = report_shell_count; struct AnomExport { float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; }; std::unordered_map anom_export; // filled by the run's own table only @@ -5814,6 +5839,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool shells, merge_friedel, possible, nthreads); for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s]; + sub("st merged walk"); CorrelationCoefficient cc_half_overall; for (const auto &m : merged) { const auto shell = shells.GetShell(m.d); @@ -5829,6 +5855,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } } + sub("st rmeas"); 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; std::vector urmeas_num(n_shells, 0.0), urmeas_den(n_shells, 0.0); @@ -5852,10 +5879,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // n_eff = (sum v)^2 / sum v^2, which is n where the sigmas are equal, and there this is the // ordinary R_meas = sum sqrt(n/(n-1)) sum|I_i - | / sum I_i. // The unweighted sums beside them give the ordinary r_meas; the weighted ones r_meas_weighted. - std::vector &rmeas = merge_rmeas; + FillVector &rmeas = merge_rmeas; rmeas.resize(n_groups); ParallelChunks(n_groups, ThreadsForWork(n_groups, nthreads), [&](int lo, int hi) { - std::fill(rmeas.begin() + lo, rmeas.begin() + hi, RmeasObs{}); + std::fill(rmeas.begin() + lo, rmeas.begin() + hi, RMEAS_EMPTY); }); bool did_gpu_rmeas = false; #ifdef JFJOCH_USE_CUDA @@ -5951,6 +5978,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } } + sub("st anomalous"); // ---- 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 @@ -6212,6 +6240,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool } } + sub("st rest"); MergeStatistics out; out.shells.resize(n_shells); for (int s = 0; s < n_shells; ++s) { @@ -7080,11 +7109,15 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool 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. + std::unique_ptr sub_mark; + auto sub = [&](const char *n) { sub_mark.reset(); sub_mark = std::make_unique(*logger, std::string("rsm/s ") + n); }; + sub("damage B"); 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. It drops nothing itself; the ranges it // segments are the ledger 4c below decides on. + sub("sweep quality"); if (refine_surfaces) MeasureSweepQuality(partial_scaled, cc, cc_n); // CC1/2 of these same fulls with the surfaces not yet folded in, for the caller's cross-pass quality @@ -7105,8 +7138,12 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool "one is judged against fits no surfaces, so this is what the two are compared on)", *cc_half_before_corrections); } + sub("decay"); if (refine_surfaces && refine_decay_b) RefineDecay(n_groups); + decay_terms_kept = {}; + decay_terms_kept_valid = false; + sub("relative B"); if (refine_surfaces && relative_b_deg > 0.0) RefineRelativeB(n_groups); // per-batch relative-B on top of the single decay slope // The three surfaces overlap, so the one fitted first gets the first claim and the others are tested @@ -7122,14 +7159,19 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // The two crystal-frame candidates come last, the smooth spherical-harmonic one before the grid: // the grid is then tested on what the smooth surface left, and where that surface is refused the // grid sees exactly what it saw before. + sub("modulation"); if (refine_surfaces && modulation_iter > 0) RefineModulation(modulation_iter, n_groups); + sub("absorption time"); if (refine_surfaces && absorption_iter > 0) RefineAbsorptionTime(absorption_iter, n_groups); + sub("absorption SH"); if (refine_surfaces && absorption_iter > 0) RefineAbsorptionSH(absorption_iter, n_groups); + sub("absorption gonio"); if (refine_surfaces && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); + sub("disposition"); // --- 4c. The disposition. What keeping each batch of the sweep costs the merged intensities // (delta-CC1/2), which frames cost too much, and what became of every frame of the run. @@ -7155,6 +7197,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool gpu_->SurfaceRelease(); #endif + sub_mark.reset(); surfaces_mark.Stop(); // --- 5. Error model + merge + statistics. --- TimingMark merge_mark(*logger, "rsm merge and statistics"); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index b6dc12a51..fd60744ec 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -546,8 +546,8 @@ private: // same reason as FullsStaging below - the merge runs several times per run and this is tens of // megabytes each time. struct MergeFields { - std::vector group; - std::vector I, sigma, corr, d; + FillVector group; + FillVector I, sigma, corr, d; void Resize(int n) { group.resize(n); I.resize(n); sigma.resize(n); corr.resize(n); d.resize(n); } }; MergeFields merge_fields; @@ -559,9 +559,9 @@ private: // scale->combine->merge chain runs several times per run and these are a few hundred megabytes // between them, so as locals every chain allocates, faults in and zeroes the lot again. struct FullsStaging { - std::vector run, peak, frame, group; - std::vector I, sigma, d, image_number, corr, var_bkg, var_per_I, capture; - std::vector on_ice, clipped; + FillVector run, peak, frame, group; + FillVector I, sigma, d, image_number, corr, var_bkg, var_per_I, capture; + FillVector on_ice, clipped; void Resize(int n) { run.resize(n); peak.resize(n); frame.resize(n); group.resize(n); I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); @@ -576,13 +576,16 @@ private: // FullsStaging - a merge runs several times per run and this is a hundred megabytes between them. // swh_typ: each half-set's precision had every observation been at the run's typical frame scale // (see MergeAndStats, cc_weight). - struct Accum { double swI = 0, sw = 0, swIh[2] = {0, 0}, swh[2] = {0, 0}, swh_typ[2] = {0, 0}; - size_t nh[2] = {0, 0}; float d = NAN; bool on_ice = false; }; - std::vector merge_acc; + // No initialisers of its own, so that the device path - whose unpack writes every field of every + // group - sizes merge_acc without zeroing it first on one thread; the host path assigns EMPTY. + struct Accum { double swI, sw, swIh[2], swh[2], swh_typ[2]; + size_t nh[2]; float d; bool on_ice; }; + static constexpr Accum ACCUM_EMPTY{0, 0, {0, 0}, {0, 0}, {0, 0}, {0, 0}, NAN, false}; + FillVector merge_acc; struct MergeAccumStaging { - std::vector swI, sw, swIh0, swIh1, swh0, swh1, swh_typ0, swh_typ1, d; - std::vector nh0, nh1, rej; - std::vector on_ice; + FillVector swI, sw, swIh0, swIh1, swh0, swh1, swh_typ0, swh_typ1, d; + FillVector nh0, nh1, rej; + FillVector on_ice; void Resize(int n) { swI.resize(n); sw.resize(n); swIh0.resize(n); swIh1.resize(n); swh0.resize(n); swh1.resize(n); swh_typ0.resize(n); swh_typ1.resize(n); @@ -600,12 +603,14 @@ private: std::vector merge_obs_key; std::vector merge_wobs; // The R_meas accumulators per group (see MergeAndStats' statistics), and the device's per-group sums. - struct RmeasObs { double sum_abs_dev = 0, sum_I = 0, w_abs_dev = 0, w_I = 0, sum_v = 0, sum_v2 = 0; - int n = 0, shell = -1; }; - std::vector merge_rmeas; + // No initialisers, as Accum: the merge resets them on all threads to RMEAS_EMPTY. + struct RmeasObs { double sum_abs_dev, sum_I, w_abs_dev, w_I, sum_v, sum_v2; + int n, shell; }; + static constexpr RmeasObs RMEAS_EMPTY{0, 0, 0, 0, 0, 0, 0, -1}; + FillVector merge_rmeas; struct RmeasStaging { - std::vector absdev, sumI, wabsdev, wsumI, sumv, sumv2; - std::vector n, nusable; + FillVector absdev, sumI, wabsdev, wsumI, sumv, sumv2; + FillVector n, nusable; void Resize(int g) { absdev.resize(g); sumI.resize(g); wabsdev.resize(g); wsumI.resize(g); sumv.resize(g); sumv2.resize(g); n.resize(g); nusable.resize(g); @@ -615,8 +620,8 @@ private: // 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; }; - std::vector asymptote_scatter; + struct GroupScatter { double sum, sum_sq, sum_var; int n; }; // no initialisers, as Accum + FillVector asymptote_scatter; // The error model's working pools (see MergeAndStats): the samples themselves, the scratch copy each // fit partitions, the subset inside the resolution cutoff, and @@ -919,6 +924,10 @@ private: std::vector g_start; }; void BuildDecayTerms(int n_groups, DecayTerms &t) const; + // The terms the radiation-damage monitor built, handed on to RefineDecay: nothing between the two + // touches corr, so they are the terms RefineDecay would build. Empty outside that stretch of RunArm. + DecayTerms decay_terms_kept; + bool decay_terms_kept_valid = false; // Solve a smooth per-batch relative-B from the per-batch normal equations for b (num_c, den_c): // data-fidelity + a second-difference (curvature) penalty, by Gauss-Seidel, each batch clamped to // +-b_max. Returns the un-anchored curve; the caller sets the gauge and the clamp it can live with.