mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions 3 Packages Projects Releases 72 Wiki Activity
Files
49dd01331b85f8649f9964a065e6d5210e5ec791
Jungfraujoch/rugnux/rugnux_cli.cpp
T
leonarski_fandClaude Opus 5 49dd01331b calibration: let the spots vote on the geometry in ring mode too 
The rings and the spots fail in different regimes, which is the whole reason to
carry both. The header is right on a well-configured instrument and is the thing
a calibration is run to check. The profile's own estimates are exact while the
error stays small and stop meaning anything beyond that - past about half a ring
spacing each ring reaches the azimuthally averaged profile as two horns rather
than one peak, and the distance search reads a list of horns as a list of rings.
The circle through the spots reads nothing from the header at all: measured on a
110 mm LaB6 exposure, --calibration spots returns the same geometry from a header
400 px out in the centre AND eight times out in distance.

So ring mode now finds spots as well - half a second - and offers what they make
of the geometry as one more starting hypothesis, fitted like the others with the
residual left to choose. It is added whole rather than as a centre alone:
GuessGeometry votes for the circle centre, clusters the radii into rings and
takes the distance from the innermost, and those two belong together. Taking only
its centre would not have helped, because the distance candidates come from a
profile averaged about the header's centre, and a ring smeared over hundreds of
pixels cannot be un-smeared by reading its bins differently.

Measured on the 110 mm exposure, whose true PONI is 765.90 at 110.03 mm. A header
40 px out in the centre now lands within 0.6 px (it landed 41 px away before).
The cases with BOTH wrong, which failed before this and equally before the beam
centre work, now come out: 250 mm / 780 px gives 110.041 mm and 766.10 px against
42.3 mm and 782.7; 250 mm / 867 px gives 110.065 and 765.48. All five datasets
are unchanged from their correct headers and the distance still recovers from any
header between 25 and 1200 mm.

Past about a hundred pixels nothing rescues ring mode, and the reason is the
profile rather than the seeding: binned about a centre that far out it shows each
ring smeared across its own sectors, so even started from the exact answer the
fit leaves rms 4.8 px and drifts. Re-binning would fix it and would need the
images read a second time; --calibration spots, which never touches the profile,
already covers it.

That regime is now visible rather than silent. The spots' beam centre is printed
beside the fitted one as a cross-check - two methods sharing no assumption, so a
reader can see at a glance whether they agree. It costs nothing, the spots having
been found already, and it separates cleanly: 0.3-0.4 px on the good runs against
201.4 px on the exposure whose profile could not be fitted at all. Reported as a
fact and not gated on, since at long distance both methods weaken together and
the honest thing is to show the number (5.2 px at 300 mm, 12.6 at 500).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27
2026-08-31 17:58:26 +02:00

