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

Stop scaling and merging what the resolution range excludes

Browse Source 
A crystal integrated to the detector corner but merged well short of it carries
observations through the whole merge that the merge then discards. On the heaviest
dataset in the rotation test set that is 63.3 M partials of which 6.4 M are ever
used: the other nine tenths are sorted, uploaded, scaled, combined, error-modelled
and post-refined before anything looks at their resolution. Ingest copied every one
of them unconditionally, and the d_min limit was first applied far downstream, in
the ASU grouping.

They are now dropped at ingest, immediately after the one big sort:

- WHOLE raw-hkl runs are dropped, on the same rawrun_d the ASU grouping already
  tests. A per-observation test is not equivalent - a run is in or out today by one
  member's d - and using a different rule here would put the two out of step.
- The drop happens AFTER the flux meter, which takes each frame's mean background
  over every reflection on it, and after the sort, so neither changes.
- The compaction runs in index order, so a frame's observations stay contiguous and
  keep their order, and every per-frame sum keeps its sequence of roundings.

Post-refinement reads the integration outcomes rather than the merge arrays, so it
still sees every reflection. That is the point: --integration-high-resolution buys
the same time by never integrating the reflections, and pays for it in the per-frame
geometry, which wants them.

Two more passes over the observation array go with it. The incident-flux divide was
15.8 % of all user cycles to read one int and divide one float across 5 GB; the
per-frame mean it needs is now accumulated by the ingest fill loop - one frame, one
thread, same order, so bit-exact - and the divide rides on the finiteness pass that
already touches that field. And the geometry post-refinement is fitted on a bounded
sample of partials, selected by a hash of the raw hkl so whole rocking events are
kept or dropped together and the sweep and the detector are thinned uniformly.

The sample size is 8 M and the reason it is not smaller is measured. Over a 126x
thinning the fitted rotation scale is flat to 2e-5 and the beam centre moves 0.03 px,
but the CELL scale breaks between 8 M and 4 M: the axis step keeps the 20 000
strongest events, so once the pool approaches that size it starts fitting weaker ones
and the second pass's cell shifts by ~0.1 %.

