mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases 67 Wiki Activity
Files
b5f5879a1d83c4c4336182dde7734cb5fdff4ee5
Jungfraujoch/rugnux/Rugnux.cpp
T
leonarski_fandClaude Opus 5 b5f5879a1d rugnux: measure the ice in the first pass, and always find its own spots 
Ice handling was gated on a measurement the run only made AFTER the images had
been processed, so the per-image pass could not use it. The flagging therefore
ran unconditionally: ice-band spots were ordered last in the --max-spots budget
and held out of the indexer seed and the geometry refinement on every crystal,
iced or not. The eleven bands are fixed geometry holding 16-26 % of the unique
reflections whether or not there is ice, so on a clean crystal that discards a
fifth of the spots - the strongest first - for nothing. Measured on a crystal
whose gate never fires, that moved the merged data by a mean of 0.85 sigma
against a run-to-run floor of 9.3e-5.

Measure it in the first pass instead. That pass already looks at ~100 images
spread over the sweep, and it already stops at the spot finder, so it sees the
azimuthal profile for the smooth channel and the unfiltered connected components
for the spot channel. Both counts SpotAnalyze takes are pre-filter, so pooling
them there is the run's own verdict, reached before anything has been discarded
and in time for the pass that acts on it. Where the sample sees no ice, the run
indexes on the ice-band spots too.

It has to be the whole sample: the spot channel is a ratio pooled over images,
because one frame carries a handful of control spots. A per-image gate is not an
alternative - two of the crystals whose indexing this rescues fire on that
channel alone, at profile scores of 1.12 and 1.22, so gating per image on the
profile score would drop exactly the cases that matter.

This also removes the first-pass spot reuse, and with it --redo-rotation-spots
and the reuse path. Finding the ~100 first-pass spots costs little, and reusing
was actively wrong here: the stored spots were found online at the acquisition's
threshold and have already had their ice-band entries ordered last and dropped
by its spot budget, so counting ice from them under-reads it by construction,
and the lattice search never saw the spot-finding settings at all. It also
removes the need for the machinery that re-found spots whenever a spot-finding
option was named, which made those options impossible to A/B.

IndexAndRefine cached index_ice_rings at construction, which happens before the
first pass; it holds a reference to the experiment, so it now reads the setting
where it uses it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-06 20:48:22 +02:00

