diff --git a/common/Reflection.h b/common/Reflection.h index 4b18846af..9f054599a 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -76,6 +76,13 @@ struct ScaledFull { float I = NAN; float sigma = NAN; float image_number = NAN; // the frame of the event's peak partial, which tells two events of one index apart + // The full's counting variance as the merge rebuilds it at an intensity I, before its error model: + // var_bkg + var_per_I * I + (capture * I)^2, all at the intensity's scale (what a later merge rebuilds + // the full from). NaN where not exported. + float var_bkg = NAN; + float var_per_I = NAN; + float capture = NAN; + float d = NAN; // resolution, A }; struct MergedReflection { diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index 8b70c6ac9..8fc2cd2c9 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -315,3 +315,12 @@ ScalingSettings& ScalingSettings::ReportResolutionRange_A(std::optional ScalingSettings::GetReportResolutionRange_A() const { return report_resolution_range; } + +ScalingSettings& ScalingSettings::SweepSegments(std::vector input) { + sweep_segments = std::move(input); + return *this; +} + +const std::vector &ScalingSettings::GetSweepSegments() const { + return sweep_segments; +} diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index 6e3d0b8e7..d1b3af5da 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -3,7 +3,9 @@ #pragma once +#include #include +#include #include "JFJochException.h" // How the high-resolution cutoff for the written reflections and the reported shell table is chosen @@ -17,8 +19,20 @@ struct ReportResolutionRange { double d_max = 0.0; }; +// One sweep of several whose frames are laid end to end in one rotation merge (the joint merge of several +// sweeps): frames [first_frame, first_frame + n_frames). Every observation of a segment's frames is a full +// already combined in its own sweep (partiality 1, filed under the frame of its peak), so the merge +// combines nothing again, recomputes nothing from the rocking geometry, and lets nothing that runs along +// the frames - the per-frame scale's smoothness above all - reach across from one sweep into the next +// (RotationScaleMerge). +struct SweepSegment { + int32_t first_frame = 0; + int32_t n_frames = 0; +}; + class ScalingSettings { bool refine_wedge = false; + std::vector sweep_segments; bool merge_friedel = true; std::optional high_resolution_limit_A; @@ -150,6 +164,8 @@ public: ScalingSettings& ResolutionCCTarget(double input); ScalingSettings& ReportShellCount(int input); ScalingSettings& ReportResolutionRange_A(std::optional input); + // The sweeps a rotation merge's frames are made of; empty (the default): one sweep of partials. + ScalingSettings& SweepSegments(std::vector input); [[nodiscard]] bool GetRefineWedge() const; @@ -189,4 +205,5 @@ public: [[nodiscard]] double GetResolutionCCTarget() const; [[nodiscard]] int GetReportShellCount() const; [[nodiscard]] std::optional GetReportResolutionRange_A() const; + [[nodiscard]] const std::vector &GetSweepSegments() const; }; diff --git a/docs/RUGNUX_CLI.md b/docs/RUGNUX_CLI.md index 6081bc20d..382cd92fc 100644 --- a/docs/RUGNUX_CLI.md +++ b/docs/RUGNUX_CLI.md @@ -30,11 +30,10 @@ rugnux -o myrun dataset_master.h5 # NXmx / EIGER / Jungfraujoch mas rugnux -o myrun sweep_00001.cbf # CBF, marCCD or SMV sweep: name one frame rugnux -o myrun dataset_master.h5 model.cif # with a model (or --model model.cif) rugnux -o myrun x_master.h5 http://host:5232/ # follow a broker's collection while it is written -rugnux -o myrun --sweeps separate a_master.h5 b_master.h5 # several sweeps, each alone +rugnux -o myrun a_master.h5 b_master.h5 # several sweeps of one crystal, merged as one ``` -Chosen when the file's goniometer axis turns. Several sweeps can be given (below, Inputs); for now -each is processed alone with `--sweeps separate`. A model is placed (rigid body, or +Chosen when the file's goniometer axis turns. Several sweeps can be given (below, Inputs). A model is placed (rigid body, or molecular replacement), refined and validated after merging; it does not change the processing decisions ([Models and reference data](RUGNUX_MODEL_AND_REFERENCE.md)). @@ -105,9 +104,19 @@ template (the name with its frame number taken out - any frame of a sweep names `rugnux a_*.cbf` is the one sweep `a`), one per rotation master. Only the inputs given are processed; other sweeps beside them are named in the log. A mix of rotation and stills masters is refused, and so are sweeps at different wavelengths, `-s`/`-e` and `--force-still` with several sweeps. -`--sweeps separate` processes each sweep alone, as if it were the only input, into -`_s_*` (k = 1, 2, ... in the order given); processing several sweeps as one dataset -(`--sweeps joint`, the default) is not implemented yet and refuses. +By default (`--sweeps joint`) the sweeps are one dataset, merged in two levels as XSCALE merges the +data sets XDS made of each sweep: each sweep is processed on its own first - indexed, integrated, its +partials combined and scaled, its space group determined - writing only its log `_s_log.txt`; +the sweeps are then put on the axes of the one with the most reflections (the change of basis, an +axis permutation or a twin law of the lattice, under which their merged intensities correlate best; +a sweep correlating below 0.5 on every axes is left out and named in the log), and their fulls are +scaled and merged together: a per-frame scale within each sweep, the sweeps' dose order read from +their relative B factors (not from the input order or the header times), a decay along that order or +one B per sweep, whichever predicts held-out frames better, then the space-group search, the +resolution cut and everything after it on the joint data, into `.*`. The detector-dependent +corrections (absorption and detector surfaces) are each sweep's own; the joint merge writes no +unmerged MTZ yet. `--sweeps separate` processes each sweep alone instead, as if it were the only +input, into `_s_*` (k = 1, 2, ... in the order given). What each format needs is in [What Rugnux reads](RUGNUX_FORMATS.md); serial XFEL stills are read without a geometry file where their own metadata places the detector @@ -142,7 +151,7 @@ With the input already on disk, the URL is ignored. | `-v, --verbose` | Show the whole log on the terminal in place of the steps; `_log.txt` is written either way | | `--simulate-live ` | Read a stored input as the collection it was while being written: image *i* becomes readable (*i*+1)·`ms` ms after the start (0: all at once). The run takes the live path a broker URL takes, with the same result; for testing the live workflow without a broker. `--mode mx` only | | `--fast-live` | With `--simulate-live` or a broker URL: run as the broker's rugnux live mode does — every decision stands where it was first taken, nothing is re-checked or run again once the sweep is in, the integration radius is fixed, and SIGINT cuts the run to what it has established. A difficult crystal can come out worse | -| `--sweeps joint\|separate` | Several rotation sweeps (see Inputs): `separate` processes each alone into `_s_*`; `joint` (the default), one dataset of all of them, is not implemented yet | +| `--sweeps joint\|separate` | Several rotation sweeps (see Inputs): `joint` (the default) merges them as one dataset into `.*`; `separate` processes each alone into `_s_*` | | `--no-gpu` | Run on the CPU only, as on a machine without a GPU, even where one is present: CUDA is not touched (nor, on macOS, the Apple GPU that Bragg integration, the model density and the molecular-replacement rotation function otherwise run on) and every GPU step takes its CPU path (`-X FFT` indexes with FFTW and `-X FFBIDX` with its CPU implementation; `-X cuFFT` is refused) | ### Mode diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 2848ca1c4..a88ab71ad 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -365,6 +365,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, // Decay + absorption correction surfaces are one master toggle (on by default; both cross-validated, // so a no-op when their systematic is absent). refine_decay_b = s.GetCorrectionSurfaces(); + segments = s.GetSweepSegments(); absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0; relative_b_deg = s.GetRelativeBDegrees(); // per-batch relative-B batch width (0 = off) @@ -6113,18 +6114,31 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // reflections are) - each with the sigma the merge weighted it by. if (export_scaled_fulls && full_stats && !for_search) { // Which fulls, in order, and then each one's record on its own - on all threads. + // The fulls of one sweep of several (uncut) also hand back the reflections the group extinguishes: + // a joint merge searches its own space group, and the screw axes are decided on exactly these. + const gemmi::GroupOps ops = Group().operations(); + auto extinguished = [&](const Obs &o) { + return export_scaled_fulls_uncut && o.group < 0 && frame_cell_ok[o.frame] && o.corr > 0.0f + && std::isfinite(o.corr) && std::isfinite(o.I) && o.d > 0.0f + && ops.is_systematically_absent({{o.h, o.k, o.l}}); + }; const std::vector kept = KeptIndices(n_full, nthreads, [&](int i) { const int g = mf.group[i]; - return g >= 0 && !rejected_obs[i] && std::isfinite(merged_I[g]) - && !(effective_d_min && acc[g].d < *effective_d_min); + if (g < 0) + return extinguished(fulls[i]); + return !rejected_obs[i] && std::isfinite(merged_I[g]) + && (export_scaled_fulls_uncut || !(effective_d_min && acc[g].d < *effective_d_min)); }); result.scaled_fulls.resize(kept.size()); ParallelChunks(static_cast(kept.size()), ThreadsForWork(kept.size(), nthreads), [&](int lo, int hi) { for (int j = lo; j < hi; ++j) { const Obs &o = fulls[kept[j]]; const float I_corr = o.I * o.corr; - result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr), - o.image_number}; + const float sigma = o.group >= 0 ? corrected_sigma(o, I_corr, o.sigma * o.corr) + : model_sigma(o, I_corr, o.sigma * o.corr); + result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, sigma, + o.image_number, o.corr * o.corr * o.var_bkg, o.corr * o.var_per_I, + o.capture, o.d}; } }); } diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index f8886d385..d1ebd935a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -81,17 +81,6 @@ std::pair TrimBlankSweepEnds(std::vector &outcomes const std::vector &tally); std::pair TrimBlankSweepEnds(std::vector &outcomes, size_t nthreads); -// One sweep of several whose frames are laid end to end in one engine (the joint merge of several -// sweeps): frames [first_frame, first_frame + n_frames) of the outcomes. Every observation of a segment's -// frames is a full already combined in its own sweep (partiality 1, filed under the frame of its peak), -// so nothing is combined again, nothing is recomputed from the rocking geometry, and nothing that runs -// along the frames - the per-frame scale's smoothness above all - reaches across from one sweep into the -// next. -struct SweepSegment { - int32_t first_frame = 0; - int32_t n_frames = 0; -}; - class RotationScaleMerge { public: // The log scale x minimising sum_f J_f (x_f - y_f)^2 + lambda * sum_f (x_f - 2 x_{f+1} + x_{f+2})^2: @@ -310,12 +299,16 @@ public: // Merge only the first n frames of the outcomes (by default all of them). Set before Ingest(). void SetFrameLimit(int n) { frame_limit = n; } - // The outcomes are several sweeps of pre-combined fulls laid end to end (SweepSegment). Set before - // Ingest(); the segments must cover the frames in order. Such an engine runs on the host: its fulls - // are a few per frame where a sweep's partials are millions, and the device pipeline is the partials'. - // Without segments (the default) nothing of this runs. + // The outcomes are several sweeps of pre-combined fulls laid end to end (SweepSegment; by default the + // experiment's ScalingSettings::GetSweepSegments). Set before Ingest(); the segments must cover the + // frames in order. Such an engine runs on the host: its fulls are a few per frame where a sweep's + // partials are millions, and the device pipeline is the partials'. Without segments nothing of this runs. void SetSweepSegments(std::vector s) { segments = std::move(s); } + // Whether the fulls handed back (SetExportScaledFulls) are all the merge kept, past its resolution cut + // too: the fulls of one sweep of several, which the joint merge cuts again on its own. + void SetExportScaledFullsUncut(bool on) { export_scaled_fulls_uncut = on; } + // The GPU the engine runs on (device 0 by default). Set before Ingest(). The device only says where // the merge runs: on cards of one model the result is the same on any of them. void SetGpuDevice(int device) { gpu_device = device; } @@ -763,6 +756,7 @@ private: std::vector &lambda, bool choose) const; int gpu_device = 0; // see SetGpuDevice bool export_scaled_fulls = false; // see SetExportScaledFulls + bool export_scaled_fulls_uncut = false; // see SetExportScaledFullsUncut bool french_wilson = true; // see SetFrenchWilson // --- helpers (each a flat pass; see the .cpp) --- diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 9f1b96385..3a5b90302 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -17,6 +17,8 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxFirstPassRescues.cpp RugnuxFirstPassAccept.cpp RugnuxScaleMerge.cpp + RugnuxJointSweeps.cpp + RugnuxJointSweeps.h RugnuxDefaults.cpp RugnuxDefaults.h RugnuxLive.cpp diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 33de9fb4c..dfbc0260c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -787,6 +787,21 @@ std::vector Rugnux::PrepassSubset(int prepass_images) const { return keep; } +ProcessResult Rugnux::RunJoint(std::vector outcomes, const UnitCell &cell, + gemmi::CrystalSystem system, char centring, RugnuxObserver *observer) { + const auto run_start = std::chrono::steady_clock::now(); + joint_outcomes_ = std::move(outcomes); + joint_cell_ = cell; + joint_system_ = system; + joint_centring_ = centring; + decide_late_ = true; + ProcessResult result = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + joint_outcomes_.reset(); + result.pass_count = 1; + result.total_time_s = std::chrono::duration(std::chrono::steady_clock::now() - run_start).count(); + return result; +} + ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) { Logger logger("Rugnux"); ProcessResult result; @@ -818,7 +833,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (config_.stride <= 0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive"); - const auto total_images_in_file = static_cast(reader_.GetNumberOfImages()); + const bool joint = joint_outcomes_.has_value(); + const auto total_images_in_file = joint ? static_cast(joint_outcomes_->size()) + : static_cast(reader_.GetNumberOfImages()); int end_image = config_.end_image; if (end_image < 0 || end_image > total_images_in_file) end_image = total_images_in_file; @@ -907,7 +924,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // After any geometry refinement, so the shadow is found about the beam centre actually used, // and before the azimuthal mapping and the output mask below, which both read pixel_mask_. - if (config_.detect_beam_stop.has_value() || config_.adaptive_integration_radius) + if (!joint && (config_.detect_beam_stop.has_value() || config_.adaptive_integration_radius)) PreScan(start_image, early_images, std::min(config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES), HotPixelFinder::MAX_SAMPLED_FRAMES), observer); @@ -1120,6 +1137,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // (see the reflection mask in BraggIntegrationEngineCPU). A fixed group used to cost real // intensity accuracy for that reason. indexer->PredictAllCentringNodes(true); + if (joint) { + indexer->GetIntegrationOutcome() = std::move(*joint_outcomes_); + joint_outcomes_.reset(); + } } const auto start_time = std::chrono::steady_clock::now(); @@ -1156,7 +1177,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b .scan_result = scan_result, .unmerged_mtz = unmerged_mtz}; // First pass of two-pass rotation indexing (full analysis only). - if (indexing && force_rotation_result_.has_value()) { + if (joint) { + // The sweeps were indexed, integrated and combined each on its own; their fulls are the frames. + result.rotation_lattice_found = true; + result.indexing_rate = 1.0f; + result.images_processed = static_cast(images_to_process); + result.consensus_cell = joint_cell_; + result.consensus_centering = joint_centring_; + end_msg.rotation_lattice = CrystalLattice(*joint_cell_); + end_msg.rotation_lattice_type = LatticeMessage{.centering = joint_centring_, .niggli_class = 0, + .crystal_system = joint_system_}; + result.rotation_lattice_type = end_msg.rotation_lattice_type; + end_msg.unit_cell = result.consensus_cell; + } else if (indexing && force_rotation_result_.has_value()) { // Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation / // search metadata / axis), not just its lattice - a lattice-only force loses that metadata and // integrates markedly worse. @@ -1194,7 +1227,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // // Only an output pass probes: the answer the run returns is always one (RunAllPasses), and no other // pass's probe is read - the pre-pass, the beam-centre arm and the probe passes leave theirs unread. - if (indexing && indexer && experiment_.IsRotationIndexing() && write_output && !geometry_prepass) { + if (!joint && indexing && indexer && experiment_.IsRotationIndexing() && write_output && !geometry_prepass) { std::vector frames; for (const int o : select_equally_spaced_image_ordinals(images_to_process, 60)) frames.push_back(o); @@ -1205,7 +1238,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::optional gpu_image_pass; // what a probe that does not fit beside it waits for if (!gpu_work_beside_) gpu_image_pass.emplace(); - ImagePass(p); + if (!joint) + ImagePass(p); gpu_image_pass.reset(); image_pass_mark.Stop(); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index f240ce0f3..2f031f189 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -114,6 +114,11 @@ struct ProcessConfig { // scaling. bool write_merged = true; + // One sweep of several merged as one (--sweeps joint): the run is made as on its own, but its written + // merge hands back all of its fulls, past the resolution cut too (MergedOutput::scaled_fulls), for the + // joint merge to scale and cut again (Rugnux::RunJoint). + bool sweep_member = false; + SpotFindingSettings spot_finding; // FullAnalysis spot finding (and Calibration by spots) // Calibration: how the powder rings are measured, and the calibrant's rings as q = 2*pi/d (1/A). @@ -1261,6 +1266,11 @@ class Rugnux { // Whether the first passes decide their near ties on all the frames: offline, on a re-check, and on // the run made again. bool decide_late_ = false; + // RunJoint's frames and cell, for the one pass it makes (RunPipeline); empty otherwise. + std::optional> joint_outcomes_; + std::optional joint_cell_; + gemmi::CrystalSystem joint_system_ = gemmi::CrystalSystem::Triclinic; + char joint_centring_ = 'P'; // Whether what the run otherwise decides on the pre-pass frames - the rotation-scale walk's refits, // the demoted metric-symmetry arm, the integration radius (RefinedPasses) - is decided on all the // frames: only on the run made again. Offline too these are decided on the pre-pass frames, and a near @@ -1656,6 +1666,14 @@ public: // Runs the configured workflow to completion or until Cancel(). Throws on setup failure. ProcessResult Run(RugnuxObserver *observer = nullptr); + // The joint merge of several sweeps (--sweeps joint): `outcomes` are the frames of all of them laid + // end to end, each holding the fulls its own sweep combined (ScalingSettings::SweepSegments says which + // frames are whose), indexed on the axes of `cell`. Nothing is read from the images: the merge, the + // space-group search and everything after it run on these as on the outcomes of one sweep. + // `system` and `centring` are those of the lattice the cell is the conventional cell of. + ProcessResult RunJoint(std::vector outcomes, const UnitCell &cell, + gemmi::CrystalSystem system, char centring, RugnuxObserver *observer = nullptr); + // The pinned host buffers the engines of a run kept for each other (CUDAMemHelpers.h), handed back to the // driver. Not part of Run: a process that outlives the run (the broker, the viewer) calls it once the // results are out, as it takes a tenth of a second the user would otherwise wait for; one that ends diff --git a/rugnux/RugnuxJointSweeps.cpp b/rugnux/RugnuxJointSweeps.cpp new file mode 100644 index 000000000..09898d4c8 --- /dev/null +++ b/rugnux/RugnuxJointSweeps.cpp @@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RugnuxJointSweeps.h" + +#include +#include + +#include "../common/CrystalLattice.h" +#include "../image_analysis/scale_merge/ReindexAmbiguity.h" + +JointSweepRelation RelateJointSweeps(const std::vector &sweeps) { + JointSweepRelation out; + out.sweeps.resize(sweeps.size()); + if (sweeps.empty()) + return out; + out.reference = 0; + for (size_t k = 1; k < sweeps.size(); ++k) + if (sweeps[k].fulls.size() > sweeps[out.reference].fulls.size()) + out.reference = static_cast(k); + const JointMemberSweep &ref = sweeps[out.reference]; + out.sweeps[out.reference] = {true, gemmi::Op::identity(), 1.0, 1.0, {}}; + for (size_t k = 0; k < sweeps.size(); ++k) { + if (static_cast(k) == out.reference) + continue; + JointSweepRelation::Sweep &s = out.sweeps[k]; + if (sweeps[k].merged.empty()) { + s.why_left_out = "it merged nothing"; + continue; + } + const ReferenceSettingChoice choice = ChooseReferenceSetting(sweeps[k].merged, sweeps[k].cell, + sweeps[k].space_group, ref.merged, ref.cell, + ref.space_group); + s.cob = choice.cob; + s.cc = choice.cc; + s.identity_cc = choice.identity_cc; + s.included = std::isfinite(choice.cc) && choice.cc >= MIN_REFERENCE_CC; + if (!s.included) + s.why_left_out = "its intensities correlate with the reference sweep's on no axes (best CC " + + std::to_string(choice.cc) + ")"; + } + return out; +} + +JointMergeInput BuildJointMergeInput(const std::vector &sweeps, + const JointSweepRelation &relation) { + JointMergeInput out; + const CrystalLattice lattice(sweeps[relation.reference].cell); + int32_t first = 0; + for (size_t k = 0; k < sweeps.size(); ++k) { + const auto &sweep = sweeps[k]; + const auto &rel = relation.sweeps[k]; + if (!rel.included || sweep.n_frames <= 0) + continue; + out.segments.push_back({first, static_cast(sweep.n_frames)}); + out.segment_sweep.push_back(static_cast(k)); + out.outcomes.resize(first + sweep.n_frames); + for (int64_t f = 0; f < sweep.n_frames; ++f) + out.outcomes[first + f].latt = lattice; + for (const ScaledFull &full : sweep.fulls) { + if (!std::isfinite(full.I) || !std::isfinite(full.var_bkg) || !(full.var_per_I > 0.0f) + || !(full.d > 0.0f)) + continue; + const auto frame = std::clamp(std::llround(full.image_number), 0, sweep.n_frames - 1); + const gemmi::Op::Miller hkl = rel.cob.apply_to_hkl({{full.h, full.k, full.l}}); + // A full as the merge takes a pre-combined observation, so that the combine rebuilds the variance + // its own sweep's merge gave it: a count c of intensity I/c at correction c = var_per_I is a + // variance c^2 var_bkg' + c I at I, which is the sweep's var_bkg + var_per_I I with + // var_bkg' = var_bkg / c^2 (the capture term folded in at the full's own intensity). + const float c = full.var_per_I; + const float capture_var = std::isfinite(full.capture) ? full.capture * full.I * full.capture * full.I : 0.0f; + Reflection r{}; + r.h = hkl[0]; r.k = hkl[1]; r.l = hkl[2]; + r.image_number = static_cast(first + frame); + r.d = full.d; + r.I = full.I / c; + r.var_bkg = (full.var_bkg + capture_var) / (c * c); + r.sigma = std::sqrt(std::max(r.var_bkg + std::max(0.0f, r.I), 1e-12f)); + r.prescaling_corr = c; + r.partiality = 1.0f; + r.zeta = 1.0f; + r.image_scale_corr = c; + r.observed = true; + out.outcomes[first + frame].reflections.push_back(r); + } + first += static_cast(sweep.n_frames); + } + return out; +} diff --git a/rugnux/RugnuxJointSweeps.h b/rugnux/RugnuxJointSweeps.h new file mode 100644 index 000000000..d66379830 --- /dev/null +++ b/rugnux/RugnuxJointSweeps.h @@ -0,0 +1,58 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../common/Reflection.h" +#include "../common/ScalingSettings.h" +#include "../common/UnitCell.h" +#include "../gemmi_gph/gemmi/symmetry.hpp" +#include "../image_analysis/IntegrationOutcome.h" + +// Several sweeps of one crystal merged as one (--sweeps joint), in two levels: each sweep is processed on +// its own - indexed, integrated, its partials combined into fulls and scaled - and the fulls of all of them +// are then scaled and merged together (RotationScaleMerge with SweepSegments), as XSCALE merges the data +// sets XDS made of each sweep. Between the two the sweeps are put on one set of axes: the indexing each +// was given is its own, and an ambiguous lattice or a different choice of axes indexes the same crystal +// differently from one sweep to the next. + +// What the joint merge takes from one sweep's own run. +struct JointMemberSweep { + std::string label; + int64_t n_frames = 0; + UnitCell cell{}; // as its reflections are written + gemmi::SpaceGroup space_group; // the group its run determined + std::vector merged; + std::vector fulls; // all its fulls, unmerged, on the same axes +}; + +// How each sweep is put on the axes of the reference - the sweep with the most fulls: the change of basis +// (an axis permutation or a twin law of the lattice, ChooseReferenceSetting) under which its merged +// intensities correlate best with the reference's, and that correlation. A sweep whose best correlation is +// below MIN_REFERENCE_CC is not the same crystal form on any axes, and is left out. +struct JointSweepRelation { + struct Sweep { + bool included = false; + gemmi::Op cob = gemmi::Op::identity(); + double cc = 0.0; // with the reference, under cob + double identity_cc = 0.0; // ... and on its own axes + std::string why_left_out; + }; + int reference = -1; + std::vector sweeps; +}; +JointSweepRelation RelateJointSweeps(const std::vector &sweeps); + +// The joint merge's input: the frames of the included sweeps laid end to end, in the order given, each +// holding its sweep's fulls on the reference's axes as pre-combined observations, and the segments that +// say whose frames are whose. +struct JointMergeInput { + std::vector outcomes; + std::vector segments; + std::vector segment_sweep; // the sweep each segment is +}; +JointMergeInput BuildJointMergeInput(const std::vector &sweeps, + const JointSweepRelation &relation); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 94e3e27cd..c29d51f29 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1688,6 +1688,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // Only a merge that may be written hands its fulls back: they are a copy as long as the // fulls, and the search merges and the pre-pass never write a .hkl. rsm->SetExportScaledFulls(!for_search && !geometry_prepass && write_output && config_.write_merged); + rsm->SetExportScaledFullsUncut(config_.sweep_member); // Amplitudes are made once, on the merge that is written, where the anisotropy tensor // their prior uses is known (see below) - not on every merge the tail makes on its way. rsm->SetFrenchWilson(false); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 1f58d54fa..14f639fcc 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -86,6 +86,7 @@ #include "RugnuxConsole.h" #include "RugnuxInputKind.h" #include "RugnuxSweeps.h" +#include "RugnuxJointSweeps.h" #ifdef __APPLE__ #include "../image_analysis/bragg_integration/BraggIntegrationEngineMetal.h" #endif @@ -123,7 +124,7 @@ never a directory; a model or reflection file among them is recognised by its co Rotation data (--mode mx, the default; chosen when the file's goniometer axis turns) rugnux -o myrun dataset_master.h5 NXmx / EIGER / Jungfraujoch master file rugnux -o myrun sweep_00001.cbf CBF, marCCD or SMV sweep: name one frame - rugnux -o myrun --sweeps separate a_master.h5 b_master.h5 several sweeps, each alone + rugnux -o myrun a_master.h5 b_master.h5 several sweeps of one crystal, merged as one rugnux -o myrun dataset_master.h5 model.cif with a model (or --model model.cif) rugnux -o myrun x_master.h5 http://host:5232/ follow a broker's collection while it is written A model is placed and refined after merging; it does not change the processing. @@ -172,7 +173,7 @@ Basic options --no-merge No scaling or merging; write the per-image _process.h5 --write-process-h5 Also write the (large) _process.h5 when merging --sweeps separate Several rotation sweeps (masters or frame templates): each alone, - into _s_*; processing them as one dataset is to come + into _s_*, instead of merged as one dataset -v, --verbose Whole log on the terminal (it is always in _log.txt) All options (about 120, for unusual data and diagnosis), each with what it does: @@ -960,6 +961,8 @@ struct RugnuxCli { std::string broker_url; std::shared_ptr model; RugnuxSweepsMode sweeps_mode = RugnuxSweepsMode::Joint; // --sweeps + // A log file for the whole run is open already (--sweeps joint): an input's own is not opened over it. + bool log_open = false; RugnuxCli(bool no_gpu_, std::chrono::steady_clock::time_point invocation_start_, RunProvenance provenance_, std::string metal_device_, std::filesystem::file_time_type files_since_, @@ -971,6 +974,8 @@ struct RugnuxCli { // The command line: RUGNUX_CONTINUE, or the exit code of a run that ends here. int Parse(int argc, char **argv); + // --sweeps joint: the sweeps one dataset (RugnuxJointSweeps.h). + int RunJointSweeps(const RugnuxSweepPlan &plan); // Each sweep of the command line, or its one input. int Run(); }; @@ -999,6 +1004,18 @@ struct RugnuxInputRun : RugnuxCli { std::optional sink; // The log file of this input, taken out of the log again once it is done. std::shared_ptr log_file_sink; + // One sweep of a joint run (--sweeps joint): processed as on its own, its result kept here and nothing + // written. + bool sweep_member = false; + std::optional member_result; + // The joint merge of the sweeps' fulls, made in place of this input's own run (Rugnux::RunJoint). + struct JointRequest { + std::vector outcomes; + UnitCell cell{}; + gemmi::CrystalSystem system = gemmi::CrystalSystem::Triclinic; + char centring = 'P'; + }; + std::optional joint_request; // Reads the input and sets up the experiment and the run's configuration: RUGNUX_CONTINUE, or the // exit code of a mode that ends here (refine, azint, calibration, scale). @@ -1764,7 +1781,7 @@ int RugnuxInputRun::PrepareRugnuxInput() { // The log file, next to the outputs, and everything logged so far into it. std::string &log_file = g_log_file; - if (!output_prefix.empty()) { + if (!output_prefix.empty() && !log_open) { log_file = output_prefix + "_log.txt"; try { auto file = std::make_shared(log_file, /*truncate=*/true); @@ -3600,7 +3617,9 @@ int RugnuxInputRun::RunPreparedInput() { ProcessResult result; try { - result = process.Run(console.get()); + result = joint_request ? process.RunJoint(std::move(joint_request->outcomes), joint_request->cell, + joint_request->system, joint_request->centring, console.get()) + : process.Run(console.get()); } catch (const NoLatticeFound &e) { // No crystal in these frames is a result, and the report says so - VERDICT= FAILED with // NO_LATTICE - before the run exits as failed. @@ -3638,6 +3657,10 @@ int RugnuxInputRun::RunPreparedInput() { throw RugnuxInputFailed{EXIT_FAILURE}; } g_active_process = nullptr; + if (sweep_member) { + member_result = std::move(result); + return EXIT_SUCCESS; + } // Where the modules are, for a later run's --module-geometry (and any NXmx reader). if (!output_prefix.empty() && result.modules_placed) { @@ -3880,11 +3903,8 @@ int RugnuxCli::Run() { logger.Error("--force-still processes one rotation sweep as stills; give the sweeps one at a time"); return EXIT_FAILURE; } - if (sweeps_mode == RugnuxSweepsMode::Joint) { - logger.Error("{} sweeps: processing them as one dataset is not implemented yet. --sweeps separate " - "processes each sweep alone, into _s_*", plan.sweeps.size()); - return EXIT_FAILURE; - } + if (sweeps_mode == RugnuxSweepsMode::Joint) + return RunJointSweeps(plan); // --sweeps separate: each sweep alone, as if it had been the one input, into _s. int failed = 0; @@ -3917,6 +3937,118 @@ int RugnuxCli::Run() { } +int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) { + // The whole run's log, from the start: each sweep writes its own beside it, and the joint merge goes on + // in this one. + if (!output_prefix.empty()) { + g_log_file = output_prefix + "_log.txt"; + try { + auto file = std::make_shared(g_log_file, /*truncate=*/true); + for (const auto &m : early_log->last_raw()) + file->log(m); + log_sinks->add_sink(file); + log_open = true; + } catch (const std::exception &e) { + logger.Warning("The log file {} cannot be written: {}", g_log_file, e.what()); + g_log_file.clear(); + } + } + // Each sweep processed on its own first, writing nothing, its fulls kept: the first level. + std::vector members; + std::vector inputs; + for (size_t k = 0; k < plan.sweeps.size(); k++) { + RugnuxInputRun sweep(*this); + sweep.input_files = {plan.sweeps[k].input}; + sweep.log_open = false; + // A sweep of a joint run writes its log, and nothing else (sweep_member). + if (!output_prefix.empty()) + sweep.output_prefix = output_prefix + "_s" + std::to_string(k + 1); + sweep.sweep_member = true; + logger.Info("Sweep {} of {}: {}, processed on its own", k + 1, plan.sweeps.size(), plan.sweeps[k].label); + int rc; + try { + rc = sweep.PrepareRugnuxInput(); + if (rc == RUGNUX_CONTINUE) { + sweep.config.sweep_member = true; + sweep.config.model.reset(); + sweep.config.force_mr = false; + sweep.config.image_sink = nullptr; + rc = sweep.RunPreparedInput(); + } + } catch (const RugnuxInputFailed &f) { + rc = f.code; + } + if (sweep.log_file_sink) + log_sinks->remove_sink(sweep.log_file_sink); + JointMemberSweep m; + m.label = plan.sweeps[k].label; + if (rc == EXIT_SUCCESS && sweep.member_result && sweep.member_result->outputs + && sweep.member_result->outputs->merged) { + const MergedOutput &merged = *sweep.member_result->outputs->merged; + m.n_frames = static_cast(sweep.member_result->images_processed); + m.cell = merged.cell; + m.space_group = merged.experiment.GetSpaceGroupOrP1(); + m.merged = merged.reflections; + m.fulls = merged.scaled_fulls; + logger.Info("Sweep {}: {} frames, {} in {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, {} fulls, {} merged", + k + 1, m.n_frames, m.space_group.xhm(), m.cell.a, m.cell.b, m.cell.c, m.cell.alpha, + m.cell.beta, m.cell.gamma, m.fulls.size(), m.merged.size()); + } else + logger.Warning("Sweep {} ({}) gave no merge and is left out of the joint one", k + 1, plan.sweeps[k].label); + members.push_back(std::move(m)); + inputs.push_back(plan.sweeps[k].input); + } + + // On one set of axes, and laid end to end: the second level's input. + const JointSweepRelation relation = RelateJointSweeps(members); + if (relation.reference < 0 || members[relation.reference].fulls.empty()) { + logger.Error("No sweep merged on its own, so there is nothing to merge jointly"); + return EXIT_FAILURE; + } + for (size_t k = 0; k < members.size(); k++) { + const auto &r = relation.sweeps[k]; + if (static_cast(k) == relation.reference) + logger.Info("Sweep {}: the reference the others are put on the axes of", k + 1); + else if (r.included) + logger.Info("Sweep {}: put on the reference's axes by {} (CC {:.3f}, {:.3f} as indexed)", k + 1, + r.cob.as_hkl().triplet('h'), r.cc, r.identity_cc); + else + logger.Warning("Sweep {} is left out of the joint merge: {}", k + 1, + r.why_left_out.empty() ? std::string("it gave no merge") : r.why_left_out); + } + JointMergeInput input = BuildJointMergeInput(members, relation); + if (input.segments.size() < 2) + logger.Warning("Only {} sweep goes into the joint merge", input.segments.size()); + + // The joint merge, in the reference sweep's setup, writing what a run writes. + const JointMemberSweep &ref = members[relation.reference]; + RugnuxInputRun joint(*this); + joint.input_files = {inputs[relation.reference]}; + if (log_open) + g_log_file = output_prefix + "_log.txt"; + int rc; + try { + rc = joint.PrepareRugnuxInput(); + if (rc != RUGNUX_CONTINUE) + return rc; + ScalingSettings s = joint.experiment.GetScalingSettings(); + s.SweepSegments(input.segments); + joint.experiment.ImportScalingSettings(s); + joint.config.rotation_postrefine_geometry = false; + joint.config.export_unmerged = false; + joint.config.export_unmerged_partials = false; + joint.config.image_sink = nullptr; + joint.joint_request = RugnuxInputRun::JointRequest{std::move(input.outcomes), ref.cell, + ref.space_group.crystal_system(), + ref.space_group.centring_type()}; + logger.Info("Joint merge of {} sweeps, {} frames", input.segments.size(), joint.joint_request->outcomes.size()); + rc = joint.RunPreparedInput(); + } catch (const RugnuxInputFailed &f) { + rc = f.code; + } + return rc; +} + // 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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d2075443c..9fbf9ff5f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -214,6 +214,7 @@ ADD_EXECUTABLE(jfjoch_test RugnuxResultFilesTest.cpp RugnuxInputKindTest.cpp RugnuxSweepsTest.cpp + RugnuxJointSweepsTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/RugnuxJointSweepsTest.cpp b/tests/RugnuxJointSweepsTest.cpp new file mode 100644 index 000000000..74ab41293 --- /dev/null +++ b/tests/RugnuxJointSweepsTest.cpp @@ -0,0 +1,108 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include + +#include "../rugnux/RugnuxJointSweeps.h" + +namespace { + +// A tetragonal P4 crystal: Laue class 4/m on a 4/mmm lattice, so (k, h, -l) indexes it a second way. +const gemmi::SpaceGroup &P4() { return *gemmi::find_spacegroup_by_name("P 4"); } +constexpr float A = 60.0f, C = 40.0f; + +std::array P4Key(int h, int k, int l) { + // the 4/m ASU representative, enough to give equivalents one intensity + std::array best{h, k, l}; + const std::array, 8> eq{{{h, k, l}, {-k, h, l}, {-h, -k, l}, {k, -h, l}, + {-h, -k, -l}, {k, -h, -l}, {h, k, -l}, {-k, h, -l}}}; + for (const auto &e: eq) best = std::max(best, e); + return best; +} + +JointMemberSweep MakeSweep(const std::map, double> &truth, const gemmi::Op &indexed_as, + int n_frames, uint32_t seed) { + std::mt19937 rng(seed); + JointMemberSweep s; + s.n_frames = n_frames; + s.cell = UnitCell{A, A, C, 90, 90, 90}; + s.space_group = P4(); + for (const auto &[hkl, I]: truth) { + const auto m = indexed_as.apply_to_hkl({{hkl[0], hkl[1], hkl[2]}}); + MergedReflection r; + r.h = m[0]; r.k = m[1]; r.l = m[2]; + r.I = static_cast(I); + r.sigma = 1.0f; + r.d = static_cast(1.0 / std::sqrt((hkl[0] * hkl[0] + hkl[1] * hkl[1]) / (A * A) + hkl[2] * hkl[2] / (C * C))); + s.merged.push_back(r); + ScaledFull f; + f.h = m[0]; f.k = m[1]; f.l = m[2]; + f.I = static_cast(I); + f.sigma = 1.0f; + f.image_number = static_cast(rng() % n_frames); + f.var_bkg = 4.0f; + f.var_per_I = 0.5f; + f.capture = 0.0f; + f.d = r.d; + s.fulls.push_back(f); + } + return s; +} + +std::map, double> Truth(uint32_t seed) { + std::mt19937 rng(seed); + std::exponential_distribution wilson(1.0); + std::map, double> t; + for (int h = -12; h <= 12; ++h) + for (int k = -12; k <= 12; ++k) + for (int l = -8; l <= 8; ++l) { + if (!(h || k || l)) continue; + const auto key = P4Key(h, k, l); + if (!t.contains(key)) t[key] = 1000.0 * wilson(rng); + } + return t; +} + +} // namespace + +TEST_CASE("RelateJointSweeps_AmbiguityRecoveredAndAStrangerLeftOut", "[Rugnux]") { + const auto truth = Truth(1); + const gemmi::Op twin = gemmi::parse_triplet("y,x,-z"); // h,k,l -> k,h,-l + std::vector sweeps; + sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2)); + sweeps.push_back(MakeSweep(truth, twin, 80, 3)); + sweeps.push_back(MakeSweep(Truth(99), gemmi::Op::identity(), 50, 4)); // another crystal + sweeps[0].fulls.push_back(sweeps[0].fulls.front()); // the reference: most fulls + + const auto relation = RelateJointSweeps(sweeps); + REQUIRE(relation.reference == 0); + CHECK(relation.sweeps[0].included); + REQUIRE(relation.sweeps[1].included); + CHECK(relation.sweeps[1].cc > 0.99); + CHECK(relation.sweeps[1].identity_cc < 0.5); + // the operator found takes sweep 2's indices back onto the reference's (up to the Laue group) + const auto back = relation.sweeps[1].cob.apply_to_hkl(twin.apply_to_hkl({{3, 1, 2}})); + CHECK(P4Key(back[0], back[1], back[2]) == P4Key(3, 1, 2)); + CHECK_FALSE(relation.sweeps[2].included); + + const auto input = BuildJointMergeInput(sweeps, relation); + REQUIRE(input.segments.size() == 2); + CHECK(input.segments[0].first_frame == 0); + CHECK(input.segments[0].n_frames == 100); + CHECK(input.segments[1].first_frame == 100); + CHECK(input.segments[1].n_frames == 80); + CHECK(input.segment_sweep == std::vector{0, 1}); + REQUIRE(input.outcomes.size() == 180); + // a full as the combine takes it: count I/c at correction c = var_per_I, variance var_bkg / c^2 + const Reflection *r = nullptr; + for (const auto &o: input.outcomes) + if (!o.reflections.empty()) { r = &o.reflections.front(); break; } + REQUIRE(r != nullptr); + CHECK(r->image_scale_corr == Catch::Approx(0.5)); + CHECK(r->var_bkg == Catch::Approx(16.0)); + CHECK(r->partiality == 1.0f); +}