mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 81 Wiki Activity
Files
a7dc45b7ea46f597cc7a77a5e3e7ca90cb18a3a7
Jungfraujoch/rugnux/RugnuxFirstPass.cpp
T
leonarski_fandClaude Opus 5.5 a7dc45b7ea rugnux: ask the long-axis rescue at both beam centres before the pooled spots choose between them 
When neither the file's nor the background-measured beam centre indexes a majority of the
validation frames, the first pass chose between them on the pooled spots each centre's lattice
puts on itself over chance, and only then asked the long-axis rescue - at the chosen centre alone.
The lattices compared there are, by construction, not the crystal's, and the comparison is close
to a coin toss. With the decisions before pass 2 taken from the first part of the sweep (part B,
5b9144410), 8v4o (P6122, c = 545 A) tipped from the file's centre to the measured one, 3.6 px out
across the spindle, on 16.7% against 15.5%; the rescue there could not see the 545 A axis, and the
run shipped a P1 sub-cell (primitive volume 0.35x). At the file's centre the same rescue on the
same early frames indexes 60/60 (checked with --beam-center-check=off on the rc175 binary).

The rescue is now a lambda, asked at both centres in that branch so each centre is judged on the
best lattice it gives, and asked as before after the beam-centre check where that branch did not
already ask it. The early-frame principle of part B is unchanged.

8v4o: P1 2.62 A, CC1/2 0.920, R_model 0.553 -> P6122 2.12 A (pass), the file's centre kept.
p.mtz md5-identical to rc175 7c9a6c7d1 (GPU) on 5nw5, 5ky6, 5reo, 7atg, 9qw8, 9rci, cuhf2,
cytidine, insu_I_x06da_low_isa, kdp_x10sa_20keV, lcystine_x10sa_20keV/25keV, lyso_x06da_half_image,
myob_x06da_powder_2, myob_x06da_sparse, nothing_1, nothing_2 and the two private sets that reach
the same branch. 9rci (doubled-c alternative under tNCS) is a different mechanism: rc175 with
--prepass-fraction 1 also returns the doubled cell.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
2026-10-08 21:11:17 +02:00

2631 lines
172 KiB
C++
Raw Blame History

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