Unmerged MTZ: build the file on several workers, byte-identical

The unmerged export (_unmerged.mtz) was serial apart from its final row
sort: gathering the partials, sorting them by (h,k,l,image), summing the
rocking events and filling the rows all ran on one thread. On a 17 M
partial rotation set that was 2.8 s (sort 1.6 s, event sums 0.6 s,
gather 0.25 s, rows 0.19 s); on a synthetic 122 M partial sweep 23 s.

- The partials are gathered into preallocated slots per outcome, left
  uninitialised so the workers first-touch the pages, and sorted with
  ParallelSort. The sort key now ends with the part's input position
  (outcome, index), so the order is total and the same on any number
  of workers. It differs from the previous std::sort order only among
  parts with equal (h,k,l,image_number); measured on four in-house
  sweeps (464-6138 such positions each), every summed full is bitwise
  the same.
- Events are summed in chunks whose bounds move forward to the next
  event start, so no event is split; the chunks are returned as pieces
  in order rather than concatenated.
- Rows are filled by index into a preallocated table; absences and the
  batch range are still decided in one serial pass. The final row
  gather after the sort is parallel too.

Measured: the 17 M set 2.82 s -> 0.86 s, the 8 M set 1.37 -> 0.49 s,
all while the P1 cross-check merge shares the pool; a synthetic 122 M
partial sweep 22.7 s -> 7.7 s at 32 workers. What remains is mostly
ParallelSort's pairwise merge levels. The file write (0.1 s for 112 MB)
is I/O and left as it is.

p.hkl, p.mtz, p_P1.mtz and p_unmerged.mtz are md5-identical to the base
on four in-house sets, on the GPU and the CPU build, and at -N 1.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-09-28 16:52:13 +02:00
1 parent 351be7de0f
commit aa6ec193e9
2 files changed
+196 -110

No files matched your search

+48
View File
@@ -5,6 +5,7 @@
#include <filesystem>
#include <string>
#include <cstring>
#include <gemmi/mtz.hpp>
@@ -137,3 +138,50 @@ TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization, QE the sensor efficiency, F
CHECK(I == Catch::Approx(raw[i] * lp[i] * qe[i] * fl[i]).epsilon(1e-5));
}
}
TEST_CASE("Unmerged MTZ: built on several workers, the same rows as on one", "[write_reflections][portable]") {
// A sweep large enough that the sort, the event sums and the row fill are all split into
// pieces, with events of several parts, events cut short by the sweep's ends, and parts that
// share (h,k,l) and image_number - the tie the sort has to break the same way every time.
auto experiment = TestExperiment();
experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
uint32_t seed = 12345;
const auto next = [&seed] { seed = seed * 1664525u + 1013904223u; return (seed >> 8) / 16777216.0f; };
std::vector<IntegrationOutcome> outcomes(200);
for (int f = 0; f < 200; ++f) {
for (int q = 0; q < 8000; ++q) {
if ((f + q) % 7 >= 3)
continue; // each reflection is seen on three consecutive images out of seven
Reflection r{};
r.h = q % 20 - 10; r.k = (q / 20) % 20 - 10; r.l = q / 400 + 1;
r.image_number = static_cast<float>(f);
r.I = 1000.0f * next();
r.sigma = 1.0f + 10.0f * next();
r.bkg = next();
r.var_bkg = next();
r.prescaling_corr = 1.0f + next();
r.partiality = 0.2f + 0.2f * next();
r.predicted_x = 1000.0f * next(); r.predicted_y = 1000.0f * next();
r.observed_x = NAN; r.observed_y = NAN;
outcomes[f].reflections.push_back(r);
if (q % 97 == 0) {
r.I = 1000.0f * next();
outcomes[f].reflections.push_back(r);
}
}
}
for (bool sum_partials : {true, false}) {
INFO("sum_partials " << sum_partials);
const gemmi::Mtz serial = UnmergedMtz(outcomes, TETRAGONAL_CELL, experiment, sum_partials, 1);
const gemmi::Mtz parallel = UnmergedMtz(outcomes, TETRAGONAL_CELL, experiment, sum_partials, 8);
REQUIRE(serial.nreflections > 100000);
REQUIRE(parallel.nreflections == serial.nreflections);
CHECK(std::memcmp(parallel.data.data(), serial.data.data(), serial.data.size() * sizeof(float)) == 0);
REQUIRE(parallel.batches.size() == serial.batches.size());
CHECK(parallel.batches.front().number == serial.batches.front().number);
CHECK(parallel.batches.back().number == serial.batches.back().number);
}
}