From 97cac203a60319074830aa8e8e84138fe019ccbf Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 19:24:40 +0200 Subject: [PATCH] rugnux: split FirstPassRotationIndexing into named steps Pure refactor, no behaviour change. The 2460-line function with ~25 lambdas becomes Rugnux::FirstPassRun (RugnuxFirstPass.h): the state the steps share (spot cache, validation frames, schemes, committed spot-list state, the speculative passes) as members, and one member function per step taking and returning the standing lattice (FirstPass). FirstPassRotationIndexing is now the ordered list of decisions: PrefetchCertainFrames, IndexSchemes (+ centre/short-axis speculation), AxisSignRescue, BeamCentreCheck (IndexAtMeasuredCentre, CentreWhenFileFails, CentreByPooledSpots, CentreWhenBothIndex, CentresDisagree), LongAxisRescue, AddEndWedge, ReleaseSpotEngines, RequireLattice, ShortAxisHypothesis, SubLatticeHypothesis, Pass1Fallback, last-resort FirstPassLadder, MeasureVerdict, RejectChanceLattice, AdoptLattice, StartLeftoverCensus, ReportTiltWalk, MeasureSpotBudget, Pass1SupercellGuard. Files: RugnuxFirstPass.cpp (order + shared machinery: spots, scoring, schemes, ladder), RugnuxFirstPassRescues.cpp, RugnuxFirstPassAccept.cpp. Bodies are moved verbatim (Rugnux members reached through rx.); every log line is unchanged. Literal duplicates folded: the two MAX_SPOT_KEYS constants, and the feed-every-scheme / run-every-indexer loops of the long-axis rescue and the short-axis fallback (now FeedSchemes/RunIndexing, the same calls in the same order). The leftover census, which outlives the pass, captures what it needs explicitly instead of [&]. Verified against the rc175 0b2dde214 binary on 16 sets (7kcn, 8v4o, 5nw5, 9rci, 6z9g, 9min, nothing_1/2, 7atg, myob/cytc_x10sa, 8a1a, kdp, insu low ISa, 6p8j, 9qw8) and -N 8 on 7kcn/8v4o/9rci: p.mtz md5-identical, logs identical as multisets apart from timings. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/CMakeLists.txt | 3 + rugnux/Rugnux.h | 2 + rugnux/RugnuxFirstPass.cpp | 3369 +++++++++-------------------- rugnux/RugnuxFirstPass.h | 265 +++ rugnux/RugnuxFirstPassAccept.cpp | 500 +++++ rugnux/RugnuxFirstPassRescues.cpp | 1045 +++++++++ 6 files changed, 2808 insertions(+), 2376 deletions(-) create mode 100644 rugnux/RugnuxFirstPass.h create mode 100644 rugnux/RugnuxFirstPassAccept.cpp create mode 100644 rugnux/RugnuxFirstPassRescues.cpp diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 6b437d83b..d21fa2c66 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -9,7 +9,10 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxPreScan.cpp RugnuxPasses.cpp RugnuxImagePass.cpp + RugnuxFirstPass.h RugnuxFirstPass.cpp + RugnuxFirstPassRescues.cpp + RugnuxFirstPassAccept.cpp RugnuxScaleMerge.cpp RugnuxDefaults.cpp RugnuxDefaults.h diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4df80b6b8..9e48b3bfe 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -892,6 +892,8 @@ class Rugnux { // The first pass of the two-pass rotation indexing: index the sweep and hand the lattice to the // indexer. True where the pass ends there - an indexing probe, whose result is its validation score. bool FirstPassRotationIndexing(PipelineLocals &p); + // Its steps and the state they share (RugnuxFirstPass.h). + struct FirstPassRun; void ImagePass(PipelineLocals &p); // Scaling and merging, the space-group search, the analyses of the merge and the merged output. // True where the pass ends early without merging: a post-refinement probe, or a starved pass. diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index f6664d617..dade24cc5 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -84,6 +84,7 @@ #include #include #include "RugnuxPipeline.h" +#include "RugnuxFirstPass.h" using namespace rugnux_internal; @@ -333,20 +334,13 @@ namespace { } } +// The first pass of a rotation run, as the ordered list of its decisions: index the first part of the +// sweep on two sampling schemes, put the standing lattice through the rescues and the hypotheses that +// may replace it, accept or refuse it on the pooled validation spots, and leave it forced on the +// indexer. Each step is a member of FirstPassRun (RugnuxFirstPass.h); the comment at each says why it +// sits where it does - the order is part of the decision. bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { auto &logger = p.logger; - auto &result = p.result; - auto &observer = p.observer; - auto &geometry_prepass = p.geometry_prepass; - auto &start_image = p.start_image; - // Every frame the first pass reads is in the first part of the sweep (Rugnux::PrepassEnd). - const int early_images = p.prepass_images; - auto &dataset = p.dataset; - auto &mapping = p.mapping; - auto &indexer_pool = p.indexer_pool; - auto &indexer = p.indexer; - auto &indexing_bound_A = p.indexing_bound_A; - auto &census_future = p.census_future; // An indexing probe over inputs a first pass has already been run on returned above, except // under verify_first_pass_memo, which runs the pass and compares. const std::vector first_pass_key = FirstPassInputKey(); @@ -356,929 +350,907 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { if (indexing_probe_only_ && memo != first_pass_memo_.end()) memo_hit = memo->evidence; - if (observer) - observer->OnPhase("Rotation indexing (first pass)"); - Step(!geometry_prepass && prepass_detector_geometry_ ? "Indexing at the refined geometry" : "Indexing"); + if (p.observer) + p.observer->OnPhase("Rotation indexing (first pass)"); + Step(!p.geometry_prepass && prepass_detector_geometry_ ? "Indexing at the refined geometry" : "Indexing"); - // Mid-exposure goniometer angle for an image ordinal (matches the per-image path). - auto rot_angle = [&](int ordinal) -> std::optional { - if (const auto g = experiment_.GetGoniometer()) - return g->GetAngle_deg(static_cast(ordinal)) + g->GetWedge_deg() / 2.0f; - return std::nullopt; - }; + FirstPassRun run(*this, p); + using FirstPass = FirstPassRun::FirstPass; + const int n_validation = static_cast(run.validation.size()); - // Spots for a first-pass image ordinal, cached (a frame is reused across schemes and the - // validation set; re-finding is expensive on the --redo-rotation-spots path). The workers below - // hand their spots back and this thread puts them in, so the map itself is still only ever - // touched from here and needs no locking. - std::map> spot_cache; - // The azimuthal mapping the first pass finds its spots against. It is `mapping` until a - // beam-centre hypothesis moves the centre, and the hypothesis's own mapping while it stands - - // the mapping is where the resolution of every pixel and the finder's radial background bins - // come from, and both move with the centre. - const AzimuthalIntegrationMapping *spot_mapping = &mapping; - // The spots of one frame, found with the caller's own analysis engine. A frame's spots are a - // pure function of that frame and the settings - the engine carries nothing from one image to - // the next, which is why the main image loop can give one to each of its workers - so it does - // not matter which engine, or which thread, finds them. - const auto find_spots = [&](MXAnalysisWithoutFPGA &analysis, AzimuthalIntegrationProfile &profile, - JFJochReaderRawImage &img, int ordinal) { - return FindFirstPassSpots(reader_, config_.spot_finding, *dataset, analysis, profile, img, - ordinal, start_image + ordinal * config_.stride, logger); - }; - - // Find the spots of every ordinal in the list that is not cached yet, on several workers, and - // put them in the cache. - // - // The cache ends up holding exactly what a serial loop would have put there: the results are - // collected into a vector indexed by position and inserted afterwards, so neither the spots of - // a frame nor the set of frames in the cache depends on how the workers interleaved. - // - // Worker t keeps its analysis engine for the whole first pass. Building one is a CUDA stream - // and a full set of device buffers - measured at 0.13 s against 4 ms for a frame - and the - // "wedge" scheme only discovers how far into the sweep it has to read as the frames arrive, - // so this is called several times and rebuilding them each time cost a quarter of the pass. - // The engine belongs to a card, not to a thread, so a later call only has to make that card - // current again; workers take a card by index so an engine always meets the one it was - // built on. - const int spot_workers = get_gpu_count() > 0 - ? std::min(config_.nthreads, std::max(8, 4 * get_gpu_count())) - : config_.nthreads; - struct SpotEngine { - // Before the analysis engine, which page-locks these bytes for its uploads and - // unregisters them when it is destroyed: it must not outlive the buffer it was given. - JFJochReaderRawImage raw_image; - std::unique_ptr analysis; - std::unique_ptr profile; - }; - std::vector engines(std::max(spot_workers, 1)); - const auto use_gpu_of_worker = [](size_t t) { - if (const int32_t gpus = get_gpu_count(); gpus > 0) - pin_gpu(static_cast(t) % gpus); - }; - int prefetch_batches = 0; - double prefetch_time_s = 0.0; - const auto prefetch_spots = [&](const std::vector &ordinals) { - // What an earlier pass already found under the same key (first_pass_spots_) is taken from - // there, at this pass's angles. - const std::vector key = SpotFindingKey(*spot_mapping); - std::map> kept; - { - std::lock_guard lock(first_pass_spots_->m); - for (const auto &[k, spots] : first_pass_spots_->entries) - if (k == key) - for (const int ordinal : ordinals) - if (const auto it = spots.find(ordinal); it != spots.end()) - kept.emplace(ordinal, it->second); - } - const DiffractionGeometry spot_geometry = experiment_.GetDiffractionGeometry(); - for (auto &[ordinal, spots] : kept) - StampSpotAngles(spots, spot_geometry); - std::vector wanted; - for (const int ordinal : ordinals) - if (!spot_cache.contains(ordinal)) { - if (const auto it = kept.find(ordinal); it != kept.end() && !config_.verify_first_pass_memo) - spot_cache.emplace(ordinal, it->second); - else - wanted.push_back(ordinal); - } - if (wanted.empty()) - return; - const auto prefetch_start = std::chrono::steady_clock::now(); - ++prefetch_batches; - - std::vector> found(wanted.size()); - std::atomic next{0}; - std::vector> futures; - const size_t nworkers = std::min(engines.size(), wanted.size()); - futures.reserve(nworkers); - for (size_t t = 0; t < nworkers; t++) - futures.emplace_back(std::async(std::launch::async, [&, t] { - use_gpu_of_worker(t); // must precede engine construction - SpotEngine &e = engines[t]; - if (!e.analysis) { - e.analysis = std::make_unique( - experiment_, *spot_mapping, pixel_mask_, *indexer, - /*enable_fused_adaptive_gpu=*/true); - e.profile = std::make_unique(*spot_mapping); - } - for (size_t i = next.fetch_add(1); i < wanted.size(); i = next.fetch_add(1)) - found[i] = find_spots(*e.analysis, *e.profile, e.raw_image, wanted[i]); - })); - for (auto &f : futures) - f.get(); - - { - std::lock_guard lock(first_pass_spots_->m); - auto &entries = first_pass_spots_->entries; - auto entry = std::find_if(entries.begin(), entries.end(), - [&](const auto &e) { return e.first == key; }); - if (entry == entries.end()) { - // A handful of keys covers a run; the oldest go first. - constexpr size_t MAX_SPOT_KEYS = 8; - if (entries.size() == MAX_SPOT_KEYS) - entries.erase(entries.begin()); - entries.emplace_back(key, std::map>()); - entry = entries.end() - 1; - } - for (size_t i = 0; i < wanted.size(); i++) { - if (const auto it = kept.find(wanted[i]); it != kept.end() && !SameSpots(it->second, found[i])) - throw MemoMismatch(fmt::format("First-pass spot memo differs from a recomputation " - "for image {}", wanted[i])); - entry->second.emplace(wanted[i], found[i]); - } - } - for (size_t i = 0; i < wanted.size(); i++) - spot_cache.emplace(wanted[i], std::move(found[i])); - prefetch_time_s += std::chrono::duration( - std::chrono::steady_clock::now() - prefetch_start).count(); - }; - - // Release the engines on workers of their own, so a stream and its device buffers are freed - // with the card they were made on current. - const auto release_spot_engines = [&] { - std::vector> futures; - futures.reserve(engines.size()); - for (size_t t = 0; t < engines.size(); t++) - futures.emplace_back(std::async(std::launch::async, [&, t] { - use_gpu_of_worker(t); - engines[t] = SpotEngine{}; - })); - for (auto &f : futures) - f.get(); - }; - - // Evaluating a beam-centre hypothesis means re-finding the spots at that centre. The spot - // POSITIONS do not move with the centre - which is why the cache is keyed by image alone - but - // the spot LIST does: every spot's d, the ice-ring flag, the resolution mask the finder - // applies, the high-resolution gap cut and the strongest-N budget are all radial, and the - // adaptive finder's ring background is binned about the centre. Scoring a trial centre against - // spots found at the file's puts the trial in a world the run never enters: measured, a centre - // scored 42/60 as a trial and 5/60 when that same value was given on the command line. - // - // A trial therefore gets its own azimuthal mapping, its own engines and its own cache, and the - // starting centre's spots are PARKED rather than thrown away - a hypothesis that is not adopted - // has to leave the rest of the run exactly as it found it. - std::unique_ptr trial_mapping; - std::map> home_spot_cache; - const auto try_beam_center = [&](float x, float y) { - if (trial_mapping) - spot_cache.clear(); - else - home_spot_cache.swap(spot_cache); - experiment_.BeamX_pxl(x).BeamY_pxl(y); - release_spot_engines(); - trial_mapping = std::make_unique(experiment_, pixel_mask_, - *azint_geometry_); - spot_mapping = trial_mapping.get(); - }; - const auto restore_beam_center = [&](float x, float y) { - if (!trial_mapping) - return; - experiment_.BeamX_pxl(x).BeamY_pxl(y); - release_spot_engines(); - trial_mapping.reset(); - spot_mapping = &mapping; - spot_cache = std::move(home_spot_cache); - home_spot_cache.clear(); - }; - - // How many validation frames a candidate global lattice actually indexes per image, scored on - // the *real* per-image path (force the candidate, then run the same indexing+refinement the - // main loop uses). A cell that merely fits the accumulated first-pass cloud (a spurious - // sub-lattice from an under-sampled scheme) indexes few frames, while the true cell indexes - // many - this is the discriminator that lets the best scheme win. - const std::vector validation = select_equally_spaced_image_ordinals( - early_images, std::min(early_images, 60)); - - const float index_tol = experiment_.GetIndexingSettings().GetTolerance(); - - // The validation frames' spots in reciprocal space, under a candidate's own refined geometry - // and axis. The same construction the rotation indexer makes on its accumulated frames, but - // over the validation set and over each frame's WHOLE spot list - the arbiter counts spots and - // wants the deep sample, where the first pass wants a bounded one. A superstructure layer is - // faintest exactly where the first pass's accumulation cap cuts. - auto validation_cloud = [&](const RotationIndexerResult &r) { - std::vector cloud; - if (!r.axis) - return cloud; - for (const int ordinal : validation) { - const auto it = spot_cache.find(ordinal); - if (it == spot_cache.end()) - continue; - const float ang = rot_angle(ordinal).value_or( - r.axis->GetAngle_deg(static_cast(ordinal)) + r.axis->GetWedge_deg() / 2.0f); - const auto rot = r.axis->GetTransformationAngle(ang); - for (const auto &sp : it->second) - cloud.emplace_back(rot * sp.ReciprocalCoord(r.geom)); - } - return cloud; - }; - - auto validation_settings = config_.spot_finding; - validation_settings.indexing = true; - validation_settings.quick_integration = false; - auto count_indexed = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int { - idx.ForceRotationIndexerResult(r); - // Fill the spot cache first. It is a plain map filled on demand, and the loop below runs the - // frames in parallel, where a lookup racing an insert is not something the map survives. - // After the first candidate this costs nothing. - prefetch_spots(validation); - // Each frame's verdict is its own - the score is just how many of them index - so they run - // together. Refining one frame's geometry is a few Ceres solves and no GPU, and it is the - // same call the main image loop already makes from every one of its workers on this same - // IndexAndRefine, writing only unit_cells[] under its own mutex. With the candidate forced, - // GetLattice() returns it and the branch that would advance the indexer's own state is never - // reached. The offline solver stops on an iteration count rather than a clock - // (IndexAndRefine::RefineGeometryIfNeeded), so a loaded machine cannot change a verdict. - std::atomic count = 0; - ParallelFor(static_cast(validation.size()), - std::min(validation.size(), config_.nthreads), [&](int i) { - DataMessage m{}; - m.number = validation[i]; - m.spots = spot_cache.at(validation[i]); - if (idx.IndexFrameOnly(m, validation_settings)) - ++count; - }); - return count.load(); - }; - - // The run's acceptance test counts SPOTS, not frames. Summing each validation frame's spots - // before any threshold lets a four-spot frame contribute four spots of evidence instead of an - // automatic "no", and removes the spots-per-frame denominator that makes a frame count measure - // how sparse the pattern is rather than how good the lattice is: with 4 spots on an average - // frame no lattice, the true one included, can clear a bar of six indexed spots per frame. - // - // `shift` displaces each frame's spots onto another validation frame's spindle angle. The - // lattice, the spots, the geometry and the per-frame orientation polish are all unchanged and - // only the claim under test is removed - that these spots are at these angles - so what a - // displaced run still puts on the lattice is what chance and the polish give for free. That is - // the bar the real count is measured against, and it is measured rather than chosen, so no - // spots-per-frame number survives anywhere in the test. - auto pooled_on_lattice = [&](IndexAndRefine &idx, const RotationIndexerResult &r, - int shift) -> std::pair { - idx.ForceRotationIndexerResult(r); - prefetch_spots(validation); - const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); - std::mutex reflections_mutex; - std::set> reflections; - std::atomic on_lattice = 0; - ParallelFor(static_cast(validation.size()), - std::min(validation.size(), config_.nthreads), [&](int i) { - DataMessage m{}; - m.number = validation[(i + shift) % validation.size()]; - m.spots = spot_cache.at(validation[i]); - on_lattice += idx.CountSpotsOnLattice(m, validation_settings); - // The spots the count took are the ones it gave a Miller index (none is 0,0,0). - std::unique_lock ul(reflections_mutex); - for (const auto &sp : m.spots) - if ((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0)) - reflections.insert({sp.h, sp.k, sp.l}); - }); - return {on_lattice.load(), static_cast(reflections.size())}; - }; - - // The pooled evidence for one candidate lattice: the validation frames' non-ice spots, how many - // of them lie on the candidate, and how many a wrong spindle angle still puts there. - using PooledEvidence = ValidationSpotEvidence; - - // Displacements are a fixed fraction of the sweep, so the null is a function of the data alone - // and the same file gives the same verdict every time. Sevenths: no crystallographic rotation - // is a seventh of a turn, so a displacement cannot coincide with a symmetry about the spindle - // and hand the null the signal's own score. The median over five of them bounds what one - // unlucky displacement can do. - auto pooled_evidence = [&](IndexAndRefine &idx, const RotationIndexerResult &r) { - PooledEvidence e; - std::tie(e.on_lattice, e.reflections) = pooled_on_lattice(idx, r, 0); - const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); - for (const int ordinal : validation) - for (const auto &sp : spot_cache.at(ordinal)) - if (index_ice_rings || !sp.ice_ring) - e.spots++; - const int n = static_cast(validation.size()); - std::vector chance; - for (int k = 1; k <= VALIDATION_NULL_DISPLACEMENTS; k++) { - const auto [spots_on, reflections_on] = pooled_on_lattice(idx, r, std::max(1, k * n / 7)); - chance.push_back(spots_on); - e.reflections_by_chance += reflections_on; - } - std::sort(chance.begin(), chance.end()); - e.by_chance = chance[chance.size() / 2]; - return e; - }; - - // This crystal's lattice is the one that puts materially more reflections on itself than a - // wrong angle does (ValidationEvidenceBeatsChance): every displaced angle's count is used, so - // the bar carries the null's own noise, and it carries no absolute count; the significance is - // the one-false-positive-in-a-thousand bar the spot-budget measurement is already set at. - auto beats_chance = [](const PooledEvidence &e) { - return ValidationEvidenceBeatsChance(e); - }; - - // Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the - // pooled spots each one puts on itself over and above what a wrong spindle angle puts there. - // Each is scored against its own null, which is what makes the comparison fair between a cell - // and its axis harmonic. - - // How deep into an image's intensity-ordered spot list this lattice is still being seen - the - // measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the - // spots rather than the frames: each frame's spots are tallied by rank and the run keeps the - // ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than - // that admits detections that are not this crystal's reflections, and the frame gate then counts - // them against the frame - which is how a larger budget can integrate FEWER images. - auto measure_budget = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int64_t { - idx.ForceRotationIndexerResult(r); - prefetch_spots(validation); - std::mutex evidence_mutex; - std::vector indexed(experiment_.GetMaxSpotCount(), 0); - std::vector counted(indexed.size(), 0); - const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings(); - ParallelFor(static_cast(validation.size()), - std::min(validation.size(), config_.nthreads), [&](int i) { - DataMessage m{}; - m.number = validation[i]; - m.spots = spot_cache.at(validation[i]); - if (!idx.IndexFrameOnly(m, validation_settings)) - return; - std::unique_lock ul(evidence_mutex); - AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted); - }); - return SpotBudgetFromEvidence(indexed, counted); - }; - - // Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to - // be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately - // (RunIndexing) so the two schemes' indexing can overlap. - // - // The feed itself stays in frame order and stops where it always did, but it keeps the cache a - // chunk of frames ahead of itself so those frames' spots are found on several workers. A chunk - // rather than the whole list because the "wedge" scheme stops on the accumulation cap, and how - // far into the sweep that is only becomes known as the frames arrive. - constexpr int SPOT_PREFETCH_CHUNK = 64; - size_t first_pass_spots_per_image = 0; // 0 = the indexer's own default - // x: the experiment the indexer is to read, experiment_ unless given. - auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector &ordinals, - const DiffractionExperiment *x = nullptr) { - auto ri = std::make_unique(x ? *x : experiment_, pool); - ri->MaxSpotsPerImage(first_pass_spots_per_image); - ri->VerifyMemo(config_.verify_first_pass_memo); - for (size_t i = 0; i < ordinals.size(); i++) { - if (cancelled_ || ri->AccumulationFull()) - break; - if (!spot_cache.contains(ordinals[i])) - prefetch_spots({ordinals.begin() + i, - ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)}); - ri->ProcessImage(ordinals[i], spot_cache.at(ordinals[i]), rot_angle(ordinals[i])); - } - return ri; - }; - - // Sampling schemes, tried in order of preference. We index each and keep the one that indexes - // the most frames; the schemes cover complementary failure modes and no single one is best for - // every crystal: - // - "spread": images spread over the first part of the rotation (PrepassEnd) - wide angular range, - // coarse stride; the long-standing default, best for typical crystals. - // - "wedge": consecutive frames from the clean start - native stride keeps each reflection's - // rocking curve continuous across frames, so the FFT can trace a long axis whose fine - // reciprocal spacing the coarse spread cannot resolve (e.g. an F-cubic cell viewed down - // [111]); bounded by the accumulation cap to the least radiation-damaged early wedge. - std::vector all_consecutive(early_images); - std::iota(all_consecutive.begin(), all_consecutive.end(), 0); - const std::vector spread_ordinals = select_equally_spaced_image_ordinals( - early_images, config_.rotation_indexing_image_count); - const std::vector>> schemes = { - {"spread", spread_ordinals}, - {"wedge", all_consecutive}, - }; - // A third sample of the same frames, asked only when nothing else indexes a majority (see - // add_end_wedge): the wedge read back from the END of the first part. Which orientations a sweep - // passes through first is chance, and the first part can open on frames that see the long axis - // badly - measured on a sparse crystal with a ~390 A axis, the spread and the wedge over its - // first 120 deg returned wrong cells (5/60 and 0/60 validation frames) at every spread density - // tried, while consecutive frames from 60 deg on indexed 60/60. Still inside [start, - // PrepassEnd), so nothing before pass 2 waits for the rest of the sweep. - const std::vector>> end_wedge = { - {"end wedge", {all_consecutive.rbegin(), all_consecutive.rend()}}, - }; - - // indexer_error: what the indexer threw, if it threw. An empty result with no error is a - // verdict on the crystal; an empty result with one is a verdict on the machine. - struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; - std::optional indexer_error; }; - - // Feed every scheme (single-threaded spot cache), ready to be indexed. - using SchemeIndexers = std::vector>; - auto feed_schemes = [&](IndexerThreadPool &pool, const DiffractionExperiment *x = nullptr) { - SchemeIndexers ris; - ris.reserve(schemes.size()); - for (const auto &[name, ordinals] : schemes) { - if (cancelled_) - break; - ris.push_back(feed_scheme(pool, ordinals, x)); - } - return ris; - }; - - // Run RunIndexing() of every fed scheme in the list at once. Each indexer reads only what it - // was fed and its own experiment, so the answers do not depend on what runs beside them. - auto run_indexing = [](const std::vector &ris) { - std::vector> index_futures; - index_futures.reserve(ris.size()); - for (RotationIndexer *rp : ris) - index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); - for (auto &f : index_futures) - f.get(); - }; - - // Keep the lattice of the indexed schemes that indexes the most validation frames on the real - // per-image path. Scored against the spot lists and settings the run stands on now, which must - // be those the schemes were fed under. - // `named` names the indexers (`schemes` unless given) and `bp` is the standing pass they are - // judged against (none unless given). - auto score_schemes = [&](const SchemeIndexers &ris, IndexAndRefine &idx, - const std::vector>> *named = nullptr, - FirstPass bp = {}) -> FirstPass { - if (!named) - named = &schemes; - for (size_t i = 0; i < ris.size(); i++) { - if (cancelled_) - break; - if (!bp.indexer_error) - bp.indexer_error = ris[i]->GetIndexerError(); - auto found = ris[i]->GetLattice(); - if (!found.has_value()) - continue; - RotationIndexerResult result = std::move(*found); - int score = count_indexed(idx, result); - const std::string &name = (*named)[i].first; - logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score, - static_cast(validation.size())); - - // Metric symmetry that costs nearly every frame is not this crystal's symmetry. The - // Bravais class is decided on the UNREFINED FFT candidate against a fixed angular - // tolerance, so a lattice that is pseudo-symmetric to a few tenths of a degree gets - // promoted a class too far; the constraint then snaps a real angle to the ideal one, - // which throws nearly every reflection of every frame out of tolerance (measured: 2/60 - // frames for a monoclinic lattice promoted to C-orthorhombic, 39/60 for its primitive - // cell). The indexer's own guard compares the two on the accumulated first-pass cloud, - // where they differ by less than a factor 2 - far too little to act on. - // The bar is a clear majority, not a margin: the unconstrained refinement holds NO cell - // parameter fixed, so it can only index at least as many frames as the constrained one, - // and on genuine symmetry it does index a few more (measured: 47->54 frames on a real - // I-centred orthorhombic crystal, 49->58 on an F-cubic one - a 10% margin demotes both). - // Only a constrained cell that fails outright while its unconstrained cell works is - // evidence of a false promotion, so demand exactly that. It is the same "fails to index - // half the frames" test the long-axis rescue below uses, and it leaves genuine symmetry - // alone; the space group is in any case decided from the intensities later. - const int majority = static_cast(validation.size()) / 2; - if (result.unconstrained && score < majority) { - RotationIndexerResult alt = *result.unconstrained; - const int alt_score = count_indexed(idx, alt); - if (alt_score > majority) { - logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained " - "cell indexes {}/{} - the metric symmetry is a false promotion, dropping it", - name, result.search_result.centering, - gemmi::crystal_system_str(result.search_result.system), score, - static_cast(validation.size()), alt_score, - static_cast(validation.size())); - result = std::move(alt); - score = alt_score; - } - } - // The demoted arm of an angle promotion (RunAllPasses): the free refinement instead. - if (lattice_angles_demoted_ && result.unconstrained - && (result.search_result.system == gemmi::CrystalSystem::Monoclinic - || result.search_result.system == gemmi::CrystalSystem::Orthorhombic)) { - RotationIndexerResult alt = *result.unconstrained; - result = std::move(alt); - score = count_indexed(idx, result); - } - - // Compare PRIMITIVE volumes: two schemes can find the same lattice in different - // settings, and a centred setting's cell is an exact integer multiple of its primitive - // one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive - // rhombohedral cell. Comparing the centred volumes makes the supercell test below fire - // on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which - // is enough to change the space group the search then picks. - const double vol = std::abs( - result.lattice.ToPrimitive(result.search_result.centering).CalcVolume()); - - // A later scheme wins if it indexes clearly more frames (>10%). - const bool clearly_more = static_cast(score) > bp.score * 1.1f + 0.5f; - - // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis - // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the - // same frame total and the count alone cannot tell them apart; the default then keeps - // whichever ran first. When the two schemes tie on frames but their cell volumes are - // related by an integer factor >=2, which of the pair is real is decided on the SPOTS by - // the arbiter below - regardless of scheme order, and in either direction. A near-integer - // ratio separates a real axis multiplication from a centering coincidence. Volumes are - // primitive (see above), so a pure setting difference is a ratio of 1 and never fires. - bool integer_subcell = false; - bool harmonic_settled = false; - double harmonic_ratio = 0.0; - if (bp.result.has_value() && vol > 1.0 && bp.vol > 1.0) { - const bool tied = static_cast(score) >= bp.score * 0.9f - 0.5f; - const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol; - const double nearest = std::round(ratio); - // Bounded at 4. The window is absolute in a unit integer spacing, so a ratio - // between two UNRELATED lattices passes it about a third of the time whatever n - // is, and at n = 124 the pair is not a cell and its harmonic at all. Measured - // over the corpus: all 36 firings at n <= 4 read 1.1 to 51.3 points BELOW the - // chance occupancy (n-1)/n, so an index-n sub-lattice really exists in each; all - // 8 above it read within 7.8 points OF chance, so none does. The largest n at - // which a real sub-lattice was ever seen is 3. - const bool integer_multiple = nearest >= 2.0 && nearest <= 4.0 - && std::abs(ratio - nearest) < 0.15; - - // Ask the data which of the pair is right instead of assuming the larger is - // spurious: which cell accounts for more of the validation frames' spots (see - // IndexedMillerIndices)? The measurement is made on their FULL spot lists, which - // reach far deeper into each frame's intensity distribution than the first pass's - // own accumulation cap - and a superstructure layer is faintest exactly where that - // cap cuts. - if (tied && !clearly_more && integer_multiple) { - const bool cand_is_larger = vol > bp.vol; - const RotationIndexerResult &larger = cand_is_larger ? result : *bp.result; - const RotationIndexerResult &smaller = cand_is_larger ? *bp.result : result; - const std::vector cloud = validation_cloud(larger); - const auto hkl = IndexedMillerIndices(larger.lattice, cloud, index_tol); - const size_t n_smaller = IndexedMillerIndices(smaller.lattice, - validation_cloud(smaller), - index_tol).size(); - // The added class on the PRIMITIVE basis: a centred basis has absences of its - // own, which would read as an empty class. Both cells are counted on their - // Niggli-reduced bases, because a tolerance on fractional Miller indices depends - // on the basis, and a skewed one makes the sub-lattice miss spots it contains. - const CrystalLattice larger_p = larger.lattice.ToPrimitive(larger.search_result.centering) - .NiggliReduce(); - const auto hkl_p = IndexedMillerIndices(larger_p, cloud, index_tol); - const auto ev = AddedClassOccupancy(hkl_p, static_cast(nearest)); - const CrystalLattice sub = SubLattice(larger_p, ev, static_cast(nearest)).NiggliReduce(); - const size_t n_sub = IndexedMillerIndices(sub, cloud, index_tol).size(); - const bool larger_is_real = hkl.size() > n_smaller && hkl_p.size() > n_sub; - logger.Info("Axis-harmonic arbiter: '{}' ({:.0f} A^3) vs '{}' ({:.0f} A^3), {:.2f}x - " - "the larger cell accounts for {} validation spots against the smaller " - "cell's {}; on its own geometry {} against its ({},{},{}) sub-lattice's " - "{}, {:.1f}% of them off it, so the class it adds is {}", - name, vol, bp.name, bp.vol, ratio, hkl.size(), n_smaller, - hkl_p.size(), ev.u, ev.v, ev.w, n_sub, 100.0 * ev.occupancy, - larger_is_real ? "occupied enough to pay for the longer axis - the " - "larger cell is real" - : "not occupied enough to pay for the longer axis - the " - "larger cell is a multiple"); - // Adopt the candidate exactly when its role matches the verdict. - integer_subcell = (cand_is_larger == larger_is_real); - // Only a win for the larger cell settles the pair. Declining it is a fall back - // to the default, and the warning below - which names the beam centre - is then - // still the most useful thing the run can say. It is the right diagnosis in - // particular for the failure this arbiter cannot see: a centre error along the - // spindle translates the derotated cloud rigidly by HALF a lattice spacing, and - // a half-offset lattice is indexable only on a doubled axis. Such a cell fits no - // index-n sublattice at all, so its added class reads the chance value (n-1)/n - - // the largest the occupancy can be - and it accounts for far fewer spots than - // the true cell, which is what this test sees it by. - harmonic_settled = larger_is_real; - } - // sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by - // sqrt(3) is the other harmonic seen on real data and its volume ratio is not a - // whole number. This is a REPORT and not a decision - it is deliberately taken - // over the whole pair, `clearly_more` included, because the harmonic is often the - // one that indexes more frames (see below). - if (integer_multiple || std::abs(ratio - std::sqrt(3.0)) < 0.15) - harmonic_ratio = ratio; - } - - // Two candidates whose primitive volumes differ by a small integer or by sqrt(3) are - // one lattice and its axis harmonic, and what decides between them is the BEAM CENTRE. - // A centre error along the spindle translates the derotated cloud rigidly, which the - // FFT amplitude cannot see at all, so the peaks stay sharp and the indexer commits - - // with confidence - to a sub-multiple. Measured on real data: a 0.12 px change of - // centre, 0.03x of what the J0 law asks across the spindle, is the whole difference - // between the deposited cell and a halved axis, and the halved one is the answer that - // indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate - // this pair, and until this line nothing in the run said what the pair meant. - // ... unless the arbiter above established the larger cell on the data, which is - // exactly the question this warning says the run cannot answer. - if (harmonic_ratio > 0.0 && !harmonic_settled) - logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) " - "differ by {:.2f}x - one is an axis harmonic of the other. Which of " - "the two is right is settled by the beam centre to a fraction of a " - "pixel, and the harmonic is often the one that indexes MORE frames, " - "so the frame counts above cannot decide it. Check the beam centre", - name, vol, bp.name, bp.vol, harmonic_ratio); - - if (!bp.result.has_value() || clearly_more || integer_subcell) { - if (integer_subcell && vol < bp.vol) - logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " - "equal frame count - adopting the smaller true cell", name, vol, - std::round(bp.vol / vol), bp.name, bp.vol); - else if (integer_subcell) - logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x supercell of '{}' (vol {:.0f}) at " - "equal frame count and accounts for more of the spots - adopting the " - "larger true cell", name, vol, std::round(vol / bp.vol), bp.name, bp.vol); - bp.score = score; - bp.vol = vol; - bp.name = name; - bp.result = std::move(result); - } - } - return bp; - }; - - // Where the first pass indexes no majority of the validation frames, the end wedge is asked too - // and judged against the standing pass as any later scheme is. Only then: on a run the schemes - // index it would be one more lattice for the axis-harmonic arbiter, which is not decided on - // frame counts. Asked once, after the axis-sign, beam-centre and long-axis rescues (below), not - // inside every pass they run: a consecutive wedge indexes part of the validation frames at the - // WRONG axis sign too (measured: 33/60 where the spread at the right sign indexes 60/60), and - // counted in the pass a rescue is gated on, it lifts that pass past the gate and the rescue - // that would have indexed every frame is never asked. - auto add_end_wedge = [&](FirstPass bp, IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass { - if (cancelled_ || bp.score > static_cast(validation.size()) / 2) - return bp; - SchemeIndexers ris; - ris.push_back(feed_scheme(pool, end_wedge[0].second)); - run_indexing({ris[0].get()}); - return score_schemes(ris, idx, &end_wedge, std::move(bp)); - }; - - // Feed the schemes, index them concurrently, and keep the best (score_schemes). - // alongside: started once the schemes are fed and before they are indexed, to run beside them. - auto pick_best = [&](IndexerThreadPool &pool, IndexAndRefine &idx, - const std::function &alongside = {}) -> FirstPass { - const SchemeIndexers ris = feed_schemes(pool); - if (alongside) - alongside(); - std::vector rp; - for (const auto &ri : ris) - rp.push_back(ri.get()); - run_indexing(rp); - return score_schemes(ris, idx); - }; - - // Powder-ring rescue - and, with it, the lean seed the rotation first pass never had. - // - // The pre-scan has already measured this sample's powder rings (see PreScan): a crystalline - // phase other than the crystal, whose reflections the indexer is otherwise handed as if they - // were this crystal's. Hexagonal ice is the only such phase whose rings can be named in - // advance, and the spot finder flags it from that fixed list; everything else - a shower of - // microcrystals, a salt out of the cryoprotectant, a phase nobody expected - has to be - // measured, and until it is, every one of its spots goes into the FFT. Measured on one such - // crystal: 876 spots a frame against the 80 a clean crystal on the same beamline gives, 70 % - // of them on rings, no scheme indexing a single validation frame, and the FFT taking a 10-30 A - // cell out of the shells. - // - // Setting the rings aside is half of it. The other half is HOW DEEP into each frame's spot - // list the first pass reads. It has always read 200 spots an image; the per-frame indexer, on - // the same problem, escalates 30 -> 80 -> the whole list and keeps the seed a lattice explains - // the largest fraction of (IndexAndRefine::DetermineLatticeAndSymmetry), because on a flooded - // frame a lean high-quality seed finds the lattice where the full list derails the search. The - // rotation pass gets the same escalation here. It is the general statement of what a - // resolution cut would do on this one crystal - the contaminant and the noise sit at the deep - // end of the intensity-ordered list wherever in resolution they happen to lie - and it needs - // nothing measured about where they are. - // - // It runs LAST: after every other rescue and after the pass-1 fallback above, on a run that - // would otherwise stop with no lattice at all. Measured over 100 datasets: placed earlier it - // pre-empted the rescues that already work - a crystal the beam-centre ladder fixes at 59/60 - // and one the pass-1 fallback fixes at 60/60 both came out WORSE, because this ladder had - // already replaced the standing pass those rescues start from. A last resort has to be last. - // - // The two compose: FilterSpotsByCount has already ordered each image's spots non-ring first - // and strongest first, so a lean cap takes the crystal's strongest reflections and nothing - // else. Decided late, the way the axis sign above is: index each rung, keep the one that - // indexes the most validation frames. A run that already indexes never reaches here, and one - // that does cannot lose by trying - where the rings are the crystal's own rows, or where the - // crystal needs its deep spots to index at all, the pass that already ran stands. - - // Called from two places, because the two questions it answers are coupled: a centre that is - // wrong and a spot list that is too deep each hide the other, exactly as the centre and the - // axis sign do above. `score_before` is what the caller's standing pass indexes and `min_gain` - // is what a rung must beat it by; nothing is changed and the settings are put back when no - // rung does. - // What the run is CURRENTLY reading each frame at. It starts at what the file asked for, and - // only an ADOPTED rung moves it. Every rung the ladder tries is put back to this, so a rung - // tried and rejected leaves nothing behind - and a later pass that finds no rung restores what - // an earlier pass adopted rather than resetting the run to the defaults. - const auto first_pass_file_d_min = config_.spot_finding.high_resolution_limit; - FirstPassState committed{false, 0, first_pass_file_d_min}; - const auto apply_state = [&](const FirstPassState &st) { - const std::vector r = st.rings ? powder_.rings_q_recipA : std::vector{}; - config_.spot_finding.measured_ring_q_recipA = r; - validation_settings.measured_ring_q_recipA = r; - config_.spot_finding.high_resolution_limit = st.d_min; - validation_settings.high_resolution_limit = st.d_min; - first_pass_spots_per_image = st.seed; - spot_cache.clear(); // the per-image spot budget is spent again under these settings - }; // The measured-centre arm of a carried beam-centre hypothesis reads each frame as deep as the // check found that centre indexing at: the arm is the hypothesis the check measured, not a new // search at the same centre that may end somewhere else. if (beam_center_alternative_state_) { - committed = *beam_center_alternative_state_; + run.committed = *beam_center_alternative_state_; beam_center_alternative_state_.reset(); - apply_state(committed); + run.ApplyState(run.committed); } - // min_score is an ABSOLUTE bar beside the gain. A caller that takes a rung only when it - // indexes a majority must not be handed a state the ladder committed on a lesser rung: - // measured, the ladder adopted at 29/60, the caller rejected it against a majority of 30, and - // nothing put the settings back - the beam-centre walk that followed then scored every - // hypothesis on the leaked lean spot lists. - const auto first_pass_ladder = - [&](int score_before, int min_gain, int min_score) -> std::optional { - if (cancelled_) - return std::nullopt; + run.PrefetchCertainFrames(); + FirstPass best = run.IndexSchemes(); + best = run.AxisSignRescue(std::move(best)); + bool long_axis_asked = false; + best = run.BeamCentreCheck(std::move(best), long_axis_asked); + // The long-axis rescue, at the centre the run now stands on. + if (!long_axis_asked) + if (auto rescued = run.LongAxisRescue(best)) + best = *rescued; + // The end wedge, at the geometry the rescues settled on. + best = run.AddEndWedge(std::move(best), *p.indexer_pool, *p.indexer); + run.ReleaseSpotEngines(); // no more frames to read; give the cards their memory back + if (cancelled_) + return false; - // The seeds, leanest first. The same ladder the per-frame indexer uses, one rung deeper - // because a rotation cloud pools many frames and can afford more. - // The last rung is EVERY spot of each frame: SIZE_MAX, not 0 - RotationIndexer reads 0 as - // "use my default", which is 200, so a 0 here made the deepest rung a duplicate of the - // 200 one and the log line "the whole spot list" a lie. - constexpr std::array SEEDS = {30, 80, 200, std::numeric_limits::max()}; - struct Rung { bool rings; size_t seed; std::optional d_min; }; - // How much of each frame a rung let the first pass see: more spots an image, and further - // in resolution, is more evidence. Rings set aside is not part of it - that is a - // measurement of the sample, not a depth. - const auto Deeper = [](const Rung &a, const Rung &b) { - if (a.seed != b.seed) return a.seed > b.seed; - const float da = a.d_min.value_or(0.0f), db = b.d_min.value_or(0.0f); - return da < db; // a smaller d_min reaches further; 0 means the detector edge - }; - // The file's OWN limit, captured before any rung committed anything. Reading it out of - // config_ here would read back the previous rung's floor once a pass has adopted one, and - // the restore path would then "put back" that instead of what the file asked for. - const auto file_limit = first_pass_file_d_min; - std::vector rungs; - // Where a floor is worth trying at all - the rings merged somewhere inside the measured - // range, so past that point nothing can be told from the crystal and a pass that reads - // there is reading the contaminant. - std::vector> floors = {file_limit}; - for (const float d : prescan_spot_quantile_d_A_) - if (!file_limit || d > *file_limit) floors.emplace_back(d); - // Floors outer, seeds inner, so rungs that share the finder-facing settings (rings, - // d_min) sit adjacent: the seed is an indexer-side cap over the SAME found spots - // (MaxSpotsPerImage over the cached list), so the loop below spot-finds once per floor - // instead of once per rung. Adoption is decided over all rungs at once, so the visit - // order carries no meaning. - for (const auto &d_min : floors) - for (const size_t seed : SEEDS) { - if (!powder_.rings_q_recipA.empty()) - rungs.push_back({true, seed, d_min}); - if (powder_.rings_q_recipA.empty() && seed != SEEDS.back()) - rungs.push_back({false, seed, d_min}); - } + run.RequireLattice(best); + best = run.ShortAxisHypothesis(std::move(best)); + best = run.SubLatticeHypothesis(std::move(best)); + if (!p.geometry_prepass && prepass_result_.has_value()) + best = run.Pass1Fallback(std::move(best)); + // The last resort, and ONLY on a run that has no usable lattice: the pass-1 fallback above + // has had its turn, and below this the run stops with an error. A run that indexes must + // not pay for it - measured, an unguarded call spent 41-76 % of the wall clock on a + // healthy dataset evaluating rungs that arithmetically could not be adopted (one crystal + // 41 s guarded against 175 s unguarded, same answer). + if (best.score < n_validation / 6) { + if (auto rescued = run.FirstPassLadder(best.score, n_validation / 6, 0)) + best = *rescued; + } - // A new floor (or ring handling) changes what the finder returns, so the spot lists are - // spent and found again; a new seed alone does not - it caps how much of the SAME lists - // the indexer reads - and re-finding 60 validation frames per rung was most of what a - // hopeless run's minutes were spent on. - const auto enter_rung = [&](size_t i) { - if (i == 0 || rungs[i].rings != rungs[i - 1].rings || rungs[i].d_min != rungs[i - 1].d_min) - apply_state({rungs[i].rings, rungs[i].seed, rungs[i].d_min}); + const FirstPassRun::Verdict verdict = run.MeasureVerdict(best, first_pass_key, memo_hit); + // A pass run only to score a lattice has its score, whatever it is: a lattice that does not + // beat chance is a result for the comparison that asked, not a reason to stop the run. + if (indexing_probe_only_) { + logger.Info("Indexing probe: scheme '{}', {}/{} validation frames, {}/{} validation spots " + "= {:.1f}% against {:.1f}% at a wrong spindle angle", best.name, best.score, + n_validation, verdict.evidence.on_lattice, verdict.evidence.spots, + 100.0 * verdict.pooled, 100.0 * verdict.pooled_chance); + return true; + } + run.RejectChanceLattice(best, verdict); + run.AdoptLattice(best, verdict); + if (!p.geometry_prepass && !postrefine_probe_only_) + run.StartLeftoverCensus(best); + if (best.result->tilt_walk) + run.ReportTiltWalk(best); + if (config_.measure_spot_budget) + run.MeasureSpotBudget(best); + if (!p.geometry_prepass && prepass_result_.has_value()) + run.Pass1SupercellGuard(best); + return false; +} + +Rugnux::FirstPassRun::FirstPassRun(Rugnux &rx, PipelineLocals &p) + : rx(rx), logger(p.logger), result(p.result), geometry_prepass(p.geometry_prepass), + start_image(p.start_image), early_images(p.prepass_images), dataset(p.dataset), mapping(p.mapping), + indexer_pool(p.indexer_pool), indexer(p.indexer), indexing_bound_A(p.indexing_bound_A), + census_future(p.census_future), + spot_workers(get_gpu_count() > 0 + ? std::min(rx.config_.nthreads, std::max(8, 4 * get_gpu_count())) + : rx.config_.nthreads), + engines(std::max(spot_workers, 1)), + validation_settings(rx.config_.spot_finding), + all_consecutive(early_images), + spread_ordinals(select_equally_spaced_image_ordinals(early_images, + rx.config_.rotation_indexing_image_count)), + first_pass_file_d_min(rx.config_.spot_finding.high_resolution_limit), + committed{false, 0, first_pass_file_d_min} { + validation_settings.indexing = true; + validation_settings.quick_integration = false; + std::iota(all_consecutive.begin(), all_consecutive.end(), 0); + schemes = { + {"spread", spread_ordinals}, + {"wedge", all_consecutive}, + }; + end_wedge = { + {"end wedge", {all_consecutive.rbegin(), all_consecutive.rend()}}, + }; +} + +// Mid-exposure goniometer angle for an image ordinal (matches the per-image path). +std::optional Rugnux::FirstPassRun::RotAngle(int ordinal) const { + if (const auto g = rx.experiment_.GetGoniometer()) + return g->GetAngle_deg(static_cast(ordinal)) + g->GetWedge_deg() / 2.0f; + return std::nullopt; +} + +// The spots of one frame, found with the caller's own analysis engine. A frame's spots are a +// pure function of that frame and the settings - the engine carries nothing from one image to +// the next, which is why the main image loop can give one to each of its workers - so it does +// not matter which engine, or which thread, finds them. +std::vector Rugnux::FirstPassRun::FindSpots(MXAnalysisWithoutFPGA &analysis, + AzimuthalIntegrationProfile &profile, + JFJochReaderRawImage &img, int ordinal) { + return FindFirstPassSpots(rx.reader_, rx.config_.spot_finding, *dataset, analysis, profile, img, + ordinal, start_image + ordinal * rx.config_.stride, logger); +} + +void Rugnux::FirstPassRun::UseGpuOfWorker(size_t t) { + if (const int32_t gpus = get_gpu_count(); gpus > 0) + pin_gpu(static_cast(t) % gpus); +} + +// Find the spots of every ordinal in the list that is not cached yet, on several workers, and +// put them in the cache. +// +// The cache ends up holding exactly what a serial loop would have put there: the results are +// collected into a vector indexed by position and inserted afterwards, so neither the spots of +// a frame nor the set of frames in the cache depends on how the workers interleaved. +// +// Worker t keeps its analysis engine for the whole first pass. Building one is a CUDA stream +// and a full set of device buffers - measured at 0.13 s against 4 ms for a frame - and the +// "wedge" scheme only discovers how far into the sweep it has to read as the frames arrive, +// so this is called several times and rebuilding them each time cost a quarter of the pass. +// The engine belongs to a card, not to a thread, so a later call only has to make that card +// current again; workers take a card by index so an engine always meets the one it was +// built on. +void Rugnux::FirstPassRun::PrefetchSpots(const std::vector &ordinals) { + // What an earlier pass already found under the same key (rx.first_pass_spots_) is taken from + // there, at this pass's angles. + const std::vector key = rx.SpotFindingKey(*spot_mapping); + std::map> kept; + { + std::lock_guard lock(rx.first_pass_spots_->m); + for (const auto &[k, spots] : rx.first_pass_spots_->entries) + if (k == key) + for (const int ordinal : ordinals) + if (const auto it = spots.find(ordinal); it != spots.end()) + kept.emplace(ordinal, it->second); + } + const DiffractionGeometry spot_geometry = rx.experiment_.GetDiffractionGeometry(); + for (auto &[ordinal, spots] : kept) + StampSpotAngles(spots, spot_geometry); + std::vector wanted; + for (const int ordinal : ordinals) + if (!spot_cache.contains(ordinal)) { + if (const auto it = kept.find(ordinal); it != kept.end() && !rx.config_.verify_first_pass_memo) + spot_cache.emplace(ordinal, it->second); else - first_pass_spots_per_image = rungs[i].seed; - }; - // Every rung is fed first and all of them are indexed at once: the rungs do not depend on - // one another, and indexing them one after the other left most of the machine idle for - // most of a hopeless run. They are then scored in order, each under its own settings again. - std::vector rung_ris(rungs.size()); - std::vector all_ris; - for (size_t i = 0; i < rungs.size() && !cancelled_; i++) { - enter_rung(i); - rung_ris[i] = feed_schemes(*indexer_pool); - for (const auto &ri : rung_ris[i]) - all_ris.push_back(ri.get()); + wanted.push_back(ordinal); } - run_indexing(all_ris); + if (wanted.empty()) + return; + const auto prefetch_start = std::chrono::steady_clock::now(); + ++prefetch_batches; - // Every rung's score, so the choice below can be made on all of them at once. - std::vector rung_score(rungs.size(), -1); - std::vector rung_pass(rungs.size()); - int won = -1; - for (size_t i = 0; i < rungs.size() && !cancelled_; i++) { - enter_rung(i); - FirstPass alt = score_schemes(rung_ris[i], *indexer); - logger.Info(" first-pass rung: {} spots/image, rings {}, stop at {} -> {}/{} frames{}", - rungs[i].seed == std::numeric_limits::max() - ? std::string("all") : std::to_string(rungs[i].seed), - rungs[i].rings ? "aside" : "kept", - rungs[i].d_min ? fmt::format("{:.2f} A", *rungs[i].d_min) - : std::string("the detector edge"), - alt.result.has_value() ? alt.score : -1, - static_cast(validation.size()), - alt.result.has_value() - ? fmt::format(", {}{} {:.1f} {:.1f} {:.1f} {:.1f} {:.1f} {:.1f}", - alt.result->search_result.centering, - gemmi::crystal_system_str(alt.result->search_result.system), - alt.result->search_result.conventional.GetUnitCell().a, - alt.result->search_result.conventional.GetUnitCell().b, - alt.result->search_result.conventional.GetUnitCell().c, - alt.result->search_result.conventional.GetUnitCell().alpha, - alt.result->search_result.conventional.GetUnitCell().beta, - alt.result->search_result.conventional.GetUnitCell().gamma) - : std::string()); - // A rung has to WIN by the count the run itself keeps a lattice on - a sixth of the - // validation frames, the bar below - not merely come out ahead. Three frames out of - // sixty is the count's noise, and adopting on it would hand the rescues that follow - // (the beam centre, the long axis) a spot list cut for no reason. Measured: a marginal - // crystal that indexed 26/60 went to 29/60 on a rung and must not take it, where the - // crystals this exists for go from 0/60 to 24 and to 32/60. - if (alt.result.has_value()) { - rung_score[i] = alt.score; - rung_pass[i] = std::move(alt); + std::vector> found(wanted.size()); + std::atomic next{0}; + std::vector> futures; + const size_t nworkers = std::min(engines.size(), wanted.size()); + futures.reserve(nworkers); + for (size_t t = 0; t < nworkers; t++) + futures.emplace_back(std::async(std::launch::async, [&, t] { + UseGpuOfWorker(t); // must precede engine construction + SpotEngine &e = engines[t]; + if (!e.analysis) { + e.analysis = std::make_unique( + rx.experiment_, *spot_mapping, rx.pixel_mask_, *indexer, + /*enable_fused_adaptive_gpu=*/true); + e.profile = std::make_unique(*spot_mapping); + } + for (size_t i = next.fetch_add(1); i < wanted.size(); i = next.fetch_add(1)) + found[i] = FindSpots(*e.analysis, *e.profile, e.raw_image, wanted[i]); + })); + for (auto &f : futures) + f.get(); + + { + std::lock_guard lock(rx.first_pass_spots_->m); + auto &entries = rx.first_pass_spots_->entries; + auto entry = std::find_if(entries.begin(), entries.end(), + [&](const auto &e) { return e.first == key; }); + if (entry == entries.end()) { + if (entries.size() == MAX_SPOT_KEYS) + entries.erase(entries.begin()); + entries.emplace_back(key, std::map>()); + entry = entries.end() - 1; + } + for (size_t i = 0; i < wanted.size(); i++) { + if (const auto it = kept.find(wanted[i]); it != kept.end() && !SameSpots(it->second, found[i])) + throw MemoMismatch(fmt::format("First-pass spot memo differs from a recomputation " + "for image {}", wanted[i])); + entry->second.emplace(wanted[i], found[i]); + } + } + for (size_t i = 0; i < wanted.size(); i++) + spot_cache.emplace(wanted[i], std::move(found[i])); + prefetch_time_s += std::chrono::duration( + std::chrono::steady_clock::now() - prefetch_start).count(); +} + +// Release the engines on workers of their own, so a stream and its device buffers are freed +// with the card they were made on current. +void Rugnux::FirstPassRun::ReleaseSpotEngines() { + std::vector> futures; + futures.reserve(engines.size()); + for (size_t t = 0; t < engines.size(); t++) + futures.emplace_back(std::async(std::launch::async, [&, t] { + UseGpuOfWorker(t); + engines[t] = SpotEngine{}; + })); + for (auto &f : futures) + f.get(); +} + +// A beam-centre hypothesis, tried at its own mapping, engines and spot cache (see trial_mapping). +void Rugnux::FirstPassRun::TryBeamCenter(float x, float y) { + if (trial_mapping) + spot_cache.clear(); + else + home_spot_cache.swap(spot_cache); + rx.experiment_.BeamX_pxl(x).BeamY_pxl(y); + ReleaseSpotEngines(); + trial_mapping = std::make_unique(rx.experiment_, rx.pixel_mask_, + *rx.azint_geometry_); + spot_mapping = trial_mapping.get(); +} + +void Rugnux::FirstPassRun::RestoreBeamCenter(float x, float y) { + if (!trial_mapping) + return; + rx.experiment_.BeamX_pxl(x).BeamY_pxl(y); + ReleaseSpotEngines(); + trial_mapping.reset(); + spot_mapping = &mapping; + spot_cache = std::move(home_spot_cache); + home_spot_cache.clear(); +} + +// The validation frames' spots in reciprocal space, under a candidate's own refined geometry +// and axis. The same construction the rotation indexer makes on its accumulated frames, but +// over the validation set and over each frame's WHOLE spot list - the arbiter counts spots and +// wants the deep sample, where the first pass wants a bounded one. A superstructure layer is +// faintest exactly where the first pass's accumulation cap cuts. +std::vector Rugnux::FirstPassRun::ValidationCloud(const RotationIndexerResult &r) { + std::vector cloud; + if (!r.axis) + return cloud; + for (const int ordinal : validation) { + const auto it = spot_cache.find(ordinal); + if (it == spot_cache.end()) + continue; + const float ang = RotAngle(ordinal).value_or( + r.axis->GetAngle_deg(static_cast(ordinal)) + r.axis->GetWedge_deg() / 2.0f); + const auto rot = r.axis->GetTransformationAngle(ang); + for (const auto &sp : it->second) + cloud.emplace_back(rot * sp.ReciprocalCoord(r.geom)); + } + return cloud; +} + +// How many validation frames a candidate indexes on the real per-image path (see validation). +int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r) { + idx.ForceRotationIndexerResult(r); + // Fill the spot cache first. It is a plain map filled on demand, and the loop below runs the + // frames in parallel, where a lookup racing an insert is not something the map survives. + // After the first candidate this costs nothing. + PrefetchSpots(validation); + // Each frame's verdict is its own - the score is just how many of them index - so they run + // together. Refining one frame's geometry is a few Ceres solves and no GPU, and it is the + // same call the main image loop already makes from every one of its workers on this same + // IndexAndRefine, writing only unit_cells[] under its own mutex. With the candidate forced, + // GetLattice() returns it and the branch that would advance the indexer's own state is never + // reached. The offline solver stops on an iteration count rather than a clock + // (IndexAndRefine::RefineGeometryIfNeeded), so a loaded machine cannot change a verdict. + std::atomic count = 0; + ParallelFor(static_cast(validation.size()), + std::min(validation.size(), rx.config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[i]; + m.spots = spot_cache.at(validation[i]); + if (idx.IndexFrameOnly(m, validation_settings)) + ++count; + }); + return count.load(); +} + +// The run's acceptance test counts SPOTS, not frames. Summing each validation frame's spots +// before any threshold lets a four-spot frame contribute four spots of evidence instead of an +// automatic "no", and removes the spots-per-frame denominator that makes a frame count measure +// how sparse the pattern is rather than how good the lattice is: with 4 spots on an average +// frame no lattice, the true one included, can clear a bar of six indexed spots per frame. +// +// `shift` displaces each frame's spots onto another validation frame's spindle angle. The +// lattice, the spots, the geometry and the per-frame orientation polish are all unchanged and +// only the claim under test is removed - that these spots are at these angles - so what a +// displaced run still puts on the lattice is what chance and the polish give for free. That is +// the bar the real count is measured against, and it is measured rather than chosen, so no +// spots-per-frame number survives anywhere in the test. +std::pair Rugnux::FirstPassRun::PooledOnLattice(IndexAndRefine &idx, + const RotationIndexerResult &r, int shift) { + idx.ForceRotationIndexerResult(r); + PrefetchSpots(validation); + const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings(); + std::mutex reflections_mutex; + std::set> reflections; + std::atomic on_lattice = 0; + ParallelFor(static_cast(validation.size()), + std::min(validation.size(), rx.config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[(i + shift) % validation.size()]; + m.spots = spot_cache.at(validation[i]); + on_lattice += idx.CountSpotsOnLattice(m, validation_settings); + // The spots the count took are the ones it gave a Miller index (none is 0,0,0). + std::unique_lock ul(reflections_mutex); + for (const auto &sp : m.spots) + if ((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0)) + reflections.insert({sp.h, sp.k, sp.l}); + }); + return {on_lattice.load(), static_cast(reflections.size())}; +} + +// Displacements are a fixed fraction of the sweep, so the null is a function of the data alone +// and the same file gives the same verdict every time. Sevenths: no crystallographic rotation +// is a seventh of a turn, so a displacement cannot coincide with a symmetry about the spindle +// and hand the null the signal's own score. The median over five of them bounds what one +// unlucky displacement can do. +// +// This crystal's lattice is the one that puts materially more reflections on itself than a +// wrong angle does (ValidationEvidenceBeatsChance): every displaced angle's count is used, so +// the bar carries the null's own noise, and it carries no absolute count; the significance is +// the one-false-positive-in-a-thousand bar the spot-budget measurement is already set at. +// +// Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the +// pooled spots each one puts on itself over and above what a wrong spindle angle puts there. +// Each is scored against its own null, which is what makes the comparison fair between a cell +// and its axis harmonic. +auto Rugnux::FirstPassRun::MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r) -> PooledEvidence { + PooledEvidence e; + std::tie(e.on_lattice, e.reflections) = PooledOnLattice(idx, r, 0); + const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings(); + for (const int ordinal : validation) + for (const auto &sp : spot_cache.at(ordinal)) + if (index_ice_rings || !sp.ice_ring) + e.spots++; + const int n = static_cast(validation.size()); + std::vector chance; + for (int k = 1; k <= VALIDATION_NULL_DISPLACEMENTS; k++) { + const auto [spots_on, reflections_on] = PooledOnLattice(idx, r, std::max(1, k * n / 7)); + chance.push_back(spots_on); + e.reflections_by_chance += reflections_on; + } + std::sort(chance.begin(), chance.end()); + e.by_chance = chance[chance.size() / 2]; + return e; +} + +// How deep into an image's intensity-ordered spot list this lattice is still being seen - the +// measured spot budget. Same frames, same per-image path as CountIndexed above, but scoring the +// spots rather than the frames: each frame's spots are tallied by rank and the run keeps the +// ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than +// that admits detections that are not this crystal's reflections, and the frame gate then counts +// them against the frame - which is how a larger budget can integrate FEWER images. +int64_t Rugnux::FirstPassRun::MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r) { + idx.ForceRotationIndexerResult(r); + PrefetchSpots(validation); + std::mutex evidence_mutex; + std::vector indexed(rx.experiment_.GetMaxSpotCount(), 0); + std::vector counted(indexed.size(), 0); + const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings(); + ParallelFor(static_cast(validation.size()), + std::min(validation.size(), rx.config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[i]; + m.spots = spot_cache.at(validation[i]); + if (!idx.IndexFrameOnly(m, validation_settings)) + return; + std::unique_lock ul(evidence_mutex); + AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted); + }); + return SpotBudgetFromEvidence(indexed, counted); +} + +// Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to +// be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately +// (RunIndexing) so the two schemes' indexing can overlap. +// +// The feed itself stays in frame order and stops where it always did, but it keeps the cache a +// chunk of frames ahead of itself so those frames' spots are found on several workers. A chunk +// rather than the whole list because the "wedge" scheme stops on the accumulation cap, and how +// far into the sweep that is only becomes known as the frames arrive. +// x: the experiment the indexer is to read, experiment_ unless given. +std::unique_ptr Rugnux::FirstPassRun::FeedScheme(IndexerThreadPool &pool, + const std::vector &ordinals, + const DiffractionExperiment *x) { + auto ri = std::make_unique(x ? *x : rx.experiment_, pool); + ri->MaxSpotsPerImage(first_pass_spots_per_image); + ri->VerifyMemo(rx.config_.verify_first_pass_memo); + for (size_t i = 0; i < ordinals.size(); i++) { + if (rx.cancelled_ || ri->AccumulationFull()) + break; + if (!spot_cache.contains(ordinals[i])) + PrefetchSpots({ordinals.begin() + i, + ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)}); + ri->ProcessImage(ordinals[i], spot_cache.at(ordinals[i]), RotAngle(ordinals[i])); + } + return ri; +} + +// Feed every scheme (single-threaded spot cache), ready to be indexed. +auto Rugnux::FirstPassRun::FeedSchemes(IndexerThreadPool &pool, const DiffractionExperiment *x) -> SchemeIndexers { + SchemeIndexers ris; + ris.reserve(schemes.size()); + for (const auto &[name, ordinals] : schemes) { + if (rx.cancelled_) + break; + ris.push_back(FeedScheme(pool, ordinals, x)); + } + return ris; +} + +// Run RunIndexing() of every fed scheme in the list at once. Each indexer reads only what it +// was fed and its own experiment, so the answers do not depend on what runs beside them. +void Rugnux::FirstPassRun::RunIndexing(const std::vector &ris) { + std::vector> index_futures; + index_futures.reserve(ris.size()); + for (RotationIndexer *rp : ris) + index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); + for (auto &f : index_futures) + f.get(); +} + +void Rugnux::FirstPassRun::RunIndexing(const SchemeIndexers &ris) { + std::vector rp; + for (const auto &ri : ris) + rp.push_back(ri.get()); + RunIndexing(rp); +} + +// Keep the lattice of the indexed schemes that indexes the most validation frames on the real +// per-image path. Scored against the spot lists and settings the run stands on now, which must +// be those the schemes were fed under. +// `named` names the indexers (`schemes` unless given) and `bp` is the standing pass they are +// judged against (none unless given). +auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefine &idx, const Schemes *named, + FirstPass bp) -> FirstPass { + if (!named) + named = &schemes; + for (size_t i = 0; i < ris.size(); i++) { + if (rx.cancelled_) + break; + if (!bp.indexer_error) + bp.indexer_error = ris[i]->GetIndexerError(); + auto found = ris[i]->GetLattice(); + if (!found.has_value()) + continue; + RotationIndexerResult result = std::move(*found); + int score = CountIndexed(idx, result); + const std::string &name = (*named)[i].first; + logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score, + static_cast(validation.size())); + + // Metric symmetry that costs nearly every frame is not this crystal's symmetry. The + // Bravais class is decided on the UNREFINED FFT candidate against a fixed angular + // tolerance, so a lattice that is pseudo-symmetric to a few tenths of a degree gets + // promoted a class too far; the constraint then snaps a real angle to the ideal one, + // which throws nearly every reflection of every frame out of tolerance (measured: 2/60 + // frames for a monoclinic lattice promoted to C-orthorhombic, 39/60 for its primitive + // cell). The indexer's own guard compares the two on the accumulated first-pass cloud, + // where they differ by less than a factor 2 - far too little to act on. + // The bar is a clear majority, not a margin: the unconstrained refinement holds NO cell + // parameter fixed, so it can only index at least as many frames as the constrained one, + // and on genuine symmetry it does index a few more (measured: 47->54 frames on a real + // I-centred orthorhombic crystal, 49->58 on an F-cubic one - a 10% margin demotes both). + // Only a constrained cell that fails outright while its unconstrained cell works is + // evidence of a false promotion, so demand exactly that. It is the same "fails to index + // half the frames" test the long-axis rescue below uses, and it leaves genuine symmetry + // alone; the space group is in any case decided from the intensities later. + const int majority = static_cast(validation.size()) / 2; + if (result.unconstrained && score < majority) { + RotationIndexerResult alt = *result.unconstrained; + const int alt_score = CountIndexed(idx, alt); + if (alt_score > majority) { + logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained " + "cell indexes {}/{} - the metric symmetry is a false promotion, dropping it", + name, result.search_result.centering, + gemmi::crystal_system_str(result.search_result.system), score, + static_cast(validation.size()), alt_score, + static_cast(validation.size())); + result = std::move(alt); + score = alt_score; } } + // The demoted arm of an angle promotion (RunAllPasses): the free refinement instead. + if (rx.lattice_angles_demoted_ && result.unconstrained + && (result.search_result.system == gemmi::CrystalSystem::Monoclinic + || result.search_result.system == gemmi::CrystalSystem::Orthorhombic)) { + RotationIndexerResult alt = *result.unconstrained; + result = std::move(alt); + score = CountIndexed(idx, result); + } - // The rung to keep. A rung has to WIN by the count the run itself keeps a lattice on - a - // sixth of the validation frames, the bar below - not merely come out ahead: three frames - // out of sixty is the count's noise, and adopting on it would hand the rescues that follow + // Compare PRIMITIVE volumes: two schemes can find the same lattice in different + // settings, and a centred setting's cell is an exact integer multiple of its primitive + // one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive + // rhombohedral cell. Comparing the centred volumes makes the supercell test below fire + // on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which + // is enough to change the space group the search then picks. + const double vol = std::abs( + result.lattice.ToPrimitive(result.search_result.centering).CalcVolume()); + + // A later scheme wins if it indexes clearly more frames (>10%). + const bool clearly_more = static_cast(score) > bp.score * 1.1f + 0.5f; + + // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis + // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the + // same frame total and the count alone cannot tell them apart; the default then keeps + // whichever ran first. When the two schemes tie on frames but their cell volumes are + // related by an integer factor >=2, which of the pair is real is decided on the SPOTS by + // the arbiter below - regardless of scheme order, and in either direction. A near-integer + // ratio separates a real axis multiplication from a centering coincidence. Volumes are + // primitive (see above), so a pure setting difference is a ratio of 1 and never fires. + bool integer_subcell = false; + bool harmonic_settled = false; + double harmonic_ratio = 0.0; + if (bp.result.has_value() && vol > 1.0 && bp.vol > 1.0) { + const bool tied = static_cast(score) >= bp.score * 0.9f - 0.5f; + const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol; + const double nearest = std::round(ratio); + // Bounded at 4. The window is absolute in a unit integer spacing, so a ratio + // between two UNRELATED lattices passes it about a third of the time whatever n + // is, and at n = 124 the pair is not a cell and its harmonic at all. Measured + // over the corpus: all 36 firings at n <= 4 read 1.1 to 51.3 points BELOW the + // chance occupancy (n-1)/n, so an index-n sub-lattice really exists in each; all + // 8 above it read within 7.8 points OF chance, so none does. The largest n at + // which a real sub-lattice was ever seen is 3. + const bool integer_multiple = nearest >= 2.0 && nearest <= 4.0 + && std::abs(ratio - nearest) < 0.15; + + // Ask the data which of the pair is right instead of assuming the larger is + // spurious: which cell accounts for more of the validation frames' spots (see + // IndexedMillerIndices)? The measurement is made on their FULL spot lists, which + // reach far deeper into each frame's intensity distribution than the first pass's + // own accumulation cap - and a superstructure layer is faintest exactly where that + // cap cuts. + if (tied && !clearly_more && integer_multiple) { + const bool cand_is_larger = vol > bp.vol; + const RotationIndexerResult &larger = cand_is_larger ? result : *bp.result; + const RotationIndexerResult &smaller = cand_is_larger ? *bp.result : result; + const std::vector cloud = ValidationCloud(larger); + const auto hkl = IndexedMillerIndices(larger.lattice, cloud, index_tol); + const size_t n_smaller = IndexedMillerIndices(smaller.lattice, + ValidationCloud(smaller), + index_tol).size(); + // The added class on the PRIMITIVE basis: a centred basis has absences of its + // own, which would read as an empty class. Both cells are counted on their + // Niggli-reduced bases, because a tolerance on fractional Miller indices depends + // on the basis, and a skewed one makes the sub-lattice miss spots it contains. + const CrystalLattice larger_p = larger.lattice.ToPrimitive(larger.search_result.centering) + .NiggliReduce(); + const auto hkl_p = IndexedMillerIndices(larger_p, cloud, index_tol); + const auto ev = AddedClassOccupancy(hkl_p, static_cast(nearest)); + const CrystalLattice sub = SubLattice(larger_p, ev, static_cast(nearest)).NiggliReduce(); + const size_t n_sub = IndexedMillerIndices(sub, cloud, index_tol).size(); + const bool larger_is_real = hkl.size() > n_smaller && hkl_p.size() > n_sub; + logger.Info("Axis-harmonic arbiter: '{}' ({:.0f} A^3) vs '{}' ({:.0f} A^3), {:.2f}x - " + "the larger cell accounts for {} validation spots against the smaller " + "cell's {}; on its own geometry {} against its ({},{},{}) sub-lattice's " + "{}, {:.1f}% of them off it, so the class it adds is {}", + name, vol, bp.name, bp.vol, ratio, hkl.size(), n_smaller, + hkl_p.size(), ev.u, ev.v, ev.w, n_sub, 100.0 * ev.occupancy, + larger_is_real ? "occupied enough to pay for the longer axis - the " + "larger cell is real" + : "not occupied enough to pay for the longer axis - the " + "larger cell is a multiple"); + // Adopt the candidate exactly when its role matches the verdict. + integer_subcell = (cand_is_larger == larger_is_real); + // Only a win for the larger cell settles the pair. Declining it is a fall back + // to the default, and the warning below - which names the beam centre - is then + // still the most useful thing the run can say. It is the right diagnosis in + // particular for the failure this arbiter cannot see: a centre error along the + // spindle translates the derotated cloud rigidly by HALF a lattice spacing, and + // a half-offset lattice is indexable only on a doubled axis. Such a cell fits no + // index-n sublattice at all, so its added class reads the chance value (n-1)/n - + // the largest the occupancy can be - and it accounts for far fewer spots than + // the true cell, which is what this test sees it by. + harmonic_settled = larger_is_real; + } + // sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by + // sqrt(3) is the other harmonic seen on real data and its volume ratio is not a + // whole number. This is a REPORT and not a decision - it is deliberately taken + // over the whole pair, `clearly_more` included, because the harmonic is often the + // one that indexes more frames (see below). + if (integer_multiple || std::abs(ratio - std::sqrt(3.0)) < 0.15) + harmonic_ratio = ratio; + } + + // Two candidates whose primitive volumes differ by a small integer or by sqrt(3) are + // one lattice and its axis harmonic, and what decides between them is the BEAM CENTRE. + // A centre error along the spindle translates the derotated cloud rigidly, which the + // FFT amplitude cannot see at all, so the peaks stay sharp and the indexer commits - + // with confidence - to a sub-multiple. Measured on real data: a 0.12 px change of + // centre, 0.03x of what the J0 law asks across the spindle, is the whole difference + // between the deposited cell and a halved axis, and the halved one is the answer that + // indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate + // this pair, and until this line nothing in the run said what the pair meant. + // ... unless the arbiter above established the larger cell on the data, which is + // exactly the question this warning says the run cannot answer. + if (harmonic_ratio > 0.0 && !harmonic_settled) + logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) " + "differ by {:.2f}x - one is an axis harmonic of the other. Which of " + "the two is right is settled by the beam centre to a fraction of a " + "pixel, and the harmonic is often the one that indexes MORE frames, " + "so the frame counts above cannot decide it. Check the beam centre", + name, vol, bp.name, bp.vol, harmonic_ratio); + + if (!bp.result.has_value() || clearly_more || integer_subcell) { + if (integer_subcell && vol < bp.vol) + logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " + "equal frame count - adopting the smaller true cell", name, vol, + std::round(bp.vol / vol), bp.name, bp.vol); + else if (integer_subcell) + logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x supercell of '{}' (vol {:.0f}) at " + "equal frame count and accounts for more of the spots - adopting the " + "larger true cell", name, vol, std::round(vol / bp.vol), bp.name, bp.vol); + bp.score = score; + bp.vol = vol; + bp.name = name; + bp.result = std::move(result); + } + } + return bp; +} + +// Feed the schemes, index them concurrently, and keep the best (ScoreSchemes). +// alongside: started once the schemes are fed and before they are indexed, to run beside them. +auto Rugnux::FirstPassRun::PickBest(IndexerThreadPool &pool, IndexAndRefine &idx, + const std::function &alongside) -> FirstPass { + const SchemeIndexers ris = FeedSchemes(pool); + if (alongside) + alongside(); + RunIndexing(ris); + return ScoreSchemes(ris, idx); +} + +void Rugnux::FirstPassRun::ApplyState(const FirstPassState &st) { + const std::vector r = st.rings ? rx.powder_.rings_q_recipA : std::vector{}; + rx.config_.spot_finding.measured_ring_q_recipA = r; + validation_settings.measured_ring_q_recipA = r; + rx.config_.spot_finding.high_resolution_limit = st.d_min; + validation_settings.high_resolution_limit = st.d_min; + first_pass_spots_per_image = st.seed; + spot_cache.clear(); // the per-image spot budget is spent again under these settings +} + +// Powder-ring rescue - and, with it, the lean seed the rotation first pass never had. +// +// The pre-scan has already measured this sample's powder rings (see PreScan): a crystalline +// phase other than the crystal, whose reflections the indexer is otherwise handed as if they +// were this crystal's. Hexagonal ice is the only such phase whose rings can be named in +// advance, and the spot finder flags it from that fixed list; everything else - a shower of +// microcrystals, a salt out of the cryoprotectant, a phase nobody expected - has to be +// measured, and until it is, every one of its spots goes into the FFT. Measured on one such +// crystal: 876 spots a frame against the 80 a clean crystal on the same beamline gives, 70 % +// of them on rings, no scheme indexing a single validation frame, and the FFT taking a 10-30 A +// cell out of the shells. +// +// Setting the rings aside is half of it. The other half is HOW DEEP into each frame's spot +// list the first pass reads. It has always read 200 spots an image; the per-frame indexer, on +// the same problem, escalates 30 -> 80 -> the whole list and keeps the seed a lattice explains +// the largest fraction of (IndexAndRefine::DetermineLatticeAndSymmetry), because on a flooded +// frame a lean high-quality seed finds the lattice where the full list derails the search. The +// rotation pass gets the same escalation here. It is the general statement of what a +// resolution cut would do on this one crystal - the contaminant and the noise sit at the deep +// end of the intensity-ordered list wherever in resolution they happen to lie - and it needs +// nothing measured about where they are. +// +// It runs LAST: after every other rescue and after the pass-1 fallback above, on a run that +// would otherwise stop with no lattice at all. Measured over 100 datasets: placed earlier it +// pre-empted the rescues that already work - a crystal the beam-centre ladder fixes at 59/60 +// and one the pass-1 fallback fixes at 60/60 both came out WORSE, because this ladder had +// already replaced the standing pass those rescues start from. A last resort has to be last. +// +// The two compose: FilterSpotsByCount has already ordered each image's spots non-ring first +// and strongest first, so a lean cap takes the crystal's strongest reflections and nothing +// else. Decided late, the way the axis sign above is: index each rung, keep the one that +// indexes the most validation frames. A run that already indexes never reaches here, and one +// that does cannot lose by trying - where the rings are the crystal's own rows, or where the +// crystal needs its deep spots to index at all, the pass that already ran stands. +// +// Called from two places, because the two questions it answers are coupled: a centre that is +// wrong and a spot list that is too deep each hide the other, exactly as the centre and the +// axis sign do above. `score_before` is what the caller's standing pass indexes and `min_gain` +// is what a rung must beat it by; nothing is changed and the settings are put back when no +// rung does. +// +// min_score is an ABSOLUTE bar beside the gain. A caller that takes a rung only when it +// indexes a majority must not be handed a state the ladder committed on a lesser rung: +// measured, the ladder adopted at 29/60, the caller rejected it against a majority of 30, and +// nothing put the settings back - the beam-centre walk that followed then scored every +// hypothesis on the leaked lean spot lists. +auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int min_score) + -> std::optional { + if (rx.cancelled_) + return std::nullopt; + + // The seeds, leanest first. The same ladder the per-frame indexer uses, one rung deeper + // because a rotation cloud pools many frames and can afford more. + // The last rung is EVERY spot of each frame: SIZE_MAX, not 0 - RotationIndexer reads 0 as + // "use my default", which is 200, so a 0 here made the deepest rung a duplicate of the + // 200 one and the log line "the whole spot list" a lie. + constexpr std::array SEEDS = {30, 80, 200, std::numeric_limits::max()}; + struct Rung { bool rings; size_t seed; std::optional d_min; }; + // How much of each frame a rung let the first pass see: more spots an image, and further + // in resolution, is more evidence. Rings set aside is not part of it - that is a + // measurement of the sample, not a depth. + const auto Deeper = [](const Rung &a, const Rung &b) { + if (a.seed != b.seed) return a.seed > b.seed; + const float da = a.d_min.value_or(0.0f), db = b.d_min.value_or(0.0f); + return da < db; // a smaller d_min reaches further; 0 means the detector edge + }; + // The file's OWN limit, captured before any rung committed anything. Reading it out of + // rx.config_ here would read back the previous rung's floor once a pass has adopted one, and + // the restore path would then "put back" that instead of what the file asked for. + const auto file_limit = first_pass_file_d_min; + std::vector rungs; + // Where a floor is worth trying at all - the rings merged somewhere inside the measured + // range, so past that point nothing can be told from the crystal and a pass that reads + // there is reading the contaminant. + std::vector> floors = {file_limit}; + for (const float d : rx.prescan_spot_quantile_d_A_) + if (!file_limit || d > *file_limit) floors.emplace_back(d); + // Floors outer, seeds inner, so rungs that share the finder-facing settings (rings, + // d_min) sit adjacent: the seed is an indexer-side cap over the SAME found spots + // (MaxSpotsPerImage over the cached list), so the loop below spot-finds once per floor + // instead of once per rung. Adoption is decided over all rungs at once, so the visit + // order carries no meaning. + for (const auto &d_min : floors) + for (const size_t seed : SEEDS) { + if (!rx.powder_.rings_q_recipA.empty()) + rungs.push_back({true, seed, d_min}); + if (rx.powder_.rings_q_recipA.empty() && seed != SEEDS.back()) + rungs.push_back({false, seed, d_min}); + } + + // A new floor (or ring handling) changes what the finder returns, so the spot lists are + // spent and found again; a new seed alone does not - it caps how much of the SAME lists + // the indexer reads - and re-finding 60 validation frames per rung was most of what a + // hopeless run's minutes were spent on. + const auto enter_rung = [&](size_t i) { + if (i == 0 || rungs[i].rings != rungs[i - 1].rings || rungs[i].d_min != rungs[i - 1].d_min) + ApplyState({rungs[i].rings, rungs[i].seed, rungs[i].d_min}); + else + first_pass_spots_per_image = rungs[i].seed; + }; + // Every rung is fed first and all of them are indexed at once: the rungs do not depend on + // one another, and indexing them one after the other left most of the machine idle for + // most of a hopeless run. They are then scored in order, each under its own settings again. + std::vector rung_ris(rungs.size()); + std::vector all_ris; + for (size_t i = 0; i < rungs.size() && !rx.cancelled_; i++) { + enter_rung(i); + rung_ris[i] = FeedSchemes(*indexer_pool); + for (const auto &ri : rung_ris[i]) + all_ris.push_back(ri.get()); + } + RunIndexing(all_ris); + + // Every rung's score, so the choice below can be made on all of them at once. + std::vector rung_score(rungs.size(), -1); + std::vector rung_pass(rungs.size()); + int won = -1; + for (size_t i = 0; i < rungs.size() && !rx.cancelled_; i++) { + enter_rung(i); + FirstPass alt = ScoreSchemes(rung_ris[i], *indexer); + logger.Info(" first-pass rung: {} spots/image, rings {}, stop at {} -> {}/{} frames{}", + rungs[i].seed == std::numeric_limits::max() + ? std::string("all") : std::to_string(rungs[i].seed), + rungs[i].rings ? "aside" : "kept", + rungs[i].d_min ? fmt::format("{:.2f} A", *rungs[i].d_min) + : std::string("the detector edge"), + alt.result.has_value() ? alt.score : -1, + static_cast(validation.size()), + alt.result.has_value() + ? fmt::format(", {}{} {:.1f} {:.1f} {:.1f} {:.1f} {:.1f} {:.1f}", + alt.result->search_result.centering, + gemmi::crystal_system_str(alt.result->search_result.system), + alt.result->search_result.conventional.GetUnitCell().a, + alt.result->search_result.conventional.GetUnitCell().b, + alt.result->search_result.conventional.GetUnitCell().c, + alt.result->search_result.conventional.GetUnitCell().alpha, + alt.result->search_result.conventional.GetUnitCell().beta, + alt.result->search_result.conventional.GetUnitCell().gamma) + : std::string()); + // A rung has to WIN by the count the run itself keeps a lattice on - a sixth of the + // validation frames, the bar below - not merely come out ahead. Three frames out of + // sixty is the count's noise, and adopting on it would hand the rescues that follow // (the beam centre, the long axis) a spot list cut for no reason. Measured: a marginal - // crystal that indexed 26/60 went to 29/60 on a rung and must not take it. - // - // Among the rungs that do win it is the highest count that is kept, and an exact tie goes - // to the rung that saw the MOST of each frame. Ties are common - on a pattern where the - // contaminant sits past a resolution, every seed behind that resolution returns the same - // lattice and the same count - and where they happen the rung that read deepest rests on - // the most evidence. Only an exact tie, not a band: the count is what tells a lattice that - // works from one that does not, and a rung three frames behind is behind. - for (size_t i = 0; i < rungs.size(); i++) { - if (rung_score[i] < 0 || rung_score[i] - score_before < min_gain - || rung_score[i] < min_score) - continue; - if (won < 0 || rung_score[i] > rung_score[won] - || (rung_score[i] == rung_score[won] && Deeper(rungs[i], rungs[won]))) - won = static_cast(i); + // crystal that indexed 26/60 went to 29/60 on a rung and must not take it, where the + // crystals this exists for go from 0/60 to 24 and to 32/60. + if (alt.result.has_value()) { + rung_score[i] = alt.score; + rung_pass[i] = std::move(alt); } + } - // ...and it must not be an axis SUB-MULTIPLE of a cell another rung found. The frame - // count the choice above is made on cannot arbitrate an axis harmonic - a cell twice as - // long has to place every spot twice as accurately to score the same, which the scheme - // comparison above says in the same words - and here the two biases compound, because - // the leanest and shallowest rung is both the one that scores highest on a harmonic and - // the one a smaller cell is easiest on. Measured on one crystal: the "30 spots/image, - // stop at 4.91 A" rung indexed 48/60 on a HALVED c axis while the deeper rungs found 2x, - // 3x and 4x of the true one; the ladder took the halved cell, and the run ended in a - // different Bravais class on a different cell from the one it reaches with no ladder at - // all - the ladder had pre-empted the beam-centre walk that finds the crystal's own. - // - // So put the winner through the same axis-harmonic arbiter the scheme comparison uses: - // against every rung whose primitive volume is a near-integer multiple of the winner's, - // which cell accounts for more of the validation frames' SPOTS. A rung that loses that - // question is a sub-multiple, and the ladder then adopts NOTHING. It is a last resort, - // and a last resort that cannot tell an axis from its multiple has not rescued anything - - // promoting the rival instead is no better, since on the measured crystal the rival is a - // multiple of the true axis in its own right. - // - // Judged on the spot list the run STANDS on rather than on either rung's cut: the two - // cells must be counted against the same spots, and the deepest list is where the class - // an axis multiple adds is visible at all (see IndexedMillerIndices). - if (won >= 0 && rung_pass[won].vol > 1.0) { - apply_state(committed); - prefetch_spots(validation); - const double won_vol = rung_pass[won].vol; - const size_t won_spots = IndexedMillerIndices(rung_pass[won].result->lattice, - validation_cloud(*rung_pass[won].result), - index_tol).size(); - for (size_t i = 0; i < rungs.size() && won >= 0; i++) { - // Only a rung the caller would have taken can demote one: the winner is chosen - // among those, and a rung that indexes too little to be adopted is too little - // evidence to overturn one either. - if (static_cast(i) == won || rung_score[i] < 0 - || rung_score[i] - score_before < min_gain || rung_score[i] < min_score) - continue; - const double ratio = rung_pass[i].vol / won_vol; - const double nearest = std::round(ratio); - // Only a LARGER rival can demote, and only in the window the arbiter is - // calibrated over (see the scheme comparison: every firing at n <= 4 read below - // the chance occupancy, none above it did). - if (nearest < 2.0 || nearest > 4.0 || std::abs(ratio - nearest) >= 0.15) - continue; - const auto hkl = IndexedMillerIndices(rung_pass[i].result->lattice, - validation_cloud(*rung_pass[i].result), - index_tol); - const auto ev = AddedClassOccupancy(hkl, static_cast(nearest)); - logger.Info("First-pass ladder: the winning rung's cell ({:.0f} A^3) and the " - "rung at {}/{} frames ({:.0f} A^3) differ by {:.2f}x - one is an " - "axis harmonic of the other. The larger accounts for {} validation " - "spots against the winner's {}, and {:.1f}% of them lie off its " - "({},{},{}) sub-lattice", - won_vol, rung_score[i], static_cast(validation.size()), - rung_pass[i].vol, ratio, hkl.size(), won_spots, - 100.0 * ev.occupancy, ev.u, ev.v, ev.w); - if (hkl.size() > won_spots) { - logger.Warning("First-pass ladder: the rung that indexes the most frames does " - "so on a {:.0f}x sub-multiple of an axis another rung found, " - "which the frame count cannot see - the ladder adopts nothing " - "and the run keeps the pass it already stands on", - nearest); - won = -1; - } + // The rung to keep. A rung has to WIN by the count the run itself keeps a lattice on - a + // sixth of the validation frames, the bar below - not merely come out ahead: three frames + // out of sixty is the count's noise, and adopting on it would hand the rescues that follow + // (the beam centre, the long axis) a spot list cut for no reason. Measured: a marginal + // crystal that indexed 26/60 went to 29/60 on a rung and must not take it. + // + // Among the rungs that do win it is the highest count that is kept, and an exact tie goes + // to the rung that saw the MOST of each frame. Ties are common - on a pattern where the + // contaminant sits past a resolution, every seed behind that resolution returns the same + // lattice and the same count - and where they happen the rung that read deepest rests on + // the most evidence. Only an exact tie, not a band: the count is what tells a lattice that + // works from one that does not, and a rung three frames behind is behind. + for (size_t i = 0; i < rungs.size(); i++) { + if (rung_score[i] < 0 || rung_score[i] - score_before < min_gain + || rung_score[i] < min_score) + continue; + if (won < 0 || rung_score[i] > rung_score[won] + || (rung_score[i] == rung_score[won] && Deeper(rungs[i], rungs[won]))) + won = static_cast(i); + } + + // ...and it must not be an axis SUB-MULTIPLE of a cell another rung found. The frame + // count the choice above is made on cannot arbitrate an axis harmonic - a cell twice as + // long has to place every spot twice as accurately to score the same, which the scheme + // comparison above says in the same words - and here the two biases compound, because + // the leanest and shallowest rung is both the one that scores highest on a harmonic and + // the one a smaller cell is easiest on. Measured on one crystal: the "30 spots/image, + // stop at 4.91 A" rung indexed 48/60 on a HALVED c axis while the deeper rungs found 2x, + // 3x and 4x of the true one; the ladder took the halved cell, and the run ended in a + // different Bravais class on a different cell from the one it reaches with no ladder at + // all - the ladder had pre-empted the beam-centre walk that finds the crystal's own. + // + // So put the winner through the same axis-harmonic arbiter the scheme comparison uses: + // against every rung whose primitive volume is a near-integer multiple of the winner's, + // which cell accounts for more of the validation frames' SPOTS. A rung that loses that + // question is a sub-multiple, and the ladder then adopts NOTHING. It is a last resort, + // and a last resort that cannot tell an axis from its multiple has not rescued anything - + // promoting the rival instead is no better, since on the measured crystal the rival is a + // multiple of the true axis in its own right. + // + // Judged on the spot list the run STANDS on rather than on either rung's cut: the two + // cells must be counted against the same spots, and the deepest list is where the class + // an axis multiple adds is visible at all (see IndexedMillerIndices). + if (won >= 0 && rung_pass[won].vol > 1.0) { + ApplyState(committed); + PrefetchSpots(validation); + const double won_vol = rung_pass[won].vol; + const size_t won_spots = IndexedMillerIndices(rung_pass[won].result->lattice, + ValidationCloud(*rung_pass[won].result), + index_tol).size(); + for (size_t i = 0; i < rungs.size() && won >= 0; i++) { + // Only a rung the caller would have taken can demote one: the winner is chosen + // among those, and a rung that indexes too little to be adopted is too little + // evidence to overturn one either. + if (static_cast(i) == won || rung_score[i] < 0 + || rung_score[i] - score_before < min_gain || rung_score[i] < min_score) + continue; + const double ratio = rung_pass[i].vol / won_vol; + const double nearest = std::round(ratio); + // Only a LARGER rival can demote, and only in the window the arbiter is + // calibrated over (see the scheme comparison: every firing at n <= 4 read below + // the chance occupancy, none above it did). + if (nearest < 2.0 || nearest > 4.0 || std::abs(ratio - nearest) >= 0.15) + continue; + const auto hkl = IndexedMillerIndices(rung_pass[i].result->lattice, + ValidationCloud(*rung_pass[i].result), + index_tol); + const auto ev = AddedClassOccupancy(hkl, static_cast(nearest)); + logger.Info("First-pass ladder: the winning rung's cell ({:.0f} A^3) and the " + "rung at {}/{} frames ({:.0f} A^3) differ by {:.2f}x - one is an " + "axis harmonic of the other. The larger accounts for {} validation " + "spots against the winner's {}, and {:.1f}% of them lie off its " + "({},{},{}) sub-lattice", + won_vol, rung_score[i], static_cast(validation.size()), + rung_pass[i].vol, ratio, hkl.size(), won_spots, + 100.0 * ev.occupancy, ev.u, ev.v, ev.w); + if (hkl.size() > won_spots) { + logger.Warning("First-pass ladder: the rung that indexes the most frames does " + "so on a {:.0f}x sub-multiple of an axis another rung found, " + "which the frame count cannot see - the ladder adopts nothing " + "and the run keeps the pass it already stands on", + nearest); + won = -1; } } + } - if (won >= 0) { + if (won >= 0) { // The rung says which spots let the FFT FIND the lattice; how deep the run then READS each // frame is asked of the lattice it found. Which floor's FFT happens to succeed is close to // a coin toss on a hard pattern - one crystal took 1.13 A on one sample of frames and @@ -1297,8 +1269,8 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { return a.value_or(0.0f) < b.value_or(0.0f); }); for (const auto &d : deeper) { - apply_state({rungs[won].rings, rungs[won].seed, d}); - const int n = count_indexed(*indexer, *rung_pass[won].result); + ApplyState({rungs[won].rings, rungs[won].seed, d}); + const int n = CountIndexed(*indexer, *rung_pass[won].result); logger.Info("First-pass depth: the winning rung's lattice indexes {}/{} validation " "frames reading to {}", n, static_cast(validation.size()), d ? fmt::format("{:.2f} A", *d) : std::string("the detector edge")); @@ -1309,14 +1281,14 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { } } committed = {rungs[won].rings, rungs[won].seed, read_to}; - apply_state(committed); + ApplyState(committed); logger.Info("First-pass seed: {}{}{} takes it from {}/{} validation frames to {}/{}", rungs[won].seed == std::numeric_limits::max() ? std::string("the whole spot list") : fmt::format("the strongest {} spots an image", rungs[won].seed), rungs[won].rings ? fmt::format(", with this sample's {} measured powder rings set aside " - "for the whole run", powder_.rings_q_recipA.size()) + "for the whole run", rx.powder_.rings_q_recipA.size()) : std::string(), read_to && read_to != file_limit ? fmt::format(", stopping spot finding at {:.2f} A", *read_to) @@ -1325,1472 +1297,117 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { rung_pass[won].score, static_cast(validation.size())); return rung_pass[won]; } - apply_state(committed); // back to what this run stands on, not to the defaults + ApplyState(committed); // back to what this run stands on, not to the defaults return std::nullopt; - }; - - // Find the spots of everything the first pass is certain to want in one go - the schemes' - // opening frames and the validation set. A batch of workers costs an analysis engine each, - // which is far more than a frame is to analyse, so one batch for the lot beats one per scheme. +} + +// Find the spots of everything the first pass is certain to want in one go - the schemes' +// opening frames and the validation set. A batch of workers costs an analysis engine each, +// which is far more than a frame is to analyse, so one batch for the lot beats one per scheme. +void Rugnux::FirstPassRun::PrefetchCertainFrames() { + std::vector upfront = spread_ordinals; + upfront.insert(upfront.end(), validation.begin(), validation.end()); + for (int i = 0; i < std::min(SPOT_PREFETCH_CHUNK, early_images); i++) + upfront.push_back(i); + std::sort(upfront.begin(), upfront.end()); + upfront.erase(std::unique(upfront.begin(), upfront.end()), upfront.end()); + PrefetchSpots(upfront); +} + +IndexingSettings Rugnux::FirstPassRun::ShortAxisSettings(const DiffractionExperiment &x) const { + auto settings = x.GetIndexingSettings(); + settings.FFT_MinUnitCell_A(IndexingSettings::fft_min_unit_cell_limit_A); + settings.IndexingThreads(2); + settings.RefineThreads(std::clamp(rx.config_.nthreads / 2, 1, MAX_REFINE_THREADS)); + return settings; +} + +// The beam-centre check's pass at the measured centre, started beside the file centre's (see +// centre_spec). +void Rugnux::FirstPassRun::StartCentreSpeculation() { + if (get_gpu_count() < 2) + return; + rx.JoinBeamCenterCapture(); + if (rx.cancelled_ || !rx.config_.beam_center_check || !rx.background_center_) + return; + auto &cs = centre_spec.emplace(); + cs.experiment = std::make_unique(rx.experiment_); + cs.experiment->BeamX_pxl(rx.background_center_->beam_x_pxl).BeamY_pxl(rx.background_center_->beam_y_pxl); + cs.mask = std::make_unique(rx.pixel_mask_); + cs.mapping = std::make_unique(*cs.experiment, *cs.mask, *rx.azint_geometry_); + cs.experiment_key = rx.ExperimentKey(*cs.experiment); + cs.key = rx.SpotFindingKey(*cs.experiment, *cs.mapping); + cs.spots_per_image = first_pass_spots_per_image; { - std::vector upfront = spread_ordinals; - upfront.insert(upfront.end(), validation.begin(), validation.end()); - for (int i = 0; i < std::min(SPOT_PREFETCH_CHUNK, early_images); i++) - upfront.push_back(i); - std::sort(upfront.begin(), upfront.end()); - upfront.erase(std::unique(upfront.begin(), upfront.end()), upfront.end()); - prefetch_spots(upfront); + // What an earlier pass already found at this centre under these settings. + std::lock_guard lock(rx.first_pass_spots_->m); + for (const auto &[k, spots] : rx.first_pass_spots_->entries) + if (k == cs.key) + cs.spots = spots; } - - // The short-axis pass below indexes the same schemes at a lower floor, after the rescues in - // between - and on nearly every run they change nothing it reads. So it is started here, on a - // copy of the experiment, beside the standard indexing it would otherwise wait for, and taken - // at its place only if what it read is still what the run has there (ExperimentKey, the spot - // list's mapping and budget); otherwise it is run again there as before. - const auto short_axis_settings = [&](const DiffractionExperiment &x) { - auto settings = x.GetIndexingSettings(); - settings.FFT_MinUnitCell_A(IndexingSettings::fft_min_unit_cell_limit_A); - settings.IndexingThreads(2); - settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); - return settings; - }; - struct ShortAxisSpeculation { - std::vector key; - const AzimuthalIntegrationMapping *mapping = nullptr; - size_t spots_per_image = 0; - std::unique_ptr experiment; - std::unique_ptr pool; - std::vector> sris; - std::vector> futures; // last, so it is waited for before the rest goes - }; - std::optional short_spec; - // The beam-centre check below runs a second first pass at the centre measured from the - // background, and it reads the first one's answer only after that pass is over. So the spots at - // the measured centre are found, and both schemes fed and indexed, here beside the file centre's - // pass - on a copy of the experiment at that centre, with workers, a mapping and a spot list of - // their own - and taken at the check's place only if what they read is still what the run has - // there (the experiment with its rotation angles, SpotFindingKey at the measured centre, the spot - // budget); otherwise that pass is run there as before. The schemes are scored there in either - // case, on the run's own indexer. - // - // Only with a second GPU: on one card the two passes share the card and the reader, and measured, - // running them side by side saved nothing. There the check runs its pass after the first, as it - // always did. With two or more its spot finding spreads over all the cards, as the first pass's - // does: the file centre's spot finding is over long before this pass's, so this one is what the - // first pass waits on. - struct CentreSpeculation { - std::vector experiment_key; - std::vector key; - size_t spots_per_image = 0; - std::unique_ptr experiment; - std::unique_ptr mask; - std::unique_ptr mapping; - std::map> spots; - SchemeIndexers ris; - std::future done; // last, so it is waited for before the rest goes - }; - std::optional centre_spec; - const auto start_centre_speculation = [&] { - if (get_gpu_count() < 2) - return; - JoinBeamCenterCapture(); - if (cancelled_ || !config_.beam_center_check || !background_center_) - return; - auto &cs = centre_spec.emplace(); - cs.experiment = std::make_unique(experiment_); - cs.experiment->BeamX_pxl(background_center_->beam_x_pxl).BeamY_pxl(background_center_->beam_y_pxl); - cs.mask = std::make_unique(pixel_mask_); - cs.mapping = std::make_unique(*cs.experiment, *cs.mask, *azint_geometry_); - cs.experiment_key = ExperimentKey(*cs.experiment); - cs.key = SpotFindingKey(*cs.experiment, *cs.mapping); - cs.spots_per_image = first_pass_spots_per_image; - { - // What an earlier pass already found at this centre under these settings. - std::lock_guard lock(first_pass_spots_->m); - for (const auto &[k, spots] : first_pass_spots_->entries) - if (k == cs.key) - cs.spots = spots; - } - // At this pass's angles, as prefetch_spots takes them. - const DiffractionGeometry spot_geometry = cs.experiment->GetDiffractionGeometry(); - for (auto &[ordinal, spots] : cs.spots) - StampSpotAngles(spots, spot_geometry); - cs.done = std::async(std::launch::async, [this, &cs, &schemes, &validation, &indexer_pool, &run_indexing, - settings = config_.spot_finding, dataset, start_image, - nworkers = engines.size()] { - SpotFindingWorkers workers(reader_, *cs.experiment, *cs.mapping, *cs.mask, settings, *dataset, - start_image, config_.stride, nworkers); - const auto find = [&](const std::vector &ordinals) { - std::vector wanted; - for (const int o : ordinals) - if (!cs.spots.contains(o)) - wanted.push_back(o); - auto found = workers.Find(wanted); - for (size_t i = 0; i < wanted.size(); i++) - cs.spots.emplace(wanted[i], std::move(found[i])); - }; - // As feed_scheme, at the measured centre. - const auto gonio = cs.experiment->GetGoniometer(); - for (const auto &[name, ordinals] : schemes) { - auto ri = std::make_unique(*cs.experiment, *indexer_pool); - ri->MaxSpotsPerImage(cs.spots_per_image); - ri->VerifyMemo(config_.verify_first_pass_memo); - for (size_t i = 0; i < ordinals.size(); i++) { - if (cancelled_ || ri->AccumulationFull()) - break; - if (!cs.spots.contains(ordinals[i])) - find({ordinals.begin() + i, - ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)}); - std::optional angle; - if (gonio) - angle = gonio->GetAngle_deg(static_cast(ordinals[i])) + gonio->GetWedge_deg() / 2.0f; - ri->ProcessImage(ordinals[i], cs.spots.at(ordinals[i]), angle); - } - cs.ris.push_back(std::move(ri)); + // At this pass's angles, as PrefetchSpots takes them. + const DiffractionGeometry spot_geometry = cs.experiment->GetDiffractionGeometry(); + for (auto &[ordinal, spots] : cs.spots) + StampSpotAngles(spots, spot_geometry); + cs.done = std::async(std::launch::async, [this, &cs, settings = rx.config_.spot_finding, dataset = dataset, + start_image = start_image, nworkers = engines.size()] { + SpotFindingWorkers workers(rx.reader_, *cs.experiment, *cs.mapping, *cs.mask, settings, *dataset, + start_image, rx.config_.stride, nworkers); + const auto find = [&](const std::vector &ordinals) { + std::vector wanted; + for (const int o : ordinals) + if (!cs.spots.contains(o)) + wanted.push_back(o); + auto found = workers.Find(wanted); + for (size_t i = 0; i < wanted.size(); i++) + cs.spots.emplace(wanted[i], std::move(found[i])); + }; + // As FeedScheme, at the measured centre. + const auto gonio = cs.experiment->GetGoniometer(); + for (const auto &[name, ordinals] : schemes) { + auto ri = std::make_unique(*cs.experiment, *indexer_pool); + ri->MaxSpotsPerImage(cs.spots_per_image); + ri->VerifyMemo(rx.config_.verify_first_pass_memo); + for (size_t i = 0; i < ordinals.size(); i++) { + if (rx.cancelled_ || ri->AccumulationFull()) + break; + if (!cs.spots.contains(ordinals[i])) + find({ordinals.begin() + i, + ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)}); + std::optional angle; + if (gonio) + angle = gonio->GetAngle_deg(static_cast(ordinals[i])) + gonio->GetWedge_deg() / 2.0f; + ri->ProcessImage(ordinals[i], cs.spots.at(ordinals[i]), angle); } - find(validation); - std::vector rp; - for (const auto &ri : cs.ris) - rp.push_back(ri.get()); - run_indexing(rp); - }); - }; - FirstPass best = pick_best(*indexer_pool, *indexer, [&] { - start_centre_speculation(); - if (cancelled_ || experiment_.GetUnitCell().has_value()) - return; - auto &sp = short_spec.emplace(); - sp.key = ExperimentKey(); - sp.mapping = spot_mapping; - sp.spots_per_image = first_pass_spots_per_image; - sp.experiment = std::make_unique(experiment_); - sp.pool = std::make_unique(short_axis_settings(*sp.experiment), - IndexerConstruction::OnFirstUse); - for (const auto &[name, ordinals] : schemes) - sp.sris.push_back(feed_scheme(*sp.pool, ordinals, sp.experiment.get())); - for (auto &ri : sp.sris) { - RotationIndexer *rp = ri.get(); - sp.futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); + cs.ris.push_back(std::move(ri)); } + find(validation); + RunIndexing(cs.ris); }); - - // Rotation-axis SIGN rescue. The sign is a convention the input cannot settle: a PILATUS - // miniCBF header carries only "Oscillation_axis omega" and no direction at all, and facilities - // disagree about the sign of the NXmx vector (measured across ten of them: most want the - // stored [-1,0,0], APS 19-ID and SPring-8 BL41XU want the opposite). Nothing downstream can - // repair it either - the default spindle fit refines the axis DIRECTION by fractions of a - // degree and cannot flip it, that being a discrete change. - // - // It is, however, decidable from the data: the wrong sign does not index at all. Measured on - // one dataset by changing nothing but the sign: 100% of frames indexed one way, 0 of 60 the - // other. So try the other sign and keep whichever indexes more frames - the same - // validation-frame count the scheme choice above already uses, so no new metric and no new - // threshold. Runs before the long-axis rescue: with the sign wrong every candidate lattice is - // wrong, and the coarse pass below would spend itself on a mis-signed cloud. - // - // Only after a poor pass, so a correctly-signed file costs nothing; and spot finding is not - // repeated, the cache being keyed by image and the spots not depending on the axis. - if (!cancelled_ && best.score < 0.5 * static_cast(validation.size())) { - if (const auto gon = experiment_.GetGoniometer(); gon.has_value() && gon->IsScanning()) { - GoniometerAxis flipped = *gon; - flipped.Axis(-gon->GetAxis()); - experiment_.Goniometer(flipped); - FirstPass alt = pick_best(*indexer_pool, *indexer); - // More frames is not enough on its own while both counts sit at the floor: 1/60 against - // 0/60 is noise, and adopting it puts every rescue below on the wrong sign (measured: a - // run whose file sign, kept, went on to index found nothing once flipped on 1/60). So - // the flipped lattice must also beat chance on the pooled spots - the test the run puts - // its final lattice to - which a sparse pattern at the right sign passes on few frames - // (2/60 frames, 89% of the validation spots against 0% by chance). - const bool alt_real = alt.result.has_value() && alt.score > best.score - && beats_chance(pooled_evidence(*indexer, *alt.result)); - if (alt_real) { - logger.Info("Rotation axis sign: the file's axis indexes {}/{} validation frames, " - "the opposite sign {}/{} - adopting the opposite sign for this run", - best.score, static_cast(validation.size()), - alt.score, static_cast(validation.size())); - best = alt; - } else { - if (alt.result.has_value() && alt.score > best.score) - logger.Info("Rotation axis sign: the opposite sign indexes {}/{} validation frames " - "against the file's {}/{}, but its lattice does not beat chance on the " - "pooled spots - the file's sign is kept", - alt.score, static_cast(validation.size()), - best.score, static_cast(validation.size())); - experiment_.Goniometer(*gon); // the file was right; put it back - } - } +} + +// The short-axis pass, started beside the standard indexing (see short_spec). +void Rugnux::FirstPassRun::StartShortAxisSpeculation() { + if (rx.cancelled_ || rx.experiment_.GetUnitCell().has_value()) + return; + auto &sp = short_spec.emplace(); + sp.key = rx.ExperimentKey(); + sp.mapping = spot_mapping; + sp.spots_per_image = first_pass_spots_per_image; + sp.experiment = std::make_unique(rx.experiment_); + sp.pool = std::make_unique(ShortAxisSettings(*sp.experiment), + IndexerConstruction::OnFirstUse); + for (const auto &[name, ordinals] : schemes) + sp.sris.push_back(FeedScheme(*sp.pool, ordinals, sp.experiment.get())); + for (auto &ri : sp.sris) { + RotationIndexer *rp = ri.get(); + sp.futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); } +} - // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced - // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let - // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the - // many high-order reflections off along the fine axis, so the true cell scores worst on the raw - // cloud). Recover the true metric with a second search whose FFT unit-cell bound is widened past - // the default, then RE-INDEX at full resolution constrained by the cell that comes back, as a - // reference (the -C path): the reference filter drops the collapsed/supercell candidates and - // refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing - // crystals never pay for it and their result is untouched. Returns the re-index where it indexes - // more validation frames than the `standing` pass it was asked about. - auto max_axis = [](const RotationIndexerResult &r) { - const auto uc = r.lattice.GetUnitCell(); - return std::max({uc.a, uc.b, uc.c}); - }; - const auto long_axis_rescue = [&](const FirstPass &standing) -> std::optional { - if (cancelled_ || !standing.result.has_value() - || standing.score >= 0.5 * static_cast(validation.size())) - return std::nullopt; - auto coarse_settings = experiment_.GetIndexingSettings(); - // The resolution setting only sizes the FFT histogram - it does not filter which spots go - // into it - so this is here to pay for the widened bound below, not to make the long axis - // easier to pick out. Measured, coarsening does the opposite: the acceptance angle widens - // (0.159 -> 0.579 deg) but the true axis falls from rank 1996 to 16361 of 16384 directions. - coarse_settings.FFT_HighResolution_A(3.5f); - // The widened bound is what the rescue turns on, and it costs nothing: the histogram length is - // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already - // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at - // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover: - // a cell whose long axis is past the default 500 A bound is unrepresentable, so the - // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a - // well-indexing crystal never pays. Stills never reach here. - constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f; - if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A) - coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A); - coarse_settings.IndexingThreads(2); - coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); - IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); - // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- - // resolution per-frame validation would be low (the coarse cell is metrically right but - // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. - std::vector> cris; - cris.reserve(schemes.size()); - for (const auto &[name, ordinals] : schemes) { - if (cancelled_) break; - cris.push_back(feed_scheme(coarse_pool, ordinals)); - } - std::vector> cfut; - cfut.reserve(cris.size()); - for (auto &ri : cris) { - RotationIndexer *rp = ri.get(); - cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); - } - for (auto &f : cfut) f.get(); - std::optional coarse_ref; - for (auto &ri : cris) { - auto r = ri->GetLattice(); - if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref))) - coarse_ref = r; - } - // A recovered axis at least 30% longer than the standard pass found signals a rescued metric. - if (!cancelled_ && coarse_ref.has_value() - && max_axis(*coarse_ref) > 1.3 * max_axis(*standing.result)) { - logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); " - "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*standing.result)); - experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell()); - // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The - // standard pool was sized at the default bound, so re-using it here would search for a - // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace. - // count_indexed forces the lattice and never touches the pool, so the existing - // IndexAndRefine is reused; only the FFT needs the wider one. - auto wide_settings = experiment_.GetIndexingSettings(); - wide_settings.IndexingThreads(2); - wide_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS)); - RaiseFFTBoundForKnownCell(wide_settings, experiment_); - FirstPass constrained; - if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) { - IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse); - constrained = pick_best(wide_pool, *indexer); - } else { - constrained = pick_best(*indexer_pool, *indexer); - } - experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo - if (constrained.result.has_value() && constrained.score > standing.score) - return constrained; - } - return std::nullopt; - }; - // Set where the beam-centre check below has already asked the rescue at the centre it keeps. - bool long_axis_asked = false; - - // The background centre, tried as a SECOND FIRST PASS. - // - // The estimate itself is free - it is read off the projection --detect-beam-stop has already - // built - and a first pass re-uses the spots it has already found, so indexing the run a - // second time at the measured centre costs about one more first pass. - // That buys the comparison an on-failure search structurally cannot make: a centre wrong - // ALONG the spindle does not fail. It keeps 96-100 % of frames indexed and quietly returns a - // 2x, 3x or sqrt(3) axis harmonic, so no failure trigger ever fires on it. Indexing both - // centres and comparing the two ANSWERS is what can see it. - // - // What must NOT arbitrate the two is the frame count. Acceptance is a fractional-Miller test, - // so a cell twice as long has to place every spot twice as accurately to score the same: - // measured on real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a - // halved axis at 100.00 %, and the WRONG answer indexed better. A rule of the form "take the - // centre that indexes more" would therefore pick wrong in exactly the case this comparison - // exists for. So only what needs no arbiter is decided here: - // - // * the file's centre does not index a majority -> nothing defends it, so try the measured - // one (and the axis sign with it, see below). A run that indexes nothing has nothing to - // lose; where the measured one is no better, the file's centre is put back; - // * the file's centre indexes and the two agree on the lattice -> keep the file's centre and - // say they agree. This is the overwhelming majority, and saying it is worth the pass: it - // is a free statement that the header is good enough for this crystal; - // * they index and disagree -> say so, with both cells and their primitive-volume ratio. - // Which of an integer-related pair is real is not decided on the frame counts, for the - // reason above; it is decided, in ONE direction only, on the pooled spots - see the - // branch itself. Naming the disagreement is useful even where nothing is decided. - // - // Where the file's centre indexes nothing and the measured one does not either, the axis-sign - // question is asked AGAIN at the measured centre before giving up. The two unknowns are - // discrete and coupled, and the rescue above asked its question at a centre that may be 73 px - // wrong: on the two depositions here whose header is the geometric centre of the detector, - // each error hides the other, both runs sit at 0/60 until both are moved, and moving both - // takes them to 60/60 on the deposited lattice. Asking the pair costs one more first pass on - // a run that has already failed twice. - JoinBeamCenterCapture(); - if (!cancelled_ && config_.beam_center_check && background_center_) { - // Consumed, not re-read: a two-pass rotation run comes back through here with the geometry - // the first pass post-refined, and a pre-scan estimate of the FILE's centre says nothing - // about that one. The pre-scan does not re-measure it either - by then the beam-stop shadow - // is already in the mask, so the projection is not rebuilt. - const BeamCenterEstimate measured = *background_center_; - background_center_.reset(); - const float header_x = experiment_.GetBeamX_pxl(), header_y = experiment_.GetBeamY_pxl(); - const float measured_x = measured.beam_x_pxl; - const float measured_y = measured.beam_y_pxl; - const float moved = std::hypot(measured_x - header_x, measured_y - header_y); - // No gate on how significant the move is. It is tempting to skip the second pass when the - // two centres are within the fit's own noise - it would save it on a quarter of the runs - // here - but that is deciding at a threshold before the evidence is in, and the evidence - // is cheap. It would also gate away exactly the case worth asking about: measured on real - // data, a 0.12 px change of centre is the whole difference between the deposited cell and - // a halved axis, and 0.12 px is well inside the sigma of any fit that could gate it. - { - const int majority = static_cast(validation.size()) / 2; - try_beam_center(measured_x, measured_y); - FirstPass alt; - if (centre_spec) - centre_spec->done.get(); - if (centre_spec && centre_spec->experiment_key == ExperimentKey() - && centre_spec->key == SpotFindingKey(*spot_mapping) - && centre_spec->spots_per_image == first_pass_spots_per_image) { - // What pick_best would have found and fed here: hand the spots to the run's cache and - // its memo as prefetch_spots does, and score the schemes on the run's own indexer. - spot_cache = centre_spec->spots; - { - std::lock_guard lock(first_pass_spots_->m); - auto &entries = first_pass_spots_->entries; - if (std::none_of(entries.begin(), entries.end(), - [&](const auto &e) { return e.first == centre_spec->key; })) { - constexpr size_t MAX_SPOT_KEYS = 8; - if (entries.size() == MAX_SPOT_KEYS) - entries.erase(entries.begin()); - entries.emplace_back(centre_spec->key, centre_spec->spots); - } - } - alt = score_schemes(centre_spec->ris, *indexer); - } else - alt = pick_best(*indexer_pool, *indexer); - centre_spec.reset(); - const bool header_indexes = best.result.has_value() && best.score > majority; - bool measured_indexes = alt.result.has_value() && alt.score > majority; - if (!header_indexes) { - // Nothing defends the file's centre: the data do not index at it, and the - // background says it is `moved` px away. So try the measured one, and if that - // does not index either, ask the AXIS-SIGN question again there. The two unknowns - // are discrete and coupled, and the sign rescue above asked its question at a - // centre that may itself be wrong - each error hides the other, and only moving - // both reveals anything. - bool joint = false; - const auto gonio = experiment_.GetGoniometer(); - if (!measured_indexes && gonio.has_value() && gonio->IsScanning()) { - GoniometerAxis flipped = *gonio; - flipped.Axis(-gonio->GetAxis()); - experiment_.Goniometer(flipped); - const FirstPass both = pick_best(*indexer_pool, *indexer); - if (both.result.has_value() && both.score > majority) { - logger.Warning("Beam centre from the background: the file's centre " - "({:.2f},{:.2f}) indexes {}/{} validation frames and the " - "measured ({:.2f},{:.2f}) {}/{}, but the measured centre " - "WITH the opposite rotation-axis sign indexes {}/{}. The " - "file has both wrong, and neither shows while the other " - "stands; the run takes the measured centre and the opposite " - "sign", header_x, header_y, best.score, - static_cast(validation.size()), measured_x, measured_y, - alt.score, static_cast(validation.size()), both.score, - static_cast(validation.size())); - alt = both; - measured_indexes = true; - joint = true; - } else { - experiment_.Goniometer(*gonio); - } - } - // Third coupled unknown: how DEEP the pass reads. A centre 170 px out and a spot - // list too deep to index at ANY centre hide each other exactly as the centre and - // the axis sign do - the measured centre is right, and it is thrown away because - // the full list indexes nothing there either. Measured on one such crystal: the - // file's centre indexed 1/60 and the measured one 0/60, so the file's was kept and - // the run died; with the ladder's lean seed the measured centre indexes 60/60, and - // the run then indexes every frame of the sweep and merges to the deposited - // resolution. So ask the ladder here too, before giving the centre up. - // - // The depth is asked at BOTH centres, and the centre moves only if it is the - // centre that paid. A lean rung at the measured centre scored against the whole - // spot list at the file's changes two things at once, and the count cannot say - // which of them the gain came from - the same reason a trial centre is scored on - // its own spots at all. Measured on one crystal: the measured centre's ladder - // reached 40/60 on a FIVE-fold axis harmonic while the file's own centre reached - // 49/60 on the deposited cell, and the run moved the centre, shipped the - // supercell, and skipped the beam-centre walk that would have found the crystal's - // own. An exact tie keeps the file's centre: nothing then says the centre is wrong. - bool lean = false; - bool deeper_at_home = false; - if (!measured_indexes) { - const FirstPassState before = committed; - restore_beam_center(header_x, header_y); - const auto home = first_pass_ladder(best.score, 1, majority + 1); - const FirstPassState home_state = committed; - committed = before; - apply_state(committed); - try_beam_center(measured_x, measured_y); - const auto away = first_pass_ladder(alt.score, 1, majority + 1); - const FirstPassState away_state = committed; - // Both centres index a majority, and on the same lattice. Then the two counts - // are two readings of one hypothesis at two geometries, and which is higher is - // the count's noise: measured on one crystal, the same two ladders gave the - // measured centre 56/60 against the file's 50/60, and 48/60 against 50/60 once - // a few thousand defective pixels were masked - and the run that then kept the - // file's centre, 12.6 px out, merged with twice the frames rejected. So the - // count decides nothing here: the file's centre stays for this pass, and the - // measured one is carried forward to be judged on the merges (RunAllPasses). - const bool same_lattice = - home && away - && home->result->search_result.system == away->result->search_result.system - && home->result->search_result.centering - == away->result->search_result.centering - && std::max(home->vol, away->vol) < 1.02 * std::min(home->vol, away->vol); - if (home && (!away || same_lattice || home->score >= away->score)) { - committed = home_state; - restore_beam_center(header_x, header_y); - apply_state(committed); - best = *home; - deeper_at_home = true; - if (same_lattice && geometry_prepass) { - beam_center_alternative_ = std::array{measured_x, measured_y}; - beam_center_alternative_state_ = away_state; - logger.Warning("Beam centre check: reading less of each frame, the " - "file's centre indexes {}/{} validation frames and the " - "measured centre {:.2f} px away {}/{}, on the same " - "lattice - the counts cannot say which centre is right, " - "so the run carries BOTH centres forward and adopts " - "whichever measures more of the crystal", - home->score, static_cast(validation.size()), moved, - away->score, static_cast(validation.size())); - } else - logger.Warning("Beam centre check: reading less of each frame takes the " - "file's centre to {}/{} validation frames against {}/{} at " - "the measured centre {:.2f} px away - it is how deep this " - "run was reading and not the centre, so the centre stays", - home->score, static_cast(validation.size()), - away ? away->score : 0, - static_cast(validation.size()), moved); - } else if (away) { - alt = *away; - measured_indexes = true; - lean = true; - } - } - if (measured_indexes) { - if (lean) - logger.Warning("Beam centre from the background: neither centre indexes the " - "whole spot list, but at the measured ({:.2f},{:.2f}) the " - "lean first pass above indexes {}/{} against the file's " - "{}/{} - the beam centre in this file is wrong, and only a " - "pass that reads less of each frame could see it", - measured_x, measured_y, alt.score, - static_cast(validation.size()), best.score, - static_cast(validation.size())); - else if (!joint) - logger.Warning("Beam centre from the background: the file's ({:.2f},{:.2f}) " - "indexes {}/{} validation frames and the measured " - "({:.2f},{:.2f}) indexes {}/{} - the beam centre in this " - "file is wrong, and the run adopts the measured one", - header_x, header_y, best.score, - static_cast(validation.size()), measured_x, measured_y, - alt.score, static_cast(validation.size())); - best = alt; - } else if (!deeper_at_home) { - // Nothing worked. Put the file's centre back, so a run that fails for some - // other reason fails at the geometry it was given rather than at a fit that has - // just been tested and found no better. - // - // Unless the frame counts are both at the floor and the POOLED spots do decide: - // a crystal whose spots are a minority of every frame indexes 0/60 at any centre, - // and then only the pooled excess over chance - each lattice measured against its - // own wrong-spindle null, so a larger cell is not credited for fitting more by - // chance - says which centre sees it. - // - // Each centre is judged on the best lattice it gives, so the long-axis rescue - // is asked at BOTH before they are compared. The lattices standing here index - // no majority anywhere, which is exactly the pass that rescue exists for, and - // the two counts of spots on two lattices that are not the crystal's are close - // to a coin toss: measured on a 545 A axis, the measured centre (3.6 px out - // across the spindle) won it on 16.7% against 15.5%, the rescue then ran - // there alone and could not see the long axis, and the run shipped a P1 - // sub-cell - while at the file's centre the same rescue indexes 60/60. - if (auto rescued = long_axis_rescue(alt)) - alt = *rescued; - PooledEvidence at_measured, at_header; - if (alt.result.has_value()) - at_measured = pooled_evidence(*indexer, *alt.result); - restore_beam_center(header_x, header_y); - if (auto rescued = long_axis_rescue(best)) - best = *rescued; - if (best.result.has_value()) - at_header = pooled_evidence(*indexer, *best.result); - long_axis_asked = true; - if (alt.result.has_value() && beats_chance(at_measured) - && ValidationEvidencePrefers(at_header, at_measured)) { - try_beam_center(measured_x, measured_y); - logger.Warning("Beam centre from the background: neither centre indexes a " - "majority of the validation frames ({}/{} at the measured, " - "{}/{} at the file's), but the pooled spots put " - "{:.1f}% on the lattice at the measured ({:.2f},{:.2f}) against " - "{:.1f}% at the file's (chance {:.1f}% and {:.1f}%) - the run " - "adopts the measured centre", - alt.score, static_cast(validation.size()), - best.score, static_cast(validation.size()), - 100.0 * at_measured.on_lattice / std::max(1, at_measured.spots), - measured_x, measured_y, - 100.0 * at_header.on_lattice / std::max(1, at_header.spots), - 100.0 * at_measured.by_chance / std::max(1, at_measured.spots), - 100.0 * at_header.by_chance / std::max(1, at_header.spots)); - best = alt; - } else { - logger.Info("Beam centre check: the file's centre ({}/{} frames) and the " - "measured one {:.2f} px away ({}/{}) - the measured centre puts no " - "more of the pooled spots on its lattice, so the file's centre is " - "kept", - best.score, static_cast(validation.size()), moved, alt.score, - static_cast(validation.size())); - } - } - } else { - restore_beam_center(header_x, header_y); - if (measured_indexes) { - // Primitive volumes, as everywhere else here: two passes can find the same - // lattice in different settings, and a centred setting's cell is an exact - // integer multiple of its primitive one. - const double ratio = alt.vol < best.vol ? best.vol / alt.vol : alt.vol / best.vol; - const auto &hc = best.result->search_result.conventional.GetUnitCell(); - const auto &mc = alt.result->search_result.conventional.GetUnitCell(); - const bool same_class = - best.result->search_result.system == alt.result->search_result.system - && best.result->search_result.centering == alt.result->search_result.centering; - // The primitive volume is what the beam centre decides. A centre error along - // the spindle makes the FFT take an axis harmonic, which changes the volume by - // an integer factor or by sqrt(3); the Bravais class can differ for a reason - // that has nothing to do with the centre, so the two are reported separately. - if (ratio < 1.02 && same_class) { - Note(fmt::format("beam centre kept: measured {:.1f} px off, same lattice", moved)); - logger.Info("Beam centre check: indexing at the measured centre " - "({:.2f},{:.2f}) returns the same lattice as the file's - " - "{}-centred {}, cells agreeing to {:.2f} % in volume. The {:.2f} px " - "difference does not decide this crystal's cell", - measured_x, measured_y, alt.result->search_result.centering, - gemmi::crystal_system_str(alt.result->search_result.system), - 100.0 * (ratio - 1.0), moved); - // Unless the two are further apart than this geometry absorbs. A centre - // wrong by about one reflection spacing indexes the same cell just as - // well, with every reflection one index out over part of the sweep, and - // only the merge shows it - so the measured centre is kept as a hypothesis - // for RunAllPasses, which tries it if the file's centre merges - // inconsistently. - if (geometry_prepass && moved > 3.0f * measured.sigma_pxl - && moved > BeamCenterNeed_pxl(experiment_)) - beam_center_same_lattice_alternative_ = - std::array{measured_x, measured_y}; - } else if (ratio < 1.02) { - // Same cell, different metric symmetry read off it. The cell is not in - // question, but the CLASS is not cosmetic either: GeometryRefiner holds - // the cell to it for the rest of the run, so a class that is one axis - // equality away from the truth is refined, predicted, integrated and - // merged with that equality imposed. - // - // Whether the centre decided it is arithmetic. The equality is held to - // LATTICE_SEARCH_DIST_TOLERANCE, relative, on a pair of reduced-cell - // axes; a centre wrong by `moved` px moves an axis of length a by - // moved*p*a/(D*lambda) relative, the FFT reading that axis off spot - // positions the centre displaces. Both axes of the pair carry their own - // error, so the centre can open or close a gap of their sum - and where - // that reaches the tolerance the higher class is not a fact about the - // crystal but a consequence of where the beam was assumed to be. - // - // Nothing available here can settle it: both centres index, and the frame - // count is the one statistic this check already refuses to arbitrate on. - // So carry the other centre out of the pass and let the two arms be judged - // on their merges (RunAllPasses). - const auto reduced = best.result->search_result.primitive_reduced.GetUnitCell(); - const double axis_pair = std::max(reduced.a + reduced.b, - reduced.b + reduced.c); - const double induced = moved * experiment_.GetPixelSize_mm() * axis_pair - / (experiment_.GetDetectorDistance_mm() - * experiment_.GetWavelength_A()); - const bool centre_can_decide = induced >= LATTICE_SEARCH_DIST_TOLERANCE; - logger.Info("Beam centre check: indexing at the measured centre " - "({:.2f},{:.2f}) returns the same CELL as the file's - the " - "primitive volumes agree to {:.2f} % - but a different metric " - "symmetry: {}-centred {} from the file's centre against " - "{}-centred {} from the measured one", measured_x, measured_y, - 100.0 * (ratio - 1.0), best.result->search_result.centering, - gemmi::crystal_system_str(best.result->search_result.system), - alt.result->search_result.centering, - gemmi::crystal_system_str(alt.result->search_result.system)); - if (centre_can_decide && geometry_prepass) { - beam_center_alternative_ = std::array{measured_x, measured_y}; - logger.Warning("Beam centre check: {:.2f} px is worth {:.2f} % on the " - "axis pair the promotion compares, against the {:.2f} % " - "the metric symmetry is decided to - the difference is " - "large enough to have decided it either way, so the run " - "carries BOTH centres forward and adopts whichever " - "merges better", moved, 100.0 * induced, - 100.0 * LATTICE_SEARCH_DIST_TOLERANCE); - } else { - logger.Info("Beam centre check: {:.2f} px is worth {:.2f} % on the axis " - "pair the promotion compares, against the {:.2f} % the " - "metric symmetry is decided to - the run keeps the file's " - "centre and the symmetry it found", moved, 100.0 * induced, - 100.0 * LATTICE_SEARCH_DIST_TOLERANCE); - } - } else { - const double nearest = std::round(ratio); - const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15) - || std::abs(ratio - std::sqrt(3.0)) < 0.15; - logger.Warning("Beam centre check: the two centres DISAGREE about the " - "lattice. The file's ({:.2f},{:.2f}) gives {}-centred {}, " - "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} (primitive volume " - "{:.0f} A^3, {}/{} frames); the measured ({:.2f},{:.2f}) gives " - "{}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} " - "(primitive volume {:.0f} A^3, {}/{} frames){}", - header_x, header_y, best.result->search_result.centering, - gemmi::crystal_system_str(best.result->search_result.system), - hc.a, hc.b, hc.c, hc.alpha, hc.beta, hc.gamma, best.vol, - best.score, static_cast(validation.size()), - measured_x, measured_y, alt.result->search_result.centering, - gemmi::crystal_system_str(alt.result->search_result.system), - mc.a, mc.b, mc.c, mc.alpha, mc.beta, mc.gamma, alt.vol, - alt.score, static_cast(validation.size()), - harmonic ? fmt::format(" - the volumes differ by {:.2f}x, so one " - "is an axis harmonic of the other and the " - "beam centre is what settles which", - ratio) - : ""); - // Deliberately no decision on the frame counts. They are the one statistic - // that must not choose: on the measured case the harmonic indexed MORE. - // - // Decided in either direction on the pooled excess over chance, the same - // measurement the failed-pass branch above adopts a centre on. Each lattice is - // scored against its own wrong-spindle null, so the larger cell is not credited - // for the spots a denser lattice catches by accident, and the smaller one is - // not credited for the spots it misses. - // - // Both directions need it. Where the measured centre holds the LARGER cell, - // nothing else in the run repairs it - the second pass refuses a lattice larger - // than pass 1's (the supercell guard below). Where it holds the SMALLER one, - // the second pass can drop to it only if its own beam refinement happens to - // walk the whole centre error, and that rescue is chaotic: measured on one - // crystal, a 0.1 % change in the pass-1 distance switched off a 6.5 px walk and - // the run kept the doubled axis. - const double nearest_n = std::round(ratio); - const bool integer_harmonic = nearest_n >= 2.0 && nearest_n <= 4.0 - && std::abs(ratio - nearest_n) < 0.15; - bool adopted_measured = false; - if (integer_harmonic) { - const PooledEvidence at_header = pooled_evidence(*indexer, *best.result); - try_beam_center(measured_x, measured_y); - const PooledEvidence at_measured = pooled_evidence(*indexer, *alt.result); - adopted_measured = beats_chance(at_measured) - && ValidationEvidencePrefers(at_header, at_measured); - logger.Info("Beam centre check: pooled validation spots on the lattice " - "{:.1f}% at the file's centre and {:.1f}% at the measured " - "(chance {:.1f}% and {:.1f}%)", - 100.0 * at_header.on_lattice / std::max(1, at_header.spots), - 100.0 * at_measured.on_lattice - / std::max(1, at_measured.spots), - 100.0 * at_header.by_chance / std::max(1, at_header.spots), - 100.0 * at_measured.by_chance - / std::max(1, at_measured.spots)); - if (adopted_measured) { - logger.Warning("Beam centre check: the spots are on the {} cell of the " - "measured centre, not on the file's - so the {:.2f}x is " - "the file's centre error and not the crystal, and the " - "run adopts the measured centre", - alt.vol > best.vol ? "larger" : "smaller", ratio); - best = alt; - } else { - restore_beam_center(header_x, header_y); - } - } - if (!adopted_measured) - logger.Warning("Beam centre check: the run keeps the file's centre. The " - "frame counts cannot arbitrate this - a cell twice as long " - "must place every spot twice as accurately to score the " - "same, so the harmonic is often the one that indexes more. " - "Where the beam centre is known, give it with " - "--beam-x/--beam-y"); - } - } else { - logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes " - "{}/{} validation frames against the file's {}/{} - the file's " - "centre stands", measured_x, measured_y, alt.score, - static_cast(validation.size()), best.score, - static_cast(validation.size())); - } - } - } - } - - // The long-axis rescue (above), at the centre the run now stands on. - if (!long_axis_asked) - if (auto rescued = long_axis_rescue(best)) - best = *rescued; - // The end wedge (above), at the geometry the rescues settled on. - best = add_end_wedge(std::move(best), *indexer_pool, *indexer); - release_spot_engines(); // no more frames to read; give the cards their memory back - - if (!cancelled_) { - // The two ways the first pass can end without a usable lattice that are about the DATA. - // They are one situation from the user's side - no crystal was found in these frames - so - // both name the cause and then give the same remedies. - constexpr const char *next_steps = - " Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise " - "--max-spots if the frames are well populated with spots, and use --mode azint to " - "write the per-image analysis for inspection."; - if (!best.result.has_value()) { - // An indexer that FAILED is a third outcome, and none of those remedies fit it: it - // never examined the frames, so nothing about the data or the processing options is - // implicated and changing them cannot help. Name the failure instead. - if (best.indexer_error) - throw JFJochException(JFJochExceptionCategory::IndexerError, - fmt::format("Rotation indexing could not run: {}. Nothing was " - "integrated. The indexer never examined the frames, " - "so this is a failure of the machine and not of the " - "data - free the resource it needs (a GPU shared with " - "another process is the usual cause) and run again.", - *best.indexer_error)); - throw NoLatticeFound(std::string("Two-pass rotation indexing found no lattice.") + next_steps); - } - // Short-axis hypothesis - the small-molecule case. - // - // LAST, after the axis-sign, beam-centre and long-axis rescues, deliberately. Those all - // run the standard floor and compare their result against `best` on volume ratios; a cell - // adopted here is a k-fold sub-cell, so running before them makes every one of those - // comparisons read an integer ratio and report an axis harmonic that is not one (measured: - // the beam-centre check announced a "3.00x harmonic ... the beam centre is what settles - // which" on a pair that differs only by the FLOOR). Here the comparison is like with like, - // and this pass sees whatever geometry those rescues settled on. - // - // fft_min_unit_cell_A (10 A) is a protein constant: shorter than any protein axis and longer - // than a small molecule's. It is applied to the FFT's PER-DIRECTION PEAK SEARCH - // (FFTIndexerCPU.cpp / FFTIndexerGPU.cu skip every histogram bin below it) and only then to - // the assembled cell (FFTIndexer::ReduceAndRefine, PostIndexingRefinement). The first of - // those is what does the damage: the transform reports ONE peak per direction, so a - // lattice row whose real-space repeat is shorter than the floor is reported at its first - // HARMONIC above it, and the cell built out of those harmonics is an exact integer - // sublattice of the true one - a real-space supercell. Measured on chemical-crystallography - // data: index 2, 3 and 6 exactly, and the one dataset whose true reduced cell clears 10 A is - // the one that comes out right. Because the floor is applied to the REDUCED cell it rejects - // the lattice and not merely a basis choice, so nothing downstream can re-describe its way - // back. - // - // Lowering the floor by default is not available: it only ever ADDS candidates, and the - // Bravais classification the candidate pool feeds is under-determined, so it moves crystals - // whose cell it never touched (measured: 2 regressions and 1 improvement over 20 protein - // datasets, all three in the metric symmetry rather than the cell). So the low floor runs as - // a SECOND hypothesis and the standard pass's answer stands unless the evidence takes it. - // - // What the evidence can be, before integration: a sublattice's per-spot residuals are - // dominated pointwise by the finer lattice it contains, so every margin statistic - squares, - // counts, trimmed sums - is biased toward the LARGER cell, and a spurious supercell indexes - // every frame its true cell does. Only the SATURATION is informative. Hence the same shape - // the scheme ratchet in pick_best already uses: an integer primitive-volume ratio plus a tie - // on validation frames, or an outright win, which runs against the bias and so is real. The - // statistic that genuinely decides an axis multiple is the occupancy of the reflection class - // the multiple adds - AbsenceEvidence in SearchSpaceGroup - and that needs integrated - // intensities, so it cannot run here. - // Not when a cell was given. RaiseFFTBoundForKnownCell has already lowered the floor for - // the standard pass against that cell, so this pass has nothing left to add - and it takes - // its settings from experiment_, which never carried the RAISED max bound, so on a cell - // longer than the default search it could only ever return a harmonic of what was asked - // for, tie on frames the way a sub-multiple always does, and be adopted over it. - if (!cancelled_ && !experiment_.GetUnitCell().has_value()) { - std::unique_ptr short_pool; - std::vector> sris; - if (short_spec && short_spec->key == ExperimentKey() && short_spec->mapping == spot_mapping - && short_spec->spots_per_image == first_pass_spots_per_image) { - for (auto &f : short_spec->futures) - f.get(); - short_pool = std::move(short_spec->pool); - sris = std::move(short_spec->sris); - } else { - short_spec.reset(); - short_pool = std::make_unique(short_axis_settings(experiment_), - IndexerConstruction::OnFirstUse); - sris.reserve(schemes.size()); - for (const auto &[name, ordinals] : schemes) { - if (cancelled_) - break; - sris.push_back(feed_scheme(*short_pool, ordinals)); - } - std::vector> short_futures; - short_futures.reserve(sris.size()); - for (auto &ri : sris) { - RotationIndexer *rp = ri.get(); - short_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); })); - } - for (auto &f : short_futures) - f.get(); - } - - for (size_t i = 0; i < sris.size() && !cancelled_; i++) { - auto found = sris[i]->GetLattice(); - if (!found.has_value()) - continue; - RotationIndexerResult alt = std::move(*found); - const std::string &name = schemes[i].first; - // Primitive volumes, for the reason pick_best gives: two searches can find the same - // lattice in different settings, and a centred setting's cell is an exact integer - // multiple of its primitive one. - const double vol = std::abs( - alt.lattice.ToPrimitive(alt.search_result.centering).CalcVolume()); - if (!(vol > 1.0) || !(best.vol > 1.0)) - continue; - const double ratio = best.vol / vol; - const double nearest = std::round(ratio); - // The gate is here rather than after the scoring because scoring is the expensive - // half - count_indexed runs the real per-image path over the validation set, seconds - // on a protein - while the volume ratio is free. On every protein measured the low - // floor returns the SAME lattice (ratio 1.00-1.01), so this costs them one extra FFT - // and nothing else. - if (!(vol < best.vol) || nearest < 2.0 || std::abs(ratio - nearest) >= 0.15) { - logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': primitive volume " - "{:.0f} A^3 against the standing cell's {:.0f} A^3, ratio {:.2f} - " - "not an integer sub-cell, keeping the standing cell", - IndexingSettings::fft_min_unit_cell_limit_A, name, vol, best.vol, - ratio); - continue; - } - - // The premise of this pass, checked rather than asserted. The pass exists because a - // lattice row shorter than the standard floor is reported at its first HARMONIC - // above it, so a cell it may rightfully adopt has a reduced axis BELOW that floor - - // and the floor is applied to the reduced cell, which is the one to measure. Where - // the alternative has no such axis the floor cannot be what the standard pass - // tripped over, and an integer volume ratio is then an ordinary axis multiple: the - // question the axis-harmonic arbiter settles just above on the occupancy of the - // class the multiple adds. That is intensity-class evidence about the very - // reflections in dispute, where the tie below is a frame count this comment already - // says a sublattice cannot lose - so where the two disagree the arbiter is right and - // this block must not overrule it. Measured on the corpus: the four small-molecule - // sets the pass is for all have a reduced axis under 9 A, and the two protein sets - // it halved have 35 and 47 A. - const auto alt_reduced = alt.lattice.ToPrimitive(alt.search_result.centering) - .NiggliReduce().GetUnitCell(); - const float alt_shortest = std::min({alt_reduced.a, alt_reduced.b, alt_reduced.c}); - const float standard_floor = experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A(); - if (alt_shortest >= standard_floor) { - logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': a {:.0f}x sub-cell " - "(vol {:.0f} A^3), but its shortest reduced axis is {:.1f} A, above " - "the {:.1f} A floor the standard pass already searched - the floor " - "is not why that cell was missed, so this is an ordinary axis " - "multiple and not this pass's to decide; keeping the standing cell", - IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol, - alt_shortest, standard_floor); - continue; - } - - const int alt_score = count_indexed(*indexer, alt); - // A pass that indexes nothing is not tied with a pass that indexes nothing. Without - // this, 0 against 0 passes the tie band and the volume ratio alone decides, which is - // the one thing this block must never do. - if (alt_score <= 0 || best.score <= 0) { - logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {}/{} validation " - "frames against {} - nothing indexes, keeping the standing cell", - IndexingSettings::fft_min_unit_cell_limit_A, name, alt_score, - static_cast(validation.size()), best.score); - continue; - } - const bool tied = static_cast(alt_score) >= best.score * 0.9f - 0.5f; - const bool clearly_more = static_cast(alt_score) > best.score * 1.1f + 0.5f; - if (!tied && !clearly_more) { - logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {:.0f}x sub-cell " - "(vol {:.0f} A^3) indexes {}/{} validation frames against {} - " - "keeping the standing cell", - IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol, - alt_score, static_cast(validation.size()), best.score); - continue; - } - const auto uc = alt.lattice.GetUnitCell(); - logger.Warning("Short-axis pass (floor {:.0f} A), scheme '{}': the standing " - "cell is a {:.0f}x supercell. Its {:.0f} A^3 primitive volume against " - "this pass's {:.0f} A^3, at {}/{} validation frames against {} - " - "adopting {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}. The shortest axis " - "is below the {:.0f} A default FFT floor, which is why the standard " - "pass could only build a multiple of it", - IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, best.vol, - vol, alt_score, static_cast(validation.size()), best.score, - uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, - experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A()); - best.score = alt_score; - best.vol = vol; - best.name = name + " (short-axis)"; - best.result = std::move(alt); - } - } - - // Sub-lattice hypothesis: the standing cell's seven index-2 sub-lattices, in its own - // orientation and geometry, each scored on the validation frames. - // - // The FFT can only offer the cells its peaks build, and two lattices seen together - a twin - // domain or a second crystal whose reciprocal lattice shares an index-2 sub-lattice with - // this one's - build the lattice that contains them both: a cell twice the size of either, - // which indexes both crystals' spots and so wins every comparison of candidates the FFT - // returns. Where the true cell is then never offered, nothing downstream can choose it. - // Offering it costs seven frame counts. - // - // The frame count only nominates: a coarser lattice asks less accuracy of every spot, so a - // sub-lattice of the true cell can index as many frames as the cell itself, or a few more - // where the cell's own score is depressed. What decides is the comparison this pass already - // makes between a cell and its axis harmonic: the share of the validation frames' pooled - // spots each puts on itself over and above what a wrong spindle angle puts there, each - // against its own null (ValidationEvidencePrefers). A real superstructure keeps its cell on - // that count, because the class it adds holds reflections; a cell built out of two crystals - // loses it, because it has to place every spot of the one it describes twice as accurately - // with an orientation fitted to both. The nominee must also index a majority of the frames, - // the bar the false-promotion test above asks of a cell that "works", and no fewer frames - // than the standing cell: a standing cell that indexes almost nothing is not this crystal's - // lattice, its sub-lattices are not hypotheses about it, and the rescues below are what - // such a run needs. - if (!cancelled_ && best.result.has_value() && best.score > 0) { - const CrystalLattice prim = best.result->lattice - .ToPrimitive(best.result->search_result.centering); - const Coord A = prim.Vec0(), B = prim.Vec1(), C = prim.Vec2(); - std::optional sub_best; - std::string counts; - for (int u = 0; u < 2; u++) - for (int v = 0; v < 2; v++) - for (int w = 0; w < 2; w++) { - if (u + v + w == 0 || cancelled_) - continue; - const Coord t = (A * static_cast(u) + B * static_cast(v) - + C * static_cast(w)) / 2.0f; - const CrystalLattice sub = (w ? CrystalLattice(A, B, t) - : v ? CrystalLattice(A, t, C) - : CrystalLattice(t, B, C)).NiggliReduce(); - RotationIndexerResult cand = *best.result; - cand.lattice = sub; - cand.search_result = LatticeSearchResult{}; - cand.search_result.primitive_reduced = sub; - cand.search_result.conventional = sub; - cand.unconstrained.reset(); - cand.extra_lattices.clear(); - const int score = count_indexed(*indexer, cand); - counts += fmt::format(" ({},{},{}) {}", u, v, w, score); - if (!sub_best || score > sub_best->score) - sub_best = FirstPass{std::move(cand), score, - fmt::format("{} (sub-lattice {}{}{})", best.name, u, v, w), - std::abs(sub.CalcVolume())}; - } - bool adopt = false; - std::string verdict = " - none indexes a majority of the frames and as many as the standing cell"; - if (sub_best && sub_best->score >= best.score - && sub_best->score > static_cast(validation.size()) / 2) { - const PooledEvidence at_standing = pooled_evidence(*indexer, *best.result); - const PooledEvidence at_sub = pooled_evidence(*indexer, *sub_best->result); - adopt = ValidationEvidencePrefers(at_standing, at_sub); - verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance) against " - "the standing cell's {}/{} ({}), {}", - at_sub.on_lattice, at_sub.spots, at_sub.by_chance, - at_standing.on_lattice, at_standing.spots, at_standing.by_chance, - adopt ? "adopting it" : "keeping the standing cell"); - } - logger.Info("Sub-lattice hypothesis: the standing cell ({:.0f} A^3) indexes {}/{} validation " - "frames; its index-2 sub-lattices (half-translation u,v,w):{}{}", - best.vol, best.score, static_cast(validation.size()), counts, verdict); - if (adopt) - best = std::move(*sub_best); - } - - // A lattice that indexes almost none of the validation frames is not this crystal's lattice. - // Forcing it anyway integrates every image against it; the few partials that survive yield no - // fulls in the 3D combine, and the run then dies inside the merge complaining about resolution - - // far from the cause. Stop here instead, while the reason is still visible. The floor is one - // frame in six, several times below the weakest crystal that still merges. It is checked AFTER - // the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue - // score. - // The second (refined-geometry) pass re-indexes DE NOVO, and the post-refined geometry - fitted - // to the lattice pass 1 already found - can tip the blind FFT, or the long-axis rescue above, - // onto another lattice: an axis harmonic of a long cell (measured: a tripled ~100 A c-axis at - // 7/60), or a related centred cell that indexes fewer frames than pass 1's lattice does at - // this very geometry (measured: a C-centred 303 A cell at 41/60 against pass 1's triclinic - // cell at 60/60). Pass 1's lattice is therefore carried here as a hypothesis of its own and - // the two are judged on the same validation frames: where the re-index found a different - // lattice and indexes fewer frames than pass 1's, or too few to be this crystal's lattice at - // all, pass 1's whole result is integrated instead. The same lattice found again is kept as - // the re-index refined it at this geometry. This is the floor-side companion of the - // supercell-collapse guard in RunAllPasses, which only sees a pass-2 lattice that is LARGER. - if (!geometry_prepass && prepass_result_.has_value()) { - RotationIndexerResult pass1_lattice = *prepass_result_; - // Pass 1's result carries pass 1's goniometer; forcing it whole would put the - // uncorrected angles back (the same repair as the supercell re-run in RunAllPasses). - if (prepass_rotation_scale_ && pass1_lattice.axis) - pass1_lattice.axis = ScaleRotation(*pass1_lattice.axis, *prepass_rotation_scale_); - // Pass 1's lattice was FITTED at pass 1's detector distance, and this pass - // integrates at the post-refined one. Re-scoring it here is not the same as - // re-fitting it: a real-space cell is measured against the distance the spots were - // seen at, so carrying it across a distance change scales the whole cell by the - // ratio of the two. Measured on a 545 A axis over 225 -> 226 mm: 0.4 %, exactly - // 225/226, shipped as the run's answer. Scale it to the distance it will be used - // at; the orientation is unaffected. - const float d_fit = prepass_result_->geom.GetDetectorDistance_mm(); - const float d_use = experiment_.GetDetectorDistance_mm(); - const bool rescaled = d_fit > 0 && d_use > 0 && d_fit != d_use; - if (rescaled) { - const float f = d_use / d_fit; - const auto &l = pass1_lattice.lattice; - pass1_lattice.lattice = CrystalLattice(l.Vec0() * f, l.Vec1() * f, l.Vec2() * f); - } - const CrystalLattice pass1_primitive = pass1_lattice.lattice - .ToPrimitive(pass1_lattice.search_result.centering).NiggliReduce(); - const double pass1_vol = std::abs(pass1_primitive.CalcVolume()); - // The same lattice whatever setting or class it is held in - a centred cell and the - // triclinic primitive of the same points are one lattice, and the re-index's symmetric - // setting must not lose to its own unconstrained primitive on a frame count: the Niggli- - // reduced primitive edges agree within 2 % (the battery's lattice-identity test). - const auto reduced_edges = [](const CrystalLattice &l) { - const auto uc = l.GetUnitCell(); - std::array e = {uc.a, uc.b, uc.c}; - std::sort(e.begin(), e.end()); - return e; - }; - const auto e1 = reduced_edges(pass1_primitive); - const auto e2 = reduced_edges(best.result->lattice - .ToPrimitive(best.result->search_result.centering).NiggliReduce()); - bool same_lattice = true; - for (int i = 0; i < 3; i++) - same_lattice = same_lattice && std::abs(e2[i] - e1[i]) <= 0.02 * e1[i]; - const int floor = static_cast(validation.size()) / 6; - const int pass1_score = (best.score < floor || !same_lattice) - ? count_indexed(*indexer, pass1_lattice) : -1; - if (best.score < floor || pass1_score > best.score) { - const auto &pc = pass1_lattice.search_result.conventional.GetUnitCell(); - const auto &bc = best.result->search_result.conventional.GetUnitCell(); - logger.Warning("Two-pass: the refined-geometry re-index indexes {}/{} validation frames " - "({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}), pass-1's lattice " - "{}/{} at this geometry ({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} " - "{:.2f}) - integrating with pass-1's lattice", - best.score, static_cast(validation.size()), - best.result->search_result.centering, - gemmi::crystal_system_str(best.result->search_result.system), - bc.a, bc.b, bc.c, bc.alpha, bc.beta, bc.gamma, - pass1_score, static_cast(validation.size()), - pass1_lattice.search_result.centering, - gemmi::crystal_system_str(pass1_lattice.search_result.system), - pc.a, pc.b, pc.c, pc.alpha, pc.beta, pc.gamma); - if (rescaled) - logger.Info("Two-pass: pass-1's cell was fitted at {:.3f} mm and is used at " - "{:.3f} mm - scaled by {:.6f} to the distance it is integrated at", - d_fit, d_use, d_use / d_fit); - best.result = std::move(pass1_lattice); - best.score = pass1_score; - best.vol = pass1_vol; - best.name = "pass-1 lattice (refined-geometry re-index found less)"; - } - } - - // The last resort, and ONLY on a run that has no usable lattice: the pass-1 fallback above - // has had its turn, and below this the run stops with an error. A run that indexes must - // not pay for it - measured, an unguarded call spent 41-76 % of the wall clock on a - // healthy dataset evaluating rungs that arithmetically could not be adopted (one crystal - // 41 s guarded against 175 s unguarded, same answer). - if (best.score < static_cast(validation.size()) / 6) { - if (auto rescued = first_pass_ladder(best.score, - static_cast(validation.size()) / 6, 0)) - best = *rescued; - } - - // Whether the run has a lattice at all is decided here, on the pooled spots rather than on - // the frame count the rescues above are steered by. A frame count asks each frame to clear - // a bar of its own and then counts the frames that did, which measures the spots per frame - // as much as the lattice; the pooled fraction asks the whole validation set at once and is - // judged against what the same lattice scores at a wrong spindle angle, so a sparse pattern - // is refused only when its spots genuinely do not know about the lattice. - const PooledEvidence evidence = pooled_evidence(*indexer, *best.result); - const double pooled = static_cast(evidence.on_lattice) - / static_cast(std::max(1, evidence.spots)); - const double pooled_chance = static_cast(evidence.by_chance) - / static_cast(std::max(1, evidence.spots)); - result.validation_evidence = evidence; - if (memo_hit && (memo_hit->spots != evidence.spots || memo_hit->on_lattice != evidence.on_lattice - || memo_hit->by_chance != evidence.by_chance)) - throw MemoMismatch(fmt::format("First-pass memo mismatch: stored {}/{}/{}, recomputed {}/{}/{}", - memo_hit->on_lattice, memo_hit->spots, memo_hit->by_chance, - evidence.on_lattice, evidence.spots, evidence.by_chance)); - // Kept only where the pass ended on the inputs it started from. - if (!geometry_prepass && FirstPassInputKey() == first_pass_key) - AddFirstPassMemo(FirstPassMemo{first_pass_key, evidence}); - // A pass run only to score a lattice has its score, whatever it is: a lattice that does not - // beat chance is a result for the comparison that asked, not a reason to stop the run. - if (indexing_probe_only_) { - logger.Info("Indexing probe: scheme '{}', {}/{} validation frames, {}/{} validation spots " - "= {:.1f}% against {:.1f}% at a wrong spindle angle", best.name, best.score, - static_cast(validation.size()), evidence.on_lattice, evidence.spots, - 100.0 * pooled, 100.0 * pooled_chance); - return true; - } - if (!beats_chance(evidence)) { - { - // Name the cell and Bravais class that was rejected. The commonest cause is a metric - // symmetry promoted one class too far - the constrained cell then misses every - // reflection by the small angle the constraint snapped away - and without the cell in - // the message there is nothing to see that from. - const auto &c = best.result->search_result.conventional.GetUnitCell(); - throw NoLatticeFound(fmt::format("Two-pass rotation indexing found only a lattice that " - "takes {:.1f}% of the validation spots against the " - "{:.1f}% a wrong spindle angle takes ({} distinct " - "reflections against {} over {} wrong angles) - it is " - "not this crystal's lattice. It was {}-centred {}, " - "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}", - 100.0 * pooled, 100.0 * pooled_chance, - evidence.reflections, evidence.reflections_by_chance, - VALIDATION_NULL_DISPLACEMENTS, - best.result->search_result.centering, - gemmi::crystal_system_str(best.result->search_result.system), - c.a, c.b, c.c, c.alpha, c.beta, c.gamma, next_steps)); - } - } - indexer->ForceRotationIndexerResult(*best.result); - // A lattice the pooled spots accept while most validation frames fail the per-frame floor - // (LATTICE_MIN_INDEXED_FRACTION) on their own, because the sweep as a whole puts fewer of its - // spots on the lattice than that floor asks of each frame: the crystal is real but its - // spots are a minority of every frame - a weak crystal inside a powder of its own - // microcrystals and ice. The floor cannot tell its frames from anything else, so integrate - // every frame from the sweep's lattice, as XDS and DIALS do, and leave the frames to the - // scaling. The frames that do clear it are the upper tail of the same population: measured - // on a crystal in a powder of itself (13% of the pooled spots on the lattice), 12/60 - // cleared it, and keeping the floor integrated 15% of the sweep. A sweep that is sparse - // only in spots per frame (a small molecule) keeps the floor. - const bool sparse_lattice = best.score < static_cast(validation.size()) / 2 - && pooled < LATTICE_MIN_INDEXED_FRACTION; - indexer->IntegrateEveryFrame(sparse_lattice); - if (sparse_lattice) - logger.Warning("Sparse lattice: {}/{} validation frames clear the per-frame floor, but the " - "pooled spots put {:.1f}% on the lattice against {:.1f}% by chance - every " - "frame is integrated from the sweep's lattice", - best.score, static_cast(validation.size()), 100.0 * pooled, - 100.0 * pooled_chance); - logger.Info("First-pass spot finding: {} frames in {} batches, {:.2f} s", - spot_cache.size(), prefetch_batches, prefetch_time_s); - logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames, " - "{}/{} validation spots = {:.1f}% against {:.1f}% at a wrong spindle angle; {} " - "distinct reflections against {} over {} wrong angles)", - best.name, best.score, static_cast(validation.size()), - evidence.on_lattice, evidence.spots, 100.0 * pooled, 100.0 * pooled_chance, - evidence.reflections, evidence.reflections_by_chance, VALIDATION_NULL_DISPLACEMENTS); - { - const auto &c = best.result->search_result.conventional.GetUnitCell(); - Note(fmt::format("lattice {} {}: {:.2f} {:.2f} {:.2f} {:.1f} {:.1f} {:.1f}", - gemmi::crystal_system_str(best.result->search_result.system), - best.result->search_result.centering, c.a, c.b, c.c, c.alpha, c.beta, c.gamma)); - Note(fmt::format("{}/{} validation frames, {:.0f}% of their spots on it", - best.score, validation.size(), 100.0 * pooled)); - } - // Leftover-lattice census, report only: index what the crystal's lattice leaves over, to name - // a split crystal, a twin domain or a second crystal in the report. Every spot of the spread - // and validation frames the lattice takes is marked; the rest (not on ice) goes to a fresh - // rotation indexer; a lattice it finds is judged on the validation leftovers against the same - // displaced-spindle-angle null as the first pass, its spots are marked in turn, and the search - // goes on for up to three lattices. Its own IndexAndRefine instances and indexer, so nothing - // the run goes on with is touched. Run on the pass that writes the result only. - // - // Its spots are found here; the rest runs in the background, beside the image loop, on - // copies of everything it reads, and its answer is collected before the pass returns. Not on - // a pass that only post-refines: nothing reads its result there. - if (!geometry_prepass && !postrefine_probe_only_) { - std::vector census_frames = spread_ordinals; - census_frames.insert(census_frames.end(), validation.begin(), validation.end()); - std::ranges::sort(census_frames); - census_frames.erase(std::unique(census_frames.begin(), census_frames.end()), census_frames.end()); - prefetch_spots(census_frames); - std::map> census_spots; - for (const int f : census_frames) - census_spots.emplace(f, spot_cache.at(f)); - census_future = std::async(std::launch::async, - [&, x = experiment_, census_frames, census_spots, - lattice = *best.result, validation, spread_ordinals, - validation_settings]() -> LeftoverLatticeCensus { - const auto census_start = std::chrono::steady_clock::now(); - const bool ice_ok = x.GetIndexingSettings().GetIndexIceRings(); - const auto &spot_cache = census_spots; - const auto rot_angle = [&](int ordinal) -> std::optional { - if (const auto g = x.GetGoniometer()) - return g->GetAngle_deg(static_cast(ordinal)) + g->GetWedge_deg() / 2.0f; - return std::nullopt; - }; - // assign[frame][spot]: -2 ice, -1 no lattice, 0 the crystal's lattice, 1.. a leftover lattice - std::map> assign; - for (const int f : census_frames) { - auto &a = assign[f]; - a.assign(spot_cache.at(f).size(), -1); - for (size_t k = 0; k < a.size(); k++) - if (!ice_ok && spot_cache.at(f)[k].ice_ring) - a[k] = -2; - } - // Mark the still unassigned spots lattice `r` takes as `id`. - auto mark = [&](const RotationIndexerResult &r, int id) { - IndexAndRefine idx(x, indexer_pool.get(), /*retain_outcomes=*/false); - idx.ForceRotationIndexerResult(r); - idx.IntegrateEveryFrame(true); - std::mutex m_assign; - ParallelFor(static_cast(census_frames.size()), - std::min(census_frames.size(), config_.nthreads), [&](int i) { - const int f = census_frames[i]; - DataMessage m{}; - m.number = f; - m.spots = spot_cache.at(f); - idx.CountSpotsOnLattice(m, validation_settings); - std::unique_lock ul(m_assign); - auto &a = assign[f]; - if (m.spots.size() != a.size()) - return; - for (size_t k = 0; k < a.size(); k++) - if (a[k] == -1 && m.spots[k].indexed && m.spots[k].lattice == 0) - a[k] = id; - }); - }; - // Leftover validation spots on `r` at the true spindle angles, and the median of five - // displaced ones - the first pass's null (pooled_evidence) on the leftovers. - auto leftover_evidence = [&](const RotationIndexerResult &r) -> std::pair { - IndexAndRefine idx(x, indexer_pool.get(), /*retain_outcomes=*/false); - idx.ForceRotationIndexerResult(r); - const int nv = static_cast(validation.size()); - auto at_shift = [&](int shift) { - std::atomic n = 0; - ParallelFor(nv, std::min(validation.size(), config_.nthreads), [&](int i) { - DataMessage m{}; - m.number = validation[(i + shift) % nv]; - const auto &sp = spot_cache.at(validation[i]); - const auto &a = assign.at(validation[i]); - for (size_t k = 0; k < sp.size(); k++) - if (a[k] == -1) - m.spots.push_back(sp[k]); - n += idx.CountSpotsOnLattice(m, validation_settings); - }); - return n.load(); - }; - std::vector chance; - for (int k = 1; k <= 5; k++) - chance.push_back(at_shift(std::max(1, k * nv / 7))); - std::ranges::sort(chance); - return {at_shift(0), chance[2]}; - }; - // Spots and summed intensity assigned `id`, and the share of them in the 2 densest of 8 - // equal blocks of the census frames (0.25 is uniform over the sweep). - struct Tally { int64_t n = 0; double intensity = 0.0; double concentration = 0.0; }; - auto tally = [&](int id) { - Tally t; - std::array block{}; - const size_t nf = census_frames.size(); - for (size_t i = 0; i < nf; i++) { - const int f = census_frames[i]; - const auto &a = assign.at(f); - for (size_t k = 0; k < a.size(); k++) - if (a[k] == id) { - t.n++; - t.intensity += spot_cache.at(f)[k].intensity; - // numpy.array_split: the first nf % 8 blocks hold one frame more - const size_t q = nf / 8, r = nf % 8; - block[i < r * (q + 1) ? i / (q + 1) : r + (i - r * (q + 1)) / std::max(1, q)]++; - } - } - std::ranges::sort(block); - t.concentration = static_cast(block[6] + block[7]) / std::max(1, t.n); - return t; - }; - auto reduced = [](const RotationIndexerResult &r) { - return LatticeSearch(r.lattice.ToPrimitive(r.search_result.centering)).primitive_reduced; - }; - - mark(lattice, 0); - int64_t n_all = 0; - double i_all = 0.0; - for (const auto &[f, a] : assign) - for (size_t k = 0; k < a.size(); k++) - if (a[k] != -2) { - n_all++; - i_all += spot_cache.at(f)[k].intensity; - } - const Tally t0 = tally(0); - LeftoverLatticeCensus census; - census.frames = static_cast(census_frames.size()); - census.primary_spots_pct = 100.0 * t0.n / std::max(1, n_all); - census.primary_intensity_pct = 100.0 * t0.intensity / std::max(1.0, i_all); - const CrystalLattice primary_reduced = reduced(lattice); - for (int id = 1; id <= 3; id++) { - RotationIndexer ri(x, *indexer_pool); - ri.VerifyMemo(config_.verify_first_pass_memo); - for (const int f : spread_ordinals) { - std::vector left; - const auto &sp = spot_cache.at(f); - const auto &a = assign.at(f); - for (size_t k = 0; k < sp.size(); k++) - if (a[k] == -1) - left.push_back(sp[k]); - ri.ProcessImage(f, left, rot_angle(f)); - } - ri.RunIndexing(); - const auto r = ri.GetLattice(); - if (!r) - break; - const auto [real, chance] = leftover_evidence(*r); - mark(*r, id); - const Tally t = tally(id); - const CrystalLattice r_reduced = reduced(*r); - LeftoverLattice l; - l.cell = r->search_result.conventional.GetUnitCell(); - l.centering = r->search_result.centering; - l.system = gemmi::crystal_system_str(r->search_result.system); - l.volume_ratio = std::abs(r_reduced.CalcVolume()) / std::max(1.0f, std::abs(primary_reduced.CalcVolume())); - l.misorientation_deg = LatticeMisorientationDeg(primary_reduced, r_reduced); - l.z = (real - chance) / std::sqrt(std::max(1.0, chance)); - l.spots_pct = 100.0 * t.n / std::max(1, n_all); - l.intensity_pct = 100.0 * t.intensity / std::max(1.0, i_all); - l.sweep_concentration = t.concentration; - l.kind = ClassifyLeftoverLattice(l.misorientation_deg, l.volume_ratio, l.sweep_concentration); - logger.Info("Leftover lattice {}: {}-centred {} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}, " - "volume ratio {:.3f}, misorientation {:.2f} deg, z {:.1f} ({} vs {} at a wrong " - "angle), {:.1f}% of non-ice spots, {:.1f}% of their intensity, sweep " - "concentration {:.2f} -> {}", - id, l.centering, l.system, l.cell.a, l.cell.b, l.cell.c, l.cell.alpha, l.cell.beta, - l.cell.gamma, l.volume_ratio, l.misorientation_deg, l.z, real, chance, l.spots_pct, - l.intensity_pct, l.sweep_concentration, LeftoverLatticeKindName(l.kind)); - if (l.z >= 5.0 && l.kind != LeftoverLatticeKind::Self) - census.lattices.push_back(l); - if (real - chance < 3.0 * std::sqrt(std::max(1.0, chance)) + 5) - break; - } - logger.Info("Leftover-lattice census: the crystal's lattice takes {:.1f}% of {} non-ice spots " - "({:.1f}% of their intensity) on {} frames; {} further lattice(s); {:.1f} s", - census.primary_spots_pct, n_all, census.primary_intensity_pct, census.frames, - census.lattices.size(), - std::chrono::duration(std::chrono::steady_clock::now() - census_start).count()); - return census; - }); - } - if (best.result->tilt_walk) { - // What the walked tilt does to the prediction: the keystone at the far corner of the - // detector, against the disc the integration sums over, and the whole-pattern shift - // it is equivalent to, which the beam centre carries instead where it is refused. - const auto &tw = *best.result->tilt_walk; - const auto &g = best.result->geom; - const double start_rot1 = tw.refused ? g.GetPoniRot1_rad() : experiment_.GetPoniRot1_rad(); - const double start_rot2 = tw.refused ? g.GetPoniRot2_rad() : experiment_.GetPoniRot2_rad(); - const double walk_rad = std::hypot(tw.tilt_rad[0] - start_rot1, tw.tilt_rad[1] - start_rot2); - const double corner_px = std::hypot( - std::max(g.GetBeamX_pxl(), experiment_.GetXPixelsNumConv() - g.GetBeamX_pxl()), - std::max(g.GetBeamY_pxl(), experiment_.GetYPixelsNumConv() - g.GetBeamY_pxl())); - const double lever_px = g.GetDetectorDistance_mm() / g.GetPixelSize_mm(); - const std::string what = fmt::format( - "Detector tilt: the free fit walked to {:.4f},{:.4f} deg, {:.2f} deg from the " - "{:.4f},{:.4f} deg this pass started at - further than a mounted detector is off " - "square by - which moves the far corner of the detector by {:.0f} px against an " - "integration disc of {:.1f} px, so the lattice was refined again at the starting " - "tilt with the beam centre free to take the {:.0f} px shift the tilt is equivalent " - "to, and the two were judged on the spots they index: {:.0f} at the walked tilt " - "against {:.0f} at the starting one", - tw.tilt_rad[0] * 180.0 / PI, tw.tilt_rad[1] * 180.0 / PI, walk_rad * 180.0 / PI, - start_rot1 * 180.0 / PI, start_rot2 * 180.0 / PI, - walk_rad * corner_px * corner_px / lever_px, - experiment_.GetBraggIntegrationSettings().GetR1(), walk_rad * lever_px, - tw.spots_walked, tw.spots_held); - if (tw.refused) - logger.Warning("{} - the spots cannot tell the two apart, the walked tilt measured " - "nothing, and the pass runs at the tilt it started at", what); - else - logger.Warning("{} - the spots confirm the walked tilt, and the pass runs at it; the " - "file misdescribes this detector's tilt", what); - } - - // The spots in hand were already cut to the budget in force, so the measurement can only - // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the - // list the first one started from rather than the list the first one settled on. - if (config_.measure_spot_budget) { - const int64_t budget = measure_budget(*indexer, *best.result); - // The measurement can honestly land under MIN_SPOT_COUNT, a floor written for a - // stills FIT, where fewer spots cannot pin a lattice down. Here the lattice is - // already in hand and merely tested - but the validation gate still needs - // ViableCellMinSpots indexed spots per frame (AnalyzeIndexing), so a shorter list - // could never index anything: keep the floor's worth of spots rather than die on - // the setter. A budget that shallow is a finding in its own right - the lattice and - // the file's geometry explain only the brightest few detections (measured on a - // header detector distance 20% off, compensated by a uniformly scaled cell that - // only matches near the beam centre) - so it is warned about, not just adopted. - const int64_t adopted = std::max(budget, MIN_SPOT_COUNT); - if (budget > 0 && adopted < experiment_.GetMaxSpotCount()) { - if (adopted != budget) - logger.Warning("Spot budget: only the strongest {} spots per image lie on " - "the lattice (measured on {} frames) - fewer than the {} the " - "frame gate can act on, so {} are kept. A budget this shallow " - "usually means the geometry in the file (detector distance, " - "beam centre) or the lattice is wrong, not that the crystal " - "is weak", budget, static_cast(validation.size()), - MIN_SPOT_COUNT, adopted); - else - logger.Info("Spot budget: keeping the strongest {} spots per image " - "(measured on {} frames; was {})", budget, - static_cast(validation.size()), experiment_.GetMaxSpotCount()); - experiment_.MaxSpotCount(adopted); - } - } - - // Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image - // with it. A markedly larger cell is the bistable supercell collapse, which was already - // rejected in Run() - but only after the whole pass had run, so it cost an entire extra pass - // on a dataset that ended up on pass-1's lattice anyway. Take pass-1's result now instead: - // the same remedy, one pass earlier. - if (!geometry_prepass && prepass_result_.has_value()) { - // Primitive volumes, like the scheme comparison above: a centred conventional cell is an - // exact integer multiple of its primitive one, so two settings of the SAME lattice differ - // by that factor and comparing the conventional cells reads a mere change of setting as a - // supercell - which then forces pass 1's setting and, with it, its lower symmetry. - const double v1 = std::abs(prepass_result_->lattice - .ToPrimitive(prepass_result_->search_result.centering).CalcVolume()); - const double v2 = std::abs(best.result->lattice - .ToPrimitive(best.result->search_result.centering).CalcVolume()); - const bool supercell = v1 > 1.0 && v2 > 1.5 * v1; - // Dropping to the primitive sub-cell of pass-1's centred lattice is the same - // disagreement seen from the other side, and the volumes above cannot see it: the two - // have the SAME primitive metric. Integrating there never predicts the reflections the - // centring implies, and the pass then merges in P1 (measured on a C-centred monoclinic - // crystal). Only this direction - a centring pass 2 finds where pass 1 had none is the - // refined geometry doing its job, and the intensity test in the search confirms it. - // Same primitive metric only: a pass 2 that lands on a primitive cell several times - // smaller is a different lattice - pass 1's supercell corrected - not a dropped centring. - const bool lost_centring = prepass_result_->search_result.centering != 'P' - && best.result->search_result.centering == 'P' - && v1 > 1.0 && std::abs(v2 / v1 - 1.0) < 0.1; - if (supercell || lost_centring) { - logger.Info("Two-pass: that lattice ({}-centred, volume {:.0f} A^3) disagrees with " - "pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead", - best.result->search_result.centering, v2, - prepass_result_->search_result.centering, v1); - RotationIndexerResult forced = *prepass_result_; - if (prepass_rotation_scale_ && forced.axis) - forced.axis = ScaleRotation(*forced.axis, *prepass_rotation_scale_); - indexer->ForceRotationIndexerResult(forced); - } - } - } - return false; +auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass { + return PickBest(*indexer_pool, *indexer, [&] { + StartCentreSpeculation(); + StartShortAxisSpeculation(); + }); } diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h new file mode 100644 index 000000000..a268c97eb --- /dev/null +++ b/rugnux/RugnuxFirstPass.h @@ -0,0 +1,265 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// Private to the Rugnux implementation files, like RugnuxPipeline.h: the first pass of a rotation +// run (Rugnux::FirstPassRotationIndexing) as the state its steps share - the spot cache, the +// validation frames, the sampling schemes, the spot-list settings the run stands on - and one +// member function per step. FirstPassRotationIndexing itself is the ordered list of the steps; the +// standing lattice (FirstPass) goes in and comes out of each of them. +// +// RugnuxFirstPass.cpp: the shared machinery (spots, scoring, schemes, the ladder) and the order. +// RugnuxFirstPassRescues.cpp: the rescues and hypotheses that may replace the standing lattice. +// RugnuxFirstPassAccept.cpp: the no-crystal test and what the run takes from the accepted lattice. + +#include +#include +#include +#include +#include +#include + +#include "RugnuxPipeline.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" + +namespace rugnux_internal { + // A candidate lattice of the first pass and what it scored on the validation frames. + // indexer_error: what the indexer threw, if it threw. An empty result with no error is a + // verdict on the crystal; an empty result with one is a verdict on the machine. + struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; + std::optional indexer_error; }; +} + +struct Rugnux::FirstPassRun { + FirstPassRun(Rugnux &rx, PipelineLocals &p); + + using FirstPass = rugnux_internal::FirstPass; + using SchemeIndexers = std::vector>; + using Schemes = std::vector>>; + // The pooled evidence for one candidate lattice: the validation frames' non-ice spots, how many + // of them lie on the candidate, and how many a wrong spindle angle still puts there. + using PooledEvidence = ValidationSpotEvidence; + + // --- The steps, in the order FirstPassRotationIndexing takes them --- + + // The spots of every frame the pass is certain to want, in one batch. + void PrefetchCertainFrames(); + // The schemes, indexed and scored, with the speculative passes started beside them. + FirstPass IndexSchemes(); + FirstPass AxisSignRescue(FirstPass best); + // The beam-centre check; sets long_axis_asked where it has asked the long-axis rescue itself. + FirstPass BeamCentreCheck(FirstPass best, bool &long_axis_asked); + std::optional LongAxisRescue(const FirstPass &standing); + FirstPass AddEndWedge(FirstPass bp, IndexerThreadPool &pool, IndexAndRefine &idx); + void ReleaseSpotEngines(); + // Throws where the pass ends with no lattice at all. + void RequireLattice(const FirstPass &best) const; + FirstPass ShortAxisHypothesis(FirstPass best); + FirstPass SubLatticeHypothesis(FirstPass best); + FirstPass Pass1Fallback(FirstPass best); + std::optional FirstPassLadder(int score_before, int min_gain, int min_score); + + // The pooled spots of the validation frames on the lattice the pass ended on, and their share. + struct Verdict { + PooledEvidence evidence; + double pooled = 0.0; + double pooled_chance = 0.0; + }; + Verdict MeasureVerdict(const FirstPass &best, const std::vector &first_pass_key, + const std::optional &memo_hit); + // Throws where the lattice does not beat chance. + void RejectChanceLattice(const FirstPass &best, const Verdict &v) const; + void AdoptLattice(const FirstPass &best, const Verdict &v); + void StartLeftoverCensus(const FirstPass &best); + void ReportTiltWalk(const FirstPass &best); + void MeasureSpotBudget(const FirstPass &best); + void Pass1SupercellGuard(const FirstPass &best); + + // --- The beam-centre check, by case --- + + // The two centres the check compares. + struct Centres { + BeamCenterEstimate measured; + float header_x = 0.0f, header_y = 0.0f; + float measured_x = 0.0f, measured_y = 0.0f; + float moved = 0.0f; + }; + FirstPass IndexAtMeasuredCentre(); + FirstPass CentreWhenFileFails(FirstPass best, FirstPass alt, const Centres &c, bool &long_axis_asked); + FirstPass CentreByPooledSpots(FirstPass best, FirstPass alt, const Centres &c, bool &long_axis_asked); + FirstPass CentreWhenBothIndex(FirstPass best, const FirstPass &alt, const Centres &c); + FirstPass CentresDisagree(FirstPass best, const FirstPass &alt, const Centres &c, double ratio); + + // --- The machinery the steps share --- + + std::optional RotAngle(int ordinal) const; + std::vector FindSpots(MXAnalysisWithoutFPGA &analysis, AzimuthalIntegrationProfile &profile, + JFJochReaderRawImage &img, int ordinal); + static void UseGpuOfWorker(size_t t); + void PrefetchSpots(const std::vector &ordinals); + void TryBeamCenter(float x, float y); + void RestoreBeamCenter(float x, float y); + std::vector ValidationCloud(const RotationIndexerResult &r); + int CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r); + std::pair PooledOnLattice(IndexAndRefine &idx, const RotationIndexerResult &r, int shift); + PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r); + int64_t MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r); + std::unique_ptr FeedScheme(IndexerThreadPool &pool, const std::vector &ordinals, + const DiffractionExperiment *x = nullptr); + SchemeIndexers FeedSchemes(IndexerThreadPool &pool, const DiffractionExperiment *x = nullptr); + static void RunIndexing(const std::vector &ris); + static void RunIndexing(const SchemeIndexers &ris); + FirstPass ScoreSchemes(const SchemeIndexers &ris, IndexAndRefine &idx, const Schemes *named = nullptr, + FirstPass bp = {}); + FirstPass PickBest(IndexerThreadPool &pool, IndexAndRefine &idx, + const std::function &alongside = {}); + void ApplyState(const FirstPassState &st); + IndexingSettings ShortAxisSettings(const DiffractionExperiment &x) const; + void StartCentreSpeculation(); + void StartShortAxisSpeculation(); + + // --- State --- + + Rugnux ℞ + Logger &logger; + ProcessResult &result; + const bool geometry_prepass; + const int start_image; + // Every frame the first pass reads is in the first part of the sweep (Rugnux::PrepassEnd). + const int early_images; + const std::shared_ptr &dataset; + AzimuthalIntegrationMapping &mapping; + std::unique_ptr &indexer_pool; + std::unique_ptr &indexer; + const float indexing_bound_A; + std::future &census_future; + + // Spots for a first-pass image ordinal, cached (a frame is reused across schemes and the + // validation set; re-finding is expensive on the --redo-rotation-spots path). The workers below + // hand their spots back and this thread puts them in, so the map itself is still only ever + // touched from here and needs no locking. + std::map> spot_cache; + // The azimuthal mapping the first pass finds its spots against. It is `mapping` until a + // beam-centre hypothesis moves the centre, and the hypothesis's own mapping while it stands - + // the mapping is where the resolution of every pixel and the finder's radial background bins + // come from, and both move with the centre. + const AzimuthalIntegrationMapping *spot_mapping = &mapping; + + // The workers PrefetchSpots finds spots on, each with an analysis engine it keeps (see there). + const int spot_workers; + struct SpotEngine { + // Before the analysis engine, which page-locks these bytes for its uploads and + // unregisters them when it is destroyed: it must not outlive the buffer it was given. + JFJochReaderRawImage raw_image; + std::unique_ptr analysis; + std::unique_ptr profile; + }; + std::vector engines; + int prefetch_batches = 0; + double prefetch_time_s = 0.0; + + // Evaluating a beam-centre hypothesis means re-finding the spots at that centre. The spot + // POSITIONS do not move with the centre - which is why the cache is keyed by image alone - but + // the spot LIST does: every spot's d, the ice-ring flag, the resolution mask the finder + // applies, the high-resolution gap cut and the strongest-N budget are all radial, and the + // adaptive finder's ring background is binned about the centre. Scoring a trial centre against + // spots found at the file's puts the trial in a world the run never enters: measured, a centre + // scored 42/60 as a trial and 5/60 when that same value was given on the command line. + // + // A trial therefore gets its own azimuthal mapping, its own engines and its own cache, and the + // starting centre's spots are PARKED rather than thrown away - a hypothesis that is not adopted + // has to leave the rest of the run exactly as it found it. + std::unique_ptr trial_mapping; + std::map> home_spot_cache; + + // How many validation frames a candidate global lattice actually indexes per image, scored on + // the *real* per-image path (force the candidate, then run the same indexing+refinement the + // main loop uses). A cell that merely fits the accumulated first-pass cloud (a spurious + // sub-lattice from an under-sampled scheme) indexes few frames, while the true cell indexes + // many - this is the discriminator that lets the best scheme win. + const std::vector validation = rugnux_internal::select_equally_spaced_image_ordinals( + early_images, std::min(early_images, 60)); + + const float index_tol = rx.experiment_.GetIndexingSettings().GetTolerance(); + + // The run's spot-finding settings, with indexing on: what a validation frame is indexed under. + SpotFindingSettings validation_settings; + + // How many spot-finding keys first_pass_spots_ holds; a handful covers a run, the oldest go first. + static constexpr size_t MAX_SPOT_KEYS = 8; + // How far ahead of itself FeedScheme keeps the spot cache. + static constexpr int SPOT_PREFETCH_CHUNK = 64; + size_t first_pass_spots_per_image = 0; // 0 = the indexer's own default + + // Sampling schemes, tried in order of preference. We index each and keep the one that indexes + // the most frames; the schemes cover complementary failure modes and no single one is best for + // every crystal: + // - "spread": images spread over the first part of the rotation (PrepassEnd) - wide angular range, + // coarse stride; the long-standing default, best for typical crystals. + // - "wedge": consecutive frames from the clean start - native stride keeps each reflection's + // rocking curve continuous across frames, so the FFT can trace a long axis whose fine + // reciprocal spacing the coarse spread cannot resolve (e.g. an F-cubic cell viewed down + // [111]); bounded by the accumulation cap to the least radiation-damaged early wedge. + std::vector all_consecutive; + std::vector spread_ordinals; + Schemes schemes; + // A third sample of the same frames, asked only when nothing else indexes a majority (see + // AddEndWedge): the wedge read back from the END of the first part. Which orientations a sweep + // passes through first is chance, and the first part can open on frames that see the long axis + // badly - measured on a sparse crystal with a ~390 A axis, the spread and the wedge over its + // first 120 deg returned wrong cells (5/60 and 0/60 validation frames) at every spread density + // tried, while consecutive frames from 60 deg on indexed 60/60. Still inside [start, + // PrepassEnd), so nothing before pass 2 waits for the rest of the sweep. + Schemes end_wedge; + + // What the run is CURRENTLY reading each frame at. It starts at what the file asked for, and + // only an ADOPTED rung moves it. Every rung the ladder tries is put back to this, so a rung + // tried and rejected leaves nothing behind - and a later pass that finds no rung restores what + // an earlier pass adopted rather than resetting the run to the defaults. + const std::optional first_pass_file_d_min; + FirstPassState committed; + + // The short-axis pass below indexes the same schemes at a lower floor, after the rescues in + // between - and on nearly every run they change nothing it reads. So it is started here, on a + // copy of the experiment, beside the standard indexing it would otherwise wait for, and taken + // at its place only if what it read is still what the run has there (ExperimentKey, the spot + // list's mapping and budget); otherwise it is run again there as before. + struct ShortAxisSpeculation { + std::vector key; + const AzimuthalIntegrationMapping *mapping = nullptr; + size_t spots_per_image = 0; + std::unique_ptr experiment; + std::unique_ptr pool; + std::vector> sris; + std::vector> futures; // last, so it is waited for before the rest goes + }; + std::optional short_spec; + // The beam-centre check below runs a second first pass at the centre measured from the + // background, and it reads the first one's answer only after that pass is over. So the spots at + // the measured centre are found, and both schemes fed and indexed, here beside the file centre's + // pass - on a copy of the experiment at that centre, with workers, a mapping and a spot list of + // their own - and taken at the check's place only if what they read is still what the run has + // there (the experiment with its rotation angles, SpotFindingKey at the measured centre, the spot + // budget); otherwise that pass is run there as before. The schemes are scored there in either + // case, on the run's own indexer. + // + // Only with a second GPU: on one card the two passes share the card and the reader, and measured, + // running them side by side saved nothing. There the check runs its pass after the first, as it + // always did. With two or more its spot finding spreads over all the cards, as the first pass's + // does: the file centre's spot finding is over long before this pass's, so this one is what the + // first pass waits on. + struct CentreSpeculation { + std::vector experiment_key; + std::vector key; + size_t spots_per_image = 0; + std::unique_ptr experiment; + std::unique_ptr mask; + std::unique_ptr mapping; + std::map> spots; + SchemeIndexers ris; + std::future done; // last, so it is waited for before the rest goes + }; + std::optional centre_spec; +}; diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp new file mode 100644 index 000000000..abca5c1bf --- /dev/null +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -0,0 +1,500 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/ModelValidation.h" +#include "../writer/WriteModel.h" +#include "../image_analysis/indexing/SpindleCuspLoss.h" +#include "../image_analysis/bragg_integration/SpotWidth.h" +#include "../image_analysis/bragg_integration/SpotFootprint.h" +#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA +#include "../image_analysis/hot_pixels/HotPixels.h" +#include "../image_analysis/SensorAbsorption.h" +#include "DiagnosticOutput.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../reader/JFJochHDF5Reader.h" +#include "../common/JFJochMath.h" +#include "../common/ParallelFor.h" +#include "../common/Logger.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../common/CUDAWrapper.h" +#include "../common/JFJochException.h" +#include "../common/time_utc.h" +#include "../writer/FileWriter.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" +#include "../image_analysis/indexing/HarmonicContamination.h" +#include "../image_analysis/geom_refinement/BeamCenterFFT.h" +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" +#include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#include "../image_analysis/spot_finding/SpotUtils.h" +#include "gemmi/twin.hpp" +#include "../image_analysis/azint/AzIntEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#endif +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../image_analysis/scale_merge/ScalingResult.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/geom_refinement/PostRefine.h" +#include "../image_analysis/lattice_search/LatticeSearch.h" +#include "../image_analysis/lattice_search/LePageLattice.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" +#include "../image_analysis/scale_merge/TranslationalNCS.h" +#include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/bragg_integration/CalcISigma.h" +#include "../common/Definitions.h" +#include "../common/CorrelationCoefficient.h" +#include +#include +#include "RugnuxPipeline.h" +#include "RugnuxFirstPass.h" + +// The end of the first pass: whether it found a lattice at all, the no-crystal test on the pooled +// validation spots, and what the run takes from the lattice it accepts. FirstPassRotationIndexing +// gives their order. + +using namespace rugnux_internal; + +namespace { + // The two ways the first pass can end without a usable lattice that are about the DATA. + // They are one situation from the user's side - no crystal was found in these frames - so + // both name the cause and then give the same remedies. + constexpr const char *next_steps = + " Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise " + "--max-spots if the frames are well populated with spots, and use --mode azint to " + "write the per-image analysis for inspection."; +} + +void Rugnux::FirstPassRun::RequireLattice(const FirstPass &best) const { + if (!best.result.has_value()) { + // An indexer that FAILED is a third outcome, and none of those remedies fit it: it + // never examined the frames, so nothing about the data or the processing options is + // implicated and changing them cannot help. Name the failure instead. + if (best.indexer_error) + throw JFJochException(JFJochExceptionCategory::IndexerError, + fmt::format("Rotation indexing could not run: {}. Nothing was " + "integrated. The indexer never examined the frames, " + "so this is a failure of the machine and not of the " + "data - free the resource it needs (a GPU shared with " + "another process is the usual cause) and run again.", + *best.indexer_error)); + throw NoLatticeFound(std::string("Two-pass rotation indexing found no lattice.") + next_steps); + } +} + +// Whether the run has a lattice at all is decided here, on the pooled spots rather than on +// the frame count the rescues above are steered by. A frame count asks each frame to clear +// a bar of its own and then counts the frames that did, which measures the spots per frame +// as much as the lattice; the pooled fraction asks the whole validation set at once and is +// judged against what the same lattice scores at a wrong spindle angle, so a sparse pattern +// is refused only when its spots genuinely do not know about the lattice. +auto Rugnux::FirstPassRun::MeasureVerdict(const FirstPass &best, const std::vector &first_pass_key, + const std::optional &memo_hit) -> Verdict { + const PooledEvidence evidence = MeasurePooled(*indexer, *best.result); + const double pooled = static_cast(evidence.on_lattice) + / static_cast(std::max(1, evidence.spots)); + const double pooled_chance = static_cast(evidence.by_chance) + / static_cast(std::max(1, evidence.spots)); + result.validation_evidence = evidence; + if (memo_hit && (memo_hit->spots != evidence.spots || memo_hit->on_lattice != evidence.on_lattice + || memo_hit->by_chance != evidence.by_chance)) + throw MemoMismatch(fmt::format("First-pass memo mismatch: stored {}/{}/{}, recomputed {}/{}/{}", + memo_hit->on_lattice, memo_hit->spots, memo_hit->by_chance, + evidence.on_lattice, evidence.spots, evidence.by_chance)); + // Kept only where the pass ended on the inputs it started from. + if (!geometry_prepass && rx.FirstPassInputKey() == first_pass_key) + rx.AddFirstPassMemo(FirstPassMemo{first_pass_key, evidence}); + return {evidence, pooled, pooled_chance}; +} + +void Rugnux::FirstPassRun::RejectChanceLattice(const FirstPass &best, const Verdict &v) const { + const PooledEvidence &evidence = v.evidence; + const double pooled = v.pooled, pooled_chance = v.pooled_chance; + if (!ValidationEvidenceBeatsChance(evidence)) { + { + // Name the cell and Bravais class that was rejected. The commonest cause is a metric + // symmetry promoted one class too far - the constrained cell then misses every + // reflection by the small angle the constraint snapped away - and without the cell in + // the message there is nothing to see that from. + const auto &c = best.result->search_result.conventional.GetUnitCell(); + throw NoLatticeFound(fmt::format("Two-pass rotation indexing found only a lattice that " + "takes {:.1f}% of the validation spots against the " + "{:.1f}% a wrong spindle angle takes ({} distinct " + "reflections against {} over {} wrong angles) - it is " + "not this crystal's lattice. It was {}-centred {}, " + "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}", + 100.0 * pooled, 100.0 * pooled_chance, + evidence.reflections, evidence.reflections_by_chance, + VALIDATION_NULL_DISPLACEMENTS, + best.result->search_result.centering, + gemmi::crystal_system_str(best.result->search_result.system), + c.a, c.b, c.c, c.alpha, c.beta, c.gamma, next_steps)); + } + } +} + +void Rugnux::FirstPassRun::AdoptLattice(const FirstPass &best, const Verdict &v) { + const double pooled = v.pooled, pooled_chance = v.pooled_chance; + const PooledEvidence &evidence = v.evidence; + indexer->ForceRotationIndexerResult(*best.result); + // A lattice the pooled spots accept while most validation frames fail the per-frame floor + // (LATTICE_MIN_INDEXED_FRACTION) on their own, because the sweep as a whole puts fewer of its + // spots on the lattice than that floor asks of each frame: the crystal is real but its + // spots are a minority of every frame - a weak crystal inside a powder of its own + // microcrystals and ice. The floor cannot tell its frames from anything else, so integrate + // every frame from the sweep's lattice, as XDS and DIALS do, and leave the frames to the + // scaling. The frames that do clear it are the upper tail of the same population: measured + // on a crystal in a powder of itself (13% of the pooled spots on the lattice), 12/60 + // cleared it, and keeping the floor integrated 15% of the sweep. A sweep that is sparse + // only in spots per frame (a small molecule) keeps the floor. + const bool sparse_lattice = best.score < static_cast(validation.size()) / 2 + && pooled < LATTICE_MIN_INDEXED_FRACTION; + indexer->IntegrateEveryFrame(sparse_lattice); + if (sparse_lattice) + logger.Warning("Sparse lattice: {}/{} validation frames clear the per-frame floor, but the " + "pooled spots put {:.1f}% on the lattice against {:.1f}% by chance - every " + "frame is integrated from the sweep's lattice", + best.score, static_cast(validation.size()), 100.0 * pooled, + 100.0 * pooled_chance); + logger.Info("First-pass spot finding: {} frames in {} batches, {:.2f} s", + spot_cache.size(), prefetch_batches, prefetch_time_s); + logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames, " + "{}/{} validation spots = {:.1f}% against {:.1f}% at a wrong spindle angle; {} " + "distinct reflections against {} over {} wrong angles)", + best.name, best.score, static_cast(validation.size()), + evidence.on_lattice, evidence.spots, 100.0 * pooled, 100.0 * pooled_chance, + evidence.reflections, evidence.reflections_by_chance, VALIDATION_NULL_DISPLACEMENTS); + { + const auto &c = best.result->search_result.conventional.GetUnitCell(); + rx.Note(fmt::format("lattice {} {}: {:.2f} {:.2f} {:.2f} {:.1f} {:.1f} {:.1f}", + gemmi::crystal_system_str(best.result->search_result.system), + best.result->search_result.centering, c.a, c.b, c.c, c.alpha, c.beta, c.gamma)); + rx.Note(fmt::format("{}/{} validation frames, {:.0f}% of their spots on it", + best.score, validation.size(), 100.0 * pooled)); + } +} + +// Leftover-lattice census, report only: index what the crystal's lattice leaves over, to name +// a split crystal, a twin domain or a second crystal in the report. Every spot of the spread +// and validation frames the lattice takes is marked; the rest (not on ice) goes to a fresh +// rotation indexer; a lattice it finds is judged on the validation leftovers against the same +// displaced-spindle-angle null as the first pass, its spots are marked in turn, and the search +// goes on for up to three lattices. Its own IndexAndRefine instances and indexer, so nothing +// the run goes on with is touched. Run on the pass that writes the result only. +// +// Its spots are found here; the rest runs in the background, beside the image loop, on +// copies of everything it reads, and its answer is collected before the pass returns. Not on +// a pass that only post-refines: nothing reads its result there. +void Rugnux::FirstPassRun::StartLeftoverCensus(const FirstPass &best) { + std::vector census_frames = spread_ordinals; + census_frames.insert(census_frames.end(), validation.begin(), validation.end()); + std::ranges::sort(census_frames); + census_frames.erase(std::unique(census_frames.begin(), census_frames.end()), census_frames.end()); + PrefetchSpots(census_frames); + std::map> census_spots; + for (const int f : census_frames) + census_spots.emplace(f, spot_cache.at(f)); + census_future = std::async(std::launch::async, + [&logger = logger, &indexer_pool = indexer_pool, &config_ = rx.config_, + x = rx.experiment_, census_frames, census_spots, lattice = *best.result, + validation = validation, spread_ordinals = spread_ordinals, + validation_settings = validation_settings]() -> LeftoverLatticeCensus { + const auto census_start = std::chrono::steady_clock::now(); + const bool ice_ok = x.GetIndexingSettings().GetIndexIceRings(); + const auto &spot_cache = census_spots; + const auto rot_angle = [&](int ordinal) -> std::optional { + if (const auto g = x.GetGoniometer()) + return g->GetAngle_deg(static_cast(ordinal)) + g->GetWedge_deg() / 2.0f; + return std::nullopt; + }; + // assign[frame][spot]: -2 ice, -1 no lattice, 0 the crystal's lattice, 1.. a leftover lattice + std::map> assign; + for (const int f : census_frames) { + auto &a = assign[f]; + a.assign(spot_cache.at(f).size(), -1); + for (size_t k = 0; k < a.size(); k++) + if (!ice_ok && spot_cache.at(f)[k].ice_ring) + a[k] = -2; + } + // Mark the still unassigned spots lattice `r` takes as `id`. + auto mark = [&](const RotationIndexerResult &r, int id) { + IndexAndRefine idx(x, indexer_pool.get(), /*retain_outcomes=*/false); + idx.ForceRotationIndexerResult(r); + idx.IntegrateEveryFrame(true); + std::mutex m_assign; + ParallelFor(static_cast(census_frames.size()), + std::min(census_frames.size(), config_.nthreads), [&](int i) { + const int f = census_frames[i]; + DataMessage m{}; + m.number = f; + m.spots = spot_cache.at(f); + idx.CountSpotsOnLattice(m, validation_settings); + std::unique_lock ul(m_assign); + auto &a = assign[f]; + if (m.spots.size() != a.size()) + return; + for (size_t k = 0; k < a.size(); k++) + if (a[k] == -1 && m.spots[k].indexed && m.spots[k].lattice == 0) + a[k] = id; + }); + }; + // Leftover validation spots on `r` at the true spindle angles, and the median of five + // displaced ones - the first pass's null (MeasurePooled) on the leftovers. + auto leftover_evidence = [&](const RotationIndexerResult &r) -> std::pair { + IndexAndRefine idx(x, indexer_pool.get(), /*retain_outcomes=*/false); + idx.ForceRotationIndexerResult(r); + const int nv = static_cast(validation.size()); + auto at_shift = [&](int shift) { + std::atomic n = 0; + ParallelFor(nv, std::min(validation.size(), config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = validation[(i + shift) % nv]; + const auto &sp = spot_cache.at(validation[i]); + const auto &a = assign.at(validation[i]); + for (size_t k = 0; k < sp.size(); k++) + if (a[k] == -1) + m.spots.push_back(sp[k]); + n += idx.CountSpotsOnLattice(m, validation_settings); + }); + return n.load(); + }; + std::vector chance; + for (int k = 1; k <= 5; k++) + chance.push_back(at_shift(std::max(1, k * nv / 7))); + std::ranges::sort(chance); + return {at_shift(0), chance[2]}; + }; + // Spots and summed intensity assigned `id`, and the share of them in the 2 densest of 8 + // equal blocks of the census frames (0.25 is uniform over the sweep). + struct Tally { int64_t n = 0; double intensity = 0.0; double concentration = 0.0; }; + auto tally = [&](int id) { + Tally t; + std::array block{}; + const size_t nf = census_frames.size(); + for (size_t i = 0; i < nf; i++) { + const int f = census_frames[i]; + const auto &a = assign.at(f); + for (size_t k = 0; k < a.size(); k++) + if (a[k] == id) { + t.n++; + t.intensity += spot_cache.at(f)[k].intensity; + // numpy.array_split: the first nf % 8 blocks hold one frame more + const size_t q = nf / 8, r = nf % 8; + block[i < r * (q + 1) ? i / (q + 1) : r + (i - r * (q + 1)) / std::max(1, q)]++; + } + } + std::ranges::sort(block); + t.concentration = static_cast(block[6] + block[7]) / std::max(1, t.n); + return t; + }; + auto reduced = [](const RotationIndexerResult &r) { + return LatticeSearch(r.lattice.ToPrimitive(r.search_result.centering)).primitive_reduced; + }; + + mark(lattice, 0); + int64_t n_all = 0; + double i_all = 0.0; + for (const auto &[f, a] : assign) + for (size_t k = 0; k < a.size(); k++) + if (a[k] != -2) { + n_all++; + i_all += spot_cache.at(f)[k].intensity; + } + const Tally t0 = tally(0); + LeftoverLatticeCensus census; + census.frames = static_cast(census_frames.size()); + census.primary_spots_pct = 100.0 * t0.n / std::max(1, n_all); + census.primary_intensity_pct = 100.0 * t0.intensity / std::max(1.0, i_all); + const CrystalLattice primary_reduced = reduced(lattice); + for (int id = 1; id <= 3; id++) { + RotationIndexer ri(x, *indexer_pool); + ri.VerifyMemo(config_.verify_first_pass_memo); + for (const int f : spread_ordinals) { + std::vector left; + const auto &sp = spot_cache.at(f); + const auto &a = assign.at(f); + for (size_t k = 0; k < sp.size(); k++) + if (a[k] == -1) + left.push_back(sp[k]); + ri.ProcessImage(f, left, rot_angle(f)); + } + ri.RunIndexing(); + const auto r = ri.GetLattice(); + if (!r) + break; + const auto [real, chance] = leftover_evidence(*r); + mark(*r, id); + const Tally t = tally(id); + const CrystalLattice r_reduced = reduced(*r); + LeftoverLattice l; + l.cell = r->search_result.conventional.GetUnitCell(); + l.centering = r->search_result.centering; + l.system = gemmi::crystal_system_str(r->search_result.system); + l.volume_ratio = std::abs(r_reduced.CalcVolume()) / std::max(1.0f, std::abs(primary_reduced.CalcVolume())); + l.misorientation_deg = LatticeMisorientationDeg(primary_reduced, r_reduced); + l.z = (real - chance) / std::sqrt(std::max(1.0, chance)); + l.spots_pct = 100.0 * t.n / std::max(1, n_all); + l.intensity_pct = 100.0 * t.intensity / std::max(1.0, i_all); + l.sweep_concentration = t.concentration; + l.kind = ClassifyLeftoverLattice(l.misorientation_deg, l.volume_ratio, l.sweep_concentration); + logger.Info("Leftover lattice {}: {}-centred {} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}, " + "volume ratio {:.3f}, misorientation {:.2f} deg, z {:.1f} ({} vs {} at a wrong " + "angle), {:.1f}% of non-ice spots, {:.1f}% of their intensity, sweep " + "concentration {:.2f} -> {}", + id, l.centering, l.system, l.cell.a, l.cell.b, l.cell.c, l.cell.alpha, l.cell.beta, + l.cell.gamma, l.volume_ratio, l.misorientation_deg, l.z, real, chance, l.spots_pct, + l.intensity_pct, l.sweep_concentration, LeftoverLatticeKindName(l.kind)); + if (l.z >= 5.0 && l.kind != LeftoverLatticeKind::Self) + census.lattices.push_back(l); + if (real - chance < 3.0 * std::sqrt(std::max(1.0, chance)) + 5) + break; + } + logger.Info("Leftover-lattice census: the crystal's lattice takes {:.1f}% of {} non-ice spots " + "({:.1f}% of their intensity) on {} frames; {} further lattice(s); {:.1f} s", + census.primary_spots_pct, n_all, census.primary_intensity_pct, census.frames, + census.lattices.size(), + std::chrono::duration(std::chrono::steady_clock::now() - census_start).count()); + return census; + }); +} + +void Rugnux::FirstPassRun::ReportTiltWalk(const FirstPass &best) { + // What the walked tilt does to the prediction: the keystone at the far corner of the + // detector, against the disc the integration sums over, and the whole-pattern shift + // it is equivalent to, which the beam centre carries instead where it is refused. + const auto &tw = *best.result->tilt_walk; + const auto &g = best.result->geom; + const double start_rot1 = tw.refused ? g.GetPoniRot1_rad() : rx.experiment_.GetPoniRot1_rad(); + const double start_rot2 = tw.refused ? g.GetPoniRot2_rad() : rx.experiment_.GetPoniRot2_rad(); + const double walk_rad = std::hypot(tw.tilt_rad[0] - start_rot1, tw.tilt_rad[1] - start_rot2); + const double corner_px = std::hypot( + std::max(g.GetBeamX_pxl(), rx.experiment_.GetXPixelsNumConv() - g.GetBeamX_pxl()), + std::max(g.GetBeamY_pxl(), rx.experiment_.GetYPixelsNumConv() - g.GetBeamY_pxl())); + const double lever_px = g.GetDetectorDistance_mm() / g.GetPixelSize_mm(); + const std::string what = fmt::format( + "Detector tilt: the free fit walked to {:.4f},{:.4f} deg, {:.2f} deg from the " + "{:.4f},{:.4f} deg this pass started at - further than a mounted detector is off " + "square by - which moves the far corner of the detector by {:.0f} px against an " + "integration disc of {:.1f} px, so the lattice was refined again at the starting " + "tilt with the beam centre free to take the {:.0f} px shift the tilt is equivalent " + "to, and the two were judged on the spots they index: {:.0f} at the walked tilt " + "against {:.0f} at the starting one", + tw.tilt_rad[0] * 180.0 / PI, tw.tilt_rad[1] * 180.0 / PI, walk_rad * 180.0 / PI, + start_rot1 * 180.0 / PI, start_rot2 * 180.0 / PI, + walk_rad * corner_px * corner_px / lever_px, + rx.experiment_.GetBraggIntegrationSettings().GetR1(), walk_rad * lever_px, + tw.spots_walked, tw.spots_held); + if (tw.refused) + logger.Warning("{} - the spots cannot tell the two apart, the walked tilt measured " + "nothing, and the pass runs at the tilt it started at", what); + else + logger.Warning("{} - the spots confirm the walked tilt, and the pass runs at it; the " + "file misdescribes this detector's tilt", what); +} + +void Rugnux::FirstPassRun::MeasureSpotBudget(const FirstPass &best) { + // The spots in hand were already cut to the budget in force, so the measurement can only + // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the + // list the first one started from rather than the list the first one settled on. + const int64_t budget = MeasureBudget(*indexer, *best.result); + // The measurement can honestly land under MIN_SPOT_COUNT, a floor written for a + // stills FIT, where fewer spots cannot pin a lattice down. Here the lattice is + // already in hand and merely tested - but the validation gate still needs + // ViableCellMinSpots indexed spots per frame (AnalyzeIndexing), so a shorter list + // could never index anything: keep the floor's worth of spots rather than die on + // the setter. A budget that shallow is a finding in its own right - the lattice and + // the file's geometry explain only the brightest few detections (measured on a + // header detector distance 20% off, compensated by a uniformly scaled cell that + // only matches near the beam centre) - so it is warned about, not just adopted. + const int64_t adopted = std::max(budget, MIN_SPOT_COUNT); + if (budget > 0 && adopted < rx.experiment_.GetMaxSpotCount()) { + if (adopted != budget) + logger.Warning("Spot budget: only the strongest {} spots per image lie on " + "the lattice (measured on {} frames) - fewer than the {} the " + "frame gate can act on, so {} are kept. A budget this shallow " + "usually means the geometry in the file (detector distance, " + "beam centre) or the lattice is wrong, not that the crystal " + "is weak", budget, static_cast(validation.size()), + MIN_SPOT_COUNT, adopted); + else + logger.Info("Spot budget: keeping the strongest {} spots per image " + "(measured on {} frames; was {})", budget, + static_cast(validation.size()), rx.experiment_.GetMaxSpotCount()); + rx.experiment_.MaxSpotCount(adopted); + } +} + +// Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image +// with it. A markedly larger cell is the bistable supercell collapse, which was already +// rejected in Run() - but only after the whole pass had run, so it cost an entire extra pass +// on a dataset that ended up on pass-1's lattice anyway. Take pass-1's result now instead: +// the same remedy, one pass earlier. +void Rugnux::FirstPassRun::Pass1SupercellGuard(const FirstPass &best) { + // Primitive volumes, like the scheme comparison above: a centred conventional cell is an + // exact integer multiple of its primitive one, so two settings of the SAME lattice differ + // by that factor and comparing the conventional cells reads a mere change of setting as a + // supercell - which then forces pass 1's setting and, with it, its lower symmetry. + const double v1 = std::abs(rx.prepass_result_->lattice + .ToPrimitive(rx.prepass_result_->search_result.centering).CalcVolume()); + const double v2 = std::abs(best.result->lattice + .ToPrimitive(best.result->search_result.centering).CalcVolume()); + const bool supercell = v1 > 1.0 && v2 > 1.5 * v1; + // Dropping to the primitive sub-cell of pass-1's centred lattice is the same + // disagreement seen from the other side, and the volumes above cannot see it: the two + // have the SAME primitive metric. Integrating there never predicts the reflections the + // centring implies, and the pass then merges in P1 (measured on a C-centred monoclinic + // crystal). Only this direction - a centring pass 2 finds where pass 1 had none is the + // refined geometry doing its job, and the intensity test in the search confirms it. + // Same primitive metric only: a pass 2 that lands on a primitive cell several times + // smaller is a different lattice - pass 1's supercell corrected - not a dropped centring. + const bool lost_centring = rx.prepass_result_->search_result.centering != 'P' + && best.result->search_result.centering == 'P' + && v1 > 1.0 && std::abs(v2 / v1 - 1.0) < 0.1; + if (supercell || lost_centring) { + logger.Info("Two-pass: that lattice ({}-centred, volume {:.0f} A^3) disagrees with " + "pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead", + best.result->search_result.centering, v2, + rx.prepass_result_->search_result.centering, v1); + RotationIndexerResult forced = *rx.prepass_result_; + if (rx.prepass_rotation_scale_ && forced.axis) + forced.axis = ScaleRotation(*forced.axis, *rx.prepass_rotation_scale_); + indexer->ForceRotationIndexerResult(forced); + } +} diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp new file mode 100644 index 000000000..b6f2872ce --- /dev/null +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -0,0 +1,1045 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "Rugnux.h" +#include "../image_analysis/structure_refinement/ModelValidation.h" +#include "../writer/WriteModel.h" +#include "../image_analysis/indexing/SpindleCuspLoss.h" +#include "../image_analysis/bragg_integration/SpotWidth.h" +#include "../image_analysis/bragg_integration/SpotFootprint.h" +#include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA +#include "../image_analysis/hot_pixels/HotPixels.h" +#include "../image_analysis/SensorAbsorption.h" +#include "DiagnosticOutput.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../reader/JFJochHDF5Reader.h" +#include "../common/JFJochMath.h" +#include "../common/ParallelFor.h" +#include "../common/Logger.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" +#include "../common/CUDAWrapper.h" +#include "../common/JFJochException.h" +#include "../common/time_utc.h" +#include "../writer/FileWriter.h" +#include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" +#include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/indexing/AnalyzeIndexing.h" +#include "../image_analysis/indexing/HarmonicContamination.h" +#include "../image_analysis/geom_refinement/BeamCenterFFT.h" +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../image_analysis/geom_refinement/GeometryRefiner.h" +#include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#include "../image_analysis/spot_finding/SpotUtils.h" +#include "gemmi/twin.hpp" +#include "../image_analysis/azint/AzIntEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorGPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#endif +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" +#include "../image_analysis/scale_merge/ResolutionCutoff.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" +#include "../image_analysis/scale_merge/ScalingResult.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/geom_refinement/PostRefine.h" +#include "../image_analysis/lattice_search/LatticeSearch.h" +#include "../image_analysis/lattice_search/LePageLattice.h" +#include "../image_analysis/scale_merge/AnisotropyAnalysis.h" +#include "../image_analysis/scale_merge/FrenchWilson.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" +#include "../image_analysis/scale_merge/TranslationalNCS.h" +#include "../image_analysis/scale_merge/HKLKey.h" +#include "../image_analysis/scale_merge/ScaleOnTheFly.h" +#include "../image_analysis/scale_merge/StillsPartialityRefine.h" +#include "../image_analysis/WriteReflections.h" +#include "../image_analysis/bragg_integration/CalcISigma.h" +#include "../common/Definitions.h" +#include "../common/CorrelationCoefficient.h" +#include +#include +#include "RugnuxPipeline.h" +#include "RugnuxFirstPass.h" + +// The steps of the first pass that may replace the standing lattice: the rescues, asked where it +// indexes too little, and the hypotheses, asked of every run. FirstPassRotationIndexing gives their +// order. + +using namespace rugnux_internal; + +namespace { + auto MaxAxis(const RotationIndexerResult &r) { + const auto uc = r.lattice.GetUnitCell(); + return std::max({uc.a, uc.b, uc.c}); + } +} + +// Rotation-axis SIGN rescue. The sign is a convention the input cannot settle: a PILATUS +// miniCBF header carries only "Oscillation_axis omega" and no direction at all, and facilities +// disagree about the sign of the NXmx vector (measured across ten of them: most want the +// stored [-1,0,0], APS 19-ID and SPring-8 BL41XU want the opposite). Nothing downstream can +// repair it either - the default spindle fit refines the axis DIRECTION by fractions of a +// degree and cannot flip it, that being a discrete change. +// +// It is, however, decidable from the data: the wrong sign does not index at all. Measured on +// one dataset by changing nothing but the sign: 100% of frames indexed one way, 0 of 60 the +// other. So try the other sign and keep whichever indexes more frames - the same +// validation-frame count the scheme choice above already uses, so no new metric and no new +// threshold. Runs before the long-axis rescue: with the sign wrong every candidate lattice is +// wrong, and the coarse pass below would spend itself on a mis-signed cloud. +// +// Only after a poor pass, so a correctly-signed file costs nothing; and spot finding is not +auto Rugnux::FirstPassRun::AxisSignRescue(FirstPass best) -> FirstPass { + if (!rx.cancelled_ && best.score < 0.5 * static_cast(validation.size())) { + if (const auto gon = rx.experiment_.GetGoniometer(); gon.has_value() && gon->IsScanning()) { + GoniometerAxis flipped = *gon; + flipped.Axis(-gon->GetAxis()); + rx.experiment_.Goniometer(flipped); + FirstPass alt = PickBest(*indexer_pool, *indexer); + // More frames is not enough on its own while both counts sit at the floor: 1/60 against + // 0/60 is noise, and adopting it puts every rescue below on the wrong sign (measured: a + // run whose file sign, kept, went on to index found nothing once flipped on 1/60). So + // the flipped lattice must also beat chance on the pooled spots - the test the run puts + // its final lattice to - which a sparse pattern at the right sign passes on few frames + // (2/60 frames, 89% of the validation spots against 0% by chance). + const bool alt_real = alt.result.has_value() && alt.score > best.score + && ValidationEvidenceBeatsChance(MeasurePooled(*indexer, *alt.result)); + if (alt_real) { + logger.Info("Rotation axis sign: the file's axis indexes {}/{} validation frames, " + "the opposite sign {}/{} - adopting the opposite sign for this run", + best.score, static_cast(validation.size()), + alt.score, static_cast(validation.size())); + best = alt; + } else { + if (alt.result.has_value() && alt.score > best.score) + logger.Info("Rotation axis sign: the opposite sign indexes {}/{} validation frames " + "against the file's {}/{}, but its lattice does not beat chance on the " + "pooled spots - the file's sign is kept", + alt.score, static_cast(validation.size()), + best.score, static_cast(validation.size())); + rx.experiment_.Goniometer(*gon); // the file was right; put it back + } + } + } + return best; +} + +// Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced +// axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let +// a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the +// many high-order reflections off along the fine axis, so the true cell scores worst on the raw +// cloud). Recover the true metric with a second search whose FFT unit-cell bound is widened past +// the default, then RE-INDEX at full resolution constrained by the cell that comes back, as a +// reference (the -C path): the reference filter drops the collapsed/supercell candidates and +// refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing +// crystals never pay for it and their result is untouched. Returns the re-index where it indexes +// more validation frames than the `standing` pass it was asked about. +auto Rugnux::FirstPassRun::LongAxisRescue(const FirstPass &standing) -> std::optional { + if (rx.cancelled_ || !standing.result.has_value() + || standing.score >= 0.5 * static_cast(validation.size())) + return std::nullopt; + auto coarse_settings = rx.experiment_.GetIndexingSettings(); + // The resolution setting only sizes the FFT histogram - it does not filter which spots go + // into it - so this is here to pay for the widened bound below, not to make the long axis + // easier to pick out. Measured, coarsening does the opposite: the acceptance angle widens + // (0.159 -> 0.579 deg) but the true axis falls from rank 1996 to 16361 of 16384 directions. + coarse_settings.FFT_HighResolution_A(3.5f); + // The widened bound is what the rescue turns on, and it costs nothing: the histogram length is + // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already + // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at + // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover: + // a cell whose long axis is past the default 500 A bound is unrepresentable, so the + // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a + // well-indexing crystal never pays. Stills never reach here. + constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f; + if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A) + coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A); + coarse_settings.IndexingThreads(2); + coarse_settings.RefineThreads(std::clamp(rx.config_.nthreads / 2, 1, MAX_REFINE_THREADS)); + IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse); + // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full- + // resolution per-frame validation would be low (the coarse cell is metrically right but + // imprecise), so do NOT score it here - it is only the reference for the constrained re-index. + const SchemeIndexers cris = FeedSchemes(coarse_pool); + RunIndexing(cris); + std::optional coarse_ref; + for (auto &ri : cris) { + auto r = ri->GetLattice(); + if (r.has_value() && (!coarse_ref.has_value() || MaxAxis(*r) > MaxAxis(*coarse_ref))) + coarse_ref = r; + } + // A recovered axis at least 30% longer than the standard pass found signals a rescued metric. + if (!rx.cancelled_ && coarse_ref.has_value() + && MaxAxis(*coarse_ref) > 1.3 * MaxAxis(*standing.result)) { + logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); " + "re-indexing constrained by that cell", MaxAxis(*coarse_ref), MaxAxis(*standing.result)); + rx.experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell()); + // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The + // standard pool was sized at the default bound, so re-using it here would search for a + // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace. + // CountIndexed forces the lattice and never touches the pool, so the existing + // IndexAndRefine is reused; only the FFT needs the wider one. + auto wide_settings = rx.experiment_.GetIndexingSettings(); + wide_settings.IndexingThreads(2); + wide_settings.RefineThreads(std::clamp(rx.config_.nthreads / 2, 1, MAX_REFINE_THREADS)); + RaiseFFTBoundForKnownCell(wide_settings, rx.experiment_); + FirstPass constrained; + if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) { + IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse); + constrained = PickBest(wide_pool, *indexer); + } else { + constrained = PickBest(*indexer_pool, *indexer); + } + rx.experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo + if (constrained.result.has_value() && constrained.score > standing.score) + return constrained; + } + return std::nullopt; +} + +// The background centre, tried as a SECOND FIRST PASS. +// +// The estimate itself is free - it is read off the projection --detect-beam-stop has already +// built - and a first pass re-uses the spots it has already found, so indexing the run a +// second time at the measured centre costs about one more first pass. +// That buys the comparison an on-failure search structurally cannot make: a centre wrong +// ALONG the spindle does not fail. It keeps 96-100 % of frames indexed and quietly returns a +// 2x, 3x or sqrt(3) axis harmonic, so no failure trigger ever fires on it. Indexing both +// centres and comparing the two ANSWERS is what can see it. +// +// What must NOT arbitrate the two is the frame count. Acceptance is a fractional-Miller test, +// so a cell twice as long has to place every spot twice as accurately to score the same: +// measured on real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a +// halved axis at 100.00 %, and the WRONG answer indexed better. A rule of the form "take the +// centre that indexes more" would therefore pick wrong in exactly the case this comparison +// exists for. So only what needs no arbiter is decided here: +// +// * the file's centre does not index a majority -> nothing defends it, so try the measured +// one (and the axis sign with it, see below). A run that indexes nothing has nothing to +// lose; where the measured one is no better, the file's centre is put back; +// * the file's centre indexes and the two agree on the lattice -> keep the file's centre and +// say they agree. This is the overwhelming majority, and saying it is worth the pass: it +// is a free statement that the header is good enough for this crystal; +// * they index and disagree -> say so, with both cells and their primitive-volume ratio. +// Which of an integer-related pair is real is not decided on the frame counts, for the +// reason above; it is decided, in ONE direction only, on the pooled spots - see the +// branch itself. Naming the disagreement is useful even where nothing is decided. +// +// Where the file's centre indexes nothing and the measured one does not either, the axis-sign +// question is asked AGAIN at the measured centre before giving up. The two unknowns are +// discrete and coupled, and the rescue above asked its question at a centre that may be 73 px +// wrong: on the two depositions here whose header is the geometric centre of the detector, +// each error hides the other, both runs sit at 0/60 until both are moved, and moving both +// takes them to 60/60 on the deposited lattice. Asking the pair costs one more first pass on +// a run that has already failed twice. +auto Rugnux::FirstPassRun::BeamCentreCheck(FirstPass best, bool &long_axis_asked) -> FirstPass { + rx.JoinBeamCenterCapture(); + if (rx.cancelled_ || !rx.config_.beam_center_check || !rx.background_center_) + return best; + // Consumed, not re-read: a two-pass rotation run comes back through here with the geometry + // the first pass post-refined, and a pre-scan estimate of the FILE's centre says nothing + // about that one. The pre-scan does not re-measure it either - by then the beam-stop shadow + // is already in the mask, so the projection is not rebuilt. + const BeamCenterEstimate measured = *rx.background_center_; + rx.background_center_.reset(); + const float header_x = rx.experiment_.GetBeamX_pxl(), header_y = rx.experiment_.GetBeamY_pxl(); + const float measured_x = measured.beam_x_pxl; + const float measured_y = measured.beam_y_pxl; + const float moved = std::hypot(measured_x - header_x, measured_y - header_y); + const Centres c{measured, header_x, header_y, measured_x, measured_y, moved}; + // No gate on how significant the move is. It is tempting to skip the second pass when the + // two centres are within the fit's own noise - it would save it on a quarter of the runs + // here - but that is deciding at a threshold before the evidence is in, and the evidence + // is cheap. It would also gate away exactly the case worth asking about: measured on real + // data, a 0.12 px change of centre is the whole difference between the deposited cell and + // a halved axis, and 0.12 px is well inside the sigma of any fit that could gate it. + const int majority = static_cast(validation.size()) / 2; + TryBeamCenter(measured_x, measured_y); + FirstPass alt = IndexAtMeasuredCentre(); + const bool header_indexes = best.result.has_value() && best.score > majority; + const bool measured_indexes = alt.result.has_value() && alt.score > majority; + if (!header_indexes) + return CentreWhenFileFails(std::move(best), std::move(alt), c, long_axis_asked); + RestoreBeamCenter(header_x, header_y); + if (measured_indexes) + return CentreWhenBothIndex(std::move(best), alt, c); + logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes " + "{}/{} validation frames against the file's {}/{} - the file's " + "centre stands", measured_x, measured_y, alt.score, + static_cast(validation.size()), best.score, + static_cast(validation.size())); + return best; +} + +// The first pass at the measured centre: the speculative one started beside the file centre's +// (StartCentreSpeculation) where it read what the run has here, otherwise run here. +auto Rugnux::FirstPassRun::IndexAtMeasuredCentre() -> FirstPass { + FirstPass alt; + if (centre_spec) + centre_spec->done.get(); + if (centre_spec && centre_spec->experiment_key == rx.ExperimentKey() + && centre_spec->key == rx.SpotFindingKey(*spot_mapping) + && centre_spec->spots_per_image == first_pass_spots_per_image) { + // What PickBest would have found and fed here: hand the spots to the run's cache and + // its memo as PrefetchSpots does, and score the schemes on the run's own indexer. + spot_cache = centre_spec->spots; + { + std::lock_guard lock(rx.first_pass_spots_->m); + auto &entries = rx.first_pass_spots_->entries; + if (std::none_of(entries.begin(), entries.end(), + [&](const auto &e) { return e.first == centre_spec->key; })) { + if (entries.size() == MAX_SPOT_KEYS) + entries.erase(entries.begin()); + entries.emplace_back(centre_spec->key, centre_spec->spots); + } + } + alt = ScoreSchemes(centre_spec->ris, *indexer); + } else + alt = PickBest(*indexer_pool, *indexer); + centre_spec.reset(); + return alt; +} + +// The file's centre does not index a majority. +auto Rugnux::FirstPassRun::CentreWhenFileFails(FirstPass best, FirstPass alt, const Centres &c, + bool &long_axis_asked) -> FirstPass { + const float header_x = c.header_x, header_y = c.header_y; + const float measured_x = c.measured_x, measured_y = c.measured_y; + const float moved = c.moved; + const int majority = static_cast(validation.size()) / 2; + bool measured_indexes = alt.result.has_value() && alt.score > majority; + // Nothing defends the file's centre: the data do not index at it, and the + // background says it is `moved` px away. So try the measured one, and if that + // does not index either, ask the AXIS-SIGN question again there. The two unknowns + // are discrete and coupled, and the sign rescue above asked its question at a + // centre that may itself be wrong - each error hides the other, and only moving + // both reveals anything. + bool joint = false; + const auto gonio = rx.experiment_.GetGoniometer(); + if (!measured_indexes && gonio.has_value() && gonio->IsScanning()) { + GoniometerAxis flipped = *gonio; + flipped.Axis(-gonio->GetAxis()); + rx.experiment_.Goniometer(flipped); + const FirstPass both = PickBest(*indexer_pool, *indexer); + if (both.result.has_value() && both.score > majority) { + logger.Warning("Beam centre from the background: the file's centre " + "({:.2f},{:.2f}) indexes {}/{} validation frames and the " + "measured ({:.2f},{:.2f}) {}/{}, but the measured centre " + "WITH the opposite rotation-axis sign indexes {}/{}. The " + "file has both wrong, and neither shows while the other " + "stands; the run takes the measured centre and the opposite " + "sign", header_x, header_y, best.score, + static_cast(validation.size()), measured_x, measured_y, + alt.score, static_cast(validation.size()), both.score, + static_cast(validation.size())); + alt = both; + measured_indexes = true; + joint = true; + } else { + rx.experiment_.Goniometer(*gonio); + } + } + // Third coupled unknown: how DEEP the pass reads. A centre 170 px out and a spot + // list too deep to index at ANY centre hide each other exactly as the centre and + // the axis sign do - the measured centre is right, and it is thrown away because + // the full list indexes nothing there either. Measured on one such crystal: the + // file's centre indexed 1/60 and the measured one 0/60, so the file's was kept and + // the run died; with the ladder's lean seed the measured centre indexes 60/60, and + // the run then indexes every frame of the sweep and merges to the deposited + // resolution. So ask the ladder here too, before giving the centre up. + // + // The depth is asked at BOTH centres, and the centre moves only if it is the + // centre that paid. A lean rung at the measured centre scored against the whole + // spot list at the file's changes two things at once, and the count cannot say + // which of them the gain came from - the same reason a trial centre is scored on + // its own spots at all. Measured on one crystal: the measured centre's ladder + // reached 40/60 on a FIVE-fold axis harmonic while the file's own centre reached + // 49/60 on the deposited cell, and the run moved the centre, shipped the + // supercell, and skipped the beam-centre walk that would have found the crystal's + // own. An exact tie keeps the file's centre: nothing then says the centre is wrong. + bool lean = false; + bool deeper_at_home = false; + if (!measured_indexes) { + const FirstPassState before = committed; + RestoreBeamCenter(header_x, header_y); + const auto home = FirstPassLadder(best.score, 1, majority + 1); + const FirstPassState home_state = committed; + committed = before; + ApplyState(committed); + TryBeamCenter(measured_x, measured_y); + const auto away = FirstPassLadder(alt.score, 1, majority + 1); + const FirstPassState away_state = committed; + // Both centres index a majority, and on the same lattice. Then the two counts + // are two readings of one hypothesis at two geometries, and which is higher is + // the count's noise: measured on one crystal, the same two ladders gave the + // measured centre 56/60 against the file's 50/60, and 48/60 against 50/60 once + // a few thousand defective pixels were masked - and the run that then kept the + // file's centre, 12.6 px out, merged with twice the frames rejected. So the + // count decides nothing here: the file's centre stays for this pass, and the + // measured one is carried forward to be judged on the merges (RunAllPasses). + const bool same_lattice = + home && away + && home->result->search_result.system == away->result->search_result.system + && home->result->search_result.centering + == away->result->search_result.centering + && std::max(home->vol, away->vol) < 1.02 * std::min(home->vol, away->vol); + if (home && (!away || same_lattice || home->score >= away->score)) { + committed = home_state; + RestoreBeamCenter(header_x, header_y); + ApplyState(committed); + best = *home; + deeper_at_home = true; + if (same_lattice && geometry_prepass) { + rx.beam_center_alternative_ = std::array{measured_x, measured_y}; + rx.beam_center_alternative_state_ = away_state; + logger.Warning("Beam centre check: reading less of each frame, the " + "file's centre indexes {}/{} validation frames and the " + "measured centre {:.2f} px away {}/{}, on the same " + "lattice - the counts cannot say which centre is right, " + "so the run carries BOTH centres forward and adopts " + "whichever measures more of the crystal", + home->score, static_cast(validation.size()), moved, + away->score, static_cast(validation.size())); + } else + logger.Warning("Beam centre check: reading less of each frame takes the " + "file's centre to {}/{} validation frames against {}/{} at " + "the measured centre {:.2f} px away - it is how deep this " + "run was reading and not the centre, so the centre stays", + home->score, static_cast(validation.size()), + away ? away->score : 0, + static_cast(validation.size()), moved); + } else if (away) { + alt = *away; + measured_indexes = true; + lean = true; + } + } + if (measured_indexes) { + if (lean) + logger.Warning("Beam centre from the background: neither centre indexes the " + "whole spot list, but at the measured ({:.2f},{:.2f}) the " + "lean first pass above indexes {}/{} against the file's " + "{}/{} - the beam centre in this file is wrong, and only a " + "pass that reads less of each frame could see it", + measured_x, measured_y, alt.score, + static_cast(validation.size()), best.score, + static_cast(validation.size())); + else if (!joint) + logger.Warning("Beam centre from the background: the file's ({:.2f},{:.2f}) " + "indexes {}/{} validation frames and the measured " + "({:.2f},{:.2f}) indexes {}/{} - the beam centre in this " + "file is wrong, and the run adopts the measured one", + header_x, header_y, best.score, + static_cast(validation.size()), measured_x, measured_y, + alt.score, static_cast(validation.size())); + best = alt; + } else if (!deeper_at_home) { + best = CentreByPooledSpots(std::move(best), std::move(alt), c, long_axis_asked); + } + return best; +} + +// Nothing worked. Put the file's centre back, so a run that fails for some +// other reason fails at the geometry it was given rather than at a fit that has +// just been tested and found no better. +// +// Unless the frame counts are both at the floor and the POOLED spots do decide: +// a crystal whose spots are a minority of every frame indexes 0/60 at any centre, +// and then only the pooled excess over chance - each lattice measured against its +// own wrong-spindle null, so a larger cell is not credited for fitting more by +// chance - says which centre sees it. +// +// Each centre is judged on the best lattice it gives, so the long-axis rescue +// is asked at BOTH before they are compared. The lattices standing here index +// no majority anywhere, which is exactly the pass that rescue exists for, and +// the two counts of spots on two lattices that are not the crystal's are close +// to a coin toss: measured on a 545 A axis, the measured centre (3.6 px out +// across the spindle) won it on 16.7% against 15.5%, the rescue then ran +// there alone and could not see the long axis, and the run shipped a P1 +// sub-cell - while at the file's centre the same rescue indexes 60/60. +auto Rugnux::FirstPassRun::CentreByPooledSpots(FirstPass best, FirstPass alt, const Centres &c, + bool &long_axis_asked) -> FirstPass { + const float header_x = c.header_x, header_y = c.header_y; + const float measured_x = c.measured_x, measured_y = c.measured_y; + const float moved = c.moved; + if (auto rescued = LongAxisRescue(alt)) + alt = *rescued; + PooledEvidence at_measured, at_header; + if (alt.result.has_value()) + at_measured = MeasurePooled(*indexer, *alt.result); + RestoreBeamCenter(header_x, header_y); + if (auto rescued = LongAxisRescue(best)) + best = *rescued; + if (best.result.has_value()) + at_header = MeasurePooled(*indexer, *best.result); + long_axis_asked = true; + if (alt.result.has_value() && ValidationEvidenceBeatsChance(at_measured) + && ValidationEvidencePrefers(at_header, at_measured)) { + TryBeamCenter(measured_x, measured_y); + logger.Warning("Beam centre from the background: neither centre indexes a " + "majority of the validation frames ({}/{} at the measured, " + "{}/{} at the file's), but the pooled spots put " + "{:.1f}% on the lattice at the measured ({:.2f},{:.2f}) against " + "{:.1f}% at the file's (chance {:.1f}% and {:.1f}%) - the run " + "adopts the measured centre", + alt.score, static_cast(validation.size()), + best.score, static_cast(validation.size()), + 100.0 * at_measured.on_lattice / std::max(1, at_measured.spots), + measured_x, measured_y, + 100.0 * at_header.on_lattice / std::max(1, at_header.spots), + 100.0 * at_measured.by_chance / std::max(1, at_measured.spots), + 100.0 * at_header.by_chance / std::max(1, at_header.spots)); + best = alt; + } else { + logger.Info("Beam centre check: the file's centre ({}/{} frames) and the " + "measured one {:.2f} px away ({}/{}) - the measured centre puts no " + "more of the pooled spots on its lattice, so the file's centre is " + "kept", + best.score, static_cast(validation.size()), moved, alt.score, + static_cast(validation.size())); + } + return best; +} + +// The file's centre indexes a majority, and so does the measured one. +auto Rugnux::FirstPassRun::CentreWhenBothIndex(FirstPass best, const FirstPass &alt, const Centres &c) -> FirstPass { + const BeamCenterEstimate &measured = c.measured; + const float measured_x = c.measured_x, measured_y = c.measured_y; + const float moved = c.moved; + // Primitive volumes, as everywhere else here: two passes can find the same + // lattice in different settings, and a centred setting's cell is an exact + // integer multiple of its primitive one. + const double ratio = alt.vol < best.vol ? best.vol / alt.vol : alt.vol / best.vol; + const bool same_class = + best.result->search_result.system == alt.result->search_result.system + && best.result->search_result.centering == alt.result->search_result.centering; + // The primitive volume is what the beam centre decides. A centre error along + // the spindle makes the FFT take an axis harmonic, which changes the volume by + // an integer factor or by sqrt(3); the Bravais class can differ for a reason + // that has nothing to do with the centre, so the two are reported separately. + if (ratio < 1.02 && same_class) { + rx.Note(fmt::format("beam centre kept: measured {:.1f} px off, same lattice", moved)); + logger.Info("Beam centre check: indexing at the measured centre " + "({:.2f},{:.2f}) returns the same lattice as the file's - " + "{}-centred {}, cells agreeing to {:.2f} % in volume. The {:.2f} px " + "difference does not decide this crystal's cell", + measured_x, measured_y, alt.result->search_result.centering, + gemmi::crystal_system_str(alt.result->search_result.system), + 100.0 * (ratio - 1.0), moved); + // Unless the two are further apart than this geometry absorbs. A centre + // wrong by about one reflection spacing indexes the same cell just as + // well, with every reflection one index out over part of the sweep, and + // only the merge shows it - so the measured centre is kept as a hypothesis + // for RunAllPasses, which tries it if the file's centre merges + // inconsistently. + if (geometry_prepass && moved > 3.0f * measured.sigma_pxl + && moved > BeamCenterNeed_pxl(rx.experiment_)) + rx.beam_center_same_lattice_alternative_ = + std::array{measured_x, measured_y}; + } else if (ratio < 1.02) { + // Same cell, different metric symmetry read off it. The cell is not in + // question, but the CLASS is not cosmetic either: GeometryRefiner holds + // the cell to it for the rest of the run, so a class that is one axis + // equality away from the truth is refined, predicted, integrated and + // merged with that equality imposed. + // + // Whether the centre decided it is arithmetic. The equality is held to + // LATTICE_SEARCH_DIST_TOLERANCE, relative, on a pair of reduced-cell + // axes; a centre wrong by `moved` px moves an axis of length a by + // moved*p*a/(D*lambda) relative, the FFT reading that axis off spot + // positions the centre displaces. Both axes of the pair carry their own + // error, so the centre can open or close a gap of their sum - and where + // that reaches the tolerance the higher class is not a fact about the + // crystal but a consequence of where the beam was assumed to be. + // + // Nothing available here can settle it: both centres index, and the frame + // count is the one statistic this check already refuses to arbitrate on. + // So carry the other centre out of the pass and let the two arms be judged + // on their merges (RunAllPasses). + const auto reduced = best.result->search_result.primitive_reduced.GetUnitCell(); + const double axis_pair = std::max(reduced.a + reduced.b, + reduced.b + reduced.c); + const double induced = moved * rx.experiment_.GetPixelSize_mm() * axis_pair + / (rx.experiment_.GetDetectorDistance_mm() + * rx.experiment_.GetWavelength_A()); + const bool centre_can_decide = induced >= LATTICE_SEARCH_DIST_TOLERANCE; + logger.Info("Beam centre check: indexing at the measured centre " + "({:.2f},{:.2f}) returns the same CELL as the file's - the " + "primitive volumes agree to {:.2f} % - but a different metric " + "symmetry: {}-centred {} from the file's centre against " + "{}-centred {} from the measured one", measured_x, measured_y, + 100.0 * (ratio - 1.0), best.result->search_result.centering, + gemmi::crystal_system_str(best.result->search_result.system), + alt.result->search_result.centering, + gemmi::crystal_system_str(alt.result->search_result.system)); + if (centre_can_decide && geometry_prepass) { + rx.beam_center_alternative_ = std::array{measured_x, measured_y}; + logger.Warning("Beam centre check: {:.2f} px is worth {:.2f} % on the " + "axis pair the promotion compares, against the {:.2f} % " + "the metric symmetry is decided to - the difference is " + "large enough to have decided it either way, so the run " + "carries BOTH centres forward and adopts whichever " + "merges better", moved, 100.0 * induced, + 100.0 * LATTICE_SEARCH_DIST_TOLERANCE); + } else { + logger.Info("Beam centre check: {:.2f} px is worth {:.2f} % on the axis " + "pair the promotion compares, against the {:.2f} % the " + "metric symmetry is decided to - the run keeps the file's " + "centre and the symmetry it found", moved, 100.0 * induced, + 100.0 * LATTICE_SEARCH_DIST_TOLERANCE); + } + } else { + best = CentresDisagree(std::move(best), alt, c, ratio); + } + return best; +} + +// The two centres index different lattices. +auto Rugnux::FirstPassRun::CentresDisagree(FirstPass best, const FirstPass &alt, const Centres &c, + double ratio) -> FirstPass { + const float header_x = c.header_x, header_y = c.header_y; + const float measured_x = c.measured_x, measured_y = c.measured_y; + const auto &hc = best.result->search_result.conventional.GetUnitCell(); + const auto &mc = alt.result->search_result.conventional.GetUnitCell(); + const double nearest = std::round(ratio); + const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15) + || std::abs(ratio - std::sqrt(3.0)) < 0.15; + logger.Warning("Beam centre check: the two centres DISAGREE about the " + "lattice. The file's ({:.2f},{:.2f}) gives {}-centred {}, " + "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} (primitive volume " + "{:.0f} A^3, {}/{} frames); the measured ({:.2f},{:.2f}) gives " + "{}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} " + "(primitive volume {:.0f} A^3, {}/{} frames){}", + header_x, header_y, best.result->search_result.centering, + gemmi::crystal_system_str(best.result->search_result.system), + hc.a, hc.b, hc.c, hc.alpha, hc.beta, hc.gamma, best.vol, + best.score, static_cast(validation.size()), + measured_x, measured_y, alt.result->search_result.centering, + gemmi::crystal_system_str(alt.result->search_result.system), + mc.a, mc.b, mc.c, mc.alpha, mc.beta, mc.gamma, alt.vol, + alt.score, static_cast(validation.size()), + harmonic ? fmt::format(" - the volumes differ by {:.2f}x, so one " + "is an axis harmonic of the other and the " + "beam centre is what settles which", + ratio) + : ""); + // Deliberately no decision on the frame counts. They are the one statistic + // that must not choose: on the measured case the harmonic indexed MORE. + // + // Decided in either direction on the pooled excess over chance, the same + // measurement the failed-pass branch above adopts a centre on. Each lattice is + // scored against its own wrong-spindle null, so the larger cell is not credited + // for the spots a denser lattice catches by accident, and the smaller one is + // not credited for the spots it misses. + // + // Both directions need it. Where the measured centre holds the LARGER cell, + // nothing else in the run repairs it - the second pass refuses a lattice larger + // than pass 1's (the supercell guard below). Where it holds the SMALLER one, + // the second pass can drop to it only if its own beam refinement happens to + // walk the whole centre error, and that rescue is chaotic: measured on one + // crystal, a 0.1 % change in the pass-1 distance switched off a 6.5 px walk and + // the run kept the doubled axis. + const double nearest_n = std::round(ratio); + const bool integer_harmonic = nearest_n >= 2.0 && nearest_n <= 4.0 + && std::abs(ratio - nearest_n) < 0.15; + bool adopted_measured = false; + if (integer_harmonic) { + const PooledEvidence at_header = MeasurePooled(*indexer, *best.result); + TryBeamCenter(measured_x, measured_y); + const PooledEvidence at_measured = MeasurePooled(*indexer, *alt.result); + adopted_measured = ValidationEvidenceBeatsChance(at_measured) + && ValidationEvidencePrefers(at_header, at_measured); + logger.Info("Beam centre check: pooled validation spots on the lattice " + "{:.1f}% at the file's centre and {:.1f}% at the measured " + "(chance {:.1f}% and {:.1f}%)", + 100.0 * at_header.on_lattice / std::max(1, at_header.spots), + 100.0 * at_measured.on_lattice + / std::max(1, at_measured.spots), + 100.0 * at_header.by_chance / std::max(1, at_header.spots), + 100.0 * at_measured.by_chance + / std::max(1, at_measured.spots)); + if (adopted_measured) { + logger.Warning("Beam centre check: the spots are on the {} cell of the " + "measured centre, not on the file's - so the {:.2f}x is " + "the file's centre error and not the crystal, and the " + "run adopts the measured centre", + alt.vol > best.vol ? "larger" : "smaller", ratio); + best = alt; + } else { + RestoreBeamCenter(header_x, header_y); + } + } + if (!adopted_measured) + logger.Warning("Beam centre check: the run keeps the file's centre. The " + "frame counts cannot arbitrate this - a cell twice as long " + "must place every spot twice as accurately to score the " + "same, so the harmonic is often the one that indexes more. " + "Where the beam centre is known, give it with " + "--beam-x/--beam-y"); + return best; +} + +// Where the first pass indexes no majority of the validation frames, the end wedge is asked too +// and judged against the standing pass as any later scheme is. Only then: on a run the schemes +// index it would be one more lattice for the axis-harmonic arbiter, which is not decided on +// frame counts. Asked once, after the axis-sign, beam-centre and long-axis rescues (below), not +// inside every pass they run: a consecutive wedge indexes part of the validation frames at the +// WRONG axis sign too (measured: 33/60 where the spread at the right sign indexes 60/60), and +// counted in the pass a rescue is gated on, it lifts that pass past the gate and the rescue +// that would have indexed every frame is never asked. +auto Rugnux::FirstPassRun::AddEndWedge(FirstPass bp, IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass { + if (rx.cancelled_ || bp.score > static_cast(validation.size()) / 2) + return bp; + SchemeIndexers ris; + ris.push_back(FeedScheme(pool, end_wedge[0].second)); + RunIndexing({ris[0].get()}); + return ScoreSchemes(ris, idx, &end_wedge, std::move(bp)); +} + +// Short-axis hypothesis - the small-molecule case. +// +// LAST, after the axis-sign, beam-centre and long-axis rescues, deliberately. Those all +// run the standard floor and compare their result against `best` on volume ratios; a cell +// adopted here is a k-fold sub-cell, so running before them makes every one of those +// comparisons read an integer ratio and report an axis harmonic that is not one (measured: +// the beam-centre check announced a "3.00x harmonic ... the beam centre is what settles +// which" on a pair that differs only by the FLOOR). Here the comparison is like with like, +// and this pass sees whatever geometry those rescues settled on. +// +// fft_min_unit_cell_A (10 A) is a protein constant: shorter than any protein axis and longer +// than a small molecule's. It is applied to the FFT's PER-DIRECTION PEAK SEARCH +// (FFTIndexerCPU.cpp / FFTIndexerGPU.cu skip every histogram bin below it) and only then to +// the assembled cell (FFTIndexer::ReduceAndRefine, PostIndexingRefinement). The first of +// those is what does the damage: the transform reports ONE peak per direction, so a +// lattice row whose real-space repeat is shorter than the floor is reported at its first +// HARMONIC above it, and the cell built out of those harmonics is an exact integer +// sublattice of the true one - a real-space supercell. Measured on chemical-crystallography +// data: index 2, 3 and 6 exactly, and the one dataset whose true reduced cell clears 10 A is +// the one that comes out right. Because the floor is applied to the REDUCED cell it rejects +// the lattice and not merely a basis choice, so nothing downstream can re-describe its way +// back. +// +// Lowering the floor by default is not available: it only ever ADDS candidates, and the +// Bravais classification the candidate pool feeds is under-determined, so it moves crystals +// whose cell it never touched (measured: 2 regressions and 1 improvement over 20 protein +// datasets, all three in the metric symmetry rather than the cell). So the low floor runs as +// a SECOND hypothesis and the standard pass's answer stands unless the evidence takes it. +// +// What the evidence can be, before integration: a sublattice's per-spot residuals are +// dominated pointwise by the finer lattice it contains, so every margin statistic - squares, +// counts, trimmed sums - is biased toward the LARGER cell, and a spurious supercell indexes +// every frame its true cell does. Only the SATURATION is informative. Hence the same shape +// the scheme ratchet in PickBest already uses: an integer primitive-volume ratio plus a tie +// on validation frames, or an outright win, which runs against the bias and so is real. The +// statistic that genuinely decides an axis multiple is the occupancy of the reflection class +// the multiple adds - AbsenceEvidence in SearchSpaceGroup - and that needs integrated +// intensities, so it cannot run here. +// Not when a cell was given. RaiseFFTBoundForKnownCell has already lowered the floor for +// the standard pass against that cell, so this pass has nothing left to add - and it takes +// its settings from rx.experiment_, which never carried the RAISED max bound, so on a cell +// longer than the default search it could only ever return a harmonic of what was asked +// for, tie on frames the way a sub-multiple always does, and be adopted over it. +auto Rugnux::FirstPassRun::ShortAxisHypothesis(FirstPass best) -> FirstPass { + if (!rx.cancelled_ && !rx.experiment_.GetUnitCell().has_value()) { + std::unique_ptr short_pool; + std::vector> sris; + if (short_spec && short_spec->key == rx.ExperimentKey() && short_spec->mapping == spot_mapping + && short_spec->spots_per_image == first_pass_spots_per_image) { + for (auto &f : short_spec->futures) + f.get(); + short_pool = std::move(short_spec->pool); + sris = std::move(short_spec->sris); + } else { + short_spec.reset(); + short_pool = std::make_unique(ShortAxisSettings(rx.experiment_), + IndexerConstruction::OnFirstUse); + sris = FeedSchemes(*short_pool); + RunIndexing(sris); + } + + for (size_t i = 0; i < sris.size() && !rx.cancelled_; i++) { + auto found = sris[i]->GetLattice(); + if (!found.has_value()) + continue; + RotationIndexerResult alt = std::move(*found); + const std::string &name = schemes[i].first; + // Primitive volumes, for the reason PickBest gives: two searches can find the same + // lattice in different settings, and a centred setting's cell is an exact integer + // multiple of its primitive one. + const double vol = std::abs( + alt.lattice.ToPrimitive(alt.search_result.centering).CalcVolume()); + if (!(vol > 1.0) || !(best.vol > 1.0)) + continue; + const double ratio = best.vol / vol; + const double nearest = std::round(ratio); + // The gate is here rather than after the scoring because scoring is the expensive + // half - CountIndexed runs the real per-image path over the validation set, seconds + // on a protein - while the volume ratio is free. On every protein measured the low + // floor returns the SAME lattice (ratio 1.00-1.01), so this costs them one extra FFT + // and nothing else. + if (!(vol < best.vol) || nearest < 2.0 || std::abs(ratio - nearest) >= 0.15) { + logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': primitive volume " + "{:.0f} A^3 against the standing cell's {:.0f} A^3, ratio {:.2f} - " + "not an integer sub-cell, keeping the standing cell", + IndexingSettings::fft_min_unit_cell_limit_A, name, vol, best.vol, + ratio); + continue; + } + + // The premise of this pass, checked rather than asserted. The pass exists because a + // lattice row shorter than the standard floor is reported at its first HARMONIC + // above it, so a cell it may rightfully adopt has a reduced axis BELOW that floor - + // and the floor is applied to the reduced cell, which is the one to measure. Where + // the alternative has no such axis the floor cannot be what the standard pass + // tripped over, and an integer volume ratio is then an ordinary axis multiple: the + // question the axis-harmonic arbiter settles just above on the occupancy of the + // class the multiple adds. That is intensity-class evidence about the very + // reflections in dispute, where the tie below is a frame count this comment already + // says a sublattice cannot lose - so where the two disagree the arbiter is right and + // this block must not overrule it. Measured on the corpus: the four small-molecule + // sets the pass is for all have a reduced axis under 9 A, and the two protein sets + // it halved have 35 and 47 A. + const auto alt_reduced = alt.lattice.ToPrimitive(alt.search_result.centering) + .NiggliReduce().GetUnitCell(); + const float alt_shortest = std::min({alt_reduced.a, alt_reduced.b, alt_reduced.c}); + const float standard_floor = rx.experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A(); + if (alt_shortest >= standard_floor) { + logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': a {:.0f}x sub-cell " + "(vol {:.0f} A^3), but its shortest reduced axis is {:.1f} A, above " + "the {:.1f} A floor the standard pass already searched - the floor " + "is not why that cell was missed, so this is an ordinary axis " + "multiple and not this pass's to decide; keeping the standing cell", + IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol, + alt_shortest, standard_floor); + continue; + } + + const int alt_score = CountIndexed(*indexer, alt); + // A pass that indexes nothing is not tied with a pass that indexes nothing. Without + // this, 0 against 0 passes the tie band and the volume ratio alone decides, which is + // the one thing this block must never do. + if (alt_score <= 0 || best.score <= 0) { + logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {}/{} validation " + "frames against {} - nothing indexes, keeping the standing cell", + IndexingSettings::fft_min_unit_cell_limit_A, name, alt_score, + static_cast(validation.size()), best.score); + continue; + } + const bool tied = static_cast(alt_score) >= best.score * 0.9f - 0.5f; + const bool clearly_more = static_cast(alt_score) > best.score * 1.1f + 0.5f; + if (!tied && !clearly_more) { + logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {:.0f}x sub-cell " + "(vol {:.0f} A^3) indexes {}/{} validation frames against {} - " + "keeping the standing cell", + IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol, + alt_score, static_cast(validation.size()), best.score); + continue; + } + const auto uc = alt.lattice.GetUnitCell(); + logger.Warning("Short-axis pass (floor {:.0f} A), scheme '{}': the standing " + "cell is a {:.0f}x supercell. Its {:.0f} A^3 primitive volume against " + "this pass's {:.0f} A^3, at {}/{} validation frames against {} - " + "adopting {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}. The shortest axis " + "is below the {:.0f} A default FFT floor, which is why the standard " + "pass could only build a multiple of it", + IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, best.vol, + vol, alt_score, static_cast(validation.size()), best.score, + uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, + rx.experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A()); + best.score = alt_score; + best.vol = vol; + best.name = name + " (short-axis)"; + best.result = std::move(alt); + } + } + return best; +} + +// Sub-lattice hypothesis: the standing cell's seven index-2 sub-lattices, in its own +// orientation and geometry, each scored on the validation frames. +// +// The FFT can only offer the cells its peaks build, and two lattices seen together - a twin +// domain or a second crystal whose reciprocal lattice shares an index-2 sub-lattice with +// this one's - build the lattice that contains them both: a cell twice the size of either, +// which indexes both crystals' spots and so wins every comparison of candidates the FFT +// returns. Where the true cell is then never offered, nothing downstream can choose it. +// Offering it costs seven frame counts. +// +// The frame count only nominates: a coarser lattice asks less accuracy of every spot, so a +// sub-lattice of the true cell can index as many frames as the cell itself, or a few more +// where the cell's own score is depressed. What decides is the comparison this pass already +// makes between a cell and its axis harmonic: the share of the validation frames' pooled +// spots each puts on itself over and above what a wrong spindle angle puts there, each +// against its own null (ValidationEvidencePrefers). A real superstructure keeps its cell on +// that count, because the class it adds holds reflections; a cell built out of two crystals +// loses it, because it has to place every spot of the one it describes twice as accurately +// with an orientation fitted to both. The nominee must also index a majority of the frames, +// the bar the false-promotion test above asks of a cell that "works", and no fewer frames +// than the standing cell: a standing cell that indexes almost nothing is not this crystal's +// lattice, its sub-lattices are not hypotheses about it, and the rescues below are what +// such a run needs. +auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { + if (!rx.cancelled_ && best.result.has_value() && best.score > 0) { + const CrystalLattice prim = best.result->lattice + .ToPrimitive(best.result->search_result.centering); + const Coord A = prim.Vec0(), B = prim.Vec1(), C = prim.Vec2(); + std::optional sub_best; + std::string counts; + for (int u = 0; u < 2; u++) + for (int v = 0; v < 2; v++) + for (int w = 0; w < 2; w++) { + if (u + v + w == 0 || rx.cancelled_) + continue; + const Coord t = (A * static_cast(u) + B * static_cast(v) + + C * static_cast(w)) / 2.0f; + const CrystalLattice sub = (w ? CrystalLattice(A, B, t) + : v ? CrystalLattice(A, t, C) + : CrystalLattice(t, B, C)).NiggliReduce(); + RotationIndexerResult cand = *best.result; + cand.lattice = sub; + cand.search_result = LatticeSearchResult{}; + cand.search_result.primitive_reduced = sub; + cand.search_result.conventional = sub; + cand.unconstrained.reset(); + cand.extra_lattices.clear(); + const int score = CountIndexed(*indexer, cand); + counts += fmt::format(" ({},{},{}) {}", u, v, w, score); + if (!sub_best || score > sub_best->score) + sub_best = FirstPass{std::move(cand), score, + fmt::format("{} (sub-lattice {}{}{})", best.name, u, v, w), + std::abs(sub.CalcVolume())}; + } + bool adopt = false; + std::string verdict = " - none indexes a majority of the frames and as many as the standing cell"; + if (sub_best && sub_best->score >= best.score + && sub_best->score > static_cast(validation.size()) / 2) { + const PooledEvidence at_standing = MeasurePooled(*indexer, *best.result); + const PooledEvidence at_sub = MeasurePooled(*indexer, *sub_best->result); + adopt = ValidationEvidencePrefers(at_standing, at_sub); + verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance) against " + "the standing cell's {}/{} ({}), {}", + at_sub.on_lattice, at_sub.spots, at_sub.by_chance, + at_standing.on_lattice, at_standing.spots, at_standing.by_chance, + adopt ? "adopting it" : "keeping the standing cell"); + } + logger.Info("Sub-lattice hypothesis: the standing cell ({:.0f} A^3) indexes {}/{} validation " + "frames; its index-2 sub-lattices (half-translation u,v,w):{}{}", + best.vol, best.score, static_cast(validation.size()), counts, verdict); + if (adopt) + best = std::move(*sub_best); + } + return best; +} + +// A lattice that indexes almost none of the validation frames is not this crystal's lattice. +// Forcing it anyway integrates every image against it; the few partials that survive yield no +// fulls in the 3D combine, and the run then dies inside the merge complaining about resolution - +// far from the cause. Stop here instead, while the reason is still visible. The floor is one +// frame in six, several times below the weakest crystal that still merges. It is checked AFTER +// the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue +// score. +// The second (refined-geometry) pass re-indexes DE NOVO, and the post-refined geometry - fitted +// to the lattice pass 1 already found - can tip the blind FFT, or the long-axis rescue above, +// onto another lattice: an axis harmonic of a long cell (measured: a tripled ~100 A c-axis at +// 7/60), or a related centred cell that indexes fewer frames than pass 1's lattice does at +// this very geometry (measured: a C-centred 303 A cell at 41/60 against pass 1's triclinic +// cell at 60/60). Pass 1's lattice is therefore carried here as a hypothesis of its own and +// the two are judged on the same validation frames: where the re-index found a different +// lattice and indexes fewer frames than pass 1's, or too few to be this crystal's lattice at +// all, pass 1's whole result is integrated instead. The same lattice found again is kept as +// the re-index refined it at this geometry. This is the floor-side companion of the +// supercell-collapse guard in RunAllPasses, which only sees a pass-2 lattice that is LARGER. +auto Rugnux::FirstPassRun::Pass1Fallback(FirstPass best) -> FirstPass { + RotationIndexerResult pass1_lattice = *rx.prepass_result_; + // Pass 1's result carries pass 1's goniometer; forcing it whole would put the + // uncorrected angles back (the same repair as the supercell re-run in RunAllPasses). + if (rx.prepass_rotation_scale_ && pass1_lattice.axis) + pass1_lattice.axis = ScaleRotation(*pass1_lattice.axis, *rx.prepass_rotation_scale_); + // Pass 1's lattice was FITTED at pass 1's detector distance, and this pass + // integrates at the post-refined one. Re-scoring it here is not the same as + // re-fitting it: a real-space cell is measured against the distance the spots were + // seen at, so carrying it across a distance change scales the whole cell by the + // ratio of the two. Measured on a 545 A axis over 225 -> 226 mm: 0.4 %, exactly + // 225/226, shipped as the run's answer. Scale it to the distance it will be used + // at; the orientation is unaffected. + const float d_fit = rx.prepass_result_->geom.GetDetectorDistance_mm(); + const float d_use = rx.experiment_.GetDetectorDistance_mm(); + const bool rescaled = d_fit > 0 && d_use > 0 && d_fit != d_use; + if (rescaled) { + const float f = d_use / d_fit; + const auto &l = pass1_lattice.lattice; + pass1_lattice.lattice = CrystalLattice(l.Vec0() * f, l.Vec1() * f, l.Vec2() * f); + } + const CrystalLattice pass1_primitive = pass1_lattice.lattice + .ToPrimitive(pass1_lattice.search_result.centering).NiggliReduce(); + const double pass1_vol = std::abs(pass1_primitive.CalcVolume()); + // The same lattice whatever setting or class it is held in - a centred cell and the + // triclinic primitive of the same points are one lattice, and the re-index's symmetric + // setting must not lose to its own unconstrained primitive on a frame count: the Niggli- + // reduced primitive edges agree within 2 % (the battery's lattice-identity test). + const auto reduced_edges = [](const CrystalLattice &l) { + const auto uc = l.GetUnitCell(); + std::array e = {uc.a, uc.b, uc.c}; + std::sort(e.begin(), e.end()); + return e; + }; + const auto e1 = reduced_edges(pass1_primitive); + const auto e2 = reduced_edges(best.result->lattice + .ToPrimitive(best.result->search_result.centering).NiggliReduce()); + bool same_lattice = true; + for (int i = 0; i < 3; i++) + same_lattice = same_lattice && std::abs(e2[i] - e1[i]) <= 0.02 * e1[i]; + const int floor = static_cast(validation.size()) / 6; + const int pass1_score = (best.score < floor || !same_lattice) + ? CountIndexed(*indexer, pass1_lattice) : -1; + if (best.score < floor || pass1_score > best.score) { + const auto &pc = pass1_lattice.search_result.conventional.GetUnitCell(); + const auto &bc = best.result->search_result.conventional.GetUnitCell(); + logger.Warning("Two-pass: the refined-geometry re-index indexes {}/{} validation frames " + "({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}), pass-1's lattice " + "{}/{} at this geometry ({}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} " + "{:.2f}) - integrating with pass-1's lattice", + best.score, static_cast(validation.size()), + best.result->search_result.centering, + gemmi::crystal_system_str(best.result->search_result.system), + bc.a, bc.b, bc.c, bc.alpha, bc.beta, bc.gamma, + pass1_score, static_cast(validation.size()), + pass1_lattice.search_result.centering, + gemmi::crystal_system_str(pass1_lattice.search_result.system), + pc.a, pc.b, pc.c, pc.alpha, pc.beta, pc.gamma); + if (rescaled) + logger.Info("Two-pass: pass-1's cell was fitted at {:.3f} mm and is used at " + "{:.3f} mm - scaled by {:.6f} to the distance it is integrated at", + d_fit, d_use, d_use / d_fit); + best.result = std::move(pass1_lattice); + best.score = pass1_score; + best.vol = pass1_vol; + best.name = "pass-1 lattice (refined-geometry re-index found less)"; + } + return best; +}