2668 lines
163 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <csignal>
#include <future>
#include <optional>
#include <getopt.h>
#include <iostream>
#include <limits>
#include <sstream>
#include <string>
#include <thread>
#include <type_traits>
#include <vector>
#include "../reader/JFJochHDF5Reader.h"
#include "../reader/JFJochCBFReader.h"
#include "../common/Logger.h"
#include "../common/Definitions.h"
#include "../common/DiffractionExperiment.h"
#include "../common/PixelMask.h"
#include "../common/print_license.h"
#include "../common/CUDAWrapper.h"
#include "../common/JFJochMath.h"
#include "../image_analysis/bragg_integration/CalcISigma.h"
#include "../image_analysis/geom_refinement/Calibrants.h"
#include "../image_analysis/indexing/IndexerThreadPool.h"
#include "../image_analysis/LoadFCalcFromMtz.h"
#include "../image_analysis/UpdateReflectionResolution.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/RfreeFlags.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "Rugnux.h"
#include "RugnuxDefaults.h"
#include "ModelValidation.h"
#include "ResultReport.h"
// Spots kept per image (the strongest ones) and handed to indexing. Offline reprocessing is not
// bound by the online spot budget, so this is rugnux's own default rather than the 250 the
// DatasetSettings constructor uses for the detector. Measured over the 37-crystal rotation battery,
// lifting it to 1000 leaves low-resolution R_meas better or equal on every crystal (16 better, 0
// worse, the rest untouched because their frames never reach the cap), with R_meas 14/4 and ISa 14/6
// in its favour and no measurable cost in wall clock. It is also what jfjoch_viewer already sends.
constexpr int64_t RUGNUX_MAX_SPOT_COUNT = 1000;
// Default background-ring high-side clip for broadband (non-zero bandwidth) data, tighter than the
// monochromatic default because a pink-beam reflection is long and its wings reach into the ring.
constexpr float BROADBAND_BACKGROUND_CLIP_NSIGMA = 3.0f;
// Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation.
constexpr double SMOOTH_G_DEFAULT_DEG = 5.0;
// Default rot3d per-batch relative-B batch width (bare --relative-b), in degrees of rotation.
constexpr double RELATIVE_B_DEFAULT_DEG = 10.0;
// What the run does (--mode). Each of these was once a flag of its own; one option keeps them
// mutually exclusive by construction.
enum class RugnuxMode { MX, Azint, Scale, Calibration };
// Azimuthal sectors used by --mode calibration --calibration rings when --azim-phi-bins was not given.
// The default of 1 is a plain radial profile, in which a ring has been averaged over every direction so
// nothing is left to say where its centre is (RingsFromAzimuthalProfile refuses below 4). Sweeping the
// count on a LaB6 exposure and on a crystal's ice rings, 32 and above agree to about 0.1 px and the
// residual is flat, while 8 sectors is visibly coarser (0.3 px away, twice the residual). Going higher
// buys nothing and runs into the 65534-bin limit on q x azimuth, so 32 is the default.
constexpr int CALIBRATION_AZIM_BINS_DEFAULT = 32;
void print_usage() {
std::cout << "Usage rugnux {<options>} <input.h5>" << std::endl;
std::cout << "Options:" << std::endl;
std::cout << " -o, --output-prefix <txt> Output file prefix (default: output). mx and scale runs always write "
"<prefix>_report.txt, the results report; an empty prefix writes nothing at all" << std::endl;
std::cout << " -N, --threads <num> Number of threads (default: all hardware threads, of which the "
"per-image loop uses at most 16 per GPU; an explicit value is used as given)" << std::endl;
std::cout << " -s, --start-image <num> Start image number (default: 0)" << std::endl;
std::cout << " -e, --end-image <num> End image number (default: all)" << std::endl;
std::cout << " -t, --stride <num> Image stride (default: 1)" << std::endl;
std::cout << " -v, --verbose Verbose output" << std::endl;
std::cout << std::endl;
std::cout << " Mode" << std::endl;
std::cout << " --mode <name> What this run does (default: mx)" << std::endl;
std::cout << " mx Full analysis - spot finding, indexing, integration and merging" << std::endl;
std::cout << " azint Only azimuthal integration (no spot finding/indexing); writes <prefix>_process.h5" << std::endl;
std::cout << " scale Only re-scale/merge the already-integrated reflections in <input> (no re-integration)" << std::endl;
std::cout << " calibration Determine the detector geometry from powder rings; writes <prefix>.poni" << std::endl;
std::cout << std::endl;
std::cout << " Calibration (--mode calibration)" << std::endl;
std::cout << " --calibrant <name> Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive). Overridden by -C: a unit cell given there IS the standard, with its rings enumerated from the cell and its absences from -S where that is given, or taken as primitive where it is not. -S supplies SYMMETRY absences only: a standard whose extinctions come from its atomic basis rather than from a symmetry element - silicon, whose 222 is absent in Fd-3m for that reason - is better named with --calibrant" << std::endl;
std::cout << " --calibration <txt> How the rings are measured: rings|spots (default: rings). rings sums the (q x azimuth) azimuthal profile over every processed image and fits the ring arcs in it; spots pools the found spots and fits those. -s/-e/--stride select the images; rings defaults --azim-phi-bins to 32" << std::endl;
std::cout << " --no-refine-tilt Do not refine the detector tilt: hold rot1/rot2 at the header value and fit only the beam centre and the distance. For a calibration handed to a program that cannot express a tilted detector, e.g. XDS" << std::endl;
std::cout << std::endl;
std::cout << " Detector mask" << std::endl;
std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl;
std::cout << " --estimate-beam-center Place the beam centre before anything is indexed, and use it in place of the header value when it is measured precisely enough. On a sweep that reaches half a turn it comes from the symmetry of the spot positions - the frames 180 deg apart are each other's mirror image, and every reflection is recorded twice - and where the sweep does not reach that far, from the isotropy of the scattered background, which needs only a few frames. It reads frames of its own, chosen as pairs half a turn apart, so it does not change the mask --detect-beam-stop finds. Ignored when --beam-x/--beam-y are given, and when a stills geometry refinement has already placed the centre from indexed spots" << std::endl;
std::cout << " --beam-center-search[=N|off] After a first pass that indexes fewer than half the validation frames, try the beam centre a pixel at a time out to N px along each detector axis and keep the first one that indexes a majority. A beam-centre error is fixed in the lab frame, so it smears the accumulated reciprocal-space cloud and the FFT takes an axis harmonic instead of the true axis; nothing downstream repairs that. ON by default (12 px); =off disables. It runs only after a pass that has already failed, so a run that indexes never pays for it" << std::endl;
std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. Nothing is committed on a run that indexes: the measured centre is adopted only where the file's indexes nothing and the measured one indexes a majority. ON by default; =off disables" << std::endl;
std::cout << " --no-fit-spindle Take the goniometer axis from the file rather than measuring the spindle's rotation about the beam from the spots. Measuring it is the DEFAULT and the file is never right: every master writes an exact lab axis and no goniometer is one. Both mirror lines of the beam-centre estimator turn with the spindle, so an axis a few tenths of a milliradian out smears the vote until a neighbouring tooth wins. Only has an effect with --estimate-beam-center" << std::endl;
std::cout << std::endl;
std::cout << " Spot finding" << std::endl;
std::cout << " --spot-sigma <num> Noise sigma level for spot finding (default: 4.0)" << std::endl;
std::cout << " --spot-threshold <num> Photon count threshold for spot finding (default: 10)" << std::endl;
std::cout << " --min-pix-per-spot <num> Minimum connected strong pixels per spot. If omitted, min-pix is chosen PER IMAGE (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed count x indexed fraction is kept. Give an explicit value to force a fixed min-pix instead." << std::endl;
std::cout << " --adaptive-spots Self-calibrating detection (DEFAULT): the strong-pixel threshold comes from each image's own per-resolution-ring noise instead of the fixed --spot-threshold, so one setting adapts across datasets (no per-dataset --spot-threshold/--spot-sigma tuning)." << std::endl;
std::cout << " --no-adaptive-spots Turn adaptive detection off and use the fixed --spot-threshold / --spot-sigma finder instead" << std::endl;
std::cout << " --spot-false-pixels <num> Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl;
std::cout << " --spot-high-resolution <num> High resolution limit for spot finding. If omitted (or 0), spot finding is not clipped in resolution and extends as far as the detector reaches" << std::endl;
std::cout << " --spot-low-resolution <num> Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data; 0 removes the limit)" << std::endl;
std::cout << " --max-spots <num> Max spot count per image, the strongest ones, handed to indexing. Default: measured from the data on rotation (as deep into each image's spot list as its spots still lie on the lattice), capped at 1000" << std::endl;
std::cout << " --detect-ice-rings[=on|off] Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling. Default: the master file's setting, or - where the file says nothing - on for rotation and off for stills" << std::endl;
std::cout << std::endl;
std::cout << " Indexing" << std::endl;
std::cout << " (A dataset with a rotation goniometer axis is processed as rotation data by default; use --force-still to override)" << std::endl;
std::cout << " --force-still Process a rotation (goniometer) dataset as independent stills (still indexing + per-image ScaleOnTheFly) instead of rotation" << std::endl;
std::cout << " -R, --two-pass-rotation[=num] Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default: 100)" << std::endl;
std::cout << " --single-pass-rotation[=num] Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl;
std::cout << " --min-indexed-spots <num> Spots a frame must have on the lattice before it counts as indexed (default: 9). The rotation first pass scores candidate lattices on that count, so on a pattern with fewer spots per frame no lattice can score and the run gives up; integration itself is not gated by it" << std::endl;
std::cout << " --force-rotation-lattice <vec> Force rotation indexer with external lattice (in Angstrom) : \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, skips first pass)" << std::endl;
std::cout << " --rotation-scale <k> Goniometer rotation scale: the stage turned k times the angle stored in the file (the commanded one). Applied to both passes; overrides the fitted correction" << std::endl;
std::cout << " --rotation-no-postrefine Disable the (default-on) two-pass rotation post-refine (post-refine detector distance/beam + cell/axis, then re-integrate; only the refined pass is written, as the canonical <prefix>_* output)" << std::endl;
std::cout << " -X, --indexing-algorithm <txt> Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl;
std::cout << " -S, --space-group <num|symbol> Space group number (92) or symbol (P43212) - for indexing and scaling" << std::endl;
std::cout << " -C, --unit-cell <cell> Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl;
std::cout << " --fft-min-unit-cell <num> Shortest cell axis the FFT search accepts, in A (default: 10). A candidate with a shorter axis is discarded, so a crystal below the default cannot be indexed at all; -C lowers it on its own to cover the cell given" << std::endl;
std::cout << " -r, --refine <txt> Geometry refinement algorithm (none|orientation|beam_and_lattice|flex); flex tries all three per image and keeps whichever indexes the most spots (alias: multi)" << std::endl;
std::cout << " --refine-geometry[=N|off] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200) then re-indexes (lifts weak-stills indexing). Default ON for stills when a reference cell is given (-C / reference MTZ); =off disables" << std::endl;
std::cout << " --index-ice-rings[=on|off] Index on the spots flagged as sitting on an ice ring too, instead of setting them aside (default: off; no effect without --detect-ice-rings)" << std::endl;
std::cout << std::endl;
std::cout << " Scaling and merging (on by default)" << std::endl;
std::cout << " --no-merge Skip scaling and merging; write only the per-image _process.h5" << std::endl;
std::cout << " --scale-fulls rot3d: after the 3D combine, refit a per-frame scale on the fulls (XDS order, Unity model). Default ON for rot3d" << std::endl;
std::cout << " --no-scale-fulls Disable the rot3d scale-fulls refit (it is on by default for rot3d)" << std::endl;
std::cout << " --write-process-h5 Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl;
std::cout << " --smooth-g[=deg] rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl;
std::cout << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl;
std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption + modulation correction surfaces fitted on the fulls after scale-fulls" << std::endl;
std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl;
std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl;
std::cout << " --scaling-high-resolution <num> High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl;
std::cout << " --scaling-low-resolution <num> Low resolution limit for scaling/merging, in A (default: 50, the value XDS configurations use; 0 = no limit). Reflections coarser than this sit behind or beside the beam stop and are measured on a background it has eaten into" << std::endl;
std::cout << " --resolution-cutoff <txt> Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl;
std::cout << " --resolution-cc-target <num> CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl;
std::cout << " --resolution-shells <num> Number of resolution shells in the reported statistics table (default: 9, as XDS)" << std::endl;
std::cout << " --ice-min-score <num> Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging and the exclusion from scaling (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl;
std::cout << " --ice-min-spot-ratio <num> Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl;
std::cout << " --min-partiality <num> Minimum partiality to accept reflection (default: 0.02)" << std::endl;
std::cout << " --capture-uncertainty <num> rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl;
std::cout << " --min-captured-fraction <num> rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl;
std::cout << " --mosaicity <num> Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
std::cout << " --prediction-mosaicity <num> Diagnostic: fix the rocking width the PREDICTION window opens to (deg), leaving the partiality on the per-image sigma_M. Prediction and partiality are one number by default, so a sigma_M that moves takes the integrated reflection population with it" << std::endl;
std::cout << " --reject-outliers <num> Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl;
std::cout << " --search-min-zeta <num> De-novo space-group search only: also search a merge of just the observations whose Lorentz geometry |zeta| reaches this, and report both answers (default: 0.85 for rotation, 0 = single search). Reflections crossing the Ewald sphere near-tangentially are measured worst and can make a real symmetry operator look like a twin law; where the two searches disagree, the merge of ALL the observations decides" << std::endl;
std::cout << " --min-image-cc <num> Per-image CC limit in percent (default: no limit)" << std::endl;
std::cout << " --scaling-iterations <num> Number of scaling iterations with no reference data (default: 3)" << std::endl;
std::cout << " -z, --reference-mtz <file> Reference MTZ file" << std::endl;
std::cout << " --reference-column <label> Reference MTZ column to use (default: auto - F-model, else IMEAN/I, else FP/FOBS/F)" << std::endl;
std::cout << " --model <file> After merging, validate vs this atomic model: R-free, 2Fo-Fc/Fo-Fc maps, and - where the merge kept the Bijvoet split - an anomalous difference map with the strongest anomalous scatterers named by the atom they sit on. PDB or mmCIF, gzipped or not, recognised by content rather than by extension" << std::endl;
std::cout << " --export-unmerged Write <prefix>_unmerged.mtz, an unmerged MTZ (POINTLESS column layout) of the integrated observations, for aimless / pointless / careless. On by default, whenever there is an output prefix. On a rotation run the partials of each reflection are summed into one full, written at the batch its rocking curve is centred on, with the summed rocking-curve fraction in FRACTIONCALC; an event that caught less of its rocking curve than --min-partiality is not written, as in the merge. Intensities carry the Lorentz-polarization factor and nothing else: the partiality is not divided out and the per-image scale is not applied at all, since those programs scale the data themselves. Written in --mode mx and --mode scale, and with --no-merge" << std::endl;
std::cout << " --no-export-unmerged Do not write <prefix>_unmerged.mtz. It is the largest file a run produces, so a run whose observations are not going to another scaling program - a regression battery, say - can skip writing it" << std::endl;
std::cout << " --export-unmerged-partials Write <prefix>_unmerged_partials.mtz, the same observations with each partial as its own row (one batch per image) flagged for the reading program to sum, instead of summed here. Off by default and independent of --export-unmerged; both can be written in one run" << std::endl;
std::cout << " --no-p1-crosscheck Do not write <prefix>_P1.mtz. That file is the same observations merged in P1 instead of the space group the run determined, so a wrong space group can be re-merged, re-solved or re-refined without processing the images again. Every rotation run that determines its own space group writes it - including one that determined P1, where it repeats the merged output - so it is there to be harvested whatever the run decided. Not written when -S fixed the space group: that run never integrated the group's centring absences, so its P1 merge would be missing whole classes of reflections. It is a full merge, not the degraded one the search runs on, and it costs one extra merge" << std::endl;
std::cout << std::endl;
std::cout << " Integration" << std::endl;
std::cout << " --bandwidth <num> Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
std::cout << " --integration-radius <r> Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
std::cout << " --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default: on for rotation, off for stills). r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl;
std::cout << " --integration-stencil <k> Push the r2..r3 background ring out by k times the beam's radial streak (bandwidth*Rpx), per reflection (default 0 = a fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially - the r1 signal box stays a circle. Needs --bandwidth: on a monochromatic beam the streak is zero and this does nothing" << std::endl;
std::cout << " --integration-high-resolution <num> High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl;
std::cout << " --max-hkl <n> Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl;
std::cout << " --background-clip <n> High-side clip of the background ring at mean + n*sqrt(mean) (default 4, or 3 when --bandwidth is set; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl;
std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl;
std::cout << " --background-trim <f> Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean. Applies whatever --bandwidth is set to" << std::endl;
std::cout << " --overlap <txt> What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl;
std::cout << " --overlap-minpk <f> Least fraction of a reflection's expected profile that must be usable for the reflection to be kept (default 0.75, XDS MINPK). Governs both ways part of a profile is lost: the fraction that must be READABLE - not masked, untrusted, in a detector gap or overloaded - in every profile mode, and, under --overlap reject, the fraction that must be cleanly the reflection's own. --integrator boxsum has no profile to renormalise with, so there a signal disk with any unreadable pixel is still discarded outright and the reject fraction is by disk AREA, which cuts harder" << std::endl;
std::cout << " --integrator <txt> Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl;
std::cout << " -q, --azim-q-spacing <num> Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;
std::cout << " --azim-min-q <num> Azimuthal-integration minimum Q (1/A)" << std::endl;
std::cout << " --azim-max-q <num> Azimuthal-integration maximum Q (1/A). If omitted, integration extends to the highest Q the detector reaches." << std::endl;
std::cout << " --azim-phi-bins <num> Number of azimuthal (phi) bins (default: 1)" << std::endl;
std::cout << " --polarization-correction <on|off> Enable/disable azimuthal polarization correction" << std::endl;
std::cout << " --solid-angle-correction <on|off> Enable/disable azimuthal solid angle correction" << std::endl;
std::cout << std::endl;
std::cout << " Geometry overrides (defaults taken from the input file)" << std::endl;
std::cout << " --beam-x <num> Beam center X (pixel)" << std::endl;
std::cout << " --beam-y <num> Beam center Y (pixel)" << std::endl;
std::cout << " --detector-distance <num> Detector distance (mm)" << std::endl;
std::cout << " --wavelength <num> Wavelength (A)" << std::endl;
std::cout << " --rot1 <num> PONI rotation 1 (rad)" << std::endl;
std::cout << " --rot2 <num> PONI rotation 2 (rad)" << std::endl;
std::cout << " --rot3 <num> PONI rotation 3, about the beam (rad)" << std::endl;
std::cout << " --detector-mirror-y Stored image is mirrored in Y vs the detector" << std::endl;
std::cout << " --detector-quarter-turns <0-3> Stored image is turned by this many 90 deg about the beam" << std::endl;
std::cout << " --polarization <num> Polarization factor" << std::endl;
}
enum {
OPT_SPOT_SIGMA = 1000,
OPT_SPOT_THRESHOLD,
OPT_MIN_PIX_PER_SPOT,
OPT_ADAPTIVE_SPOTS,
OPT_NO_ADAPTIVE_SPOTS,
OPT_SPOT_FALSE_PIXELS,
OPT_SPOT_RESOLUTION,
OPT_SPOT_LOW_RESOLUTION,
OPT_MAX_SPOTS,
OPT_MIN_PARTIALITY,
OPT_MIN_IMAGE_CC,
OPT_SEARCH_MIN_ZETA,
OPT_SCALING_ITERATIONS,
OPT_SCALING_HIGH_RESOLUTION,
OPT_SCALING_LOW_RESOLUTION,
OPT_RESOLUTION_CUTOFF,
OPT_RESOLUTION_CC_TARGET,
OPT_RESOLUTION_SHELLS,
OPT_SINGLE_PASS_ROTATION,
OPT_FORCE_ROTATION_LATTICE,
OPT_MIN_INDEXED_SPOTS,
OPT_ROTATION_NO_POSTREFINE,
OPT_ROTATION_SCALE,
OPT_BACKGROUND_CLIP,
OPT_BACKGROUND_RADIAL,
OPT_REFINE_GEOMETRY,
OPT_DETECT_BEAM_STOP,
OPT_BEAM_CENTER_SEARCH,
OPT_BEAM_CENTER_CHECK,
OPT_ESTIMATE_BEAM_CENTER,
OPT_FIT_SPINDLE,
OPT_NO_FIT_SPINDLE,
OPT_BANDWIDTH,
OPT_INTEGRATION_RADIUS,
OPT_ADAPTIVE_INTEGRATION_RADIUS,
OPT_INTEGRATION_STENCIL,
OPT_BACKGROUND_TRIM,
OPT_OVERLAP,
OPT_OVERLAP_MINPK,
OPT_MAX_HKL,
OPT_INTEGRATION_HIGH_RES,
OPT_REJECT_OUTLIERS,
OPT_REFERENCE_COLUMN,
OPT_MODEL,
OPT_EXPORT_UNMERGED,
OPT_NO_EXPORT_UNMERGED,
OPT_EXPORT_UNMERGED_PARTIALS,
OPT_DUMP_OBSERVATIONS,
OPT_INTEGRATOR,
OPT_SIMPLE_STILLS,
OPT_SCALE_FULLS,
OPT_CAPTURE_UNCERTAINTY,
OPT_MIN_CAPTURED_FRACTION,
OPT_MOSAICITY,
OPT_PREDICTION_MOSAICITY,
OPT_SMOOTH_G,
OPT_RELATIVE_B,
OPT_NO_SCALING_CORRECTIONS,
OPT_NO_EXPECTED_VARIANCE_MERGE,
OPT_DETECT_ICE_RINGS,
OPT_INDEX_ICE_RINGS,
OPT_ICE_MIN_SCORE,
OPT_ICE_MIN_SPOT_RATIO,
OPT_NO_SCALE_FULLS,
OPT_WRITE_PROCESS_H5,
OPT_FORCE_STILL,
OPT_AZIM_MIN_Q,
OPT_AZIM_MAX_Q,
OPT_AZIM_PHI_BINS,
OPT_MODE,
OPT_CALIBRANT,
OPT_CALIBRATION,
OPT_NO_MERGE,
OPT_NO_P1_CROSSCHECK,
OPT_POLARIZATION_CORRECTION,
OPT_SOLID_ANGLE_CORRECTION,
OPT_BEAM_X,
OPT_BEAM_Y,
OPT_DETECTOR_DISTANCE,
OPT_WAVELENGTH,
OPT_ROT1,
OPT_ROT2,
OPT_ROT3,
OPT_DETECTOR_MIRROR_Y,
OPT_DETECTOR_QUARTER_TURNS,
OPT_FFT_MIN_UNIT_CELL,
OPT_POLARIZATION,
OPT_NO_REFINE_TILT
};
static option long_options[] = {
{"verbose", no_argument, nullptr, 'v'},
{"output-prefix", required_argument, nullptr, 'o'},
{"threads", required_argument, nullptr, 'N'},
{"start-image", required_argument, nullptr, 's'},
{"end-image", required_argument, nullptr, 'e'},
{"stride", required_argument, nullptr, 't'},
{"indexing-algorithm", required_argument, nullptr, 'X'},
{"unit-cell", required_argument, nullptr, 'C'},
{"min-indexed-spots", required_argument, nullptr, OPT_MIN_INDEXED_SPOTS},
{"reference-mtz", required_argument, nullptr, 'z'},
{"reference-column", required_argument, nullptr, OPT_REFERENCE_COLUMN},
{"model", required_argument, nullptr, OPT_MODEL},
{"export-unmerged", no_argument, nullptr, OPT_EXPORT_UNMERGED},
{"no-export-unmerged", no_argument, nullptr, OPT_NO_EXPORT_UNMERGED},
{"export-unmerged-partials", no_argument, nullptr, OPT_EXPORT_UNMERGED_PARTIALS},
{"dump-observations", required_argument, nullptr, OPT_DUMP_OBSERVATIONS},
{"space-group", required_argument, nullptr, 'S'},
{"anomalous", no_argument, nullptr, 'A'},
{"mode", required_argument, nullptr, OPT_MODE},
{"calibrant", required_argument, nullptr, OPT_CALIBRANT},
{"calibration", required_argument, nullptr, OPT_CALIBRATION},
{"no-refine-tilt", no_argument, nullptr, OPT_NO_REFINE_TILT},
{"no-merge", no_argument, nullptr, OPT_NO_MERGE},
{"no-p1-crosscheck", no_argument, nullptr, OPT_NO_P1_CROSSCHECK},
{"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS},
{"no-scale-fulls", no_argument, nullptr, OPT_NO_SCALE_FULLS},
{"write-process-h5", no_argument, nullptr, OPT_WRITE_PROCESS_H5},
{"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G},
{"relative-b", optional_argument, nullptr, OPT_RELATIVE_B},
{"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS},
{"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE},
{"refine", required_argument, nullptr, 'r'},
{"two-pass-rotation", optional_argument, nullptr, 'R'},
{"single-pass-rotation", optional_argument, nullptr, OPT_SINGLE_PASS_ROTATION},
{"force-still", no_argument, nullptr, OPT_FORCE_STILL},
{"azim-q-spacing", required_argument, nullptr, 'q'},
{"azim-min-q", required_argument, nullptr, OPT_AZIM_MIN_Q},
{"azim-max-q", required_argument, nullptr, OPT_AZIM_MAX_Q},
{"azim-phi-bins", required_argument, nullptr, OPT_AZIM_PHI_BINS},
{"polarization-correction", required_argument, nullptr, OPT_POLARIZATION_CORRECTION},
{"solid-angle-correction", required_argument, nullptr, OPT_SOLID_ANGLE_CORRECTION},
{"beam-x", required_argument, nullptr, OPT_BEAM_X},
{"beam-y", required_argument, nullptr, OPT_BEAM_Y},
{"detector-distance", required_argument, nullptr, OPT_DETECTOR_DISTANCE},
{"wavelength", required_argument, nullptr, OPT_WAVELENGTH},
{"rot1", required_argument, nullptr, OPT_ROT1},
{"rot2", required_argument, nullptr, OPT_ROT2},
{"rot3", required_argument, nullptr, OPT_ROT3},
{"detector-mirror-y", no_argument, nullptr, OPT_DETECTOR_MIRROR_Y},
{"detector-quarter-turns", required_argument, nullptr, OPT_DETECTOR_QUARTER_TURNS},
{"fft-min-unit-cell", required_argument, nullptr, OPT_FFT_MIN_UNIT_CELL},
{"polarization", required_argument, nullptr, OPT_POLARIZATION},
{"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE},
{"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE},
{"rotation-scale", required_argument, nullptr, OPT_ROTATION_SCALE},
{"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY},
{"detect-beam-stop", optional_argument, nullptr, OPT_DETECT_BEAM_STOP},
{"beam-center-search", optional_argument, nullptr, OPT_BEAM_CENTER_SEARCH},
{"beam-center-check", optional_argument, nullptr, OPT_BEAM_CENTER_CHECK},
{"estimate-beam-center", no_argument, nullptr, OPT_ESTIMATE_BEAM_CENTER},
{"fit-spindle", no_argument, nullptr, OPT_FIT_SPINDLE},
{"no-fit-spindle", no_argument, nullptr, OPT_NO_FIT_SPINDLE},
{"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
{"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
{"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT},
{"adaptive-spots", no_argument, nullptr, OPT_ADAPTIVE_SPOTS},
{"no-adaptive-spots", no_argument, nullptr, OPT_NO_ADAPTIVE_SPOTS},
{"spot-false-pixels", required_argument, nullptr, OPT_SPOT_FALSE_PIXELS},
{"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION},
{"spot-low-resolution", required_argument, nullptr, OPT_SPOT_LOW_RESOLUTION},
{"max-spots", required_argument, nullptr, OPT_MAX_SPOTS},
{"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY},
{"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY},
{"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION},
{"mosaicity", required_argument, nullptr, OPT_MOSAICITY},
{"prediction-mosaicity", required_argument, nullptr, OPT_PREDICTION_MOSAICITY},
{"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC},
{"search-min-zeta", required_argument, nullptr, OPT_SEARCH_MIN_ZETA},
{"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
{"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
{"scaling-low-resolution", required_argument, nullptr, OPT_SCALING_LOW_RESOLUTION},
{"background-clip", required_argument, nullptr, OPT_BACKGROUND_CLIP},
{"background-radial", optional_argument, nullptr, OPT_BACKGROUND_RADIAL},
{"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF},
{"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET},
{"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS},
{"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
{"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
{"adaptive-integration-radius", optional_argument, nullptr, OPT_ADAPTIVE_INTEGRATION_RADIUS},
{"integration-stencil", required_argument, nullptr, OPT_INTEGRATION_STENCIL},
{"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
{"overlap", required_argument, nullptr, OPT_OVERLAP},
{"overlap-minpk", required_argument, nullptr, OPT_OVERLAP_MINPK},
{"max-hkl", required_argument, nullptr, OPT_MAX_HKL},
{"integration-high-resolution", required_argument, nullptr, OPT_INTEGRATION_HIGH_RES},
{"integrator", required_argument, nullptr, OPT_INTEGRATOR},
{"simple-stills", no_argument, nullptr, OPT_SIMPLE_STILLS},
{"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
{"index-ice-rings", optional_argument, nullptr, OPT_INDEX_ICE_RINGS},
{"ice-min-score", required_argument, nullptr, OPT_ICE_MIN_SCORE},
{"ice-min-spot-ratio", required_argument, nullptr, OPT_ICE_MIN_SPOT_RATIO},
{"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
{nullptr, 0, nullptr, 0}
};
void trim_in_place(std::string &t) {
size_t b = 0;
while (b < t.size() && std::isspace(static_cast<unsigned char>(t[b]))) b++;
size_t e = t.size();
while (e > b && std::isspace(static_cast<unsigned char>(t[e - 1]))) e--;
t = t.substr(b, e - b);
};
bool parse_float_strict(const std::string &t, float &out) {
try {
size_t idx = 0;
out = std::stof(t, &idx);
return idx == t.size();
} catch (...) {
return false;
}
};
// Parse a required numeric option argument, optionally bounded to [min_value, max_value], or print a
// clear error and exit. getopt hands option arguments over as raw C strings; atoi()/atof() silently
// return 0 on non-numeric input (so a typo like "--min-pix-per-spot 2O" becomes 2 or 0) and std::sto*
// throws, which would terminate the program. This rejects non-numeric input, trailing garbage
// ("1.5foo"), integer overflow, and out-of-range values. T may be integral or floating-point; the
// bounds default to the full representable range (i.e. unbounded).
template <typename T>
T parse_number_arg(const char *arg, const char *option_name, Logger &logger,
T min_value = std::numeric_limits<T>::lowest(),
T max_value = std::numeric_limits<T>::max()) {
std::string s = arg ? arg : "";
trim_in_place(s);
T value{};
bool parsed = false;
if (!s.empty()) {
try {
size_t idx = 0;
if constexpr (std::is_integral_v<T>) {
const long long v = std::stoll(s, &idx);
value = static_cast<T>(v);
parsed = (idx == s.size()) && (static_cast<long long>(value) == v); // no overflow
} else {
value = static_cast<T>(std::stod(s, &idx));
parsed = (idx == s.size());
}
} catch (...) {}
}
if (!parsed) {
logger.Error("Invalid numeric value for {}: '{}'", option_name, arg ? arg : "<null>");
print_usage();
exit(EXIT_FAILURE);
}
if (value < min_value || value > max_value) {
logger.Error("Value for {} out of range: {} (expected {} to {})",
option_name, value, min_value, max_value);
print_usage();
exit(EXIT_FAILURE);
}
return value;
}
// Thin unbounded wrappers for the existing floating-point call sites.
double parse_double_arg(const char *arg, const char *option_name, Logger &logger) {
return parse_number_arg<double>(arg, option_name, logger);
}
float parse_float_arg(const char *arg, const char *option_name, Logger &logger) {
return parse_number_arg<float>(arg, option_name, logger);
}
bool parse_on_off(const char *arg, bool &out) {
std::string s = arg ? arg : "";
std::transform(s.begin(), s.end(), s.begin(),
[](unsigned char c) { return static_cast<char>(std::tolower(c)); });
if (s == "on" || s == "1" || s == "true" || s == "yes") {
out = true;
return true;
}
if (s == "off" || s == "0" || s == "false" || s == "no") {
out = false;
return true;
}
return false;
}
std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
if (!arg)
return std::nullopt;
std::string s(arg);
trim_in_place(s);
if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
s = s.substr(1, s.size() - 2);
trim_in_place(s);
}
std::vector<std::string> parts;
parts.reserve(6);
size_t start = 0;
while (true) {
size_t pos = s.find(',', start);
if (pos == std::string::npos) {
parts.push_back(s.substr(start));
break;
}
parts.push_back(s.substr(start, pos - start));
start = pos + 1;
}
if (parts.size() != 6)
return std::nullopt;
UnitCell uc{};
if (!parse_float_strict(parts[0], uc.a)) return std::nullopt;
if (!parse_float_strict(parts[1], uc.b)) return std::nullopt;
if (!parse_float_strict(parts[2], uc.c)) return std::nullopt;
if (!parse_float_strict(parts[3], uc.alpha)) return std::nullopt;
if (!parse_float_strict(parts[4], uc.beta)) return std::nullopt;
if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt;
return uc;
}
std::optional<CrystalLattice> parse_lattice_arg(const char *arg) {
if (!arg)
return std::nullopt;
std::string s(arg);
trim_in_place(s);
if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
s = s.substr(1, s.size() - 2);
trim_in_place(s);
}
std::vector<std::string> parts;
parts.reserve(9);
size_t start = 0;
while (true) {
size_t pos = s.find(',', start);
if (pos == std::string::npos) {
parts.push_back(s.substr(start));
break;
}
parts.push_back(s.substr(start, pos - start));
start = pos + 1;
}
if (parts.size() != 9)
return std::nullopt;
std::vector<float> vals(9);
for (int i = 0; i < 9; i++) {
if (!parse_float_strict(parts[i], vals[i]))
return std::nullopt;
}
return CrystalLattice(vals);
}
// Shared offline-output settings for the _process.h5 / reflection writer, used by both the --scale
// path and the full-analysis path (each then sets its own space group and images-per-trigger).
void configure_offline_output(DiffractionExperiment &experiment, const std::string &output_prefix) {
ApplyRugnuxExperimentDefaults(experiment); // analysis policy shared with the viewer
experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4);
// Offline CLI: the operator chose the output path, so allow an absolute -o (the multi-user guard
// that FilePrefix() applies is only for remotely-supplied prefixes in the broker/writer).
experiment.FilePrefixTrusted(output_prefix);
experiment.Mode(DetectorMode::Standard); // full image analysis
experiment.PixelSigned(true);
experiment.OverwriteExistingFiles(true);
experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
experiment.NumTriggers(1);
}
namespace {
std::atomic<Rugnux *> g_active_process{nullptr};
void handle_sigint(int) {
if (auto *p = g_active_process.load())
p->Cancel();
}
// argv as one line, echoed at the top of the run and recorded in the results report. An argument
// that would not survive being pasted back into a shell is single-quoted, so a file prefix with a
// space in it comes back as the one argument it was.
std::string JoinCommandLine(int argc, char **argv) {
std::string out;
for (int i = 0; i < argc; i++) {
if (i > 0)
out += ' ';
const std::string arg = argv[i];
if (arg.find_first_of(" \t\n'\"\\$`*?()[]{}<>|&;#~") == std::string::npos) {
out += arg;
continue;
}
out += '\'';
for (char c : arg) {
if (c == '\'')
out += "'\\''";
else
out += c;
}
out += '\'';
}
return out;
}
}
// The body of main. Settings setters and the pipeline itself throw JFJochException on input the parser
// cannot reject on its own (a resolution limit of zero, a distance of zero, a polarisation above one,
// an unreadable file), so main wraps this and reports rather than letting the exception terminate the
// process with no diagnostic and exit code 134.
static int RunRugnux(int argc, char **argv) {
// The whole invocation, reported at the end and written into the results report. Started before
// anything else so it covers reading the file and setting up as well as the processing.
const auto invocation_start = std::chrono::steady_clock::now();
RunProvenance provenance;
provenance.command_line = JoinCommandLine(argc, argv);
provenance.gpu_count = get_gpu_count();
provenance.gpu_description = get_gpu_description();
std::cout << provenance.command_line << std::endl << std::endl;
// Said before the run rather than after it: the GPUs are what makes rugnux fast, and a machine
// that turns out to have none - a driver mismatch, a CUDA_VISIBLE_DEVICES left over from another
// job - is worth knowing about while there is still time to stop and fix it.
if (provenance.gpu_count > 0)
std::cout << fmt::format("GPU: {} ({} visible)",
provenance.gpu_description, provenance.gpu_count) << std::endl;
else
std::cout << "GPU: none visible - running on the CPU, which is much slower"
<< std::endl;
RegisterHDF5Filter();
print_license("rugnux");
Logger logger("rugnux");
std::string input_file;
std::string output_prefix = "output";
int nthreads = 0; // 0 = auto: resolved to all hardware threads after parsing (see below)
int start_image = 0;
int end_image = -1; // -1 indicates process until end
int image_stride = 1;
bool verbose = false;
RugnuxMode mode = RugnuxMode::MX; // --mode
std::string calibrant = "lab6"; // --calibrant
CalibrationMethod calibration_method = CalibrationMethod::Rings; // --calibration
bool calibration_refine_tilt = true; // --no-refine-tilt
bool rotation_indexing = false;
bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling)
bool two_pass_rotation = true;
// Set by any spot-finding option. The two-pass rotation first pass reuses the spots stored in the
// file when it has them, so those options would otherwise not reach the pass that determines the
// lattice - the setting would appear to do nothing at all on rotation data.
bool rotation_postrefine_geometry = true; // default on; --rotation-no-postrefine disables it
std::optional<float> rotation_scale; // --rotation-scale: asserted stage calibration
int rotation_indexing_image_count = 100;
std::optional<float> rotation_indexing_range;
std::optional<int64_t> min_indexed_spots; // --min-indexed-spots: spots on the lattice for a frame to count as indexed
bool run_scaling = true; // merge is on by default; --no-merge turns it off
bool write_p1_crosscheck = true; // <prefix>_P1.mtz is written by default; --no-p1-crosscheck declines it
std::optional<bool> scale_fulls_arg; // --scale-fulls / --no-scale-fulls; default on for rot3d
bool write_process_h5_flag = false; // --write-process-h5; also write _process.h5 when merging
std::optional<bool> detect_ice_rings; // --detect-ice-rings[=on|off]; unset => use the dataset (file) value
bool index_ice_rings = false; // --index-ice-rings[=on|off]; index on the ice-band spots too
std::optional<double> ice_min_score_arg; // --ice-min-score: ice-presence gate on the measured score
std::optional<double> ice_min_spot_ratio_arg; // --ice-min-spot-ratio: the same gate on the spot channel
std::optional<float> fft_min_unit_cell_A; // --fft-min-unit-cell: shortest axis the FFT search accepts
std::optional<float> min_q, max_q, q_spacing; // azimuthal integration range / -q spacing (1/A)
std::optional<int32_t> azimuthal_bins; // --azimuthal-bins
std::optional<bool> polarization_correction; // --polarization-correction (azimuthal integration)
std::optional<bool> solid_angle_correction; // --solid-angle-correction (azimuthal integration)
// Geometry overrides (default: keep the value stored in the input file)
std::optional<float> beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor;
bool detector_mirror_y = false;
int64_t detector_quarter_turns = 0;
std::optional<double> smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise
std::optional<double> relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given
bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces
bool no_expected_variance_merge = false; // --no-expected-variance-merge: restore observed-sigma stills merge weighting
bool anomalous_mode = false;
const gemmi::SpaceGroup *space_group = nullptr;
std::optional<UnitCell> fixed_reference_unit_cell;
std::optional<int64_t> max_spot_count_override;
float sigma_spot_finding = 4.0;
int64_t photon_count_threshold_spot_finding = 10;
std::optional<int64_t> min_pix_per_spot; // unset -> adaptive per image; a value -> fixed min-pix
std::optional<bool> adaptive_spots; // unset -> on, for both workflows
float false_pixels_per_frame = 100.0f;
std::string ref_mtz;
std::string ref_column;
std::string model_pdb; // --model: PDB to validate merged intensities against (R-free + maps)
bool export_unmerged = true; // on by default; --no-export-unmerged turns it off
bool export_unmerged_partials = false; // --export-unmerged-partials: the same, one row per image
std::string dump_observations; // diagnostic: dump unmerged -P rot3d fulls to this path
double min_partiality = 0.02;
std::optional<double> min_captured_fraction_arg; // explicit --min-captured-fraction; default depends on rotation
std::optional<double> capture_uncertainty_arg; // explicit --capture-uncertainty; default depends on rot3d
std::optional<double> forced_mosaicity_arg; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed
std::optional<double> forced_prediction_mosaicity_arg; // diagnostic: fix the PREDICTION mosaicity (deg), leaving the partiality on the per-image value
double min_image_cc = 0.0;
std::optional<double> search_min_zeta_arg; // --search-min-zeta; rotation default below
int64_t scaling_iter = 3;
std::optional<CrystalLattice> forced_rotation_lattice;
std::optional<double> background_clip_arg; // --background-clip: background-ring high-side sigma clip
bool background_radial_given = false; // --background-radial seen at all (unset => auto)
std::optional<bool> background_radial_arg; // when given: set = force on/off, unset = auto
std::optional<int> detect_beam_stop = 60; // --detect-beam-stop[=N|off]; on by default
bool estimate_beam_center = false; // --estimate-beam-center
int beam_center_search = 12; // --beam-center-search[=N|off]; on by default
bool beam_center_check = true; // --beam-center-check[=off]; on by default
bool fit_spindle = true; // --fit-spindle / --no-fit-spindle
std::optional<int> refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass
bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on
std::optional<float> bandwidth_fwhm; // relative FWHM of dlambda/lambda
IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
std::optional<float> d_min_spot_finding; // unset -> as far as the detector reaches
// The low-resolution limit to apply: starts at the SpotFindingSettings default rather than restating
// it, and --spot-low-resolution 0 resets it to unset, i.e. no limit at that end.
std::optional<float> d_max_spot_finding = SpotFindingSettings{}.low_resolution_limit;
std::optional<float> d_min_scale_merge;
std::optional<float> d_max_scale_merge; // --scaling-low-resolution; 0 removes the default limit
std::optional<ResolutionCutoffMethod> resolution_cutoff_method; // --resolution-cutoff cc-logistic|off
std::optional<double> resolution_cc_target; // --resolution-cc-target
std::optional<int> report_shell_count; // --resolution-shells
std::optional<std::string> integration_radius_arg;
// --adaptive-integration-radius: unset means the per-workflow default below (on for rotation).
std::optional<bool> adaptive_integration_radius_arg;
std::optional<std::string> integration_stencil_arg; // --integration-stencil: ring elongation, in sigma
std::optional<double> background_trim_arg; // --background-trim: background-ring trimmed-mean fraction
std::optional<OverlapMode> overlap_arg; // --overlap: treatment of shared signal pixels
std::optional<double> overlap_minpk_arg; // --overlap-minpk: XDS-like MINPK threshold
std::optional<int64_t> max_hkl_arg; // --max-hkl: half-width of the predicted hkl box
std::optional<double> integration_d_min_arg; // --integration-high-resolution; unset = detector reach
std::optional<IntegratorMode> integrator_mode; // --integrator boxsum|gaussian|empirical
bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement
std::optional<double> outlier_reject_nsigma; // merge per-observation outlier rejection
if (argc == 1) {
print_usage();
exit(EXIT_FAILURE);
}
int opt;
int option_index = 0;
const char *short_opts = "vo:N:s:e:t:R::X:C:z:FAS:r:q:";
while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) {
switch (opt) {
case 'o':
output_prefix = optarg;
break;
case 'v':
verbose = true;
break;
case 'N':
nthreads = atoi(optarg);
break;
case 's':
start_image = atoi(optarg);
break;
case 'e':
end_image = atoi(optarg);
break;
case 't':
image_stride = atoi(optarg);
break;
case 'R':
if (rotation_indexing) {
logger.Error("Rotation indexing already enabled");
exit(EXIT_FAILURE);
}
rotation_indexing = true;
two_pass_rotation = true;
if (optarg)
rotation_indexing_image_count = atoi(optarg);
break;
case OPT_SINGLE_PASS_ROTATION:
if (rotation_indexing) {
logger.Error("Rotation indexing already enabled");
exit(EXIT_FAILURE);
}
rotation_indexing = true;
two_pass_rotation = false;
if (optarg)
rotation_indexing_range = atof(optarg);
break;
case OPT_ROTATION_NO_POSTREFINE:
rotation_postrefine_geometry = false;
break;
case OPT_ROTATION_SCALE:
rotation_scale = parse_number_arg<float>(optarg, "--rotation-scale", logger, 0.9f, 1.1f);
break;
case OPT_MIN_INDEXED_SPOTS:
min_indexed_spots = parse_number_arg<int64_t>(optarg, "--min-indexed-spots", logger);
break;
case OPT_DETECT_BEAM_STOP:
// Frames projected to find the shadow. The default is what the detection was validated
// on; fewer leaves the background too sparsely counted to tell a shadow from noise.
if (optarg && std::string(optarg) == "off") {
detect_beam_stop = std::nullopt;
break;
}
detect_beam_stop = optarg
? parse_number_arg<int>(optarg, "--detect-beam-stop", logger, 1, 1000000)
: 60;
break;
case OPT_BEAM_CENTER_SEARCH:
// How far the ladder reaches, in pixels, in each direction. The step is one pixel, flat:
// a header beam centre is written in pixels and is wrong by pixels. Deriving the step from
// the J0 law was tried and is wrong, because the only cell available at that point is the
// one the FAILED pass returned - on a dataset whose failed cell was a small spurious
// sub-cell the formula asked for a 6 px step, which steps clean over the lobe it is
// looking for.
if (optarg && std::string(optarg) == "off") {
beam_center_search = 0;
break;
}
beam_center_search = optarg
? parse_number_arg<int>(optarg, "--beam-center-search", logger, 1, 1000)
: 12;
break;
case OPT_BEAM_CENTER_CHECK:
beam_center_check = !optarg || std::string(optarg) != "off";
break;
case OPT_ESTIMATE_BEAM_CENTER:
estimate_beam_center = true;
break;
case OPT_FIT_SPINDLE:
fit_spindle = true;
estimate_beam_center = true;
break;
case OPT_NO_FIT_SPINDLE:
fit_spindle = false;
break;
case OPT_REFINE_GEOMETRY: {
if (optarg && std::string(optarg) == "off") {
refine_geometry = std::nullopt;
refine_geometry_disabled = true; // opt out of the stills default-on
break;
}
// Positive frame count fed to the bundle adjust; upper-bounded so refine_frames * 50 in
// RefineStillsGeometry cannot overflow int (1e6 is already far past any real dataset).
refine_geometry = optarg
? parse_number_arg<int>(optarg, "--refine-geometry", logger, 1, 1000000)
: 200;
break;
}
case OPT_FORCE_ROTATION_LATTICE: {
if (rotation_indexing) {
logger.Error("Rotation indexing already enabled");
exit(EXIT_FAILURE);
}
rotation_indexing = true;
auto latt = parse_lattice_arg(optarg);
if (!latt.has_value()) {
logger.Error(
"Invalid rotation lattice. Expected: \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, comma-separated). Got: {}",
optarg ? optarg : "<null>");
print_usage();
exit(EXIT_FAILURE);
}
forced_rotation_lattice = latt;
auto uc = latt->GetUnitCell();
logger.Info(
"Forced rotation lattice set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
break;
}
case 'X': {
std::string alg = optarg ? optarg : "";
std::transform(alg.begin(), alg.end(), alg.begin(),
[](unsigned char c) { return static_cast<char>(std::tolower(c)); });
if (alg == "ffbidx")
indexing_algorithm = IndexingAlgorithmEnum::FFBIDX;
else if (alg == "fft")
indexing_algorithm = IndexingAlgorithmEnum::FFT;
else if (alg == "fftw")
indexing_algorithm = IndexingAlgorithmEnum::FFTW;
else if (alg == "auto")
indexing_algorithm = IndexingAlgorithmEnum::Auto;
else if (alg == "none")
indexing_algorithm = IndexingAlgorithmEnum::None;
else {
logger.Error("Invalid indexing algorithm: {}", alg);
print_usage();
exit(EXIT_FAILURE);
}
break;
}
case 'r': {
std::string alg = optarg ? optarg : "";
std::transform(alg.begin(), alg.end(), alg.begin(),
[](unsigned char c) { return static_cast<char>(std::tolower(c)); });
if (alg == "none")
refinement_algorithm = GeomRefinementAlgorithmEnum::None;
else if (alg == "beam_and_lattice")
refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
else if (alg == "orientation")
refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly;
else if (alg == "flex" || alg == "multi") // "multi" kept as a back-compat alias
refinement_algorithm = GeomRefinementAlgorithmEnum::Flex;
else {
logger.Error("Invalid geom refinement algorithm: {}", alg);
print_usage();
exit(EXIT_FAILURE);
}
break;
}
case 'C': {
auto uc = parse_unit_cell_arg(optarg);
if (!uc.has_value()) {
logger.Error(
"Invalid unit cell. Expected: \"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}",
optarg ? optarg : "<null>");
print_usage();
exit(EXIT_FAILURE);
}
fixed_reference_unit_cell = uc;
logger.Info(
"Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma);
break;
}
case 'z':
ref_mtz = optarg;
break;
case OPT_REFERENCE_COLUMN:
ref_column = optarg;
break;
case OPT_MODEL:
model_pdb = optarg;
break;
case OPT_EXPORT_UNMERGED:
export_unmerged = true;
break;
case OPT_NO_EXPORT_UNMERGED:
export_unmerged = false;
break;
case OPT_EXPORT_UNMERGED_PARTIALS:
export_unmerged_partials = true;
break;
case OPT_DUMP_OBSERVATIONS:
dump_observations = optarg;
break;
case 'F':
indexing_algorithm = IndexingAlgorithmEnum::FFT;
break;
case 'A':
anomalous_mode = true;
break;
case 'S': {
// Accept a space-group number ("92") or a Hermann-Mauguin symbol ("P43212", "P 43 21 2").
// The symbol is kept as the group it names, not reduced to its number: a number can only
// ever name the reference setting, so -S P1121 would silently become P 1 21 1.
char *end = nullptr;
const long as_number = strtol(optarg, &end, 10);
space_group = (end != optarg && *end == '\0')
? gemmi::find_spacegroup_by_number(static_cast<int>(as_number))
: gemmi::find_spacegroup_by_name(optarg);
if (space_group == nullptr) {
logger.Error("Unknown space group '{}' (use a number like 92 or a symbol like P43212)", optarg);
exit(EXIT_FAILURE);
}
break;
}
case OPT_SPOT_SIGMA:
sigma_spot_finding = parse_number_arg<float>(optarg, "--spot-sigma", logger, 1.0f);
logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding);
break;
case OPT_SPOT_THRESHOLD:
photon_count_threshold_spot_finding = parse_number_arg<int64_t>(optarg, "--spot-threshold", logger, 0);
logger.Info("Photon-count threshold level for spot finding set to {:d}",
photon_count_threshold_spot_finding);
break;
case OPT_MIN_PIX_PER_SPOT:
// Giving an explicit min-pix opts out of the per-image adaptive selection.
min_pix_per_spot = parse_number_arg<int64_t>(optarg, "--min-pix-per-spot", logger, 1);
logger.Info("Minimum pixels per spot fixed at {:d} (adaptive per-image min-pix off)", *min_pix_per_spot);
break;
case OPT_ADAPTIVE_SPOTS:
adaptive_spots = true;
logger.Info("Adaptive (self-calibrating) spot detection enabled");
break;
case OPT_NO_ADAPTIVE_SPOTS:
adaptive_spots = false;
logger.Info("Adaptive spot detection off: using the fixed --spot-threshold / --spot-sigma finder");
break;
case OPT_SPOT_FALSE_PIXELS:
false_pixels_per_frame = parse_number_arg<float>(optarg, "--spot-false-pixels", logger, 1.0f);
adaptive_spots = true;
logger.Info("Adaptive spot detection: expected false pixels/frame set to {:.0f}", false_pixels_per_frame);
break;
case OPT_SPOT_LOW_RESOLUTION: {
// Same shape as --spot-high-resolution below: 0 means "no limit at that end", expressed as
// the unset optional the rest of the code understands.
const auto d_max = parse_number_arg<float>(optarg, "--spot-low-resolution", logger, 0.0f);
if (d_max > 0.0f) {
d_max_spot_finding = d_max;
logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max);
} else {
d_max_spot_finding.reset();
logger.Info("No low resolution limit for spot finding");
}
break;
}
case OPT_SPOT_RESOLUTION: {
// 0 has always meant "no limit" for this setting; keep that, but express it as the unset
// optional the rest of the code understands. Passing the 0 through instead reached
// ResolutionShells (via the spot plot), which rejects a zero d_min and threw away every image.
const auto d_min = parse_number_arg<float>(optarg, "--spot-high-resolution", logger, 0.0f);
if (d_min > 0.0f) {
d_min_spot_finding = d_min;
logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min);
} else {
d_min_spot_finding.reset();
logger.Info("No high resolution limit for spot finding: as far as the detector reaches");
}
break;
}
case OPT_MAX_SPOTS:
max_spot_count_override = parse_number_arg<int64_t>(optarg, "--max-spots", logger, 1);
break;
case OPT_MODE: {
const std::string m = optarg ? optarg : "";
if (m == "mx") mode = RugnuxMode::MX;
else if (m == "azint") mode = RugnuxMode::Azint;
else if (m == "scale") mode = RugnuxMode::Scale;
else if (m == "calibration") mode = RugnuxMode::Calibration;
else {
logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration)", m);
return 1;
}
break;
}
case OPT_CALIBRANT:
calibrant = optarg ? optarg : "";
if (CalibrantRings(calibrant).empty()) {
logger.Error("Unknown --calibrant: {} (expected {})", calibrant, CalibrantNameList());
return 1;
}
break;
case OPT_CALIBRATION:
if (strcmp(optarg, "rings") == 0) calibration_method = CalibrationMethod::Rings;
else if (strcmp(optarg, "spots") == 0) calibration_method = CalibrationMethod::Spots;
else {
logger.Error("Invalid --calibration: {} (expected rings|spots)", optarg);
return 1;
}
break;
case OPT_NO_REFINE_TILT:
calibration_refine_tilt = false;
break;
case OPT_NO_MERGE:
run_scaling = false;
break;
case OPT_NO_P1_CROSSCHECK:
write_p1_crosscheck = false;
break;
case OPT_SCALE_FULLS:
scale_fulls_arg = true;
break;
case OPT_NO_SCALE_FULLS:
scale_fulls_arg = false;
break;
case OPT_DETECT_ICE_RINGS:
if (optarg == nullptr || strcmp(optarg, "on") == 0)
detect_ice_rings = true;
else if (strcmp(optarg, "off") == 0)
detect_ice_rings = false;
else {
logger.Error("Invalid --detect-ice-rings value: {} (expected on|off)", optarg);
print_usage();
exit(EXIT_FAILURE);
}
break;
case OPT_INDEX_ICE_RINGS:
if (optarg == nullptr || strcmp(optarg, "on") == 0)
index_ice_rings = true;
else if (strcmp(optarg, "off") == 0)
index_ice_rings = false;
else {
logger.Error("Invalid --index-ice-rings value: {} (expected on|off)", optarg);
print_usage();
exit(EXIT_FAILURE);
}
break;
case OPT_ICE_MIN_SCORE:
ice_min_score_arg = parse_double_arg(optarg, "--ice-min-score", logger);
break;
case OPT_ICE_MIN_SPOT_RATIO:
ice_min_spot_ratio_arg = parse_double_arg(optarg, "--ice-min-spot-ratio", logger);
break;
case OPT_WRITE_PROCESS_H5:
write_process_h5_flag = true;
break;
case OPT_SMOOTH_G:
smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG;
break;
case OPT_RELATIVE_B:
relative_b_deg_arg = optarg ? parse_double_arg(optarg, "--relative-b", logger) : RELATIVE_B_DEFAULT_DEG;
break;
case OPT_NO_EXPECTED_VARIANCE_MERGE:
no_expected_variance_merge = true;
break;
case OPT_NO_SCALING_CORRECTIONS:
no_scaling_corrections = true;
break;
case OPT_MIN_PARTIALITY:
min_partiality = parse_double_arg(optarg, "--min-partiality", logger);
break;
case OPT_CAPTURE_UNCERTAINTY:
capture_uncertainty_arg = parse_double_arg(optarg, "--capture-uncertainty", logger);
break;
case OPT_MIN_CAPTURED_FRACTION:
min_captured_fraction_arg = parse_double_arg(optarg, "--min-captured-fraction", logger);
break;
case OPT_MOSAICITY:
forced_mosaicity_arg = parse_double_arg(optarg, "--mosaicity", logger);
break;
case OPT_PREDICTION_MOSAICITY:
forced_prediction_mosaicity_arg = parse_double_arg(optarg, "--prediction-mosaicity", logger);
break;
case OPT_INTEGRATION_RADIUS:
integration_radius_arg = optarg;
break;
case OPT_ADAPTIVE_INTEGRATION_RADIUS:
adaptive_integration_radius_arg = !optarg || std::string(optarg) != "off";
break;
case OPT_INTEGRATION_STENCIL:
integration_stencil_arg = optarg;
break;
case OPT_BACKGROUND_TRIM:
background_trim_arg = parse_double_arg(optarg, "--background-trim", logger);
break;
case OPT_OVERLAP:
if (strcmp(optarg, "off") == 0) overlap_arg = OverlapMode::Off;
else if (strcmp(optarg, "reject") == 0) overlap_arg = OverlapMode::Reject;
else if (strcmp(optarg, "exclude") == 0) overlap_arg = OverlapMode::Exclude;
else { logger.Error("--overlap expects off|reject|exclude"); return 1; }
break;
case OPT_OVERLAP_MINPK:
overlap_minpk_arg = parse_double_arg(optarg, "--overlap-minpk", logger);
break;
case OPT_MAX_HKL:
max_hkl_arg = parse_number_arg<int64_t>(optarg, "--max-hkl", logger, 1, 511);
break;
case OPT_INTEGRATION_HIGH_RES:
integration_d_min_arg = parse_double_arg(optarg, "--integration-high-resolution", logger);
break;
case OPT_BACKGROUND_RADIAL:
background_radial_given = true;
if (optarg == nullptr || strcmp(optarg, "on") == 0)
background_radial_arg = true;
else if (strcmp(optarg, "off") == 0)
background_radial_arg = false;
else if (strcmp(optarg, "auto") == 0)
background_radial_arg = std::nullopt;
else {
logger.Error("Invalid --background-radial value: {} (expected on|off|auto)", optarg);
return 1;
}
break;
case OPT_BACKGROUND_CLIP:
background_clip_arg = parse_double_arg(optarg, "--background-clip", logger);
break;
case OPT_INTEGRATOR:
if (strcmp(optarg, "boxsum") == 0) integrator_mode = IntegratorMode::BoxSum;
else if (strcmp(optarg, "gaussian") == 0) integrator_mode = IntegratorMode::ProfileGaussian;
else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical;
else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; }
break;
case OPT_SIMPLE_STILLS:
simple_stills_flag = true;
break;
case OPT_REJECT_OUTLIERS:
outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger);
break;
case OPT_MIN_IMAGE_CC:
min_image_cc = parse_double_arg(optarg, "--min-image-cc", logger);
break;
case OPT_SEARCH_MIN_ZETA:
search_min_zeta_arg = parse_double_arg(optarg, "--search-min-zeta", logger);
break;
case OPT_SCALING_HIGH_RESOLUTION:
d_min_scale_merge = parse_number_arg<float>(optarg, "--scaling-high-resolution", logger,
0.1f, 1000.0f);
break;
case OPT_SCALING_LOW_RESOLUTION:
d_max_scale_merge = parse_number_arg<float>(optarg, "--scaling-low-resolution", logger,
0.0f, 100000.0f);
break;
case OPT_RESOLUTION_CUTOFF:
if (strcmp(optarg, "cc-logistic") == 0)
resolution_cutoff_method = ResolutionCutoffMethod::CCHalfLogistic;
else if (strcmp(optarg, "off") == 0)
resolution_cutoff_method = ResolutionCutoffMethod::Off;
else {
logger.Error("Invalid --resolution-cutoff value: {} (expected cc-logistic|off)", optarg);
print_usage();
exit(EXIT_FAILURE);
}
break;
case OPT_RESOLUTION_CC_TARGET:
resolution_cc_target = parse_double_arg(optarg, "--resolution-cc-target", logger);
break;
case OPT_RESOLUTION_SHELLS:
report_shell_count = atoi(optarg);
if (report_shell_count.value() < 1) {
logger.Error("Invalid --resolution-shells value: {} (must be >= 1)", report_shell_count.value());
exit(EXIT_FAILURE);
}
break;
case OPT_FORCE_STILL:
force_still = true;
break;
case 'q':
q_spacing = atof(optarg);
break;
case OPT_AZIM_MIN_Q:
min_q = atof(optarg);
break;
case OPT_AZIM_MAX_Q:
max_q = atof(optarg);
break;
case OPT_AZIM_PHI_BINS:
azimuthal_bins = atoi(optarg);
break;
case OPT_POLARIZATION_CORRECTION: {
bool value;
if (!parse_on_off(optarg, value)) {
logger.Error("Invalid polarization correction value (expected on|off): {}", optarg);
exit(EXIT_FAILURE);
}
polarization_correction = value;
break;
}
case OPT_SOLID_ANGLE_CORRECTION: {
bool value;
if (!parse_on_off(optarg, value)) {
logger.Error("Invalid solid angle correction value (expected on|off): {}", optarg);
exit(EXIT_FAILURE);
}
solid_angle_correction = value;
break;
}
case OPT_BEAM_X: beam_x = parse_float_arg(optarg, "--beam-x", logger); break;
case OPT_BEAM_Y: beam_y = parse_float_arg(optarg, "--beam-y", logger); break;
case OPT_DETECTOR_DISTANCE: detector_distance_mm = parse_float_arg(optarg, "--detector-distance", logger); break;
case OPT_WAVELENGTH: {
// Guard > 0: wavelength is used as a divisor (WVL_1A_IN_KEV / wavelength) below, and a 0
// would produce a non-finite incident energy that throws unguarded and aborts the process.
float w = parse_float_arg(optarg, "--wavelength", logger);
if (!(w > 0.0f)) {
logger.Error("Invalid wavelength (must be > 0 A): {}", optarg);
exit(EXIT_FAILURE);
}
wavelength_A = w;
break;
}
case OPT_ROT1: rot1_rad = parse_float_arg(optarg, "--rot1", logger); break;
case OPT_ROT2: rot2_rad = parse_float_arg(optarg, "--rot2", logger); break;
case OPT_ROT3: rot3_rad = parse_float_arg(optarg, "--rot3", logger); break;
case OPT_DETECTOR_MIRROR_Y: detector_mirror_y = true; break;
case OPT_DETECTOR_QUARTER_TURNS:
detector_quarter_turns = atoi(optarg);
if (detector_quarter_turns < 0 || detector_quarter_turns > 3) {
logger.Error("--detector-quarter-turns must be 0, 1, 2 or 3");
return 1;
}
break;
case OPT_FFT_MIN_UNIT_CELL:
fft_min_unit_cell_A = parse_float_arg(optarg, "--fft-min-unit-cell", logger); break;
case OPT_POLARIZATION: polarization_factor = parse_float_arg(optarg, "--polarization", logger); break;
case OPT_SCALING_ITERATIONS:
scaling_iter = atoi(optarg);
if (scaling_iter <= 0) {
logger.Error("Invalid scaling iteration count: {}", scaling_iter);
exit(EXIT_FAILURE);
}
break;
case OPT_BANDWIDTH:
bandwidth_fwhm = atof(optarg);
if (!(bandwidth_fwhm.value() >= 0.0f)) {
logger.Error("Invalid bandwidth: {}", optarg);
exit(EXIT_FAILURE);
}
break;
default:
print_usage();
exit(EXIT_FAILURE);
}
}
if (optind != argc - 1) {
logger.Error("Input file not specified");
print_usage();
exit(EXIT_FAILURE);
}
input_file = argv[optind];
logger.Verbose(verbose);
// -N defaults to 0 = "use all hardware threads"; resolve it to a concrete count here so every mode
// behaves the same. The scale/merge engines expand 0 on their own, but the per-image processing
// loop (Rugnux) spawns exactly nthreads workers, so passing 0 there would spawn none and process
// nothing - hence resolving it centrally rather than relying on each consumer.
const bool nthreads_auto = nthreads <= 0;
if (nthreads <= 0) {
unsigned int hw = std::thread::hardware_concurrency();
nthreads = hw > 0 ? static_cast<int>(hw) : 1;
}
if (space_group != nullptr)
logger.Info("Using space group {} (number {})", space_group->hm, space_group->number);
// Off the critical path, so cuFFT's one-time initialisation overlaps the file open, the mask read
// and the beam-stop pre-scan instead of landing inside the first pass.
auto cufft_warmup = std::async(std::launch::async, WarmUpCuFFT);
// 1. Read Input File. Either an HDF5 master, or one PILATUS miniCBF frame - naming any frame of a
// sweep reads the whole sweep, which is how these are archived and how a user thinks of them.
JFJochHDF5Reader hdf5_reader;
JFJochCBFReader cbf_reader;
JFJochReader *reader_ptr = nullptr;
const bool input_is_cbf = JFJochCBFReader::CanRead(input_file);
try {
if (input_is_cbf) {
cbf_reader.ReadFiles(input_file);
reader_ptr = &cbf_reader;
logger.Info("Read {} CBF images of the sweep containing {}",
cbf_reader.GetNumberOfImages(), input_file);
} else {
hdf5_reader.ReadFile(input_file);
reader_ptr = &hdf5_reader;
}
} catch (const std::exception &e) {
logger.Error("Error reading input file: {}", e.what());
exit(EXIT_FAILURE);
}
JFJochReader &reader = *reader_ptr;
const auto dataset = reader.GetDataset();
if (!dataset) {
logger.Error("No experiment dataset found in the input file");
exit(EXIT_FAILURE);
}
if (rotation_indexing_image_count <= 0) {
logger.Error("Invalid number of rotation indexing images: {}", rotation_indexing_image_count);
exit(EXIT_FAILURE);
}
logger.Info("Loaded dataset from {}", input_file);
std::vector<MergedReflection> reference_data;
bool reference_has_free_flags = false;
if (!ref_mtz.empty()) {
try {
const auto reference = LoadReferenceMtz(
ref_mtz, ref_column.empty() ? std::nullopt : std::optional<std::string>(ref_column));
reference_data = reference.reflections;
reference_has_free_flags = reference.has_free_flags;
logger.Info("Loaded {} reference reflections from {} (column {}{}{})",
reference_data.size(), ref_mtz, reference.used_column,
reference.squared ? ", squared to intensity" : "",
reference.default_column ? ", auto-selected" : ", user-specified");
if (reference.has_free_flags)
logger.Info("Reference carries R-free flags (column {}): {} of {} free; the merged "
"reflections will inherit this test set",
reference.free_column, reference.n_free, reference_data.size());
if (reference.d_max > 0.0)
logger.Info("Reference resolution range {:.2f} - {:.2f} A", reference.d_max, reference.d_min);
if (reference.cell.has_value())
logger.Info("Reference unit cell: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
reference.cell->a, reference.cell->b, reference.cell->c,
reference.cell->alpha, reference.cell->beta, reference.cell->gamma);
if (!reference.space_group_name.empty())
logger.Info("Reference space group: {} (number {})",
reference.space_group_name, reference.space_group_number.value_or(0));
// Check the reference against the cell that will actually drive the merge. --scale merges
// in the cell stored in the input file (as the former jfjoch_scale did); the -C override
// only takes effect on the full-analysis path, which otherwise determines its cell later by
// indexing (unknown here, so nothing can be checked yet).
const std::optional<UnitCell> data_cell =
mode == RugnuxMode::Scale ? dataset->experiment.GetUnitCell() : fixed_reference_unit_cell;
const auto warning = ReferenceConsistencyWarning(
reference, data_cell,
space_group != nullptr ? std::optional<int>(space_group->number) : std::nullopt);
if (!warning.empty())
logger.Warning("{}", warning);
// A reference MTZ fixes the space group and unit cell, unless -S / -C override them.
// (-S with the wrong enantiomorph, or -C with a different cell, is allowed - the explicit
// flag always wins.) The cell is a soft reference: indexing may drift within tolerance.
if (space_group == nullptr && reference.space_group_number.has_value()) {
space_group = gemmi::find_spacegroup_by_number(static_cast<int>(*reference.space_group_number));
logger.Info("Fixing space group from reference MTZ: {} ({})",
reference.space_group_name, *reference.space_group_number);
}
if (!fixed_reference_unit_cell.has_value() && reference.cell.has_value()) {
fixed_reference_unit_cell = reference.cell;
logger.Info("Fixing reference unit cell from reference MTZ (indexing may drift within tolerance)");
}
} catch (const std::exception &e) {
logger.Error("Error reading reference MTZ {}: {}", ref_mtz, e.what());
exit(EXIT_FAILURE);
}
}
// --mode scale: re-scale and merge the already-integrated reflections stored in the input file,
// without re-running spot finding or integration (folded in from the former rugnux_scale tool).
if (mode == RugnuxMode::Scale) {
// Re-scaling reads reflections a previous run integrated, which only a _process.h5 holds.
if (input_is_cbf) {
logger.Error("--mode scale needs the integrated reflections in a _process.h5; "
"a CBF sweep holds none. Run --mode mx on it first.");
exit(EXIT_FAILURE);
}
const auto total_images = static_cast<int>(reader.GetNumberOfImages());
const int last_image = (end_image < 0 || end_image >= total_images) ? total_images - 1 : end_image;
auto reflections = hdf5_reader.ReadReflections(start_image, last_image);
DiffractionExperiment experiment(dataset->experiment);
configure_offline_output(experiment, output_prefix);
// The reflections in the file are already indexed, so the cell and space group they were
// integrated in are the file's to supply here - but an explicit -S / -C still wins.
if (space_group != nullptr)
experiment.SetSpaceGroup(*space_group);
if (fixed_reference_unit_cell.has_value())
experiment.SetUnitCell(fixed_reference_unit_cell);
// A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling.
IndexingSettings indexing_settings;
indexing_settings.RotationIndexing(experiment.GetGoniometer().has_value() && !force_still);
// --detect-ice-rings, applied here as well as on the full path below: this block returns
// before that one runs, so without it the flag is silently ignored by --scale. Same
// precedence as there - command line, then the file, then the geometry's default.
if (detect_ice_rings.has_value())
experiment.DetectIceRings(detect_ice_rings.value());
else if (!dataset->file_detect_ice_rings.has_value())
experiment.DetectIceRings(indexing_settings.GetRotationIndexing());
experiment.ImportIndexingSettings(indexing_settings);
// Start from the same defaults the full pipeline uses, so --scale reproduces the merge that wrote
// the _process.h5 rather than a weaker model of its own. Every CLI override below then applies on
// top, exactly as in the full-analysis path.
const bool rot = experiment.GetGoniometer().has_value() && !force_still;
ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rot);
if (d_min_scale_merge)
scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
if (d_max_scale_merge)
scaling_settings.LowResolutionLimit_A(*d_max_scale_merge > 0.0f
? std::optional<double>(*d_max_scale_merge) : std::nullopt);
if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
scaling_settings.MergeFriedel(!anomalous_mode);
scaling_settings.MinPartiality(min_partiality);
scaling_settings.MinCapturedFraction(
min_captured_fraction_arg.value_or(scaling_settings.GetMinCapturedFraction()));
scaling_settings.CaptureUncertaintyCoeff(
capture_uncertainty_arg.value_or(scaling_settings.GetCaptureUncertaintyCoeff()));
scaling_settings.ForcedMosaicity(forced_mosaicity_arg);
scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is percent; the setting is a fraction
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
if (ice_min_score_arg)
scaling_settings.IceMinScore(static_cast<float>(*ice_min_score_arg));
if (ice_min_spot_ratio_arg)
scaling_settings.IceMinSpotRatio(static_cast<float>(*ice_min_spot_ratio_arg));
scaling_settings.OutlierRejectNsigma(
outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma()));
scaling_settings.ScaleFulls(scale_fulls_arg.value_or(scaling_settings.GetScaleFulls()));
scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees()));
scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off
if (no_scaling_corrections)
scaling_settings.CorrectionSurfaces(false);
experiment.ImportScalingSettings(scaling_settings);
if (!experiment.GetUnitCell()) {
logger.Error("Experiment unit cell not found, cannot update reflection resolution");
exit(EXIT_FAILURE);
}
// The stored reflections have to be in the setting of the cell they are about to be merged in,
// or the merge's per-frame cell gate rejects every observation and the run dies far downstream
// complaining about resolution shells. The reader has already applied /entry/MX/reindexMatrix,
// so a file that carries one arrives consistent; one written before that dataset existed, by a
// run whose space-group search re-seated the lattice after the images had gone to file, does
// not - and the change of basis cannot be recovered from it, because the stored lattices carry
// each frame's own goniometer rotation. Say so here, where both cells can be named, instead of
// merging indices in a frame they do not belong to.
{
const auto &merge_cell = experiment.GetUnitCell().value();
const float dist_tol = experiment.GetIndexingSettings().GetUnitCellDistTolerance();
const float ang_tol = experiment.GetIndexingSettings().GetUnitCellAngleTolerance_deg();
size_t n_indexed = 0, n_agree = 0;
std::optional<UnitCell> file_cell;
for (const auto &o : reflections)
if (!o.reflections.empty()) {
++n_indexed;
const auto cell = o.latt.GetUnitCell();
if (cell.is_close(merge_cell, dist_tol, ang_tol))
++n_agree;
else if (!file_cell)
file_cell = cell;
}
if (n_indexed > 0 && n_agree == 0) {
logger.Error("The per-image lattices stored in this file are not in the setting of the "
"cell it is being merged in, so no observation can be used. Merging in "
"a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} gamma={:.2f}, but all "
"{} indexed images are in a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} "
"beta={:.2f} gamma={:.2f}{}. Merge in the setting the file is in with "
"-S 1 -C {:.4f},{:.4f},{:.4f},{:.4f},{:.4f},{:.4f}, or reprocess the "
"images with --mode mx.",
merge_cell.a, merge_cell.b, merge_cell.c, merge_cell.alpha,
merge_cell.beta, merge_cell.gamma, n_indexed,
file_cell->a, file_cell->b, file_cell->c, file_cell->alpha,
file_cell->beta, file_cell->gamma,
dataset->reindex_matrix
? "" : " (the file carries no /entry/MX/reindexMatrix, so it was "
"written before that was recorded)",
file_cell->a, file_cell->b, file_cell->c, file_cell->alpha,
file_cell->beta, file_cell->gamma);
exit(EXIT_FAILURE);
}
}
auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections);
logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images);
experiment.ImagesPerTrigger(refl_stats.n_images);
// Ice-ring handling, as the full pipeline does it (Rugnux.cpp): flag reflections on a
// hexagonal-ice powder ring so scaling skips them while the merge keeps them. The flag is not
// stored per reflection, so it has to be recomputed here from the resolution just assigned -
// otherwise --scale re-scales a dataset the writing run had scaled without those reflections,
// and the per-image scales come out of a different fit than the ones in the file.
const float ice_width = SpotFindingSettings().ice_ring_width_Q_recipA;
// ...and gated the same way, on the per-image ice score the writing run stored in the file, so
// --scale reaches the same verdict on the same data as the pipeline that produced it.
double ice_sum = 0.0;
size_t ice_n = 0;
for (const float s : dataset->ice_ring_score)
if (std::isfinite(s)) {
ice_sum += s;
++ice_n;
}
const float ice_min_score = experiment.GetScalingSettings().GetIceMinScore();
// ...and the spot channel, pooled over the run exactly as the full pipeline pools it.
double ring_sum = 0.0, ctrl_sum = 0.0;
for (const float v : dataset->spot_count_ice_rings)
if (std::isfinite(v)) ring_sum += v;
for (const float v : dataset->spot_count_ice_control)
if (std::isfinite(v)) ctrl_sum += v;
// Empty control + spots on the rings = the strongest ice evidence, not its absence.
const double ice_spot_ratio = ctrl_sum > 0.0 ? ring_sum / ctrl_sum : (ring_sum > 0.0 ? 1.0e3 : 0.0);
const float ice_min_spot_ratio = experiment.GetScalingSettings().GetIceMinSpotRatio();
const bool ice_present = (ice_n == 0 || ice_sum / static_cast<double>(ice_n) >= ice_min_score)
|| (ice_min_spot_ratio > 0.0f && ice_spot_ratio >= ice_min_spot_ratio);
if (experiment.IsDetectIceRings() && !ice_present) {
logger.Info("Ice-ring handling: measured ice score {:.2f} and spot ratio {:.2f} below the "
"gates ({:.2f} / {:.2f}), no ice detected - ice-ring handling skipped entirely",
ice_sum / static_cast<double>(ice_n), ice_spot_ratio, ice_min_score,
ice_min_spot_ratio);
} else if (experiment.IsDetectIceRings()) {
size_t total = 0, flagged = 0;
for (auto &outcome : reflections) {
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: flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
"excluded from scaling, kept for merging", flagged, total, ice_width);
}
const auto scale_start = std::chrono::steady_clock::now();
std::vector<MergedReflection> merged_reflections;
MergeStatistics merged_statistics;
double error_model_isa = 0.0;
double error_model_isa_asymptotic = 0.0;
std::optional<double> resolution_fit_A; // CC1/2 crossing, without the one-shell extension
double error_model_a = 0.0;
double error_model_b = 0.0;
// Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine ->
// merge, including the default-on decay + absorption correction surfaces. It does not support
// external-reference scaling or wedge refinement.
// Everything else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly.
const bool is_rotation = experiment.IsRotationIndexing();
if (is_rotation) {
if (!reference_data.empty()
|| experiment.GetRefineRotationWedgeInScaling()
|| experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Rotation scaling/merging (RotationScaleMerge) does not support reference "
"scaling or wedge refinement");
RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(),
scaling_iter, nthreads, logger);
rsm.Ingest();
// No pass to compare against here - --mode scale merges the stored reflections once - so
// decline the extra merge the pre-correction CC1/2 would cost.
auto r = rsm.Run(false, /*full_stats=*/true, /*measure_cc_before_corrections=*/false);
merged_reflections = std::move(r.merged);
merged_statistics = std::move(r.statistics);
error_model_isa = r.isa;
error_model_isa_asymptotic = r.isa_asymptotic;
error_model_a = r.error_model_a;
error_model_b = r.error_model_b;
resolution_fit_A = r.resolution_fit_A;
} else {
// Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports
// CCref and provides the R-free test set, but is NOT 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. The per-image scale G is the exact one-pass solution, so one pass
// (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills).
ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads);
// Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
// orientation tilt against the merge and recompute each reflection's partiality + scale correction
// (no re-integration), then merge with the improved corrections.
if (experiment.GetScalingSettings().GetStillsPartialityRefine()) {
StillsPartialityRefine refiner(experiment);
const double mean_tilt = refiner.Run(reflections, nthreads);
logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
}
MergeOnTheFly merge_engine(experiment);
merge_engine.ReferenceCell(experiment.GetUnitCell());
// --min-image-cc has to hold for the merge itself, not only for the reported statistics.
merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0);
// Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full
// pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity
// model and produce much worse stills intensities (no (b*I)^2 systematic term at all).
merge_engine.RefineErrorModel(reflections);
if (merge_engine.ErrorModelActive())
{
// XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is
// already the whole-range 1/b, there being no asymptote tier on this path.
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}",
em.a, em.b, em.isa, merge_engine.ErrorModelChi2());
}
for (size_t i = 0; i < reflections.size(); ++i)
merge_engine.AddImage(reflections[i], static_cast<int64_t>(i));
merged_reflections = merge_engine.ExportReflections();
// Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual
// --scaling-high-resolution wins, otherwise trim the written reflections + reported shells
// to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.)
const auto &cut_ss = experiment.GetScalingSettings();
// The offline --scale path re-scales a stored _process.h5 and is never a P1 search merge.
const std::optional<double> effective_d_min = ApplyResolutionCutoff(
merged_reflections, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(),
cut_ss.GetResolutionCCTarget(), /*for_search=*/false, logger, &resolution_fit_A);
merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data,
effective_d_min);
{ // stills: no asymptote tier, and its ISa already IS the whole-range 1/b
const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
error_model_isa = em.isa;
error_model_a = em.a;
error_model_b = em.b;
}
}
logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)",
std::chrono::duration<double>(std::chrono::steady_clock::now() - scale_start).count(),
merged_reflections.size());
// Inherit the campaign's shared R-free test set from the reference MTZ (overriding the
// per-hkl hash the merge assigned), so every dataset flags the same free reflections.
if (reference_has_free_flags && !reference_data.empty() && !merged_reflections.empty()) {
const size_t matched = ApplyReferenceFreeFlags(merged_reflections, experiment.GetSpaceGroupOrP1(),
reference_data);
logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)",
matched, merged_reflections.size());
}
// Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B), as the full pipeline reports
// it (Rugnux.cpp). Diagnostic only - it is not fed back into scaling; it lands in the printed
// statistics, the mmCIF and the log. Without it --mode scale reports WILSON_B= nan.
{
const GlobalWilsonB wilson = CalcGlobalWilsonB(merged_reflections);
merged_statistics.wilson_b = wilson.b;
merged_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);
}
std::cout << merged_statistics;
// Space-group determination lives in the full rugnux pipeline; --scale only consumes a space
// group (from the file or -S) and merges in it.
const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
if (!fixed_space_group)
logger.Warning("No space group in the input file or on the command line - merged in P1. "
"Re-run rugnux (which determines and stores the space group) or pass "
"-S to scale and merge in the correct symmetry.");
const auto &twin_sg_opt = experiment.GetGemmiSpaceGroup();
const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr;
const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg);
std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl;
// Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - report-only,
// nothing is corrected. The stored reflections have just been re-scaled above, so they carry
// this merge's own per-image scale.
if (experiment.GetUnitCell()) {
AnisotropyRunInfo aniso_run;
if (merged_statistics.sweep_quality.measured && merged_statistics.sweep_quality.sweep_deg > 0.0f)
aniso_run.observed_rotation_deg = merged_statistics.sweep_quality.sweep_deg;
aniso_run.dose_term_in_scale_model = experiment.GetScalingSettings().GetCorrectionSurfaces();
aniso_run.radiation_damage_relative_b = merged_statistics.radiation_damage_delta_b;
merged_statistics.anisotropy = AnalyzeAnisotropy(merged_reflections,
ScaledObservations(reflections, is_rotation, twin_sg),
*experiment.GetUnitCell(), twin_sg, aniso_run);
std::cout << AnisotropyToText(merged_statistics.anisotropy) << std::endl;
}
// Before the reflection files, as in the full pipeline: the model settles the enantiomorph and,
// with no reference MTZ, the alternative indexing, and the files are written in what it settled.
// A model that was asked for and could not be used has to reach the report, or the run ends
// looking exactly like one that was never given --model.
std::string model_validation_failure;
std::optional<ModelValidationResult> model_validation;
if (!output_prefix.empty() && !model_pdb.empty()) {
const auto data_sg = experiment.GetGemmiSpaceGroup();
// With a reference MTZ the merohedral indexing was already resolved (stills per-image
// scaling); only probe indexing by R-free when model-only, with no reference.
const auto validation =
ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb,
output_prefix, logger,
data_sg ? &*data_sg : nullptr,
/*probe_indexing_ambiguity=*/reference_data.empty());
model_validation = validation;
if (!validation.failure_reason.empty())
model_validation_failure = validation.failure_reason;
if (data_sg.has_value()) {
const gemmi::SpaceGroup *adopted =
AdoptModelFrame(validation, merged_reflections, *data_sg,
experiment.GetScalingSettings().GetMergeFriedel(), logger);
if (adopted->number != data_sg->number)
experiment.SetSpaceGroup(*adopted);
// The unmerged export below writes the integrated observations, so it follows the same
// relabeling as the merged set - only the indexing moves a reflection; the enantiomorph
// is a label - or the two files would describe different indexings of one dataset.
const gemmi::Op &op = validation.indexing_op;
if (!(op == gemmi::Op::identity()))
for (auto &io : reflections)
for (auto &r : io.reflections) {
const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}});
r.h = h[0]; r.k = h[1]; r.l = h[2];
}
}
}
// Unmerged observations, from the integrated observations rather than the merged ones: the
// partiality and the per-image scale are left for the reading program, which fits a scale
// model of its own. Without an output prefix there is nowhere to put it, and the name would
// come out as a bare "_unmerged.mtz".
if (export_unmerged && !output_prefix.empty()) {
const std::string path = output_prefix + "_unmerged.mtz";
WriteUnmergedMtzReflections(reflections, *experiment.GetUnitCell(), experiment, true, path);
logger.Info("Unmerged observations written to {}", path);
}
if (export_unmerged_partials && !output_prefix.empty()) {
const std::string path = output_prefix + "_unmerged_partials.mtz";
WriteUnmergedMtzReflections(reflections, *experiment.GetUnitCell(), experiment, false, path);
logger.Info("Unmerged partials written to {}", path);
}
if (!output_prefix.empty()) {
const ErrorModelReport em_report{
error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?",
error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", error_model_isa_asymptotic)
: std::string(),
error_model_a > 0 ? fmt::format("{:.3f}", error_model_a) : std::string(),
error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string()};
WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics,
em_report, twinning, output_prefix, static_cast<size_t>(nthreads));
}
// --mode scale re-merges stored reflections, so it determines a merging result and gets the
// same report. It has no ProcessResult of its own; fill the parts that exist here.
ProcessResult scale_result;
scale_result.images_processed = reader.GetNumberOfImages();
scale_result.consensus_cell = experiment.GetUnitCell();
scale_result.space_group = experiment.GetGemmiSpaceGroup();
scale_result.used_beam_x_pxl = experiment.GetBeamX_pxl();
scale_result.used_beam_y_pxl = experiment.GetBeamY_pxl();
scale_result.used_distance_mm = experiment.GetDetectorDistance_mm();
scale_result.has_merge_statistics = true;
scale_result.merge_statistics = merged_statistics;
{ std::ostringstream s; s << merged_statistics; scale_result.merge_statistics_text = s.str(); }
scale_result.error_model_isa = error_model_isa;
scale_result.error_model_isa_asymptotic = error_model_isa_asymptotic;
scale_result.error_model_a = error_model_a;
scale_result.error_model_b = error_model_b;
scale_result.resolution_fit_A = resolution_fit_A;
scale_result.has_reference = !reference_data.empty();
scale_result.twinning = twinning;
scale_result.model_validation = model_validation;
if (!model_validation_failure.empty())
scale_result.warnings.push_back("Model validation did not run: " + model_validation_failure);
provenance.wall_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - invocation_start).count();
std::cout << fmt::format("Total wall time: {:.2f} s", provenance.wall_time_s) << std::endl;
WriteResultReport(output_prefix, input_file, experiment, scale_result, logger, provenance);
return 0;
}
uint64_t total_images_in_file = reader.GetNumberOfImages();
if (end_image < 0 || end_image > total_images_in_file)
end_image = total_images_in_file;
if (image_stride < 0) {
logger.Error("Image stride cannot be negative");
exit(EXIT_FAILURE);
}
if (image_stride == 0) {
logger.Error("Image stride cannot be zero");
exit(EXIT_FAILURE);
}
// Round UP: the worker loop walks ordinals until start + ordinal * stride reaches end_image, so a
// stride that does not divide the range evenly still processes a final partial step. Rounding down
// here would size the per-image arrays one short of the ordinals that index them.
int images_to_process = (end_image - start_image + image_stride - 1) / image_stride;
if (images_to_process <= 0) {
logger.Warning("No images to process (Start: {}, End: {} Stride: {}, Total: {})", start_image, end_image,
image_stride, total_images_in_file);
return 0;
}
// 2. Setup Experiment & Components
DiffractionExperiment experiment(dataset->experiment);
// Geometry overrides (default: keep the value stored in the input file). Applied before the
// azimuthal-integration settings are derived, which depend on the geometry.
if (beam_x) experiment.BeamX_pxl(beam_x.value());
if (beam_y) experiment.BeamY_pxl(beam_y.value());
// A centre the user typed in is the answer, not a starting point.
if (beam_x || beam_y) estimate_beam_center = false;
if (detector_distance_mm) experiment.DetectorDistance_mm(detector_distance_mm.value());
if (wavelength_A) experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A.value());
if (rot1_rad) experiment.PoniRot1_rad(rot1_rad.value());
if (rot2_rad) experiment.PoniRot2_rad(rot2_rad.value());
if (rot3_rad) experiment.PoniRot3_rad(rot3_rad.value());
// The discrete part of the detector orientation, which the file does not state for a detector this
// system did not assemble. Never derived from rot3 - it says how the image is stored, and a fitted
// angle must not change that.
if (detector_mirror_y || detector_quarter_turns != 0)
experiment.Detector().ImageOrientation(
DetectorOrientation(detector_mirror_y, detector_quarter_turns));
// --polarization is applied after configure_offline_output below, which sets the rugnux default.
// Calibrating from the run-summed profile needs the profile to be binned in azimuth - see
// CALIBRATION_AZIM_BINS_DEFAULT. An explicit --azim-phi-bins always wins.
if (mode == RugnuxMode::Calibration && calibration_method == CalibrationMethod::Rings
&& !azimuthal_bins.has_value())
azimuthal_bins = CALIBRATION_AZIM_BINS_DEFAULT;
// Azimuthal integration (default q-spacing 0.01 1/A, from AzimuthalIntegrationSettings): the profile
// resolves the narrow ice rings for the ice-ring score. Shared by every mode that reads images.
// -q / --azim-* / correction flags override; defaults come from the input file.
{
AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings();
if (min_q || max_q)
azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()),
max_q ? max_q : azint_settings.GetRequestedHighQ_recipA());
if (q_spacing)
azint_settings.QSpacing_recipA(q_spacing.value());
if (azimuthal_bins)
azint_settings.AzimuthalBinCount(azimuthal_bins.value());
if (polarization_correction)
azint_settings.PolarizationCorrection(polarization_correction.value());
if (solid_angle_correction)
azint_settings.SolidAngleCorrection(solid_angle_correction.value());
experiment.ImportAzimuthalIntegrationSettings(azint_settings);
logger.Info("Azimuthal integration: Q [{:.4f}, {:.4f}] 1/A, spacing {:.4f}, {} Q x {} azimuthal bins",
azint_settings.GetLowQ_recipA(), azint_settings.GetHighQ_recipA(),
azint_settings.GetQSpacing_recipA(), azint_settings.GetQBinCount(),
azint_settings.GetAzimuthalBinCount());
}
// --mode azint: azimuthal integration only (no spot finding / indexing / scaling). Rugnux reads
// the geometry and azimuthal-integration settings configured above off the experiment.
if (mode == RugnuxMode::Azint) {
// As for --mode calibration below: the polarization factor is a rugnux analysis default and it
// scales the profile this mode exists to produce, so apply the defaults and then let an explicit
// --polarization override them. Without this the flag was parsed and then dropped on this path.
ApplyRugnuxExperimentDefaults(experiment);
if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value());
ProcessConfig config;
config.mode = ProcessMode::AzimuthalIntegration;
config.start_image = start_image;
config.end_image = end_image;
config.stride = image_stride;
config.nthreads = nthreads;
config.nthreads_auto = nthreads_auto;
config.output_prefix = output_prefix;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
Rugnux process(reader, experiment, *dataset->pixel_mask, config);
g_active_process = &process;
std::signal(SIGINT, handle_sigint);
ProcessResult result;
try {
result = process.Run();
} catch (const std::exception &e) {
logger.Error("Processing failed: {}", e.what());
exit(EXIT_FAILURE);
}
g_active_process = nullptr;
std::cout << fmt::format("Processing time: {:.2f} s", result.total_time_s) << std::endl;
std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl;
std::cout << fmt::format("Total throughput: {:.2f} MB/s", result.throughput_MBs) << std::endl;
if (result.cancelled)
logger.Warning("Processing was cancelled after {} images", result.images_processed);
return 0;
}
// --mode calibration: fit the detector geometry to the calibrant's powder rings, measured over
// every selected image, and write it as a PONI file.
if (mode == RugnuxMode::Calibration) {
// The same analysis defaults every other mode gets (configure_offline_output applies them for
// the MX path further down). The polarization factor is one of them and it enters the summed
// azimuthal profile the rings are measured in, so without this the CLI and the viewer's
// "Analyze dataset", which always applies them, would fit slightly different ring positions.
ApplyRugnuxExperimentDefaults(experiment);
if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value());
ProcessConfig config;
config.mode = ProcessMode::Calibration;
config.calibration_method = calibration_method;
config.calibration_refine_tilt = calibration_refine_tilt;
// -C wins over --calibrant: a cell given on the command line IS the standard, and the built-in
// table is only the convenience for the five that have names. -S, where it is given, supplies
// the absences - without it the cell is taken as primitive, which for a centred standard would
// open the ring list with a reflection that is not there and scale the whole calibration.
std::string calibrant_label = calibrant;
if (fixed_reference_unit_cell.has_value()) {
config.calibrant_ring_q = space_group != nullptr
? CalculateXtalRings(*fixed_reference_unit_cell, *space_group)
: CalculateXtalRings(*fixed_reference_unit_cell);
calibrant_label = fmt::format("cell {:.4f},{:.4f},{:.4f},{:.3f},{:.3f},{:.3f} ({})",
fixed_reference_unit_cell->a, fixed_reference_unit_cell->b,
fixed_reference_unit_cell->c, fixed_reference_unit_cell->alpha,
fixed_reference_unit_cell->beta, fixed_reference_unit_cell->gamma,
space_group != nullptr ? space_group->hm : "assumed primitive");
if (config.calibrant_ring_q.empty()) {
logger.Error("The unit cell given with -C produces no rings");
return 1;
}
} else {
config.calibrant_ring_q = CalibrantRings(calibrant);
}
config.start_image = start_image;
config.end_image = end_image;
config.stride = image_stride;
config.nthreads = nthreads;
config.nthreads_auto = nthreads_auto;
config.output_prefix = output_prefix;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
config.write_process_h5 = false; // the .poni below is the output of this mode
// Ice-ring handling has no place in a calibration run and does active harm. Flagged spots are
// sorted LAST by the spot budget (FilterSpotsByCount), so they are the first thrown away - which
// for --calibrant ice discards exactly the spots being calibrated on, and for every other
// standard discards whichever of its rings happen to fall in the fixed hexagonal bands.
experiment.DetectIceRings(false);
// Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions
// and nothing else, and the calibrant is a powder with no lattice to index.
config.spot_finding.enable = true;
config.spot_finding.indexing = false;
config.spot_finding.signal_to_noise_threshold = sigma_spot_finding;
config.spot_finding.photon_count_threshold = photon_count_threshold_spot_finding;
config.spot_finding.min_pix_per_spot = min_pix_per_spot.value_or(2);
config.spot_finding.adaptive_threshold = adaptive_spots.value_or(true);
config.spot_finding.high_resolution_limit = d_min_spot_finding;
config.spot_finding.false_pixels_per_frame = false_pixels_per_frame;
config.spot_finding.low_resolution_limit = d_max_spot_finding;
experiment.MaxSpotCount(max_spot_count_override.value_or(RUGNUX_MAX_SPOT_COUNT));
logger.Info("Powder calibration: calibrant {} ({} rings), method {}, {} azimuthal bins, tilt {}",
calibrant_label, config.calibrant_ring_q.size(),
calibration_method == CalibrationMethod::Rings ? "rings" : "spots",
experiment.GetAzimuthalIntegrationSettings().GetAzimuthalBinCount(),
calibration_refine_tilt ? "refined" : "fixed");
Rugnux process(reader, experiment, *dataset->pixel_mask, config);
g_active_process = &process;
std::signal(SIGINT, handle_sigint);
ProcessResult result;
try {
result = process.Run();
} catch (const std::exception &e) {
logger.Error("Calibration failed: {}", e.what());
return 1;
}
g_active_process = nullptr;
if (!result.calibration.has_value()) {
logger.Error("Calibration produced no geometry");
return 1;
}
const auto &cal = *result.calibration;
const auto &g = cal.geometry;
const DiffractionGeometry header = experiment.GetDiffractionGeometry();
const float pxl_mm = g.GetPixelSize_mm();
constexpr double RAD_TO_DEG = 180.0 / PI;
std::cout << fmt::format("Calibration: {} ring points, rms {:.3f} px "
"(beam centre sigma {:.3f} px)",
cal.ring_points, cal.rms_radial_pxl, cal.beam_sigma_pxl) << std::endl;
// What the rings said the distance was before the fit ran. The fit only converges from a few
// per cent out, so this is what decides whether the run had a chance at all - and a large gap
// from the header is the answer to the question the calibration was run to ask.
if (cal.seed_distance_mm > 0.0f)
std::cout << fmt::format("Ring seed: distance {:.4f} mm from the rings themselves "
"({:+.4f} mm from the header)",
cal.seed_distance_mm,
cal.seed_distance_mm - cal.header_distance_mm) << std::endl;
else
std::cout << "Ring seed: not measurable - the profile shows fewer than two rings"
<< std::endl;
std::cout << fmt::format("PONI x: {:.3f} px = {:.4f} mm ({:+.3f} px from the header)",
g.GetBeamX_pxl(), g.GetBeamX_pxl() * pxl_mm,
g.GetBeamX_pxl() - header.GetBeamX_pxl()) << std::endl;
std::cout << fmt::format("PONI y: {:.3f} px = {:.4f} mm ({:+.3f} px from the header)",
g.GetBeamY_pxl(), g.GetBeamY_pxl() * pxl_mm,
g.GetBeamY_pxl() - header.GetBeamY_pxl()) << std::endl;
std::cout << fmt::format("Rot1: {:+.6f} rad = {:+.4f} deg ({:+.6f} rad from the header)",
g.GetPoniRot1_rad(), g.GetPoniRot1_rad() * RAD_TO_DEG,
g.GetPoniRot1_rad() - header.GetPoniRot1_rad()) << std::endl;
std::cout << fmt::format("Rot2: {:+.6f} rad = {:+.4f} deg ({:+.6f} rad from the header)",
g.GetPoniRot2_rad(), g.GetPoniRot2_rad() * RAD_TO_DEG,
g.GetPoniRot2_rad() - header.GetPoniRot2_rad()) << std::endl;
std::cout << fmt::format("Distance: {:.4f} mm ({:+.4f} mm from the header)",
g.GetDetectorDistance_mm(),
g.GetDetectorDistance_mm() - header.GetDetectorDistance_mm()) << std::endl;
// The PONI is the point of normal incidence, which is what the .poni file wants; the direct beam
// is where the beam actually lands, which is the number every other program reports as the beam
// centre. They part company by distance*tan(rot) as soon as the detector is tilted, so print both.
const auto [beam_x, beam_y] = g.GetDirectBeam_pxl();
std::cout << fmt::format("Direct beam: {:.3f}, {:.3f} px", beam_x, beam_y) << std::endl;
// Where the spots put the beam, independently of everything above. The circle through them
// reads nothing from the header, so it holds in exactly the regime the summed profile does not:
// a profile binned about a badly wrong centre shows each ring smeared across its sectors, and no
// amount of seeding recovers that - the fit then converges, prettily, on the wrong geometry.
// Printed as a fact rather than gated on, because the reader can see at a glance whether two
// methods that share no assumption agree.
if (cal.spots_available)
std::cout << fmt::format("Spot cross-check: the circle through the spots puts the beam at "
"{:.1f}, {:.1f} px - {:.1f} px from the fitted centre",
cal.spots_beam_x_pxl, cal.spots_beam_y_pxl,
cal.spots_disagreement_pxl) << std::endl;
// What the fit says about itself. The scatter line above is about the MEASUREMENTS; this is
// about the PARAMETERS, and the two disagree exactly where the calibration is worth doubting:
// a fit with few rings can sit tightly on the points it has while leaving the tilt free to
// trade tens of pixels of beam centre for itself. The correlations name that trade - rot1
// against the beam in x, rot2 against y - and approach 1 as the two stop being separable.
if (const auto &u = cal.uncertainty; u.valid) {
std::cout << fmt::format("Fit sigma: PONI {:.3f}, {:.3f} px distance {:.4f} mm",
u.sigma_beam_x_pxl, u.sigma_beam_y_pxl, u.sigma_distance_mm)
<< std::endl;
if (u.sigma_rot1_rad > 0.0 || u.sigma_rot2_rad > 0.0)
std::cout << fmt::format(" rot1 {:.4f} deg, rot2 {:.4f} deg "
"correlation with PONI {:+.3f} / {:+.3f}",
u.sigma_rot1_rad * RAD_TO_DEG, u.sigma_rot2_rad * RAD_TO_DEG,
u.corr_beam_x_rot1, u.corr_beam_y_rot2)
<< std::endl;
else
std::cout << " tilt held fixed" << std::endl;
} else {
std::cout << "Fit sigma: not available - the fit is degenerate at its solution"
<< std::endl;
}
// Whether the tilt is a measurement or an artefact of a fit that could not tell it from the beam
// centre. Reported either way: a declined tilt changes what the PONI means, and a kept one is
// not a certified one - the estimator is limited by systematics rather than by this sigma.
if (!calibration_refine_tilt) {
std::cout << "Tilt: held at the header value (--no-refine-tilt)" << std::endl;
} else if (cal.tilt_refined) {
std::cout << fmt::format("Tilt: refined; it stands {:.1f}x its own sigma",
cal.tilt_significance) << std::endl;
} else if (cal.tilt_significance > 0.0f) {
std::cout << fmt::format("Tilt: DECLINED and pinned - it came out at only "
"{:.1f}x its own sigma, which is not a measurement of a tilt but "
"of a shift of the beam centre", cal.tilt_significance) << std::endl;
} else {
// Distinct from the case above and worth saying so: there, the tilt was measured and found
// too small to believe; here the fit could not measure it at all, because the rings it had
// do not separate a tilt from a beam-centre shift in the first place.
std::cout << "Tilt: DECLINED and pinned - these rings cannot separate a tilt "
"from a shift of the beam centre at all" << std::endl;
}
const std::string poni_path = output_prefix + ".poni";
try {
WritePoniFile(poni_path, experiment, g);
} catch (const std::exception &e) {
logger.Error("{}", e.what());
return 1;
}
logger.Info("Wrote {}", poni_path);
return 0;
}
configure_offline_output(experiment, output_prefix);
// configure_offline_output applies the rugnux analysis defaults, one of which is the polarization
// factor, so an explicit --polarization has to land after it or it is silently overwritten.
if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value());
ClearStoredCrystal(experiment); // shared with the viewer; -S / -C below override it
if (space_group != nullptr)
experiment.SetSpaceGroup(*space_group);
experiment.ImagesPerTrigger(images_to_process);
// Re-determine the unit cell from scratch: discard any cell stored in the input file so
// indexing is not biased by it. A stale or wrong stored cell otherwise resolves the indexing
// algorithm to FFBIDX and drives it to the wrong lattice (e.g. a non-cubic cell for a cubic
// crystal). A user-supplied -C cell still takes effect (clears to nullopt when absent).
experiment.SetUnitCell(fixed_reference_unit_cell);
// --refine-geometry defaults ON for stills whenever a reference cell is available (-C or a
// reference MTZ): that is exactly when the geometry bundle-adjust can act (it anchors on a known
// cell) and where it helps weak/sparse stills. It is a no-op for rotation (which has its own
// two-pass) and for de-novo stills (no cell yet), so auto-enabling it only where it does something
// avoids spurious "skipping" warnings. Explicit --refine-geometry[=N] still forces it on;
// --refine-geometry=off opts out.
if (!refine_geometry.has_value() && !refine_geometry_disabled) {
const bool is_stills = !(experiment.GetGoniometer().has_value() && !force_still);
if (is_stills && experiment.GetUnitCell().has_value())
refine_geometry = 200;
}
experiment.MaxSpotCount(max_spot_count_override.value_or(RUGNUX_MAX_SPOT_COUNT));
if (max_spot_count_override.has_value())
logger.Info("Max spot count overridden to {}", max_spot_count_override.value());
// X-ray bandwidth: CLI overrides the value carried in the dataset; otherwise
// keep whatever the dataset provided (0 / none -> monochromatic).
if (bandwidth_fwhm)
experiment.BandwidthFWHM(bandwidth_fwhm);
if (experiment.GetBandwidthFWHM())
logger.Info("X-ray bandwidth FWHM set to {:.4f}", experiment.GetBandwidthFWHM().value());
// Rotation vs stills. A dataset collected on a rotation goniometer is processed as rotation data
// (two-pass indexing) by default; --force-still forces per-frame stills. The rotation flags
// (-R / --single-pass-rotation / --force-rotation-lattice) still request rotation explicitly and
// choose the pass/lattice; at this point they show up as rotation_indexing already being set.
// A stationary axis records where the head was, not that the crystal turned - processing such a
// dataset as rotation would leave the indexer with a zero stride and never index anything.
const bool has_goniometer = experiment.GetGoniometer().has_value()
&& experiment.GetGoniometer()->IsScanning();
if (force_still) {
if (rotation_indexing) {
logger.Error("--force-still conflicts with -R / --single-pass-rotation / --force-rotation-lattice");
exit(EXIT_FAILURE);
}
if (has_goniometer)
logger.Info("--force-still: treating the rotation dataset as independent stills");
} else if (!rotation_indexing && has_goniometer) {
rotation_indexing = true;
two_pass_rotation = true;
logger.Info("Dataset has a rotation goniometer axis: processing as rotation data (two-pass "
"indexing). Use --force-still to treat it as stills.");
}
// A stride cuts every rocking curve. The angles still come out right - the goniometer is shifted so
// one ordinal is one stride - but a reflection's partials are then sampled every k-th frame, so the
// combine sees a fraction of each event, the captured fraction collapses, and the partiality divides
// by a width the sweep never delivered. There is no reading of a strided rotation sweep worth having,
// so it is refused rather than processed into a plausible wrong answer.
if (rotation_indexing && image_stride != 1) {
logger.Error("-t/--stride {} cannot be used on rotation data: skipping frames cuts every "
"reflection's rocking curve, so the combined fulls and their partiality would be "
"measured over frames the sweep never recorded. Process the sweep whole, or select "
"a contiguous range with -s/-e. --mode azint and --force-still still take a stride.",
image_stride);
exit(EXIT_FAILURE);
}
// Scaling and merging are on by default (run_scaling initialised true); --no-merge turns them off
// for both rotation and stills, in which case only the per-image _process.h5 is written.
// Configure Indexing
IndexingSettings indexing_settings;
indexing_settings.Algorithm(indexing_algorithm);
if (min_indexed_spots) indexing_settings.ViableCellMinSpots(min_indexed_spots.value());
indexing_settings.RotationIndexing(rotation_indexing);
if (rotation_indexing_range.has_value())
indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value());
indexing_settings.GeomRefinementAlgorithm(refinement_algorithm);
indexing_settings.IndexIceRings(index_ice_rings);
if (fft_min_unit_cell_A) indexing_settings.FFT_MinUnitCell_A(fft_min_unit_cell_A.value());
experiment.ImportIndexingSettings(indexing_settings);
// --detect-ice-rings[=on|off] overrides the value carried in from the dataset (HDF5MetadataSource
// sets DetectIceRings from the master file's detect_ice_rings key); with no flag the dataset stands.
// Where the file says nothing at all, the default is the geometry's: on for rotation, off for
// stills. A rotation sweep sits on the same rings for the whole run, so ice there is a coherent
// systematic worth handling, and the ice-presence gate keeps it inert on a clean crystal; a serial
// stills run has too few spots per image to spend any of them on flagging.
if (detect_ice_rings.has_value())
experiment.DetectIceRings(detect_ice_rings.value());
else if (!dataset->file_detect_ice_rings.has_value())
experiment.DetectIceRings(rotation_indexing);
// Scale-fulls refits the per-frame scale on the rotation combined fulls; on by default for rotation
// data (where it lifts ISa substantially) and off for stills. --no-scale-fulls overrides.
const bool scale_fulls = scale_fulls_arg.value_or(rotation_indexing);
ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rotation_indexing);
scaling_settings.ScaleFulls(scale_fulls);
scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(scaling_settings.GetSmoothGDegrees()));
scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off
if (no_scaling_corrections)
scaling_settings.CorrectionSurfaces(false);
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
if (ice_min_score_arg)
scaling_settings.IceMinScore(static_cast<float>(*ice_min_score_arg));
if (ice_min_spot_ratio_arg)
scaling_settings.IceMinSpotRatio(static_cast<float>(*ice_min_spot_ratio_arg));
if (d_min_scale_merge)
scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
if (d_max_scale_merge)
scaling_settings.LowResolutionLimit_A(*d_max_scale_merge > 0.0f
? std::optional<double>(*d_max_scale_merge) : std::nullopt);
if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
scaling_settings.MergeFriedel(!anomalous_mode);
scaling_settings.MinPartiality(min_partiality);
// Drop edge-of-sweep truncated fulls (rocking curve captured < this fraction) from the rot3d combine.
// Defaults ON (0.7) for rotation - removes the low-capture fulls that inflate low-res R-meas and
// slightly bias accuracy; off for non-rot3d (no combine). 0.7 (rather than 0.5) also strips the
// partiality-extrapolated fulls that dominate the intensity second moment on weakly-diffracting
// crystals, so the de-novo space-group search is no longer starved by the error-model I/sigma floor
// (e.g. a weakly-diffracting F-cubic or hexagonal crystal recovers its true space group instead of
// P1). An explicit --min-captured-fraction wins.
scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or(scaling_settings.GetMinCapturedFraction()));
// Capture-aware systematic sigma defaults ON (1.0) for the rot3d combine - it down-weights the
// over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous
// substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins.
scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(scaling_settings.GetCaptureUncertaintyCoeff()));
scaling_settings.ForcedMosaicity(forced_mosaicity_arg);
scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is in percent; the setting is a fraction
// Rotation default: run the de-novo space-group search a second time on a merge of only the
// well-measured observations and keep whichever found more symmetry. It cannot lose symmetry - see
// Rugnux.cpp - so the cut being imperfect only means the second opinion contributes nothing.
scaling_settings.SearchMinZeta(search_min_zeta_arg.value_or(scaling_settings.GetSearchMinZeta()));
scaling_settings.OutlierRejectNsigma(
outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma()));
experiment.ImportScalingSettings(scaling_settings);
// Integration radii: r1 (signal box), r2/r3 (background annulus).
if (integration_radius_arg) {
std::vector<float> rr;
std::stringstream ss(*integration_radius_arg);
std::string tok;
while (std::getline(ss, tok, ',')) {
trim_in_place(tok);
if (!tok.empty())
rr.push_back(parse_number_arg<float>(tok.c_str(), "--integration-radius", logger,
0.1f, 1000.0f));
}
float r1, r2, r3;
if (rr.size() == 1) { r1 = rr[0]; r2 = r1 + 2.0f; r3 = r1 + 4.0f; }
else if (rr.size() == 3) { r1 = rr[0]; r2 = rr[1]; r3 = rr[2]; }
else { logger.Error("--integration-radius expects r1 or r1,r2,r3"); return 1; }
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.R1(r1).R2(r2).R3(r3);
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Integration radii set to r1={:.1f} r2={:.1f} r3={:.1f}", r1, r2, r3);
} else if (!rotation_indexing) {
// Stills spots span a range of crystal orientations captured in a single shot, so they land
// wider on the detector than the r1=4 monochromatic-rotation default assumes. A larger signal
// box lets the profile-fit integrator capture the whole spot while its profile weighting keeps
// the extra background from adding noise (a plain box-sum degrades with it). Measured R-free
// gains on serial stills. An explicit --integration-radius always wins.
//
// r3 = 14 rather than 12: the reported sigma carries the ring mean's own error as
// (dI/dbkg)^2 * bkg/n_ring, and the ring only holds 244 pixels at r3 = 12 against the
// rotation default's 408, so that term is a larger share of the variance here than
// anywhere else. 14 takes the ring to about 410 pixels. What normally stops a ring being
// widened is crowding - a neighbour inside the annulus is excluded from it - and stills
// barely crowd: under a percent of stills reflections have any neighbour within 13 px,
// against most of them on rotation.
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.R1(6.0f).R2(8.0f).R3(14.0f);
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=14.0 (override with --integration-radius)");
}
// Adaptive signal radius: measured in the pre-scan and applied there (Rugnux::PreScan). Radii the
// user set by hand always win, and stills already default to a wide box, so the rule is on for
// rotation only.
const bool adaptive_integration_radius =
!integration_radius_arg
&& adaptive_integration_radius_arg.value_or(rotation_indexing);
if (forced_prediction_mosaicity_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.ForcedPredictionMosaicity_deg(static_cast<float>(*forced_prediction_mosaicity_arg));
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Prediction window fixed at a rocking width of {:.4f} deg; the partiality still uses "
"the per-image sigma_M", *forced_prediction_mosaicity_arg);
}
if (integration_stencil_arg) {
const float k = parse_number_arg<float>(integration_stencil_arg->c_str(), "--integration-stencil",
logger, 0.0f, 10.0f);
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.StencilKSigma(k);
experiment.ImportBraggIntegrationSettings(bis);
if (!(k > 0.0f))
logger.Info("Background ring left circular (--integration-stencil 0)");
else if (!(experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f))
logger.Info("--integration-stencil {:.2f} has no effect without --bandwidth: "
"the ring is elongated by the beam's radial streak, which is zero here", k);
else
logger.Info("Background ring elongated by {:.2f} times the beam's radial streak, per reflection", k);
}
if (integrator_mode) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.Integrator(*integrator_mode);
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Integrator set to {}", *integrator_mode == IntegratorMode::BoxSum ? "box-sum"
: *integrator_mode == IntegratorMode::ProfileGaussian ? "profile (gaussian)"
: "profile (empirical)");
}
if (integration_d_min_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
// 0 spells "no limit" for the sibling resolution options, so it has to mean the same here.
bis.DMinLimit_A(*integration_d_min_arg > 0.0
? std::optional<float>(static_cast<float>(*integration_d_min_arg))
: std::nullopt);
experiment.ImportBraggIntegrationSettings(bis);
}
if (max_hkl_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.MaxHKL(static_cast<int>(*max_hkl_arg));
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Predicting reflections with |h|,|k|,|l| <= {} (overriding the per-crystal bound)", *max_hkl_arg);
}
if (background_trim_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.BackgroundTrimFraction(static_cast<float>(*background_trim_arg));
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip",
*background_trim_arg);
}
if (background_clip_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.BackgroundClipNSigma(static_cast<float>(*background_clip_arg));
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg);
} else if (!background_trim_arg && experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f) {
// Broadband (pink-beam / DMM) reflections are long and their wings reach into the background
// ring, so the default clip is tighter there. It is only a default - --background-clip and
// --background-trim override it as they do for any other data.
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.BackgroundClipNSigma(BROADBAND_BACKGROUND_CLIP_NSIGMA);
experiment.ImportBraggIntegrationSettings(bis);
logger.Info("Background ring: high-side clip at {:.1f} sigma (broadband default)",
BROADBAND_BACKGROUND_CLIP_NSIGMA);
}
if (overlap_arg || overlap_minpk_arg) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
if (overlap_arg)
bis.Overlap(*overlap_arg);
if (overlap_minpk_arg)
bis.OverlapMinPeak(static_cast<float>(*overlap_minpk_arg));
experiment.ImportBraggIntegrationSettings(bis);
const auto ovl = bis.GetOverlap();
logger.Info("Overlapping signal regions: {}", ovl == OverlapMode::Reject
? fmt::format("reject below {:.2f} of the profile", bis.GetOverlapMinPeak())
: (ovl == OverlapMode::Exclude ? std::string("exclude the shared pixels")
: std::string("off")));
}
if (background_radial_given) {
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
bis.BackgroundRadialCorrection(background_radial_arg);
experiment.ImportBraggIntegrationSettings(bis);
}
// Say what was decided whether or not the flag was given. In auto the correction is applied per
// image from that image's ice score, so a log that only spoke when asked left no record of whether
// it ran - and a run whose numbers differ from another's gave no way to see why.
{
const auto radial = experiment.GetBraggIntegrationSettings().GetBackgroundRadialCorrection();
logger.Info("Background ring: radial curvature correction {}",
radial.has_value() ? (*radial ? "on" : "off")
: "auto (per image, where the ice score sees a smooth ring)");
}
SpotFindingSettings spot_settings;
spot_settings.enable = true;
spot_settings.indexing = true;
spot_settings.signal_to_noise_threshold = sigma_spot_finding;
spot_settings.photon_count_threshold = photon_count_threshold_spot_finding;
// Detection defaults differ by workflow; each is overridden by its flag, which always wins.
// - min-pix: choosing it per image (unset) only means something where each frame is indexed on its
// own. Rotation indexing builds ONE lattice from all frames, so it keeps the fixed value.
// - adaptive detection: on by default for both workflows.
// The high-resolution limit is NOT one of them: unset means "as far as the detector reaches" for
// rotation as well as stills. Rotation used to keep 1.5 A on the strength of an indexing-rate
// measurement, but over the 33-crystal battery that limit changes nothing on 29 crystals, changes
// no space-group decision at all, and on the crystals where it does bite it is the LIMIT that is
// worse: the one crystal that loses appreciable indexing rate without it (99.5 -> 94.2%) comes back
// with better R_meas, better high-resolution CC1/2 and better ISa. Fewer frames, better data.
spot_settings.min_pix_per_spot = min_pix_per_spot;
if (rotation_indexing && !spot_settings.min_pix_per_spot.has_value())
spot_settings.min_pix_per_spot = 2;
spot_settings.adaptive_threshold = adaptive_spots.value_or(true);
spot_settings.high_resolution_limit = d_min_spot_finding;
spot_settings.false_pixels_per_frame = false_pixels_per_frame;
spot_settings.low_resolution_limit = d_max_spot_finding;
// Validate the assembled spot-finding settings the same way the online receivers do (broker and
// receiver call this same function). It enforces the cross-field constraints that per-argument
// bounds cannot express - in particular that the low-resolution limit is coarser than the
// high-resolution limit, so --spot-low-resolution below the high-res cut no longer silently
// rejects every pixel.
try {
DiffractionExperiment::CheckDataProcessingSettings(spot_settings);
} catch (const std::exception &e) {
logger.Error("Invalid spot-finding settings: {}", e.what());
return 1;
}
// Run the shared full-analysis workflow (rotation indexing + scaling/merging live in
// Rugnux; the experiment above carries all algorithm settings).
ProcessConfig config;
config.mode = ProcessMode::FullAnalysis;
config.start_image = start_image;
config.end_image = end_image;
config.stride = image_stride;
config.nthreads = nthreads;
config.nthreads_auto = nthreads_auto;
config.output_prefix = output_prefix;
config.spot_finding = spot_settings;
config.rotation_indexing = rotation_indexing;
config.two_pass_rotation = two_pass_rotation;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_search_pxl = beam_center_search;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
config.adaptive_integration_radius = adaptive_integration_radius;
config.rotation_postrefine_geometry = rotation_postrefine_geometry;
// Measure the spot budget from the data unless the user pinned it.
config.measure_spot_budget = !max_spot_count_override.has_value();
config.rotation_scale = rotation_scale;
config.rotation_indexing_image_count = rotation_indexing_image_count;
config.forced_rotation_lattice = forced_rotation_lattice;
config.refine_geometry = refine_geometry;
config.run_scaling = run_scaling;
config.scaling_iter = scaling_iter;
config.reference_data = reference_data;
config.reference_has_free_flags = reference_has_free_flags;
config.observation_dump_path = dump_observations;
config.export_unmerged = export_unmerged;
config.export_unmerged_partials = export_unmerged_partials;
config.write_p1_crosscheck = write_p1_crosscheck;
config.model_path = model_pdb;
// When merging, the merged reflections (.mtz/.cif) are the wanted output; skip the large
// _process.h5 unless explicitly requested. Without merging, the _process.h5 is the only output.
config.write_process_h5 = run_scaling ? write_process_h5_flag : true;
Rugnux process(reader, experiment, *dataset->pixel_mask, config);
g_active_process = &process;
std::signal(SIGINT, handle_sigint);
ProcessResult result;
try {
result = process.Run();
} catch (const std::exception &e) {
logger.Error("Processing failed: {}", e.what());
exit(EXIT_FAILURE);
}
g_active_process = nullptr;
// The results report, next to the .mtz/.cif/.hkl. Written on every run with an output prefix -
// including --no-merge, which still determined an indexing and geometry result worth recording.
// The wall time is the whole invocation, which is more than result.total_time_s: that one starts
// inside Rugnux::Run, so it counts neither opening the file nor setting up the analysis.
provenance.wall_time_s = std::chrono::duration<double>(
std::chrono::steady_clock::now() - invocation_start).count();
WriteResultReport(output_prefix, input_file, experiment, result, logger, provenance);
// The space-group search is rendered here (not in the library) so the viewer does not emit it on
// stdout and the CLI owns the format.
if (result.space_group_search.has_value())
std::cout << std::endl << SearchSpaceGroupResultToText(*result.space_group_search) << std::endl;
if (!result.merge_statistics_text.empty())
std::cout << std::endl << result.merge_statistics_text << std::endl;
// Report statistics. The time is the whole run - a rotation two-pass integrates everything twice,
// and the pre-scan runs before either pass, so the per-pass number understated it by 3-4x. The rate
// and throughput stay per-pass: they say how fast rugnux moves through images, which running a
// second pass does not change.
std::cout << fmt::format("Processing time: {:.2f} s", result.total_time_s) << std::endl;
std::cout << fmt::format("Total wall time: {:.2f} s", provenance.wall_time_s) << std::endl;
if (result.pass_count > 1)
std::cout << fmt::format(" last pass: {:.2f} s (of {} passes)",
result.processing_time_s, result.pass_count) << std::endl;
std::cout << fmt::format("Frame rate: {:.2f} Hz", result.frame_rate_hz) << std::endl;
std::cout << fmt::format("Total throughput:{:.2f} MB/s", result.throughput_MBs) << std::endl;
if (result.indexing_rate.has_value())
std::cout << fmt::format("Indexing rate: {:.2f}%", result.indexing_rate.value() * 100.0) << std::endl;
// Final one-line summary of the adopted crystal (whether de-novo determined or fixed with -S),
// so it is not buried in the space-group-search block (which is de-novo only) or only in the mmCIF.
// Only when something actually indexed: with a zero indexing rate the cell is whatever the lattice
// search happened to return and no reflection was measured on it, so printing it as the run's answer
// states a result the data do not support.
const bool anything_indexed = result.indexing_rate.value_or(0.0f) > 0.0f;
if (result.space_group.has_value() && anything_indexed) {
std::string line = fmt::format("Space group: {} (No. {})", result.space_group->short_name(),
result.space_group->number);
// Name every group the data cannot separate, not just the representative. Some pairs share
// their whole absence pattern - an enantiomorphic pair (P4_1 vs P4_3), or I23 vs I2_13 and
// I222 vs I2_12_12_1, where the screw condition h00: h=2n is already implied by the
// I-centering - so the choice among them is a convention, not a measurement. The search
// reports the representative as the lowest space-group number; saying so here keeps the
// summary from claiming a decision the diffraction did not make.
// Not the adopted group itself: where --model has taken the enantiomorph, the group now
// reported is the one the search listed as the alternative, and naming it as its own
// alternative would read as a decision that was not made.
bool named_alternative = false;
if (result.space_group_search.has_value())
for (const auto &alt : result.space_group_search->alternatives)
if (alt.number != result.space_group->number) {
line += fmt::format(" or {} (No. {})", alt.short_name(), alt.number);
named_alternative = true;
}
if (named_alternative)
line += " - indistinguishable from these data";
std::cout << line << std::endl;
}
if (result.consensus_cell.has_value() && anything_indexed) {
const auto &c = result.consensus_cell.value();
std::cout << fmt::format("Unit cell: a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
c.a, c.b, c.c, c.alpha, c.beta, c.gamma) << std::endl;
}
if (result.indexing_rate.has_value() && !anything_indexed)
std::cout << "No image indexed - no crystal lattice was determined from this dataset" << std::endl;
// Each stage timer measures wall time inside one worker, so it counts the time that worker spent
// BLOCKED on a contended resource - above all the GPU - as well as its own work. Dividing that by
// the worker count, as this used to, assumes every worker was busy for the whole loop; the loop
// in fact runs at a fraction of the workers it was given (a third of them on a large detector,
// less on a small one), so the printed cost came out low by that same factor.
//
// Report the share instead. A stage's fraction of a worker's own per-image time is what that
// stage is responsible for, whatever the contention was, and spending it against the loop's wall
// time per image gives a number that is attributable and that sums to the loop. The worker mean
// itself is printed at the end, because the gap between it and the wall is the waiting, and that
// is worth seeing rather than dividing away.
const auto &t = result.mean_processing_time;
const double loop_ms = (result.images_processed > 0 && result.image_loop_time_s > 0.0)
? result.image_loop_time_s * 1e3 / static_cast<double>(result.images_processed) : 0.0;
const double worker_mean_ms = std::isfinite(t.processing) ? t.processing * 1e3 : 0.0;
auto stage = [&](const char *name, float mean_s) {
// A stage that never ran has no mean at all - the per-image indexing and scaling timers are
// never fed on the two-pass rotation path, where the lattice is forced and the merge happens
// outside the image loop. Say nothing rather than printing nan.
return (std::isfinite(mean_s) && worker_mean_ms > 0.0)
? fmt::format(" {} {:.2f}", name, mean_s * 1e3 / worker_mean_ms * loop_ms) : std::string();
};
std::cout << fmt::format("Per-image cost (ms of loop wall, {} workers):",
result.image_loop_threads > 0 ? result.image_loop_threads : nthreads)
<< stage("decompress", t.compression) << stage("preprocess", t.preprocessing)
<< stage("azint", t.azint) << stage("spot-finding", t.spot_finding)
<< stage("indexing", t.indexing) << stage("refinement", t.refinement)
<< stage("indexing-analysis", t.indexing_analysis) << stage("prediction", t.bragg_prediction)
<< stage("integration", t.integration) << stage("scaling", t.image_scale)
<< stage("total", t.processing)
<< fmt::format(" [{:.1f} ms inside a worker; the difference is waiting]", worker_mean_ms)
<< std::endl;
// The stage timers only cover the per-image loop. On a rotation run the first-pass indexing and the
// scaling/merging sit outside it and can be a large share of the run, so report the loop against the
// whole run instead of leaving the difference unexplained. Both are the last pass only: a two-pass
// rotation run does all of this twice.
if (result.images_processed > 0 && result.image_loop_time_s > 0.0) {
const double loop_ms = result.image_loop_time_s * 1e3 / static_cast<double>(result.images_processed);
const double outside_s = result.processing_time_s - result.image_loop_time_s;
std::cout << fmt::format("Per-image wall: {:.2f} ms in the image loop ({:.2f} s); "
"{:.2f} s outside it (first-pass indexing, scaling/merging) [last pass]",
loop_ms, result.image_loop_time_s, std::max(0.0, outside_s)) << std::endl;
}
if (result.cancelled)
logger.Warning("Processing was cancelled after {} images", result.images_processed);
return EXIT_SUCCESS;
}
int main(int argc, char **argv) {
try {
return RunRugnux(argc, argv);
} catch (const std::exception &e) {
Logger("rugnux").Error("{}", e.what());
return EXIT_FAILURE;
}
}
Reference in New Issue View Git Blame Copy Permalink