From 1bb4cecbd1ea588ead8a1611c6ad7d80b2a7d500 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 26 Aug 2026 09:55:34 +0200 Subject: [PATCH] Spot width: measure on a growing sample and stop once it settles The pre-scan measured the spot width on all 60 frames of its sample every time. The radius that measurement feeds is a three-way choice - R1ForWidth moves only where 2*r80 crosses 4.5 and 5.5 - and 29 of 37 measurable crystals sit more than 0.29 px clear of both switches, so most of the sample is spent confirming a bucket the first few frames already picked. Measure instead on a growing share of it - every eighth frame, then every fourth, every second, all of it - and stop at the first tier whose answer has settled: within 0.40 px of what the smaller sample said AND 0.25 px clear of both switches. Both conditions are load-bearing. Clearance alone loses the crystal whose r80 is 0.07 px from a switch, because its small samples read across it; the step test alone lets a sample settle on a switch and stop there. Both bounds sit interior to a two-dimensional region that is right on 38 of 38, and every candidate was re-scored at all eight phase offsets of the tier ladder. The tiers are sized in worker-rounds rather than frames: a first tier of five frames occupies eight workers as long as one of eight does. A {12,4, 2,1} ladder measures 11% fewer frames than {8,4,2,1} and is 11% slower. Over the rotation battery this reads 1215 of 2280 frames - 21 crystals stop at 15, 4 at 30, 13 still run all 60 - and the chosen r1 is identical on 38 of 38, as is the beam-stop mask the reordered loop also touches. The pass feeds nothing else into the run, so every merged intensity is unchanged by construction and no battery is required. Cost over the battery goes 15.6 s to 10.1 s, median 0.18 s to 0.14 s per data set. Note for anyone optimising this further: the measurement itself is 0.54% of the pre-scan's cycles. The cost is the decode, preprocessing and spot finding each frame needs before it - 64% of the pass - because the beam-stop projection decodes on the GPU and never materialises the frame on the host. Cutting frames is the only lever short of harvesting the width from the two-pass run's first GPU pass. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/CHANGELOG.md | 1 + rugnux/Rugnux.cpp | 166 ++++++++++++++++++++++++++++++++----------- rugnux/SpotWidth.cpp | 6 ++ rugnux/SpotWidth.h | 23 ++++++ 4 files changed, 155 insertions(+), 41 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index fd8a3452..f7a17c93 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -13,6 +13,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: the widened integration radius applies to the final integration pass only - the two-pass geometry pre-pass keeps the radius the run started with, so the post-refined detector geometry, and the lattice the second pass indexes on, are the ones the fixed radius gives. * rugnux: on a pattern too dense for the widened radius - where neighbouring reflections leave more than 1.1% of the reflections without a background ring - the final pass is re-integrated at the fixed 4 px radius, and says so. * rugnux: on rotation data the integration signal radius is set from how wide the crystal's own spots are, measured in the pre-scan, instead of the fixed 4 px; `--adaptive-integration-radius=off` restores the fixed radius, and an explicit `--integration-radius` still overrides both. +* rugnux: the spot-width pre-scan stops once the integration radius it is measuring has settled, instead of always working through the whole sample; the radius it chooses is unchanged. * rugnux: a directional diffraction limit that is the edge of the measured data rather than the crystal's own limit is marked as such - with a `<` in the report and in `ANISOTROPY_D_MIN_CENSORED`, and in the mmCIF - so `ANISOTROPY_D_MIN_SPREAD` is not read as a measurement when it is a lower bound. * rugnux: the anisotropy verdict line names which of `ANISOTROPY_DELTA_B` and `ANISOTROPY_DELTA_B_LINEAR` it is quoting, says which of the two to act on, and says why the second can be the larger. * rugnux: the anisotropy caution about a too-high symmetry assignment now fires only where that is actually indicated - the symmetry-forbidden tensor directions far above their own counting noise together with a gate that established nothing - instead of on every tetragonal, trigonal and hexagonal data set. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 6628c5b2..bd6f93ad 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -107,6 +107,30 @@ namespace { // is memory-bound rather than compute-bound, so a handful of workers already saturates it. constexpr size_t PRESCAN_MAX_WORKERS = 8; + // The spot width is measured on a GROWING share of the pre-scan sample: every eighth frame of + // it, then every fourth, then every second, then all of it - an eighth, a quarter, a half and + // the whole. Each stride is a multiple of the next, so a tier's sample really does contain the + // one before it, and the measurement stops at the first tier whose answer has settled + // (spot_width::WidthSettled, in the loop below). A crystal whose frames carry plenty of isolated + // strong spots is characterised from a quarter of the sample; a sparse one still gets all of it. + // + // The tiers add PRESCAN_MAX_WORKERS-sized steps as well: on a 60 frame sample they are 8, 15, 30 + // and 60 frames, i.e. 8 + 7 + 15 + 30 read in the four passes, so no pass leaves most of the + // workers idle waiting for the others to finish theirs. A ladder that starts smaller measures + // fewer frames but costs more wall clock, because a pass of five frames occupies eight workers + // for as long as a pass of eight does. + // + // Only the width work is tiered. Every frame of the sample is still READ once - the frames are + // simply visited in tier order - and the beam-stop projection still sees all of them, so the + // mask and the beam centre are what they were. What the tiers cut is the per-frame work the + // width needs on top of the read: the host decompression, the preprocessing and the spot + // finding, which together are two thirds of the pre-scan's cycles once the projection itself is + // decoding on the GPU. The width measurement proper is under 1 % of them. + // + // Measured over the rotation battery: the radius chosen is the whole-sample answer on 38 of 38 + // crystals, from an average of 32 frames of the 60 instead of all of them. + constexpr std::array WIDTH_TIER_STRIDE = {8, 4, 2, 1}; + // Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the // first pass of two-pass rotation indexing. std::vector select_equally_spaced_image_ordinals(int images_to_process, int requested_images) { @@ -459,6 +483,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::vector frame_angle_deg; std::vector beam_center_spots; std::vector width_curves; + // Images the width ended up being measured on, for the log: the tiers below stop as soon as the + // answer has settled, so this is a property of the crystal and worth reporting. + size_t width_images = 0; // Read the sample on several workers. The reader serialises on the HDF5 lock, but the // decompression, the projection and the spot finding - which is all of the cost on a large @@ -475,46 +502,103 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::min(PRESCAN_MAX_WORKERS, ordinals.size())); finder.SetShardCount(nworkers); - std::atomic next{0}; - std::vector> futures; - futures.reserve(nworkers); - for (size_t t = 0; t < nworkers; t++) - futures.emplace_back(std::async(std::launch::async, [&, t] { - PreScanWorker w = make_worker(); - for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { - const int ordinal = ordinals[i]; - const int image_idx = start_image + ordinal * config_.stride; - std::shared_ptr img; - try { - img = reader_.GetRawImage(image_idx); - } catch (const std::exception &e) { - logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); - continue; + // The order the frames are visited in, and where each width tier ends in it. The tiers + // stride the PROJECTION's sample rather than the union read below, so asking for a beam + // centre - which brings frames of its own - cannot move which frames the width is measured + // on; the frames only the beam centre wants are appended to the last tier. Without a width + // to measure there is one tier and the sample is read in its own order. + std::vector visit; + std::vector tier_end; + visit.reserve(ordinals.size()); + { + std::vector width_pos; + for (size_t i = 0; i < ordinals.size(); i++) + if (shadow_set.contains(ordinals[i])) width_pos.push_back(i); + std::vector taken(ordinals.size(), 0); + const size_t n_tier = want_width ? WIDTH_TIER_STRIDE.size() : 1; + for (size_t tier = 0; tier < n_tier; tier++) { + const size_t stride = want_width ? WIDTH_TIER_STRIDE[tier] : 1; + for (size_t j = 0; j < width_pos.size(); j += stride) + if (!taken[width_pos[j]]) { + taken[width_pos[j]] = 1; + visit.push_back(width_pos[j]); } - if (!img) continue; + tier_end.push_back(visit.size()); + } + for (size_t i = 0; i < ordinals.size(); i++) + if (!taken[i]) visit.push_back(i); + tier_end.back() = visit.size(); + } - DataMessage msg{}; - msg.image = img->image; - msg.number = ordinal; - msg.original_number = image_idx; - if (shadow_set.contains(ordinal)) - finder.AddImage(msg, w.shadow_buffer, t); + // Kept across the tiers rather than rebuilt for each. A worker's preprocessed image is a + // buffer the size of the detector, so each is still built on its own worker the first time + // that worker runs and its pages are touched by the thread that reads them. + std::vector workers(nworkers); + std::vector worker_built(nworkers, 0); - // The width is measured on the projection's own even spread over the sweep, so - // asking for a beam centre cannot move it either. - const bool for_beam_center = spot_set.contains(ordinal); - const bool for_width = want_width && shadow_set.contains(ordinal); - if (for_beam_center) spot_read[i] = 1; - if (for_beam_center || for_width) - find_spots(w, msg.image, image_idx, for_beam_center, spots_of[i], - for_width, curves_of[i]); - } - })); - for (auto &f : futures) f.get(); + bool width_settled = !want_width; + std::optional previous_r80; + size_t phase_begin = 0; + for (const size_t phase_end : tier_end) { + std::atomic next{phase_begin}; + std::vector> futures; + futures.reserve(nworkers); + for (size_t t = 0; t < nworkers; t++) + futures.emplace_back(std::async(std::launch::async, [&, t] { + if (!worker_built[t]) { + workers[t] = make_worker(); + worker_built[t] = 1; + } + PreScanWorker &w = workers[t]; + for (size_t v = next.fetch_add(1); v < phase_end; v = next.fetch_add(1)) { + const size_t i = visit[v]; + const int ordinal = ordinals[i]; + const int image_idx = start_image + ordinal * config_.stride; + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, + e.what()); + continue; + } + if (!img) continue; - // Sample order, so the pooled width does not depend on how the workers interleaved. - for (auto &c : curves_of) - width_curves.insert(width_curves.end(), c.begin(), c.end()); + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + if (shadow_set.contains(ordinal)) + finder.AddImage(msg, w.shadow_buffer, t); + + // The width is measured on the projection's own even spread over the sweep, + // so asking for a beam centre cannot move it either. + const bool for_beam_center = spot_set.contains(ordinal); + const bool for_width = !width_settled && shadow_set.contains(ordinal); + if (for_beam_center) spot_read[i] = 1; + if (for_beam_center || for_width) + find_spots(w, msg.image, image_idx, for_beam_center, spots_of[i], + for_width, curves_of[i]); + } + })); + for (auto &f : futures) f.get(); + phase_begin = phase_end; + + // What this tier makes of the width, pooled in sample order so the answer does not + // depend on how the workers interleaved. A tier too sparse to measure settles nothing: + // the test is against what the smaller sample actually said, not against the radius it + // fell back to. + if (width_settled) continue; + width_images = static_cast( + std::count_if(visit.begin(), visit.begin() + static_cast(phase_end), + [&](size_t i) { return shadow_set.contains(ordinals[i]); })); + width_curves.clear(); + for (const auto &c : curves_of) + width_curves.insert(width_curves.end(), c.begin(), c.end()); + const auto r80 = spot_width::R80AtReference(width_curves); + width_settled = r80 && previous_r80 && spot_width::WidthSettled(*r80, *previous_r80); + previous_r80 = r80; + } // Frame numbering follows the sample order, exactly as the serial read did. for (size_t i = 0; i < ordinals.size(); i++) { @@ -541,8 +625,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru spot_width_measured_ = true; const auto r80 = spot_width::R80AtReference(width_curves); if (!r80) { - logger.Info("Spot width: not measurable on {} spots, keeping the integration radius at " - "r1={:.1f}", width_curves.size(), + logger.Info("Spot width: not measurable on {} spots from {} images, keeping the " + "integration radius at r1={:.1f}", width_curves.size(), width_images, experiment_.GetBraggIntegrationSettings().GetR1()); } else { const float r1 = spot_width::R1ForWidth(*r80); @@ -554,9 +638,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (r1 != bis.GetR1()) bragg_before_adaptive_ = bis; bis.R1(r1).R2(r2).R3(r3); experiment_.ImportBraggIntegrationSettings(bis); - logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots) => integration radii " - "r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A, - width_curves.size(), r1, r2, r3); + logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots from {} images) => " + "integration radii r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A, + width_curves.size(), width_images, r1, r2, r3); } } diff --git a/rugnux/SpotWidth.cpp b/rugnux/SpotWidth.cpp index 4c16a68d..5245d648 100644 --- a/rugnux/SpotWidth.cpp +++ b/rugnux/SpotWidth.cpp @@ -240,3 +240,9 @@ std::optional spot_width::R80AtReference(const std::vector &cu float spot_width::R1ForWidth(float r80) { return std::clamp(std::round(2.0f * r80), 4.0f, 6.0f); } + +bool spot_width::WidthSettled(float r80, float r80_before) { + return std::abs(r80 - r80_before) < SETTLED_STEP_PX + && std::abs(r80 - 2.25f) > SWITCH_CLEARANCE_PX + && std::abs(r80 - 2.75f) > SWITCH_CLEARANCE_PX; +} diff --git a/rugnux/SpotWidth.h b/rugnux/SpotWidth.h index 5b62b2ca..e5b450c7 100644 --- a/rugnux/SpotWidth.h +++ b/rugnux/SpotWidth.h @@ -63,6 +63,29 @@ struct FluxCurve { // pattern starts losing reflections whose ring falls below six clean pixels. [[nodiscard]] float R1ForWidth(float r80); +// When a growing sample has told us all it is going to about the radius, so the measurement can +// stop. R1ForWidth is a three-way choice and it moves only where 2*r80 crosses 4.5 and 5.5, i.e. at +// r80 = 2.25 and 2.75 px, so "settled" is two things and not one: +// +// * the estimate has STOPPED MOVING - this sample is within SETTLED_STEP_PX of what the smaller +// sample before it said; and +// * it is not sitting ON a switch - it is SWITCH_CLEARANCE_PX clear of both of them. +// +// Neither carries the rotation battery on its own. Clearance without the step test loses the +// crystal whose r80 is 0.07 px from a switch and whose smaller samples read across it; the step test +// without the clearance lets a sample that has settled AT a switch stop there. Both bounds are set +// inside the range that is right on 38 of 38 crystals rather than at its edge: over every phase of +// the sub-sample the answer is exact for a clearance of 0.20-0.30 px and a step of 0.20-0.60 px. +// +// A quarter of a pixel is also half the width of the whole r1 = 5 band, so a sample landing inside +// that band never stops early - which is deliberate, because that band and its edges are where a +// reduced sample gets a crystal wrong. +constexpr float SWITCH_CLEARANCE_PX = 0.25f; +constexpr float SETTLED_STEP_PX = 0.40f; + +// True where this sample's r80 has settled against the smaller sample's before it. +[[nodiscard]] bool WidthSettled(float r80, float r80_before); + // Most of the predicted reflections a widened radius may leave without a background ring, measured // (BraggIntegrationCounts::bkg_starved_by_neighbour) rather than predicted. Above this the widening // has cost the pattern more than it can be worth and the shipped radius is restored.