mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 68 Wiki Activity
Files
642cbc7271ea2cdbca9b4fb6c2834e3368a29350
Jungfraujoch/rugnux/Rugnux.cpp
T
jungfrauandClaude Opus 5 642cbc7271 Stop the rotation pre-pass once it has what the second pass needs, and say where time went 
Three changes to the rotation two-pass and what it reports.

THE PRE-PASS STOPS EARLY. The first pass exists to measure five things for the second:
the post-refined detector distance and beam, the goniometer rotation scale, the
frame-order-smoothed mosaicity, the space group, and its own indexing result for the
supercell guard. All five are settled by the time the group is adopted.

It then carried on and did a full production run anyway - a fourth scale/combine/merge
in the adopted group, the correction surfaces, the error model, the resolution cutoff,
the twinning analysis, the Wilson B, the statistics table - and wrote MTZ, CIF, HKL and
a report to <prefix>_01_*. The second pass remakes every bit of that at the refined
geometry a few seconds later, which is the result anybody reads. Measured over the
rotation test set, running it anyway was 10 % of the battery, and up to 21 % of a
single merge-dominated crystal. It now stops at the space-group decision and writes
nothing; the _01_* files go with it, since they existed to compare the two passes and
nothing in the pipeline or the test harness reads them.

One value had to move for that. Pass 2 asks whether pass 1 reached its space group by
PROMOTING the point group, and read that off result.twinning - a struct only filled
once the final merge has run. It is now recorded where it is decided.

EVERY RUN NOW REPORTS WHERE ITS TIME WENT. The phase marks that drove the GUI progress
label are now also timed, and the run ends with a table of wall time per phase summed
over both passes, with the mean number of cores each phase kept busy beside it (CPU
time over wall time). That second column is the one that matters: a phase with a large
share and one core is a single thread doing all the work, which is a different problem
from a phase that is simply large. It costs two clock reads per phase, and it means a
slow dataset on a machine with no profiler on it can still be diagnosed from its log.
The cores column needs process CPU time and is zero where the platform has no
getrusage; the wall column works everywhere.

AND THE PER-IMAGE COST LINE WAS WRONG. It divided the per-stage means, which are per
worker, by the thread count - but the image loop is capped at four workers per GPU, so
with 2 GPUs at -N 48 those are 8 and 48 and every stage was reported 6x too small. It
divides by the loop's own worker count now, and names it, so the line says what it
means. This was the one output anybody tuning the program would read first.

Battery 6m52s -> 6m17s, space group 21/24 with the same three disagreements, no
failures, and no merged result outside the spread the binary already has against itself
between two runs. (One weak, 82%-complete crystal moves several points of R_meas
between any two runs of any binary, this one included - measured at 39.7 / 31.7 / 35.5 /
35.8 / 40.1 across five battery runs of four different builds.)

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2026-08-18 01:12:21 -04:00

