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/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index ede0116af..8f41d08a4 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -384,6 +384,15 @@ J. Beilsten-Edmands, G. Winter, R. Gildea et al., "Scaling diffraction data in t package: algorithms and new approaches for multi-crystal scaling" (2020), Acta Cryst. D76, 385-399 [doi:10.1107/S2059798320003198](https://doi.org/10.1107/S2059798320003198). +**Several sweeps merged as one** - each sweep's partials combined into fulls on their own, and the +fulls of all the sweeps then scaled together, every sweep with its own frame scales and a relative B +against the joint merge, is how XSCALE merges several data sets (W. Kabsch, "Integration, scaling, +space-group assignment and post-refinement" (2010), above). Reading the radiation dose a sweep was taken +at from its relative B, which orders the sweeps by dose from the data rather than from the order they +were given in, follows K. Diederichs, "Some aspects of quantitative analysis and correction of radiation +damage" (2006), Acta Cryst. D62, 96-101 +[doi:10.1107/S0907444905031537](https://doi.org/10.1107/S0907444905031537). + **Uncertainty conventions** follow the IUCr Commission on Crystallographic Nomenclature: D. Schwarzenbach, S. C. Abrahams, H. D. Flack et al., "Statistical descriptors in crystallography: Report of the IUCr Subcommittee on Statistical Descriptors" (1989), Acta Cryst. A45, 63-75 diff --git a/docs/RUGNUX_CLI.md b/docs/RUGNUX_CLI.md index 5d5834b33..da955336f 100644 --- a/docs/RUGNUX_CLI.md +++ b/docs/RUGNUX_CLI.md @@ -30,9 +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 a_master.h5 b_master.h5 # several sweeps of one crystal, merged as one ``` -Chosen when the file's goniometer axis turns. One sweep per run. 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)). @@ -91,14 +92,34 @@ spot finding or integration ([Re-scaling and re-merging](RUGNUX_ADVANCED.md#re-s ## Inputs Inputs are files, never a directory, read in the order given (a shell glob gives them sorted by -name). Several inputs are read together only as one stills dataset in `--mode mx`: all the files of +name). Several inputs are one stills dataset in `--mode mx`, or several rotation sweeps (below): all the files of one dataset - `rugnux run*.h5` (SACLA), `*.cxi`, `*.xtc` (LCLS) - or several Jungfraujoch/NXmx masters of stills, read as one run: their images one after another, merged together, with the indexing rate reported per master as well. Such masters must share the detector, wavelength, -distance and beam centre. Several inputs must be one format; a rotation sweep (name one frame of a -CBF/marCCD/SMV sweep), any `--mode` other than `mx`, and following a broker take one. 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 +distance and beam centre. Several inputs must be one format; any `--mode` other than `mx`, and +following a broker, take one. + +Several rotation sweeps are told apart from the files: one sweep per CBF/marCCD/SMV file-name +template (the name with its frame number taken out - any frame of a sweep names the whole sweep, so +`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. +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 ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)). An atomic model and a reference reflection file can be given among the inputs, in any order, @@ -130,6 +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): `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 42b4e2b73..57a8962da 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) @@ -531,11 +532,28 @@ void RotationScaleMerge::IngestHost() { if (frame_limit >= 0) n_frames = std::min(n_frames, frame_limit); partials_released = false; resident_ingest = false; + frame_segment.clear(); + if (!segments.empty()) { + frame_segment.assign(n_frames, -1); + int32_t next = 0; + for (size_t k = 0; k < segments.size(); ++k) { + const SweepSegment &seg = segments[k]; + if (seg.first_frame != next || seg.n_frames <= 0 || seg.first_frame + seg.n_frames > n_frames) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "RotationScaleMerge: the sweep segments do not cover the frames in order"); + for (int32_t f = seg.first_frame; f < seg.first_frame + seg.n_frames; ++f) + frame_segment[f] = static_cast(k); + next = seg.first_frame + seg.n_frames; + } + if (next != n_frames) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "RotationScaleMerge: the sweep segments do not cover the frames in order"); + } #ifdef JFJOCH_USE_CUDA // Probed before the observations are built: whether the build lays down the full Obs array or // only the ingest arrays depends on whether the device pipeline will be active. gpu_ = std::make_unique(gpu_device); - gpu_active_ = gpu_->Available() && !scale_on_cpu_; + gpu_active_ = gpu_->Available() && !scale_on_cpu_ && segments.empty(); resident_ingest = gpu_active_ && !observation_dump; #endif @@ -605,7 +623,12 @@ void RotationScaleMerge::IngestHost() { } }); - MeasureIncidentFlux(mean_bkg); + // Fulls of several sweeps carry each sweep's own background: not one flux meter, so the per-frame + // scale alone takes the flux. + if (segments.empty()) + MeasureIncidentFlux(mean_bkg); + else + frame_flux.assign(n_frames, 1.0); // The occupied range of h, which is the buckets the sort works on. int b_lo = H_ALL, b_hi = -1; @@ -728,12 +751,19 @@ void RotationScaleMerge::IngestHost() { } }); } - BuildInRangeObservations(keys, static_cast(total), SmoothedFrameGeometry()); + // Pre-combined fulls have no rocking geometry left to smooth or to recompute a partiality from. + BuildInRangeObservations(keys, static_cast(total), + segments.empty() ? SmoothedFrameGeometry() : SmoothedGeometry{}); rawrun_group.assign(rawrun_start.size(), -1); logger->Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)", n_partials_obs, n_frames, rawrun_start.size()); - SmoothMosaicityAndPartiality(); + if (segments.empty()) { + SmoothMosaicityAndPartiality(); + } else { + ingest_zeta.reset(); + ingest_d.reset(); + } ParallelChunks(static_cast(partials.size()), ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) @@ -2174,6 +2204,7 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop( has_fulls[f] = frame_obs_count[f] > 0 ? 1 : 0; } double best_step = INFINITY, lambda = 0.0; + std::vector segment_lambda(segments.size(), 0.0); bool lambda_kept = false; std::vector steps, g_best; AndersonMixing anderson(SCALING_ANDERSON_DEPTH); @@ -2188,9 +2219,16 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop( if (fitted) y[f] += g_fit[f] > 0.0 ? std::max(std::log(g_fit[f] / g[f]), -1.0) : -1.0; } - if (!lambda_kept) - lambda = ChooseLogScaleSmoothing(y, J, std::max(1, rocking_event_frames_at_start), nthreads); - const std::vector x = SmoothLogScale(y, J, lambda); + std::vector x; + if (segments.empty()) { + if (!lambda_kept) + lambda = ChooseLogScaleSmoothing(y, J, std::max(1, rocking_event_frames_at_start), nthreads); + x = SmoothLogScale(y, J, lambda); + } else { + x = SmoothLogScaleBySegment(y, J, segment_lambda, !lambda_kept); + std::vector l = segment_lambda; + lambda = median_of(l); + } for (int f = 0; f < n_frames; ++f) g_new[f] = std::exp(x[f]); const double g_ref = TypicalFrameScale(g_new, n_obs); @@ -2342,8 +2380,14 @@ void RotationScaleMerge::RefineDecay(int n_groups) { return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && (o.d > 0.0f) && o.partiality >= min_partiality; }; + // Where an observation sits on the dose axis: its frame, or with several sweeps its place in the dose + // order the data put them in (MeasureSegmentDose). + auto dose = [&](const auto &o) { + return segments.empty() ? static_cast(o.image_number) + : frame_dose[o.frame] + (o.image_number - static_cast(o.frame)); + }; auto decay_factor = [&](const auto &o, double slope) { - return std::exp(slope * (o.image_number - fcenter) * s2_of(o.d)); + return std::exp(slope * (dose(o) - fcenter) * s2_of(o.d)); }; // The fields the seven passes below read (three slope fits, four held-out scores), pulled out of the @@ -2358,6 +2402,8 @@ void RotationScaleMerge::RefineDecay(int n_groups) { BuildDecayTerms(n_groups, t); decay_terms_kept = {}; decay_terms_kept_valid = false; + if (!segments.empty()) + MeasureSegmentDose(t, n_groups); const int nt = static_cast(ThreadsForWork(t.term.size(), nthreads)); std::vector sw(n_groups), swI(n_groups); @@ -2385,7 +2431,7 @@ void RotationScaleMerge::RefineDecay(int n_groups) { const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return; const double w = (Is / sc) * (Is / sc); - const double x = (o.image_number - fcenter) * s2_of(o.d), y = std::log(Iref / Is); + const double x = (dose(o) - fcenter) * s2_of(o.d), y = std::log(Iref / Is); a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y; }, [](Lin &a, const Lin &b) { for (int j = 0; j < 5; ++j) a[j] += b[j]; }); const double Sw = S[0], Sx = S[1], Sy = S[2], Sxx = S[3], Sxy = S[4]; @@ -2435,6 +2481,125 @@ void RotationScaleMerge::RefineDecay(int n_groups) { // sit at the same s^2 but different frames, so the correction spreads them) - net harmful. slope = // 2 dB/dframe, so the total relative-B change across the run is |slope/2 * n_frames|. constexpr double DECAY_MIN_DELTA_B = 2.0; // A^2, minimum total relative-B over the run to engage + if (!segments.empty()) { + // Two hypotheses of how damage accrued over several sweeps, decided on held-out frames: a B rising + // linearly along the dose order (frame_dose), or one B per sweep - a step at each sweep boundary, + // flat within a sweep. Both are regressions of y = ln(I_ref / I_obs) on s^2 = 1/(4 d^2), the ramp + // on x = (dose - centre) s^2, the steps on x = s^2 with a slope per sweep, at the weights the slope + // fit above uses and within each frame (below). Each is fitted on the even frames and its + // weighted squared residual taken on the odd ones, against their own reference, and the other way + // round; the hypothesis that predicts the held-out frames better is applied, where it beats no + // decay at all by the margin a correction has to clear on one sweep (and the physical floor). + // The held-out frames' residual is read at weights fixed by the observations, so a correction + // cannot win by reshaping sigma. + using Sums = std::array; // Sw, Sx, Sy, Sxx, Sxy + struct Line { double a = 0.0, b = 0.0; }; + auto line_of = [](const Sums &S) { + Line l; + const double var = S[0] > 0.0 ? S[3] - S[1] * S[1] / S[0] : 0.0; + if (var > 0.0) l.b = (S[4] - S[1] * S[2] / S[0]) / var; + if (S[0] > 0.0) l.a = (S[2] - l.b * S[1]) / S[0]; + return l; + }; + auto reference = [&](int parity) { + std::fill(sw.begin(), sw.end(), 0.0); + std::fill(swI.begin(), swI.end(), 0.0); + for (const DecayTerm &o : t.term) { + if (parity >= 0 && (o.frame & 1) != parity) continue; + const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + sw[o.group] += w; + swI[o.group] += w * static_cast(o.I) * o.corr; + } + }; + // Every term of `parity` with its y and weight against the current reference. + auto for_terms = [&](int parity, const auto &fn) { + for (const DecayTerm &o : t.term) { + if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue; + const double Iref = swI[o.group] / sw[o.group]; + const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + fn(o, std::log(Iref / Is), (Is / sc) * (Is / sc)); + } + }; + auto x_ramp = [&](const DecayTerm &o) { return (dose(o) - fcenter) * s2_of(o.d); }; + // The terms of `parity` with y and both x taken about their frame's weighted mean: each frame has a + // scale of its own (the per-frame G), which took up its part of the B averaged over its + // reflections, so only how y varies with resolution WITHIN a frame is decay. + struct Point { int32_t seg; double w, y, xr, xs; }; + auto points = [&](int parity) { + reference(parity); + std::vector pts; + std::vector frame; + std::vector> mean(n_frames, std::array{}); // Sw, Sy, Sxr, Sxs + for_terms(parity, [&](const DecayTerm &o, double y, double w) { + pts.push_back({frame_segment[o.frame], w, y, x_ramp(o), s2_of(o.d)}); + frame.push_back(o.frame); + auto &m = mean[o.frame]; + m[0] += w; m[1] += w * y; m[2] += w * pts.back().xr; m[3] += w * pts.back().xs; + }); + for (size_t k = 0; k < pts.size(); ++k) { + const auto &m = mean[frame[k]]; + pts[k].y -= m[1] / m[0]; pts[k].xr -= m[2] / m[0]; pts[k].xs -= m[3] / m[0]; + } + return pts; + }; + struct Fit { double ramp = 0.0; std::vector step; }; + auto fit = [&](int parity) { + Sums r{}; + std::vector st(segments.size(), Sums{}); + for (const Point &q : points(parity)) { + r[0] += q.w; r[1] += q.w * q.xr; r[2] += q.w * q.y; r[3] += q.w * q.xr * q.xr; r[4] += q.w * q.xr * q.y; + Sums &a = st[q.seg]; + a[0] += q.w; a[1] += q.w * q.xs; a[2] += q.w * q.y; a[3] += q.w * q.xs * q.xs; a[4] += q.w * q.xs * q.y; + } + Fit out{line_of(r).b, {}}; + for (const Sums &a : st) out.step.push_back(line_of(a).b); + return out; + }; + // Held-out weighted squared residuals of no decay, the ramp and the steps. + auto score = [&](int test, const Fit &m, std::array &ssr) { + for (const Point &q : points(test)) { + const double e1 = q.y - m.ramp * q.xr, e2 = q.y - m.step[q.seg] * q.xs; + ssr[0] += q.w * q.y * q.y; ssr[1] += q.w * e1 * e1; ssr[2] += q.w * e2 * e2; + } + }; + std::array ssr{}; + const Fit even = fit(0), odd = fit(1); + score(1, even, ssr); + score(0, odd, ssr); + const Fit all = fit(-1); + + const double slope = all.ramp; // the slope decay_factor applies + std::vector b(segments.size()); + double b_mid = 0.0, b_lo = INFINITY, b_hi = -INFINITY; + for (size_t k = 0; k < segments.size(); ++k) { + b[k] = 0.5 * all.step[k]; + b_mid += b[k] * segments[k].n_frames / static_cast(n_frames); + b_lo = std::min(b_lo, b[k]); + b_hi = std::max(b_hi, b[k]); + } + const double delta_lin = std::fabs(0.5 * slope * n_frames), delta_step = b_hi - b_lo; + const bool lin_ok = delta_lin >= DECAY_MIN_DELTA_B && ssr[1] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0]; + const bool step_ok = delta_step >= DECAY_MIN_DELTA_B && ssr[2] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0]; + const bool use_step = step_ok && (!lin_ok || ssr[2] < ssr[1]); + const bool use_lin = lin_ok && !use_step; + segment_decay_model = use_step ? "per-sweep B" : use_lin ? "dose-linear" : "none"; + const double rel = ssr[0] > 0.0 ? 100.0 / ssr[0] : 0.0; + logger->Info("Decay over {} sweeps, held-out residual against no decay: dose-linear (dB = {:.2f} A^2 over " + "the dose order) {:.1f}%, one B per sweep (spread {:.2f} A^2) {:.1f}% => {}", segments.size(), + delta_lin, ssr[1] * rel, delta_step, ssr[2] * rel, segment_decay_model); + if (use_step || use_lin) + ParallelChunks(static_cast(fulls.size()), ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) { + Obs &o = fulls[i]; + if (!usable(o)) continue; + const double f = use_step ? std::exp(2.0 * (b[frame_segment[o.frame]] - b_mid) * s2_of(o.d)) + : decay_factor(o, slope); + o.corr = static_cast(o.corr * f); + } + }); + return; + } const double slope = fit_slope(-1); const double total_delta_B = std::fabs(0.5 * slope * n_frames); if (total_delta_B < DECAY_MIN_DELTA_B) { @@ -2460,6 +2625,92 @@ void RotationScaleMerge::RefineDecay(int n_groups) { 0.5 * slope * n_frames, 0.5 * slope); } +std::vector RotationScaleMerge::SegmentB(const DecayTerms &t, int n_groups, int parity) const { + // Per segment, the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against + // s^2 = 1/(4 d^2), which is 2 B, with an intercept of its own (the scale is the per-frame G's), over + // the frames of `parity` (< 0: all) against a reference made of the same frames. NaN where a segment + // has nothing to fit. + const size_t ns = segments.size(); + std::vector sw(n_groups, 0.0), swI(n_groups, 0.0); + for (const DecayTerm &o : t.term) { + if (parity >= 0 && (o.frame & 1) != parity) continue; + const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + sw[o.group] += w; + swI[o.group] += w * static_cast(o.I) * o.corr; + } + std::vector> S(ns, std::array{}); // Sw, Sx, Sy, Sxx, Sxy + for (const DecayTerm &o : t.term) { + if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue; + const double Iref = swI[o.group] / sw[o.group]; + const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; + if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; + const double w = (Is / sc) * (Is / sc); + const double x = 1.0 / (4.0 * static_cast(o.d) * o.d), y = std::log(Iref / Is); + auto &a = S[frame_segment[o.frame]]; + a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y; + } + std::vector b(ns, NAN); + for (size_t k = 0; k < ns; ++k) { + const auto &a = S[k]; + const double var = a[0] > 0.0 ? a[3] - a[1] * a[1] / a[0] : 0.0; + if (var > 0.0) + b[k] = 0.5 * (a[4] - a[1] * a[2] / a[0]) / var; + } + return b; +} + +void RotationScaleMerge::MeasureSegmentDose(const DecayTerms &t, int n_groups) { + // A relative B per data set against the joint merge, as XSCALE fits it: Kabsch (2010) Acta Cryst. + // D66, 133-144; the dose read from it: Diederichs (2006) Acta Cryst. D62, 96-101. + // Radiation damage raises the relative B of a sweep collected later, so the B factors of the + // segments against their joint merge give the order the crystal took the dose in - read from the + // data, not from the order the sweeps were given in or the times their headers state. Per segment, + // the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against s^2 = 1/(4 d^2) is 2 B, + // with an intercept of its own (the scale is the per-frame G's). The least damaged segment is B = 0, + // and its frames start the dose axis. + const size_t ns = segments.size(); + const std::vector b = SegmentB(t, n_groups, -1); + // A segment with nothing to measure keeps its place in the input order, after the measured ones. + std::vector order(ns); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), [&](int p, int q) { + if (std::isfinite(b[p]) != std::isfinite(b[q])) return std::isfinite(b[p]); + return std::isfinite(b[p]) && b[p] < b[q]; + }); + const double b0 = std::isfinite(b[order[0]]) ? b[order[0]] : 0.0; + frame_dose.assign(n_frames, 0.0); + double offset = 0.0; + std::string text; + for (size_t rank = 0; rank < ns; ++rank) { + const int k = order[rank]; + const SweepSegment &seg = segments[k]; + for (int32_t f = seg.first_frame; f < seg.first_frame + seg.n_frames; ++f) + frame_dose[f] = offset + (f - seg.first_frame); + offset += seg.n_frames; + segment_results[k].relative_b = b[k] - b0; + segment_results[k].dose_rank = static_cast(rank); + text += fmt::format("{}segment {} (B {:+.2f} A^2)", rank == 0 ? "" : ", ", k + 1, b[k] - b0); + } + logger->Info("Dose order of the {} sweep segments, from their relative B against the joint merge: {}", + ns, text); +} + +std::vector RotationScaleMerge::SmoothLogScaleBySegment(const std::vector &y, + const std::vector &J, + std::vector &lambda, bool choose) const { + std::vector x(y.size()); + for (size_t k = 0; k < segments.size(); ++k) { + const auto lo = y.begin() + segments[k].first_frame, hi = lo + segments[k].n_frames; + const auto jlo = J.begin() + segments[k].first_frame; + const std::vector ys(lo, hi), Js(jlo, jlo + segments[k].n_frames); + if (choose) + lambda[k] = ChooseLogScaleSmoothing(ys, Js, std::max(1, rocking_event_frames_at_start), nthreads); + const std::vector xs = SmoothLogScale(ys, Js, lambda[k]); + std::copy(xs.begin(), xs.end(), x.begin() + segments[k].first_frame); + } + return x; +} + // A^2, the bound the Gauss-Seidel solve clamps each batch's relative-B to. It is a numerical guard against // divergence, never an answer: a run of identical values in the solved curve means the solve reached the // bound, not that the crystal did anything. The CORRECTION keeps the tight bound it was tuned with; the @@ -4716,12 +4967,15 @@ void RotationScaleMerge::Combine() { size_t i = 0; while (i < ev.size()) { size_t kk = i + 1; - float last_frame = partials[ev[i]].r->image_number; - while (kk < ev.size()) { - const float frame = partials[ev[kk]].r->image_number; - if (frame - last_frame > max_frame_gap) break; - last_frame = frame; - ++kk; + // A pre-combined full is a whole event already (SweepSegment). + if (segments.empty()) { + float last_frame = partials[ev[i]].r->image_number; + while (kk < ev.size()) { + const float frame = partials[ev[kk]].r->image_number; + if (frame - last_frame > max_frame_gap) break; + last_frame = frame; + ++kk; + } } double pooled_bkg = 0.0; @@ -5860,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}; } }); } @@ -6697,6 +6964,12 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st const gemmi::SpaceGroup *group) { const ScopedGpuPin pin(gpu_device); // see IngestHost run_group = group; + // Pre-combined fulls: there are no partials to fit a scale on, and every event is one frame. + if (!segments.empty()) { + scale_partials_choice = false; + rocking_event_frames_at_start = 1; + rocking_events_at_start = n_partials_obs; + } struct ClearGroup { const gemmi::SpaceGroup *&g; ~ClearGroup() { g = nullptr; } } clear_group{run_group}; // Whether a frame's scale can be fitted on its partials depends on whether the partials tell a change // of scale from an error of the partiality model. The count of rocking events per frame is the prior @@ -7047,7 +7320,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // kernels on it, and the restore re-uploads. // --- 2a. On the de-novo search pass only, drop observations whose Lorentz geometry is poor. --- - if (for_search && search_min_zeta > 0.0) { + if (for_search && search_min_zeta > 0.0 && segments.empty()) { // a full has no rocking geometry int64_t n_dropped = 0; bool zeta_on_gpu = false; #ifdef JFJOCH_USE_CUDA @@ -7357,20 +7630,27 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // is a correction to the merged INTENSITIES, and a pass whose intensities are thrown away has // nowhere to put one - the fitted gain lives on the fulls, which the next pass rebuilds. --- const bool refine_surfaces = !for_search && full_stats; + // What runs along the frames of ONE sweep - the damage monitor, the sweep ledger, the per-batch B and + // the surfaces, each read in one goniometer and one detector setting - has no meaning across the + // segments of several, and is left out there. The decay is fitted along the dose order instead. + const bool along_frames = segments.empty(); const bool corrections = - refine_surfaces && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0); + refine_surfaces && (refine_decay_b || (along_frames && (absorption_iter > 0 || modulation_iter > 0 + || relative_b_deg > 0.0))); + segment_results.assign(segments.size(), Result::Segment{}); + segment_decay_model.clear(); // Radiation-damage monitor: measure the per-batch relative-B on the scaled fulls BEFORE any correction // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass. std::unique_ptr sub_mark; auto sub = [&](const char *n) { sub_mark.reset(); sub_mark = std::make_unique(*logger, std::string("rsm/s ") + n); }; sub("damage B"); - if (refine_surfaces) + if (refine_surfaces && along_frames) MeasureRadiationDamageB(n_groups); // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is // already out of it) and the per-frame CC computed above. It drops nothing itself; the ranges it // segments are the ledger 4c below decides on. sub("sweep quality"); - if (refine_surfaces) + if (refine_surfaces && along_frames) MeasureSweepQuality(partial_scaled, cc, cc_n); // CC1/2 of these same fulls with the surfaces not yet folded in, for the caller's cross-pass quality // guard: the pass this one is judged against fits no surfaces of its own, so it is the uncorrected @@ -7396,7 +7676,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool decay_terms_kept = {}; decay_terms_kept_valid = false; sub("relative B"); - if (refine_surfaces && relative_b_deg > 0.0) + if (refine_surfaces && along_frames && relative_b_deg > 0.0) RefineRelativeB(n_groups); // per-batch relative-B on top of the single decay slope // The three surfaces overlap, so the one fitted first gets the first claim and the others are tested // on what it left - the order decides what is adopted. Detector first: every frame of the sweep @@ -7412,16 +7692,16 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // the grid is then tested on what the smooth surface left, and where that surface is refused the // grid sees exactly what it saw before. sub("modulation"); - if (refine_surfaces && modulation_iter > 0) + if (refine_surfaces && along_frames && modulation_iter > 0) RefineModulation(modulation_iter, n_groups); sub("absorption time"); - if (refine_surfaces && absorption_iter > 0) + if (refine_surfaces && along_frames && absorption_iter > 0) RefineAbsorptionTime(absorption_iter, n_groups); sub("absorption SH"); - if (refine_surfaces && absorption_iter > 0) + if (refine_surfaces && along_frames && absorption_iter > 0) RefineAbsorptionSH(absorption_iter, n_groups); sub("absorption gonio"); - if (refine_surfaces && absorption_iter > 0) + if (refine_surfaces && along_frames && absorption_iter > 0) RefineAbsorption(absorption_iter, n_groups); sub("disposition"); @@ -7462,6 +7742,17 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // its CC1/2 is the answer without a second merge. A caller that did not ask keeps the NaN. if (measure_cc_before_corrections) r.cc_half_before_corrections = cc_half_before_corrections.value_or(r.statistics.overall.cc_half); + // Each segment's scale: the median of its frames' scales, at the run's gauge. + for (size_t k = 0; k < segments.size(); ++k) { + std::vector g; + for (int32_t f = segments[k].first_frame; f < segments[k].first_frame + segments[k].n_frames; ++f) + if (fulls_frame_count[f] > 0 && std::isfinite(g_full[f]) && g_full[f] > 0.0) + g.push_back(g_full[f]); + if (!g.empty()) + segment_results[k].scale = median_of(g); + } + r.segments = segment_results; + r.segment_decay_model = segment_decay_model; // The fulls and the merge's working arrays are rebuilt by every arm before they are read, so they go // with it: kept, they sat beside the next arm's combine and through everything the caller did between // merges - on a large cell in P1, the fulls and about 330 bytes per ASU group, gigabytes. Faulting them diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 6291bdf1e..d1ebd935a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -139,6 +140,16 @@ public: // The fulls the merge kept, unmerged (see ScaledFull): what passed the outlier tests, inside the // resolution cut. Filled only when SetExportScaledFulls asked for it, and never on a search merge. std::vector scaled_fulls; + // Per sweep segment (SetSweepSegments), in segment order; empty without segments. + struct Segment { + double scale = NAN; // median per-frame scale of the segment's fulls (run gauge) + double relative_b = NAN; // A^2, against the joint merge, the least damaged segment at 0 + int dose_rank = -1; // the segment's place in the dose order those B factors decide + }; + std::vector segments; + // How the decay over the segments was modelled - "dose-linear" (along the dose order), "per-sweep B" + // (one B a segment) or "none" - as the held-out frames decided; empty without segments or decay. + std::string segment_decay_model; }; // experiment: read live (its space group is changed by the caller between Run() calls). @@ -288,6 +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; 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; } @@ -724,8 +745,18 @@ private: std::vector last_cc_n; std::vector last_frame_scaled; int frame_limit = -1; // see SetFrameLimit + std::vector segments; // see SetSweepSegments + std::vector frame_segment; // each frame's segment (Ingest; empty without segments) + std::vector segment_results; // the last arm's, copied into its Result + // Each frame's place on the dose axis the decay is fitted along: the frames of the segments laid end + // to end in the order the segments' relative B factors put them (MeasureSegmentDose). + std::vector frame_dose; + // SmoothLogScale / ChooseLogScaleSmoothing within each segment on its own; lambda per segment. + std::vector SmoothLogScaleBySegment(const std::vector &y, const std::vector &J, + 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) --- @@ -958,6 +989,12 @@ private: std::vector g_start; }; void BuildDecayTerms(int n_groups, DecayTerms &t) const; + // Per segment: the relative B against the merged reference, and the dose order it puts the segments + // in (segment_results, frame_dose). + void MeasureSegmentDose(const DecayTerms &t, int n_groups); + // Each segment's B against a reference made of the same frames: those of `parity` (< 0: all). + std::vector SegmentB(const DecayTerms &t, int n_groups, int parity) const; + std::string segment_decay_model; // the last decay's choice (Result::segment_decay_model) // The terms the radiation-damage monitor built, handed on to RefineDecay: nothing between the two // touches corr, so they are the terms RefineDecay would build. Empty outside that stretch of RunArm. DecayTerms decay_terms_kept; diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 396d610f5..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 @@ -40,7 +42,8 @@ TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochPrevie # The files of a run, written by the programs that run one (rugnux, the viewer) from what JFJochRugnux # returns: the library itself writes none. -ADD_LIBRARY(JFJochRugnuxOutput STATIC RugnuxOutputFiles.cpp RugnuxOutputFiles.h RugnuxInputKind.cpp RugnuxInputKind.h) +ADD_LIBRARY(JFJochRugnuxOutput STATIC RugnuxOutputFiles.cpp RugnuxOutputFiles.h RugnuxInputKind.cpp RugnuxInputKind.h + RugnuxSweeps.cpp RugnuxSweeps.h) TARGET_LINK_LIBRARIES(JFJochRugnuxOutput JFJochRugnux JFJochWriter JFJochModelWriter ZLIB::ZLIB) # rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 90fa50985..b1dc52b18 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -823,6 +823,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; @@ -854,7 +869,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; @@ -943,7 +960,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); @@ -1156,6 +1173,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(); @@ -1192,7 +1213,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. @@ -1230,7 +1263,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); @@ -1241,7 +1274,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 d5defa3b7..21d77bf6d 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -115,6 +115,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). @@ -1289,6 +1294,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 @@ -1691,6 +1701,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 4d031458d..8153d745f 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1693,6 +1693,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/RugnuxSweeps.cpp b/rugnux/RugnuxSweeps.cpp new file mode 100644 index 000000000..1f88a5de4 --- /dev/null +++ b/rugnux/RugnuxSweeps.cpp @@ -0,0 +1,140 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RugnuxSweeps.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace { + // How far two sweeps may state the wavelength apart and still be one dataset: JFJochMultiReader's + // tolerance for the inputs of one stills run. + constexpr double WAVELENGTH_TOLERANCE = 1e-3; // relative + + std::string Lower(std::string s) { + std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); }); + return s; + } +} + +std::string RugnuxFrameTemplate(const std::string &path) { + const std::filesystem::path p(path); + std::string name = p.filename().string(); + size_t end = name.size(); + while (end > 0 && !std::isdigit(static_cast(name[end - 1]))) + end--; + if (end == 0) + return {}; + size_t start = end; + while (start > 0 && std::isdigit(static_cast(name[start - 1]))) + start--; + name.replace(start, end - start, std::string(end - start, '#')); + return (p.parent_path() / name).string(); +} + +RugnuxSweepPlan PlanRugnuxSweeps(const std::vector &inputs, + const std::function &describe_file, + const std::function &describe_sweep) { + RugnuxSweepPlan plan; + if (inputs.size() < 2) + return plan; + + std::vector files; + for (const auto &f: inputs) + files.push_back(describe_file(f)); + const bool any_foreign = std::any_of(files.begin(), files.end(), [](const auto &d) { return d.foreign; }); + const size_t n_frames = std::count_if(files.begin(), files.end(), + [](const auto &d) { return !d.frame_template.empty(); }); + // Inputs of the stills readers, or a mix of frames and masters: read as before, which takes the + // former and refuses the latter. + if (any_foreign || (n_frames != 0 && n_frames != inputs.size())) + return plan; + + if (n_frames == inputs.size()) { + std::set seen; + for (size_t i = 0; i < inputs.size(); i++) + if (seen.insert(files[i].frame_template).second) + plan.sweeps.push_back({inputs[i], files[i].frame_template, true}); + if (plan.sweeps.size() == 1) + return plan; + } else { + // A first input of stills is a stills run, read together as before - and refused there if any + // other input turns. A rotation first input asks every other one. + if (!describe_sweep(inputs[0]).rotation) + return plan; + std::vector rotation{inputs[0]}, stills; + for (size_t i = 1; i < inputs.size(); i++) + (describe_sweep(inputs[i]).rotation ? rotation : stills).push_back(inputs[i]); + if (!stills.empty()) + throw std::runtime_error(fmt::format( + "{} is a rotation sweep and {} a stills collection: rotation sweeps and stills are not " + "processed together; give them in separate runs", rotation.front(), stills.front())); + for (const auto &f: inputs) + plan.sweeps.push_back({f, f, false}); + } + + std::vector wavelength; + for (const auto &s: plan.sweeps) + wavelength.push_back(describe_sweep(s.input).wavelength_A); + for (size_t k = 1; k < plan.sweeps.size(); k++) + if (std::fabs(wavelength[k] - wavelength[0]) > WAVELENGTH_TOLERANCE * wavelength[0]) + throw std::runtime_error(fmt::format( + "{} is at {:.5f} A and {} at {:.5f} A: sweeps at different wavelengths are not processed " + "together; give them in separate runs", plan.sweeps[0].label, wavelength[0], + plan.sweeps[k].label, wavelength[k])); + plan.wavelength_A = wavelength[0]; + return plan; +} + +std::vector RugnuxSiblingSweeps(const RugnuxSweepPlan &plan, const std::vector &extensions) { + namespace fs = std::filesystem; + std::set given, dirs_frames, dirs_masters; + for (const auto &s: plan.sweeps) { + const fs::path p(s.input); + const std::string dir = p.parent_path().empty() ? "." : p.parent_path().string(); + if (s.frames) { + given.insert(fs::path(s.label).filename().string()); + dirs_frames.insert(dir); + } else { + given.insert(p.filename().string()); + dirs_masters.insert(dir); + } + } + auto has_extension = [&](std::string name) { + name = Lower(name); + if (name.size() > 3 && name.ends_with(".gz")) + name.resize(name.size() - 3); + return std::any_of(extensions.begin(), extensions.end(), [&](const auto &e) { return name.ends_with(e); }); + }; + + // template (or master) name -> files, in name order + std::map found; + std::error_code ec; + for (const auto &dir: dirs_frames) + for (const auto &e: fs::directory_iterator(dir, ec)) { + const std::string name = e.path().filename().string(); + if (!e.is_regular_file() || !has_extension(name)) + continue; + const std::string t = fs::path(RugnuxFrameTemplate(e.path().string())).filename().string(); + if (!t.empty() && !given.contains(t)) + found[(fs::path(dir) / t).string()]++; + } + for (const auto &dir: dirs_masters) + for (const auto &e: fs::directory_iterator(dir, ec)) { + const std::string name = e.path().filename().string(); + if (e.is_regular_file() && name.ends_with("_master.h5") && !given.contains(name)) + found[(fs::path(dir) / name).string()] = 0; + } + + std::vector out; + for (const auto &[name, n]: found) + out.push_back(n > 0 ? fmt::format("{} ({} frame{})", name, n, n == 1 ? "" : "s") : name); + return out; +} diff --git a/rugnux/RugnuxSweeps.h b/rugnux/RugnuxSweeps.h new file mode 100644 index 000000000..eec15b1a7 --- /dev/null +++ b/rugnux/RugnuxSweeps.h @@ -0,0 +1,56 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +// The rotation sweeps a rugnux command line names, decided from the files in the order given: +// - frame files (CBF, marCCD, SMV): one sweep per file-name template - the name with its last run of +// digits, the frame number, taken out. Any frame of a sweep names the whole sweep, so several frames +// of one template are one sweep; +// - masters: one sweep per rotation master. Several stills masters stay one stills run, as before, and +// a mix of rotation and stills masters is refused. +// Sweeps at different wavelengths are refused. Nothing is planned for a single input, or for inputs only +// the stills readers take (XFEL HDF5, CXI, XTC): those are read as they always were. + +// What the planner needs to know about one input. `frame_template` is empty for anything but a frame +// file; `rotation` and `wavelength_A` are asked only of the input that stands for a sweep. +struct RugnuxInputDescription { + std::string frame_template; + bool foreign = false; // read only by the stills readers (XFEL HDF5, CXI, XTC) + bool rotation = false; + double wavelength_A = 0.0; +}; + +struct RugnuxSweep { + std::string input; // the file the sweep is read from: its first frame given, or its master + std::string label; // the frame template or the master, for messages + bool frames = false; // a sweep of frame files (label is its template) +}; + +struct RugnuxSweepPlan { + // Empty: the inputs are read together as before (one input, a stills run, or a mix the readers refuse). + std::vector sweeps; + double wavelength_A = 0.0; // of every sweep, where there are several +}; + +// `describe_file` says what an input is from its name and contents without reading the sweep (its +// frame template, or whether it is only for the stills readers); `describe_sweep` reads the input that +// stands for a sweep for its rotation and wavelength. Throws std::runtime_error, with a message for the +// user, on a mix of rotation sweeps and stills, or on sweeps at different wavelengths. +RugnuxSweepPlan PlanRugnuxSweeps(const std::vector &inputs, + const std::function &describe_file, + const std::function &describe_sweep); + +// The template of a frame file name, with its directory: the last run of digits in the name replaced by +// '#' characters, one per digit ("/d/x_1_0001.cbf" -> "/d/x_1_####.cbf"). Empty for a name with no digits. +std::string RugnuxFrameTemplate(const std::string &path); + +// The sweeps beside the given ones that were not given: frame templates of `extensions` (lower case, with +// the dot; a trailing ".gz" is looked through) in the directories of the frame sweeps, and other +// "*_master.h5" files beside the given masters - one line each, for the log. Only the inputs given are +// processed; these are named so that a sweep left out by accident shows. +std::vector RugnuxSiblingSweeps(const RugnuxSweepPlan &plan, const std::vector &extensions); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 0e2156f43..d2c4424b3 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -86,6 +86,8 @@ #include "RugnuxOutputFiles.h" #include "RugnuxConsole.h" #include "RugnuxInputKind.h" +#include "RugnuxSweeps.h" +#include "RugnuxJointSweeps.h" #ifdef __APPLE__ #include "../image_analysis/bragg_integration/BraggIntegrationEngineMetal.h" #endif @@ -101,6 +103,9 @@ constexpr double RELATIVE_B_DEFAULT_DEG = 10.0; // mutually exclusive by construction. enum class RugnuxMode { MX, Azint, Scale, Calibration, Refine }; +// --sweeps: several rotation sweeps as one dataset (joint), or each alone into _s (separate). +enum class RugnuxSweepsMode { Joint, Separate }; + // 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 @@ -120,6 +125,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 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. @@ -167,6 +173,8 @@ Basic options --no-gpu CPU only, even where a GPU is present --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_*, 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: @@ -280,7 +288,8 @@ enum { OPT_POLARIZATION, OPT_NO_REFINE_TILT, OPT_REFINE_SCORE_ONLY, - OPT_REFINE_D_MIN + OPT_REFINE_D_MIN, + OPT_SWEEPS }; static option long_options[] = { @@ -304,6 +313,7 @@ static option long_options[] = { {"no-refine-solvent", no_argument, nullptr, OPT_NO_REFINE_SOLVENT}, {"refine-score-only", no_argument, nullptr, OPT_REFINE_SCORE_ONLY}, {"refine-d-min", required_argument, nullptr, OPT_REFINE_D_MIN}, + {"sweeps", required_argument, nullptr, OPT_SWEEPS}, {"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}, @@ -805,70 +815,27 @@ extern "C" void fftwf_assertion_failed(const char *s, int line, const char *file std::abort(); } -// 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, bool no_gpu) { - // cuFFT plans are reused through the run and destroyed with it, ahead of the device reset in main. - FFTPlanReuse fft_plans; - // 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(); -#ifdef __APPLE__ - // Without CUDA, the Bragg integration's integer pixel work runs on the Apple GPU (MXAnalysisWithoutFPGA). - const std::string metal_device = BraggIntegrationEngineMetal::DeviceName(); -#else - const std::string metal_device; -#endif - provenance.gpu_description = get_gpu_description(); +namespace { + // One input's run ended: what exit(EXIT_FAILURE) was where a command line ran one input only. A run of + // several sweeps goes on with the next sweep; RunRugnux returns the code. + struct RugnuxInputFailed { int code; }; + // What Parse and PrepareRugnuxInput return when the run goes on. + constexpr int RUGNUX_CONTINUE = -1; +} - // The terminal shows the run as steps (RugnuxConsole) and every Logger writes to the log file - // instead; with -v the whole log is on the terminal as well. -v is looked for here, before the - // options are parsed, because the parser already logs. +// The command line as parsed, with what RunRugnux set up before it was parsed. The options are members so +// that the code reading and running one input (RugnuxInputRun) reads them by name, as it did when it was +// one function, and so that each sweep of a run of several is run from a copy of its own. +struct RugnuxCli { + bool no_gpu; bool verbose = false; - for (int i = 1; i < argc; i++) - if (std::strcmp(argv[i], "-v") == 0 || std::strcmp(argv[i], "--verbose") == 0) - verbose = true; - console = std::make_shared(verbose); - std::atexit(CloseConsoleOnExit); - // Files with an older time than this are not this run's (the margin covers coarse file times). - const auto files_since = std::filesystem::file_time_type::clock::now() - std::chrono::seconds(2); - auto log_sinks = std::make_shared(); - // What is logged before the log file can be opened - its name comes from -o - waits here. - auto early_log = std::make_shared(10000); - if (verbose) - log_sinks->add_sink(std::make_shared()); - else - log_sinks->add_sink(MakeConsoleSink(console)); - log_sinks->add_sink(early_log); - Logger::RouteAll(log_sinks); - // The HDF5 library prints its own error stack on stderr; every failure it reports reaches the - // reader as an error code anyway, which is turned into a message there. - H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr); - - PrintBanner(*console); - - Logger logger("rugnux"); -#ifdef JFJOCH_USE_CUDA - const char *build = "CUDA build"; -#else - const char *build = "CPU build"; -#endif - logger.Info("rugnux Version {}{} - {}", jfjoch_version(), - jfjoch_git_sha1().empty() ? std::string() - : fmt::format(" (git {} {})", jfjoch_git_sha1().substr(0, 6), jfjoch_git_date()), build); - logger.Info("Build flags: {}", jfjoch_build_cxx_flags().empty() ? "none" : jfjoch_build_cxx_flags()); - logger.Info("Build: {}", jfjoch_build_info(true)); - logger.Info("Command line: {}", provenance.command_line); - logger.Info("GPU: {}", no_gpu ? std::string("not used (--no-gpu) - running on the CPU") - : provenance.gpu_count > 0 ? fmt::format("{} ({} visible)", provenance.gpu_description, provenance.gpu_count) - : !metal_device.empty() ? fmt::format("{} through Metal, for the Bragg integration only - the rest runs on the CPU", metal_device) - : std::string("none visible - running on the CPU")); - RegisterHDF5Filter(); + std::chrono::steady_clock::time_point invocation_start; + RunProvenance provenance; + std::string metal_device; + std::filesystem::file_time_type files_since; + std::shared_ptr log_sinks; + std::shared_ptr early_log; + Logger logger{"rugnux"}; // made after RunRugnux has routed every Logger to the log std::string input_file; std::string output_prefix = "output"; @@ -991,6 +958,79 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement std::optional outlier_reject_nsigma; // merge per-observation outlier rejection + std::vector input_files; + 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_, + std::shared_ptr log_sinks_, + std::shared_ptr early_log_) + : no_gpu(no_gpu_), invocation_start(invocation_start_), provenance(std::move(provenance_)), + metal_device(std::move(metal_device_)), files_since(files_since_), log_sinks(std::move(log_sinks_)), + early_log(std::move(early_log_)) {} + + // 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(); +}; + +// One input - one sweep, or one stills run - read, set up and processed. Copied from the command line it +// belongs to, with its own input and output prefix. +struct RugnuxInputRun : RugnuxCli { + explicit RugnuxInputRun(const RugnuxCli &cli) : RugnuxCli(cli) {} + + // What PrepareRugnuxInput sets up and RunPreparedInput runs with. + std::future cufft_warmup, context_warmup; + SimulatedLive hdf5_reader; + SimulatedLive cbf_reader; + SimulatedLive marccd_reader; + SimulatedLive smv_reader; + SimulatedLive xfel_reader; + SimulatedLive xtc_reader; + SimulatedLive sfrm_reader; + JFJochStartMessageReader live_reader; + std::unique_ptr multi_reader; + std::unique_ptr broker_feed; + JFJochReader *reader_ptr = nullptr; + std::shared_ptr dataset; + OutputProvenance output_provenance; + DiffractionExperiment experiment; + ProcessConfig config; + 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). + int PrepareRugnuxInput(); + // Runs the input PrepareRugnuxInput set up and writes its files. + int RunPreparedInput(); + int Process() { + const int rc = PrepareRugnuxInput(); + return rc != RUGNUX_CONTINUE ? rc : RunPreparedInput(); + } +}; + +int RugnuxCli::Parse(int argc, char **argv) { if (argc == 1) { print_usage(); exit(EXIT_FAILURE); @@ -1198,6 +1238,16 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { case OPT_REFINE_D_MIN: refine_d_min = parse_number_arg(optarg, "--refine-d-min", logger, 0.1, 100.0); break; + case OPT_SWEEPS: + if (std::strcmp(optarg, "joint") == 0) + sweeps_mode = RugnuxSweepsMode::Joint; + else if (std::strcmp(optarg, "separate") == 0) + sweeps_mode = RugnuxSweepsMode::Separate; + else { + logger.Error("--sweeps expects joint or separate, not {}", optarg); + exit(EXIT_FAILURE); + } + break; case OPT_EXPORT_UNMERGED: export_unmerged = true; break; @@ -1625,8 +1675,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // The inputs, then a broker URL if the last word is one: told apart by the scheme, so that several // stills inputs can be given. - std::vector input_files(argv + optind, argv + argc); - std::string broker_url; + input_files.assign(argv + optind, argv + argc); + broker_url.clear(); if (!input_files.empty() && (input_files.back().starts_with("http://") || input_files.back().starts_with("https://"))) { broker_url = input_files.back(); @@ -1696,7 +1746,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { exit(EXIT_FAILURE); } // Read before anything runs: a file that holds no model ends the run here rather than after it. - std::shared_ptr model; + model.reset(); if (const std::string path = mr_model.empty() ? model_pdb : mr_model; !path.empty()) { try { model = std::make_shared(ReadRugnuxModelFile(path)); @@ -1712,6 +1762,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { exit(EXIT_FAILURE); } + return RUGNUX_CONTINUE; +} + +int RugnuxInputRun::PrepareRugnuxInput() { input_file = input_files[0]; // -N defaults to 0 = "use every CPU the run may use" (AvailableCpus); resolve it to a concrete count here so every mode @@ -1729,13 +1783,14 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // 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); for (const auto &m : early_log->last_raw()) file->log(m); log_sinks->add_sink(file); + log_file_sink = file; } catch (const std::exception &e) { logger.Warning("The log file {} cannot be written: {}", log_file, e.what()); log_file.clear(); @@ -1768,7 +1823,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { data = ReadRugnuxReflectionFile(input_file); } catch (const std::exception &e) { logger.Error("Cannot read the reflections {}: {}", input_file, e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } console->OnStep(refine_score_only ? "Scoring the model" : "Refining the model"); RefinementSettings settings; @@ -1815,27 +1870,16 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // 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. It also times a plan, which // shows whether the CUDA JIT cache works. - auto cufft_warmup = std::async(std::launch::async, [&logger] { LogCudaJitCache(logger, WarmUpCuFFT()); }); + cufft_warmup = std::async(std::launch::async, [this] { LogCudaJitCache(logger, WarmUpCuFFT()); }); // The same for the CUDA contexts of the other cards (the warm-up above creates card 0's), which // the pre-scan does not use: created one after another at start-up, they cost a fifth of a second // per card before the run began. - auto context_warmup = std::async(std::launch::async, create_gpu_contexts, 1); + context_warmup = std::async(std::launch::async, create_gpu_contexts, 1); // 1. Read Input File. An HDF5 master, one PILATUS miniCBF frame, or one marCCD frame - naming // any frame of a sweep reads the whole sweep, which is how these are archived and how a user // thinks of them. console->OnStep("Reading input"); - SimulatedLive hdf5_reader; - SimulatedLive cbf_reader; - SimulatedLive marccd_reader; - SimulatedLive smv_reader; - SimulatedLive xfel_reader; - SimulatedLive xtc_reader; - SimulatedLive sfrm_reader; - JFJochStartMessageReader live_reader; - std::unique_ptr multi_reader; - std::unique_ptr broker_feed; - JFJochReader *reader_ptr = nullptr; // Which format this is. Asked INSIDE the try: these open and read the file, and an input the program // cannot read has to come out as the message below, not as an unhandled exception with nothing // printed at all. @@ -1843,13 +1887,13 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { bool input_is_cxi = false, input_is_xtc = false, input_is_panel_table = false, input_is_sfrm = false; try { // Every input is the format the first one is. - const auto all_inputs = [&input_files](auto can_read) { + const auto all_inputs = [this](auto can_read) { return std::all_of(input_files.begin(), input_files.end(), can_read); }; const auto refuse_mixed = [&]() { logger.Error("The inputs are not all of one format; several inputs must be files of the same " "kind - Jungfraujoch/NXmx masters, SACLA HDF5, CXI or XTC"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; }; const auto is_foreign = [](const std::string &f) { return JFJochXFELReader::CanReadPanelTable(f) || JFJochXFELReader::CanReadCXI(f) || JFJochXTCReader::CanRead(f); @@ -1861,7 +1905,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { || JFJochSMVReader::CanRead(input_file) || JFJochSFRMReader::CanRead(input_file)) { logger.Error("{} is one frame of a sweep, which names the whole sweep: give just one of its " "frames", input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (!all_inputs([&is_foreign](const std::string &f) { return !is_foreign(f) && !JFJochCBFReader::CanRead(f) && !JFJochMarCCDReader::CanRead(f) @@ -1874,7 +1918,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { r->ReadFile(f); } catch (const std::exception &e) { logger.Error("{} is not a Jungfraujoch/NXmx master that can be read: {}", f, e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } inputs.push_back(std::move(r)); } @@ -1893,7 +1937,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (!extra_masks.empty() && !input_is_crystfel && !input_is_panel_table && !input_is_cxi) { logger.Error("--mask-good and --mask-bad are for HDF5/CXI stills (SACLA, CXI or --crystfel-geometry " "input); {} is not one", input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } xfel_reader.ExtraMasks(extra_masks); if (input_is_crystfel) { @@ -1901,7 +1945,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { reader_ptr = &xfel_reader; if (!xfel_reader.HasPhotonEnergy() && !wavelength_A) { logger.Error("Nothing read states the photon energy of {}; give it with --wavelength", input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } logger.Info("Read {} frames from {} files with the CrystFEL geometry {}", xfel_reader.GetNumberOfImages(), input_files.size(), crystfel_geometry); @@ -1915,7 +1959,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { logger.Info("Read {} frames from {} CXI files", xfel_reader.GetNumberOfImages(), input_files.size()); } else if (const auto refusal = JFJochXFELReader::CXIRefusal(input_file)) { logger.Error("{}", *refusal); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } input_is_cbf = !reader_ptr && JFJochCBFReader::CanRead(input_file); input_is_marccd = !reader_ptr && !input_is_cbf && JFJochMarCCDReader::CanRead(input_file); @@ -1930,7 +1974,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { && !input_is_cbf) { logger.Error("--gain is for input in detector units (SACLA, CXI, XTC, SMV, marCCD, SFRM, a CCD imgCIF or " "--crystfel-geometry input); {} is read in the units it states", input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (reader_ptr) { // read above @@ -1972,7 +2016,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { logger.Error("The XTC does not state the detector distance{}{}: give {} (the beam centre is " "measured; with -C the distance is too)", xtc_reader.HasPixelSize() ? "" : ", the pixel size", xtc_reader.HasWavelength() ? "" : ", the wavelength", options); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } reader_ptr = &xtc_reader; logger.Info("Read {} XTC still images from {} streams", xtc_reader.GetNumberOfImages(), @@ -1981,13 +2025,13 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // The collection is still being written: follow it from the broker. if (mode != RugnuxMode::MX) { logger.Error("Following a collection from the broker is for --mode mx only"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } const char *token = std::getenv("JUNGFRAUJOCH_HTTP_TOKEN"); if (token == nullptr) { logger.Error("{} does not exist yet, and following it from the broker needs the collection's token " "in the environment variable JUNGFRAUJOCH_HTTP_TOKEN", input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } broker_feed = std::make_unique(broker_url, token); // The broker's event stream starts with the start message: at once on a JUNGFRAU, when the @@ -1999,11 +2043,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { const std::string given = std::filesystem::path(input_file).lexically_normal().generic_string(); if (given != master && !given.ends_with("/" + master)) { logger.Error("The broker at {} is writing {}, not {}", broker_feed->Host(), master, input_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (!start.goniometer.has_value() || !start.goniometer->IsScanning()) { logger.Error("Following a collection from the broker is for rotation data, and {} is not", master); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } live_reader.Open(start, input_file); broker_feed->Attach(live_reader); @@ -2018,7 +2062,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } } catch (const std::exception &e) { logger.Error("Error reading input file: {}", e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } JFJochReader &reader = *reader_ptr; @@ -2085,20 +2129,19 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { ? fmt::format(", {:.3f} of them to a photon; intensities and sigmas in photons", adu_per_photon) : std::string(", one to a photon")); - const auto dataset = reader.GetDataset(); + dataset = reader.GetDataset(); if (!dataset) { logger.Error("No experiment dataset found in the input file"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (rotation_indexing_image_count <= 0) { logger.Error("Invalid number of rotation indexing images: {}", rotation_indexing_image_count); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } logger.Info("Loaded dataset from {}", input_file); // Recorded in the reflection files. - OutputProvenance output_provenance; output_provenance.command_line = provenance.command_line; DescribeInput(input_file, output_provenance); output_provenance.collection_date = dataset->arm_date; @@ -2192,7 +2235,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } } catch (const std::exception &e) { logger.Error("Error reading reference {}: {}", ref_mtz, e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } } @@ -2203,7 +2246,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { const auto flags = LoadReferenceMtz(free_flags_file, std::nullopt); if (!flags.has_free_flags) { logger.Error("{} carries no R-free flags (no FreeR_flag-like column, no _refln.status)", free_flags_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } free_flags_data = flags.reflections; logger.Info("R-free flags from {} (column {}): {} of {} free, taken by the merge where these data match " @@ -2211,7 +2254,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { flags.reflections.size()); } catch (const std::exception &e) { logger.Error("Error reading {}: {}", free_flags_file, e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } } @@ -2231,7 +2274,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { || input_is_xtc) { logger.Error("--mode scale needs the integrated reflections in a _process.h5; " "a raw image sweep holds none. Run --mode mx on it first."); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } const auto total_images = static_cast(reader.GetNumberOfImages()); const int last_image = (end_image < 0 || end_image >= total_images) ? total_images - 1 : end_image; @@ -2300,7 +2343,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (!experiment.GetUnitCell()) { logger.Error("Experiment unit cell not found, cannot update reflection resolution"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{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 @@ -2342,7 +2385,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "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); + throw RugnuxInputFailed{EXIT_FAILURE}; } } @@ -2366,7 +2409,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "has, or reprocess the images with --mode mx.", merge_sg->xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, violation, MAX_METRIC_VIOLATION); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } } @@ -2804,12 +2847,12 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (image_stride < 0) { logger.Error("Image stride cannot be negative"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (image_stride == 0) { logger.Error("Image stride cannot be zero"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } // Round UP: the worker loop walks ordinals until start + ordinal * stride reaches end_image, so a @@ -2824,7 +2867,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } // 2. Setup Experiment & Components - DiffractionExperiment experiment(dataset->experiment); + experiment = DiffractionExperiment(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. @@ -2860,20 +2903,20 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } if (modules.empty()) { logger.Error("--module-geometry {} describes fewer than two modules", module_geometry_file); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } for (const auto &m: modules.Get()) if (m.x0 < 0 || m.y0 < 0 || m.x0 + m.width > experiment.GetXPixelsNum() || m.y0 + m.height > experiment.GetYPixelsNum()) { logger.Error("--module-geometry {}: a module lies outside the {} x {} image", module_geometry_file, experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } experiment.Modules(modules); logger.Info("--module-geometry: {} modules from {}", modules.size(), module_geometry_file); } catch (const std::exception &e) { logger.Error("Cannot read --module-geometry {}: {}", module_geometry_file, e.what()); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } } if (!experiment.GetModules().empty()) @@ -2942,7 +2985,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { result = process.Run(console.get()); } catch (const std::exception &e) { logger.Error("Processing failed: {}", FailureText(e)); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } g_active_process = nullptr; @@ -3228,7 +3271,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { if (force_still) { if (rotation_indexing) { logger.Error("--force-still conflicts with -R / --single-pass-rotation / --force-rotation-lattice"); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } if (has_goniometer) logger.Info("--force-still: treating the rotation dataset as independent stills"); @@ -3250,7 +3293,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { "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); + throw RugnuxInputFailed{EXIT_FAILURE}; } // Scaling and merging are on by default (run_scaling initialised true); --no-merge turns them off @@ -3510,7 +3553,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // 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.provenance = output_provenance; config.start_image = start_image; @@ -3563,7 +3605,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.refine_solvent = refine_solvent; // 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. - std::optional sink; if (!output_prefix.empty() && (run_scaling ? write_process_h5_flag : true)) config.image_sink = &sink.emplace(output_prefix); config.finalist_ledger = finalist_ledger_flag; @@ -3579,6 +3620,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { WarmUpModelValidation(model, static_cast(config.nthreads), logger); simulated_start = std::chrono::steady_clock::now(); + return RUGNUX_CONTINUE; +} + +int RugnuxInputRun::RunPreparedInput() { + JFJochReader &reader = *reader_ptr; Rugnux process(reader, experiment, *dataset->pixel_mask, config); g_active_process = &process; @@ -3586,7 +3632,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { 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. @@ -3617,13 +3665,17 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { provenance); console->Fail(); PrintSummary(*console, doc, output_prefix, files_since, provenance.wall_time_s); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } catch (const std::exception &e) { logger.Error("Processing failed: {}", FailureText(e)); sink.reset(); - exit(EXIT_FAILURE); + 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) { @@ -3641,7 +3693,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { logger); } catch (const std::exception &e) { logger.Error("Processing failed: {}", FailureText(e)); - exit(EXIT_FAILURE); + throw RugnuxInputFailed{EXIT_FAILURE}; } } @@ -3785,6 +3837,304 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { return EXIT_SUCCESS; } +int RugnuxCli::Run() { + // The sweeps the inputs are (RugnuxSweeps.h). One input, a stills run and anything the readers refuse + // are read as they were. + const auto describe_file = [this](const std::string &f) { + RugnuxInputDescription d; + if (!crystfel_geometry.empty()) + d.foreign = true; + else if (JFJochCBFReader::CanRead(f) || JFJochMarCCDReader::CanRead(f) || JFJochSMVReader::CanRead(f)) + d.frame_template = RugnuxFrameTemplate(f); + else + d.foreign = JFJochXFELReader::CanReadPanelTable(f) || JFJochXFELReader::CanReadCXI(f) + || JFJochXTCReader::CanRead(f); + return d; + }; + // What a sweep is read as, for its rotation and wavelength: the reader the run itself would take. + const auto read_dataset = [](const std::string &f) -> std::shared_ptr { + if (JFJochCBFReader::CanRead(f)) { + JFJochCBFReader r; + r.ReadFiles(f); + return r.GetDataset(); + } + if (JFJochMarCCDReader::CanRead(f)) { + JFJochMarCCDReader r; + r.ReadFiles(f); + return r.GetDataset(); + } + if (JFJochSMVReader::CanRead(f)) { + JFJochSMVReader r; + r.ReadFiles(f); + return r.GetDataset(); + } + JFJochHDF5Reader r; + r.ReadFile(f); + return r.GetDataset(); + }; + const auto describe_sweep = [&read_dataset](const std::string &f) { + std::shared_ptr dataset; + try { + dataset = read_dataset(f); + } catch (const std::exception &e) { + throw std::runtime_error(fmt::format("{} cannot be read: {}", f, e.what())); + } + RugnuxInputDescription d; + if (dataset) { + const auto g = dataset->experiment.GetGoniometer(); + d.rotation = g.has_value() && g->GetIncrement_deg() != 0.0f; + d.wavelength_A = dataset->experiment.GetWavelength_A(); + } + return d; + }; + RugnuxSweepPlan plan; + try { + plan = PlanRugnuxSweeps(input_files, describe_file, describe_sweep); + } catch (const std::exception &e) { + logger.Error("{}", e.what()); + return EXIT_FAILURE; + } + + if (plan.sweeps.size() < 2) { + // Frames of one sweep: any of them names it. + if (plan.sweeps.size() == 1) + input_files = {plan.sweeps[0].input}; + try { + return RugnuxInputRun(*this).Process(); + } catch (const RugnuxInputFailed &f) { + return f.code; + } + } + + for (size_t k = 0; k < plan.sweeps.size(); k++) + logger.Info("Sweep {}: {}", k + 1, plan.sweeps[k].label); + for (const auto &s: RugnuxSiblingSweeps(plan, {".cbf", ".img", ".mccd", ".osc"})) + logger.Info("Also beside these inputs, not given and so not processed: {}", s); + if (start_image != 0 || end_image >= 0) { + logger.Error("-s / -e select images of one sweep; with several sweeps give each sweep's frames instead"); + return EXIT_FAILURE; + } + if (force_still) { + 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) + return RunJointSweeps(plan); + + // --sweeps separate: each sweep alone, as if it had been the one input, into _s. + int failed = 0; + for (size_t k = 0; k < plan.sweeps.size(); k++) { + RugnuxInputRun sweep(*this); + sweep.input_files = {plan.sweeps[k].input}; + // Each sweep's report times that sweep; the first also the start-up before it. + if (k > 0) + sweep.invocation_start = std::chrono::steady_clock::now(); + if (!output_prefix.empty()) + sweep.output_prefix = output_prefix + "_s" + std::to_string(k + 1); + logger.Info("Sweep {} of {}: {} into {}_*", k + 1, plan.sweeps.size(), plan.sweeps[k].label, + sweep.output_prefix); + int rc; + try { + rc = sweep.Process(); + } catch (const RugnuxInputFailed &f) { + rc = f.code; + } + if (sweep.log_file_sink) + log_sinks->remove_sink(sweep.log_file_sink); + if (rc != EXIT_SUCCESS) { + logger.Error("Sweep {} ({}) failed", k + 1, plan.sweeps[k].label); + failed++; + } + } + if (failed > 0) + logger.Error("{} of {} sweeps failed", failed, plan.sweeps.size()); + return failed > 0 ? EXIT_FAILURE : EXIT_SUCCESS; +} + + +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 +// process with no diagnostic and exit code 134. +static int RunRugnux(int argc, char **argv, bool no_gpu) { + // cuFFT plans are reused through the run and destroyed with it, ahead of the device reset in main. + FFTPlanReuse fft_plans; + // 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(); +#ifdef __APPLE__ + // Without CUDA, the Bragg integration's integer pixel work runs on the Apple GPU (MXAnalysisWithoutFPGA). + const std::string metal_device = BraggIntegrationEngineMetal::DeviceName(); +#else + const std::string metal_device; +#endif + provenance.gpu_description = get_gpu_description(); + + // The terminal shows the run as steps (RugnuxConsole) and every Logger writes to the log file + // instead; with -v the whole log is on the terminal as well. -v is looked for here, before the + // options are parsed, because the parser already logs. + bool verbose = false; + for (int i = 1; i < argc; i++) + if (std::strcmp(argv[i], "-v") == 0 || std::strcmp(argv[i], "--verbose") == 0) + verbose = true; + console = std::make_shared(verbose); + std::atexit(CloseConsoleOnExit); + // Files with an older time than this are not this run's (the margin covers coarse file times). + const auto files_since = std::filesystem::file_time_type::clock::now() - std::chrono::seconds(2); + auto log_sinks = std::make_shared(); + // What is logged before the log file can be opened - its name comes from -o - waits here. + auto early_log = std::make_shared(10000); + if (verbose) + log_sinks->add_sink(std::make_shared()); + else + log_sinks->add_sink(MakeConsoleSink(console)); + log_sinks->add_sink(early_log); + Logger::RouteAll(log_sinks); + // The HDF5 library prints its own error stack on stderr; every failure it reports reaches the + // reader as an error code anyway, which is turned into a message there. + H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr); + + PrintBanner(*console); + + Logger logger("rugnux"); +#ifdef JFJOCH_USE_CUDA + const char *build = "CUDA build"; +#else + const char *build = "CPU build"; +#endif + logger.Info("rugnux Version {}{} - {}", jfjoch_version(), + jfjoch_git_sha1().empty() ? std::string() + : fmt::format(" (git {} {})", jfjoch_git_sha1().substr(0, 6), jfjoch_git_date()), build); + logger.Info("Build flags: {}", jfjoch_build_cxx_flags().empty() ? "none" : jfjoch_build_cxx_flags()); + logger.Info("Build: {}", jfjoch_build_info(true)); + logger.Info("Command line: {}", provenance.command_line); + logger.Info("GPU: {}", no_gpu ? std::string("not used (--no-gpu) - running on the CPU") + : provenance.gpu_count > 0 ? fmt::format("{} ({} visible)", provenance.gpu_description, provenance.gpu_count) + : !metal_device.empty() ? fmt::format("{} through Metal, for the Bragg integration only - the rest runs on the CPU", metal_device) + : std::string("none visible - running on the CPU")); + RegisterHDF5Filter(); + + RugnuxCli cli(no_gpu, invocation_start, provenance, metal_device, files_since, log_sinks, early_log); + if (const int rc = cli.Parse(argc, argv); rc != RUGNUX_CONTINUE) + return rc; + return cli.Run(); +} + int main(int argc, char **argv) { try { // Before anything creates a CUDA context, which fixes the flags for the device. A thread diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c220dd46c..4e7b73abe 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -141,6 +141,7 @@ ADD_EXECUTABLE(jfjoch_test PostRefineTest.cpp CellEsdTest.cpp RotationScaleWalkTest.cpp + RotationScaleMergeSegmentsTest.cpp AndersonMixingTest.cpp RugnuxLargeTest.cpp TestData.h @@ -214,6 +215,8 @@ ADD_EXECUTABLE(jfjoch_test WriteReflectionsTest.cpp RugnuxResultFilesTest.cpp RugnuxInputKindTest.cpp + RugnuxSweepsTest.cpp + RugnuxJointSweepsTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/RotationScaleMergeSegmentsTest.cpp b/tests/RotationScaleMergeSegmentsTest.cpp new file mode 100644 index 000000000..f20d05f49 --- /dev/null +++ b/tests/RotationScaleMergeSegmentsTest.cpp @@ -0,0 +1,231 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include +#include +#include + +#include "../common/CUDAWrapper.h" +#include "../common/Logger.h" +#include "../image_analysis/scale_merge/RotationScaleMerge.h" + +namespace { + +// A synthetic crystal measured as pre-combined fulls on several sweeps: a cubic 20 A cell, Wilson +// intensities at B = 15 A^2 to 1.6 A, each sweep at its own scale and relative B, Poisson-like noise. +struct SweepSpec { + int n_frames; + double scale; + double relative_b; // at the sweep's first frame + double relative_b_end = NAN; // at its last; NaN: the same B throughout +}; + +struct Synthetic { + std::vector frames; + std::map, double> I_true; // by the Friedel-merged index +}; + +constexpr double CELL = 20.0; + +std::array FriedelKey(int h, int k, int l) { + const std::array a{h, k, l}, b{-h, -k, -l}; + return std::max(a, b); +} + +Synthetic MakeSweeps(const std::vector &sweeps, uint32_t seed) { + std::mt19937 rng(seed); + std::vector> hkl; + for (int h = -12; h <= 12; ++h) + for (int k = -12; k <= 12; ++k) + for (int l = -12; l <= 12; ++l) { + const double d = CELL / std::sqrt(static_cast(h * h + k * k + l * l)); + if ((h || k || l) && d >= 1.6) hkl.push_back({h, k, l}); + } + Synthetic out; + std::exponential_distribution wilson(1.0); + for (const auto &m: hkl) { + const auto key = FriedelKey(m[0], m[1], m[2]); + if (!out.I_true.contains(key)) { + const double s2 = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]) / (4.0 * CELL * CELL); + out.I_true[key] = 2000.0 * std::exp(-2.0 * 15.0 * s2) * (0.05 + wilson(rng)); + } + } + std::uniform_int_distribution pick(0, hkl.size() - 1); + std::normal_distribution noise(0.0, 1.0); + int frame = 0; + for (const auto &sweep: sweeps) + for (int f = 0; f < sweep.n_frames; ++f, ++frame) { + IntegrationOutcome o; + o.latt = CrystalLattice(UnitCell{static_cast(CELL), static_cast(CELL), + static_cast(CELL), 90, 90, 90}); + for (int j = 0; j < 400; ++j) { + const auto &m = hkl[pick(rng)]; + const double s2 = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]) / (4.0 * CELL * CELL); + const double b = std::isfinite(sweep.relative_b_end) + ? sweep.relative_b + (sweep.relative_b_end - sweep.relative_b) * f / (sweep.n_frames - 1) + : sweep.relative_b; + const double expected = out.I_true.at(FriedelKey(m[0], m[1], m[2])) * sweep.scale + * std::exp(-2.0 * b * s2); + const double sigma = std::sqrt(expected + 20.0); + Reflection r{}; + r.h = m[0]; r.k = m[1]; r.l = m[2]; + r.image_number = static_cast(frame); + r.d = static_cast(1.0 / std::sqrt(4.0 * s2)); + r.I = static_cast(expected + sigma * noise(rng)); + r.sigma = static_cast(sigma); + // The combine rebuilds a full's variance as var_bkg + I: what is not the intensity's own + r.var_bkg = static_cast(sigma * sigma - std::max(0.0, static_cast(r.I))); + r.bkg = 20.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + r.zeta = 1.0f; + r.image_scale_corr = 1.0f; + r.predicted_x = 500.0f; + r.predicted_y = 500.0f; + o.reflections.push_back(r); + } + out.frames.push_back(std::move(o)); + } + return out; +} + +DiffractionExperiment SyntheticExperiment(bool correction_surfaces) { + DiffractionExperiment x; + x.IncidentEnergy_keV(12.4f); + x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), {})); + x.SetUnitCell(UnitCell{static_cast(CELL), static_cast(CELL), static_cast(CELL), 90, 90, 90}); + ScalingSettings s; + s.ScaleFulls(true).CorrectionSurfaces(correction_surfaces).ResolutionCutoff(ResolutionCutoffMethod::Off); + x.ImportScalingSettings(s); + return x; +} + +std::vector Segments(const std::vector &sweeps) { + std::vector out; + int32_t first = 0; + for (const auto &s: sweeps) { + out.push_back({first, s.n_frames}); + first += s.n_frames; + } + return out; +} + +RotationScaleMerge::Result Merge(Synthetic &data, const DiffractionExperiment &x, + const std::vector &segments, bool export_fulls = false, + bool *on_gpu = nullptr) { + Logger logger("RotationScaleMergeSegmentsTest"); + const UnitCell cell{static_cast(CELL), static_cast(CELL), static_cast(CELL), 90, 90, 90}; + RotationScaleMerge engine(x, data.frames, cell, 20, 4, logger); + if (!segments.empty()) + engine.SetSweepSegments(segments); + engine.SetExportScaledFulls(export_fulls); + engine.Ingest(); + if (on_gpu) + *on_gpu = engine.OnGPU(); + return engine.Run(false, true, false); +} + +// The median over a frame range of the scaled fulls' ratio to the true intensity (the strong ones). +double MedianRatio(const RotationScaleMerge::Result &r, const Synthetic &data, int f0, int f1) { + std::vector ratio; + for (const auto &o: r.scaled_fulls) + if (o.image_number >= static_cast(f0) && o.image_number < static_cast(f1)) { + const double t = data.I_true.at(FriedelKey(o.h, o.k, o.l)); + if (t > 500.0) + ratio.push_back(o.I / t); + } + REQUIRE(ratio.size() > 20); + std::nth_element(ratio.begin(), ratio.begin() + ratio.size() / 2, ratio.end()); + return ratio[ratio.size() / 2]; +} + +} // namespace + +TEST_CASE("RotationScaleMergeSegments_PerSweepScaleRecovered", "[RotationScale]") { + const std::vector sweeps{{60, 1.0, 0.0}, {60, 0.5, 0.0}}; + auto data = MakeSweeps(sweeps, 11); + const auto r = Merge(data, SyntheticExperiment(false), Segments(sweeps), true); + REQUIRE(r.segments.size() == 2); + CHECK(r.segments[1].scale / r.segments[0].scale == Catch::Approx(0.5).epsilon(0.03)); + CHECK(r.statistics.overall.cc_half > 0.99); + + // No coupling across the boundary: the frames on either side of it are put on the common scale as + // well as the middle of their sweeps are. + const double middle = MedianRatio(r, data, 20, 40); + CHECK(MedianRatio(r, data, 57, 60) / middle == Catch::Approx(1.0).epsilon(0.03)); + CHECK(MedianRatio(r, data, 60, 63) / middle == Catch::Approx(1.0).epsilon(0.03)); + CHECK(MedianRatio(r, data, 80, 100) / middle == Catch::Approx(1.0).epsilon(0.03)); +} + +TEST_CASE("RotationScaleMergeSegments_DoseOrderFromTheRelativeB", "[RotationScale]") { + // Given first, the sweep taken second: its B is 8 A^2 higher. The order comes from the data. + const std::vector sweeps{{60, 1.0, 8.0}, {60, 1.0, 0.0}, {60, 1.0, 4.0}}; + auto data = MakeSweeps(sweeps, 12); + const auto r = Merge(data, SyntheticExperiment(true), Segments(sweeps)); + REQUIRE(r.segments.size() == 3); + CHECK(r.segments[1].dose_rank == 0); + CHECK(r.segments[2].dose_rank == 1); + CHECK(r.segments[0].dose_rank == 2); + CHECK(r.segments[1].relative_b == Catch::Approx(0.0).margin(1e-9)); + CHECK(r.segments[0].relative_b == Catch::Approx(8.0).margin(1.0)); + CHECK(r.segments[2].relative_b == Catch::Approx(4.0).margin(1.0)); + // Flat within each sweep: one B per sweep predicts the held-out frames better than a dose-linear ramp + CHECK(r.segment_decay_model == "per-sweep B"); +} + +TEST_CASE("RotationScaleMergeSegments_DoseLinearDecayAcrossSweeps", "[RotationScale]") { + // One crystal damaged steadily: the sweep given first was taken second, its B going on from where the + // other's ended. Along the dose order the B rises in a straight line, which the ramp describes. + const std::vector sweeps{{60, 1.0, 6.0, 12.0}, {60, 1.0, 0.0, 6.0}}; + auto data = MakeSweeps(sweeps, 15); + const auto r = Merge(data, SyntheticExperiment(true), Segments(sweeps)); + REQUIRE(r.segments.size() == 2); + CHECK(r.segments[1].dose_rank == 0); + CHECK(r.segments[0].dose_rank == 1); + CHECK(r.segments[0].relative_b == Catch::Approx(6.0).margin(1.5)); + CHECK(r.segment_decay_model == "dose-linear"); +} + +TEST_CASE("RotationScaleMergeSegments_SegmentsMustCoverTheFrames", "[RotationScale]") { + const std::vector sweeps{{10, 1.0, 0.0}, {10, 1.0, 0.0}}; + auto data = MakeSweeps(sweeps, 13); + const auto x = SyntheticExperiment(false); + CHECK_THROWS_AS(Merge(data, x, {{0, 10}, {11, 9}}), JFJochException); + CHECK_THROWS_AS(Merge(data, x, {{0, 10}}), JFJochException); +} + +TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[RotationScale]") { + // The same synthetic sweep merged without segments gives, to the bit, what the engine gave before + // segments existed: FNV-1a over the merged I and sigma, recorded on the engine of rc181 e56b19abb - + // on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host + // path to the bit only where every reduction does; this one carries its own, so it is checked against + // its own record). + const std::vector sweeps{{40, 1.0, 0.0}, {40, 0.7, 3.0}}; + auto data = MakeSweeps(sweeps, 14); + bool on_gpu = false; + const auto r = Merge(data, SyntheticExperiment(true), {}, false, &on_gpu); + uint64_t hash = 1469598103934665603ULL; + auto mix = [&hash](float v) { + uint32_t u; + std::memcpy(&u, &v, sizeof u); + for (int b = 0; b < 4; ++b) { + hash ^= (u >> (8 * b)) & 0xffu; + hash *= 1099511628211ULL; + } + }; + for (const auto &m: r.merged) { + mix(m.I); + mix(m.sigma); + } + CHECK(r.segments.empty()); + CHECK(r.merged.size() == 4106); + if (!on_gpu) + CHECK(hash == 3376995616949734726ULL); + else if (get_gpu_description().find("RTX 5080") != std::string::npos) + CHECK(hash == 6594155396455850264ULL); +} 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); +} diff --git a/tests/RugnuxSweepsTest.cpp b/tests/RugnuxSweepsTest.cpp new file mode 100644 index 000000000..a46541304 --- /dev/null +++ b/tests/RugnuxSweepsTest.cpp @@ -0,0 +1,132 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include + +#include "../rugnux/RugnuxSweeps.h" + +namespace { + +// Inputs described by name alone: "