Measured on the heaviest crystal, three A/B pairs with the order alternated:
68.4 s -> 37.3 s wall, 530 s -> 221 s of CPU. Whole battery 8m07s -> 7m15s, space
group 21/24 with the same three disagreements as before, no failures. Bit-identical
is not available on the GPU path - the resident reductions and the fulls emit order
depend on array length - so what is shown is that every difference sits inside the
spread the unmodified binary has against itself between two runs.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
jungfrau
2026-08-17 18:04:55 -04:00
co-authored by Claude Opus 5
parent 20f1827550
commit 2308bbad8c
5 changed files with 178 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CHANGELOG.md
+1
View File
@@ -1,6 +1,7 @@
# Changelog
## 1.0.0
### Unreleased
* **rugnux: scaling and merging are much faster on crystals that integrate far beyond the resolution they merge at.** Observations outside the scaling resolution range are dropped as they are ingested instead of being scaled, combined and post-refined first, and the geometry post-refinement fits on a bounded sample of them. On a large-cell dataset merging at less than half the resolution its detector reaches, a run drops from 68 s to 37 s; merged statistics are unchanged.
* rugnux: the detector-frame modulation correction is fitted on a grid spanning the detector rather than the reflections that happen to be present, so whether it is applied no longer depends on how far integration reached. On a crystal where it was being refused, merged R_meas improves by 4 percentage points.
* rugnux: the beam-stop pre-scan reads and accumulates its frames on several threads instead of one. On a 16M-pixel rotation dataset it drops from 7.9 s to 4.8 s; the shadow it finds is unchanged.
image_analysis/geom_refinement/PostRefine.cpp
+28 -4
View File
@@ -109,18 +109,40 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
const Coord ax = axis.GetAxis();
const Coord Astar = reference_latt.Astar(), Bstar = reference_latt.Bstar(), Cstar = reference_latt.Cstar();
// This refinement fits five numbers - a uniform cell scale, the rotation-axis direction, the
// detector distance and the beam centre - and the two solver steps below already cap themselves at
// 20 000 rocking events and 20 000 spot positions. Gathering, sorting and event-splitting every
// partial of a large dataset to feed those caps is work nothing reads. So thin the input by RAW
// HKL: a whole rocking event is then kept or dropped together, so no phi_obs centroid is biased;
// every part of the sweep and every part of the detector face is thinned by the same factor, which
// is what the distance/beam step and the rotation-scale scan are sensitive to; and the subset
// depends on the indices alone, so it is the same one on every run.
const size_t nthreads = std::max(1, settings.num_threads);
const int n_out = static_cast<int>(outcomes.size());
size_t n_integrated = 0;
for (const auto &o : outcomes) n_integrated += o.reflections.size();
const size_t stride = std::max<size_t>(1, n_integrated / std::max<size_t>(1, settings.max_partials));
// Mixed differently from the cross-validation split further down (is_val), so which hkls are
// sampled and which fold they land in are independent.
const auto in_subset = [stride](int h, int k, int l) {
if (stride == 1) return true;
unsigned u = static_cast<unsigned>(h) * 2246822519u + static_cast<unsigned>(k) * 3266489917u
+ static_cast<unsigned>(l) * 668265263u;
u ^= u >> 15; u *= 2654435761u; u ^= u >> 16;
return u % stride == 0;
};
// Count first, then fill. Growing one vector by push_back over tens of millions of
// reflections copies the whole thing every time it doubles - several gigabytes of pure
// copying - and the counts are cheap to take. Each outcome then owns a slice, so the fill
// runs on all threads and lands in the order the serial loop produced.
const size_t nthreads = std::max(1, settings.num_threads);
const int n_out = static_cast<int>(outcomes.size());
std::vector<size_t> pts_offset(n_out + 1, 0);
ParallelChunks(n_out, ThreadsForWork(static_cast<size_t>(n_out), nthreads), [&](int lo, int hi) {
for (int o = lo; o < hi; o++) {
size_t keep = 0;
for (const auto &r : outcomes[o].reflections)
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f)
if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f
&& in_subset(r.h, r.k, r.l))
keep++;
pts_offset[o + 1] = keep;
}
@@ -134,6 +156,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
size_t at = pts_offset[o];
for (const auto &r : outcomes[o].reflections) {
if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue;
if (!in_subset(r.h, r.k, r.l)) continue;
const double mid_deg = axis.GetAngle_deg(r.image_number) + wedge_half;
const double ox = std::isfinite(r.observed_x) ? r.observed_x : NAN;
const double oy = std::isfinite(r.observed_y) ? r.observed_y : NAN;
@@ -142,7 +165,8 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
}
}
});
logger.Info("Post-refine: {} partials gathered", pts.size());
logger.Info("Post-refine: {} partials gathered (1 raw hkl in {} of {} integrated)", pts.size(),
stride, n_integrated);
if (pts.size() < static_cast<size_t>(settings.min_events)) return result;
// Bucket by h, then sort the buckets. h is the leading key, so the sorted array is the
image_analysis/geom_refinement/PostRefine.h
+4
View File
@@ -51,6 +51,10 @@ struct PostRefineSettings {
// separate cross-validated steps. The only supported refinement mode.
double excitation_weight = 1.0; // weight of the phi/excitation residual vs the positional one
int min_events = 50;
// Cap on how many partials are gathered, reached by thinning whole raw hkls out of the set (see the
// .cpp). Five geometry parameters do not need tens of millions of observations, and the two solver
// steps already cap themselves an order of magnitude below this.
size_t max_partials = 8000000;
int num_threads = 1;
};
image_analysis/scale_merge/RotationScaleMerge.cpp
+131 -25
View File
@@ -184,8 +184,8 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
mosaicity_deg = s.GetDefaultMosaicity();
}
// Divide the incident flux out of the scale before anything is fitted, using the per-frame mean of the
// integrator's background estimate as the flux meter.
// Measure the incident flux per frame from the per-frame mean of the integrator's background estimate,
// so it can be divided out of the scale before anything is fitted.
//
// The beam is not constant. On a large-bandwidth source it oscillates by ~10% with a period of a few
// frames, and the per-frame scale G cannot follow that: G is smoothed over --smooth-g degrees (5 by
@@ -200,21 +200,13 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
// changing with the goniometer angle also move the background. Those are all SLOW, and what they put in
// here the fitted G takes back out, because smoothing leaves G its low frequencies. Only the fast part -
// the part no smoothed per-frame scale can represent anyway - is taken on trust.
void RotationScaleMerge::DivideOutIncidentFlux() {
//
// Neither half of this costs a pass over the observations of its own: the per-frame mean background is
// accumulated by the ingest fill loop, which already reads every reflection, and the division into rlp
// rides on the finiteness pass that follows it.
void RotationScaleMerge::MeasureIncidentFlux(const std::vector<double> &mean_bkg) {
frame_flux.assign(n_frames, 1.0);
// A frame at a time, not an observation at a time: each frame's mean is a running sum over its
// own contiguous range, so giving a thread whole frames keeps every sum in one thread and in the
// order it had. Splitting by observation would cut a frame across two threads and the sum would
// have to be recombined, which is a different sequence of roundings.
std::vector<double> mean_bkg(n_frames, NAN);
ParallelFor(n_frames, nthreads, [&](int f) {
double sum = 0.0;
int n = 0;
for (int i = frame_start[f]; i < frame_start[f] + frame_count[f]; ++i)
if (std::isfinite(partials[i].bkg) && partials[i].bkg > 0.0f) { sum += partials[i].bkg; ++n; }
if (n > 0) mean_bkg[f] = sum / n;
});
// Collected in frame order, serially: it is one value per frame, and the median only depends on
// the set, but keeping the order removes the question.
std::vector<double> finite;
@@ -231,13 +223,6 @@ void RotationScaleMerge::DivideOutIncidentFlux() {
const double median = finite[mid];
for (int f = 0; f < n_frames; ++f)
if (std::isfinite(mean_bkg[f])) frame_flux[f] = mean_bkg[f] / median;
ParallelChunks(static_cast<int>(partials.size()), nthreads, [&](int lo, int hi) {
for (int i = lo; i < hi; ++i) {
auto &o = partials[i];
o.rlp = static_cast<float>(o.rlp / frame_flux[o.frame]);
}
});
}
void RotationScaleMerge::Ingest() {
@@ -269,14 +254,21 @@ void RotationScaleMerge::Ingest() {
acc += frame_count[o];
}
}
// The per-frame mean background - the incident-flux meter, see MeasureIncidentFlux - is accumulated
// here rather than in a pass of its own. One frame is still summed by one thread and in its own
// order, which is what keeps the mean bit-exact whatever the thread count.
std::vector<double> mean_bkg(n_frames, NAN);
partials.resize(total);
ParallelFor(n_frames, ThreadsForWork(total, nthreads), [&](int o) {
if (reference_cell) {
const auto cell = partials_out[o].latt.GetUnitCell();
frame_cell_ok[o] = cell.is_close(*reference_cell, dist_tol, ang_tol) ? 1 : 0;
}
double bkg_sum = 0.0;
int bkg_n = 0;
int32_t at = frame_start[o];
for (const auto &r : partials_out[o].reflections) {
if (std::isfinite(r.bkg) && r.bkg > 0.0f) { bkg_sum += r.bkg; ++bkg_n; }
Obs &obs = partials[at++];
obs.h = r.h; obs.k = r.k; obs.l = r.l;
obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; obs.rlp = r.rlp;
@@ -288,9 +280,10 @@ void RotationScaleMerge::Ingest() {
obs.corr = r.image_scale_corr;
obs.group = -1;
}
if (bkg_n > 0) mean_bkg[o] = bkg_sum / bkg_n;
});
DivideOutIncidentFlux();
MeasureIncidentFlux(mean_bkg);
// Per-obs AcceptReflection finiteness (immutable) - lets ComputeAsuGroups stamp the ASU-group id per
// obs from a flat 1-byte array instead of re-reading the fat Obs struct for every space group.
@@ -316,7 +309,11 @@ void RotationScaleMerge::Ingest() {
int lmin = INT_HI, lmax = INT_LO;
const int lo = t * chunk, hi = std::min(n_obs, lo + chunk);
for (int i = lo; i < hi; ++i) {
const auto &o = partials[i];
auto &o = partials[i];
// Divide the incident flux out of rlp here: this pass already has the observation
// in hand, and rlp is the numerator of corr, so the fits below see only what is
// left. Before the finiteness test, which is where it has always been.
o.rlp = static_cast<float>(o.rlp / frame_flux[o.frame]);
finite_ok[i] = (std::isfinite(o.I) && std::isfinite(o.rlp) && o.rlp != 0.0f
&& std::isfinite(o.sigma) && o.sigma > 0.0f) ? 1 : 0;
h_of[i] = o.h;
@@ -464,9 +461,10 @@ void RotationScaleMerge::Ingest() {
}
});
}
DropOutOfRangeObservations();
rawrun_group.assign(rawrun_start.size(), -1);
logger.Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)",
total, n_frames, rawrun_start.size());
partials.size(), n_frames, rawrun_start.size());
SmoothMosaicityAndPartiality();
@@ -506,6 +504,114 @@ void RotationScaleMerge::Ingest() {
#endif
}
// Compact `partials` (and the perm / raw-hkl tables that index it) down to the observations whose
// resolution is inside the requested range. See the header for why this cannot change any result: a raw
// hkl outside the range is group -1 in every pass, and dropping WHOLE runs on the same per-hkl resolution
// ComputeAsuGroups tests leaves rawrun_d, the surviving observations and their order untouched.
//
// It has to come after the incident flux has been measured (the flux meter is the frame's mean background
// over ALL its reflections) and after the sort (which is what defines a raw-hkl run in the first place).
void RotationScaleMerge::DropOutOfRangeObservations() {
if (!d_min_limit && !d_max_limit)
return;
const int n_run = static_cast<int>(rawrun_start.size());
const int n_obs = static_cast<int>(partials.size());
std::vector<uint8_t> keep_run(n_run, 0);
std::vector<uint8_t> keep(n_obs, 0);
ParallelChunks(n_run, nthreads, [&](int lo, int hi) {
for (int r = lo; r < hi; ++r) {
const float d = rawrun_d[r];
if (!std::isfinite(d) || d <= 0.0f) continue;
if (d_min_limit && d < *d_min_limit) continue;
if (d_max_limit && d > *d_max_limit) continue;
keep_run[r] = 1;
const int p0 = rawrun_start[r], p1 = p0 + rawrun_count[r];
for (int p = p0; p < p1; ++p) keep[perm[p]] = 1; // distinct runs own disjoint observations
}
});
// New index of every kept observation. Compacting in index order keeps each frame's observations
// contiguous and in the order they had, so the per-frame ranges still describe them and every
// per-frame sum keeps its sequence of roundings.
const int nt = static_cast<int>(std::clamp<size_t>(nthreads, 1, std::max(1, n_obs)));
const int chunk = (n_obs + nt - 1) / nt;
std::vector<int32_t> base(nt + 1, 0);
ParallelChunks(nt, nthreads, [&](int tlo, int thi) {
for (int t = tlo; t < thi; ++t) {
int32_t c = 0;
for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) c += keep[i];
base[t + 1] = c;
}
});
for (int t = 0; t < nt; ++t) base[t + 1] += base[t];
const int32_t n_keep = base[nt];
if (n_keep == n_obs)
return;
std::vector<int32_t> new_idx(n_obs);
std::vector<Obs> kept_obs(n_keep);
std::vector<uint8_t> kept_finite(n_keep);
ParallelChunks(nt, nthreads, [&](int tlo, int thi) {
for (int t = tlo; t < thi; ++t) {
int32_t at = base[t];
for (int i = t * chunk; i < std::min(n_obs, (t + 1) * chunk); ++i) {
if (!keep[i]) continue;
new_idx[i] = at;
kept_obs[at] = partials[i];
kept_finite[at] = finite_ok[i];
++at;
}
}
});
partials.swap(kept_obs);
finite_ok.swap(kept_finite);
kept_obs.clear(); kept_obs.shrink_to_fit();
frame_start.assign(n_frames, 0);
frame_count.assign(n_frames, 0);
for (int i = 0; i < n_keep; ) {
const int f = partials[i].frame;
int j = i;
while (j < n_keep && partials[j].frame == f) ++j;
frame_start[f] = i;
frame_count[f] = j - i;
i = j;
}
// The surviving runs, in the same order, each owning the same perm entries remapped to the new
// observation indices.
int n_keep_run = 0;
std::vector<int32_t> dst_run(n_run, -1), dst_start(n_run, 0);
{
int32_t at = 0;
for (int r = 0; r < n_run; ++r)
if (keep_run[r]) { dst_run[r] = n_keep_run++; dst_start[r] = at; at += rawrun_count[r]; }
}
std::vector<int32_t> new_perm(n_keep), ns(n_keep_run), nc(n_keep_run),
nh(n_keep_run), nk(n_keep_run), nl(n_keep_run);
std::vector<float> nd(n_keep_run);
ParallelChunks(n_run, nthreads, [&](int lo, int hi) {
for (int r = lo; r < hi; ++r) {
const int rr = dst_run[r];
if (rr < 0) continue;
const int32_t at = dst_start[r];
ns[rr] = at; nc[rr] = rawrun_count[r];
nh[rr] = rawrun_h[r]; nk[rr] = rawrun_k[r]; nl[rr] = rawrun_l[r]; nd[rr] = rawrun_d[r];
for (int32_t p = 0; p < rawrun_count[r]; ++p)
new_perm[at + p] = new_idx[perm[rawrun_start[r] + p]];
}
});
perm.swap(new_perm);
rawrun_start.swap(ns); rawrun_count.swap(nc);
rawrun_h.swap(nh); rawrun_k.swap(nk); rawrun_l.swap(nl); rawrun_d.swap(nd);
logger.Info("RotationScaleMerge: dropped {} of {} observations ({} of {} distinct hkl) outside the "
"requested resolution range - nothing downstream could have merged them",
n_obs - n_keep, n_obs, n_run - n_keep_run, n_run);
}
namespace {
// Rotate v about a unit axis by angle (Rodrigues).
Coord RotateAbout(const Coord &v, const Coord &axis, double angle) {
image_analysis/scale_merge/RotationScaleMerge.h
+14 -4
View File
@@ -201,10 +201,20 @@ private:
#endif
// --- helpers (each a flat pass; see the .cpp) ---
// Measure the incident flux per frame from the mean background under its reflections and fold it into
// rlp, so corr = rlp / (partiality * G) divides it out and G fits only the residual. See the .cpp for
// why a background is a usable flux meter and what it costs when it is not.
void DivideOutIncidentFlux();
// Turn the per-frame mean background under the reflections (accumulated by the ingest fill loop) into
// the per-frame incident flux, which the finiteness pass then folds into rlp so that
// corr = rlp / (partiality * G) divides it out and G fits only the residual. See the .cpp for why a
// background is a usable flux meter and what it costs when it is not.
void MeasureIncidentFlux(const std::vector<double> &mean_bkg);
// Drop the observations whose resolution can never be in range: --scaling-high/low-resolution are the
// coarsest limits any Run() uses (the space-group search only ever RAISES d_min), and an out-of-range
// raw hkl gets group -1 in every pass, which keeps it out of the scaling reference, the per-frame fit,
// the combine, the merge and the error model alike. On a crystal that integrates to the detector
// corner and merges well short of it that is most of the array. Whole raw-hkl RUNS are dropped, on the
// same per-hkl resolution ComputeAsuGroups tests, so what survives - and the order of every sum formed
// over it - is exactly what it was. No-op when no manual limit was given.
void DropOutOfRangeObservations();
// Compute the dense ASU-group id for the current space group by grouping the (pre-sorted) raw-hkl
// runs by their ASU key - one gemmi ASU reduction per distinct raw hkl, not per observation. Fills