Merge branch 'partb-decisions' into rc175: near ties decided early, re-checked late after the run, the run made again on an overturn

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-10-11 03:36:02 +02:00
commit 65fde24fb6
5 files changed
+453 -124

No files matched your search

+77 -4
View File
@@ -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<Rugnux>(*this);
again->near_tie_rechecks_ = std::make_shared<NearTieRechecks>();
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<double>(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<Rugnux> 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<Rugnux> 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<Rugnux>(*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).
+38
View File
@@ -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<Rugnux> 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<bool> carried = false; // a first pass decided a near tie early
std::vector<Pending> pending;
};
std::shared_ptr<NearTieRechecks> near_tie_rechecks_ = std::make_shared<NearTieRechecks>();
// 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<Rugnux> 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
+228 -94
View File
@@ -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<int> LateValidationFrames(int early, int images) {
std::vector<int> 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<int>(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<std::array<float, 2>> &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<int>(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<Coord> Rugnux::FirstPassRun::ValidationCloud(const RotationIndexerResult &r) {
return ValidationCloud(r, validation);
}
std::vector<Coord> Rugnux::FirstPassRun::ValidationCloud(const RotationIndexerResult &r,
const std::vector<int> &frames) {
std::vector<Coord> 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<Coord> 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<int> &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<int> count = 0;
ParallelFor(static_cast<int>(validation.size()),
std::min<size_t>(validation.size(), rx.config_.nthreads), [&](int i) {
ParallelFor(static_cast<int>(frames.size()),
std::min<size_t>(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<int64_t, int64_t> Rugnux::FirstPassRun::PooledOnLattice(IndexAndRefine &idx,
const RotationIndexerResult &r, int shift) {
const RotationIndexerResult &r, int shift,
const std::vector<int> &frames) {
idx.ForceRotationIndexerResult(r);
PrefetchSpots(validation);
PrefetchSpots(frames);
const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings();
std::mutex reflections_mutex;
std::set<std::array<int64_t, 3>> reflections;
std::atomic<int64_t> on_lattice = 0;
ParallelFor(static_cast<int>(validation.size()),
std::min<size_t>(validation.size(), rx.config_.nthreads), [&](int i) {
ParallelFor(static_cast<int>(frames.size()),
std::min<size_t>(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<int64_t, int64_t> 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<int> &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<int>(validation.size());
const int n = static_cast<int>(frames.size());
std::vector<int64_t> 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<float>(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<float>(score) > bp.score * 1.1f + 0.5f,
static_cast<float>(count) > bp_count * 1.1f + 0.5f);
}
const bool clearly_more = static_cast<float>(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<float>(score) >= bp.score * 0.9f - 0.5f;
const bool tied = static_cast<float>(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<Coord> 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<int>(nearest));
const CrystalLattice sub = SubLattice(larger_p, ev, static_cast<int>(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<int64_t>(hkl.size()),
"smaller cell", static_cast<int64_t>(n_smaller));
NearTieCounts("axis-harmonic arbiter (own sub-lattice)", "larger cell",
static_cast<int64_t>(hkl_p.size()), "sub-lattice", static_cast<int64_t>(n_sub));
struct Counts { size_t larger = 0, smaller = 0, larger_p = 0, sub = 0; AddedClassEvidence ev; };
const auto count_spots = [&](const std::vector<int> &frames) {
Counts c;
const std::vector<Coord> 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<int>(nearest));
const CrystalLattice sub = SubLattice(larger_p, c.ev, static_cast<int>(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<int64_t>(c.larger), "smaller cell",
static_cast<int64_t>(c.smaller));
const bool tie_sub = NearTieCounts("axis-harmonic arbiter (own sub-lattice)", "larger cell",
static_cast<int64_t>(c.larger_p), "sub-lattice",
static_cast<int64_t>(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<bool> 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<int>(i);
if (top >= 0) {
// Near tie (log only) on the gain.
const double n = static_cast<double>(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<int>(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<int>(i);
if (top >= 0) {
const double n = static_cast<double>(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<int>(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<int>(i);
if (rival >= 0)
NearTieFrames("ladder rung choice", "winner", rung_score[won], "other lattice", rung_score[rival]);
}
if (static_cast<int>(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<int>(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<int>(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<int>(validation.size()) / 6)
return;
return false;
const double n = static_cast<double>(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<double>(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<double>(na - nb), std::sqrt(static_cast<double>(na + nb)), "Poisson");
return NearTie(decision, fmt::format("{} {}", a, na), fmt::format("{} {}", b, nb),
static_cast<double>(na - nb), std::sqrt(static_cast<double>(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<double>(std::max<int64_t>(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<int> Rugnux::FirstPassRun::AllValidationFrames() const {
std::vector<int> 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<double>(validation.size());
return score + late > gate * (n + static_cast<double>(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");
}
+37 -10
View File
@@ -122,10 +122,15 @@ struct Rugnux::FirstPassRun {
void PrefetchSpots(const std::vector<int> &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<Coord> ValidationCloud(const RotationIndexerResult &r);
std::vector<Coord> ValidationCloud(const RotationIndexerResult &r, const std::vector<int> &frames);
int CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r);
std::pair<int64_t, int64_t> PooledOnLattice(IndexAndRefine &idx, const RotationIndexerResult &r, int shift);
int CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector<int> &frames);
std::pair<int64_t, int64_t> PooledOnLattice(IndexAndRefine &idx, const RotationIndexerResult &r, int shift,
const std::vector<int> &frames);
PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r);
PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector<int> &frames);
int64_t MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r);
std::unique_ptr<RotationIndexer> FeedScheme(IndexerThreadPool &pool, const std::vector<int> &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<int> 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 &rx;
@@ -237,6 +261,9 @@ struct Rugnux::FirstPassRun {
// many - this is the discriminator that lets the best scheme win.
const std::vector<int> 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<int> late_validation;
const float index_tol = rx.experiment_.GetIndexingSettings().GetTolerance();
+73 -16
View File
@@ -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<int> 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<int>(validation.size()) / 2);
if (sub_best && sub_best->score >= best.score
&& sub_best->score > static_cast<int>(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<int>(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<int> 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 "