leonarski_fandClaude Opus 5.5
097668e674
The two-pass supercell guards (in the first-pass indexing and after the second pass) revert to pass 1's lattice whenever the second pass's primitive cell is markedly larger. On a crystal with a doubled axis pass 1's schemes did not tie on frames, so pass 1 never asked the axis-harmonic arbiter and kept the halved cell; pass 2's schemes tied, its arbiter established the larger cell on the spots (17048 against 16099 validation spots), and the guard then threw it away. The scheme comparison now records when the lattice it keeps is the larger cell of an axis-harmonic pair the arbiter settled, and both guards leave such a lattice in place. Where the arbiter did not settle the larger cell - every other set the guard fires on - nothing changes. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi
1417 lines
80 KiB
C++
1417 lines
80 KiB
C++
// 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"
|
|
#include "RugnuxFirstPass.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, for work that runs beside the first pass: spread
|
|
// over the cards other than card 0, which the first pass itself keeps busy (worker t on card
|
|
// 1 + t % (cards - 1); card 0 only on a one-card machine). 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)
|
|
: reader(reader), x(x), mapping(mapping), mask(mask), settings(settings), dataset(dataset),
|
|
start_image(start_image), stride(stride), 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, thread t on the card its engine lives on.
|
|
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] {
|
|
if (const int32_t gpus = get_gpu_count(); gpus > 1)
|
|
pin_gpu(1 + static_cast<int32_t>(t) % (gpus - 1));
|
|
else if (gpus == 1)
|
|
pin_gpu(0);
|
|
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;
|
|
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.
|
|
//
|
|
// That argument holds only on ONE cloud. Each scheme's cell comes with the geometry it refined,
|
|
// and the spots are placed in reciprocal space by that geometry, so counted on their own clouds
|
|
// the two cells also compare two geometries - and a scheme that refined on frames spread over the
|
|
// sweep places spots better than one that refined on a wedge, whichever cell it found (measured:
|
|
// the 'spread' cell won the count at 2x, 3x and on either side of the pair). So the larger cell
|
|
// is also asked against its OWN sub-lattice (SubLattice), on its own cloud, where only the added
|
|
// class can tell the two apart.
|
|
//
|
|
// 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.u + out.v + out.w == 0) {
|
|
out.occupancy = occ;
|
|
out.u = u; out.v = v; out.w = w;
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// The index-n sub-lattice the form (u,v,w) describes: the reflections with u h + v k + w l = 0
|
|
// (mod n), i.e. the lattice with the translation (u a + v b + w c) / n added - the cell the larger
|
|
// one would be if the class it adds were empty. One axis i whose coefficient is prime to n is
|
|
// replaced by that translation, scaled so its coefficient becomes 1.
|
|
CrystalLattice SubLattice(const CrystalLattice &latt, const AddedClassEvidence &ev, int n) {
|
|
const int c[3] = {ev.u, ev.v, ev.w};
|
|
const Coord e[3] = {latt.Vec0(), latt.Vec1(), latt.Vec2()};
|
|
for (int i = 0; i < 3; i++) {
|
|
if (std::gcd(c[i], n) != 1)
|
|
continue;
|
|
int k = 1;
|
|
while ((k * c[i]) % n != 1)
|
|
k++;
|
|
Coord t = e[i] / static_cast<float>(n);
|
|
for (int j = 0; j < 3; j++)
|
|
if (j != i)
|
|
t += e[j] * static_cast<float>((k * c[j]) % n) / static_cast<float>(n);
|
|
Coord v[3] = {e[0], e[1], e[2]};
|
|
v[i] = t;
|
|
return {v[0], v[1], v[2]};
|
|
}
|
|
return latt;
|
|
}
|
|
}
|
|
|
|
// The first pass of a rotation run, as the ordered list of its decisions: index the first part of the
|
|
// sweep on two sampling schemes, put the standing lattice through the rescues and the hypotheses that
|
|
// may replace it, accept or refuse it on the pooled validation spots, and leave it forced on the
|
|
// indexer. Each step is a member of FirstPassRun (RugnuxFirstPass.h); the comment at each says why it
|
|
// sits where it does - the order is part of the decision.
|
|
bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) {
|
|
auto &logger = p.logger;
|
|
// 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 (p.observer)
|
|
p.observer->OnPhase("Rotation indexing (first pass)");
|
|
Step(!p.geometry_prepass && prepass_detector_geometry_ ? "Indexing at the refined geometry" : "Indexing");
|
|
|
|
FirstPassRun run(*this, p);
|
|
using FirstPass = FirstPassRun::FirstPass;
|
|
const int n_validation = static_cast<int>(run.validation.size());
|
|
|
|
// 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_) {
|
|
run.committed = *beam_center_alternative_state_;
|
|
beam_center_alternative_state_.reset();
|
|
run.ApplyState(run.committed);
|
|
}
|
|
|
|
run.PrefetchCertainFrames();
|
|
FirstPass best = run.IndexSchemes();
|
|
best = run.AxisSignRescue(std::move(best));
|
|
bool long_axis_asked = false;
|
|
best = run.BeamCentreCheck(std::move(best), long_axis_asked);
|
|
// The long-axis rescue, at the centre the run now stands on.
|
|
if (!long_axis_asked)
|
|
if (auto rescued = run.LongAxisRescue(best))
|
|
best = *rescued;
|
|
// The end wedge, at the geometry the rescues settled on.
|
|
best = run.AddEndWedge(std::move(best), *p.indexer_pool, *p.indexer);
|
|
run.ReleaseSpotEngines(); // no more frames to read; give the cards their memory back
|
|
if (cancelled_)
|
|
return false;
|
|
|
|
run.RequireLattice(best);
|
|
best = run.ShortAxisHypothesis(std::move(best));
|
|
best = run.SubLatticeHypothesis(std::move(best));
|
|
if (!p.geometry_prepass && prepass_result_.has_value())
|
|
best = run.Pass1Fallback(std::move(best));
|
|
// The last resort, and ONLY on a run that has no usable lattice: the pass-1 fallback above
|
|
// has had its turn, and below this the run stops with an error. A run that indexes must
|
|
// not pay for it - measured, an unguarded call spent 41-76 % of the wall clock on a
|
|
// healthy dataset evaluating rungs that arithmetically could not be adopted (one crystal
|
|
// 41 s guarded against 175 s unguarded, same answer).
|
|
if (best.score < n_validation / 6) {
|
|
if (auto rescued = run.FirstPassLadder(best.score, n_validation / 6, 0))
|
|
best = *rescued;
|
|
}
|
|
|
|
const FirstPassRun::Verdict verdict = run.MeasureVerdict(best, first_pass_key, memo_hit);
|
|
// A pass run only to score a lattice has its score, whatever it is: a lattice that does not
|
|
// beat chance is a result for the comparison that asked, not a reason to stop the run.
|
|
if (indexing_probe_only_) {
|
|
logger.Info("Indexing probe: scheme '{}', {}/{} validation frames, {}/{} validation spots "
|
|
"= {:.1f}% against {:.1f}% at a wrong spindle angle", best.name, best.score,
|
|
n_validation, verdict.evidence.on_lattice, verdict.evidence.spots,
|
|
100.0 * verdict.pooled, 100.0 * verdict.pooled_chance);
|
|
return true;
|
|
}
|
|
run.RejectChanceLattice(best, verdict);
|
|
run.AdoptLattice(best, verdict);
|
|
if (!p.geometry_prepass && !postrefine_probe_only_)
|
|
run.StartLeftoverCensus(best);
|
|
if (best.result->tilt_walk)
|
|
run.ReportTiltWalk(best);
|
|
if (config_.measure_spot_budget)
|
|
run.MeasureSpotBudget(best);
|
|
if (!p.geometry_prepass && prepass_result_.has_value())
|
|
run.Pass1SupercellGuard(best);
|
|
return false;
|
|
}
|
|
|
|
Rugnux::FirstPassRun::FirstPassRun(Rugnux &rx, PipelineLocals &p)
|
|
: rx(rx), logger(p.logger), result(p.result), geometry_prepass(p.geometry_prepass),
|
|
start_image(p.start_image), early_images(p.prepass_images), dataset(p.dataset), mapping(p.mapping),
|
|
indexer_pool(p.indexer_pool), indexer(p.indexer), indexing_bound_A(p.indexing_bound_A),
|
|
census_future(p.census_future),
|
|
spot_workers(get_gpu_count() > 0
|
|
? std::min(rx.config_.nthreads, std::max(8, 4 * get_gpu_count()))
|
|
: rx.config_.nthreads),
|
|
engines(std::max(spot_workers, 1)),
|
|
validation_settings(rx.config_.spot_finding),
|
|
all_consecutive(early_images),
|
|
spread_ordinals(select_equally_spaced_image_ordinals(early_images,
|
|
rx.config_.rotation_indexing_image_count)),
|
|
first_pass_file_d_min(rx.config_.spot_finding.high_resolution_limit),
|
|
committed{false, 0, first_pass_file_d_min} {
|
|
validation_settings.indexing = true;
|
|
validation_settings.quick_integration = false;
|
|
std::iota(all_consecutive.begin(), all_consecutive.end(), 0);
|
|
schemes = {
|
|
{"spread", spread_ordinals},
|
|
{"wedge", all_consecutive},
|
|
};
|
|
end_wedge = {
|
|
{"end wedge", {all_consecutive.rbegin(), all_consecutive.rend()}},
|
|
};
|
|
}
|
|
|
|
// Mid-exposure goniometer angle for an image ordinal (matches the per-image path).
|
|
std::optional<float> Rugnux::FirstPassRun::RotAngle(int ordinal) const {
|
|
if (const auto g = rx.experiment_.GetGoniometer())
|
|
return g->GetAngle_deg(static_cast<float>(ordinal)) + g->GetWedge_deg() / 2.0f;
|
|
return std::nullopt;
|
|
}
|
|
|
|
// The spots of one frame, found with the caller's own analysis engine. A frame's spots are a
|
|
// pure function of that frame and the settings - the engine carries nothing from one image to
|
|
// the next, which is why the main image loop can give one to each of its workers - so it does
|
|
// not matter which engine, or which thread, finds them.
|
|
std::vector<SpotToSave> Rugnux::FirstPassRun::FindSpots(MXAnalysisWithoutFPGA &analysis,
|
|
AzimuthalIntegrationProfile &profile,
|
|
JFJochReaderRawImage &img, int ordinal) {
|
|
return FindFirstPassSpots(rx.reader_, rx.config_.spot_finding, *dataset, analysis, profile, img,
|
|
ordinal, start_image + ordinal * rx.config_.stride, logger);
|
|
}
|
|
|
|
void Rugnux::FirstPassRun::UseGpuOfWorker(size_t t) {
|
|
if (const int32_t gpus = get_gpu_count(); gpus > 0)
|
|
pin_gpu(static_cast<int32_t>(t) % gpus);
|
|
}
|
|
|
|
// Find the spots of every ordinal in the list that is not cached yet, on several workers, and
|
|
// put them in the cache.
|
|
//
|
|
// The cache ends up holding exactly what a serial loop would have put there: the results are
|
|
// collected into a vector indexed by position and inserted afterwards, so neither the spots of
|
|
// a frame nor the set of frames in the cache depends on how the workers interleaved.
|
|
//
|
|
// Worker t keeps its analysis engine for the whole first pass. Building one is a CUDA stream
|
|
// and a full set of device buffers - measured at 0.13 s against 4 ms for a frame - and the
|
|
// "wedge" scheme only discovers how far into the sweep it has to read as the frames arrive,
|
|
// so this is called several times and rebuilding them each time cost a quarter of the pass.
|
|
// The engine belongs to a card, not to a thread, so a later call only has to make that card
|
|
// current again; workers take a card by index so an engine always meets the one it was
|
|
// built on.
|
|
void Rugnux::FirstPassRun::PrefetchSpots(const std::vector<int> &ordinals) {
|
|
// What an earlier pass already found under the same key (rx.first_pass_spots_) is taken from
|
|
// there, at this pass's angles.
|
|
const std::vector<double> key = rx.SpotFindingKey(*spot_mapping);
|
|
std::map<int, std::vector<SpotToSave>> kept;
|
|
{
|
|
std::lock_guard lock(rx.first_pass_spots_->m);
|
|
for (const auto &[k, spots] : rx.first_pass_spots_->entries)
|
|
if (k == key)
|
|
for (const int ordinal : ordinals)
|
|
if (const auto it = spots.find(ordinal); it != spots.end())
|
|
kept.emplace(ordinal, it->second);
|
|
}
|
|
const DiffractionGeometry spot_geometry = rx.experiment_.GetDiffractionGeometry();
|
|
for (auto &[ordinal, spots] : kept)
|
|
StampSpotAngles(spots, spot_geometry);
|
|
std::vector<int> wanted;
|
|
for (const int ordinal : ordinals)
|
|
if (!spot_cache.contains(ordinal)) {
|
|
if (const auto it = kept.find(ordinal); it != kept.end() && !rx.config_.verify_first_pass_memo)
|
|
spot_cache.emplace(ordinal, it->second);
|
|
else
|
|
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] {
|
|
UseGpuOfWorker(t); // must precede engine construction
|
|
SpotEngine &e = engines[t];
|
|
if (!e.analysis) {
|
|
e.analysis = std::make_unique<MXAnalysisWithoutFPGA>(
|
|
rx.experiment_, *spot_mapping, rx.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] = FindSpots(*e.analysis, *e.profile, e.raw_image, wanted[i]);
|
|
}));
|
|
for (auto &f : futures)
|
|
f.get();
|
|
|
|
{
|
|
std::lock_guard lock(rx.first_pass_spots_->m);
|
|
auto &entries = rx.first_pass_spots_->entries;
|
|
auto entry = std::find_if(entries.begin(), entries.end(),
|
|
[&](const auto &e) { return e.first == key; });
|
|
if (entry == entries.end()) {
|
|
if (entries.size() == MAX_SPOT_KEYS)
|
|
entries.erase(entries.begin());
|
|
entries.emplace_back(key, std::map<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.
|
|
void Rugnux::FirstPassRun::ReleaseSpotEngines() {
|
|
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] {
|
|
UseGpuOfWorker(t);
|
|
engines[t] = SpotEngine{};
|
|
}));
|
|
for (auto &f : futures)
|
|
f.get();
|
|
}
|
|
|
|
// A beam-centre hypothesis, tried at its own mapping, engines and spot cache (see trial_mapping).
|
|
void Rugnux::FirstPassRun::TryBeamCenter(float x, float y) {
|
|
if (trial_mapping)
|
|
spot_cache.clear();
|
|
else
|
|
home_spot_cache.swap(spot_cache);
|
|
rx.experiment_.BeamX_pxl(x).BeamY_pxl(y);
|
|
ReleaseSpotEngines();
|
|
trial_mapping = std::make_unique<AzimuthalIntegrationMapping>(rx.experiment_, rx.pixel_mask_,
|
|
*rx.azint_geometry_);
|
|
spot_mapping = trial_mapping.get();
|
|
}
|
|
|
|
void Rugnux::FirstPassRun::RestoreBeamCenter(float x, float y) {
|
|
if (!trial_mapping)
|
|
return;
|
|
rx.experiment_.BeamX_pxl(x).BeamY_pxl(y);
|
|
ReleaseSpotEngines();
|
|
trial_mapping.reset();
|
|
spot_mapping = &mapping;
|
|
spot_cache = std::move(home_spot_cache);
|
|
home_spot_cache.clear();
|
|
}
|
|
|
|
// The validation frames' spots in reciprocal space, under a candidate's own refined geometry
|
|
// and axis. The same construction the rotation indexer makes on its accumulated frames, but
|
|
// over the validation set and over each frame's WHOLE spot list - the arbiter counts spots and
|
|
// wants the deep sample, where the first pass wants a bounded one. A superstructure layer is
|
|
// faintest exactly where the first pass's accumulation cap cuts.
|
|
std::vector<Coord> Rugnux::FirstPassRun::ValidationCloud(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 = RotAngle(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;
|
|
}
|
|
|
|
// How many validation frames a candidate indexes on the real per-image path (see validation).
|
|
int Rugnux::FirstPassRun::CountIndexed(IndexAndRefine &idx, const RotationIndexerResult &r) {
|
|
idx.ForceRotationIndexerResult(r);
|
|
// Fill the spot cache first. It is a plain map filled on demand, and the loop below runs the
|
|
// frames in parallel, where a lookup racing an insert is not something the map survives.
|
|
// After the first candidate this costs nothing.
|
|
PrefetchSpots(validation);
|
|
// Each frame's verdict is its own - the score is just how many of them index - so they run
|
|
// together. Refining one frame's geometry is a few Ceres solves and no GPU, and it is the
|
|
// same call the main image loop already makes from every one of its workers on this same
|
|
// IndexAndRefine, writing only unit_cells[] under its own mutex. With the candidate forced,
|
|
// GetLattice() returns it and the branch that would advance the indexer's own state is never
|
|
// reached. The offline solver stops on an iteration count rather than a clock
|
|
// (IndexAndRefine::RefineGeometryIfNeeded), so a loaded machine cannot change a verdict.
|
|
std::atomic<int> count = 0;
|
|
ParallelFor(static_cast<int>(validation.size()),
|
|
std::min<size_t>(validation.size(), rx.config_.nthreads), [&](int i) {
|
|
DataMessage m{};
|
|
m.number = validation[i];
|
|
m.spots = spot_cache.at(validation[i]);
|
|
if (idx.IndexFrameOnly(m, validation_settings))
|
|
++count;
|
|
});
|
|
return count.load();
|
|
}
|
|
|
|
// The run's acceptance test counts SPOTS, not frames. Summing each validation frame's spots
|
|
// before any threshold lets a four-spot frame contribute four spots of evidence instead of an
|
|
// automatic "no", and removes the spots-per-frame denominator that makes a frame count measure
|
|
// how sparse the pattern is rather than how good the lattice is: with 4 spots on an average
|
|
// frame no lattice, the true one included, can clear a bar of six indexed spots per frame.
|
|
//
|
|
// `shift` displaces each frame's spots onto another validation frame's spindle angle. The
|
|
// lattice, the spots, the geometry and the per-frame orientation polish are all unchanged and
|
|
// only the claim under test is removed - that these spots are at these angles - so what a
|
|
// displaced run still puts on the lattice is what chance and the polish give for free. That is
|
|
// the bar the real count is measured against, and it is measured rather than chosen, so no
|
|
// spots-per-frame number survives anywhere in the test.
|
|
std::pair<int64_t, int64_t> Rugnux::FirstPassRun::PooledOnLattice(IndexAndRefine &idx,
|
|
const RotationIndexerResult &r, int shift) {
|
|
idx.ForceRotationIndexerResult(r);
|
|
PrefetchSpots(validation);
|
|
const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings();
|
|
std::mutex reflections_mutex;
|
|
std::set<std::array<int64_t, 3>> reflections;
|
|
std::atomic<int64_t> on_lattice = 0;
|
|
ParallelFor(static_cast<int>(validation.size()),
|
|
std::min<size_t>(validation.size(), rx.config_.nthreads), [&](int i) {
|
|
DataMessage m{};
|
|
m.number = validation[(i + shift) % validation.size()];
|
|
m.spots = spot_cache.at(validation[i]);
|
|
on_lattice += idx.CountSpotsOnLattice(m, validation_settings);
|
|
// The spots the count took are the ones it gave a Miller index (none is 0,0,0).
|
|
std::unique_lock ul(reflections_mutex);
|
|
for (const auto &sp : m.spots)
|
|
if ((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0))
|
|
reflections.insert({sp.h, sp.k, sp.l});
|
|
});
|
|
return {on_lattice.load(), static_cast<int64_t>(reflections.size())};
|
|
}
|
|
|
|
// Displacements are a fixed fraction of the sweep, so the null is a function of the data alone
|
|
// and the same file gives the same verdict every time. Sevenths: no crystallographic rotation
|
|
// is a seventh of a turn, so a displacement cannot coincide with a symmetry about the spindle
|
|
// and hand the null the signal's own score. The median over five of them bounds what one
|
|
// unlucky displacement can do.
|
|
//
|
|
// This crystal's lattice is the one that puts materially more reflections on itself than a
|
|
// wrong angle does (ValidationEvidenceBeatsChance): every displaced angle's count is used, so
|
|
// the bar carries the null's own noise, and it carries no absolute count; the significance is
|
|
// the one-false-positive-in-a-thousand bar the spot-budget measurement is already set at.
|
|
//
|
|
// Two lattices are COMPARED on the same quantity, ValidationEvidencePrefers: the share of the
|
|
// pooled spots each one puts on itself over and above what a wrong spindle angle puts there.
|
|
// Each is scored against its own null, which is what makes the comparison fair between a cell
|
|
// and its axis harmonic.
|
|
auto Rugnux::FirstPassRun::MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r) -> PooledEvidence {
|
|
PooledEvidence e;
|
|
std::tie(e.on_lattice, e.reflections) = PooledOnLattice(idx, r, 0);
|
|
const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings();
|
|
for (const int ordinal : validation)
|
|
for (const auto &sp : spot_cache.at(ordinal))
|
|
if (index_ice_rings || !sp.ice_ring)
|
|
e.spots++;
|
|
const int n = static_cast<int>(validation.size());
|
|
std::vector<int64_t> chance;
|
|
for (int k = 1; k <= VALIDATION_NULL_DISPLACEMENTS; k++) {
|
|
const auto [spots_on, reflections_on] = PooledOnLattice(idx, r, std::max(1, k * n / 7));
|
|
chance.push_back(spots_on);
|
|
e.reflections_by_chance += reflections_on;
|
|
}
|
|
std::sort(chance.begin(), chance.end());
|
|
e.by_chance = chance[chance.size() / 2];
|
|
return e;
|
|
}
|
|
|
|
// How deep into an image's intensity-ordered spot list this lattice is still being seen - the
|
|
// measured spot budget. Same frames, same per-image path as CountIndexed above, but scoring the
|
|
// spots rather than the frames: each frame's spots are tallied by rank and the run keeps the
|
|
// ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than
|
|
// that admits detections that are not this crystal's reflections, and the frame gate then counts
|
|
// them against the frame - which is how a larger budget can integrate FEWER images.
|
|
int64_t Rugnux::FirstPassRun::MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r) {
|
|
idx.ForceRotationIndexerResult(r);
|
|
PrefetchSpots(validation);
|
|
std::mutex evidence_mutex;
|
|
std::vector<int64_t> indexed(rx.experiment_.GetMaxSpotCount(), 0);
|
|
std::vector<int64_t> counted(indexed.size(), 0);
|
|
const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings();
|
|
ParallelFor(static_cast<int>(validation.size()),
|
|
std::min<size_t>(validation.size(), rx.config_.nthreads), [&](int i) {
|
|
DataMessage m{};
|
|
m.number = validation[i];
|
|
m.spots = spot_cache.at(validation[i]);
|
|
if (!idx.IndexFrameOnly(m, validation_settings))
|
|
return;
|
|
std::unique_lock ul(evidence_mutex);
|
|
AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted);
|
|
});
|
|
return SpotBudgetFromEvidence(indexed, counted);
|
|
}
|
|
|
|
// Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to
|
|
// be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately
|
|
// (RunIndexing) so the two schemes' indexing can overlap.
|
|
//
|
|
// The feed itself stays in frame order and stops where it always did, but it keeps the cache a
|
|
// chunk of frames ahead of itself so those frames' spots are found on several workers. A chunk
|
|
// rather than the whole list because the "wedge" scheme stops on the accumulation cap, and how
|
|
// far into the sweep that is only becomes known as the frames arrive.
|
|
// x: the experiment the indexer is to read, experiment_ unless given.
|
|
std::unique_ptr<RotationIndexer> Rugnux::FirstPassRun::FeedScheme(IndexerThreadPool &pool,
|
|
const std::vector<int> &ordinals,
|
|
const DiffractionExperiment *x) {
|
|
auto ri = std::make_unique<RotationIndexer>(x ? *x : rx.experiment_, pool);
|
|
ri->MaxSpotsPerImage(first_pass_spots_per_image);
|
|
ri->VerifyMemo(rx.config_.verify_first_pass_memo);
|
|
for (size_t i = 0; i < ordinals.size(); i++) {
|
|
if (rx.cancelled_ || ri->AccumulationFull())
|
|
break;
|
|
if (!spot_cache.contains(ordinals[i]))
|
|
PrefetchSpots({ordinals.begin() + i,
|
|
ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)});
|
|
ri->ProcessImage(ordinals[i], spot_cache.at(ordinals[i]), RotAngle(ordinals[i]));
|
|
}
|
|
return ri;
|
|
}
|
|
|
|
// Feed every scheme (single-threaded spot cache), ready to be indexed.
|
|
auto Rugnux::FirstPassRun::FeedSchemes(IndexerThreadPool &pool, const DiffractionExperiment *x) -> SchemeIndexers {
|
|
SchemeIndexers ris;
|
|
ris.reserve(schemes.size());
|
|
for (const auto &[name, ordinals] : schemes) {
|
|
if (rx.cancelled_)
|
|
break;
|
|
ris.push_back(FeedScheme(pool, ordinals, x));
|
|
}
|
|
return ris;
|
|
}
|
|
|
|
// Run RunIndexing() of every fed scheme in the list at once. Each indexer reads only what it
|
|
// was fed and its own experiment, so the answers do not depend on what runs beside them.
|
|
void Rugnux::FirstPassRun::RunIndexing(const std::vector<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();
|
|
}
|
|
|
|
void Rugnux::FirstPassRun::RunIndexing(const SchemeIndexers &ris) {
|
|
std::vector<RotationIndexer *> rp;
|
|
for (const auto &ri : ris)
|
|
rp.push_back(ri.get());
|
|
RunIndexing(rp);
|
|
}
|
|
|
|
// Keep the lattice of the indexed schemes that indexes the most validation frames on the real
|
|
// per-image path. Scored against the spot lists and settings the run stands on now, which must
|
|
// be those the schemes were fed under.
|
|
// `named` names the indexers (`schemes` unless given) and `bp` is the standing pass they are
|
|
// judged against (none unless given).
|
|
auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefine &idx, const Schemes *named,
|
|
FirstPass bp) -> FirstPass {
|
|
if (!named)
|
|
named = &schemes;
|
|
for (size_t i = 0; i < ris.size(); i++) {
|
|
if (rx.cancelled_)
|
|
break;
|
|
if (!bp.indexer_error)
|
|
bp.indexer_error = ris[i]->GetIndexerError();
|
|
auto found = ris[i]->GetLattice();
|
|
if (!found.has_value())
|
|
continue;
|
|
RotationIndexerResult result = std::move(*found);
|
|
int score = CountIndexed(idx, result);
|
|
const std::string &name = (*named)[i].first;
|
|
logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score,
|
|
static_cast<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 = CountIndexed(idx, alt);
|
|
if (alt_score > majority) {
|
|
logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained "
|
|
"cell indexes {}/{} - the metric symmetry is a false promotion, dropping it",
|
|
name, result.search_result.centering,
|
|
gemmi::crystal_system_str(result.search_result.system), score,
|
|
static_cast<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 (rx.lattice_angles_demoted_ && result.unconstrained
|
|
&& (result.search_result.system == gemmi::CrystalSystem::Monoclinic
|
|
|| result.search_result.system == gemmi::CrystalSystem::Orthorhombic)) {
|
|
RotationIndexerResult alt = *result.unconstrained;
|
|
result = std::move(alt);
|
|
score = CountIndexed(idx, result);
|
|
}
|
|
|
|
// 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 std::vector<Coord> cloud = ValidationCloud(larger);
|
|
const auto hkl = IndexedMillerIndices(larger.lattice, cloud, index_tol);
|
|
const size_t n_smaller = IndexedMillerIndices(smaller.lattice,
|
|
ValidationCloud(smaller),
|
|
index_tol).size();
|
|
// The added class on the PRIMITIVE basis: a centred basis has absences of its
|
|
// own, which would read as an empty class. Both cells are counted on their
|
|
// Niggli-reduced bases, because a tolerance on fractional Miller indices depends
|
|
// on the basis, and a skewed one makes the sub-lattice miss spots it contains.
|
|
const CrystalLattice larger_p = larger.lattice.ToPrimitive(larger.search_result.centering)
|
|
.NiggliReduce();
|
|
const auto hkl_p = IndexedMillerIndices(larger_p, cloud, index_tol);
|
|
const auto ev = AddedClassOccupancy(hkl_p, static_cast<int>(nearest));
|
|
const CrystalLattice sub = SubLattice(larger_p, ev, static_cast<int>(nearest)).NiggliReduce();
|
|
const size_t n_sub = IndexedMillerIndices(sub, cloud, index_tol).size();
|
|
const bool larger_is_real = hkl.size() > n_smaller && hkl_p.size() > n_sub;
|
|
logger.Info("Axis-harmonic arbiter: '{}' ({:.0f} A^3) vs '{}' ({:.0f} A^3), {:.2f}x - "
|
|
"the larger cell accounts for {} validation spots against the smaller "
|
|
"cell's {}; on its own geometry {} against its ({},{},{}) sub-lattice's "
|
|
"{}, {:.1f}% of them off it, so the class it adds is {}",
|
|
name, vol, bp.name, bp.vol, ratio, hkl.size(), n_smaller,
|
|
hkl_p.size(), ev.u, ev.v, ev.w, n_sub, 100.0 * ev.occupancy,
|
|
larger_is_real ? "occupied enough to pay for the longer axis - the "
|
|
"larger cell is real"
|
|
: "not occupied enough to pay for the longer axis - the "
|
|
"larger cell is a multiple");
|
|
// Adopt the candidate exactly when its role matches the verdict.
|
|
integer_subcell = (cand_is_larger == larger_is_real);
|
|
// Only a win for the larger cell settles the pair. Declining it is a fall back
|
|
// to the default, and the warning below - which names the beam centre - is then
|
|
// still the most useful thing the run can say. It is the right diagnosis in
|
|
// particular for the failure this arbiter cannot see: a centre error along the
|
|
// spindle translates the derotated cloud rigidly by HALF a lattice spacing, and
|
|
// a half-offset lattice is indexable only on a doubled axis. Such a cell fits no
|
|
// index-n sublattice at all, so its added class reads the chance value (n-1)/n -
|
|
// the largest the occupancy can be - and it accounts for far fewer spots than
|
|
// the true cell, which is what this test sees it by.
|
|
harmonic_settled = larger_is_real;
|
|
}
|
|
// sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by
|
|
// sqrt(3) is the other harmonic seen on real data and its volume ratio is not a
|
|
// whole number. This is a REPORT and not a decision - it is deliberately taken
|
|
// over the whole pair, `clearly_more` included, because the harmonic is often the
|
|
// one that indexes more frames (see below).
|
|
if (integer_multiple || std::abs(ratio - std::sqrt(3.0)) < 0.15)
|
|
harmonic_ratio = ratio;
|
|
}
|
|
|
|
// Two candidates whose primitive volumes differ by a small integer or by sqrt(3) are
|
|
// one lattice and its axis harmonic, and what decides between them is the BEAM CENTRE.
|
|
// A centre error along the spindle translates the derotated cloud rigidly, which the
|
|
// FFT amplitude cannot see at all, so the peaks stay sharp and the indexer commits -
|
|
// with confidence - to a sub-multiple. Measured on real data: a 0.12 px change of
|
|
// centre, 0.03x of what the J0 law asks across the spindle, is the whole difference
|
|
// between the deposited cell and a halved axis, and the halved one is the answer that
|
|
// indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate
|
|
// this pair, and until this line nothing in the run said what the pair meant.
|
|
// ... unless the arbiter above established the larger cell on the data, which is
|
|
// exactly the question this warning says the run cannot answer.
|
|
if (harmonic_ratio > 0.0 && !harmonic_settled)
|
|
logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) "
|
|
"differ by {:.2f}x - one is an axis harmonic of the other. Which of "
|
|
"the two is right is settled by the beam centre to a fraction of a "
|
|
"pixel, and the harmonic is often the one that indexes MORE frames, "
|
|
"so the frame counts above cannot decide it. Check the beam centre",
|
|
name, vol, bp.name, bp.vol, harmonic_ratio);
|
|
|
|
if (!bp.result.has_value() || clearly_more || integer_subcell) {
|
|
bp.larger_settled = harmonic_settled && vol > bp.vol;
|
|
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);
|
|
} else if (harmonic_settled) {
|
|
bp.larger_settled = true;
|
|
}
|
|
}
|
|
return bp;
|
|
}
|
|
|
|
// Feed the schemes, index them concurrently, and keep the best (ScoreSchemes).
|
|
// alongside: started once the schemes are fed and before they are indexed, to run beside them.
|
|
auto Rugnux::FirstPassRun::PickBest(IndexerThreadPool &pool, IndexAndRefine &idx,
|
|
const std::function<void()> &alongside) -> FirstPass {
|
|
const SchemeIndexers ris = FeedSchemes(pool);
|
|
if (alongside)
|
|
alongside();
|
|
RunIndexing(ris);
|
|
return ScoreSchemes(ris, idx);
|
|
}
|
|
|
|
void Rugnux::FirstPassRun::ApplyState(const FirstPassState &st) {
|
|
const std::vector<float> r = st.rings ? rx.powder_.rings_q_recipA : std::vector<float>{};
|
|
rx.config_.spot_finding.measured_ring_q_recipA = r;
|
|
validation_settings.measured_ring_q_recipA = r;
|
|
rx.config_.spot_finding.high_resolution_limit = st.d_min;
|
|
validation_settings.high_resolution_limit = st.d_min;
|
|
first_pass_spots_per_image = st.seed;
|
|
spot_cache.clear(); // the per-image spot budget is spent again under these settings
|
|
}
|
|
|
|
// Powder-ring rescue - and, with it, the lean seed the rotation first pass never had.
|
|
//
|
|
// The pre-scan has already measured this sample's powder rings (see PreScan): a crystalline
|
|
// phase other than the crystal, whose reflections the indexer is otherwise handed as if they
|
|
// were this crystal's. Hexagonal ice is the only such phase whose rings can be named in
|
|
// advance, and the spot finder flags it from that fixed list; everything else - a shower of
|
|
// microcrystals, a salt out of the cryoprotectant, a phase nobody expected - has to be
|
|
// measured, and until it is, every one of its spots goes into the FFT. Measured on one such
|
|
// crystal: 876 spots a frame against the 80 a clean crystal on the same beamline gives, 70 %
|
|
// of them on rings, no scheme indexing a single validation frame, and the FFT taking a 10-30 A
|
|
// cell out of the shells.
|
|
//
|
|
// Setting the rings aside is half of it. The other half is HOW DEEP into each frame's spot
|
|
// list the first pass reads. It has always read 200 spots an image; the per-frame indexer, on
|
|
// the same problem, escalates 30 -> 80 -> the whole list and keeps the seed a lattice explains
|
|
// the largest fraction of (IndexAndRefine::DetermineLatticeAndSymmetry), because on a flooded
|
|
// frame a lean high-quality seed finds the lattice where the full list derails the search. The
|
|
// rotation pass gets the same escalation here. It is the general statement of what a
|
|
// resolution cut would do on this one crystal - the contaminant and the noise sit at the deep
|
|
// end of the intensity-ordered list wherever in resolution they happen to lie - and it needs
|
|
// nothing measured about where they are.
|
|
//
|
|
// It runs LAST: after every other rescue and after the pass-1 fallback above, on a run that
|
|
// would otherwise stop with no lattice at all. Measured over 100 datasets: placed earlier it
|
|
// pre-empted the rescues that already work - a crystal the beam-centre ladder fixes at 59/60
|
|
// and one the pass-1 fallback fixes at 60/60 both came out WORSE, because this ladder had
|
|
// already replaced the standing pass those rescues start from. A last resort has to be last.
|
|
//
|
|
// The two compose: FilterSpotsByCount has already ordered each image's spots non-ring first
|
|
// and strongest first, so a lean cap takes the crystal's strongest reflections and nothing
|
|
// else. Decided late, the way the axis sign above is: index each rung, keep the one that
|
|
// indexes the most validation frames. A run that already indexes never reaches here, and one
|
|
// that does cannot lose by trying - where the rings are the crystal's own rows, or where the
|
|
// crystal needs its deep spots to index at all, the pass that already ran stands.
|
|
//
|
|
// Called from two places, because the two questions it answers are coupled: a centre that is
|
|
// wrong and a spot list that is too deep each hide the other, exactly as the centre and the
|
|
// axis sign do above. `score_before` is what the caller's standing pass indexes and `min_gain`
|
|
// is what a rung must beat it by; nothing is changed and the settings are put back when no
|
|
// rung does.
|
|
//
|
|
// min_score is an ABSOLUTE bar beside the gain. A caller that takes a rung only when it
|
|
// indexes a majority must not be handed a state the ladder committed on a lesser rung:
|
|
// measured, the ladder adopted at 29/60, the caller rejected it against a majority of 30, and
|
|
// nothing put the settings back - the beam-centre walk that followed then scored every
|
|
// hypothesis on the leaked lean spot lists.
|
|
auto Rugnux::FirstPassRun::FirstPassLadder(int score_before, int min_gain, int min_score)
|
|
-> std::optional<FirstPass> {
|
|
if (rx.cancelled_)
|
|
return std::nullopt;
|
|
|
|
// The seeds, leanest first. The same ladder the per-frame indexer uses, one rung deeper
|
|
// because a rotation cloud pools many frames and can afford more.
|
|
// The last rung is EVERY spot of each frame: SIZE_MAX, not 0 - RotationIndexer reads 0 as
|
|
// "use my default", which is 200, so a 0 here made the deepest rung a duplicate of the
|
|
// 200 one and the log line "the whole spot list" a lie.
|
|
constexpr std::array<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
|
|
// rx.config_ here would read back the previous rung's floor once a pass has adopted one, and
|
|
// the restore path would then "put back" that instead of what the file asked for.
|
|
const auto file_limit = first_pass_file_d_min;
|
|
std::vector<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 : rx.prescan_spot_quantile_d_A_)
|
|
if (!file_limit || d > *file_limit) floors.emplace_back(d);
|
|
// Floors outer, seeds inner, so rungs that share the finder-facing settings (rings,
|
|
// d_min) sit adjacent: the seed is an indexer-side cap over the SAME found spots
|
|
// (MaxSpotsPerImage over the cached list), so the loop below spot-finds once per floor
|
|
// instead of once per rung. Adoption is decided over all rungs at once, so the visit
|
|
// order carries no meaning.
|
|
for (const auto &d_min : floors)
|
|
for (const size_t seed : SEEDS) {
|
|
if (!rx.powder_.rings_q_recipA.empty())
|
|
rungs.push_back({true, seed, d_min});
|
|
if (rx.powder_.rings_q_recipA.empty() && seed != SEEDS.back())
|
|
rungs.push_back({false, seed, d_min});
|
|
}
|
|
|
|
// A new floor (or ring handling) changes what the finder returns, so the spot lists are
|
|
// spent and found again; a new seed alone does not - it caps how much of the SAME lists
|
|
// the indexer reads - and re-finding 60 validation frames per rung was most of what a
|
|
// hopeless run's minutes were spent on.
|
|
const auto enter_rung = [&](size_t i) {
|
|
if (i == 0 || rungs[i].rings != rungs[i - 1].rings || rungs[i].d_min != rungs[i - 1].d_min)
|
|
ApplyState({rungs[i].rings, rungs[i].seed, rungs[i].d_min});
|
|
else
|
|
first_pass_spots_per_image = rungs[i].seed;
|
|
};
|
|
// Every rung is fed first and all of them are indexed at once: the rungs do not depend on
|
|
// one another, and indexing them one after the other left most of the machine idle for
|
|
// most of a hopeless run. They are then scored in order, each under its own settings again.
|
|
std::vector<SchemeIndexers> rung_ris(rungs.size());
|
|
std::vector<RotationIndexer *> all_ris;
|
|
for (size_t i = 0; i < rungs.size() && !rx.cancelled_; i++) {
|
|
enter_rung(i);
|
|
rung_ris[i] = FeedSchemes(*indexer_pool);
|
|
for (const auto &ri : rung_ris[i])
|
|
all_ris.push_back(ri.get());
|
|
}
|
|
RunIndexing(all_ris);
|
|
|
|
// Every rung's score, so the choice below can be made on all of them at once.
|
|
std::vector<int> rung_score(rungs.size(), -1);
|
|
std::vector<FirstPass> rung_pass(rungs.size());
|
|
int won = -1;
|
|
for (size_t i = 0; i < rungs.size() && !rx.cancelled_; i++) {
|
|
enter_rung(i);
|
|
FirstPass alt = ScoreSchemes(rung_ris[i], *indexer);
|
|
logger.Info(" first-pass rung: {} spots/image, rings {}, stop at {} -> {}/{} frames{}",
|
|
rungs[i].seed == std::numeric_limits<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) {
|
|
ApplyState(committed);
|
|
PrefetchSpots(validation);
|
|
const double won_vol = rung_pass[won].vol;
|
|
const size_t won_spots = IndexedMillerIndices(rung_pass[won].result->lattice,
|
|
ValidationCloud(*rung_pass[won].result),
|
|
index_tol).size();
|
|
for (size_t i = 0; i < rungs.size() && won >= 0; i++) {
|
|
// Only a rung the caller would have taken can demote one: the winner is chosen
|
|
// among those, and a rung that indexes too little to be adopted is too little
|
|
// evidence to overturn one either.
|
|
if (static_cast<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,
|
|
ValidationCloud(*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) {
|
|
ApplyState({rungs[won].rings, rungs[won].seed, d});
|
|
const int n = CountIndexed(*indexer, *rung_pass[won].result);
|
|
logger.Info("First-pass depth: the winning rung's lattice indexes {}/{} validation "
|
|
"frames reading to {}", n, static_cast<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};
|
|
ApplyState(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", rx.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];
|
|
}
|
|
ApplyState(committed); // back to what this run stands on, not to the defaults
|
|
return std::nullopt;
|
|
}
|
|
|
|
// Find the spots of everything the first pass is certain to want in one go - the schemes'
|
|
// opening frames and the validation set. A batch of workers costs an analysis engine each,
|
|
// which is far more than a frame is to analyse, so one batch for the lot beats one per scheme.
|
|
void Rugnux::FirstPassRun::PrefetchCertainFrames() {
|
|
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());
|
|
PrefetchSpots(upfront);
|
|
}
|
|
|
|
IndexingSettings Rugnux::FirstPassRun::ShortAxisSettings(const DiffractionExperiment &x) const {
|
|
auto settings = x.GetIndexingSettings();
|
|
settings.FFT_MinUnitCell_A(IndexingSettings::fft_min_unit_cell_limit_A);
|
|
settings.IndexingThreads(2);
|
|
settings.RefineThreads(std::clamp(rx.config_.nthreads / 2, 1, MAX_REFINE_THREADS));
|
|
return settings;
|
|
}
|
|
|
|
// The beam-centre check's pass at the measured centre, started beside the file centre's (see
|
|
// centre_spec).
|
|
void Rugnux::FirstPassRun::StartCentreSpeculation() {
|
|
if (get_gpu_count() < 2)
|
|
return;
|
|
rx.JoinBeamCenterCapture();
|
|
if (rx.cancelled_ || !rx.config_.beam_center_check || !rx.background_center_)
|
|
return;
|
|
auto &cs = centre_spec.emplace();
|
|
cs.experiment = std::make_unique<DiffractionExperiment>(rx.experiment_);
|
|
cs.experiment->BeamX_pxl(rx.background_center_->beam_x_pxl).BeamY_pxl(rx.background_center_->beam_y_pxl);
|
|
cs.mask = std::make_unique<PixelMask>(rx.pixel_mask_);
|
|
cs.mapping = std::make_unique<AzimuthalIntegrationMapping>(*cs.experiment, *cs.mask, *rx.azint_geometry_);
|
|
cs.experiment_key = rx.ExperimentKey(*cs.experiment);
|
|
cs.key = rx.SpotFindingKey(*cs.experiment, *cs.mapping);
|
|
cs.spots_per_image = first_pass_spots_per_image;
|
|
{
|
|
// What an earlier pass already found at this centre under these settings.
|
|
std::lock_guard lock(rx.first_pass_spots_->m);
|
|
for (const auto &[k, spots] : rx.first_pass_spots_->entries)
|
|
if (k == cs.key)
|
|
cs.spots = spots;
|
|
}
|
|
// At this pass's angles, as PrefetchSpots takes them.
|
|
const DiffractionGeometry spot_geometry = cs.experiment->GetDiffractionGeometry();
|
|
for (auto &[ordinal, spots] : cs.spots)
|
|
StampSpotAngles(spots, spot_geometry);
|
|
cs.done = std::async(std::launch::async, [this, &cs, settings = rx.config_.spot_finding, dataset = dataset,
|
|
start_image = start_image, nworkers = engines.size()] {
|
|
SpotFindingWorkers workers(rx.reader_, *cs.experiment, *cs.mapping, *cs.mask, settings, *dataset,
|
|
start_image, rx.config_.stride, nworkers);
|
|
const auto find = [&](const std::vector<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 FeedScheme, 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(rx.config_.verify_first_pass_memo);
|
|
for (size_t i = 0; i < ordinals.size(); i++) {
|
|
if (rx.cancelled_ || ri->AccumulationFull())
|
|
break;
|
|
if (!cs.spots.contains(ordinals[i]))
|
|
find({ordinals.begin() + i,
|
|
ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)});
|
|
std::optional<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);
|
|
RunIndexing(cs.ris);
|
|
});
|
|
}
|
|
|
|
// The short-axis pass, started beside the standard indexing (see short_spec).
|
|
void Rugnux::FirstPassRun::StartShortAxisSpeculation() {
|
|
if (rx.cancelled_ || rx.experiment_.GetUnitCell().has_value())
|
|
return;
|
|
auto &sp = short_spec.emplace();
|
|
sp.key = rx.ExperimentKey();
|
|
sp.mapping = spot_mapping;
|
|
sp.spots_per_image = first_pass_spots_per_image;
|
|
sp.experiment = std::make_unique<DiffractionExperiment>(rx.experiment_);
|
|
sp.pool = std::make_unique<IndexerThreadPool>(ShortAxisSettings(*sp.experiment),
|
|
IndexerConstruction::OnFirstUse);
|
|
for (const auto &[name, ordinals] : schemes)
|
|
sp.sris.push_back(FeedScheme(*sp.pool, ordinals, sp.experiment.get()));
|
|
for (auto &ri : sp.sris) {
|
|
RotationIndexer *rp = ri.get();
|
|
sp.futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
|
|
}
|
|
}
|
|
|
|
auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass {
|
|
return PickBest(*indexer_pool, *indexer, [&] {
|
|
StartCentreSpeculation();
|
|
StartShortAxisSpeculation();
|
|
});
|
|
}
|