mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity

Build the merge's per-frame quantities once per frame, and count what a shell can hold only when it is reported

Browse Source 
Three passes over the ingested observations were doing more than they needed.

The smoothed-geometry pass rebuilt a CrystalLattice and its three reciprocal vectors for
every observation, each one a cross product and a cell volume, for a value that depends only
on which frame the observation came from. On a large sweep that is tens of millions of
constructions against a couple of thousand distinct answers.

The completeness column counts how many unique reflections a shell could hold. It is read off
a merge that gets written out, never off the ones the space-group search runs on the way
there - and those are the expensive ones to count, because the search merges in P1, where the
list is the whole hemisphere rather than an asymmetric unit of it.

The keep flags were written over the whole observation array as 1 and then immediately over
it again as 0 whenever a resolution limit is set, which the default low-resolution limit
always does. They are filled once now. The merged reflections also get their capacity up
front rather than doubling their way to it several times per pass.

Merged intensities are unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
jungfrau
2026-08-24 05:32:12 -04:00
co-authored by Claude Opus 5
parent 0d07141b5c
commit be75001833
1 changed files with 26 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+26 -6
View File
@@ -472,9 +472,10 @@ void RotationScaleMerge::BuildInRangeObservations(const std::vector<SortKey> &ke
const bool limited = d_min_limit.has_value() || d_max_limit.has_value();
std::vector<uint8_t> keep_run(n_run, 1);
std::vector<uint8_t> keep(n_obs, 1);
// Filled once, with the value the branch below actually wants: a resolution-limited run used to
// write the whole array as 1 and then immediately write it again as 0.
std::vector<uint8_t> keep(n_obs, limited ? 0 : 1);
if (limited) {
keep.assign(n_obs, 0);
ParallelChunks(n_run, nthreads, [&](int lo, int hi) {
for (int r = lo; r < hi; ++r) {
const float d = rawrun_d[r];
@@ -733,6 +734,16 @@ void RotationScaleMerge::SmoothGeometry() {
const Coord m3 = (m1 % m2).Normalize();
const float m2_S0 = m2 * S0, m3_S0 = m3 * S0;
// The reciprocal basis follows the frame's smoothed lattice, not the observation, so build it once
// per frame. Each of Astar/Bstar/Cstar is a cross product and a cell volume, and there are tens of
// millions of observations against a couple of thousand frames.
std::vector<std::array<Coord, 3> > frame_recip(n_frames);
for (int f = 0; f < n_frames; ++f) {
if (!ok[f]) continue;
const CrystalLattice sl(smoothed[f][0], smoothed[f][1], smoothed[f][2]);
frame_recip[f] = {sl.Astar(), sl.Bstar(), sl.Cstar()};
}
std::atomic<int64_t> changed{0};
ParallelChunks(static_cast<int>(partials.size()), nthreads, [&](int lo, int hi) {
int64_t local = 0;
@@ -740,9 +751,9 @@ void RotationScaleMerge::SmoothGeometry() {
auto &o = partials[i];
if (o.frame < 0 || o.frame >= n_frames || !ok[o.frame])
continue;
const CrystalLattice sl(smoothed[o.frame][0], smoothed[o.frame][1], smoothed[o.frame][2]);
const Coord p0 = sl.Astar() * static_cast<float>(o.h) + sl.Bstar() * static_cast<float>(o.k)
+ sl.Cstar() * static_cast<float>(o.l);
const std::array<Coord, 3> &r = frame_recip[o.frame];
const Coord p0 = r[0] * static_cast<float>(o.h) + r[1] * static_cast<float>(o.k)
+ r[2] * static_cast<float>(o.l);
const float phi_near = o.delta_phi * static_cast<float>(PI) / 180.0f;
const auto phi = SolvePhi(p0, S0, m1, m2, m3, m2_S0, m3_S0, phi_near);
if (!phi || !std::isfinite(*phi))
@@ -2857,6 +2868,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
auto run_merge = [&]() {
acc.assign(n_groups, Accum{});
result.merged.clear();
// One reflection per group at most, and a merge runs several times per pass. Without this the
// vector doubles its way up to a few hundred megabytes on a large cell, copying everything it
// already holds each time.
result.merged.reserve(n_groups);
merged_I.assign(n_groups, NAN);
reject_count = 0;
rejected_obs.assign(fulls.size(), 0);
@@ -3023,7 +3038,12 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
};
std::vector<ShellAcc> sa(n_shells);
std::vector<int> possible(n_shells, 0);
if (reference_cell)
// Completeness is a reported quantity: it is read off the merge a pass writes out, never off the
// ones the space-group search runs on the way there. Those are also the expensive ones to count -
// the search merges in P1, where the list of possible reflections is the whole hemisphere rather
// than an asymmetric unit of it - so enumerating it there is the largest single piece of work in
// the merge that nothing goes on to read.
if (reference_cell && !for_search)
PossiblePerShell(x.GetSpaceGroupNumber().value_or(1), *reference_cell, d_min_pad, d_max_pad,
shells, merge_friedel, possible);
for (int s = 0; s < n_shells; ++s) sa[s].possible = possible[s];