diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d9523860e..105b8138f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -388,19 +388,82 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // first pass being invisible in it. Time the whole thing here, where every pass is inside. const auto run_start = std::chrono::steady_clock::now(); DiffractionExperiment header_experiment = experiment_; - ProcessResult result = RunAllPasses(observer); + // The run as it starts, to be made again where a near tie decided early ends elsewhere on the rest + // of the sweep (decide_late_). + const auto again = std::make_shared(*this); + again->near_tie_rechecks_ = std::make_shared(); + again->decide_late_ = true; + Rugnux *run = this; + ProcessResult result; + bool lost = false; + try { + result = RunAllPasses(observer); + } catch (const NoLatticeFound &) { + // A run that found no lattice after deciding a near tie early is made again deciding late. + if (!near_tie_rechecks_->carried) + throw; + lost = true; + } + if (lost || (!result.cancelled && NearTieOverturned())) { + JoinBeamCenterCapture(); + Logger logger("Rugnux"); + logger.Info("A near tie decided early ends elsewhere on the rest of the sweep - the run is made again, " + "deciding near ties on all the frames"); + result = again->RunAllPasses(observer); + run = again.get(); + } result.total_time_s = std::chrono::duration(std::chrono::steady_clock::now() - run_start).count(); // For the report's beam-centre warning: what the background measured, and what this crystal asks // of the centre - at the file's geometry, as the pre-scan's check asks it, but on the cell the run // determined rather than the long axis that check has to assume. - JoinBeamCenterCapture(); - result.background_beam_centre = measured_beam_center_; + run->JoinBeamCenterCapture(); + result.background_beam_centre = run->measured_beam_center_; if (result.consensus_cell) header_experiment.SetUnitCell(*result.consensus_cell); result.beam_centre_need_pxl = BeamCenterNeed_pxl(header_experiment); return result; } +void Rugnux::StartNearTieRecheck(std::shared_ptr run, bool geometry_prepass) { + // The pass again, up to its first pass, as the probe it would be: the same inputs, no memo to + // answer from, nothing started beside it. + run->decide_late_ = true; + run->near_tie_recheck_ = true; + run->indexing_probe_only_ = true; + run->first_pass_memo_.clear(); + run->first_pass_outcome_.clear(); + run->speculative_probe_.reset(); + run->speculate_geometry_probe_ = false; + run->rotation_scale_check_.reset(); + std::lock_guard lock(near_tie_rechecks_->m); + near_tie_rechecks_->pending.push_back({std::move(run), geometry_prepass, first_pass_outcome_}); +} + +bool Rugnux::NearTieOverturned() { + // After the run: by then every frame is in, live too, and nothing else holds the card. + Logger logger("Rugnux"); + TimingMark mark(logger, "near-tie re-checks"); + std::lock_guard lock(near_tie_rechecks_->m); + int elsewhere = 0; + for (auto &[run, geometry_prepass, outcome] : near_tie_rechecks_->pending) { + try { + run->RunPipeline(nullptr, /*write_output=*/false, geometry_prepass); + } catch (const NoLatticeFound &) { + run->first_pass_outcome_.clear(); // the pass found one; on all the frames it finds none + } + if (run->first_pass_outcome_ != outcome) { + logger.Info("Near-tie re-check: the first pass ends elsewhere on all the frames - {} against {}", + run->first_pass_outcome_, outcome); + elsewhere++; + } + } + if (!near_tie_rechecks_->pending.empty()) + logger.Info("Near ties re-checked on the rest of the sweep: {} first pass(es), {} ending elsewhere", + near_tie_rechecks_->pending.size(), elsewhere); + near_tie_rechecks_->pending.clear(); + return elsewhere > 0; +} + void Rugnux::AddFirstPassMemo(const FirstPassMemo &memo) { for (const auto &m : first_pass_memo_) if (m.key == memo.key) @@ -558,6 +621,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ASKED for. A pass that determined one leaves it on experiment_, and the pass after it - the // re-runs in Run(), or a second Run() on the same object - would otherwise inherit that answer and // skip the search that is the whole point of running again. + // The run as this pass finds it, for the re-check of a near tie its first pass may decide early. + std::shared_ptr recheck; + if (!decide_late_ && config_.mode == ProcessMode::FullAnalysis && config_.rotation_indexing + && config_.two_pass_rotation && !force_rotation_result_.has_value() + && !config_.forced_rotation_lattice.has_value() && experiment_.IsRotationIndexing()) + recheck = std::make_shared(*this); + experiment_.SetSpaceGroup(user_fixed_sg_); const auto dataset = reader_.GetDataset(); @@ -866,7 +936,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Rotation indexer lattice forced externally - skipping first pass"); } else if (full && config_.rotation_indexing && config_.two_pass_rotation) { TimingMark mark(logger, "first-pass indexing"); - if (FirstPassRotationIndexing(p)) + const bool probe_done = FirstPassRotationIndexing(p); + if (recheck && near_tie_carried_) + StartNearTieRecheck(std::move(recheck), geometry_prepass); + if (probe_done) return result; } // The rotation-scale walk beside this pass starts once its frames are indexed (see RunAllPasses). diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0bb4d5da6..24b901afc 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -938,6 +938,44 @@ class Rugnux { // lattice on the validation frames (ProcessResult::validation_evidence), before a single image is // integrated. The rotation-scale walk (see RunAllPasses) compares angle models on that score. bool indexing_probe_only_ = false; + + // --- Near ties: decided early, re-checked late (see FirstPassRun::DecideLate) --- + // + // A first pass acts on its validation frames at once. Where it decided a near tie that the rest + // of the sweep could re-check, the copy of the run taken as its pass began is kept, and once the + // run is over - every frame is in by then, live too - that pass is run again on it as a probe, + // deciding late. Only where a re-check's first pass ends elsewhere than the pass's did is the + // run made again, every first pass deciding late. Offline the answer is then the same as deciding + // late throughout; live, no near tie holds the run up while the collection goes on. + // + // Whether the first passes decide their near ties on all the frames: on a re-check, and on the run + // made again. + bool decide_late_ = false; + // Whether this copy of the run is a re-check (for the log). + bool near_tie_recheck_ = false; + // Set by the first pass now running where it decided a near tie the later frames can re-check. + bool near_tie_carried_ = false; + // What the last first pass ended on - the lattice, the path to it, the centre and the spot-finding + // settings it stands on, the centres it carries forward - for its re-check to compare with. + std::string first_pass_outcome_; + // The re-checks of the whole run, shared with the copies the run makes of itself: a probe on a copy + // leaves its re-check here too. + struct NearTieRechecks { + struct Pending { + std::shared_ptr run; // the run as the pass found it + bool geometry_prepass; + std::string outcome; // what its first pass ended on + }; + std::mutex m; + std::atomic carried = false; // a first pass decided a near tie early + std::vector pending; + }; + std::shared_ptr near_tie_rechecks_ = std::make_shared(); + // Leaves the first pass the pass just run made to be re-checked; `run` is the copy of this run taken + // as the pass began. + void StartNearTieRecheck(std::shared_ptr run, bool geometry_prepass); + // Runs the re-checks. Whether one ended elsewhere than its pass. + bool NearTieOverturned(); // The validation evidence of the last canonical-pass or probe first pass, keyed by the inputs that // pass read (FirstPassInputKey). The first pass is deterministic in them, so an indexing probe // whose key matches takes the stored evidence instead of indexing again - the geometry walk's diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 1b75b11c9..7f8281b00 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -225,6 +225,14 @@ namespace { s.phi = rotation ? rotation->GetAngle_deg(s.image) + rotation->GetWedge_deg() / 2.0f : 0.0f; } + // As many frames as the validation set, equally spaced over the rest of the sweep [early, images). + std::vector LateValidationFrames(int early, int images) { + std::vector late = select_equally_spaced_image_ordinals(images - early, std::min(images - early, 60)); + for (int &ordinal : late) + ordinal += early; + return late; + } + // === Is the larger of two integer-related cells real, or a spurious axis multiple? === // // The validation-FRAME count cannot answer it: a spurious multiple indexes every frame its true @@ -361,6 +369,7 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { p.observer->OnPhase("Rotation indexing (first pass)"); Step(!p.geometry_prepass && prepass_detector_geometry_ ? "Indexing at the refined geometry" : "Indexing"); + near_tie_carried_ = false; FirstPassRun run(*this, p); using FirstPass = FirstPassRun::FirstPass; const int n_validation = static_cast(run.validation.size()); @@ -407,6 +416,25 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { } const FirstPassRun::Verdict verdict = run.MeasureVerdict(best, first_pass_key, memo_hit); + // What the pass ended on, for the re-check of a near tie it decided early to compare with. The + // numbers are printed exactly (shortest round trip), so equal means equal. + { + const auto vec = [](const Coord &c) { return fmt::format("{} {} {}", c.x, c.y, c.z); }; + const auto centre = [](const std::optional> &c) { + return c ? fmt::format("{} {}", (*c)[0], (*c)[1]) : std::string("-"); + }; + const auto state = [](const FirstPassState &st) { + return fmt::format("{} {} {}", st.rings, st.seed, st.d_min.value_or(0.0f)); + }; + const CrystalLattice &l = best.result->lattice; + first_pass_outcome_ = fmt::format( + "'{}' {} {} | {} | {} | {} {} {} | {} {} {} | {} | {} | {} | {}", best.name, best.score, + best.larger_settled, vec(l.Vec0()), vec(l.Vec1()), vec(l.Vec2()), best.result->search_result.centering, + static_cast(best.result->search_result.system), experiment_.GetBeamX_pxl(), + experiment_.GetBeamY_pxl(), experiment_.GetDetectorDistance_mm(), state(run.committed), + centre(beam_center_alternative_), centre(beam_center_same_lattice_alternative_), + beam_center_alternative_state_ ? state(*beam_center_alternative_state_) : std::string("-")); + } // 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_) { @@ -438,6 +466,7 @@ Rugnux::FirstPassRun::FirstPassRun(Rugnux &rx, PipelineLocals &p) ? std::min(rx.config_.nthreads, std::max(8, 4 * get_gpu_count())) : rx.config_.nthreads), engines(std::max(spot_workers, 1)), + late_validation(LateValidationFrames(p.prepass_images, p.images_to_process)), validation_settings(rx.config_.spot_finding), all_consecutive(early_images), spread_ordinals(select_equally_spaced_image_ordinals(early_images, @@ -611,10 +640,15 @@ void Rugnux::FirstPassRun::RestoreBeamCenter(float x, float y) { // 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) { + return ValidationCloud(r, validation); +} + +std::vector Rugnux::FirstPassRun::ValidationCloud(const RotationIndexerResult &r, + const std::vector &frames) { std::vector cloud; if (!r.axis) return cloud; - for (const int ordinal : validation) { + for (const int ordinal : frames) { const auto it = spot_cache.find(ordinal); if (it == spot_cache.end()) continue; @@ -629,11 +663,16 @@ std::vector Rugnux::FirstPassRun::ValidationCloud(const RotationIndexerRe // How many validation frames a candidate indexes on the real per-image path (see validation). int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r) { + return CountIndexed(idx, r, validation); +} + +int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r, + const std::vector &frames) { 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); + PrefetchSpots(frames); // 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 @@ -642,11 +681,11 @@ int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexe // 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) { + ParallelFor(static_cast(frames.size()), + std::min(frames.size(), rx.config_.nthreads), [&](int i) { DataMessage m{}; - m.number = validation[i]; - m.spots = spot_cache.at(validation[i]); + m.number = frames[i]; + m.spots = spot_cache.at(frames[i]); if (idx.IndexFrameOnly(m, validation_settings)) ++count; }); @@ -666,18 +705,19 @@ int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexe // 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) { + const RotationIndexerResult &r, int shift, + const std::vector &frames) { idx.ForceRotationIndexerResult(r); - PrefetchSpots(validation); + PrefetchSpots(frames); 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) { + ParallelFor(static_cast(frames.size()), + std::min(frames.size(), rx.config_.nthreads), [&](int i) { DataMessage m{}; - m.number = validation[(i + shift) % validation.size()]; - m.spots = spot_cache.at(validation[i]); + m.number = frames[(i + shift) % frames.size()]; + m.spots = spot_cache.at(frames[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); @@ -704,17 +744,22 @@ std::pair Rugnux::FirstPassRun::PooledOnLattice(IndexAndRefine // 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 { + return MeasurePooled(idx, r, validation); +} + +auto Rugnux::FirstPassRun::MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r, + const std::vector &frames) -> PooledEvidence { PooledEvidence e; - std::tie(e.on_lattice, e.reflections) = PooledOnLattice(idx, r, 0); + std::tie(e.on_lattice, e.reflections) = PooledOnLattice(idx, r, 0, frames); const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings(); - for (const int ordinal : validation) + for (const int ordinal : frames) for (const auto &sp : spot_cache.at(ordinal)) if (index_ice_rings || !sp.ice_ring) e.spots++; - const int n = static_cast(validation.size()); + const int n = static_cast(frames.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)); + const auto [spots_on, reflections_on] = PooledOnLattice(idx, r, std::max(1, k * n / 7), frames); chance.push_back(spots_on); e.reflections_by_chance += reflections_on; } @@ -846,9 +891,19 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin if (result.unconstrained && score < majority) { RotationIndexerResult alt = *result.unconstrained; const int alt_score = CountIndexed(idx, alt); - NearTieGate("metric-symmetry demotion", fmt::format("scheme '{}' unconstrained", name), - alt_score, majority); - if (alt_score > majority) { + bool works = alt_score > majority; + if (NearTieGate("metric-symmetry demotion", fmt::format("scheme '{}' unconstrained", name), + alt_score, majority) + && DecideLate()) { + const int late = CountIndexed(idx, alt, late_validation); + const bool early = works; + works = LateOverGate(alt_score, late, majority); + LogDecidedLate("metric-symmetry demotion", + fmt::format("scheme '{}' unconstrained {}/{} of the later frames, {}", name, late, + late_validation.size(), works ? "demoting" : "keeping the symmetry"), + early, works); + } + if (works) { 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, @@ -877,11 +932,24 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin 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; - if (bp.result.has_value() && !SameLattice(result, *bp.result)) - NearTieFrames("scheme choice", fmt::format("'{}'", name), score, - fmt::format("'{}'", bp.name), bp.score); + // A later scheme wins if it indexes clearly more frames (>10%). On a near tie between two + // lattices, of all the frames: the counts below are then the two over both sets. + int count = score, bp_count = bp.score; + if (bp.result.has_value() && !SameLattice(result, *bp.result) + && NearTieFrames("scheme choice", fmt::format("'{}'", name), score, + fmt::format("'{}'", bp.name), bp.score) + && DecideLate()) { + const int late = CountIndexed(idx, result, late_validation); + const int bp_late = CountIndexed(idx, *bp.result, late_validation); + count += late; + bp_count += bp_late; + LogDecidedLate("scheme choice", fmt::format("'{}' {}/{} of the later frames, '{}' {}/{}", name, late, + late_validation.size(), bp.name, bp_late, + late_validation.size()), + static_cast(score) > bp.score * 1.1f + 0.5f, + static_cast(count) > bp_count * 1.1f + 0.5f); + } + const bool clearly_more = static_cast(count) > bp_count * 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 @@ -895,7 +963,7 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin 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 bool tied = static_cast(count) >= bp_count * 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 @@ -918,32 +986,51 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin 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; - NearTieCounts("axis-harmonic arbiter", "larger cell", static_cast(hkl.size()), - "smaller cell", static_cast(n_smaller)); - NearTieCounts("axis-harmonic arbiter (own sub-lattice)", "larger cell", - static_cast(hkl_p.size()), "sub-lattice", static_cast(n_sub)); + struct Counts { size_t larger = 0, smaller = 0, larger_p = 0, sub = 0; AddedClassEvidence ev; }; + const auto count_spots = [&](const std::vector &frames) { + Counts c; + const std::vector cloud = ValidationCloud(larger, frames); + c.larger = IndexedMillerIndices(larger.lattice, cloud, index_tol).size(); + c.smaller = IndexedMillerIndices(smaller.lattice, ValidationCloud(smaller, frames), + index_tol).size(); + const auto hkl_p = IndexedMillerIndices(larger_p, cloud, index_tol); + c.larger_p = hkl_p.size(); + c.ev = AddedClassOccupancy(hkl_p, static_cast(nearest)); + const CrystalLattice sub = SubLattice(larger_p, c.ev, static_cast(nearest)).NiggliReduce(); + c.sub = IndexedMillerIndices(sub, cloud, index_tol).size(); + return c; + }; + Counts c = count_spots(validation); + // Both lines are said, whichever is the tie. + const bool tie_smaller = NearTieCounts("axis-harmonic arbiter", "larger cell", + static_cast(c.larger), "smaller cell", + static_cast(c.smaller)); + const bool tie_sub = NearTieCounts("axis-harmonic arbiter (own sub-lattice)", "larger cell", + static_cast(c.larger_p), "sub-lattice", + static_cast(c.sub)); + if ((tie_smaller || tie_sub) && DecideLate()) { + const bool early = c.larger > c.smaller && c.larger_p > c.sub; + PrefetchSpots(late_validation); + c = count_spots(AllValidationFrames()); + LogDecidedLate("axis-harmonic arbiter", + fmt::format("larger cell {} spots against the smaller cell's {}, {} against " + "its sub-lattice's {}", c.larger, c.smaller, c.larger_p, c.sub), + early, c.larger > c.smaller && c.larger_p > c.sub); + } + const auto &ev = c.ev; + const bool larger_is_real = c.larger > c.smaller && c.larger_p > c.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, + name, vol, bp.name, bp.vol, ratio, c.larger, c.smaller, + c.larger_p, ev.u, ev.v, ev.w, c.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 " @@ -1193,46 +1280,57 @@ auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int m // 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. + // + // The absolute bar is a share of the validation frames, and where the rung that indexes the + // most is a near tie at it, the bar is asked of it on all the frames, at its own settings. + std::vector clears_bar(rungs.size()); + for (size_t i = 0; i < rungs.size(); i++) + clears_bar[i] = rung_score[i] >= min_score; + int top = -1; + for (size_t i = 0; i < rungs.size(); i++) + if (rung_score[i] >= 0 && (top < 0 || rung_score[i] > rung_score[top])) + top = static_cast(i); + if (top >= 0) { + // Near tie (log only) on the gain. + const double n = static_cast(validation.size()); + const auto var = [n](int k) { + const double p = (k + 0.5) / (n + 1.0); + return n * p * (1.0 - p); + }; + NearTie("ladder adoption", fmt::format("best rung gains {}/{}", rung_score[top] - score_before, + validation.size()), + fmt::format("min gain {}", min_gain), rung_score[top] - score_before - min_gain, + std::sqrt(var(rung_score[top]) + var(std::max(0, score_before))), + "binomial on validation frames"); + if (min_score > 0 && NearTieGate("ladder adoption (bar)", "best rung", rung_score[top], min_score) + && DecideLate()) { + ApplyState({rungs[top].rings, rungs[top].seed, rungs[top].d_min}); + const int late = CountIndexed(*indexer, *rung_pass[top].result, late_validation); + // At least min_score of the validation frames is more than min_score - 1 of them. + clears_bar[top] = LateOverGate(rung_score[top], late, min_score - 1); + LogDecidedLate("ladder adoption (bar)", fmt::format("best rung {}/{} of the later frames, {}", late, + late_validation.size(), + clears_bar[top] ? "clears it" : "does not clear it"), + rung_score[top] >= min_score, clears_bar[top]); + } + } 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) + if (rung_score[i] < 0 || rung_score[i] - score_before < min_gain || !clears_bar[i]) 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); } - - // Near ties (log only): the rung that indexes the most against the bars it has to clear, and - // against the best rung that found a different lattice. - { - int top = -1; + // Near tie (log only): the winner against the best rung that found a different lattice. + if (won >= 0) { + int rival = -1; for (size_t i = 0; i < rungs.size(); i++) - if (rung_score[i] >= 0 && (top < 0 || rung_score[i] > rung_score[top])) - top = static_cast(i); - if (top >= 0) { - const double n = static_cast(validation.size()); - const auto var = [n](int k) { - const double p = (k + 0.5) / (n + 1.0); - return n * p * (1.0 - p); - }; - NearTie("ladder adoption", fmt::format("best rung gains {}/{}", rung_score[top] - score_before, - validation.size()), - fmt::format("min gain {}", min_gain), rung_score[top] - score_before - min_gain, - std::sqrt(var(rung_score[top]) + var(std::max(0, score_before))), - "binomial on validation frames"); - if (min_score > 0) - NearTieGate("ladder adoption (bar)", "best rung", rung_score[top], min_score); - } - if (won >= 0) { - int rival = -1; - for (size_t i = 0; i < rungs.size(); i++) - if (static_cast(i) != won && rung_score[i] >= 0 && rung_score[i] - score_before >= min_gain - && rung_score[i] >= min_score && !SameLattice(*rung_pass[i].result, *rung_pass[won].result) - && (rival < 0 || rung_score[i] > rung_score[rival])) - rival = static_cast(i); - if (rival >= 0) - NearTieFrames("ladder rung choice", "winner", rung_score[won], "other lattice", rung_score[rival]); - } + if (static_cast(i) != won && rung_score[i] >= 0 && rung_score[i] - score_before >= min_gain + && clears_bar[i] && !SameLattice(*rung_pass[i].result, *rung_pass[won].result) + && (rival < 0 || rung_score[i] > rung_score[rival])) + rival = static_cast(i); + if (rival >= 0) + NearTieFrames("ladder rung choice", "winner", rung_score[won], "other lattice", rung_score[rival]); } // ...and it must not be an axis SUB-MULTIPLE of a cell another rung found. The frame @@ -1269,7 +1367,7 @@ auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int m // 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) + || rung_score[i] - score_before < min_gain || !clears_bar[i]) continue; const double ratio = rung_pass[i].vol / won_vol; const double nearest = std::round(ratio); @@ -1587,51 +1685,54 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass // Which first pass a near-tie line comes from: the geometry pre-pass, the canonical pass, or a probe // that only scores a lattice. std::string Rugnux::FirstPassRun::PassTag() const { - return geometry_prepass ? "pass 1" : rx.indexing_probe_only_ ? "probe" : "pass 2"; + const char *pass = geometry_prepass ? "pass 1" : rx.indexing_probe_only_ ? "probe" : "pass 2"; + return rx.near_tie_recheck_ ? std::string(pass) + " re-check" : pass; } -void Rugnux::FirstPassRun::NearTie(const std::string &decision, const std::string &a, const std::string &b, +bool Rugnux::FirstPassRun::NearTie(const std::string &decision, const std::string &a, const std::string &b, double margin, double noise, const char *model) const { - if (std::abs(margin) <= NEAR_TIE_SIGMA * noise) - logger.Info("near tie: {} {} vs {} margin {:.3g} (noise {:.3g}, {}; {}) - would carry both", - decision, a, b, margin, noise, model, PassTag()); + if (std::abs(margin) > NEAR_TIE_SIGMA * noise) + return false; + logger.Info("near tie: {} {} vs {} margin {:.3g} (noise {:.3g}, {}; {}) - would carry both", + decision, a, b, margin, noise, model, PassTag()); + return true; } std::string Rugnux::FirstPassRun::Frames(int score) const { return score < 0 ? std::string("no lattice") : fmt::format("{}/{}", score, validation.size()); } -void Rugnux::FirstPassRun::NearTieFrames(const std::string &decision, const std::string &a, int ka, +bool Rugnux::FirstPassRun::NearTieFrames(const std::string &decision, const std::string &a, int ka, const std::string &b, int kb) const { // Two lattices that both index fewer frames than the bar the run keeps a lattice on (a sixth of // the validation frames, see FirstPassRotationIndexing) are not two hypotheses about the crystal. if (std::max(ka, kb) < static_cast(validation.size()) / 6) - return; + return false; const double n = static_cast(validation.size()); const auto var = [n](int k) { const double p = (std::max(k, 0) + 0.5) / (n + 1.0); return n * p * (1.0 - p); }; - NearTie(decision, fmt::format("{} {}", a, Frames(ka)), fmt::format("{} {}", b, Frames(kb)), - std::max(ka, 0) - std::max(kb, 0), std::sqrt(var(ka) + var(kb)), "binomial on validation frames"); + return NearTie(decision, fmt::format("{} {}", a, Frames(ka)), fmt::format("{} {}", b, Frames(kb)), + std::max(ka, 0) - std::max(kb, 0), std::sqrt(var(ka) + var(kb)), "binomial on validation frames"); } -void Rugnux::FirstPassRun::NearTieGate(const std::string &decision, const std::string &what, int score, +bool Rugnux::FirstPassRun::NearTieGate(const std::string &decision, const std::string &what, int score, double gate) const { const double n = static_cast(validation.size()); const double p = (std::max(score, 0) + 0.5) / (n + 1.0); - NearTie(decision, fmt::format("{} {}", what, Frames(score)), - fmt::format("gate {:.1f}", gate), std::max(score, 0) - gate, std::sqrt(n * p * (1.0 - p)), - "binomial on validation frames"); + return NearTie(decision, fmt::format("{} {}", what, Frames(score)), + fmt::format("gate {:.1f}", gate), std::max(score, 0) - gate, std::sqrt(n * p * (1.0 - p)), + "binomial on validation frames"); } -void Rugnux::FirstPassRun::NearTieCounts(const std::string &decision, const std::string &a, int64_t na, +bool Rugnux::FirstPassRun::NearTieCounts(const std::string &decision, const std::string &a, int64_t na, const std::string &b, int64_t nb) const { - NearTie(decision, fmt::format("{} {}", a, na), fmt::format("{} {}", b, nb), - static_cast(na - nb), std::sqrt(static_cast(na + nb)), "Poisson"); + return NearTie(decision, fmt::format("{} {}", a, na), fmt::format("{} {}", b, nb), + static_cast(na - nb), std::sqrt(static_cast(na + nb)), "Poisson"); } -void Rugnux::FirstPassRun::NearTiePooled(const std::string &decision, const std::string &a, +bool Rugnux::FirstPassRun::NearTiePooled(const std::string &decision, const std::string &a, const PooledEvidence &ea, const std::string &b, const PooledEvidence &eb) const { const auto excess = [](const PooledEvidence &e) { @@ -1639,6 +1740,39 @@ void Rugnux::FirstPassRun::NearTiePooled(const std::string &decision, const std: / static_cast(std::max(1, e.spots)); }; const auto [difference, sigma] = ValidationEvidenceDifference(ea, eb); - NearTie(decision, fmt::format("{} {:.1f}%", a, excess(ea)), fmt::format("{} {:.1f}%", b, excess(eb)), - -difference, sigma, "binomial on pooled validation spots"); + return NearTie(decision, fmt::format("{} {:.1f}%", a, excess(ea)), fmt::format("{} {:.1f}%", b, excess(eb)), + -difference, sigma, "binomial on pooled validation spots"); +} + +bool Rugnux::FirstPassRun::DecideLate() const { + if (late_validation.empty()) + return false; + if (rx.decide_late_) + return true; + rx.near_tie_carried_ = true; + rx.near_tie_rechecks_->carried = true; + return false; +} + +std::vector Rugnux::FirstPassRun::AllValidationFrames() const { + std::vector all = validation; + all.insert(all.end(), late_validation.begin(), late_validation.end()); + return all; +} + +bool Rugnux::FirstPassRun::LateOverGate(int score, int late, double gate) const { + // A lattice that indexes none of the later frames at all is not judged on them: where the first + // part indexes, that is a crystal that has left the beam or stopped diffracting, and frames + // with no crystal on them say nothing about the share of the crystal's frames a cell indexes. + if (late == 0) + return score > gate; + const double n = static_cast(validation.size()); + return score + late > gate * (n + static_cast(late_validation.size())) / n; +} + +void Rugnux::FirstPassRun::LogDecidedLate(const std::string &decision, const std::string &what, bool early, + bool late) const { + logger.Info("near tie decided late: {} - {} ({} validation frames and {} from the rest of the sweep; {}) - {}", + decision, what, validation.size(), late_validation.size(), PassTag(), + early == late ? "as on the validation frames alone" : "overturns the validation frames alone"); } diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index b91870ff0..4d3ddfb35 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -122,10 +122,15 @@ struct Rugnux::FirstPassRun { void PrefetchSpots(const std::vector &ordinals); void TryBeamCenter(float x, float y); void RestoreBeamCenter(float x, float y); + // Each on the validation frames, or on the frames given. std::vector ValidationCloud(const RotationIndexerResult &r); + std::vector ValidationCloud(const RotationIndexerResult &r, const std::vector &frames); int CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r); - std::pair PooledOnLattice(IndexAndRefine &idx, const RotationIndexerResult &r, int shift); + int CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector &frames); + std::pair PooledOnLattice(IndexAndRefine &idx, const RotationIndexerResult &r, int shift, + const std::vector &frames); PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r); + PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector &frames); int64_t MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r); std::unique_ptr FeedScheme(IndexerThreadPool &pool, const std::vector &ordinals, const DiffractionExperiment *x = nullptr); @@ -141,31 +146,50 @@ struct Rugnux::FirstPassRun { void StartCentreSpeculation(); void StartShortAxisSpeculation(); - // --- Near ties (log only; nothing here changes a decision) --- + // --- Near ties --- // // A decision whose margin lies within NEAR_TIE_SIGMA standard deviations of its own noise was - // taken on noise: the line names the alternative that would have to be carried forward beside - // it. The noise is the comparison's own - binomial for a count of validation frames (on + // taken on noise: the line names the alternative that has to be carried forward beside it. The + // noise is the comparison's own - binomial for a count of validation frames (on // n = validation.size(), p = (k + 1/2) / (n + 1)), Poisson for counts of spots or reflections, // and for two pooled spot shares the binomial standard error ValidationEvidencePrefers uses. + // Each returns whether the decision is a near tie. static constexpr double NEAR_TIE_SIGMA = 2.0; std::string PassTag() const; // A frame count as the log gives it: "k/N", or "no lattice" for a pass that found none. std::string Frames(int score) const; - void NearTie(const std::string &decision, const std::string &a, const std::string &b, + bool NearTie(const std::string &decision, const std::string &a, const std::string &b, double margin, double noise, const char *model) const; - // Two frame counts, a against b; nothing where neither reaches the run's floor of a sixth. - void NearTieFrames(const std::string &decision, const std::string &a, int ka, + // Two frame counts, a against b; never a tie where neither reaches the run's floor of a sixth. + bool NearTieFrames(const std::string &decision, const std::string &a, int ka, const std::string &b, int kb) const; // A frame count against the gate it is compared with. - void NearTieGate(const std::string &decision, const std::string &what, int score, double gate) const; + bool NearTieGate(const std::string &decision, const std::string &what, int score, double gate) const; // Two counts of spots or reflections. - void NearTieCounts(const std::string &decision, const std::string &a, int64_t na, + bool NearTieCounts(const std::string &decision, const std::string &a, int64_t na, const std::string &b, int64_t nb) const; // Two lattices' pooled validation evidence, a the standing one. - void NearTiePooled(const std::string &decision, const std::string &a, const PooledEvidence &ea, + bool NearTiePooled(const std::string &decision, const std::string &a, const PooledEvidence &ea, const std::string &b, const PooledEvidence &eb) const; + // A near tie is decided early and checked late. The pass acts on the validation frames at once; + // the comparison made again, unchanged, on the validation frames and late_validation together - + // a gate set at a share of the validation frames at the same share of all of them - is the + // re-check, which Run() makes once the rest of the sweep is in (Rugnux::decide_late_). A + // decision that is not a near tie never reads later frames; nor does one on a sweep the first + // part already covers, where there is nothing later to read (late_validation, below). + // + // Whether this near tie is decided on all the frames now: on the re-check, and on the run Run() + // makes again where a re-check ended elsewhere. Otherwise it notes that the pass carried one. + bool DecideLate() const; + // The validation frames followed by late_validation. + std::vector AllValidationFrames() const; + // A frame count on the validation frames, plus `late` on late_validation, over the gate set at + // `gate` of the validation frames: the same share of all of them. + bool LateOverGate(int score, int late, double gate) const; + // `early` and `late` are the decision on the validation frames alone and on all of them. + void LogDecidedLate(const std::string &decision, const std::string &what, bool early, bool late) const; + // --- State --- Rugnux ℞ @@ -237,6 +261,9 @@ struct Rugnux::FirstPassRun { // 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)); + // As many frames of the rest of the sweep, [early_images, images_to_process), equally spaced over + // it: read only to decide a near tie (see there). + const std::vector late_validation; const float index_tol = rx.experiment_.GetIndexingSettings().GetTolerance(); diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp index 3c1d8f8cc..765007839 100644 --- a/rugnux/RugnuxFirstPassRescues.cpp +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -219,9 +219,23 @@ auto Rugnux::FirstPassRun::LongAxisRescue(const FirstPass &standing) -> std::opt constrained = PickBest(*indexer_pool, *indexer); } rx.experiment_.SetUnitCell(std::nullopt); // leave the space-group / cell determination de-novo - if (constrained.result.has_value()) - NearTieFrames("long-axis rescue", "re-index", constrained.score, "standing", standing.score); - if (constrained.result.has_value() && constrained.score > std::max(standing.score, 0)) + if (!constrained.result.has_value()) + return std::nullopt; + // A standing pass that found no lattice counts as 0 frames. + bool more = constrained.score > std::max(standing.score, 0); + if (standing.result.has_value() && !SameLattice(*constrained.result, *standing.result) + && NearTieFrames("long-axis rescue", "re-index", constrained.score, "standing", standing.score) + && DecideLate()) { + const int late = CountIndexed(*indexer, *constrained.result, late_validation); + const int standing_late = CountIndexed(*indexer, *standing.result, late_validation); + const bool early = more; + more = constrained.score + late > standing.score + standing_late; + LogDecidedLate("long-axis rescue", fmt::format("re-index {}/{} of the later frames, standing {}/{}", late, + late_validation.size(), standing_late, + late_validation.size()), + early, more); + } + if (more) return constrained; } return std::nullopt; @@ -499,8 +513,21 @@ auto Rugnux::FirstPassRun::CentreByPooledSpots(FirstPass best, FirstPass alt, co if (best.result.has_value()) at_header = MeasurePooled(*indexer, *best.result); long_axis_asked = true; - if (alt.result.has_value() && best.result.has_value()) - NearTiePooled("beam centre (pooled spots)", "file centre", at_header, "measured centre", at_measured); + if (alt.result.has_value() && best.result.has_value() + && NearTiePooled("beam centre (pooled spots)", "file centre", at_header, "measured centre", at_measured) + && DecideLate()) { + // Each centre on its own spots, as above. + const bool early = ValidationEvidenceBeatsChance(at_measured) + && ValidationEvidencePrefers(at_header, at_measured); + const std::vector all = AllValidationFrames(); + at_header = MeasurePooled(*indexer, *best.result, all); + TryBeamCenter(measured_x, measured_y); + at_measured = MeasurePooled(*indexer, *alt.result, all); + RestoreBeamCenter(header_x, header_y); + LogDecidedLate("beam centre (pooled spots)", "the pooled spots counted again on all the frames at both centres", + early, ValidationEvidenceBeatsChance(at_measured) + && ValidationEvidencePrefers(at_header, at_measured)); + } if (alt.result.has_value() && ValidationEvidenceBeatsChance(at_measured) && ValidationEvidencePrefers(at_header, at_measured)) { TryBeamCenter(measured_x, measured_y); @@ -941,15 +968,35 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { } bool adopt = false; std::string verdict = " - none indexes a majority of the frames and as many as the standing cell"; - if (sub_best) - NearTieGate("sub-lattice majority gate", "best sub-lattice", sub_best->score, - static_cast(validation.size()) / 2); - 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); + const int majority = static_cast(validation.size()) / 2; + bool nominated = sub_best && sub_best->score >= best.score && sub_best->score > majority; + // Carried only where the gate is what decides: a sub-lattice that indexes fewer frames than + // the standing cell is not nominated whatever the gate says. + if (sub_best && NearTieGate("sub-lattice majority gate", "best sub-lattice", sub_best->score, majority) + && sub_best->score >= best.score && DecideLate()) { + const int late = CountIndexed(*indexer, *sub_best->result, late_validation); + const int standing_late = CountIndexed(*indexer, *best.result, late_validation); + const bool early = nominated; + nominated = sub_best->score + late >= best.score + standing_late + && LateOverGate(sub_best->score, late, majority); + LogDecidedLate("sub-lattice majority gate", + fmt::format("best sub-lattice {}/{} of the later frames, the standing cell {}/{}", late, + late_validation.size(), standing_late, late_validation.size()), + early, nominated); + } + if (nominated) { + PooledEvidence at_standing = MeasurePooled(*indexer, *best.result); + PooledEvidence at_sub = MeasurePooled(*indexer, *sub_best->result); + if (NearTiePooled("sub-lattice hypothesis", "standing cell", at_standing, "sub-lattice", at_sub) + && DecideLate()) { + const bool early = ValidationEvidencePrefers(at_standing, at_sub); + const std::vector all = AllValidationFrames(); + at_standing = MeasurePooled(*indexer, *best.result, all); + at_sub = MeasurePooled(*indexer, *sub_best->result, all); + LogDecidedLate("sub-lattice hypothesis", "the pooled spots counted again on all the frames", + early, ValidationEvidencePrefers(at_standing, at_sub)); + } adopt = ValidationEvidencePrefers(at_standing, at_sub); - NearTiePooled("sub-lattice hypothesis", "standing cell", at_standing, "sub-lattice", 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, @@ -1027,9 +1074,19 @@ auto Rugnux::FirstPassRun::Pass1Fallback(FirstPass best) -> FirstPass { const int pass1_score = (best.score < floor || !same_lattice) ? CountIndexed(*indexer, pass1_lattice) : -1; NearTieGate("pass-1 fallback floor", "re-index", best.score, floor); - if (pass1_score >= 0 && !same_lattice) - NearTieFrames("pass-1 fallback", "re-index", best.score, "pass-1 lattice", pass1_score); - if (best.score < floor || pass1_score > best.score) { + bool pass1_more = pass1_score > best.score; + if (pass1_score >= 0 && !same_lattice + && NearTieFrames("pass-1 fallback", "re-index", best.score, "pass-1 lattice", pass1_score) + && DecideLate()) { + const int late = CountIndexed(*indexer, *best.result, late_validation); + const int pass1_late = CountIndexed(*indexer, pass1_lattice, late_validation); + const bool early = pass1_more; + pass1_more = pass1_score + pass1_late > best.score + late; + LogDecidedLate("pass-1 fallback", fmt::format("re-index {}/{} of the later frames, pass-1 lattice {}/{}", late, + late_validation.size(), pass1_late, late_validation.size()), + early, pass1_more); + } + if (best.score < floor || pass1_more) { 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 "