1989 lines
126 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "Rugnux.h"
#include "ModelValidation.h"
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <functional>
#include <future>
#include <mutex>
#include <numeric>
#include <set>
#include <sstream>
#include "../reader/JFJochHDF5Reader.h"
#include "../common/JFJochMath.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/IndexAndRefine.h"
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/indexing/IndexerThreadPool.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/IceRingMask.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 {
// 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;
}
}
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)) {}
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());
if (observer)
observer->OnPhase("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);
}
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 std::string base_prefix = config_.output_prefix;
const auto gonio_snapshot = experiment_.GetGoniometer();
prepass_detector_geometry_.reset();
prepass_result_.reset();
force_rotation_result_.reset();
prepass_sg_reindexed_ = false;
prepass_merge_sg_.reset();
prepass_promoted_point_group_ = false;
logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry) -> {}_01_*",
base_prefix);
// Keep an empty prefix empty (the "compute stats, persist nothing" mode): the "_01" suffix would
// make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted.
config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01";
auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true);
if (cancelled_) { config_.output_prefix = base_prefix; 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]);
}
// 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 - EXCEPT a reindex-derived group, whose conventional
// setting the de-novo primitive frame does not share, which is left to re-search.
experiment_.SpaceGroupNumber(std::nullopt);
prepass_merge_sg_ = prepass_sg_reindexed_ ? std::nullopt : pass1.space_group_number;
prepass_promoted_point_group_ = pass1.twinning.laue_class_was_chosen_by_promotion;
logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
base_prefix);
config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix)
auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
// 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_;
pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
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]);
experiment_.SpaceGroupNumber(std::nullopt);
config_.output_prefix = base_prefix;
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;
return redo;
}
}
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;
}
}
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;
const int images_to_process = (end_image - start_image) / 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);
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);
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;
start_message.hdf5_source_data = reader_.GetHDF5DataSource(start_image, images_to_process);
std::unique_ptr<FileWriter> writer;
if (write_files && config_.write_process_h5)
writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
/*trusted_path=*/true);
logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]",
images_to_process, start_image, end_image, config_.stride, config_.nthreads,
full ? "full analysis" : "azimuthal integration");
if (observer)
observer->OnPhase(full ? "Full analysis" : "Azimuthal integration");
// Full-analysis shared engines.
std::unique_ptr<IndexerThreadPool> indexer_pool;
std::unique_ptr<IndexAndRefine> indexer;
if (full) {
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).
// Ice measured on the first-pass sample (see where they are filled): pooled here so the verdict
// below, which runs after the first pass and before the per-image loop, can read them.
std::vector<float> sample_ice_scores;
double sample_ice_ring_spots = 0.0, sample_ice_control_spots = 0.0;
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) {
if (observer)
observer->OnPhase("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). Accessed only
// from single-threaded sections (the scheme feed and the validation loop), so no locking.
std::map<int, std::vector<SpotToSave>> spot_cache;
// One analysis engine for the whole first pass. Constructing it allocates a CUDA stream and a
// full set of GPU buffers (preprocessing, spot finding, azimuthal integration, Bragg
// integration) - far more work than analysing a frame - so building it per image would pay for
// hundreds of them, serially. This section is single-threaded, so one engine is enough.
auto analysis = std::make_unique<MXAnalysisWithoutFPGA>(experiment_, mapping, pixel_mask_, *indexer,
/*enable_fused_adaptive_gpu=*/true);
auto profile = std::make_unique<AzimuthalIntegrationProfile>(mapping);
auto get_spots = [&](int ordinal) -> const std::vector<SpotToSave> & {
auto it = spot_cache.find(ordinal);
if (it != spot_cache.end())
return it->second;
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);
if (m.ice_ring_score)
sample_ice_scores.push_back(*m.ice_ring_score);
if (m.spot_count_ice_rings)
sample_ice_ring_spots += static_cast<double>(*m.spot_count_ice_rings);
if (m.spot_count_ice_control)
sample_ice_control_spots += static_cast<double>(*m.spot_count_ice_control);
spots = std::move(m.spots);
}
} catch (const std::exception &e) {
logger.Warning("First-pass spot read failed for image {}: {}", image_idx, e.what());
}
return spot_cache.emplace(ordinal, std::move(spots)).first->second;
};
// 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);
int count = 0;
for (const int ordinal : validation) {
DataMessage m{};
m.number = ordinal;
m.spots = get_spots(ordinal);
if (idx.IndexFrameOnly(m, validation_settings))
count++;
}
return count;
};
// Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to
// be indexed. Spots are pulled from the single-threaded cache here; the FFT + refinement runs
// separately (RunIndexing) so the two schemes' indexing can overlap.
auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector<int> &ordinals) {
auto ri = std::make_unique<RotationIndexer>(experiment_, pool);
for (const int ordinal : ordinals) {
if (cancelled_ || ri->AccumulationFull())
break;
ri->ProcessImage(ordinal, get_spots(ordinal), rot_angle(ordinal));
}
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;
const auto result = ris[i]->GetLattice();
if (!result.has_value())
continue;
const int score = count_indexed(idx, *result);
// Compare PRIMITIVE volumes: two schemes can find the same lattice in different
// settings, and a centred setting's cell is an exact integer multiple of its primitive
// one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive
// rhombohedral cell. Comparing the centred volumes makes the supercell test below fire
// on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which
// is enough to change the space group the search then picks.
const double vol = std::abs(
result->lattice.ToPrimitive(result->search_result.centering).CalcVolume());
const std::string &name = schemes[i].first;
logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score,
static_cast<int>(validation.size()));
// 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 = result;
}
}
return bp;
};
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;
}
}
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 --azint-only 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)
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.{}",
best.score, static_cast<int>(validation.size()), next_steps));
indexer->ForceRotationIndexerResult(*best.result);
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_);
}
}
}
}
// Ice-presence verdict, from the first-pass sample and before the pass that acts on it. The eleven
// hexagonal bands are fixed geometry and hold 16-26 % of the unique reflections whether or not the
// crystal has ice, so setting their spots aside on a clean crystal discards a fifth of them for
// nothing - and de-prioritising them in the --max-spots budget throws away the strongest first.
// Where the sample sees no ice, index on all the spots instead. Only the SAMPLE can decide this:
// the spot channel is a ratio pooled over images, so it is meaningless on any single frame.
if (full && !sample_ice_scores.empty() && experiment_.IsDetectIceRings()
&& !experiment_.GetIndexingSettings().GetIndexIceRings()) {
const float min_score = experiment_.GetScalingSettings().GetIceMinScore();
const float min_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio();
const float score = std::accumulate(sample_ice_scores.begin(), sample_ice_scores.end(), 0.0f)
/ static_cast<float>(sample_ice_scores.size());
const bool smooth = score >= min_score;
const bool spotty = min_ratio > 0.0f && sample_ice_control_spots > 0.0
&& sample_ice_ring_spots / sample_ice_control_spots >= min_ratio;
if (!smooth && !spotty) {
IndexingSettings is = experiment_.GetIndexingSettings();
is.IndexIceRings(true);
experiment_.ImportIndexingSettings(is);
logger.Info("First-pass sample of {} images sees no ice (score {:.2f} < {:.2f}, spot ratio "
"{:.2f} < {:.2f}): indexing on the ice-band spots too",
sample_ice_scores.size(), score, min_score,
sample_ice_control_spots > 0.0
? sample_ice_ring_spots / sample_ice_control_spots : NAN, min_ratio);
} else {
logger.Info("First-pass sample of {} images sees ice (score {:.2f}, spot ratio {:.2f}): "
"ice-band spots set aside for indexing",
sample_ice_scores.size(), score,
sample_ice_control_spots > 0.0
? sample_ice_ring_spots / sample_ice_control_spots : NAN);
}
}
// 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;
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) writer->Write(msg);
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();
plots.Add(msg, profile);
if (writer) writer->Write(msg);
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);
}
};
if (observer)
observer->OnPhase("Processing images");
std::function<void()> worker = full ? std::function<void()>(full_worker)
: std::function<void()>(azint_worker);
std::vector<std::future<void> > futures;
futures.reserve(config_.nthreads);
const auto image_loop_start = std::chrono::steady_clock::now();
for (int i = 0; i < config_.nthreads; ++i)
futures.push_back(std::async(std::launch::async, worker));
for (auto &f: futures)
f.get();
// 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();
// 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) {
if (observer) observer->OnPhase(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 whole mechanism - this flagging and the merge-time CC1/2 ring mask below -
// 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.
// Ice rings masked from the merge by the CC1/2 test after the final merge; empty until then,
// at which point a re-merge applies them.
std::vector<char> masked_ice_rings;
// 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_.spot_finding.ice_ring_width_Q_recipA,
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, masked_ice_rings);
result.error_model_isa = r.isa;
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);
merge_engine.MaskIceRings(masked_ice_rings, config_.spot_finding.ice_ring_width_Q_recipA);
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())
logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(),
merge_engine.ErrorModelB(),
merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0,
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);
result.error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
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);
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");
}
}
}
}
// 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, but only when a good low-res region
// exists above the cutoff. If even the lowest-res shell (b == 0) is below the
// threshold - as happens for very weak data whose I/sigma is floored by the
// error model although CC1/2 is still high - there is no clean/noise boundary,
// and cutting here would discard essentially all reflections and leave the
// re-merge with no resolution range. Leave d_min_search at 0 (no cut) instead.
if (b > 0) d_min_search = rs[lo].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;
// Set when the space group was determined via a change of basis (the intensity-based centred-lattice
// test): its conventional setting does not match the indexer's primitive frame, so the two-pass must
// NOT reuse it - the second pass re-searches (and re-runs the same test) instead.
bool sg_reindexed = false;
// 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.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), and take whichever search found MORE symmetry.
// The two disagree in one direction only: every crystal on which the Lorentz filter changes
// the answer LOSES an operator, because discarding 40-80% of the observations starves the
// operator correlations - it never invents one. So the higher of the two is the safe pick,
// and it costs nothing on crystals where the filter is irrelevant. 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 filtered merge is used ONLY to rescue a point group the full merge failed to
// confirm. Everything else - and 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";
}
if (order_of(sg_search) <= order_of(alt)) {
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);
sg_search = alt;
} else {
logger.Info("Space-group search: Lorentz-filtered merge supports {} where the "
"all-observation one supports {} - taking the higher symmetry",
sg_search.point_group_hm, alt.point_group_hm);
// Keep the point group the filter rescued, but re-decide the ABSENCES on the merge
// that still has all the observations in it. Screws and centering live in the weak
// reflections, which is precisely what a |zeta| cut throws away - reading them off
// the filtered merge is how four crystals lost their screw axes (P2_1 read as P2).
if (sg_search.point_group_representative.has_value()) {
SearchSpaceGroupOptions abs_opts = sg_opts;
abs_opts.fixed_point_group = *sg_search.point_group_representative;
sg_search = SearchSpaceGroup(sm_all.merged, abs_opts);
}
}
// Either way 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 exists
// only to answer the point-group question.
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) {
sg_reindexed = true;
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_.spot_finding.ice_ring_width_Q_recipA,
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;
phase("Re-scaling in space group " + sg.short_name());
sm = scale_and_merge(sg.short_name(), false);
}
result.space_group_search = sg_search;
} 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). EXCEPTION: a
// reindex-derived space group (the centred-lattice test) is in a conventional setting the indexer's
// primitive frame does not share, so reusing it would break the second pass's indexing - remember it so
// the caller clears it and lets the second pass re-search. Unlike a throwaway pre-pass, this one now
// CONTINUES to the write below so its own (header-geometry) result is saved as the first pass.
if (geometry_prepass)
prepass_sg_reindexed_ = sg_reindexed;
// 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_.spot_finding.ice_ring_width_Q_recipA,
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);
}
}
// Ice-ring CC1/2 mask: a hexagonal-ice ring whose merged half-set CC1/2 collapses well below its
// resolution shoulders has been decorrelated by ice (its Bragg is unrecoverable) - drop it and
// re-merge. Weak/absent rings track their neighbours and stay, so completeness on clean crystals
// is untouched. Uses the ring band (+/- ice width in q=2pi/d) vs the shoulders either side.
if (experiment_.IsDetectIceRings() && experiment_.GetScalingSettings().GetIceRingMergeMask()
&& ice_present && !sm.merged.empty()) {
auto mask = FindDecorrelatedIceRings(sm.merged, config_.spot_finding.ice_ring_width_Q_recipA,
logger);
if (!mask.empty()) {
masked_ice_rings = std::move(mask);
const auto final_sg = experiment_.GetGemmiSpaceGroup();
sm = scale_and_merge(final_sg ? final_sg->short_name() : "P1", 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";
} 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";
} 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;
{
// 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. A large
// magnitude flags radiation damage; positive is the typical direction (high-res fades with dose).
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,
std::fabs(sm.statistics.radiation_damage_delta_b) > 5.0
? " (significant dose-dependent scaling: radiation damage)" : "");
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 << fmt::format(" {:.1f}", bb);
}
logger.Info("{}", os.str());
}
}
if (result.consensus_cell && write_files && config_.write_merged) {
phase("Writing reflections");
WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics,
result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?",
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";
}
if (observer)
observer->OnPhase(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);
return result;
}
Reference in New Issue View Git Blame Copy Permalink