From 34d5bb563138055c31c3cd9ceb482f4871eb3ab5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 14 Oct 2026 10:09:23 +0200 Subject: [PATCH] stills merge 2.3-2.5x faster, output unchanged: parallel observation passes, hash accumulator, no amplitudes on scaling references The stills scale/merge after the image loop was almost entirely serial: - every scaling-reference merge (MergeAll: the per-image scaling reference and the two partiality post-refinement references) ran French-Wilson single-threaded on its export, about 3 s each, and assigned free-R flags - neither is read by a reference, which uses I only. MergeAll now exports intensities alone (MergeOnTheFly::ExportIntensities). - the merge accumulated into a std::map; it is now an unordered_map exported in key order, which is the std::map order. - AddImage, the error-model collection and the R_meas pass of MergeStats did the per-observation work (systematic-absence test, ASU key, filters, error-model sigma, outlier test, merged-I lookup) in one serial loop. They now reduce a block of images in parallel and consume it in image order (ReduceImagesInOrder), so every floating-point sum, half-set choice and hash-map insertion is made in the serial order. Measured on a 25k-frame XFEL stills run (box loaded, load 9-14): merge P1 19.0 -> 7.1 s, final-group merge 27.7 -> 13.2 s; on a second run 23.4 -> 7.7 s and 29.7 -> 13.6 s. mtz/unmerged mtz data, hkl, cif and report identical on both runs, mx and --mode scale. New test MergeOnTheFly_AddImagesIsTheSerialMerge checks the parallel merge against the serial loop across a block boundary. Co-Authored-By: Claude Opus 5.5 (1M context) --- image_analysis/scale_merge/Merge.cpp | 253 ++++++++++++------ image_analysis/scale_merge/Merge.h | 28 +- .../scale_merge/StillsPartialityRefine.cpp | 2 +- rugnux/RugnuxScaleMerge.cpp | 11 +- rugnux/rugnux_cli.cpp | 6 +- tests/MergeScaleTest.cpp | 65 +++++ 6 files changed, 269 insertions(+), 96 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 54ec8d9eb..2c3bf37f1 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -15,6 +15,7 @@ #include "../../common/CorrelationCoefficient.h" #include "../../common/ResolutionShells.h" #include "../../common/Definitions.h" +#include "../../common/ParallelFor.h" #include "HKLKey.h" #include "RfreeFlags.h" #include "FrenchWilson.h" @@ -32,6 +33,27 @@ namespace { z = z ^ (z >> 31); return static_cast(z & 1ULL); } + + // Calls reduce(i, out) for the images in parallel, a block of them at a time, and then consume(i, out) + // for each image of the block in image order - the order that every floating-point sum, every + // half-set choice and every hash-map insertion of the merge depends on. The result is therefore the + // serial loop's, bit for bit, at any thread count. + template + void ReduceImagesInOrder(size_t n_images, size_t nthreads, Reduce reduce, Consume consume) { + constexpr size_t BLOCK_IMAGES = 4096; + std::vector> out(std::min(BLOCK_IMAGES, n_images)); + for (size_t first = 0; first < n_images; first += BLOCK_IMAGES) { + const int n = static_cast(std::min(BLOCK_IMAGES, n_images - first)); + ParallelChunks(n, nthreads, [&](int lo, int hi) { + for (int j = lo; j < hi; ++j) { + out[j].clear(); + reduce(first + j, out[j]); + } + }); + for (int j = 0; j < n; ++j) + consume(first + j, out[j]); + } + } } MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) @@ -57,14 +79,8 @@ MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional &cell) return *this; } -void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id) { - std::unique_lock ul(merged_mutex); - - if (Mask(outcome)) - return; - - const int prefer = HalfForImage(image_id); - +size_t MergeOnTheFly::Contributions(const IntegrationOutcome &outcome, std::vector &out) const { + size_t rejected = 0; for (const auto &r: outcome.reflections) { if (generator.IsSystematicallyAbsent(r)) continue; @@ -94,24 +110,36 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id const auto mit = reject_median_I.find(hkl_key); if (mit != reject_median_I.end() && std::fabs(I_corr - mit->second) > reject_nsigma * sigma_corr) { - ++reject_count; + ++rejected; continue; } } - auto it = accumulator.find(hkl_key); - if (it == accumulator.end()) - it = accumulator.emplace(hkl_key, MergeAccum{ - .h = hkl.plus ? hkl.h : -hkl.h, - .k = hkl.plus ? hkl.k : -hkl.k, - .l = hkl.plus ? hkl.l : -hkl.l, - }).first; - const float w = 1.0f / (sigma_corr * sigma_corr); - const float wI = w * I_corr; + out.push_back(Contribution{ + .key = hkl_key, + .h = hkl.plus ? hkl.h : -hkl.h, + .k = hkl.plus ? hkl.k : -hkl.k, + .l = hkl.plus ? hkl.l : -hkl.l, + .w = w, + .wI = w * I_corr, + .d = r.d, + .on_ice_ring = r.on_ice_ring + }); + } + return rejected; +} - it->second.sum_wI += wI; - it->second.sum_w += w; +void MergeOnTheFly::Accumulate(const std::vector &contributions, int64_t image_id) { + const int prefer = HalfForImage(image_id); + + for (const auto &c: contributions) { + auto it = accumulator.find(c.key); + if (it == accumulator.end()) + it = accumulator.emplace(c.key, MergeAccum{.h = c.h, .k = c.k, .l = c.l}).first; + + it->second.sum_wI += c.wI; + it->second.sum_w += c.w; // CC1/2 correlates the two half-set means, so a reflection whose observations all land in one // half has no second mean and drops out of the statistic entirely. Taking the image's own @@ -119,18 +147,44 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id // reflection floor(n/2) against ceil(n/2), which is the split cctbx and XDS use. const int half = it->second.n_half[prefer] > it->second.n_half[1 - prefer] ? 1 - prefer : prefer; - it->second.sum_wI_half[half] += wI; - it->second.sum_w_half[half] += w; + it->second.sum_wI_half[half] += c.wI; + it->second.sum_w_half[half] += c.w; it->second.n_half[half]++; - if (r.on_ice_ring) + if (c.on_ice_ring) it->second.on_ice_ring = true; - if (!std::isfinite(it->second.d) && std::isfinite(r.d) && r.d > 0.0f) - it->second.d = r.d; + if (!std::isfinite(it->second.d) && std::isfinite(c.d) && c.d > 0.0f) + it->second.d = c.d; } } +void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id) { + std::unique_lock ul(merged_mutex); + + if (Mask(outcome)) + return; + std::vector contributions; + reject_count += Contributions(outcome, contributions); + Accumulate(contributions, image_id); +} + +void MergeOnTheFly::AddImages(const std::vector &outcomes) { + std::unique_lock ul(merged_mutex); + + std::vector rejected(outcomes.size(), 0); + ReduceImagesInOrder( + outcomes.size(), nthreads, + [&](size_t i, std::vector &out) { + if (!Mask(outcomes[i])) + rejected[i] = Contributions(outcomes[i], out); + }, + [&](size_t i, const std::vector &out) { + reject_count += rejected[i]; + Accumulate(out, static_cast(i)); + }); +} + float MergeOnTheFly::CorrectedSigma(float I_corr, float sigma_corr, float image_scale_corr, float var_bkg, uint64_t hkl_key) const { @@ -188,29 +242,36 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc // --- 1. Collect accepted, scaled observations grouped by symmetry-equivalent hkl, // applying exactly the filters AddImage uses. --- struct Obs { float I, sigma; }; + struct KeyedObs { uint64_t key; Obs obs; }; std::unordered_map> groups; - for (const auto &outcome: outcomes) { - if (Mask(outcome)) - continue; - for (const auto &r: outcome.reflections) { - if (generator.IsSystematicallyAbsent(r)) - continue; - if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) - continue; - if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) - continue; - if (exclude_ice_rings && r.on_ice_ring) - continue; - if (r.partiality < min_partiality) - continue; - const float I_corr = r.I * r.image_scale_corr; - const float sigma_corr = r.sigma * r.image_scale_corr; - if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) - continue; - groups[generator(r).pack()].push_back({I_corr, sigma_corr}); - } - } + ReduceImagesInOrder( + outcomes.size(), nthreads, + [&](size_t i, std::vector &out) { + if (Mask(outcomes[i])) + return; + for (const auto &r: outcomes[i].reflections) { + if (generator.IsSystematicallyAbsent(r)) + continue; + if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) + continue; + if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) + continue; + if (exclude_ice_rings && r.on_ice_ring) + continue; + if (r.partiality < min_partiality) + continue; + const float I_corr = r.I * r.image_scale_corr; + const float sigma_corr = r.sigma * r.image_scale_corr; + if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) + continue; + out.push_back({generator(r).pack(), {I_corr, sigma_corr}}); + } + }, + [&](size_t, const std::vector &out) { + for (const auto &o: out) + groups[o.key].push_back(o.obs); + }); // --- 2. One global pool of (sigma^2, ^2, bias-corrected squared deviation). For an // observation in a group of n, the residual from the inverse-variance mean has @@ -337,7 +398,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc error_model_chi2 = chi2.empty() ? 0.0 : median(chi2) / CHI2_1_MEDIAN; } -bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome) { +bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome) const { if (reference_cell) { auto cell = outcome.latt.GetUnitCell(); if (!cell.is_close(*reference_cell, @@ -355,12 +416,20 @@ bool MergeOnTheFly::Mask(const IntegrationOutcome &outcome) { return false; } -std::vector MergeOnTheFly::ExportReflections() { +std::vector MergeOnTheFly::ExportIntensities() { std::unique_lock ul(merged_mutex); std::vector out; + std::vector> by_key; + by_key.reserve(accumulator.size()); + for (const auto &[key, accum]: accumulator) + by_key.emplace_back(key, &accum); + std::sort(by_key.begin(), by_key.end(), + [](const auto &a, const auto &b) { return a.first < b.first; }); + out.reserve(accumulator.size()); - for (const auto &accum: accumulator | std::views::values) { + for (const auto &[key, accum_ptr]: by_key) { + const MergeAccum &accum = *accum_ptr; if (accum.sum_w <= 0.0) continue; @@ -388,18 +457,22 @@ std::vector MergeOnTheFly::ExportReflections() { out.emplace_back(mr); } + return out; +} +std::vector MergeOnTheFly::ExportReflections() { + auto out = ExportIntensities(); AssignRfreeFlags(out, space_group, scaling_settings.GetRfreeFraction(), 500, reference_cell); ApplyFrenchWilson(out, space_group); return out; } std::vector MergeAll(const DiffractionExperiment &x, - const std::vector &integration_outcome) { + const std::vector &integration_outcome, + size_t nthreads) { MergeOnTheFly merge(x); - for (size_t i = 0; i < integration_outcome.size(); ++i) - merge.AddImage(integration_outcome[i], static_cast(i)); - return merge.ExportReflections(); + merge.Threads(nthreads).AddImages(integration_outcome); + return merge.ExportIntensities(); } struct ShellAccum { @@ -630,43 +703,51 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m std::unordered_map rmeas_obs; rmeas_obs.reserve(merged.size()); - for (int i = 0; i < integration_outcome.size(); ++i) { - if (Mask(integration_outcome[i])) - continue; + struct ShellObs { uint64_t key; float I_corr; float merged_I; int shell; bool merged; }; + ReduceImagesInOrder( + integration_outcome.size(), nthreads, + [&](size_t i, std::vector &out) { + if (Mask(integration_outcome[i])) + return; + for (const auto &r: integration_outcome[i].reflections) { + if (generator.IsSystematicallyAbsent(r)) + continue; + if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) + continue; + if (!AcceptReflection(r, d_min_limit_A, d_max_limit_A)) + continue; + if (r.partiality < min_partiality) + continue; - for (const auto &r: integration_outcome[i].reflections) { - if (generator.IsSystematicallyAbsent(r)) - continue; - if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) - continue; - if (!AcceptReflection(r, d_min_limit_A, d_max_limit_A)) - continue; - if (r.partiality < min_partiality) - continue; + const float I_corr = r.I * r.image_scale_corr; + const float sigma_corr = r.sigma * r.image_scale_corr; + if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) + continue; - const float I_corr = r.I * r.image_scale_corr; - const float sigma_corr = r.sigma * r.image_scale_corr; - if (!std::isfinite(I_corr) || !std::isfinite(sigma_corr) || sigma_corr <= 0.0f) - continue; - - const auto shell = shells.GetShell(r.d); - if (!shell.has_value()) - continue; - const int s = *shell; - if (s >= 0 && s < n_shells) { - acc[s].total_obs++; - const auto key = generator(r).pack(); - const auto mit = merged_I.find(key); - if (mit != merged_I.end()) { - auto &ra = rmeas_obs[key]; - ra.sum_abs_dev += std::abs(static_cast(I_corr) - mit->second); - ra.sum_I += I_corr; - ra.n++; - ra.shell = s; + const auto shell = shells.GetShell(r.d); + if (!shell.has_value()) + continue; + const int s = *shell; + if (s >= 0 && s < n_shells) { + const auto key = generator(r).pack(); + const auto mit = merged_I.find(key); + out.push_back({key, I_corr, mit != merged_I.end() ? mit->second : NAN, s, + mit != merged_I.end()}); } } - } - } + }, + [&](size_t, const std::vector &out) { + for (const auto &o: out) { + acc[o.shell].total_obs++; + if (o.merged) { + auto &ra = rmeas_obs[o.key]; + ra.sum_abs_dev += std::abs(static_cast(o.I_corr) - o.merged_I); + ra.sum_I += o.I_corr; + ra.n++; + ra.shell = o.shell; + } + } + }); // R_meas per shell: sum over reflections of sqrt(n/(n-1)) * sum_i|I_i - |, over sum of I_i. // R_merge and R_pim beside it, with the factor 1 and sqrt(1/(n-1)) - R_pim following Weiss (2001) diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index b8f2b797e..7cca243a3 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -283,10 +283,11 @@ class MergeOnTheFly { // merged intensities feed the space-group search and the error model - those model fits must not // see the ice-contaminated intensities. The final in-symmetry merge keeps them (for completeness). bool exclude_ice_rings = false; + size_t nthreads = 1; HKLKeyGenerator generator; - std::map accumulator; + std::unordered_map accumulator; // exported in key order // Global error model (XDS form): sigma_corr^2 = a*sigma^2 + (b*)^2. a rescales the // (under-estimated) counting variance; the (b*)^2 term adds the intensity- @@ -315,12 +316,26 @@ class MergeOnTheFly { std::unordered_map reject_median_I; size_t reject_count = 0; - bool Mask(const IntegrationOutcome &outcome); + bool Mask(const IntegrationOutcome &outcome) const; + + // One accepted observation, reduced to what the accumulator adds for it. + struct Contribution { + uint64_t key; + int32_t h, k, l; // the reflection's own index, as the accumulator records it + float w, wI, d; + bool on_ice_ring; + }; + // The observations of one image that reach the merge; returns how many the outlier test dropped. + size_t Contributions(const IntegrationOutcome &outcome, std::vector &out) const; + void Accumulate(const std::vector &contributions, int64_t image_id); public: MergeOnTheFly(const DiffractionExperiment &x); MergeOnTheFly& ReferenceCell(const std::optional &cell); MergeOnTheFly& ExcludeIceRings(bool input) { exclude_ice_rings = input; return *this; } MergeOnTheFly& FilterByImageCC(bool input) { filter_by_image_cc = input; return *this; } + // Workers for the passes over every observation (AddImages, RefineErrorModel, MergeStats); the + // result does not depend on the number. + MergeOnTheFly& Threads(size_t input) { nthreads = std::max(1, input); return *this; } // Fit the global error model from the spread of symmetry-equivalent observations. // Call once before merging; AddImage then applies it. @@ -341,6 +356,9 @@ public: // because it holds every observation at once; a streaming accumulator cannot, and there is no // second, device-side walk of the stills observations for this to disagree with. void AddImage(const IntegrationOutcome& outcome, int64_t image_id); + // AddImage for every outcome, the image id being its index: the same merge, bit for bit, with the + // per-observation work spread over the workers and the accumulation kept in image order. + void AddImages(const std::vector &outcomes); // d_min_override, when set, is the effective high-resolution limit for the shell table (used for // the automatic resolution cutoff computed by the caller); otherwise the manual @@ -355,10 +373,14 @@ public: std::optional declared_range = std::nullopt); std::vector ExportReflections(); + // The merged intensities alone - no free-R flags, no French-Wilson amplitudes - which is all a + // scaling reference reads. + std::vector ExportIntensities(); }; std::vector MergeAll(const DiffractionExperiment &x, - const std::vector &reflections); + const std::vector &reflections, + size_t nthreads); // Pearson CC between one image's corrected intensities (I * image_scale_corr) and a reference set of // full intensities, over the reflections that would enter the merge (non-ice, within the resolution diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 44910f6cb..12265a510 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -358,7 +358,7 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. - const std::vector merged = MergeAll(experiment_, outcomes); + const std::vector merged = MergeAll(experiment_, outcomes, nthreads); std::map reference; for (const auto &m: merged) reference[hkl_key_generator_(m)] = m.I; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index a2b5f44a4..e42731bfa 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1250,7 +1250,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // degrade weak stills (measured CC1/2 collapse at the default 3 iters). { phase("Scaling images (" + label + ")"); - auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome()); + TimingMark scale_mark(logger, "stills scaling images " + label); + auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), config_.nthreads); indexer->ScaleAllImages(merge_result); } // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal @@ -1258,6 +1259,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // correction (no re-integration). Never for the P1 search pass. if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) { phase("Partiality post-refine (" + label + ")"); + TimingMark refine_mark(logger, "stills partiality post-refine " + label); StillsPartialityRefine refiner(experiment_); const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads); logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); @@ -1272,7 +1274,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const std::vector &merge_input = indexer->GetIntegrationOutcome(); phase("Merging"); + TimingMark merge_mark(logger, "stills error model and merge " + label); MergeOnTheFly merge_engine(experiment_); + merge_engine.Threads(config_.nthreads); // For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged // intensities and the error model, so the space-group search sees clean data; the final // in-symmetry merge keeps them so completeness is not lost. @@ -1296,11 +1300,11 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}", em.a, em.b, em.isa, merge_engine.ErrorModelChi2()); } - for (size_t i = 0; i < merge_input.size(); ++i) - merge_engine.AddImage(merge_input[i], static_cast(i)); + merge_engine.AddImages(merge_input); ScaleMergeResult out; out.merged = merge_engine.ExportReflections(); + merge_mark.Stop(); // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full @@ -1313,6 +1317,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { cut_ss.GetResolutionCCTarget(), for_search, logger, &result.resolution_fit_A); phase("Computing statistics"); + TimingMark stats_mark(logger, "stills statistics " + label); out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min); // The same kept reflections binned over the reference range (--report-resolution); the // stills error model is fitted over the whole range already, so no ISa of its own. diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 714bd73d1..5527130e9 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2103,7 +2103,7 @@ static int RunRugnux(int argc, char **argv) { // against a foreign dataset injects cross-dataset systematics and is a worse reference than // the data's own merge. The per-image scale G is the exact one-pass solution, so one pass // (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills). - ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads); + ScaleOnTheFly(experiment, MergeAll(experiment, reflections, nthreads)).Scale(reflections, nthreads); // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal // orientation tilt against the merge and recompute each reflection's partiality + scale correction // (no re-integration), then merge with the improved corrections. @@ -2113,6 +2113,7 @@ static int RunRugnux(int argc, char **argv) { logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); } MergeOnTheFly merge_engine(experiment); + merge_engine.Threads(nthreads); merge_engine.ReferenceCell(experiment.GetUnitCell()); // --min-image-cc has to hold for the merge itself, not only for the reported statistics. merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0); @@ -2128,8 +2129,7 @@ static int RunRugnux(int argc, char **argv) { logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}", em.a, em.b, em.isa, merge_engine.ErrorModelChi2()); } - for (size_t i = 0; i < reflections.size(); ++i) - merge_engine.AddImage(reflections[i], static_cast(i)); + merge_engine.AddImages(reflections); merged_reflections = merge_engine.ExportReflections(); // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index aab3d5e9d..ed548a294 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -3,6 +3,7 @@ #include +#include #include #include "../image_analysis/bragg_integration/CalcISigma.h" #include "../image_analysis/scale_merge/ErrorModel.h" @@ -585,3 +586,67 @@ TEST_CASE("Global Wilson B over a widened window, with no upper cap", "[wilson]" CHECK(1.0 / (3.9 * 3.9) - 1.0 / (low_res.d_low_A * low_res.d_low_A) == Catch::Approx(1.0 / 9.0 - 1.0 / 16.0).epsilon(0.01)); } + +// The parallel AddImages is the serial AddImage loop, bit for bit: the per-observation work is spread +// over the workers, but every sum and every half-set choice is still made in image order. More images +// than one of its blocks, so the block boundary is crossed. +TEST_CASE("MergeOnTheFly_AddImagesIsTheSerialMerge") { + DiffractionExperiment x; + x.SetSpaceGroup(TestSpaceGroup()); + + std::mt19937 rng(7); + std::uniform_int_distribution index(-6, 6); + std::uniform_real_distribution unit(0.0f, 1.0f); + std::vector outcomes(5000); + for (auto &o: outcomes) { + for (int j = 0; j < 20; ++j) { + Reflection r{}; + r.h = index(rng); + r.k = index(rng); + r.l = index(rng); + if (r.h == 0 && r.k == 0 && r.l == 0) + continue; + r.d = 20.0f / std::sqrt(static_cast(r.h * r.h + r.k * r.k + r.l * r.l)); + r.I = 100.0f * unit(rng) - 5.0f; + r.sigma = 1.0f + 5.0f * unit(rng); + r.var_bkg = r.sigma * r.sigma * unit(rng); + r.prescaling_corr = 1.0f; + r.partiality = 0.5f + 0.5f * unit(rng); + r.image_scale_corr = 0.5f + unit(rng); + r.on_ice_ring = unit(rng) < 0.05f; + o.reflections.push_back(r); + } + } + + MergeOnTheFly serial(x); + serial.RefineErrorModel(outcomes); + for (size_t i = 0; i < outcomes.size(); ++i) + serial.AddImage(outcomes[i], static_cast(i)); + const auto serial_merged = serial.ExportReflections(); + + MergeOnTheFly parallel(x); + parallel.Threads(8); + parallel.RefineErrorModel(outcomes); + parallel.AddImages(outcomes); + const auto parallel_merged = parallel.ExportReflections(); + + REQUIRE(serial.ErrorModelActive()); + CHECK(serial.ErrorModelA() == parallel.ErrorModelA()); + CHECK(serial.ErrorModelB() == parallel.ErrorModelB()); + REQUIRE(serial_merged.size() == parallel_merged.size()); + for (size_t i = 0; i < serial_merged.size(); ++i) { + CHECK(serial_merged[i].h == parallel_merged[i].h); + CHECK(serial_merged[i].k == parallel_merged[i].k); + CHECK(serial_merged[i].l == parallel_merged[i].l); + CHECK(serial_merged[i].I == parallel_merged[i].I); + CHECK(serial_merged[i].sigma == parallel_merged[i].sigma); + for (int half = 0; half < 2; ++half) // NaN where a half is empty, on both sides + CHECK(std::memcmp(&serial_merged[i].I_half[half], ¶llel_merged[i].I_half[half], sizeof(float)) == 0); + } + + const auto serial_stats = serial.MergeStats(serial_merged, outcomes); + const auto parallel_stats = parallel.MergeStats(parallel_merged, outcomes); + CHECK(serial_stats.overall.total_observations == parallel_stats.overall.total_observations); + CHECK(serial_stats.overall.cc_half == parallel_stats.overall.cc_half); + CHECK(serial_stats.overall.r_meas == parallel_stats.overall.r_meas); +}