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.