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) <noreply@anthropic.com>
2026-10-14 10:09:23 +02:00
parent 70276a69fe
commit 34d5bb5631
6 changed files with 269 additions and 96 deletions
@@ -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<int>(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 <typename T, typename Reduce, typename Consume>
+    void ReduceImagesInOrder(size_t n_images, size_t nthreads, Reduce reduce, Consume consume) {
+        constexpr size_t BLOCK_IMAGES = 4096;
+        std::vector<std::vector<T>> out(std::min(BLOCK_IMAGES, n_images));
+        for (size_t first = 0; first < n_images; first += BLOCK_IMAGES) {
+            const int n = static_cast<int>(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<UnitCell> &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<Contribution> &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<Contribution> &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<Contribution> contributions;
+    reject_count += Contributions(outcome, contributions);
+    Accumulate(contributions, image_id);
+}
+
+void MergeOnTheFly::AddImages(const std::vector<IntegrationOutcome> &outcomes) {
+    std::unique_lock ul(merged_mutex);
+
+    std::vector<size_t> rejected(outcomes.size(), 0);
+    ReduceImagesInOrder<Contribution>(
+        outcomes.size(), nthreads,
+        [&](size_t i, std::vector<Contribution> &out) {
+            if (!Mask(outcomes[i]))
+                rejected[i] = Contributions(outcomes[i], out);
+        },
+        [&](size_t i, const std::vector<Contribution> &out) {
+            reject_count += rejected[i];
+            Accumulate(out, static_cast<int64_t>(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<IntegrationOutcome> &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<uint64_t, std::vector<Obs>> 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<KeyedObs>(
+        outcomes.size(), nthreads,
+        [&](size_t i, std::vector<KeyedObs> &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<KeyedObs> &out) {
+            for (const auto &o: out)
+                groups[o.key].push_back(o.obs);
+        });

    // --- 2. One global pool of (sigma^2, <I>^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<IntegrationOutcome> &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<MergedReflection> MergeOnTheFly::ExportReflections() {
+std::vector<MergedReflection> MergeOnTheFly::ExportIntensities() {
    std::unique_lock ul(merged_mutex);

    std::vector<MergedReflection> out;
+    std::vector<std::pair<uint64_t, const MergeAccum *>> 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<MergedReflection> MergeOnTheFly::ExportReflections() {

        out.emplace_back(mr);
    }
+    return out;
+}

+std::vector<MergedReflection> MergeOnTheFly::ExportReflections() {
+    auto out = ExportIntensities();
    AssignRfreeFlags(out, space_group, scaling_settings.GetRfreeFraction(), 500, reference_cell);
    ApplyFrenchWilson(out, space_group);
    return out;
 }

 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
-                                       const std::vector<IntegrationOutcome> &integration_outcome) {
+                                       const std::vector<IntegrationOutcome> &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<int64_t>(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<MergedReflection> &m
    std::unordered_map<uint64_t, RmeasObs> 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<ShellObs>(
+        integration_outcome.size(), nthreads,
+        [&](size_t i, std::vector<ShellObs> &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<double>(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<ShellObs> &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<double>(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 - <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)
@@ -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<uint64_t, MergeAccum> accumulator;
+    std::unordered_map<uint64_t, MergeAccum> accumulator;   // exported in key order

    // Global error model (XDS form): sigma_corr^2 = a*sigma^2 + (b*<I>)^2. a rescales the
    // (under-estimated) counting variance; the (b*<I>)^2 term adds the intensity-
@@ -315,12 +316,26 @@ class MergeOnTheFly {
    std::unordered_map<uint64_t, float> 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<Contribution> &out) const;
+    void Accumulate(const std::vector<Contribution> &contributions, int64_t image_id);
 public:
    MergeOnTheFly(const DiffractionExperiment &x);
    MergeOnTheFly& ReferenceCell(const std::optional<UnitCell> &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<size_t>(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<IntegrationOutcome> &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<ReportResolutionRange> declared_range = std::nullopt);

    std::vector<MergedReflection> ExportReflections();
+    // The merged intensities alone - no free-R flags, no French-Wilson amplitudes - which is all a
+    // scaling reference reads.
+    std::vector<MergedReflection> ExportIntensities();
 };

 std::vector<MergedReflection> MergeAll(const DiffractionExperiment &x,
-                                       const std::vector<IntegrationOutcome> &reflections);
+                                       const std::vector<IntegrationOutcome> &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
@@ -358,7 +358,7 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si

    for (int outer = 0; outer < settings_.outer_iterations; ++outer) {
        // Reference full intensities from the current corrections.
-        const std::vector<MergedReflection> merged = MergeAll(experiment_, outcomes);
+        const std::vector<MergedReflection> merged = MergeAll(experiment_, outcomes, nthreads);
        std::map<HKLKey, double> reference;
        for (const auto &m: merged)
            reference[hkl_key_generator_(m)] = m.I;
@@ -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<IntegrationOutcome> &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<int64_t>(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.
@@ -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<int64_t>(i));
+            merge_engine.AddImages(reflections);
            merged_reflections = merge_engine.ExportReflections();

            // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual
@@ -3,6 +3,7 @@

 #include <catch2/catch_all.hpp>

+#include <cstring>
 #include <random>
 #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<int> index(-6, 6);
+    std::uniform_real_distribution<float> unit(0.0f, 1.0f);
+    std::vector<IntegrationOutcome> 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<float>(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<int64_t>(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], &parallel_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);
+}