2849 lines
177 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#ifdef __linux__
#include <sys/resource.h>
#endif
#include "Rugnux.h"
#include "ModelValidation.h"
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <condition_variable>
#include <deque>
#include <functional>
#include <future>
#include <mutex>
#include <numeric>
#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/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/geom_refinement/BeamCenterFromBackground.h"
#include "../image_analysis/geom_refinement/BeamCenterFromSpots.h"
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/indexing/IndexerThreadPool.h"
#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
#include "../image_analysis/azint/AzIntEngineCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
#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/scale_merge/TwinningAnalysis.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>
namespace {
// How precisely the beam centre has to be placed before the estimate is used at all. An estimate
// that fails it is treated as no estimate: the spot symmetry falls through to the background and
// the background to the header. A fit's own scatter will not serve as the gate - the ways either
// estimator goes badly wrong all leave a small one - so the spot symmetry reports instead how
// far its answer moves when the search is started elsewhere, which does separate them. On the 38
// rotation regression crystals that stays below 0.35 px against errors of at most 0.74; over a
// ring of injected geometries it rises past 1.0 px on 96 % of the answers that are 2 px or more
// out, and on none of the 38 as they are.
constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f;
// Images the projection is built from when the beam centre is wanted but the beam-stop pre-pass
// is off; with it on, that pre-pass's image count is used and one projection serves both.
constexpr int BEAM_CENTER_PROJECTION_IMAGES = 60;
// Spots per pre-scan frame kept for the beam centre, the strongest first.
constexpr size_t BEAM_CENTER_SPOTS_PER_IMAGE = 500;
// Frames read for the spot symmetry alone when the first sample could not measure it. Both
// estimators are counting statistics over PAIRS of frames - the Friedel one over frames half a
// turn apart, the second crossing over any two frames a reflection is recorded on - so a sweep
// whose frames carry few spots can only be measured by reading more of them. On the two
// rotation regression crystals that decline, this is what turns the decline into an answer
// 0.04 px from the post-refined centre; the sample the beam stop is built from is left alone.
constexpr int BEAM_CENTER_SPARSE_IMAGES = 400;
// Frames the spot symmetry reads, as a multiple of the beam stop's. Its two estimators are
// counts of matched PAIRS and the pairs grow faster than the frames do, so this is the one knob
// that buys independence from WHICH frames were drawn - and that dependence is the estimator's
// largest remaining failure. Asking the same code for one frame more or less re-draws the whole
// sample: over five draws of 38 sweeps, 60 frames leaves three answers more than 1 px out, none
// of which reports an unusual sigma, and 120 leaves none in 186. It is also what makes the end
// margin below free: at 60 frames excluding the ends moved three centres, at 120 it moves none.
constexpr int BEAM_CENTER_IMAGE_FACTOR = 2;
// Images left out at each end of the sweep. Shutter synchronisation disturbs the first and last
// frame or two of a run, and a partly exposed frame is dark over the whole detector - which is
// exactly the kind of error a projection cannot average away and a background comparison reads
// as signal. The margin is capped at a tenth of the sweep so a short run still has a sample.
constexpr int PRESCAN_END_MARGIN_IMAGES = 5;
// Workers reading the pre-scan sample. Each owns a shard of the beam-stop projection so no two
// threads touch the same accumulator, and a shard costs 20 bytes per pixel - 362 MB on a 16M
// detector - so this is capped well below the worker count of the run proper. The accumulation
// is memory-bound rather than compute-bound, so a handful of workers already saturates it.
constexpr size_t PRESCAN_MAX_WORKERS = 8;
// Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the
// first pass of two-pass rotation indexing.
std::vector<int> select_equally_spaced_image_ordinals(int images_to_process, int requested_images) {
std::vector<int> ret;
if (images_to_process <= 0 || requested_images <= 0)
return ret;
const int n = std::min(images_to_process, requested_images);
if (n == 1) {
ret.push_back(0);
return ret;
}
std::set<int> unique_ordinals;
for (int i = 0; i < n; i++)
unique_ordinals.insert(static_cast<int>(
std::llround(static_cast<double>(i) * static_cast<double>(images_to_process - 1) /
static_cast<double>(n - 1))));
ret.assign(unique_ordinals.begin(), unique_ordinals.end());
return ret;
}
// A sample arranged as pairs of images half a turn apart: half the ordinals spread over the
// part of the sweep that has a partner 180 deg later, and those partners. That is what the
// Friedel mirror needs. It is NOT the plain selection rearranged - even on a full turn the two
// sets are all but disjoint, sharing 2 of 60 ordinals on a 3600 image sweep - so the caller
// keeps them apart and gives each consumer the one it needs. Empty when the sweep is shorter
// than half a turn, which is the caller's signal that the spot symmetry cannot be measured here.
std::vector<int> select_half_turn_paired_ordinals(int images_to_process, int requested_images,
float degrees_per_ordinal) {
if (!(degrees_per_ordinal > 0.0f))
return {};
const int half_turn = static_cast<int>(std::lround(180.0 / degrees_per_ordinal));
if (half_turn < 1 || images_to_process - half_turn < 1)
return {};
std::set<int> ordinals;
for (const int ordinal : select_equally_spaced_image_ordinals(images_to_process - half_turn,
requested_images / 2)) {
ordinals.insert(ordinal);
ordinals.insert(ordinal + half_turn);
}
return {ordinals.begin(), ordinals.end()};
}
// Apply a goniometer rotation SCALE: the stage turned k times the angle the file records. What is
// wrong is the SWEEP, not where it began, so the per-frame increment and the per-frame oscillation
// width both scale by k while the starting angle is left alone - the stage reached that position
// before the sweep started. A constant offset in phi is in any case exactly degenerate with a
// rotation of the crystal orientation about the spindle, which indexing refines away, so anchoring
// the stretched sweep at image 0 costs nothing and keeps the first frame at the angle the file gives
// it. Every angle in the pipeline comes from this object, so scaling it is the whole correction.
GoniometerAxis ScaleRotation(const GoniometerAxis &g, float k) {
GoniometerAxis scaled(g.GetName(), g.GetStart_deg(), g.GetIncrement_deg() * k,
g.GetAxis(), g.GetHelicalStep());
if (const auto wedge = g.GetScreeningWedge())
scaled.ScreeningWedge(*wedge * k);
return scaled;
}
// Writing an image to the process file takes the global HDF5 mutex, which is the same mutex every
// worker thread needs to fetch its next image. The write itself is short, but with a worker per
// hardware thread they spend longer queueing for the lock than the writing takes: on a 16 Mpx
// rotation dataset the per-image loop runs 7.7 s with 8 workers and 13.1 s with 48. Give the file
// to one thread and let the workers post to it, and the lock has one taker again.
//
// A DataMessage does not own its pixels - it points into the reader's buffer - so the raw image is
// parked in the queue beside it. Without that the worker frees the pixels on its next iteration and
// the writer reads whatever landed there.
class ProcessFileWriter {
public:
ProcessFileWriter(FileWriter &writer, size_t capacity)
: writer_(writer), capacity_(capacity),
thread_([this] { Run(); }) {}
// Post one image. Blocks while the queue is full, which is what stops a run whose analysis
// outpaces its writer from holding every image it has ever processed in memory.
void Post(const DataMessage &msg, std::shared_ptr<JFJochReaderRawImage> img) {
std::unique_lock lock(m_);
space_.wait(lock, [this] { return queue_.size() < capacity_ || failed_; });
if (failed_)
return; // the error is rethrown by Finish(); dropping the rest is deliberate
queue_.push_back({msg, std::move(img)});
work_.notify_one();
}
// Drain, join, and rethrow whatever the writer thread hit. Must be called before the
// FileWriter is used for anything else - the end message, Finalize().
void Finish() {
{
std::lock_guard lock(m_);
done_ = true;
}
work_.notify_one();
if (thread_.joinable())
thread_.join();
if (error_)
std::rethrow_exception(error_);
}
~ProcessFileWriter() {
{
std::lock_guard lock(m_);
done_ = true;
failed_ = true; // an unwinding run should not wait for the backlog to be written
}
work_.notify_one();
space_.notify_all();
if (thread_.joinable())
thread_.join();
}
private:
struct Job {
DataMessage msg;
std::shared_ptr<JFJochReaderRawImage> img;
};
void Run() {
while (true) {
Job job;
{
std::unique_lock lock(m_);
work_.wait(lock, [this] { return !queue_.empty() || done_; });
if (queue_.empty())
return;
job = std::move(queue_.front());
queue_.pop_front();
}
space_.notify_one();
try {
writer_.Write(job.msg);
} catch (...) {
std::lock_guard lock(m_);
if (!error_)
error_ = std::current_exception();
failed_ = true;
space_.notify_all();
return;
}
}
}
FileWriter &writer_;
const size_t capacity_;
std::mutex m_;
std::condition_variable work_;
std::condition_variable space_;
std::deque<Job> queue_;
bool done_ = false;
bool failed_ = false;
std::exception_ptr error_;
std::thread thread_;
};
}
Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment,
PixelMask pixel_mask, ProcessConfig config)
: reader_(reader), experiment_(std::move(experiment)),
pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) {
// Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that
// already carries one starts clear; the user mask (bit 8) is left as it was loaded.
pixel_mask_.ClearBeamStopMask(experiment_);
// A manually asserted stage calibration applies to everything, before anything reads an angle.
if (config_.rotation_scale.has_value())
if (const auto g = experiment_.GetGoniometer()) {
Logger("Rugnux").Info("Goniometer rotation scale {:.5f} applied from the command line",
*config_.rotation_scale);
experiment_.Goniometer(ScaleRotation(*g, *config_.rotation_scale));
}
}
void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) {
Logger logger("Rugnux");
// Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find
// nothing (its pixels are masked by now) and so clear the mask the first pass established;
// the beam centre must not be re-read either, because by then the geometry is post-refined.
const auto &current = pixel_mask_.GetMask();
const bool have_shadow = std::any_of(current.begin(), current.end(),
[](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; });
const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow;
const bool want_beam_center = config_.estimate_beam_center && !beam_center_placed_;
if (!want_shadow && !want_beam_center)
return;
// The two consumers want different frames and each gets its own set. The shadow is built from
// the even spread over the sweep it has always been built from, so that asking for a beam
// centre cannot move the mask; the spot symmetry needs frames paired half a turn apart, where
// the sweep reaches that far. The two sets barely overlap, so the loop below walks their union
// and every frame is read once.
//
// Both samples are taken from the middle of the sweep, both ends left out: an even spread over
// [0, images_to_process) always picks the very first and the very last image, and those are the
// ones a shutter out of step spoils - the last carries under 90 % of the sample's median counts
// on 12 of 38 regression sweeps and is essentially unexposed on 5.
const auto goniometer = experiment_.GetGoniometer();
const int margin = std::min(PRESCAN_END_MARGIN_IMAGES, images_to_process / 10);
auto shadow_sample = select_equally_spaced_image_ordinals(images_to_process - 2 * margin,
frame_count);
auto spot_sample = want_beam_center && goniometer
? select_half_turn_paired_ordinals(
images_to_process - 2 * margin,
BEAM_CENTER_IMAGE_FACTOR * frame_count,
// Degrees per ORDINAL: the shifted goniometer's increment already
// is that, and rotation data cannot be strided anyway.
std::abs(goniometer->GetIncrement_deg()))
: std::vector<int>();
for (int &ordinal : shadow_sample)
ordinal += margin;
for (int &ordinal : spot_sample)
ordinal += margin;
const bool want_spot_symmetry = !spot_sample.empty();
const std::set<int> shadow_set(shadow_sample.begin(), shadow_sample.end());
const std::set<int> spot_set(spot_sample.begin(), spot_sample.end());
std::set<int> sample = shadow_set;
sample.insert(spot_set.begin(), spot_set.end());
if (sample.empty())
return;
// Say which of the two this is reading for. It runs before the first image of the run proper is
// processed, so without this the job sits at its initial status for the whole pre-scan - which is
// tens of seconds once the beam centre asks for its own frames on top of the projection.
MarkPhase(observer, want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop")
: "Beam centre");
ShadowFinder finder(experiment_, pixel_mask_);
// The spot symmetry is read in the same pass as the projection rather than in a pre-pass of its
// own. The image is preprocessed a second time for it: the shadow is accumulated from raw counts
// and the spot finder works on the converted image.
std::unique_ptr<AzimuthalIntegrationMapping> prescan_mapping;
if (want_spot_symmetry)
prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(experiment_, pixel_mask_);
// Everything reading one frame needs that cannot be shared between threads: its own scratch for
// each of the two decompressions, its own preprocessed image and its own spot finder.
struct PreScanWorker {
std::vector<uint8_t> shadow_buffer;
std::vector<uint8_t> decompression_buffer;
std::unique_ptr<ImagePreprocessorCPU> preprocessor;
std::unique_ptr<ImagePreprocessorBuffer> preprocessed;
std::unique_ptr<ImageSpotFinder> spot_finder;
};
const auto make_worker = [&] {
PreScanWorker w;
if (want_spot_symmetry) {
w.preprocessor = std::make_unique<ImagePreprocessorCPU>(experiment_, pixel_mask_);
w.preprocessed = std::make_unique<ImagePreprocessorBuffer>(experiment_.GetPixelsNum());
if (config_.spot_finding.adaptive_threshold)
w.spot_finder = std::make_unique<AdaptiveSpotFinderCPU>(*prescan_mapping);
else
w.spot_finder = std::make_unique<ImageSpotFinderCPU>(experiment_.GetXPixelsNumConv(),
experiment_.GetYPixelsNumConv());
}
return w;
};
// The spots of one image. Taken as a step of its own because the sparse-sweep pass below reads
// images the shadow does not. The frame they belong to is stamped by the caller, which is what
// keeps the pool independent of the order the workers finish in.
const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx,
std::vector<BeamCenterSpot> &out) {
try {
w.preprocessor->Analyze(*w.preprocessed,
image.GetUncompressedPtr(w.decompression_buffer), image.GetMode());
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to preprocess image {}: {}", image_idx, e.what());
return;
}
auto spots = w.spot_finder->Run(*w.preprocessed, config_.spot_finding);
// The strongest of a crowded frame: the symmetry is over-determined either way, and the
// matching is quadratic in the spots of one frame.
if (spots.size() > BEAM_CENTER_SPOTS_PER_IMAGE) {
std::partial_sort(spots.begin(), spots.begin() + BEAM_CENTER_SPOTS_PER_IMAGE, spots.end(),
[](const DiffractionSpot &a, const DiffractionSpot &b) {
return a.Count() > b.Count(); });
spots.resize(BEAM_CENTER_SPOTS_PER_IMAGE);
}
for (const auto &spot : spots) {
const Coord centroid = spot.RawCoord();
out.push_back({centroid.x, centroid.y, static_cast<float>(spot.Count()), 0});
}
};
std::vector<float> frame_angle_deg;
std::vector<BeamCenterSpot> beam_center_spots;
// Read the sample on several workers. The reader serialises on the HDF5 lock, but the
// decompression, the projection and the spot finding - which is all of the cost on a large
// detector - run in parallel. Each worker accumulates into a shard of its own, so nothing is
// locked while an image is added, and the per-frame results are stitched together in sample
// order below so the beam centre sees the same input however the workers interleaved.
{
const std::vector<int> ordinals(sample.begin(), sample.end());
std::vector<std::vector<BeamCenterSpot>> spots_of(ordinals.size());
// A frame joins the pool only if it could be read, as it did when this was a serial loop.
std::vector<char> spot_read(ordinals.size(), 0);
const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
std::min(PRESCAN_MAX_WORKERS, ordinals.size()));
finder.SetShardCount(nworkers);
std::atomic<size_t> next{0};
std::vector<std::future<void>> futures;
futures.reserve(nworkers);
for (size_t t = 0; t < nworkers; t++)
futures.emplace_back(std::async(std::launch::async, [&, t] {
PreScanWorker w = make_worker();
for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) {
const int ordinal = ordinals[i];
const int image_idx = start_image + ordinal * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (shadow_set.contains(ordinal))
finder.AddImage(msg, w.shadow_buffer, t);
if (spot_set.contains(ordinal)) {
spot_read[i] = 1;
find_spots(w, msg.image, image_idx, spots_of[i]);
}
}
}));
for (auto &f : futures) f.get();
// Frame numbering follows the sample order, exactly as the serial read did.
for (size_t i = 0; i < ordinals.size(); i++) {
if (!spot_read[i]) continue;
const int frame = static_cast<int>(frame_angle_deg.size());
// The goniometer has been shifted so that a LOCAL ordinal maps to the angle of the
// original image it came from (see RunPipeline), so the angle is asked for by ordinal,
// as everywhere else.
frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(ordinals[i])));
for (auto &spot : spots_of[i]) {
spot.frame = frame;
beam_center_spots.push_back(spot);
}
}
}
if (finder.GetFrameCount() == 0) {
logger.Warning("Pre-scan: no image could be read. Skipping.");
return;
}
if (want_shadow) {
const auto shadow = finder.GetMask(config_.nthreads);
const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u);
pixel_mask_.LoadBeamStopMask(experiment_, shadow);
logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",
shadowed, 100.0 * static_cast<double>(shadowed) / static_cast<double>(shadow.size()),
finder.GetFrameCount());
}
// The beam centre comes after, so that the shadow is in the mask by the time it is placed and
// the holder arm - the largest azimuthal asymmetry on the detector - is out of the way of the
// background the centre is read from.
if (!want_beam_center)
return;
beam_center_placed_ = true;
// The spot symmetry first where the sweep supports it: it is exact geometry rather than a
// property of the scattering, and it measures the sets the background cannot - a flat radial
// profile at long wavelength leaves that fit nothing to bite on. The background answers
// wherever the sweep does not reach half a turn, which is every still and every screening
// wedge, so the two cover between them what neither covers alone.
std::string source = "spot symmetry";
std::optional<BeamCenterEstimate> estimate;
SpindleEstimate spindle;
if (want_spot_symmetry)
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
// Where it does not come out - no answer at all, or one that moves when the search is started
// elsewhere - read more of the sweep and ask again. Both estimators are counts of matched spot
// PAIRS, the Friedel one over frames half a turn apart and the second crossing over any two
// frames one reflection is recorded on, so a sweep whose frames carry few spots, or whose pairs
// match rarely, is measurable only from more of them. More data first and a weaker estimator
// second: this runs before the fall-through, and only what survives it falls through. The extra
// frames go to the spot finder alone, so the beam-stop projection keeps the images it was
// validated on.
if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL)) {
// What has already been read FOR THE SPOTS, which is not the whole sample: the frames the
// shadow was built from carry none of them, and testing against the union would silently
// skip the ones this pass exists to read.
std::vector<int> extra;
for (int ordinal : select_half_turn_paired_ordinals(
images_to_process - 2 * margin, BEAM_CENTER_SPARSE_IMAGES,
std::abs(goniometer->GetIncrement_deg()))) {
ordinal += margin;
if (!spot_set.contains(ordinal))
extra.push_back(ordinal);
}
if (!extra.empty()) {
logger.Info("Beam centre: the spot symmetry does not come out on {} frames "
"({} spots); reading {} more frames",
frame_angle_deg.size(), beam_center_spots.size(), extra.size());
// Its own phase: this reads several times the frames the pass above did, and is the one
// step of the pre-scan a user is likely to sit through wondering what is happening.
MarkPhase(observer, "Beam centre (reading more frames)");
// Read on workers as above - this pass is several times the size of the first one - and
// append in `extra` order so the pool does not depend on how they interleaved.
std::vector<std::vector<BeamCenterSpot>> extra_spots(extra.size());
std::vector<char> extra_read(extra.size(), 0);
const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
std::min(PRESCAN_MAX_WORKERS, extra.size()));
std::atomic<size_t> next{0};
std::vector<std::future<void>> futures;
futures.reserve(nworkers);
for (size_t t = 0; t < nworkers; t++)
futures.emplace_back(std::async(std::launch::async, [&] {
PreScanWorker w = make_worker();
for (size_t i = next.fetch_add(1); i < extra.size(); i = next.fetch_add(1)) {
const int image_idx = start_image + extra[i] * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (img) {
extra_read[i] = 1;
find_spots(w, img->image, image_idx, extra_spots[i]);
}
}
}));
for (auto &f : futures) f.get();
for (size_t i = 0; i < extra.size(); i++) {
if (!extra_read[i]) continue;
const int frame = static_cast<int>(frame_angle_deg.size());
frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(extra[i])));
for (auto &spot : extra_spots[i]) {
spot.frame = frame;
beam_center_spots.push_back(spot);
}
}
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg,
beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
}
}
if (config_.fit_spindle && std::isfinite(spindle.azimuth_rad))
logger.Info("Spindle: azimuth about the beam {:.3f} mrad, tip towards it {:.3f} mrad "
"(Friedel vote {:.0f} -> {:.0f})",
1e3 * spindle.azimuth_rad, 1e3 * spindle.tip_rad,
spindle.vote_excess_nominal, spindle.vote_excess);
if (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL) {
source = "background";
estimate = FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection());
}
if (!estimate) {
logger.Info("Beam centre: not measurable, keeping ({:.2f},{:.2f})",
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
return;
}
const float moved = std::hypot(estimate->beam_x_pxl - experiment_.GetBeamX_pxl(),
estimate->beam_y_pxl - experiment_.GetBeamY_pxl());
const bool commit = estimate->sigma_pxl <= MAX_BEAM_CENTER_SIGMA_PXL;
logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, "
"sigma {:.2f} px => {}",
source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl,
commit ? "COMMIT" : "reject (kept header)");
if (commit)
experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl);
}
void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process,
RugnuxObserver *observer) {
Logger logger("Rugnux");
// Rotation has its own two-pass; this is stills only.
if (experiment_.IsRotationIndexing()) {
logger.Warning("--refine-geometry is stills-only; rotation data uses two-pass indexing. Ignoring.");
return;
}
// The bundle anchors a known cell; without one there is nothing to refine against.
const auto cell = experiment_.GetUnitCell();
if (!cell.has_value()) {
logger.Warning("--refine-geometry needs a reference unit cell (-C / -S / reference MTZ). Skipping.");
return;
}
const int refine_frames = std::max(1, config_.refine_geometry.value());
MarkPhase(observer, "Geometry refinement (first pass)");
const auto dataset = reader_.GetDataset();
// Index a spread sample large enough to yield ~refine_frames strong frames even at a low hit rate.
const int sample_budget = std::min(images_to_process, std::max(refine_frames * 50, 8000));
const std::vector<int> sample = select_equally_spaced_image_ordinals(images_to_process, sample_budget);
// First-pass engines, built from the current (nominal) geometry: an index-only pass (no
// integration) to obtain each frame's spots + assigned HKL + orientation.
AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);
IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse);
IndexAndRefine indexer(experiment_, &pool, /*retain_outcomes=*/false);
auto pass_settings = config_.spot_finding;
pass_settings.enable = true;
pass_settings.indexing = true;
pass_settings.quick_integration = false;
// A frame is worth bundling only if it indexed this many spots; the bundle then takes the
// refine_frames strongest of those.
constexpr int MIN_STRONG_SPOTS = 20;
// Stop sampling once there are comfortably more strong frames than the bundle will use, instead of
// draining the whole budget. The budget is sized for a low-hit-rate dataset (refine_frames * 50, at
// least 8000), so on data that indexes well it meant re-indexing the ENTIRE run to keep 200 frames -
// 99% of this pass was sampling. Four times the bundle leaves the "strongest N" selection a real
// pool to choose from.
const int strong_target = refine_frames * 4;
// The sample is cut into a FIXED number of interleaved stripes - stripe s is sample positions
// s, s+STRIPES, s+2*STRIPES, ... - and each stripe stops once it has contributed its share of the
// strong frames. Two things follow, both of which the previous shared-cursor-with-a-shared-stop
// arrangement got wrong. Every stripe is spread over the WHOLE run, so stopping early no longer
// means fitting the geometry to the beginning of it. And which frames get examined depends only
// on the data: it is not a race between the workers, and it does not change with -N, because a
// stripe is processed identically whichever worker happens to claim it.
constexpr int STRIPES = 32;
const int stripe_strong_target = std::max(1, (strong_target + STRIPES - 1) / STRIPES);
std::vector<GeomRefineFrame> frames;
std::mutex frames_mutex;
std::atomic<int> next_stripe = 0;
std::atomic<int> examined = 0;
auto worker = [&]() {
pin_gpu(); // round-robin per worker thread; must precede engine construction
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
for (int stripe = next_stripe.fetch_add(1); stripe < STRIPES && !cancelled_;
stripe = next_stripe.fetch_add(1)) {
int stripe_strong = 0;
for (int idx = stripe; idx < static_cast<int>(sample.size()) && !cancelled_
&& stripe_strong < stripe_strong_target; idx += STRIPES) {
examined.fetch_add(1, std::memory_order_relaxed);
const int ordinal = sample[idx];
const int image_idx = start_image + ordinal * config_.stride;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Warning("Geometry refinement: failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
try {
analysis.Analyze(msg, profile, pass_settings);
} catch (const std::exception &e) {
continue;
}
if (!msg.indexing_result.value_or(false) || !msg.indexing_lattice.has_value())
continue;
GeomRefineFrame f;
f.lattice = *msg.indexing_lattice;
f.ordinal = ordinal;
for (const auto &s : msg.spots)
if (s.indexed && s.lattice == 0)
f.spots.push_back(GeomRefineSpot{s.x, s.y,
static_cast<int32_t>(s.h), static_cast<int32_t>(s.k), static_cast<int32_t>(s.l)});
if (f.spots.size() >= 6) {
if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
stripe_strong++;
const std::unique_lock ul(frames_mutex);
frames.push_back(std::move(f));
}
}
}
};
std::vector<std::future<void>> futures;
futures.reserve(config_.nthreads);
for (int i = 0; i < config_.nthreads; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &fut : futures)
fut.get();
if (cancelled_)
return;
// Keep the strongest frames (most indexed spots) for the bundle - the true cell indexes many
// spots per frame, and orientation diversity comes for free from independent serial stills.
std::vector<GeomRefineFrame> strong;
for (auto &f : frames)
if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
strong.push_back(std::move(f));
// frames is in the order the workers happened to finish, so the spot count alone does not order
// it: without the ordinal tie-break, equally strong frames would swap places between runs and a
// different bundle would be refined.
std::sort(strong.begin(), strong.end(),
[](const GeomRefineFrame &a, const GeomRefineFrame &b) {
if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size();
return a.ordinal < b.ordinal;
});
if (static_cast<int>(strong.size()) > refine_frames)
strong.resize(refine_frames);
logger.Info("Geometry refinement: indexed {} of {} sampled frames (budget {}), bundling {} strong "
"frames (>= {} spots)",
frames.size(), std::min<int>(examined.load(), static_cast<int>(sample.size())),
sample.size(), strong.size(), MIN_STRONG_SPOTS);
GeometryRefinerSettings gr_settings;
gr_settings.crystal_system = experiment_.GetCrystalSystem();
gr_settings.num_threads = config_.nthreads;
const GeometryRefinerResult r = RefineGlobalGeometry(
experiment_.GetDiffractionGeometry(), *cell, strong, gr_settings);
if (!r.ok) {
logger.Warning("Geometry refinement did not run (too few strong frames / did not converge); "
"keeping the input geometry");
return;
}
logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm",
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px,
experiment_.GetDetectorDistance_mm(), r.distance_mm);
logger.Info("Geometry refinement: cell a/b/c {:.3f}/{:.3f}/{:.3f} -> {:.3f}/{:.3f}/{:.3f} A "
"(median residual {:.3f} px, {} frames, {} spots)",
cell->a, cell->b, cell->c, r.cell.a, r.cell.b, r.cell.c,
r.median_residual_px, r.frames_used, r.spots_used);
experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm);
experiment_.SetUnitCell(r.cell);
// The beam centre now comes from indexed spots, which the background cannot better.
beam_center_placed_ = true;
}
ProcessResult Rugnux::Run(RugnuxObserver *observer) {
// Rotation two-pass geometry post-refinement: the first pass integrates at the header geometry and
// post-refines the detector geometry (distance + beam from the observed spot positions, cell scale + axis
// from phi_obs); the second pass re-indexes and re-integrates with the refined geometry. BOTH passes are
// written out - the second (refined-geometry) pass is the CANONICAL result and takes the plain "<prefix>_*"
// name, while the first (header-geometry) pass is kept as "<prefix>_01_*" for comparison - so the two can be
// compared and the ice-ring flags / geometry stored in each master reflect the geometry that pass used.
if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry
&& experiment_.IsRotationIndexing()) {
Logger logger("Rugnux");
const auto gonio_snapshot = experiment_.GetGoniometer();
prepass_detector_geometry_.reset();
prepass_rotation_scale_.reset();
prepass_result_.reset();
force_rotation_result_.reset();
prepass_merge_sg_.reset();
prepass_promoted_point_group_ = false;
logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry), measuring "
"the geometry the second pass will integrate at");
auto pass1 = RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/true);
pass1.pass_number = 1;
pass1.pass_count = 2;
if (cancelled_) return pass1;
if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift
// Apply the post-refined detector geometry for the second pass, keeping the header geometry so
// the run can go back to it if the refined pass turns out worse (see the quality guard below).
const std::array<float, 3> header_geometry{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
experiment_.GetDetectorDistance_mm()};
if (prepass_detector_geometry_) {
const auto &g = *prepass_detector_geometry_;
experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]);
}
// ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles
// in the file are the commanded ones, so a stage that ran fast is a geometry error like any other.
if (prepass_rotation_scale_ && gonio_snapshot) {
experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_));
// The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has
// just stopped using, and the override can only ever raise the second pass's own estimate (it
// takes the larger of the two). Carrying it over would hold the second pass at the rocking
// width the uncorrected angles produced - the correction half-applied. Drop it and let the
// second pass fit its own.
prepass_mosaicity_.clear();
logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass "
"(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_,
gonio_snapshot->GetIncrement_deg(),
gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_);
}
// Space group: the second pass RE-INDEXES DE NOVO (clear the group here) so the indexer's pseudo-
// symmetry safeguards recover the true cell - reusing pass-1's group in the indexer forces build_sr
// onto a wrong / doubled cell (a huge oblique cell collapses; a pseudo-centred cell doubles). Pass-1's
// group is instead reinstated for the MERGE ONLY (RunPipeline, before RotationScaleMerge), so the
// symmetry stays fixed with no re-search / no flip.
//
// This holds for a group the centred-lattice test determined too, even though its conventional
// setting is not the one the de-novo primitive frame comes back in. That is precisely what the
// primitive->conventional reindex below does, and when the reindex declines - pass 2's metric is
// not the group's - the centring check right after it stops the pass being adopted at all. Letting
// such a group be re-searched instead put the second pass's determination back in play at a
// slightly different geometry, which on a lattice whose centring is pseudo-symmetric to a couple of
// tenths of a percent is a coin toss: measured on a C-centred monoclinic crystal, pass 1 confirmed
// the centring and pass 2 dropped it, and the run merged in P1.
experiment_.SpaceGroupNumber(std::nullopt);
prepass_merge_sg_ = pass1.space_group_number;
prepass_promoted_point_group_ = pass1.promoted_point_group;
logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
config_.output_prefix);
auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
// The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever
// result is returned - it is a result of the run, not of the pass that happened to fit it.
pass2.post_refine = pass1.post_refine;
pass2.pass_number = 2;
pass2.pass_count = 2;
// Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a
// periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell
// basin onto a spurious SUPERCELL (a large multiple of pass-1's cell volume), and the space group
// reused from pass-1 is then stamped onto that wrong-shaped cell - catastrophic. Pass 1 (header
// geometry, un-perturbed) found the correct lattice, so re-run the second pass with pass-1's FULL
// indexing result forced (its refined lattice + orientation + rotation axis + geometry - not just the
// cell, which loses the per-frame orientation and integrates worse): the re-run then matches de-novo
// quality but at the refined geometry. Only fires when pass 2 is markedly LARGER than pass 1 - a
// genuine de-novo demotion to a smaller primitive (the pseudo-symmetry case the second pass re-indexes
// de novo to catch) has pass2 <= pass1 and is left untouched.
if (!cancelled_ && prepass_result_ && pass1.consensus_cell && pass2.consensus_cell) {
const double v1 = gemmi::UnitCell(*pass1.consensus_cell).volume;
const double v2 = gemmi::UnitCell(*pass2.consensus_cell).volume;
if (v1 > 0.0 && v2 > 1.5 * v1) {
logger.Info("Two-pass: second-pass cell volume {:.0f} A^3 is {:.2f}x pass-1 ({:.0f} A^3) - a "
"spurious supercell; re-running the second pass with pass-1's result forced",
v2, v2 / v1, v1);
force_rotation_result_ = *prepass_result_;
// Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles
// from there - forcing it whole would put the uncorrected angles back and silently undo
// the rotation scale for the re-run.
if (prepass_rotation_scale_ && force_rotation_result_->axis)
force_rotation_result_->axis =
ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_);
pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
pass2.post_refine = pass1.post_refine;
pass2.pass_number = 3;
pass2.pass_count = 3;
pass2.pass_decision = fmt::format(
"post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo "
"cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious supercell",
v2, v2 / v1, v1);
force_rotation_result_.reset();
}
}
// The space group was DECIDED in pass 1 and only reused by pass 2, which therefore has no search
// of its own to report. Carry pass 1's over, or the evidence behind the decision (operator CCs,
// absence counts, the refused higher symmetry) is never shown on the default rotation path.
if (!pass2.space_group_search.has_value())
pass2.space_group_search = pass1.space_group_search;
// Pass 2 is normally the better answer, which is why it is the canonical output - but it is not
// guaranteed to be, and until now it was adopted whatever it produced. Two ways it can be wrong:
// it merges more unique reflections than the cell it settled on can hold (completeness above
// 100% is arithmetically impossible and means the cell is wrong), or its CC1/2 collapses
// relative to pass 1. Measured on a large-cell crystal: pass 1 195538 unique at 92.6% and
// CC1/2 0.98, pass 2 134667 at "117%" and CC1/2 0.62. Both bounds are set where only a failure
// reaches them, so a normal run - where pass 2 is a little better - keeps pass 2.
if (!cancelled_ && pass1.has_merge_statistics && pass2.has_merge_statistics) {
const auto completeness = [](const ProcessResult &r) {
const auto &o = r.merge_statistics.overall;
return o.possible_unique_reflections > 0
? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0;
};
const double compl2 = completeness(pass2);
const double cc1 = pass1.merge_statistics.overall.cc_half;
const double cc2 = pass2.merge_statistics.overall.cc_half;
constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more
constexpr double MAX_CC_HALF_LOSS = 0.05;
if (pass2.lattice_conflicts_with_prepass_sg
|| compl2 > MAX_CREDIBLE_COMPLETENESS || cc2 < cc1 - MAX_CC_HALF_LOSS) {
logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass "
"(completeness {:.1f}% vs {:.1f}%, CC1/2 {:.3f} vs {:.3f}) - going back to "
"the header geometry. The refined geometry did not help this crystal.",
compl2, completeness(pass1), cc2, cc1);
// The refined pass has already written the canonical files, so re-run at the header
// geometry to replace them - the same remedy the supercell collapse above uses, and it
// only costs a pass on a crystal that was going to be wrong otherwise.
experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1])
.DetectorDistance_mm(header_geometry[2]);
if (prepass_rotation_scale_ && gonio_snapshot)
experiment_.Goniometer(*gonio_snapshot); // and the header rotation angles with it
experiment_.SpaceGroupNumber(std::nullopt);
auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
if (!redo.space_group_search.has_value())
redo.space_group_search = pass1.space_group_search;
redo.post_refine = pass1.post_refine;
redo.pass_number = pass2.pass_count + 1;
redo.pass_count = pass2.pass_count + 1;
redo.pass_decision = fmt::format(
"header geometry re-adopted: the post-refined pass was worse "
"(completeness {:.1f}% vs {:.1f}%, CC1/2 {:.3f} vs {:.3f})",
compl2, completeness(pass1), cc2, cc1);
return redo;
}
if (pass2.pass_decision.empty())
pass2.pass_decision = fmt::format(
"post-refined geometry adopted (completeness {:.1f}% vs {:.1f}%, "
"CC1/2 {:.3f} vs {:.3f})", compl2, completeness(pass1), cc2, cc1);
}
if (pass2.pass_decision.empty())
pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_)
? "post-refined geometry adopted"
: "the post-refinement committed no geometry change, so this pass reproduces the first";
return pass2;
}
return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
}
namespace {
// Angle between the crystal's symmetry axes and the spindle. A rotation sweep never records the
// reflections whose reciprocal vector lies close to the spindle - the blind cusp - and normally
// symmetry fills it in from an equivalent elsewhere in reciprocal space. It cannot when a symmetry
// axis IS the spindle: the cusp then maps onto itself, every reflection in it is equivalent only to
// other reflections in it, and the hole stays empty however long the sweep runs. Returns the
// smallest angle in degrees between any proper rotation axis of the space group and the goniometer
// axis, with the order of the axis that achieves it.
std::optional<std::pair<double, int>> ClosestSymmetryAxisToSpindle(const gemmi::SpaceGroup &sg,
const CrystalLattice &lattice,
const Coord &spindle) {
const double spindle_len = spindle.Length();
if (spindle_len < 1e-9)
return std::nullopt;
std::optional<std::pair<double, int>> best;
for (const auto &op : sg.operations().derive_symmorphic().sym_ops) {
if (op.rot == gemmi::Op::identity().rot)
continue;
// Order of the rotation, then project onto its invariant direction by summing its powers:
// (1/n) sum_k W^k annihilates everything except the axis.
const gemmi::Op w{op.rot, {0, 0, 0}, op.notation};
int order = 0;
gemmi::Op cur = gemmi::Op::identity();
std::array<std::array<double, 3>, 3> sum{};
for (int k = 0; k < 6; ++k) {
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
sum[i][j] += static_cast<double>(cur.rot[i][j]) / gemmi::Op::DEN;
cur = cur.combine(w).wrap();
++order;
if (cur.rot == gemmi::Op::identity().rot)
break;
}
if (order < 2)
continue;
// Any non-degenerate column of the projector is the axis in fractional direct coordinates.
for (int c = 0; c < 3; ++c) {
const Coord axis = lattice.Vec0() * static_cast<float>(sum[0][c])
+ lattice.Vec1() * static_cast<float>(sum[1][c])
+ lattice.Vec2() * static_cast<float>(sum[2][c]);
const double len = axis.Length();
if (len < 1e-6 * lattice.Vec0().Length())
continue;
const double cosang = std::abs((axis * spindle) / (len * spindle_len));
const double ang = std::acos(std::min(1.0, cosang)) * 180.0 / PI;
if (!best.has_value() || ang < best->first)
best = std::make_pair(ang, order);
break;
}
}
return best;
}
}
namespace {
// CPU time this process has used across all its threads, in seconds. Compared against wall time it
// says how many cores a phase kept busy, which is the difference between a phase that is slow
// because it has work to do and one that is slow because it is running on a single thread.
double process_cpu_seconds() {
#ifdef __linux__
rusage ru{};
if (getrusage(RUSAGE_SELF, &ru) != 0)
return 0.0;
const auto secs = [](const timeval &t) { return t.tv_sec + t.tv_usec * 1e-6; };
return secs(ru.ru_utime) + secs(ru.ru_stime);
#else
return 0.0;
#endif
}
}
void Rugnux::ClosePhase() {
if (phase_current_.empty())
return;
const double wall = std::chrono::duration<double>(std::chrono::steady_clock::now() - phase_start_).count();
const double cpu = process_cpu_seconds() - phase_start_cpu_s_;
// A phase entered again - the two passes of a rotation run enter most of them twice - adds to the
// entry it already has, so the report has one line per phase rather than one per visit.
auto it = std::find_if(phase_timings_.begin(), phase_timings_.end(),
[this](const PhaseTiming &p) { return p.name == phase_current_; });
if (it == phase_timings_.end())
phase_timings_.push_back({phase_current_, wall, cpu});
else {
it->wall_s += wall;
it->cpu_s += cpu;
}
phase_current_.clear();
}
void Rugnux::MarkPhase(RugnuxObserver *observer, const std::string &phase) {
ClosePhase();
phase_current_ = phase;
phase_start_ = std::chrono::steady_clock::now();
phase_start_cpu_s_ = process_cpu_seconds();
if (observer)
observer->OnPhase(phase);
}
ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) {
Logger logger("Rugnux");
ProcessResult result;
const auto dataset = reader_.GetDataset();
if (!dataset)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No experiment dataset found in the input file");
if (config_.stride <= 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive");
const auto total_images_in_file = static_cast<int>(reader_.GetNumberOfImages());
int end_image = config_.end_image;
if (end_image < 0 || end_image > total_images_in_file)
end_image = total_images_in_file;
const int start_image = config_.start_image;
// Round UP: the worker loop below walks ordinals until start + ordinal * stride reaches end_image,
// so a stride that does not divide the range evenly still processes a final partial step. Rounding
// down here would size the per-image arrays one short of the ordinals that index them.
const int images_to_process = (end_image - start_image + config_.stride - 1) / config_.stride;
if (images_to_process <= 0) {
logger.Warning("No images to process (start {}, end {}, stride {}, total {})",
start_image, end_image, config_.stride, total_images_in_file);
return result;
}
const bool full = (config_.mode == ProcessMode::FullAnalysis);
// Powder calibration by spots finds spots and nothing else, so it runs the same per-image engine as
// the full analysis (spot finding lives inside MXAnalysisWithoutFPGA) with indexing switched off in
// its settings; by rings it only needs the azimuthal profile, which is the azint worker exactly.
const bool calibration = (config_.mode == ProcessMode::Calibration);
const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots;
const bool per_image_analysis = full || calibration_spots;
const bool write_files = write_output && !config_.output_prefix.empty();
// Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization,
// space group, unit cell, ...) are configured on experiment_ by the caller.
experiment_.BitDepthImage(32).PixelSigned(true);
experiment_.Mode(DetectorMode::Standard);
experiment_.OverwriteExistingFiles(true);
// Offline processing: the output prefix is a local path the operator chose (rugnux -o or the
// viewer's "next to the input file"), so use the trusted setter that permits an absolute path.
experiment_.FilePrefixTrusted(config_.output_prefix.empty() ? "output" : config_.output_prefix);
experiment_.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
experiment_.ImagesPerTrigger(images_to_process);
experiment_.NumTriggers(1);
if (full)
experiment_.Compression(CompressionAlgorithm::BSHUF_LZ4);
// The pipeline indexes images 0..N-1 within this run; if we process a sub-range/strided
// selection, shift the goniometer so local index i maps to the angle of original image
// start+i*stride (keeping the per-image rotation wedge), otherwise rotation angles would be
// wrong for any start_image != 0.
if (const auto g = experiment_.GetGoniometer();
g.has_value() && (start_image != 0 || config_.stride != 1)) {
const float incr = g->GetIncrement_deg();
GoniometerAxis shifted(g->GetName(),
g->GetStart_deg() + incr * static_cast<float>(start_image),
incr * static_cast<float>(config_.stride),
g->GetAxis(), g->GetHelicalStep());
shifted.ScreeningWedge(g->GetScreeningWedge().value_or(incr));
experiment_.Goniometer(shifted);
}
// Stills global geometry-refinement first pass: bundle-adjust the shared beam/distance/cell from a
// sample of strongly-indexed frames and apply it to experiment_, so everything below (azimuthal
// mapping, the main index+integrate pass, scaling) uses the refined geometry. Rotation is skipped
// inside (it has its own two-pass).
if (full && config_.refine_geometry.has_value())
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
// After any geometry refinement, so the shadow is found about the beam centre actually used,
// and before the azimuthal mapping and the output mask below, which both read pixel_mask_.
if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center)
PreScan(start_image, images_to_process,
config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES), observer);
AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);
JFJochReceiverPlots plots;
plots.Setup(experiment_, mapping);
// Output file (NXmxIntegrated master that links back to the original images).
StartMessage start_message;
experiment_.FillMessage(start_message);
start_message.arm_date = dataset->arm_date;
start_message.az_int_bin_to_q = mapping.GetBinToQ();
start_message.az_int_bin_to_two_theta = mapping.GetBinToTwoTheta();
start_message.az_int_q_bin_count = mapping.GetQBinCount();
start_message.az_int_phi_bin_count = mapping.GetAzimuthalBinCount();
if (mapping.GetAzimuthalBinCount() > 1)
start_message.az_int_bin_to_phi = mapping.GetBinToPhi();
start_message.pixel_mask["default"] = pixel_mask_.GetMask(experiment_);
if (full) {
start_message.rois = experiment_.ROI().ExportMetadata();
if (!experiment_.ROI().empty())
start_message.roi_map = experiment_.ExportROIMap();
start_message.max_spot_count = experiment_.GetMaxSpotCount();
}
start_message.master_suffix = "process";
start_message.file_format = FileWriterFormat::NXmxIntegrated;
start_message.write_master_file = true;
start_message.write_images = false;
// With the stride: image i of the process file is source image start_image + i * stride, so the
// linked pictures line up with the per-image analysis written beside them.
start_message.hdf5_source_data =
reader_.GetHDF5DataSource(start_image, images_to_process, config_.stride);
std::unique_ptr<FileWriter> writer;
std::unique_ptr<ProcessFileWriter> writer_queue;
if (write_files && config_.write_process_h5) {
writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
/*trusted_path=*/true);
// Deep enough that a worker never waits for the writer in the normal case, shallow enough that
// the backlog is bounded by the worker count rather than by the length of the run.
writer_queue = std::make_unique<ProcessFileWriter>(*writer, 4 * std::max(1, config_.nthreads));
}
const char *mode_name = full ? "full analysis"
: calibration ? "powder calibration" : "azimuthal integration";
logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]",
images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name);
MarkPhase(observer, full ? "Setting up" : calibration ? "Powder calibration" : "Azimuthal integration");
// Full-analysis shared engines.
std::unique_ptr<IndexerThreadPool> indexer_pool;
std::unique_ptr<IndexAndRefine> indexer;
if (per_image_analysis) {
auto indexing_settings = experiment_.GetIndexingSettings();
// Two-pass rotation fires the FFT indexer only twice - the two first-pass schemes, run in
// parallel (the main per-image loop then uses the forced lattice, no FFT). Size the pool to
// those two so a rotation run doesn't pay to initialise (cuFFT plan + GPU alloc) indexers it
// never uses. Stills fire the FFT once per image across all workers, so keep the default there.
if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) {
indexing_settings.IndexingThreads(2);
// Two indexer threads leave the rest of the machine idle for the whole first pass, which is
// a third of the run. Split each scheme's candidate-cell refinement over -N/2 threads so the
// two together use -N. The split is by candidate cell, so the result does not depend on it.
// Clamped to what IndexingSettings accepts - a big host would otherwise ask for more than
// the setting's limit and throw before a single image is read.
indexing_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
}
indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
indexer = std::make_unique<IndexAndRefine>(experiment_, indexer_pool.get());
if (!config_.reference_data.empty())
indexer->ReferenceIntensities(config_.reference_data);
// Second pass of the rotation two-pass: predict with the smoothed mosaicity fitted in the pre-pass.
if (!geometry_prepass && !prepass_mosaicity_.empty())
indexer->SetPredictionMosaicityOverride(prepass_mosaicity_);
}
const auto start_time = std::chrono::steady_clock::now();
// First pass of two-pass rotation indexing (full analysis only).
if (full && force_rotation_result_.has_value()) {
// Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation /
// search metadata / axis), not just its lattice - a lattice-only force loses that metadata and
// integrates markedly worse.
indexer->ForceRotationIndexerResult(*force_rotation_result_);
logger.Info("Rotation indexer result forced from pass 1 (lattice + geometry) - skipping first pass");
} else if (full && config_.forced_rotation_lattice.has_value()) {
indexer->ForceRotationIndexerLattice(*config_.forced_rotation_lattice);
logger.Info("Rotation indexer lattice forced externally - skipping first pass");
} else if (full && config_.rotation_indexing && config_.two_pass_rotation) {
MarkPhase(observer, "Rotation indexing (first pass)");
// 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 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,
int ordinal) {
const int image_idx = start_image + ordinal * config_.stride;
std::vector<SpotToSave> spots;
try {
if (auto img = reader_.GetRawImage(image_idx)) {
DataMessage m{};
m.number = ordinal;
m.original_number = image_idx;
auto first_pass = config_.spot_finding;
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) {
logger.Warning("First-pass spot read failed for image {}: {}", image_idx, e.what());
}
return spots;
};
// 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 {
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)
set_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) {
std::vector<int> wanted;
for (const int ordinal : ordinals)
if (!spot_cache.contains(ordinal))
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_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
e.profile = std::make_unique<AzimuthalIntegrationProfile>(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, wanted[i]);
}));
for (auto &f : futures)
f.get();
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();
};
// 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(
images_to_process, std::min(images_to_process, 60));
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()),
ThreadsForWork(validation.size(), config_.nthreads, 1), [&](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();
};
// 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;
auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector<int> &ordinals) {
auto ri = std::make_unique<RotationIndexer>(experiment_, pool);
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 whole rotation - full 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(images_to_process);
std::iota(all_consecutive.begin(), all_consecutive.end(), 0);
const std::vector<int> spread_ordinals = select_equally_spaced_image_ordinals(
images_to_process, config_.rotation_indexing_image_count);
const std::vector<std::pair<std::string, std::vector<int>>> schemes = {
{"spread", spread_ordinals},
{"wedge", all_consecutive},
};
struct FirstPass { std::optional<RotationIndexerResult> result; int score = -1; std::string name; double vol = 0.0; };
// Feed both schemes (single-threaded spot cache), run their RunIndexing() concurrently, then keep
// the lattice that indexes the most validation frames on the real per-image path.
auto pick_best = [&](IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass {
std::vector<std::unique_ptr<RotationIndexer>> ris;
ris.reserve(schemes.size());
for (const auto &[name, ordinals] : schemes) {
if (cancelled_)
break;
ris.push_back(feed_scheme(pool, ordinals));
}
std::vector<std::future<void>> index_futures;
index_futures.reserve(ris.size());
for (auto &ri : ris) {
RotationIndexer *rp = ri.get();
index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
}
for (auto &f : index_futures)
f.get();
FirstPass bp;
for (size_t i = 0; i < ris.size(); i++) {
if (cancelled_)
break;
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;
}
}
// 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, the larger cell is that spurious supercell and the
// smaller is the true reduced cell: take it, regardless of scheme order. 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;
if (bp.result.has_value() && !clearly_more && 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);
const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15;
integer_subcell = tied && integer_multiple && vol < bp.vol;
}
if (!bp.result.has_value() || clearly_more || integer_subcell) {
if (integer_subcell)
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);
bp.score = score;
bp.vol = vol;
bp.name = name;
bp.result = std::move(result);
}
}
return bp;
};
// 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, images_to_process); i++)
upfront.push_back(i);
std::sort(upfront.begin(), upfront.end());
upfront.erase(std::unique(upfront.begin(), upfront.end()), upfront.end());
prefetch_spots(upfront);
}
FirstPass best = pick_best(*indexer_pool, *indexer);
// 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 COARSE-resolution pass - only low-order reflections,
// where the fine axis stays robust - then RE-INDEX at full resolution constrained by that cell 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.
auto max_axis = [](const RotationIndexerResult &r) {
const auto uc = r.lattice.GetUnitCell();
return std::max({uc.a, uc.b, uc.c});
};
if (!cancelled_ && best.result.has_value()
&& best.score < 0.5 * static_cast<double>(validation.size())) {
auto coarse_settings = experiment_.GetIndexingSettings();
coarse_settings.FFT_HighResolution_A(3.5f); // low-order reflections only -> robust long axis
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(*best.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(*best.result));
experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell());
const FirstPass 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 > best.score)
best = constrained;
}
}
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. 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())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
std::string("Two-pass rotation indexing found no lattice.") + next_steps);
// 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.
if (best.score < static_cast<int>(validation.size()) / 6) {
// 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 JFJochException(JFJochExceptionCategory::InputParameterInvalid,
fmt::format("Two-pass rotation indexing found only a lattice that "
"indexes {}/{} validation frames - too few for it to be "
"this crystal's lattice. It was {}-centred {}, "
"{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}",
best.score, static_cast<int>(validation.size()),
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);
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)",
best.name, best.score, static_cast<int>(validation.size()));
// 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, and a centring the group
// carried over from pass 1 cannot describe is the same disagreement seen from the symmetry side.
// Both were already rejected - the first in Run(), the second before the merge - but only after
// the whole pass had run, so each 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. A triclinic
// de-novo cell is left alone, since that is the demotion the merge reindexes into the group's
// conventional setting.
if (!geometry_prepass && prepass_result_.has_value()) {
const auto *reuse_sg = prepass_merge_sg_.has_value()
? gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_)) : nullptr;
// 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;
const bool wrong_centering = reuse_sg != nullptr
&& best.result->search_result.system != gemmi::CrystalSystem::Triclinic
&& best.result->search_result.centering != reuse_sg->centring_type();
if (supercell || wrong_centering) {
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);
indexer->ForceRotationIndexerResult(*prepass_result_);
}
}
}
}
// Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads
// are serialized by the global hdf5_mutex; the analysis runs in parallel.
std::atomic<int> next_ordinal = 0;
std::atomic<int> finished_count = 0;
// Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than
// shifting everything after it up. The extent of the written datasets is therefore the highest
// ordinal reached, not the number that succeeded - counting successes makes the file one image
// short for every failure, and the image it loses is the LAST one, not the one that failed.
std::atomic<int> max_written_ordinal = -1;
auto note_written = [&max_written_ordinal](int ordinal) {
int prev = max_written_ordinal.load(std::memory_order_relaxed);
while (prev < ordinal
&& !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) {
}
};
std::atomic<uint64_t> total_uncompressed_bytes = 0;
// Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
// the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
// a fixed budget instead of everything it found - a powder ring is over-determined either way.
constexpr size_t CALIBRATION_MAX_SPOTS = 20000;
const size_t calibration_spots_per_image =
std::max<size_t>(1, CALIBRATION_MAX_SPOTS / static_cast<size_t>(images_to_process));
std::mutex calibration_spot_m;
std::vector<SpotToSave> calibration_spot_list;
auto azint_worker = [&]() {
std::vector<uint8_t> decompression_buffer;
ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_);
ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum());
AzIntEngineCPU azint(mapping);
AzimuthalIntegrationProfile profile(mapping);
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= end_image) break;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
total_uncompressed_bytes += msg.image.GetUncompressedSize();
const auto t0 = std::chrono::steady_clock::now();
try {
const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer);
preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode());
azint.Run(buffer, profile);
} catch (const std::exception &e) {
logger.Error("Error integrating image {}: {}", image_idx, e.what());
continue;
}
msg.azint_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.processing_time_s = msg.azint_time_s;
msg.az_int_profile = profile.GetResult();
msg.az_int_profile_count = profile.GetPixelCount();
msg.az_int_profile_std = profile.GetStd();
msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings());
msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
config_.spot_finding.ice_ring_width_Q_recipA);
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
};
auto full_worker = [&]() {
pin_gpu(); // round-robin per worker thread; must precede engine construction
MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
AzimuthalIntegrationProfile profile(mapping);
while (!cancelled_) {
const int ordinal = next_ordinal.fetch_add(1);
const int image_idx = start_image + ordinal * config_.stride;
if (image_idx >= end_image) break;
std::shared_ptr<JFJochReaderRawImage> img;
try {
img = reader_.GetRawImage(image_idx);
} catch (const std::exception &e) {
logger.Error("Failed to load image {}: {}", image_idx, e.what());
continue;
}
if (!img) continue;
DataMessage msg{};
msg.image = img->image;
msg.number = ordinal;
msg.original_number = image_idx;
if (dataset->efficiency.size() > image_idx)
msg.image_collection_efficiency = dataset->efficiency[image_idx];
total_uncompressed_bytes += msg.image.GetUncompressedSize();
const auto t0 = std::chrono::steady_clock::now();
try {
analysis.Analyze(msg, profile, config_.spot_finding);
} catch (const std::exception &e) {
logger.Error("Error analyzing image {}: {}", image_idx, e.what());
continue;
}
msg.processing_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
msg.run_number = experiment_.GetRunNumber();
msg.run_name = experiment_.GetRunName();
if (calibration_spots) {
std::lock_guard lock(calibration_spot_m);
for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i)
calibration_spot_list.push_back(msg.spots[i]);
}
plots.Add(msg, profile);
if (writer_queue) writer_queue->Post(msg, img);
note_written(ordinal);
if (observer) observer->OnImageProcessed(msg);
const int done = finished_count.fetch_add(1) + 1;
if (observer) observer->OnProgress(done, images_to_process);
}
};
MarkPhase(observer, "Processing images");
std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
: std::function<void()>(azint_worker);
// Workers for this loop only. It is GPU-bound - every worker builds its own analysis engine, tens
// of megabytes of device and pinned memory whose allocation implicitly synchronises the device,
// and they all queue on the same few cards - so past a handful per GPU another worker adds setup
// and contention and no throughput. Measured on a four-GPU node: the loop is flat from 8 workers
// up and slower at 48 than at 12. Every other phase still gets the full thread count, because
// each one spawns its own workers.
const int loop_workers = get_gpu_count() > 0
? std::min(config_.nthreads, std::max(8, 4 * get_gpu_count()))
: config_.nthreads;
// The per-stage means below are per WORKER, so anything reporting them has to divide by this and
// not by the thread count - on a machine with GPUs they are not the same number.
result.loop_workers = loop_workers;
std::vector<std::future<void> > futures;
futures.reserve(loop_workers);
const auto image_loop_start = std::chrono::steady_clock::now();
for (int i = 0; i < loop_workers; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &f: futures)
f.get();
// Everything the workers posted has to be on disk before the end message is written and the file
// finalized, and a write that failed - out of space, above all - has to surface here rather than
// leave a truncated file behind.
if (writer_queue)
writer_queue->Finish();
// Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it,
// so without this there is nothing to compare them against and no way to see how much of a run is
// spent outside it - on the first-pass indexing and on scaling/merging.
result.image_loop_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - image_loop_start).count();
result.cancelled = cancelled_;
result.images_processed = finished_count.load();
// Every image failing is a total failure, not a run that produced nothing: it used to be reported
// only as per-image log lines while the process still exited 0 with no output file.
if (!cancelled_ && images_to_process > 0 && result.images_processed == 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"No image could be analyzed - see the per-image errors above");
result.mean_processing_time = plots.GetMeanProcessingTime();
result.indexing_rate = plots.GetIndexingRate();
// Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in
// plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and
// writes the PONI file; nothing here is stored in the _process.h5.
if (calibration && !cancelled_) {
MarkPhase(observer, "Powder calibration");
if (calibration_spots) {
logger.Info("Powder calibration from {} pooled spots", calibration_spot_list.size());
result.calibration = CalibrateFromSpots(calibration_spot_list,
experiment_.GetDiffractionGeometry(),
config_.calibrant_ring_q);
} else {
const auto *summed = plots.GetAzIntProfileObject();
if (!summed)
throw JFJochException(JFJochExceptionCategory::CalibrationError,
"No azimuthal profile was accumulated");
result.calibration = CalibrateFromProfile(summed->GetResult(), mapping,
experiment_.GetDiffractionGeometry(),
config_.calibrant_ring_q);
}
}
// End message (also written to the file).
EndMessage end_msg;
// The EXTENT of the per-image datasets (see note_written), not the success count: those two
// differ by exactly the number of frames that failed, and the writer sizes its virtual datasets
// from this.
end_msg.max_image_number = static_cast<uint64_t>(max_written_ordinal.load() + 1);
end_msg.images_collected_count = result.images_processed;
end_msg.images_sent_to_write_count = result.images_processed;
end_msg.end_date = time_UTC(std::chrono::system_clock::now());
end_msg.run_number = experiment_.GetRunNumber();
end_msg.run_name = experiment_.GetRunName();
end_msg.bkg_estimate = plots.GetBkgEstimate();
end_msg.ice_ring_score = plots.GetIceRingScoreArray();
end_msg.ice_ring_score_mean = plots.GetIceRingScore();
end_msg.az_int_result["dataset"] = plots.GetAzIntProfile();
end_msg.indexing_rate = result.indexing_rate;
if (full && !cancelled_) {
if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) {
end_msg.rotation_lattice = rot->lattice;
// Write the refined geometry (not the nominal StartMessage values) into the _process.h5.
end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl();
end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl();
end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad();
end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad();
end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad();
if (rot->axis)
end_msg.refined_rotation_axis = rot->axis->GetAxis();
end_msg.rotation_lattice_type = LatticeMessage{
.centering = rot->search_result.centering,
.niggli_class = rot->search_result.niggli_class,
.crystal_system = rot->search_result.system
};
result.rotation_lattice_found = true;
}
result.consensus_cell = indexer->GetConsensusUnitCell();
end_msg.unit_cell = result.consensus_cell;
}
// Scaling and merging (full analysis only).
if (full && !cancelled_ && result.indexing_rate.has_value() && result.indexing_rate > 0
&& (config_.run_scaling || !config_.reference_data.empty())) {
// Scaling/merging is a long post-pass; report each sub-step as a phase so the GUI progress
// bar reflects what is happening instead of freezing on one label.
auto phase = [&](const std::string &p) { MarkPhase(observer, p); };
phase("Scaling and merging");
// Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder
// ring. Their integrated intensity is contaminated by the strong, variable ice background, so
// they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the
// fulls refit) skips them, while the combine, the merge and the statistics keep them - dropping
// them outright guts low/mid-resolution completeness on crystals that merge them fine.
// The eleven hexagonal-ice bands cover 16-26% of the unique reflections at typical resolutions
// whether or not the crystal has any ice, so handling ice unconditionally taxes clean data for
// nothing. Gate the flagging on the run's own measured ice strength (1 = no ice; the strongest
// ring over the smooth radial background, measured on the peak-excluded profile so Bragg peaks
// cannot pose as ice). The score is logged either way, so a run that gates off still reports why.
// Ice is looked for on TWO channels, because it arrives in two forms. Fine-grained ice is a
// smooth powder ring and shows up in the radial profile (ice_ring_score). Ice in large
// crystallites arrives as discrete reflections instead, leaves the profile flat, and is only
// visible as a pile-up of found SPOTS on the ring positions - measured against the same q width
// of ice-free flanks beside them, so the fall-off of spot density with resolution cancels. On
// the rotation battery the two channels barely overlap: the crystals with smooth ice read
// 2.1-2.4 / ~1.0, the ones whose spots are mostly ice read ~1.1 / 3.8-14.6, and a clean crystal
// reads 1.04 on both. Either channel alone is blind to half the cases.
const auto ice_score = plots.GetIceRingScore();
const auto ice_spot_ratio = plots.GetIceRingSpotRatio();
const float ice_min_score = experiment_.GetScalingSettings().GetIceMinScore();
const float ice_min_spot_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio();
const bool ice_present = (!ice_score.has_value() || *ice_score >= ice_min_score)
|| (ice_min_spot_ratio > 0.0f && ice_spot_ratio.has_value()
&& *ice_spot_ratio >= ice_min_spot_ratio);
if (experiment_.IsDetectIceRings() && !ice_present) {
logger.Info("Ice-ring handling: measured ice score {:.2f} < {:.2f} and spot ratio {:.2f} < "
"{:.2f}, no ice detected - ice-ring handling skipped entirely",
*ice_score, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio);
} else if (experiment_.IsDetectIceRings()) {
const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA;
size_t total = 0, flagged = 0;
for (auto &outcome : indexer->GetIntegrationOutcome()) {
for (auto &r : outcome.reflections) {
++total;
r.on_ice_ring = IsOnIceRing(r.d, ice_width);
if (r.on_ice_ring)
++flagged;
}
}
logger.Info("Ice-ring handling: ice score {:.2f} (gate {:.2f}), spot ratio {:.2f} (gate "
"{:.2f}); flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
"excluded from scaling, kept for merging",
ice_score.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN),
ice_min_spot_ratio, flagged, total, ice_width);
}
// Scale the images and merge. Factored so it can run twice: first in P1 to give the
// space-group search a merged dataset, then again in the determined space group so the
// scaling sees symmetry equivalents and the final statistics are in the right symmetry.
// (With a user-fixed space group it simply runs once, already in that symmetry.)
struct ScaleMergeResult {
std::vector<MergedReflection> merged;
MergeStatistics statistics;
};
// The reference path computes each image's G once (per-image scaling against the
// reference); the scaling loop below is skipped, so G is stable across the two passes.
// Smoothing it more than once would compound the correction, so do it only on the first
// pass. The no-reference path recomputes G from scratch each pass and re-smooths correctly.
// Rotation self-scaling + 3D combine + merge is done by the dedicated RotationScaleMerge (a single
// allocate-once engine that recomputes partiality from the fitted mosaicity, combines the per-frame
// partials, scales the fulls and merges; the whole hot path runs on the GPU when one is present).
// Ingested once here (after ice flagging) and reused across both space-group passes; a fixed
// (forced) mosaicity is handled by the recompute. The decay + absorption correction surfaces run as
// post-scale-fulls stages (on by default, see ScalingSettings::CorrectionSurfaces); external-
// reference scaling, the stills -B (per-image B-factor) and wedge refinement are unsupported here.
// Two-pass second pass: if the de-novo indexer DEMOTED to a triclinic (primitive) cell while the reused
// merge space group is higher-symmetry (a huge oblique cell whose constrained centred refine was
// ill-posed, so the pseudo-symmetry guard kept the primitive), the reflections are in the primitive
// frame but the merge expects the group's CONVENTIONAL setting - reindex them into it (the same
// primitive->conventional change of basis the space-group search applies to a demoted result), else the
// merge folds the wrong equivalents. No-op when the indexer already produced the group's setting.
if (prepass_merge_sg_ && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()
&& end_msg.rotation_lattice_type->crystal_system == gemmi::CrystalSystem::Triclinic) {
const auto *msg_sg = gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_));
if (msg_sg && msg_sg->crystal_system() != gemmi::CrystalSystem::Triclinic) {
const auto cand = LatticeSearch(*end_msg.rotation_lattice);
// Exact integer reindex prim->conv: P[i][j] = conv_real[i] . rot_reciprocal[j] (as the
// centred-lattice test), robust to the refined cell not being exactly LatticeSearch's Niggli cell.
const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
end_msg.rotation_lattice->Cstar()};
gemmi::Mat33 reindex;
double reindex_res = 0.0;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
const double v = cv[i] * rs[j];
reindex.a[i][j] = std::round(v);
reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
}
// Only when LatticeSearch recovers the reused group's own metric (system + centring): otherwise
// the reused group and the indexed metric disagree and reindexing would be wrong - leave it.
if (cand.system == msg_sg->crystal_system() && cand.centering == msg_sg->centring_type()
&& reindex_res < 0.1) {
for (auto &io : indexer->GetIntegrationOutcome()) {
io.latt = io.latt.Multiply(reindex);
for (auto &r : io.reflections) {
const double h = r.h, k = r.k, l = r.l;
r.h = static_cast<int32_t>(std::lround(reindex.a[0][0]*h + reindex.a[0][1]*k + reindex.a[0][2]*l));
r.k = static_cast<int32_t>(std::lround(reindex.a[1][0]*h + reindex.a[1][1]*k + reindex.a[1][2]*l));
r.l = static_cast<int32_t>(std::lround(reindex.a[2][0]*h + reindex.a[2][1]*k + reindex.a[2][2]*l));
}
}
result.consensus_cell = cand.conventional.GetUnitCell();
end_msg.unit_cell = result.consensus_cell;
end_msg.rotation_lattice = end_msg.rotation_lattice->Multiply(reindex);
end_msg.rotation_lattice_type = LatticeMessage{ .centering = cand.centering,
.niggli_class = end_msg.rotation_lattice_type->niggli_class, .crystal_system = cand.system };
logger.Info("Two-pass: reindexed the de-novo primitive cell into the space-group {} "
"conventional setting for the merge", static_cast<int>(*prepass_merge_sg_));
}
}
}
// Two-pass second pass only: the indexer above ran DE NOVO (Run() cleared the space group so its
// pseudo-symmetry safeguards find the true cell); now reinstate pass-1's space group for the merge,
// exactly as a user -S would, so the symmetry is fixed without a re-search. No-op on pass 1 / single pass.
//
// Only when the lattice pass 2 actually found HAS that group's centring. Re-indexing de novo can
// land in a different setting from pass 1 - most often the primitive sub-cell of a centred
// lattice, which is what the reindex above exists to undo, and which it declines to do when the
// metric does not match. Stamping a centred group onto a primitive cell is not a small error:
// the centring absence rule then removes half the reflections that genuinely exist, so the merge
// reports more unique reflections than its own cell can hold (measured: a C2 group on a cell of
// exactly half the C-centred volume, reported as 117% complete with CC1/2 0.62 against the first
// pass's 92.6% and 0.98). When they disagree, let pass 2 search for itself rather than trust a
// group that belongs to a different lattice.
if (prepass_merge_sg_ && end_msg.rotation_lattice_type.has_value()) {
const auto *reuse_sg = gemmi::find_spacegroup_by_number(static_cast<int>(*prepass_merge_sg_));
if (reuse_sg && end_msg.rotation_lattice_type->centering != reuse_sg->centring_type()) {
logger.Warning("Two-pass: the second pass indexed a {}-centred lattice, but the first pass "
"determined {}, whose lattice is {}-centred - the group cannot describe this "
"lattice, so this pass will not be adopted",
end_msg.rotation_lattice_type->centering, reuse_sg->xhm(),
reuse_sg->centring_type());
result.lattice_conflicts_with_prepass_sg = true;
}
}
if (prepass_merge_sg_)
experiment_.SpaceGroupNumber(*prepass_merge_sg_);
const auto &rot_ss = experiment_.GetScalingSettings();
const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge
std::optional<RotationScaleMerge> rsm;
if (is_rotation) {
if (experiment_.GetRefineRotationWedgeInScaling()
|| rot_ss.GetRotationWedgeForScaling().has_value())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Rotation scaling/merging (RotationScaleMerge) does not support "
"wedge refinement");
// A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and
// resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays
// self-consistent.
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
}
auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult {
if (rsm) {
phase("Scale/combine/merge (" + label + ")");
auto r = rsm->Run(for_search);
result.error_model_isa = r.isa;
result.error_model_isa_asymptotic = r.isa_asymptotic;
result.error_model_a = r.error_model_a;
result.error_model_b = r.error_model_b;
return ScaleMergeResult{std::move(r.merged), std::move(r.statistics)};
}
// Stills (rotation goes through RotationScaleMerge above): self-scale each image against the
// running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even
// with an external reference: the reference is used only to break the per-image indexing
// ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor (scaling each
// image against a foreign dataset injects cross-dataset systematics and is a worse reference
// than the data's own merge). One pass: the per-image scale G is the exact one-pass solution,
// so iterating would only rebuild the reference from the freshly-scaled (noisy) data and
// degrade weak stills (measured CC1/2 collapse at the default 3 iters).
{
phase("Scaling images (" + label + ")");
auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome());
indexer->ScaleAllImages(merge_result);
}
// Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
// orientation tilt against the running merge and recompute each reflection's partiality + scale
// correction (no re-integration). Never for the P1 search pass.
if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) {
phase("Partiality post-refine (" + label + ")");
StillsPartialityRefine refiner(experiment_);
const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads);
logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
// The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally
// when the image has no prior CC - which is exactly the state the guard above leaves. So
// the guard was protecting nothing: measured on a jet dataset, of 367 images it left
// unscaled only 8 were still unscaled in the output, and 53 reached the merge at or below
// a fiftieth of the run median, the worst at a 4525th. Whatever produced the final scale,
// it is the final scale that has to be credible.
ScaleOnTheFly::RejectCollapsedScales(indexer->GetIntegrationOutcome());
}
const std::vector<IntegrationOutcome> &merge_input = indexer->GetIntegrationOutcome();
phase("Merging");
MergeOnTheFly merge_engine(experiment_);
// For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged
// intensities and the error model, so the space-group search sees clean data; the final
// in-symmetry merge keeps them so completeness is not lost.
merge_engine.ExcludeIceRings(for_search);
if (result.consensus_cell.has_value())
merge_engine.ReferenceCell(*result.consensus_cell);
// Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells)
// most indexed crystals do not correlate with the true structure; --min-image-cc drops the
// crystals whose per-image CC to the reference is below the limit, so the merge keeps only
// the coherent (real) ones. Set before the error model so the model, the merge and the
// reported statistics are all fitted over the same set of images. The P1 search pass keeps
// every image: its job is to find the symmetry, not to produce final intensities.
merge_engine.FilterByImageCC(!for_search
&& experiment_.GetScalingSettings().GetMinCCForImage() > 0.0);
merge_engine.RefineErrorModel(merge_input);
if (merge_engine.ErrorModelActive())
{
// XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is
// already the whole-range 1/b, there being no asymptote tier on this path.
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}",
em.a, em.b, em.isa, merge_engine.ErrorModelChi2());
}
for (size_t i = 0; i < merge_input.size(); ++i)
merge_engine.AddImage(merge_input[i], static_cast<int64_t>(i));
ScaleMergeResult out;
out.merged = merge_engine.ExportReflections();
// Automatic high-resolution cutoff (post-merge): trim the written reflections + reported
// shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full
// range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A
// manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the
// space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.)
const auto &cut_ss = experiment_.GetScalingSettings();
const std::optional<double> effective_d_min = ApplyResolutionCutoff(
out.merged, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(),
cut_ss.GetResolutionCCTarget(), for_search, logger);
phase("Computing statistics");
out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min);
{ // stills: no asymptote tier, and its ISa already IS the whole-range 1/b
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
result.error_model_isa = em.isa;
result.error_model_a = em.a;
result.error_model_b = em.b;
}
logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label);
return out;
};
// First pass: P1 when searching, or directly the user-fixed space group.
const auto initial_sg = experiment_.GetGemmiSpaceGroup();
auto sm = scale_and_merge(initial_sg ? initial_sg->short_name() : "P1", !initial_sg.has_value());
// Rotation two-pass geometry pre-pass: this first scale/merge has now fitted a frame-order-smoothed
// mosaicity and written it back onto the per-frame outcomes. Capture it for the second pass's Bragg
// prediction (so prediction uses the smoothed value, not a fresh per-image estimate) and post-refine +
// apply the detector geometry from this integration. The pre-pass then CONTINUES through the space-
// group search (below) so the second pass can reuse the same space group and lattice - re-searching at
// the slightly changed geometry could flip a borderline determination and mix the two passes.
if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) {
const auto &outcomes = indexer->GetIntegrationOutcome();
prepass_mosaicity_.assign(outcomes.size(), NAN);
for (size_t o = 0; o < outcomes.size(); ++o)
if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg;
if (result.consensus_cell.has_value()) {
if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) {
prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback
PostRefineSettings prs;
prs.crystal_system = rot->search_result.system;
prs.num_threads = static_cast<int>(config_.nthreads);
prs.refine_geometry = true;
const auto pr = PostRefineRotationGeometry(outcomes, *rot->axis,
experiment_.GetDiffractionGeometry(), rot->lattice, prs, logger);
result.post_refine = pr;
if (pr.ok && pr.detector_refined) {
// Store, don't apply: this pass continues to write its own (header-geometry) output;
// Run() applies the refined geometry before the second pass.
prepass_detector_geometry_ = std::array<float, 3>{
static_cast<float>(pr.beam_x_after_px),
static_cast<float>(pr.beam_y_after_px),
static_cast<float>(pr.distance_after_mm)};
logger.Info("Two-pass: post-refined detector geometry (distance {:.3f} mm, "
"beam {:.2f},{:.2f} px) will drive the second integration pass",
pr.distance_after_mm, pr.beam_x_after_px, pr.beam_y_after_px);
} else {
logger.Info("Two-pass: geometry post-refine committed no detector change "
"- the second pass reproduces the first");
}
// DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the
// scale is its own cross-validated fit and the crystals that have a stage fault are
// exactly the ones whose cell and detector steps do NOT pass, because the angle error
// is what their residual is made of. A calibration fault is rare (36 of 37 rotation
// datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset
// would silently damage it, so the asymmetry is deliberate: committed only when the
// fit is both cross-validated and outside the tolerance. A manual --rotation-scale is
// already on the goniometer and is left alone.
if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value())
prepass_rotation_scale_ = static_cast<float>(pr.rotation_scale);
}
}
}
// For a de-novo search, gate the P1 scaling+merge (the pass that feeds the space-group search) at
// <I/sigma> >= 1 over thin resolution shells: noise-dominated high-res shells otherwise mislead
// the determination, where a single artifact (a zinger the box sum integrates whole, or a
// weak-reflection profile-fit runaway) has an astronomical resolution-normalised E and wrecks the
// intensity-correlation / second-moment statistics. The final in-symmetry merge keeps the full
// range (only this P1 search pass is cut); a manual --scaling-high-resolution, if coarser, wins.
double d_min_search = 0.0;
if (!initial_sg.has_value()) {
std::vector<std::pair<float, float>> rs; // (d, I/sigma) over the full P1 merge
rs.reserve(sm.merged.size());
for (const auto &m : sm.merged)
if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0f &&
std::isfinite(m.d) && m.d > 0.0f)
rs.emplace_back(m.d, m.I / m.sigma);
if (rs.size() >= 400) {
std::sort(rs.begin(), rs.end(),
[](const auto &a, const auto &b) { return a.first > b.first; }); // low -> high res
const int bins = std::clamp(static_cast<int>(rs.size() / 100), 1, 40);
const size_t per = (rs.size() + bins - 1) / static_cast<size_t>(bins);
for (size_t b = 0; b * per < rs.size(); ++b) {
const size_t lo = b * per, hi = std::min(rs.size(), lo + per);
double sum = 0.0;
for (size_t j = lo; j < hi; ++j) sum += rs[j].second;
if (sum / static_cast<double>(hi - lo) < 1.0) {
// Cut the noise-dominated high-res shells. When even the lowest-res shell fails,
// keep that shell alone rather than abandoning the cut. The bound is absolute
// while the merged I/sigma it tests is capped by the merge's own ISa, so a merge
// whose ISa has fallen below 1 cannot satisfy it in ANY shell - and that is
// exactly the merge the search must not be handed whole, since the noise which
// drove ISa down is what the cut exists to remove. Abandoning it feeds back:
// more noise into the search, a weaker operator statistic, and a lost symmetry.
if (b > 0)
d_min_search = rs[lo].first;
else if (per >= 400)
d_min_search = rs[per - 1].first;
break;
}
}
}
const double d_manual = experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0);
if (rsm && d_min_search > d_manual) {
logger.Info("Space-group search: re-scaling the P1 pass at {:.2f} A (<I/sigma> >= 1)", d_min_search);
rsm->SetDMinLimit(d_min_search);
sm = scale_and_merge("P1, space-group search", true);
rsm->SetDMinLimit(d_manual > 0.0 ? std::optional<double>(d_manual) : std::nullopt);
}
}
std::ostringstream stats_text;
// True when the point group below was chosen BY THE SEARCH (rather than given by the user) and is
// above triclinic - i.e. a promotion was made, which is what makes a later "the Laue class is
// holohedral so there is no twin law" conclusion circular.
// Pass 2 reuses pass-1's group instead of searching, so the promotion it needs to report was
// made there; a first / single pass has no carried-over group and decides this below.
bool promoted_point_group = prepass_merge_sg_.has_value() && prepass_promoted_point_group_;
if (!experiment_.GetGemmiSpaceGroup().has_value()) {
SearchSpaceGroupOptions sg_opts;
sg_opts.nthreads = static_cast<size_t>(std::max(1, config_.nthreads));
sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel();
sg_opts.d_min_limit_A = std::max<double>(
d_min_search, experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0));
// Constrain the search to subgroups of the lattice (metric) symmetry found by rotation
// indexing. Centering is not constrained here - it is determined from the absences.
if (end_msg.rotation_lattice_type.has_value())
sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system;
auto sg_search = SearchSpaceGroup(sm.merged, sg_opts);
// Second opinion from a merge that keeps only well-measured observations (see
// RotationScaleMerge::search_min_zeta). Where the two disagree the all-observation merge
// decides. The filter earns its place because a reflection that crosses the Ewald sphere
// near-tangentially is measured worst, and a handful of bad measurements can make a real
// operator look like a twin law; but it discards 40-80% of the observations, which starves
// the operator correlations. Both failure directions are real - the filtered arm has been
// seen to lose an operator the full merge confirms, and, at geometries away from the ones
// this comparison was first checked on, to confirm one the full merge refuses. Only one of
// the two has every observation behind it, so the filtered merge reports and the
// all-observation merge decides. Both searches keep their vetoes, so a merohedral twin
// stays refused in each (verified on a real twin at the three frame ranges where it
// over-promotes: refused in both arms every time).
if (rsm && rsm->GetSearchMinZeta() > 0.0 && !sg_search.point_group_hm.empty()) {
const double zeta = rsm->GetSearchMinZeta();
rsm->SetSearchMinZeta(0.0);
auto sm_all = scale_and_merge("P1, all observations", true);
rsm->SetSearchMinZeta(zeta);
const auto alt = SearchSpaceGroup(sm_all.merged, sg_opts);
// The order of the point group each arm confirmed. Read from the search, not from the
// space group it went on to pick: Stage B leaves best_space_group unset when no
// candidate has enough absences to be eligible, which is common on the filtered arm
// precisely because the filter removes the weak reflections the absences live in. An
// arm that confirmed 422 would then have scored 0 and lost to a P2 arm - a demotion,
// reported as "taking the higher symmetry".
const auto order_of = [](const SearchSpaceGroupResult &r) { return r.point_group_order; };
logger.Info("Space-group search: all-observation merge -> {} (order {}), "
"Lorentz-filtered -> {} (order {})",
alt.point_group_hm.empty() ? "?" : alt.point_group_hm, order_of(alt),
sg_search.point_group_hm.empty() ? "?" : sg_search.point_group_hm,
order_of(sg_search));
// The screw/centering determination in particular comes from the merge with all the
// observations, because systematic absences are decided by the WEAK reflections and the
// filter throws most of them away. Preferring the filtered arm on a tie cost four
// crystals their screw axes (P2(1) read as P2) for exactly that reason, with the point
// group and every intensity statistic identical.
// Checked while both arms are still intact - below, the all-observation arm is adopted
// by value and the two become the same result.
// Which point group, by its symmorphic representative rather than by name - gemmi calls
// both P321 and P312 "32", so two arms that disagree about WHICH 2-folds are real read
// as agreeing. 0 when no point group was chosen.
const auto pg_of = [](const SearchSpaceGroupResult &r) {
return r.point_group_representative.has_value() ? r.point_group_representative->number : 0;
};
if (order_of(alt) > 0 && order_of(alt) == order_of(sg_search)
&& pg_of(alt) != pg_of(sg_search)) {
// Same order, different symmetry: the two merges disagree about WHICH operators are
// real, and neither is higher, so there is nothing to prefer. Say so instead of
// picking silently - the two imply different molecular-replacement searches, and
// trying both is cheap next to processing the data again.
const std::string a = alt.best_space_group.has_value()
? alt.best_space_group->short_name() : alt.point_group_hm;
const std::string b = sg_search.best_space_group.has_value()
? sg_search.best_space_group->short_name() : sg_search.point_group_hm;
const std::string msg = fmt::format(
"Space group is AMBIGUOUS between {} and {} (point groups {} and {}, the same "
"order): the merge of all observations and the merge of only the well-measured "
"ones each support one of them, and neither is the higher symmetry, so the data "
"do not decide. Processing continues in {} (the all-observation choice) - try "
"BOTH in molecular replacement, or re-run with -S {} to force the other one.",
a, b, alt.point_group_hm, sg_search.point_group_hm, a, b);
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
}
if (order_of(alt) > order_of(sg_search))
logger.Info("Space-group search: all-observation merge supports {} where the "
"Lorentz-filtered one supports {} - taking the higher symmetry",
alt.point_group_hm, sg_search.point_group_hm);
else if (order_of(alt) < order_of(sg_search))
logger.Info("Space-group search: Lorentz-filtered merge supports {} where the "
"all-observation one supports {} - keeping the all-observation "
"answer, which has every observation behind it",
sg_search.point_group_hm, alt.point_group_hm);
sg_search = alt;
// The all-observation merge is the one carried forward: it is what the
// pre-promotion twinning numbers below are measured on, and the filtered merge is the
// second opinion reported above.
sm = std::move(sm_all);
}
// Twinning evidence measured BEFORE any promotion, on the P1/subgroup merge the search was
// given. Reported alongside the post-adoption analysis: once a higher Laue class has been
// adopted, its own twinning test can only say "no twin law exists within this class", which
// is circular - the promotion is what a twin would have caused. These numbers are not.
const auto pre_promotion_twinning = AnalyzeTwinning(sm.merged, nullptr);
if (pre_promotion_twinning.twinning_suspected)
logger.Warning("Twinning indicators BEFORE the space-group decision (subgroup merge): "
"<|L|> = {:.3f}, second moment = {:.3f} (untwinned 0.500 / 2.00) - if the "
"search promoted the point group, treat the promotion with suspicion",
pre_promotion_twinning.mean_abs_l, pre_promotion_twinning.second_moment);
promoted_point_group = !sg_search.point_group_hm.empty() && sg_search.point_group_hm != "1";
if (!sg_search.refused_point_group_hm.empty())
logger.Warning("Higher symmetry {} was confirmed by the operator correlations but refused: "
"{}. Processing in the lower symmetry (the recoverable direction).",
sg_search.refused_point_group_hm, sg_search.refused_reason);
// Miller-index reindex under a change of basis: (h,k,l)_conv = reindex * (h,k,l)_prim.
const auto reindex_hkl = [](auto &r, const gemmi::Mat33 &m) {
const double h = r.h, k = r.k, l = r.l;
r.h = static_cast<int32_t>(std::lround(m.a[0][0] * h + m.a[0][1] * k + m.a[0][2] * l));
r.k = static_cast<int32_t>(std::lround(m.a[1][0] * h + m.a[1][1] * k + m.a[1][2] * l));
r.l = static_cast<int32_t>(std::lround(m.a[2][0] * h + m.a[2][1] * k + m.a[2][2] * l));
};
// Intensity-based centred-lattice test (XDS CORRECT / POINTLESS style). When the indexer
// demoted to a primitive/triclinic cell but the metric is pseudo-symmetric, the point-group
// search above tested no operators (its candidate set is bounded by the indexer's lattice).
// A genuinely centred lattice and a merely pseudo-symmetric triclinic one are indistinguishable
// from spot POSITIONS - the 2-fold is an INTENSITY property - so decide it here, after P1
// integration: reindex the P1-merged data to the metric candidate's conventional setting and
// re-run the search. The I(h) vs I(Rh) correlation confirms a real centred lattice (high CC ->
// commit the centring) and rejects a pseudo-symmetric triclinic (low CC -> stays P1, no reindex).
std::optional<gemmi::Mat33> commit_reindex;
std::optional<UnitCell> commit_cell;
std::optional<CrystalLattice> commit_lattice;
std::optional<LatticeMessage> commit_lattice_type;
// Only a demoted (triclinic) indexing result reaches here with a primitive cell - that is both
// the case the point-group search could not test (a triclinic holohedry offers no operators)
// and the case where LatticeSearch's Niggli reduction is valid (it assumes a primitive input).
const bool searched_trivial = !sg_search.best_space_group || sg_search.best_space_group->number == 1;
const bool demoted_triclinic = end_msg.rotation_lattice_type.has_value()
&& end_msg.rotation_lattice_type->crystal_system == gemmi::CrystalSystem::Triclinic;
if (searched_trivial && demoted_triclinic && end_msg.rotation_lattice.has_value()) {
// The metric only generates the candidate; the intensity CC gate below makes the decision
// (a false pseudo-symmetry is rejected by the correlation).
const auto cand = LatticeSearch(*end_msg.rotation_lattice);
// Exact integer reindex from the indexed setting straight to the candidate conventional
// setting: P[i][j] = conv_real[i] . rot_reciprocal[j]. Derived from the two lattices, so it
// holds even if the refined indexed cell is not exactly LatticeSearch's Niggli cell (whose
// primitive_reduced->conventional reindex would otherwise be applied in the wrong setting).
const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
end_msg.rotation_lattice->Cstar()};
gemmi::Mat33 reindex;
double reindex_res = 0.0;
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j) {
const double v = cv[i] * rs[j];
reindex.a[i][j] = std::round(v);
reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
}
if (cand.system != gemmi::CrystalSystem::Triclinic && reindex_res < 0.1) {
std::vector<MergedReflection> merged_c = sm.merged;
for (auto &m : merged_c)
reindex_hkl(m, reindex);
SearchSpaceGroupOptions o2 = sg_opts;
o2.lattice_system = cand.system;
const auto s2 = SearchSpaceGroup(merged_c, o2);
if (s2.best_space_group.has_value() && s2.best_space_group->number > 1) {
// The centring cannot be confirmed from absences here (integrated in the primitive
// cell, so the centring-absent reflections do not exist) - it is metric-determined.
// Among the point-group-equivalent candidates the intensities cannot separate
// (e.g. P2 / P21 / C2), pick the one matching the metric candidate's centring.
gemmi::SpaceGroup chosen = *s2.best_space_group;
if (chosen.centring_type() != cand.centering)
for (const auto &alt : s2.alternatives)
if (alt.centring_type() == cand.centering) { chosen = alt; break; }
// Commit only if the adopted centring matches the metric candidate whose conventional
// cell we reindex into - otherwise the space group and the written cell disagree.
if (chosen.centring_type() == cand.centering) {
sg_search = s2;
// Swap the metric-matching candidate into the answer and leave the one it
// displaced among the alternatives, so the reported "or ..." list stays the
// set of groups the data cannot separate, each named once.
if (chosen.number != s2.best_space_group->number) {
sg_search.alternatives.push_back(*s2.best_space_group);
std::erase_if(sg_search.alternatives, [&](const gemmi::SpaceGroup &a) {
return a.number == chosen.number;
});
}
sg_search.best_space_group = chosen;
commit_reindex = reindex;
commit_cell = cand.conventional.GetUnitCell();
commit_lattice = end_msg.rotation_lattice->Multiply(reindex);
commit_lattice_type = LatticeMessage{ .centering = cand.centering,
.niggli_class = end_msg.rotation_lattice_type->niggli_class,
.crystal_system = cand.system };
}
}
}
}
// Adopt the determined space group and re-scale/merge in it, so scaling uses symmetry
// equivalents and the statistics come out in the right symmetry. P1 stands when nothing
// is determined.
if (sg_search.best_space_group.has_value()) {
const auto &sg = *sg_search.best_space_group;
logger.Info("Adopting space group {} (number {})", sg.short_name(), sg.number);
// A reindex is committed only when the intensity re-test confirmed a higher symmetry in a
// centred setting: bring the integrated reflections + cell into that setting and re-ingest
// the rotation merge, so the final scaling/merging folds the equivalents correctly.
if (commit_reindex) {
for (auto &io : indexer->GetIntegrationOutcome()) {
io.latt = io.latt.Multiply(*commit_reindex);
for (auto &r : io.reflections)
reindex_hkl(r, *commit_reindex);
}
result.consensus_cell = commit_cell;
end_msg.unit_cell = commit_cell;
// Keep the lattice metadata (hence the master-file UB matrix and indexed-lattice
// vectors) in the same conventional setting as the cell, reflections and space group.
end_msg.rotation_lattice = commit_lattice;
end_msg.rotation_lattice_type = commit_lattice_type;
if (rsm) {
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
}
}
experiment_.SpaceGroupNumber(sg.number);
end_msg.space_group_number = sg.number;
result.space_group_number = sg.number;
// The pre-pass wanted the GROUP, which is now determined. Merging in it is what the
// second pass does, at the refined geometry, and that merge is the one anybody reads.
if (!geometry_prepass) {
phase("Re-scaling in space group " + sg.short_name());
sm = scale_and_merge(sg.short_name(), false);
}
}
result.space_group_search = sg_search;
result.promoted_point_group = promoted_point_group;
} else {
// A space group was fixed by the user; surface it so the viewer/CLI can still show it.
result.space_group_number = experiment_.GetSpaceGroupNumber();
}
// Rotation two-pass geometry pre-pass: the determined space group is now fixed on experiment_, so the
// second pass REUSES it instead of re-searching at the slightly changed geometry (which could flip a
// borderline determination and mix the two passes). The lattice is deliberately NOT forced - the second
// pass re-indexes at the refined geometry so the CELL comes out self-consistent with it (forcing the
// pre-pass lattice would pin the nominal cell and undo the geometry refinement).
//
// And that is everything the pre-pass exists to measure: the post-refined detector geometry, the
// goniometer rotation scale, the frame-order-smoothed mosaicity, the space group and (for the
// supercell guard) its indexing result. Everything below produces the merged intensities and the
// report and files that go with them, and the second pass makes all of it again at the refined
// geometry - so for the pre-pass it is a merge nobody reads. Measured over the rotation test set,
// running it anyway was 10% of the battery.
if (!geometry_prepass) {
// Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and
// the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group),
// pick the reindexing that best correlates with the reference intensities and re-merge in it. The
// twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged).
// Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here.
if (rsm && !config_.reference_data.empty() && result.consensus_cell
&& experiment_.GetSpaceGroupNumber().has_value()) {
const int sg_num = static_cast<int>(*experiment_.GetSpaceGroupNumber());
const auto choice = ChooseReindex(
sm.merged, *result.consensus_cell, sg_num,
[&](const std::vector<MergedReflection> &m) {
return ReferenceIntensityCC(m, config_.reference_data, sg_num);
});
if (!choice.is_identity) {
logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference "
"(CC {:.3f}, vs {:.3f} unchanged)", choice.score, choice.identity_score);
for (auto &io : indexer->GetIntegrationOutcome())
for (auto &r : io.reflections) {
const gemmi::Op::Miller h = choice.op.apply_to_hkl({{r.h, r.k, r.l}});
r.h = h[0]; r.k = h[1]; r.l = h[2];
}
rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
static_cast<int>(config_.scaling_iter),
config_.nthreads, logger, config_.observation_dump_path);
rsm->Ingest();
phase("Re-merging in the reference indexing");
sm = scale_and_merge(experiment_.GetGemmiSpaceGroup()->short_name(), false);
}
}
const auto twin_sg_number = experiment_.GetSpaceGroupNumber();
const gemmi::SpaceGroup *twin_sg = twin_sg_number
? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr;
result.twinning = AnalyzeTwinning(sm.merged, twin_sg);
// Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the
// search itself made, so the text cannot claim "no twin law exists" on its own say-so.
result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group;
stats_text << TwinningAnalysisToText(result.twinning) << "\n";
// Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the
// crystal was mounted, which the user can change), warned about when the two nearly coincide.
if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) {
if (const auto gonio = experiment_.GetGoniometer()) {
const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice,
gonio->GetAxis());
if (closest.has_value()) {
// Practical bound, not a derived one: the blind cone's half-angle is the maximum
// Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg
// case shows the loss while a 16.2 deg one is 99.7% complete.
constexpr double WARN_DEG = 15.0;
const auto [angle, order] = *closest;
if (angle < WARN_DEG) {
const std::string msg = fmt::format(
"The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep "
"never records the reflections whose reciprocal vector lies within the Bragg angle "
"of the spindle, and symmetry normally supplies them from an equivalent elsewhere; "
"it cannot here, because that axis maps the blind region onto itself. Those "
"reflections stay missing however long the sweep runs. The loss is confined to "
"that cone rather than spread over the data, so overall completeness may still "
"look reasonable - check it near the spindle direction. A second sweep on a "
"different axis, or re-mounting, recovers them.",
order, angle);
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
} else {
stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at "
<< std::fixed << std::setprecision(1) << angle << " deg\n\n";
}
}
}
}
// Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds
// the Laue symmetry the crystal can be validly indexed in several hands related by twin-law
// operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral
// metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The
// reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform.
if (result.consensus_cell && twin_sg_number) {
const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell,
static_cast<int>(*twin_sg_number), 2.0);
if (!twin_ops.empty()) {
std::string laws;
for (const auto &op : twin_ops) {
std::string t = op.triplet();
std::replace(t.begin(), t.end(), 'x', 'h');
std::replace(t.begin(), t.end(), 'y', 'k');
std::replace(t.begin(), t.end(), 'z', 'l');
laws += (laws.empty() ? "" : " ; ") + t;
}
if (config_.reference_data.empty()) {
const std::string msg = fmt::format(
"Indexing ambiguity: this cell / space group admits alternative indexing "
"(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only "
"break this against an external reference. WITHOUT one the merge mixes the hands and "
"CC1/2 is degraded - supply a reference (-z reference.mtz or --model) to resolve it.",
laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals");
logger.Warning("{}", msg);
stats_text << " !! " << msg << "\n\n";
result.warnings.push_back(msg);
} else {
logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the "
"supplied reference.", laws);
}
}
}
// Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not
// fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log.
{
const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged);
sm.statistics.wilson_b = wilson.b;
sm.statistics.wilson_b_correlation = wilson.correlation;
if (std::isfinite(wilson.b) && wilson.b > 0.0)
logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)",
wilson.b, wilson.correlation, wilson.n_shells);
}
stats_text << sm.statistics;
result.merge_statistics_text = stats_text.str();
result.has_merge_statistics = true;
result.merge_statistics = sm.statistics;
// Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where
// the image is in no flagged range, plus the vocabulary the codes index. Only filled when the
// diagnostic ran, so absent datasets mean "not looked for" rather than "all clean".
if (sm.statistics.sweep_quality.measured) {
end_msg.sweep_quality.assign(end_msg.max_image_number, 0);
for (const auto &r : sm.statistics.sweep_quality.ranges)
for (int64_t i = std::max<int64_t>(0, r.first_image);
i <= r.last_image && i < static_cast<int64_t>(end_msg.sweep_quality.size()); ++i)
end_msg.sweep_quality[i] = static_cast<uint8_t>(r.reason) + 1;
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
end_msg.sweep_quality_reasons.emplace_back(
SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
}
{
// Stride rather than take the head: the merged list is ordered by hkl, so the first N
// reflections are one corner of reciprocal space and would not show the intensity range.
constexpr size_t MAX_POINTS = 4000;
const size_t n = sm.merged.size();
const size_t step = std::max<size_t>(1, n / MAX_POINTS);
result.merged_i_sigma.reserve(std::min(n, MAX_POINTS + 1));
for (size_t i = 0; i < n; i += step)
if (std::isfinite(sm.merged[i].I) && std::isfinite(sm.merged[i].sigma)
&& sm.merged[i].sigma > 0.0f)
result.merged_i_sigma.emplace_back(sm.merged[i].I, sm.merged[i].sigma);
}
result.has_reference = !config_.reference_data.empty();
// Inherit the campaign's shared R-free test set from the reference MTZ (overriding the
// per-hkl hash that the merge assigned), so every dataset flags the same free reflections.
if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) {
const auto sg = experiment_.GetSpaceGroupNumber().value_or(1);
const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast<int32_t>(sg),
config_.reference_data);
logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)",
matched, sm.merged.size());
}
// Radiation-damage report (rotation): the per-image scale CC-to-merge and mosaicity across the
// sweep (the per-image scaling the rotation merge already fits, binned by frame = dose). A
// per-image CC that falls, and/or a mosaicity that rises, with frame number is the classic
// radiation-damage signature - a data-quality-vs-dose read complementary to the fitted decay
// correction. The full per-image table is written to <prefix>_scaling.txt for detail.
if (experiment_.IsRotationIndexing()) {
const auto &outs = indexer->GetIntegrationOutcome();
const int nf = static_cast<int>(outs.size());
constexpr int nb = 10;
std::array<double, nb> cc_sum{}, mos_sum{};
std::array<int, nb> cc_n{}, mos_n{};
for (int f = 0; f < nf; ++f) {
const int b = nf > 1 ? std::min(nb - 1, f * nb / nf) : 0;
if (outs[f].image_scale_cc && std::isfinite(*outs[f].image_scale_cc)) {
cc_sum[b] += *outs[f].image_scale_cc; cc_n[b]++;
}
if (outs[f].mosaicity_deg && std::isfinite(*outs[f].mosaicity_deg)) {
mos_sum[b] += *outs[f].mosaicity_deg; mos_n[b]++;
}
}
int tot = 0;
for (int n : cc_n) tot += n;
if (tot >= 20) {
std::ostringstream os;
os << fmt::format("Radiation-damage report (per-image scale over {} frames):\n", nf);
os << " dose <CC to merge> <mosaicity deg>\n";
for (int b = 0; b < nb; ++b) {
if (cc_n[b] == 0 && mos_n[b] == 0) continue;
os << fmt::format(" {:3d}-{:3d}% {:>7} {:>7}\n", b * 10, (b + 1) * 10,
cc_n[b] ? fmt::format("{:.3f}", cc_sum[b] / cc_n[b]) : std::string("-"),
mos_n[b] ? fmt::format("{:.3f}", mos_sum[b] / mos_n[b]) : std::string("-"));
}
if (cc_n[0] > 0 && cc_n[nb - 1] > 0) {
const double cc0 = cc_sum[0] / cc_n[0], cc9 = cc_sum[nb - 1] / cc_n[nb - 1];
os << fmt::format(" => per-image CC to merge {:.3f} (first 10%) -> {:.3f} (last 10%){}",
cc0, cc9,
cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : "");
}
// Relative B-factor across the sweep (measured before any decay correction, against the
// low-dose start): the resolution-dependent complement to the CC/mosaicity read above.
// Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose.
// A curve no straight line describes gets no first->last number at all, and a batch whose
// data could not measure one prints "-" rather than a value.
if (!sm.statistics.radiation_damage_b_batch.empty()) {
if (std::isfinite(sm.statistics.radiation_damage_delta_b))
os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}",
sm.statistics.radiation_damage_delta_b,
sm.statistics.radiation_damage_delta_b > 5.0
? " (significant dose-dependent scaling: radiation damage)" : "");
else
os << "\n => relative B-factor: the per-batch curve below is not a trend, so no "
"first->last number describes it - dose does not come back, so whatever moved "
"here was not dose; see the sweep-quality report";
os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):",
sm.statistics.radiation_damage_batch_deg);
for (float bb : sm.statistics.radiation_damage_b_batch)
os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -"));
}
logger.Info("{}", os.str());
result.radiation_damage_text = os.str();
}
// Sweep-quality report: the stretches of the sweep over which the crystal delivered much
// less than the rest of the run, and what each one looks like. Nothing is excluded because
// of it - the frames still carry signal, and this is a message for the beamline, not a
// filter. Frame numbers are processed-image ordinals, inclusive, as in <prefix>_image.dat.
const auto &sq = sm.statistics.sweep_quality;
if (sq.measured) {
std::ostringstream os;
os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, "
"which varied {:.2f}x, already divided out):\n",
sq.sweep_deg, sq.flux_peak_to_trough);
if (sq.ranges.empty()) {
os << " no stretch of the sweep is materially worse than the run";
} else {
os << " frames rotation diagnosis severity scale CC indexed\n";
for (const auto &r : sq.ranges)
os << fmt::format(" {:<17s} {:6.1f} deg {:<16s} {:5.2f} {:5.2f} {:5.2f} {:4.0f}%\n",
fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg,
SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale,
r.mean_relative_cc, 100.0 * r.indexed_fraction);
os << " => severity is the fraction of the run's typical diffracting power missing "
"over the range";
if (sq.modulation_peak_to_trough >= 1.05f)
os << fmt::format("\n => once-per-revolution modulation of the per-image scale: "
"{:.1f}x peak to trough", sq.modulation_peak_to_trough);
}
logger.Info("{}", os.str());
}
}
if (result.consensus_cell && write_files && config_.write_merged) {
phase("Writing reflections");
const ErrorModelReport em_report{
result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?",
result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic)
: std::string(),
result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(),
result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()};
WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics,
em_report, result.twinning, config_.output_prefix);
// Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS
// comparison. The offline self-scaling result is otherwise not exposed (process.h5's
// per-image arrays are only filled on the online per-image path). Sourced from the
// partials, which carry the first-pass per-image scale.
ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix);
}
if (result.consensus_cell && write_files && !config_.model_path.empty()) {
phase("Validating against model");
const auto data_sg = experiment_.GetSpaceGroupNumber();
// With a reference MTZ the merohedral indexing was already resolved against it (rotation
// merge / stills scaling), so trust that; only probe indexing by R-free when model-only.
ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path,
config_.output_prefix, logger,
data_sg ? std::optional<int>(static_cast<int>(*data_sg)) : std::nullopt,
/*probe_indexing_ambiguity=*/config_.reference_data.empty());
}
}
}
// When merging is on, the _process.h5 is skipped because the merged reflections are the wanted
// output and the per-image file is large. If nothing indexed there are no merged reflections, so
// the run would finish successfully having written no file at all. Write one now: the per-image
// messages went past while there was no writer, but the end message still carries the per-image
// scalars (spot counts, background, resolution estimate) and the dataset metadata - which is what
// anyone asking "why did nothing index" needs, and it is small precisely because nothing indexed.
if (!writer && full && write_files && !cancelled_ && result.indexing_rate.value_or(0.0f) <= 0.0f) {
logger.Warning("No image indexed, so there are no merged reflections to write - writing "
"{}_process.h5 with the per-image analysis instead", config_.output_prefix);
writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
/*trusted_path=*/true);
}
if (writer) {
writer->WriteHDF5(end_msg);
writer->Finalize();
result.written_master_path = config_.output_prefix + "_process.h5";
}
MarkPhase(observer, cancelled_ ? "Cancelled" : "Done");
// Total wall time for the whole run, including scaling/merging/writing above (not just the
// per-image pass), so the reported rate reflects everything rugnux did.
result.processing_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - start_time).count();
if (result.processing_time_s > 0.0) {
result.frame_rate_hz = static_cast<double>(result.images_processed) / result.processing_time_s;
result.throughput_MBs = static_cast<double>(total_uncompressed_bytes) / (result.processing_time_s * 1e6);
}
logger.Info("{} {} images in {:.2f} s ({:.2f} Hz)", cancelled_ ? "Cancelled after" : "Processed",
result.images_processed, result.processing_time_s, result.frame_rate_hz);
// The geometry this pass integrated at, for the report. Taken from experiment_ rather than from the
// caller's copy: the second pass of the rotation two-pass runs at the post-refined geometry.
result.used_beam_x_pxl = experiment_.GetBeamX_pxl();
result.used_beam_y_pxl = experiment_.GetBeamY_pxl();
result.used_distance_mm = experiment_.GetDetectorDistance_mm();
// Close the phase still running and hand over the breakdown. The timings accumulate in the object
// across passes, so the last pass's result carries the whole run rather than its own share of it.
ClosePhase();
result.phase_timings = phase_timings_;
return result;
}
Reference in New Issue View Git Blame Copy Permalink