From e5c0066129e82f447081fe44656b09bb56bedb27 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 11 Aug 2026 18:03:07 +0200 Subject: [PATCH] rugnux: write a results report next to the reflections Everything a run determines went to stdout and nowhere else. The space group and the evidence behind it, the error model, the post-refine commit-or-reject decisions and their held-out residuals, the two-pass adopt-or-roll-back, the resolution cut, the merging statistics - all of it scrolled past interleaved with progress lines and was gone. A user who was not watching had no record, and nothing could read it. `rugnux` had no log file at all; the `rugnux.log` in the regression harness is that harness capturing stdout. Write `_report.txt` alongside the .cif/.mtz/.hkl, always, with no option to ask for it. It holds what the run DETERMINED; timing, rates, per-image progress and engine chatter stay on stdout, where they belong. Every line rugnux logs was classified result-or-process against the regression corpus to decide what crosses over. The format follows XDS's CORRECT.LP, which has been read by people and parsed by other programs for twenty years: `KEY= value` assignment lines a script greps one at a time, fixed-width tables with stable headers and a total row, `WARNING:` sentences in plain English, section banners. REPORT_VERSION says when that interface last changed. It is assembled from results the pipeline already computed, so an unconditional file costs nothing, and a failure to write it is logged and swallowed - a run that produced good reflections must not be lost to a side file. One thing CORRECT.LP does not have to solve: a rotation run integrates twice and writes both passes, so every report says which pass it describes and why that pass was adopted. `--no-merge` gets a report too, saying MERGE= NOT_PERFORMED rather than leaving a reader to infer it from absent sections. An empty output prefix still writes nothing. Co-Authored-By: Claude Opus 5 (1M context) --- common/JFJochMessages.h | 8 + docs/CHANGELOG.md | 1 + docs/HDF5.md | 15 ++ docs/RUGNUX.md | 96 +++++++++ image_analysis/scale_merge/Merge.cpp | 22 +++ image_analysis/scale_merge/Merge.h | 44 +++++ reader/HDF5MetadataSource.cpp | 8 + reader/JFJochReaderDataset.h | 7 + rugnux/CMakeLists.txt | 2 + rugnux/ResultReport.cpp | 283 +++++++++++++++++++++++++++ rugnux/ResultReport.h | 33 ++++ rugnux/Rugnux.cpp | 54 +++++ rugnux/Rugnux.h | 26 +++ rugnux/rugnux_cli.cpp | 28 ++- tests/CMakeLists.txt | 1 + tests/ResultReportTest.cpp | 160 +++++++++++++++ writer/HDF5NXmx.cpp | 8 + 17 files changed, 795 insertions(+), 1 deletion(-) create mode 100644 rugnux/ResultReport.cpp create mode 100644 rugnux/ResultReport.h create mode 100644 tests/ResultReportTest.cpp diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 1c1787ac..9a8989e1 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -399,6 +399,14 @@ struct EndMessage { // Run mean of the above, the single "how icy was this dataset" number (1 = no ice). The // bkg_estimate scalar's counterpart; written to /entry/MX/iceRingScoreMean. std::optional ice_ring_score_mean; + + // Per-image sweep-quality code: 0 = the image falls in no flagged range, otherwise 1 + the + // SweepQualityReason of the range it falls in (image_analysis/scale_merge/Merge.h). Written to + // /entry/MX/sweepQuality, with the vocabulary alongside it in /entry/MX/sweepQualityReasons. + // Empty when the sweep-quality diagnostic did not run, so nothing could be looked for - which is + // not the same as every image being clean. Offline (rugnux) only; the broker does not merge. + std::vector sweep_quality; + std::vector sweep_quality_reasons; }; struct MetadataMessage { diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 0145dbbd..379090b5 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: Every `mx` and `scale` run now writes `_report.txt`, a results report modelled on XDS's `CORRECT.LP` — `KEY= value` lines, fixed-width tables and `WARNING:` sentences covering indexing, geometry post-refinement, the space-group decision, merging, twinning, radiation damage and the stretches of the sweep over which the crystal delivered much less than the rest of the run. * Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. * Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. * Scaling: the rotation merge weights each combined full by its variance rebuilt at the reflection's mean intensity rather than by the full's own sigma, as the stills merge already did. diff --git a/docs/HDF5.md b/docs/HDF5.md index f9f2deaf..fad40276 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -262,6 +262,7 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i | `imageScaleFactor` | | on-the-fly per-image scale factor *g* | | `imageScaleCC` | | on-the-fly scaling correlation coefficient | | `imageScaleMosaicity` | deg | scaling-model mosaicity | +| `sweepQuality` | | why this image's stretch of the sweep was flagged — see below | **Per-image lattices:** `latticeIndexed` `[n_images, 9]` (Å) — the real-space lattice (flattened 3×3); `latticeIndexedExtra` `[n_images, max_extra_lattices, 9]` (Å) — additional orientation @@ -281,6 +282,20 @@ variants. | `iceRingScoreMean` | ratio | mean `iceRingScore` over the run — the single "how icy was this dataset" number (1 = no ice) | | `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* | +**Sweep quality.** `sweepQuality` `[n_images]` (`uint8`) says why the stretch of the sweep this +image belongs to was flagged as delivering much less than the rest of the run: **0** means it was +not, and any other value is a **1-based index into `sweepQualityReasons`**, a string vector written +beside it that carries the whole vocabulary, so the codes can be read without this source. The +vocabulary is closed and stable — a code is never renamed and never reused — and currently reads +`no_diffraction`, `crystal_out_of_beam`, `weak_diffraction`, `loss_of_centring`, `radiation_damage`; +[the rugnux documentation](RUGNUX.md#sweep-quality-and-the-reason-vocabulary) defines what each one +means. Both datasets are **absent** unless the sweep-quality diagnostic ran, which needs scaling and +merging; their absence therefore means "not looked for", *not* "every image clean". Written by the +offline `rugnux` path only — the broker does not merge — and not carried on the CBOR stream, in the +same way as the other offline-only fields (`space_group_number`, the refined geometry). Nothing is +excluded from processing on the strength of it. The condensed, dataset-wide form of the same finding +is in `_report.txt`. + CrystFEL can read the spots directly with: ``` diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index bce98848..71e661f1 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -60,6 +60,9 @@ reusing them would hide the spot-finding settings from the lattice search. All three carry the **refined unit cell** (from rotation indexing) and the **space group determined from systematic absences** (constrained to the indexed lattice symmetry). No-reference scaling additionally emits per-iteration `_iterN_scale.dat`. +- `_report.txt` — the **results report**: what the run determined, in a form both a person and + a beamline script can read. Always written, next to the files above. See + [The results report](#the-results-report) below. Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the console. By default the written resolution is trimmed automatically where CC1/2 falls off @@ -103,6 +106,99 @@ expect. Two properties worth knowing before using it: untouched, but the absolute scale is not meaningful. This matters only if you intend to compare magnitudes with another file; SHELXC and ANODE use ratios alone. +## The results report + +`_report.txt` records **what the run determined**, next to the reflection files. It is +written on every `--mode mx` and `--mode scale` run that has an output prefix — there is no option +to enable or disable it. Two cases follow from that: + +- An **empty output prefix** (`-o ""`, the "compute the statistics, persist nothing" mode) writes + nothing, the report included. +- **`--no-merge`** still writes a report. It determined an indexing and a geometry result, and those + are recorded; the merging section then says `MERGE= NOT_PERFORMED` rather than being omitted, so + the absence is a statement and not something a reader has to infer. + +The report is never allowed to fail a run: if it cannot be written (unwritable path, full disk) the +failure is logged as a warning and the run finishes normally. + +### Format + +The model is XDS's `CORRECT.LP`: prose and tables a crystallographer reads top to bottom, with a +structure a script can consume without parsing prose. + +- **`KEY= value` assignment lines.** Every number worth extracting is one, so a consumer gets it with + a single `grep '^ISA= '` and never has to read a sentence. Key names are stable. +- **Fixed-width tables** with a stable header row for anything that is genuinely tabular — the + resolution shells, the space-group candidates, the sweep-quality ranges. +- **`WARNING:` lines**, one per finding, in plain English: `WARNING: Frames 500-600 out of beam + (10.1 deg, scale 0.12 and CC 0.30 of the run, 2% scaled)`. `grep '^WARNING:'` finds every one. +- **Section banners** (`***…***` around a numbered title) delimiting the blocks. + +`REPORT_VERSION=` is the format's own version. Key names, table columns and the reason vocabulary +below are an interface other software may depend on: they do not change without that number moving. + +Sections, in order: `1. DATA SET`, `2. INDEXING`, `3. GEOMETRY POST-REFINEMENT` (rotation only), +`4. SPACE GROUP DETERMINATION`, `5. SCALING AND MERGING`, `6. TWINNING`, `7. RADIATION DAMAGE`, +`8. SWEEP QUALITY`, `9. WARNINGS`. + +**Which pass.** A rotation run integrates twice — once at the geometry in the input file +(`_01.*`), then again at the post-refined geometry (`.*`) — and can integrate a +third time if a guard rejects the second pass. There is **one** report, for the pass that became the +canonical output, and `PASS=` / `PASS_DECISION=` in section 1 say which pass that is and on what +evidence, so no number in the file is ambiguous about which geometry produced it. + +**Not in the report:** timing, frame rates, thread counts, per-image progress and library banners. +Those are process, not result, and stay on stdout. + +### Sweep quality and the reason vocabulary + +Section 8 lists the stretches of the sweep over which the crystal delivered much less than the rest +of the run — the feedback a beamline control system needs to tell an operator that a crystal should +be recentred or recollected. Nothing is excluded on the strength of it; the frames still carry +signal, and this is a message for the beamline, not a filter. + +``` +SWEEP_QUALITY_STATUS= COMPUTED +SWEEP_QUALITY_COUNT= 1 +SWEEP_QUALITY_REASONS= no_diffraction crystal_out_of_beam weak_diffraction loss_of_centring radiation_damage +SWEEP_ROTATION= 360.0 +FLUX_PEAK_TO_TROUGH= 1.03 +SCALE_MODULATION_PEAK_TO_TROUGH= 1.00 + + FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED + ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- + 500 600 101 10.1 crystal_out_of_beam 0.83 0.12 0.30 0.02 + ----------- ----------- --------- -------- -------------------- -------- ------ ------ -------- +``` + +`SWEEP_QUALITY_STATUS` distinguishes **`COMPUTED`** (the diagnostic ran; a count of 0 means the sweep +was clean throughout) from **`NOT_COMPUTED`** (it did not run — no scaling and merging, or stills +data). A consumer must not read a missing table or a zero count as "clean" without checking it. +`SWEEP_QUALITY_REASONS` lists the whole vocabulary this version can emit, so an unknown code is +distinguishable from a missing one. + +| Reason code | Meaning | +|-------------|---------| +| `no_diffraction` | The range recorded essentially no diffraction from the indexed lattice. | +| `crystal_out_of_beam` | Frames were lost: over the range a per-image scale could be fitted far less often than over the run. | +| `weak_diffraction` | The frames all still index, but with much less intensity — the cause was not determined. | +| `loss_of_centring` | One cycle of modulation per revolution: the crystal is off the rotation axis. | +| `radiation_damage` | The range runs to the end of a sweep whose quality was already decaying. | + +The vocabulary is **closed and stable**: a code is never renamed, and never reused for a different +meaning. New codes are only ever added, and adding one moves `REPORT_VERSION`. + +The columns are: `FIRST_IMAGE`/`LAST_IMAGE` — inclusive, in processed-image ordinals (the numbering +of `_image.dat` and of every other per-image array `rugnux` writes; with `-s`/`--stride` the +source image is `start + ordinal * stride`); `ROTATION` — the width of the range in degrees; +`SEVERITY` — the fraction of the run's typical diffracting power missing over the range, 0 (as good +as the run) to 1 (nothing at all); `SCALE` and `CC` — the range's mean per-image scale and +CC-to-merge relative to the run median; `INDEXED` — the fraction of the range's frames that were +scaled at all. Every range also appears as a `WARNING:` sentence in section 9. + +The same finding is written **per image** into the `_process.h5` as `/entry/MX/sweepQuality`, when +one is written — see [HDF5](HDF5.md#41-entrymx--spot-finding-and-indexing-cxi-style). + ## Validating against a model (`rugnux --model`) Given a PDB atomic model of the same structure, `--model model.pdb` scales the model structure diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 7ad95166..d865b884 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -703,6 +703,28 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m return out; } +const char *SweepQualityReasonCode(SweepQualityReason reason) { + switch (reason) { + case SweepQualityReason::NoDiffraction: return "no_diffraction"; + case SweepQualityReason::CrystalOutOfBeam: return "crystal_out_of_beam"; + case SweepQualityReason::WeakDiffraction: return "weak_diffraction"; + case SweepQualityReason::LossOfCentring: return "loss_of_centring"; + case SweepQualityReason::RadiationDamage: return "radiation_damage"; + } + return "unknown"; +} + +const char *SweepQualityReasonText(SweepQualityReason reason) { + switch (reason) { + case SweepQualityReason::NoDiffraction: return "no diffraction"; + case SweepQualityReason::CrystalOutOfBeam: return "out of beam"; + case SweepQualityReason::WeakDiffraction: return "weak diffraction"; + case SweepQualityReason::LossOfCentring: return "loss of centring"; + case SweepQualityReason::RadiationDamage: return "radiation damage"; + } + return "unknown"; +} + std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { double completeness = in.possible_unique_reflections > 0 ? static_cast(in.unique_reflections) / in.possible_unique_reflections * 100.0 : 0.0; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 8bd626f5..0074937d 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -40,6 +40,46 @@ struct MergeStatisticsShell { double abs_diff_over_sigma_anomalous = NAN; }; +// Why a stretch of the sweep came out much weaker than the rest of the run. Report-only: nothing is +// excluded on the strength of it. The codes are the field's own words - "crystal rotating out of the +// beam" (HKL-2000 manual), "loss of centring during crystal rotation" (autoPROC). +enum class SweepQualityReason { + NoDiffraction, // the range recorded essentially no diffraction from the indexed lattice + CrystalOutOfBeam, // frames were lost: the range gets a per-image scale far less often than the run + WeakDiffraction, // the frames all still index, but with much less intensity - cause not determined + LossOfCentring, // one cycle of modulation per revolution: the crystal is off the rotation axis + RadiationDamage // the range runs to the end of a sweep whose quality was already decaying +}; + +const char *SweepQualityReasonCode(SweepQualityReason reason); // machine-readable, e.g. "out_of_beam" +const char *SweepQualityReasonText(SweepQualityReason reason); // for a sentence, e.g. "out of beam" + +struct SweepQualityRange { + // Inclusive, in processed-image ordinals - the numbering of _image.dat and of every other + // per-image array rugnux writes. With -s/--stride the source image is start + ordinal * stride. + int first_image = 0; + int last_image = 0; + SweepQualityReason reason = SweepQualityReason::WeakDiffraction; + // Fraction of the run's typical diffracting power missing over the range: 0 = as good as the run, + // 1 = nothing at all. The rest are supporting numbers, relative to the run unless stated. + float severity = 0.0f; + float rotation_deg = 0.0f; // width of the range + float mean_relative_scale = 1.0f; // / run median + float mean_relative_cc = 1.0f; // / run median + float indexed_fraction = 1.0f; // frames in the range that got a per-image scale at all + float relative_b = NAN; // mean of the radiation-damage monitor's per-batch B over the range +}; + +// Sweep-quality diagnostic (rotation only). `measured` separates "the run is clean" from "this did not +// run": the range list is empty in the first case and in the second alike. +struct SweepQuality { + bool measured = false; + float sweep_deg = 0.0f; + float flux_peak_to_trough = 1.0f; // the incident-flux proxy, over the whole run + float modulation_peak_to_trough = 1.0f; // depth of a DIAGNOSED once-per-revolution modulation (1 = none) + std::vector ranges; +}; + struct MergeStatistics { std::vector shells; MergeStatisticsShell overall; @@ -58,6 +98,10 @@ struct MergeStatistics { double radiation_damage_delta_b = NAN; std::vector radiation_damage_b_batch; double radiation_damage_batch_deg = 0.0; + + // Stretches of the sweep over which the crystal delivered much less than the rest of the run + // (MeasureSweepQuality). Report-only - no observation is dropped because of it. + SweepQuality sweep_quality; }; diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index b047b8d8..5d40edb5 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -389,6 +389,14 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen dataset->image_scale_factor = master_file->ReadOptVector("/entry/MX/imageScaleFactor"); dataset->image_scale_cc = master_file->ReadOptVector("/entry/MX/imageScaleCC"); dataset->integrated_reflections = master_file->ReadOptVector("/entry/MX/integratedReflections"); + dataset->sweep_quality = master_file->ReadOptVector("/entry/MX/sweepQuality"); + if (master_file->Exists("/entry/MX/sweepQualityReasons")) { + const auto dim = master_file->GetDimension("/entry/MX/sweepQualityReasons"); + for (size_t i = 0; i < (dim.empty() ? 0 : dim[0]); i++) + dataset->sweep_quality_reasons.push_back( + master_file->ReadElement("/entry/MX/sweepQualityReasons", i) + .value_or("")); + } } if (master_file->Exists("/entry/image")) dataset->max_value = master_file->ReadOptVector("/entry/image/max_value"); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index de2ebf49..b03c04e1 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -57,6 +57,13 @@ struct JFJochReaderDataset { std::vector image_scale_cc; std::vector max_value; + // Per-image sweep-quality code from /entry/MX/sweepQuality: 0 = the image is in no flagged range, + // otherwise a 1-based index into sweep_quality_reasons (the vocabulary stored beside it, so the + // codes read back without this source). Both empty when the file carries no sweep quality, which + // means the run never looked - not that every image was clean. + std::vector sweep_quality; + std::vector sweep_quality_reasons; + // Maps this dataset's image index -> the original/collected image number it came from. // Empty means identity (image i == original image i). Lets a dataset be a subset (or strided // selection) of the truly collected images: reprocessing snapshots over a sub-range, and (in diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index d9328a29..6b22bb74 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -12,6 +12,8 @@ ADD_LIBRARY(Rugnux STATIC RugnuxCommandLine.h ModelValidation.cpp ModelValidation.h + ResultReport.cpp + ResultReport.h ) TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi) diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp new file mode 100644 index 00000000..c541e300 --- /dev/null +++ b/rugnux/ResultReport.cpp @@ -0,0 +1,283 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include + +#include "../common/GitInfo.h" +#include "../common/time_utc.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" + +#include "ResultReport.h" + +namespace { + // The version of this file format. Bumped when a key is renamed or removed, a table column moves, + // or a reason code changes meaning - a consumer can gate on it. + constexpr int REPORT_VERSION = 1; + + const char *BANNER = " ******************************************************************************"; + + void Section(std::ostream &os, const std::string &title) { + os << "\n" << BANNER << "\n " << title << "\n" << BANNER << "\n\n"; + } + + // Every number a consumer might want is written as one of these, so it is one grep away. + template void Key(std::ostream &os, const char *key, const T &value) { + os << key << "= " << value << "\n"; + } + + std::string CellString(const UnitCell &c) { + return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c, + c.alpha, c.beta, c.gamma); + } +} + +std::string RenderResultReport(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result) { + std::ostringstream os; + const bool rotation = experiment.IsRotationIndexing(); + const bool merged = result.has_merge_statistics; + std::vector warnings = result.warnings; + + os << BANNER << "\n" + << " RUGNUX PROCESSING REPORT\n" + << BANNER << "\n\n" + << " What this run determined, written next to its other output. The `KEY= value` lines and\n" + << " the tables below are a stable interface - a script greps them, and REPORT_VERSION says\n" + << " when that interface last changed. Timing, rates and per-image progress are not here;\n" + << " they are on stdout.\n\n"; + + Key(os, "REPORT_VERSION", REPORT_VERSION); + Key(os, "RUGNUX_VERSION", jfjoch_version()); + if (!jfjoch_git_sha1().empty()) + Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date()); + Key(os, "DATE", time_UTC(std::chrono::system_clock::now())); + Key(os, "INPUT_FILE", input_file); + Key(os, "OUTPUT_PREFIX", output_prefix); + + // ---------------------------------------------------------------- 1. DATA SET + Section(os, "1. DATA SET"); + Key(os, "EXPERIMENT_TYPE", rotation ? "ROTATION" : "STILLS"); + Key(os, "IMAGES_PROCESSED", result.images_processed); + Key(os, "WAVELENGTH", fmt::format("{:.5f}", experiment.GetWavelength_A())); + if (const auto gonio = experiment.GetGoniometer()) { + Key(os, "OSCILLATION_RANGE", fmt::format("{:.4f}", gonio->GetIncrement_deg())); + Key(os, "STARTING_ANGLE", fmt::format("{:.3f}", gonio->GetStart_deg())); + const auto ax = gonio->GetAxis(); + Key(os, "ROTATION_AXIS", fmt::format("{:.6f} {:.6f} {:.6f}", ax.x, ax.y, ax.z)); + } + Key(os, "DETECTOR_DISTANCE", fmt::format("{:.3f}", result.used_distance_mm)); + Key(os, "BEAM_CENTRE", fmt::format("{:.2f} {:.2f}", result.used_beam_x_pxl, result.used_beam_y_pxl)); + os << "\n" + << " The distance and beam centre above are the ones this result was integrated at, which on\n" + << " a rotation run is the post-refined geometry rather than the values in the input file.\n"; + + if (result.pass_count > 1) { + os << "\n"; + Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count)); + Key(os, "PASS_DECISION", result.pass_decision); + os << "\n" + << " A rotation run integrates twice: once at the geometry in the input file, then again at\n" + << " the post-refined geometry. Every number in this report describes the pass named above,\n" + << " whose files are " << output_prefix << ".*; the header-geometry pass was written to\n" + << " " << output_prefix << "_01.* and is kept only for comparison.\n"; + } + + // ---------------------------------------------------------------- 2. INDEXING + Section(os, "2. INDEXING"); + if (result.indexing_rate.has_value()) + Key(os, "INDEXING_RATE", fmt::format("{:.4f}", result.indexing_rate.value())); + Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE")); + if (result.consensus_cell.has_value()) + Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell)); + if (result.space_group_number.has_value()) + Key(os, "SPACE_GROUP_NUMBER", result.space_group_number.value()); + if (result.indexing_rate.value_or(0.0f) <= 0.0f) + warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset"); + + // ---------------------------------------------- 3. GEOMETRY POST-REFINEMENT + if (result.post_refine.has_value()) { + const auto &pr = *result.post_refine; + Section(os, "3. GEOMETRY POST-REFINEMENT"); + os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n" + << " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n" + << " step is committed only if it improves a held-out residual.\n\n"; + Key(os, "POSTREFINE_EVENTS_USED", pr.events_used); + Key(os, "POSTREFINE_OBS_USED", pr.obs_used); + Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE"); + Key(os, "POSTREFINE_DETECTOR_COMMITTED", pr.detector_refined ? "TRUE" : "FALSE"); + Key(os, "POSTREFINE_DISTANCE", fmt::format("{:.3f} -> {:.3f}", pr.distance_before_mm, + pr.distance_after_mm)); + Key(os, "POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}", + pr.beam_x_before_px, pr.beam_y_before_px, + pr.beam_x_after_px, pr.beam_y_after_px)); + Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale)); + os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n" + << " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It is\n" + << " reported only; nothing here corrects for it.\n"; + if (pr.rotation_scale_suspect) + warnings.emplace_back(fmt::format( + "The goniometer turned by a factor {:.5f} of the angles stored in the file - the " + "stage rotation looks mis-calibrated by {:+.2f}%", + pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0))); + } + + // ---------------------------------------------- 4. SPACE GROUP DETERMINATION + Section(os, "4. SPACE GROUP DETERMINATION"); + if (result.space_group_search.has_value()) { + Key(os, "SPACE_GROUP_SEARCH", "DE_NOVO"); + os << "\n" << SearchSpaceGroupResultToText(*result.space_group_search) << "\n"; + } else if (result.space_group_number.has_value()) { + Key(os, "SPACE_GROUP_SEARCH", "FIXED"); + os << "\n The space group was given, not determined here.\n"; + } else { + Key(os, "SPACE_GROUP_SEARCH", "NONE"); + os << "\n No space group was determined.\n"; + } + + // ---------------------------------------------------- 5. SCALING AND MERGING + Section(os, "5. SCALING AND MERGING"); + if (!merged) { + Key(os, "MERGE", "NOT_PERFORMED"); + os << "\n No scaling or merging was performed on this run, so there are no merging statistics, no\n" + << " error model, and no sweep-quality diagnosis below. The integrated reflections are in\n" + << " " << output_prefix << "_process.h5.\n"; + } else { + const auto &o = result.merge_statistics.overall; + Key(os, "MERGE", "PERFORMED"); + Key(os, "INCLUDE_RESOLUTION_RANGE", fmt::format("{:.3f} {:.3f}", o.d_max, o.d_min)); + Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE"); + Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections); + Key(os, "TOTAL_OBSERVATIONS", o.total_observations); + Key(os, "COMPLETENESS", o.possible_unique_reflections > 0 + ? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections) + : std::string("nan")); + Key(os, "MULTIPLICITY", o.unique_reflections > 0 + ? fmt::format("{:.2f}", static_cast(o.total_observations) / o.unique_reflections) + : std::string("nan")); + Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma)); + Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas)); + Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half)); + Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous)); + Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b)); + // The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP. + Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a)); + Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b)); + Key(os, "ISA", fmt::format("{:.2f}", result.error_model_isa)); + if (result.error_model_isa_asymptotic > 0.0) + Key(os, "ISA_ASYMPTOTIC", fmt::format("{:.2f}", result.error_model_isa_asymptotic)); + Key(os, "REFERENCE_DATA_USED", result.has_reference ? "TRUE" : "FALSE"); + // The shell table straight off the statistics rather than result.merge_statistics_text: that + // string also carries the twinning analysis and the advisories, which have sections of their own. + os << "\n ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n" + << " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n" + << " is the strong-reflection tier only.\n\n" + << result.merge_statistics; + } + + // --------------------------------------------------------------- 6. TWINNING + if (merged && result.twinning.l_test_pairs > 0) { + Section(os, "6. TWINNING"); + Key(os, "TWINNING_SUSPECTED", result.twinning.twinning_suspected ? "TRUE" : "FALSE"); + Key(os, "L_TEST_MEAN_ABS_L", fmt::format("{:.4f}", result.twinning.mean_abs_l)); + Key(os, "L_TEST_MEAN_L_SQUARED", fmt::format("{:.4f}", result.twinning.mean_l_squared)); + Key(os, "SECOND_MOMENT_I", fmt::format("{:.4f}", result.twinning.second_moment)); + Key(os, "ESTIMATED_TWIN_FRACTION", fmt::format("{:.3f}", result.twinning.estimated_twin_fraction)); + os << "\n" << TwinningAnalysisToText(result.twinning) << "\n"; + if (result.twinning.twinning_suspected) + warnings.emplace_back(fmt::format( + "Twinning is indicated (<|L|> = {:.3f}, /^2 = {:.3f}, estimated twin " + "fraction {:.2f}) - refine against the merged data with care", + result.twinning.mean_abs_l, result.twinning.second_moment, + result.twinning.estimated_twin_fraction)); + } + + // ------------------------------------------------------- 7. RADIATION DAMAGE + if (!result.radiation_damage_text.empty()) { + Section(os, "7. RADIATION DAMAGE"); + Key(os, "RADIATION_DAMAGE_RELATIVE_B", fmt::format("{:.2f}", + result.merge_statistics.radiation_damage_delta_b)); + os << "\n" << result.radiation_damage_text << "\n"; + } + + // ------------------------------------------------------ 8. SWEEP QUALITY + const auto &sq = result.merge_statistics.sweep_quality; + Section(os, "8. SWEEP QUALITY"); + os << " Stretches of the sweep over which the crystal delivered much less than the rest of the run.\n" + << " REASON comes from a closed vocabulary, listed below so a consumer can tell an unknown code\n" + << " from a missing one. SEVERITY is the fraction of the run's typical diffracting power missing\n" + << " over the range (0 = as good as the run, 1 = nothing at all); SCALE and CC are the range's\n" + << " mean per-image scale and CC-to-merge relative to the run median; INDEXED is the fraction of\n" + << " the range's frames that were scaled at all. Nothing is excluded on the strength of this.\n\n"; + Key(os, "SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED"); + Key(os, "SWEEP_QUALITY_COUNT", sq.ranges.size()); + { + std::string codes; + for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) + codes += (codes.empty() ? "" : " ") + + std::string(SweepQualityReasonCode(static_cast(r))); + Key(os, "SWEEP_QUALITY_REASONS", codes); + } + if (sq.measured) { + Key(os, "SWEEP_ROTATION", fmt::format("{:.1f}", sq.sweep_deg)); + Key(os, "FLUX_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.flux_peak_to_trough)); + Key(os, "SCALE_MODULATION_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.modulation_peak_to_trough)); + } + os << "\n" + << " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED\n" + << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n"; + for (const auto &r : sq.ranges) { + os << fmt::format(" {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f}\n", + r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg, + SweepQualityReasonCode(r.reason), r.severity, r.mean_relative_scale, + r.mean_relative_cc, r.indexed_fraction); + warnings.push_back(fmt::format( + "Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)", + r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg, + r.mean_relative_scale, r.mean_relative_cc, 100.0 * r.indexed_fraction)); + } + os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n"; + + // --------------------------------------------------------------- 9. WARNINGS + if (result.cancelled) + warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report " + "describes an incomplete run", result.images_processed)); + Section(os, "9. WARNINGS"); + os << " Everything that needs a person's attention, one line each, marked so a script can find\n" + << " them with a single grep for \"WARNING:\".\n\n"; + Key(os, "WARNING_COUNT", warnings.size()); + os << "\n"; + for (const auto &w : warnings) + os << "WARNING: " << w << "\n"; + if (warnings.empty()) + os << " (none)\n"; + + os << "\n" << BANNER << "\n END OF REPORT\n" << BANNER << "\n"; + return os.str(); +} + +void WriteResultReport(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + Logger &logger) { + if (output_prefix.empty()) + return; // "compute the statistics, persist nothing" + + const std::string filename = output_prefix + "_report.txt"; + // The report is unconditional, so it must never be the reason a run fails: a run that produced a + // good .mtz must survive an unwritable path or a full disk. Report the failure and carry on. + try { + std::ofstream file(filename); + file.exceptions(std::ios::failbit | std::ios::badbit); + file << RenderResultReport(output_prefix, input_file, experiment, result); + } catch (const std::exception &e) { + logger.Warning("Could not write the results report {}: {}", filename, e.what()); + } +} diff --git a/rugnux/ResultReport.h b/rugnux/ResultReport.h new file mode 100644 index 00000000..be4424c6 --- /dev/null +++ b/rugnux/ResultReport.h @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../common/DiffractionExperiment.h" +#include "../common/Logger.h" +#include "Rugnux.h" + +// _report.txt - what the run DETERMINED, written next to the .mtz/.cif/.hkl on every run +// that has an output prefix. Modelled on XDS's CORRECT.LP: prose and tables a crystallographer +// reads top to bottom, with `KEY= value` assignment lines and fixed-width tables a script greps +// without parsing prose. Timing, rates and progress are deliberately absent - those are stdout. +// +// The format is an interface: key names, table columns and the reason vocabulary are stable, and +// REPORT_VERSION is bumped if they ever change. See docs/RUGNUX.md. + +// Renders the report. Exposed for testing; RunPipeline results are the only input. +std::string RenderResultReport(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result); + +// Renders and writes _report.txt. Does nothing when the prefix is empty (the +// "compute statistics, persist nothing" mode). Never throws: a run that produced good reflections +// must not be lost because a side file could not be written, so a failure is logged and swallowed. +void WriteResultReport(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + Logger &logger); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index c24ec72c..0c7599cf 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -332,6 +332,8 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted. config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01"; auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true); + pass1.pass_number = 1; + pass1.pass_count = 2; if (cancelled_) { config_.output_prefix = base_prefix; return pass1; } if (gonio_snapshot) experiment_.Goniometer(*gonio_snapshot); // undo the pre-pass goniometer shift @@ -357,6 +359,11 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { base_prefix); config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix) auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + // The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever + // result is returned - it is a result of the run, not of the pass that happened to fit it. + pass2.post_refine = pass1.post_refine; + pass2.pass_number = 2; + pass2.pass_count = 2; // Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a // periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell @@ -377,6 +384,13 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { v2, v2 / v1, v1); force_rotation_result_ = *prepass_result_; pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); + pass2.post_refine = pass1.post_refine; + pass2.pass_number = 3; + pass2.pass_count = 3; + pass2.pass_decision = fmt::format( + "post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo " + "cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious supercell", + v2, v2 / v1, v1); force_rotation_result_.reset(); } } @@ -420,9 +434,24 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); if (!redo.space_group_search.has_value()) redo.space_group_search = pass1.space_group_search; + redo.post_refine = pass1.post_refine; + redo.pass_number = pass2.pass_count + 1; + redo.pass_count = pass2.pass_count + 1; + redo.pass_decision = fmt::format( + "header geometry re-adopted: the post-refined pass was worse " + "(completeness {:.1f}% vs {:.1f}%, CC1/2 {:.3f} vs {:.3f})", + compl2, completeness(pass1), cc2, cc1); return redo; } + if (pass2.pass_decision.empty()) + pass2.pass_decision = fmt::format( + "post-refined geometry adopted (completeness {:.1f}% vs {:.1f}%, " + "CC1/2 {:.3f} vs {:.3f})", compl2, completeness(pass1), cc2, cc1); } + if (pass2.pass_decision.empty()) + pass2.pass_decision = prepass_detector_geometry_ + ? "post-refined geometry adopted" + : "the post-refinement committed no geometry change, so this pass reproduces the first"; return pass2; } return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); @@ -1498,6 +1527,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b prs.refine_geometry = true; const auto pr = PostRefineRotationGeometry(outcomes, *rot->axis, experiment_.GetDiffractionGeometry(), rot->lattice, prs, logger); + result.post_refine = pr; if (pr.ok && pr.detector_refined) { // Store, don't apply: this pass continues to write its own (header-geometry) output; // Run() applies the refined geometry before the second pass. @@ -1645,6 +1675,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b a, b, alt.point_group_hm, sg_search.point_group_hm, a, b); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; + result.warnings.push_back(msg); } if (order_of(sg_search) <= order_of(alt)) { if (order_of(alt) > order_of(sg_search)) @@ -1891,6 +1922,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b order, angle); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; + result.warnings.push_back(msg); } else { stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " << std::fixed << std::setprecision(1) << angle << " deg\n\n"; @@ -1925,6 +1957,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals"); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; + result.warnings.push_back(msg); } else { logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved against the " "supplied reference.", laws); @@ -1947,6 +1980,20 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.merge_statistics_text = stats_text.str(); result.has_merge_statistics = true; result.merge_statistics = sm.statistics; + + // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where + // the image is in no flagged range, plus the vocabulary the codes index. Only filled when the + // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean". + if (sm.statistics.sweep_quality.measured) { + end_msg.sweep_quality.assign(end_msg.max_image_number, 0); + for (const auto &r : sm.statistics.sweep_quality.ranges) + for (int64_t i = std::max(0, r.first_image); + i <= r.last_image && i < static_cast(end_msg.sweep_quality.size()); ++i) + end_msg.sweep_quality[i] = static_cast(r.reason) + 1; + for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) + end_msg.sweep_quality_reasons.emplace_back( + SweepQualityReasonCode(static_cast(r))); + } { // Stride rather than take the head: the merged list is ordered by hkl, so the first N // reflections are one corner of reciprocal space and would not show the intensity range. @@ -2023,6 +2070,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b os << fmt::format(" {:.1f}", bb); } logger.Info("{}", os.str()); + result.radiation_damage_text = os.str(); } } @@ -2087,5 +2135,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("{} {} images in {:.2f} s ({:.2f} Hz)", cancelled_ ? "Cancelled after" : "Processed", result.images_processed, result.processing_time_s, result.frame_rate_hz); + + // The geometry this pass integrated at, for the report. Taken from experiment_ rather than from the + // caller's copy: the second pass of the rotation two-pass runs at the post-refined geometry. + result.used_beam_x_pxl = experiment_.GetBeamX_pxl(); + result.used_beam_y_pxl = experiment_.GetBeamY_pxl(); + result.used_distance_mm = experiment_.GetDetectorDistance_mm(); return result; } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index cb5291ec..3ea8cd43 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -19,6 +19,7 @@ #include "../image_analysis/scale_merge/Merge.h" // MergeStatistics #include "../image_analysis/scale_merge/TwinningAnalysis.h" // TwinningAnalysisResult #include "../image_analysis/scale_merge/SearchSpaceGroup.h" // SearchSpaceGroupResult +#include "../image_analysis/geom_refinement/PostRefine.h" // PostRefineResult #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult #include "RugnuxCalibration.h" // CalibrationMethod, CalibrationResult @@ -158,6 +159,31 @@ struct ProcessResult { // populated only when a search ran (no space group was fixed). The library no longer renders it to // text - the CLI formats it for stdout, the viewer draws it as a table. std::optional space_group_search; + + // ---- Material for the results report (_report.txt, see rugnux/ResultReport.h) ---- + + // The detector geometry this pass integrated at - not necessarily the header geometry, since the + // rotation two-pass post-refines it between the passes. + float used_beam_x_pxl = NAN; + float used_beam_y_pxl = NAN; + float used_distance_mm = NAN; + + // Rotation two-pass geometry post-refinement, as measured by the first pass. Empty on a run that + // did not post-refine (stills, single pass, or the refinement did not converge). + std::optional post_refine; + + // Which pass this result is (1-based) and how many ran, plus one line saying which pass became the + // canonical output and on what evidence. pass_count == 1 outside the rotation two-pass. + int pass_number = 1; + int pass_count = 1; + std::string pass_decision; + + // The radiation-damage report (rotation), the same text the log prints. Empty when not measured. + std::string radiation_damage_text; + + // Conditions that need a person's attention, one plain sentence each (an ambiguous space group, a + // symmetry axis on the spindle, an indexing ambiguity). The same messages the log warns about. + std::vector warnings; }; // Callbacks for progress and live results. Methods may be called from worker threads, so an diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index d7161066..f2b8022f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -38,6 +38,7 @@ #include "Rugnux.h" #include "RugnuxDefaults.h" #include "ModelValidation.h" +#include "ResultReport.h" // Spots kept per image (the strongest ones) and handed to indexing. Offline reprocessing is not // bound by the online spot budget, so this is rugnux's own default rather than the 250 the @@ -72,7 +73,8 @@ constexpr int CALIBRATION_AZIM_BINS_DEFAULT = 32; void print_usage() { std::cout << "Usage rugnux {} " << std::endl; std::cout << "Options:" << std::endl; - std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; + std::cout << " -o, --output-prefix Output file prefix (default: output). mx and scale runs always write " + "_report.txt, the results report; an empty prefix writes nothing at all" << std::endl; std::cout << " -N, --threads Number of threads (default: all hardware threads)" << std::endl; std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; std::cout << " -e, --end-image End image number (default: all)" << std::endl; @@ -1499,6 +1501,26 @@ static int RunRugnux(int argc, char **argv) { data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, /*probe_indexing_ambiguity=*/reference_data.empty()); } + + // --mode scale re-merges stored reflections, so it determines a merging result and gets the + // same report. It has no ProcessResult of its own; fill the parts that exist here. + ProcessResult scale_result; + scale_result.images_processed = reader.GetNumberOfImages(); + scale_result.consensus_cell = experiment.GetUnitCell(); + scale_result.space_group_number = experiment.GetSpaceGroupNumber(); + scale_result.used_beam_x_pxl = experiment.GetBeamX_pxl(); + scale_result.used_beam_y_pxl = experiment.GetBeamY_pxl(); + scale_result.used_distance_mm = experiment.GetDetectorDistance_mm(); + scale_result.has_merge_statistics = true; + scale_result.merge_statistics = merged_statistics; + { std::ostringstream s; s << merged_statistics; scale_result.merge_statistics_text = s.str(); } + scale_result.error_model_isa = error_model_isa; + scale_result.error_model_isa_asymptotic = error_model_isa_asymptotic; + scale_result.error_model_a = error_model_a; + scale_result.error_model_b = error_model_b; + scale_result.has_reference = !reference_data.empty(); + scale_result.twinning = twinning; + WriteResultReport(output_prefix, input_file, experiment, scale_result, logger); return 0; } @@ -2035,6 +2057,10 @@ static int RunRugnux(int argc, char **argv) { } g_active_process = nullptr; + // The results report, next to the .mtz/.cif/.hkl. Written on every run with an output prefix - + // including --no-merge, which still determined an indexing and geometry result worth recording. + WriteResultReport(output_prefix, input_file, experiment, result, logger); + // The space-group search is rendered here (not in the library) so the viewer does not emit it on // stdout and the CLI owns the format. if (result.space_group_search.has_value()) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 07f004fd..65eb3eb4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -51,6 +51,7 @@ ADD_EXECUTABLE(jfjoch_test ZMQMetadataSocketTest.cpp JFJochReaderTest.cpp RugnuxTest.cpp + ResultReportTest.cpp RugnuxLargeTest.cpp TestData.h MovingAverageTest.cpp diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp new file mode 100644 index 00000000..9462eb5a --- /dev/null +++ b/tests/ResultReportTest.cpp @@ -0,0 +1,160 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../common/DiffractionExperiment.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../reader/JFJochHDF5Reader.h" +#include "../rugnux/ResultReport.h" +#include "../writer/FileWriter.h" + +// Two synthetic ranges, one of each shape the report has to handle. +namespace { + SweepQuality TestSweepQuality() { + SweepQuality out; + out.measured = true; + out.sweep_deg = 60.0f; + out.ranges.push_back(SweepQualityRange{.first_image = 100, .last_image = 149, + .reason = SweepQualityReason::CrystalOutOfBeam, + .severity = 0.83f, .rotation_deg = 5.0f, + .mean_relative_scale = 0.12f, .mean_relative_cc = 0.30f, + .indexed_fraction = 0.015f}); + out.ranges.push_back(SweepQualityRange{.first_image = 400, .last_image = 499, + .reason = SweepQualityReason::RadiationDamage, + .severity = 0.41f, .rotation_deg = 10.0f, + .mean_relative_scale = 0.55f, .mean_relative_cc = 0.80f, + .indexed_fraction = 0.62f}); + return out; + } + + std::vector ReasonVocabulary() { + std::vector out; + for (int r = 0; r <= static_cast(SweepQualityReason::RadiationDamage); ++r) + out.emplace_back(SweepQualityReasonCode(static_cast(r))); + return out; + } +} + +TEST_CASE("SweepQuality_ReasonVocabulary", "[Diagnostics]") { + // The codes are an interface - they are written verbatim into _report.txt and into HDF5, + // where the per-image codes index this list from 1. A rename or a reorder silently breaks every + // consumer, so pin both the spelling and the order here. + const auto codes = ReasonVocabulary(); + REQUIRE(codes.size() == 5); + CHECK(codes[0] == "no_diffraction"); + CHECK(codes[1] == "crystal_out_of_beam"); + CHECK(codes[2] == "weak_diffraction"); + CHECK(codes[3] == "loss_of_centring"); + CHECK(codes[4] == "radiation_damage"); +} + +TEST_CASE("ResultReport_Render", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + x.ImagesPerTrigger(600); + + ProcessResult result; + result.images_processed = 600; + result.indexing_rate = 0.87f; + result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; + result.space_group_number = 96; + result.used_beam_x_pxl = 766.62f; + result.used_beam_y_pxl = 846.87f; + result.used_distance_mm = 243.53f; + result.has_merge_statistics = true; + result.merge_statistics.sweep_quality = TestSweepQuality(); + + const auto text = RenderResultReport("prefix", "in.h5", x, result); + + // The stable keys a consumer greps for. + CHECK(text.find("\nREPORT_VERSION= 1\n") != std::string::npos); + CHECK(text.find("\nOUTPUT_PREFIX= prefix\n") != std::string::npos); + CHECK(text.find("\nINDEXING_RATE= 0.8700\n") != std::string::npos); + CHECK(text.find("\nSPACE_GROUP_NUMBER= 96\n") != std::string::npos); + CHECK(text.find("\nSWEEP_QUALITY_STATUS= COMPUTED\n") != std::string::npos); + CHECK(text.find("\nSWEEP_QUALITY_COUNT= 2\n") != std::string::npos); + CHECK(text.find("\nSWEEP_QUALITY_REASONS= no_diffraction crystal_out_of_beam weak_diffraction " + "loss_of_centring radiation_damage\n") != std::string::npos); + + // One table row per range, with the reason code verbatim. + CHECK(text.find(" 100 149 50 5.0 crystal_out_of_beam ") + != std::string::npos); + CHECK(text.find(" 400 499 100 10.0 radiation_damage ") + != std::string::npos); + + // ... and one plain-English WARNING line per range, greppable by the marker alone. + CHECK(text.find("\nWARNING_COUNT= 2\n") != std::string::npos); + CHECK(text.find("\nWARNING: Frames 100-149 out of beam (5.0 deg,") != std::string::npos); + CHECK(text.find("\nWARNING: Frames 400-499 radiation damage (10.0 deg,") != std::string::npos); +} + +TEST_CASE("ResultReport_RenderEmpty", "[Diagnostics]") { + // A clean run and a run that never looked must be distinguishable: both have a count of 0, and + // only the STATUS key separates them. This is the property a consumer relies on. + DiffractionExperiment x(DetJF(1)); + + ProcessResult clean; + clean.has_merge_statistics = true; + clean.merge_statistics.sweep_quality.measured = true; + const auto clean_text = RenderResultReport("p", "in.h5", x, clean); + CHECK(clean_text.find("\nSWEEP_QUALITY_STATUS= COMPUTED\n") != std::string::npos); + CHECK(clean_text.find("\nSWEEP_QUALITY_COUNT= 0\n") != std::string::npos); + CHECK(clean_text.find("FIRST_IMAGE LAST_IMAGE") != std::string::npos); + + ProcessResult not_merged; + const auto not_merged_text = RenderResultReport("p", "in.h5", x, not_merged); + CHECK(not_merged_text.find("\nSWEEP_QUALITY_STATUS= NOT_COMPUTED\n") != std::string::npos); + CHECK(not_merged_text.find("\nMERGE= NOT_PERFORMED\n") != std::string::npos); + CHECK(not_merged_text.find("FIRST_IMAGE LAST_IMAGE") != std::string::npos); +} + +TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") { + // The per-image codes have to survive the writer and come back out of the reader. Without an + // assertion here the field can ship as all-zeros without anyone noticing. + DiffractionExperiment x(DetJF(1)); + x.ImagesPerTrigger(6).Compression(CompressionAlgorithm::NO_COMPRESSION) + .FilePrefix("sweep_quality_roundtrip"); + x.SetFileWriterFormat(FileWriterFormat::NXmxIntegrated).OverwriteExistingFiles(true); + + // Images 2-3 out of beam, image 5 dead from radiation damage; the rest in no flagged range. + const std::vector expected{0, 0, + static_cast(SweepQualityReason::CrystalOutOfBeam) + 1, + static_cast(SweepQualityReason::CrystalOutOfBeam) + 1, + 0, + static_cast(SweepQualityReason::RadiationDamage) + 1}; + + { + RegisterHDF5Filter(); + StartMessage start_message; + x.FillMessage(start_message); + + EndMessage end_message; + end_message.max_image_number = x.GetImageNum(); + end_message.sweep_quality = expected; + end_message.sweep_quality_reasons = ReasonVocabulary(); + + FileWriter writer(start_message); + std::vector image(x.GetPixelsNum(), 42); + for (int i = 0; i < x.GetImageNum(); i++) { + DataMessage message{}; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + REQUIRE_NOTHROW(writer.Write(message)); + } + writer.WriteHDF5(end_message); + writer.Finalize(); + } + + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("sweep_quality_roundtrip_master.h5")); + auto dataset = reader.GetDataset(); + CHECK(dataset->sweep_quality == expected); + // The vocabulary travels with the codes, so a consumer can name them without this source. + CHECK(dataset->sweep_quality_reasons == ReasonVocabulary()); + } + + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); + remove("sweep_quality_roundtrip_master.h5"); +} diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index fd779752..5ab2a398 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1050,4 +1050,12 @@ void NXmx::EndResultVectors(const EndMessage &end) { SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageScaleMosaicity", end.image_scale_mosaicity, "deg"); if (!end.niggli_class.empty()) SaveVectorIfMissing(*hdf5_file, "/entry/MX/niggliClass", end.niggli_class); + // Per-image sweep-quality code, with the vocabulary next to it so the codes can be read without + // this source: sweepQuality[i] == 0 means the image is in no flagged range, otherwise it indexes + // sweepQualityReasons from 1. Absent when the diagnostic did not run. + if (!end.sweep_quality.empty() && !end.sweep_quality_reasons.empty()) { + SaveVectorIfMissing(*hdf5_file, "/entry/MX/sweepQuality", end.sweep_quality); + if (!hdf5_file->Exists("/entry/MX/sweepQualityReasons")) + hdf5_file->SaveVector("/entry/MX/sweepQualityReasons", end.sweep_quality_reasons); + } }