From 35607057d9da72e869178153f0407cb09c059672 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 1 Jul 2026 07:25:34 +0200 Subject: [PATCH] CC1/2 merge: assign half-sets by deterministic per-image hash MergeOnTheFly::AddImage picked each image's CC1/2 half from a shared mt19937 drawn in call order (and before Mask), so the split depended on iteration/thread order and on how many images were masked. The class is mutex-guarded for concurrent "on-the-fly" use, so any parallel merge would make CC1/2 non-reproducible - a latent race. Assign the half as a splitmix64 hash of the image's stable index instead, computed after Mask. The split is now reproducible run-to-run, independent of AddImage order, parallel-safe, and decoupled from masking. Callers pass the outcome's vector index as the image id. Verified: lyso_ref two-pass -M -P rot3d gives identical CC1/2 across runs (overall 99.6%, P41212); hash split is balanced ~50/50. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/scale_merge/Merge.cpp | 22 ++++++++++++++++++---- image_analysis/scale_merge/Merge.h | 10 ++++------ process/JFJochProcess.cpp | 4 ++-- tools/jfjoch_scale.cpp | 4 ++-- 4 files changed, 26 insertions(+), 14 deletions(-) diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 3f9c1466..bbab6f4d 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -17,6 +17,19 @@ #include "../../common/ResolutionShells.h" #include "HKLKey.h" +namespace { + // Deterministic CC1/2 half-set assignment: a splitmix64 bit-mix of the image's stable index. + // A pure function of image identity (not a draw from a shared RNG in call order) keeps the split + // reproducible run-to-run, independent of AddImage call order, and safe under concurrent merging. + int HalfForImage(int64_t image_id) { + uint64_t z = static_cast(image_id) + 0x9e3779b97f4a7c15ULL; + z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL; + z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL; + z = z ^ (z >> 31); + return static_cast(z & 1ULL); + } +} + MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) : space_group_number(x.GetSpaceGroupNumber().value_or(1)), scaling_settings(x.GetScalingSettings()), @@ -41,13 +54,14 @@ MergeOnTheFly &MergeOnTheFly::ReferenceCell(const std::optional &cell) return *this; } -void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, bool cc_mask) { +void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id, bool cc_mask) { std::unique_lock ul(merged_mutex); - const int half = half_dist(rng); if (Mask(outcome, cc_mask)) return; + const int half = HalfForImage(image_id); + for (const auto &r: outcome.reflections) { if (generator.IsSystematicallyAbsent(r)) continue; @@ -481,8 +495,8 @@ std::vector MergeAll(const DiffractionExperiment &x, const std::vector &integration_outcome, bool mask) { MergeOnTheFly merge(x); - for (const auto &outcome: integration_outcome) - merge.AddImage(outcome, mask); + for (size_t i = 0; i < integration_outcome.size(); ++i) + merge.AddImage(integration_outcome[i], static_cast(i), mask); return merge.ExportReflections(); } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 1c2aa265..b6095ffb 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -70,11 +70,6 @@ class MergeOnTheFly { HKLKeyGenerator generator; - // To select images for half-datasets to calculate CC1/2, I use a random number generator with a fixed seed. - // This makes sure that images are selected randomly, but in a fully reproducible manner (at least for the same binary) - std::mt19937 rng{123456789u}; - std::bernoulli_distribution half_dist{0.5}; - std::map accumulator; // Global error model (XDS form): sigma_corr^2 = a*sigma^2 + (b*)^2. a rescales the @@ -118,7 +113,10 @@ public: // Outlier rejection (driven by ScalingSettings::GetOutlierRejectNsigma) reports its count. [[nodiscard]] size_t RejectedCount() const { return reject_count; } - void AddImage(const IntegrationOutcome& outcome, bool cc_mask = false); + // image_id is the image's stable identity (its index in the outcomes vector). The CC1/2 half-set + // is a deterministic hash of it, so the split is reproducible run-to-run and independent of the + // order (or threading) of AddImage calls - not a draw from a shared RNG in call order. + void AddImage(const IntegrationOutcome& outcome, int64_t image_id, bool cc_mask = false); // Per-crystal CC1/2-delta rejection (CrystFEL deltaCChalf): returns a per-image flag // marking images whose removal would raise CC1/2 by a low-side outlier amount diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index e206ea18..9a6f424d 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -616,8 +616,8 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { merge_engine.ErrorModelB(), merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0, merge_engine.ErrorModelChi2()); - for (const auto &outcome: merge_input) - merge_engine.AddImage(outcome); + for (size_t i = 0; i < merge_input.size(); ++i) + merge_engine.AddImage(merge_input[i], static_cast(i)); ScaleMergeResult out; out.merged = merge_engine.ExportReflections(); diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 1f992e8a..ef0404e8 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -377,8 +377,8 @@ int main(int argc, char **argv) { MergeOnTheFly merge_engine(experiment); merge_engine.ReferenceCell(experiment.GetUnitCell()); - for (auto &i : merge_input) - merge_engine.AddImage(i); + for (size_t i = 0; i < merge_input.size(); ++i) + merge_engine.AddImage(merge_input[i], static_cast(i)); auto merged_reflections = merge_engine.ExportReflections(); auto merged_statistics = merge_engine.MergeStats(merged_reflections, merge_input, reference_data);