diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 01b8c789..aad3bf52 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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. diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index e00bbeac..21948c46 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -109,18 +109,40 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(outcomes.size()); + size_t n_integrated = 0; + for (const auto &o : outcomes) n_integrated += o.reflections.size(); + const size_t stride = std::max(1, n_integrated / std::max(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(h) * 2246822519u + static_cast(k) * 3266489917u + + static_cast(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(outcomes.size()); std::vector pts_offset(n_out + 1, 0); ParallelChunks(n_out, ThreadsForWork(static_cast(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(settings.min_events)) return result; // Bucket by h, then sort the buckets. h is the leading key, so the sorted array is the diff --git a/image_analysis/geom_refinement/PostRefine.h b/image_analysis/geom_refinement/PostRefine.h index 09799d8f..6f2024c3 100644 --- a/image_analysis/geom_refinement/PostRefine.h +++ b/image_analysis/geom_refinement/PostRefine.h @@ -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; }; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 2687ec2f..ee8ee46a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -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 &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 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 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(partials.size()), nthreads, [&](int lo, int hi) { - for (int i = lo; i < hi; ++i) { - auto &o = partials[i]; - o.rlp = static_cast(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 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(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(rawrun_start.size()); + const int n_obs = static_cast(partials.size()); + + std::vector keep_run(n_run, 0); + std::vector 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(std::clamp(nthreads, 1, std::max(1, n_obs))); + const int chunk = (n_obs + nt - 1) / nt; + std::vector 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 new_idx(n_obs); + std::vector kept_obs(n_keep); + std::vector 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 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 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 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) { diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 5801f254..20f9aa09 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -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 &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