From 5f731db81dbf514c9ee6dbb44c33b49eb317d302 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 17 Oct 2026 16:51:20 +0200 Subject: [PATCH] rugnux: the library returns its outputs in memory; the callers write the files The rugnux library wrote its results itself: _process.h5 through a FileWriter, the merged/twin/_mr reflection files, the unmerged MTZs, the maps and their MTZ, the placed models, _detector.jpg and _plot.txt, and on a doubled cell it renamed files on disk. It now writes none of them. - MtzTable (common/MtzTable.h): an MTZ file's header, columns and batch headers with no file behind it. It replaces MergedReflectionsMessage and MapCoefficientsMessage in the END message, with one CBOR encoding. WriteMtz (writer/WriteRugnuxFiles.h) is the one serialiser, used by rugnux and jfjoch_writer alike; it refuses a table whose columns differ in length, name a missing dataset, carry a malformed batch header or name a space group gemmi does not know. MergedMtzTable and UnmergedMtzTable build the tables the old writers built. - ProcessResult::outputs (RugnuxOutputs) holds, for the output pass only: the answer's merge (MTZ table plus what the .cif/.hkl are written from), the twin hypothesis' and the model's choice's merges, the unmerged and unmerged-partials tables, the maps (coefficients always; the grids only with ProcessConfig::keep_map_grids), the placed models, the detector JPEG and the per-image plot table, and the start message of the no-lattice _process.h5. Map RMS is taken from the grid in memory, as the CCP4 header stated it, whether or not the map is kept. - The per-image analysis goes to a caller-supplied RugnuxImageSink; rugnux and the viewer give a FileWriter-based ProcessFileSink. The geometry pre-pass no longer makes a _01_process.h5. - The doubled cell keeps both output sets: the answer's in outputs, the other's in DoubledCell::other_outputs with its suffix; the report names files by suffix after the caller's prefix. - The output prefix, write_process_h5 and the observation dump path are gone from ProcessConfig; --dump-observations is a callback. - What used to be gated on having files to write (model validation and MR, the twin-hypothesis merges, the finalist ledger, the scaled fulls, the unmerged tables) is gated on the output pass instead, so a caller with no prefix gets the same results. - RugnuxOutputFiles (new library JFJochRugnuxOutput, linked by rugnux, the viewer and the tests) writes every file from the result with the names and formats rugnux has always used. The viewer writes its job's files the same way. Co-Authored-By: Claude Opus 5.5 --- common/CMakeLists.txt | 1 + common/JFJochMessages.h | 34 +--- common/MtzTable.h | 48 +++++ frame_serialize/CBORStream2Deserializer.cpp | 108 ++++++----- frame_serialize/CBORStream2Serializer.cpp | 87 +++++---- image_analysis/WriteReflections.cpp | 167 +++++++++-------- image_analysis/WriteReflections.h | 40 ++-- .../scale_merge/RotationScaleMerge.cpp | 23 +-- .../scale_merge/RotationScaleMerge.h | 5 +- .../structure_refinement/ModelValidation.cpp | 163 +++++++--------- .../structure_refinement/ModelValidation.h | 29 ++- receiver/JFJochReceiverRugnux.cpp | 4 +- rugnux/CMakeLists.txt | 9 +- rugnux/DiagnosticOutput.cpp | 18 +- rugnux/DiagnosticOutput.h | 9 +- rugnux/ResultReport.cpp | 45 ++--- rugnux/ResultReport.h | 11 -- rugnux/Rugnux.cpp | 122 +++++------- rugnux/Rugnux.h | 129 +++++++++---- rugnux/RugnuxCalibration.cpp | 3 - rugnux/RugnuxCommandLine.cpp | 8 +- rugnux/RugnuxCommandLine.h | 4 + rugnux/RugnuxFirstPass.cpp | 3 - rugnux/RugnuxFirstPassAccept.cpp | 3 - rugnux/RugnuxFirstPassRescues.cpp | 3 - rugnux/RugnuxImagePass.cpp | 5 +- rugnux/RugnuxLive.cpp | 4 +- rugnux/RugnuxLive.h | 8 +- rugnux/RugnuxOutputFiles.cpp | 176 ++++++++++++++++++ rugnux/RugnuxOutputFiles.h | 58 ++++++ rugnux/RugnuxPasses.cpp | 64 ++----- rugnux/RugnuxPipeline.h | 23 +-- rugnux/RugnuxPreScan.cpp | 3 - rugnux/RugnuxScaleMerge.cpp | 165 ++++++++-------- rugnux/RugnuxStills.cpp | 3 - rugnux/rugnux_cli.cpp | 45 ++++- tests/CMakeLists.txt | 2 +- tests/DiagnosticOutputTest.cpp | 11 +- tests/ModelValidationTest.cpp | 41 ++-- tests/RugnuxResultFilesTest.cpp | 77 +++++--- tests/RugnuxTest.cpp | 25 ++- tests/WriteReflectionsTest.cpp | 29 +-- viewer/CMakeLists.txt | 2 +- viewer/JFJochProcessController.cpp | 23 ++- viewer/JFJochProcessController.h | 7 +- viewer/windows/JFJochProcessingJobsWindow.cpp | 23 +-- writer/RugnuxResultBlocks.cpp | 97 +++------- writer/RugnuxResultBlocks.h | 3 +- writer/WriteRugnuxFiles.cpp | 146 +++++---------- writer/WriteRugnuxFiles.h | 12 +- 50 files changed, 1163 insertions(+), 965 deletions(-) create mode 100644 common/MtzTable.h create mode 100644 rugnux/RugnuxOutputFiles.cpp create mode 100644 rugnux/RugnuxOutputFiles.h diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index f8a8ff5c9..e447b5509 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -161,6 +161,7 @@ ADD_LIBRARY(JFJochCommon STATIC ScalingSettings.h UnitCell.cpp UnitCell.h + MtzTable.h CorrelationCoefficient.cpp CorrelationCoefficient.h JFJochReceiverPlots.cpp diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index bcb4dcd1e..449929d2f 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -16,6 +16,7 @@ #include "CompressedImage.h" #include "SpotToSave.h" #include "UnitCell.h" +#include "MtzTable.h" #include "GoniometerAxis.h" #include "DetectorTransformation.h" #include "GridScanSettings.h" @@ -368,35 +369,6 @@ struct StartMessage { std::optional smargon_position; }; -// The merged reflections as rugnux writes them to .mtz (WriteMergedMtz, -// writer/WriteRugnuxFiles.h): one row per unique reflection in the CCP4 anomalous layout, one vector -// per MTZ column. The anomalous columns (i_plus ... sigma_f_minus) are all empty where the merge kept -// no Bijvoet split, and the file then has none. free_flag is the FreeR_flag value, 0 for the test set. -struct MergedReflectionsMessage { - UnitCell unit_cell; - std::string space_group_name; // gemmi xhm(): the name keeps the setting, a number does not - float wavelength_A = 0; - std::string crystal_name; - std::vector history; // the MTZ HISTORY records after the one naming the program - std::vector h, k, l; - std::vector i_mean, sigma_i_mean; - std::vector i_plus, sigma_i_plus, i_minus, sigma_i_minus; - std::vector f, sigma_f; - std::vector f_plus, sigma_f_plus, f_minus, sigma_f_minus; - std::vector free_flag; -}; - -// The sigma_A-weighted map coefficients rugnux writes to _maps.mtz (WriteMapCoefficientsMtz), -// one vector per MTZ column. Phases in degrees, [0, 360). free_flag is FREE: 0 for the test set. -struct MapCoefficientsMessage { - UnitCell unit_cell; - std::string space_group_name; - float wavelength_A = 0; - std::vector h, k, l; - std::vector fp, fc, phic, fwt, phwt, delfwt, phdelwt, fom; - std::vector free_flag; -}; - // The atomic model as placed against the data, one vector per atom field in file order, for // WritePlacedModel to write _model.cif. A new residue starts wherever the model, chain or // residue key changes. unit_cell and space_group_name are what WritePlacedModel is given (the data's); @@ -537,8 +509,8 @@ struct EndMessage { // The merged result of an analysis run inside the broker, for the writer to put beside the master // file as .mtz, _maps.mtz and _model.cif. Absent where nothing was merged. - std::optional merged; - std::optional map_coefficients; + std::optional merged; // MergedMtzTable (image_analysis/WriteReflections.h) + std::optional map_coefficients; // ModelMaps::coefficients (ModelValidation.h) std::optional placed_model; }; diff --git a/common/MtzTable.h b/common/MtzTable.h new file mode 100644 index 000000000..c6dbf7b2e --- /dev/null +++ b/common/MtzTable.h @@ -0,0 +1,48 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "UnitCell.h" + +// What an MTZ file holds - its header, its columns and its batch headers - with no file behind it: the +// merged reflections, the map coefficients and the unmerged observations rugnux makes, for whatever +// writes them (WriteMtz, writer/WriteRugnuxFiles.h) or carries them (the END message). Every dataset +// has the table's cell; a dataset's id is its position in `datasets`. +struct MtzTable { + struct Dataset { + std::string project_name; + std::string crystal_name; + std::string dataset_name; + float wavelength_A = 0.0f; // 0 = not set + }; + struct Column { + std::string label; + char type = 'R'; // the MTZ column type: H, J, Q, F, I, ... + int32_t dataset_id = 0; + std::vector values; // one per row + }; + // The header of one image of unmerged data: the integer and real arrays as the MTZ format lays them + // out (29 and 156 long), and the names of its goniostat axes. + struct Batch { + int32_t number = 0; + std::string title; + std::vector ints; + std::vector floats; + std::vector axes; + }; + + std::string title; + UnitCell unit_cell{}; + std::string space_group_name; // extended Hermann-Mauguin (xhm), which keeps the setting + std::vector history; + std::array sort_order{}; + std::vector datasets; + std::vector columns; + std::vector batches; +}; diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 15d07c0e3..2341bbc3f 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1510,64 +1510,84 @@ namespace { cbor_value_advance(&value); } - MergedReflectionsMessage ProcessMergedElement(CborValue &value) { - MergedReflectionsMessage m; + MtzTable::Dataset ProcessMtzDatasetElement(CborValue &value) { + MtzTable::Dataset d; CborValue map_value; cborErr(cbor_value_enter_container(&value, &map_value)); while (!cbor_value_at_end(&map_value)) { auto key = GetCBORString(map_value); - if (key == "unit_cell") m.unit_cell = ProcessUnitCellElement(map_value); - else if (key == "space_group_name") m.space_group_name = GetCBORString(map_value); - else if (key == "wavelength") m.wavelength_A = GetCBORFloat(map_value); - else if (key == "crystal_name") m.crystal_name = GetCBORString(map_value); - else if (key == "history") GetCBORStringArray(map_value, m.history); - else if (key == "h") GetCBORInt32Array(map_value, m.h); - else if (key == "k") GetCBORInt32Array(map_value, m.k); - else if (key == "l") GetCBORInt32Array(map_value, m.l); - else if (key == "i_mean") GetCBORFloatArray(map_value, m.i_mean); - else if (key == "sigma_i_mean") GetCBORFloatArray(map_value, m.sigma_i_mean); - else if (key == "i_plus") GetCBORFloatArray(map_value, m.i_plus); - else if (key == "sigma_i_plus") GetCBORFloatArray(map_value, m.sigma_i_plus); - else if (key == "i_minus") GetCBORFloatArray(map_value, m.i_minus); - else if (key == "sigma_i_minus") GetCBORFloatArray(map_value, m.sigma_i_minus); - else if (key == "f") GetCBORFloatArray(map_value, m.f); - else if (key == "sigma_f") GetCBORFloatArray(map_value, m.sigma_f); - else if (key == "f_plus") GetCBORFloatArray(map_value, m.f_plus); - else if (key == "sigma_f_plus") GetCBORFloatArray(map_value, m.sigma_f_plus); - else if (key == "f_minus") GetCBORFloatArray(map_value, m.f_minus); - else if (key == "sigma_f_minus") GetCBORFloatArray(map_value, m.sigma_f_minus); - else if (key == "free_flag") GetCBORUInt8Array(map_value, m.free_flag); + if (key == "project_name") d.project_name = GetCBORString(map_value); + else if (key == "crystal_name") d.crystal_name = GetCBORString(map_value); + else if (key == "dataset_name") d.dataset_name = GetCBORString(map_value); + else if (key == "wavelength") d.wavelength_A = GetCBORFloat(map_value); else SkipCBORValue(map_value); } cborErr(cbor_value_leave_container(&value, &map_value)); - return m; + return d; } - MapCoefficientsMessage ProcessMapCoefficientsElement(CborValue &value) { - MapCoefficientsMessage m; + MtzTable::Column ProcessMtzColumnElement(CborValue &value) { + MtzTable::Column c; CborValue map_value; cborErr(cbor_value_enter_container(&value, &map_value)); while (!cbor_value_at_end(&map_value)) { auto key = GetCBORString(map_value); - if (key == "unit_cell") m.unit_cell = ProcessUnitCellElement(map_value); - else if (key == "space_group_name") m.space_group_name = GetCBORString(map_value); - else if (key == "wavelength") m.wavelength_A = GetCBORFloat(map_value); - else if (key == "h") GetCBORInt32Array(map_value, m.h); - else if (key == "k") GetCBORInt32Array(map_value, m.k); - else if (key == "l") GetCBORInt32Array(map_value, m.l); - else if (key == "fp") GetCBORFloatArray(map_value, m.fp); - else if (key == "fc") GetCBORFloatArray(map_value, m.fc); - else if (key == "phic") GetCBORFloatArray(map_value, m.phic); - else if (key == "fwt") GetCBORFloatArray(map_value, m.fwt); - else if (key == "phwt") GetCBORFloatArray(map_value, m.phwt); - else if (key == "delfwt") GetCBORFloatArray(map_value, m.delfwt); - else if (key == "phdelwt") GetCBORFloatArray(map_value, m.phdelwt); - else if (key == "fom") GetCBORFloatArray(map_value, m.fom); - else if (key == "free_flag") GetCBORUInt8Array(map_value, m.free_flag); + if (key == "label") c.label = GetCBORString(map_value); + else if (key == "type") c.type = GetCBORString(map_value).at(0); + else if (key == "dataset_id") c.dataset_id = static_cast(GetCBORInt(map_value)); + else if (key == "values") GetCBORFloatArray(map_value, c.values); else SkipCBORValue(map_value); } cborErr(cbor_value_leave_container(&value, &map_value)); - return m; + return c; + } + + MtzTable::Batch ProcessMtzBatchElement(CborValue &value) { + MtzTable::Batch b; + CborValue map_value; + cborErr(cbor_value_enter_container(&value, &map_value)); + while (!cbor_value_at_end(&map_value)) { + auto key = GetCBORString(map_value); + if (key == "number") b.number = static_cast(GetCBORInt(map_value)); + else if (key == "title") b.title = GetCBORString(map_value); + else if (key == "ints") GetCBORInt32Array(map_value, b.ints); + else if (key == "floats") GetCBORFloatArray(map_value, b.floats); + else if (key == "axes") GetCBORStringArray(map_value, b.axes); + else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return b; + } + + MtzTable ProcessMtzTableElement(CborValue &value) { + MtzTable t; + CborValue map_value; + cborErr(cbor_value_enter_container(&value, &map_value)); + while (!cbor_value_at_end(&map_value)) { + auto key = GetCBORString(map_value); + if (key == "title") t.title = GetCBORString(map_value); + else if (key == "unit_cell") t.unit_cell = ProcessUnitCellElement(map_value); + else if (key == "space_group_name") t.space_group_name = GetCBORString(map_value); + else if (key == "history") GetCBORStringArray(map_value, t.history); + else if (key == "sort_order") { + std::vector v; + GetCBORInt32Array(map_value, v); + if (v.size() != t.sort_order.size()) + throw JFJochException(JFJochExceptionCategory::CBORError, "MTZ sort order of wrong length"); + std::copy(v.begin(), v.end(), t.sort_order.begin()); + } else if (key == "datasets" || key == "columns" || key == "batches") { + CborValue array_value; + cborErr(cbor_value_enter_container(&map_value, &array_value)); + while (!cbor_value_at_end(&array_value)) { + if (key == "datasets") t.datasets.push_back(ProcessMtzDatasetElement(array_value)); + else if (key == "columns") t.columns.push_back(ProcessMtzColumnElement(array_value)); + else t.batches.push_back(ProcessMtzBatchElement(array_value)); + } + cborErr(cbor_value_leave_container(&map_value, &array_value)); + } else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return t; } PlacedModelMessage ProcessPlacedModelElement(CborValue &value) { @@ -1699,9 +1719,9 @@ namespace { else if (key == "pixel_sum") GetCBORInt64Array(value, message.pixel_sum); else if (key == "merged") - message.merged = ProcessMergedElement(value); + message.merged = ProcessMtzTableElement(value); else if (key == "map_coefficients") - message.map_coefficients = ProcessMapCoefficientsElement(value); + message.map_coefficients = ProcessMtzTableElement(value); else if (key == "placed_model") message.placed_model = ProcessPlacedModelElement(value); else if (key == "rotation_lattice") { diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 82ee03e51..eaf732fd1 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -521,53 +521,52 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector(v.data()), v.size() * sizeof(double))); } -inline void CBOR_ENC(CborEncoder &encoder, const char* key, const MergedReflectionsMessage &m) { - CborEncoder mapEncoder; +inline void CBOR_ENC(CborEncoder &encoder, const char* key, const MtzTable &t) { + CborEncoder mapEncoder, arrayEncoder, elementEncoder; cborErr(cbor_encode_text_stringz(&encoder, key)); cborErr(cbor_encoder_create_map(&encoder, &mapEncoder, CborIndefiniteLength)); - CBOR_ENC(mapEncoder, "unit_cell", m.unit_cell); - CBOR_ENC(mapEncoder, "space_group_name", m.space_group_name); - CBOR_ENC(mapEncoder, "wavelength", m.wavelength_A); - CBOR_ENC(mapEncoder, "crystal_name", m.crystal_name); - CBOR_ENC(mapEncoder, "history", m.history); - CBOR_ENC(mapEncoder, "h", m.h); - CBOR_ENC(mapEncoder, "k", m.k); - CBOR_ENC(mapEncoder, "l", m.l); - CBOR_ENC(mapEncoder, "i_mean", m.i_mean); - CBOR_ENC(mapEncoder, "sigma_i_mean", m.sigma_i_mean); - CBOR_ENC(mapEncoder, "i_plus", m.i_plus); - CBOR_ENC(mapEncoder, "sigma_i_plus", m.sigma_i_plus); - CBOR_ENC(mapEncoder, "i_minus", m.i_minus); - CBOR_ENC(mapEncoder, "sigma_i_minus", m.sigma_i_minus); - CBOR_ENC(mapEncoder, "f", m.f); - CBOR_ENC(mapEncoder, "sigma_f", m.sigma_f); - CBOR_ENC(mapEncoder, "f_plus", m.f_plus); - CBOR_ENC(mapEncoder, "sigma_f_plus", m.sigma_f_plus); - CBOR_ENC(mapEncoder, "f_minus", m.f_minus); - CBOR_ENC(mapEncoder, "sigma_f_minus", m.sigma_f_minus); - CBOR_ENC(mapEncoder, "free_flag", m.free_flag); - cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); -} + CBOR_ENC(mapEncoder, "title", t.title); + CBOR_ENC(mapEncoder, "unit_cell", t.unit_cell); + CBOR_ENC(mapEncoder, "space_group_name", t.space_group_name); + CBOR_ENC(mapEncoder, "history", t.history); + CBOR_ENC(mapEncoder, "sort_order", std::vector(t.sort_order.begin(), t.sort_order.end())); -inline void CBOR_ENC(CborEncoder &encoder, const char* key, const MapCoefficientsMessage &m) { - CborEncoder mapEncoder; - cborErr(cbor_encode_text_stringz(&encoder, key)); - cborErr(cbor_encoder_create_map(&encoder, &mapEncoder, CborIndefiniteLength)); - CBOR_ENC(mapEncoder, "unit_cell", m.unit_cell); - CBOR_ENC(mapEncoder, "space_group_name", m.space_group_name); - CBOR_ENC(mapEncoder, "wavelength", m.wavelength_A); - CBOR_ENC(mapEncoder, "h", m.h); - CBOR_ENC(mapEncoder, "k", m.k); - CBOR_ENC(mapEncoder, "l", m.l); - CBOR_ENC(mapEncoder, "fp", m.fp); - CBOR_ENC(mapEncoder, "fc", m.fc); - CBOR_ENC(mapEncoder, "phic", m.phic); - CBOR_ENC(mapEncoder, "fwt", m.fwt); - CBOR_ENC(mapEncoder, "phwt", m.phwt); - CBOR_ENC(mapEncoder, "delfwt", m.delfwt); - CBOR_ENC(mapEncoder, "phdelwt", m.phdelwt); - CBOR_ENC(mapEncoder, "fom", m.fom); - CBOR_ENC(mapEncoder, "free_flag", m.free_flag); + cborErr(cbor_encode_text_stringz(&mapEncoder, "datasets")); + cborErr(cbor_encoder_create_array(&mapEncoder, &arrayEncoder, t.datasets.size())); + for (const auto &d : t.datasets) { + cborErr(cbor_encoder_create_map(&arrayEncoder, &elementEncoder, 4)); + CBOR_ENC(elementEncoder, "project_name", d.project_name); + CBOR_ENC(elementEncoder, "crystal_name", d.crystal_name); + CBOR_ENC(elementEncoder, "dataset_name", d.dataset_name); + CBOR_ENC(elementEncoder, "wavelength", d.wavelength_A); + cborErr(cbor_encoder_close_container(&arrayEncoder, &elementEncoder)); + } + cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); + + cborErr(cbor_encode_text_stringz(&mapEncoder, "columns")); + cborErr(cbor_encoder_create_array(&mapEncoder, &arrayEncoder, t.columns.size())); + for (const auto &c : t.columns) { + cborErr(cbor_encoder_create_map(&arrayEncoder, &elementEncoder, 4)); + CBOR_ENC(elementEncoder, "label", c.label); + CBOR_ENC(elementEncoder, "type", std::string(1, c.type)); + CBOR_ENC(elementEncoder, "dataset_id", static_cast(c.dataset_id)); + CBOR_ENC(elementEncoder, "values", c.values); + cborErr(cbor_encoder_close_container(&arrayEncoder, &elementEncoder)); + } + cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); + + cborErr(cbor_encode_text_stringz(&mapEncoder, "batches")); + cborErr(cbor_encoder_create_array(&mapEncoder, &arrayEncoder, t.batches.size())); + for (const auto &b : t.batches) { + cborErr(cbor_encoder_create_map(&arrayEncoder, &elementEncoder, 5)); + CBOR_ENC(elementEncoder, "number", static_cast(b.number)); + CBOR_ENC(elementEncoder, "title", b.title); + CBOR_ENC(elementEncoder, "ints", b.ints); + CBOR_ENC(elementEncoder, "floats", b.floats); + CBOR_ENC(elementEncoder, "axes", b.axes); + cborErr(cbor_encoder_close_container(&arrayEncoder, &elementEncoder)); + } + cborErr(cbor_encoder_close_container(&mapEncoder, &arrayEncoder)); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); } diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 8cf24c2a5..53e4f69af 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -607,39 +607,48 @@ void WriteMmcifReflections(const std::vector &reflections, out.close(); } -MergedReflectionsMessage MergedReflectionsForWriter(const std::vector &reflections, - const UnitCell &unitCell, - const DiffractionExperiment &experiment, - const OutputProvenance &provenance) { - MergedReflectionsMessage m; +MtzTable MergedMtzTable(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const OutputProvenance &provenance) { + MtzTable m; m.unit_cell = unitCell; m.space_group_name = experiment.GetSpaceGroupOrP1().xhm(); - m.wavelength_A = experiment.GetWavelength_A(); - m.crystal_name = experiment.GetSampleName(); + // Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention). + m.title = "Rugnux merged reflections"; m.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); AddProvenanceHistory(m.history, provenance); + // HKL_base first, so the data dataset is id 1: a dataset written as id 0 is the reserved base dataset + // as far as CCP4's mtzlib is concerned, and the wavelength on it is dropped - mtzinfo then reports the + // Cu K-alpha default 1.54187 A, which is what a program that takes f'/f'' from the file would phase on. + m.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f}, + {"native", experiment.GetSampleName(), "native", experiment.GetWavelength_A()}}; + // One row per reflection in the CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for + // the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly. bool has_anom = true; const std::vector out_rows = BuildMergedRows(reflections, experiment, has_anom); - for (const auto& row : out_rows) { - m.h.push_back(row.h); - m.k.push_back(row.k); - m.l.push_back(row.l); - m.i_mean.push_back(row.Imean); - m.sigma_i_mean.push_back(row.sImean); - m.f.push_back(row.Fmean); - m.sigma_f.push_back(row.sFmean); - if (has_anom) { - m.i_plus.push_back(row.Ip); - m.sigma_i_plus.push_back(row.sIp); - m.i_minus.push_back(row.Im); - m.sigma_i_minus.push_back(row.sIm); - m.f_plus.push_back(row.Fp); - m.sigma_f_plus.push_back(row.sFp); - m.f_minus.push_back(row.Fm); - m.sigma_f_minus.push_back(row.sFm); - } - m.free_flag.push_back(static_cast(row.rfree)); + std::vector> columns{{"H", 'H'}, {"K", 'H'}, {"L", 'H'}, + {"IMEAN", 'J'}, {"SIGIMEAN", 'Q'}}; + if (has_anom) + columns.insert(columns.end(), {{"I(+)", 'K'}, {"SIGI(+)", 'M'}, {"I(-)", 'K'}, {"SIGI(-)", 'M'}}); + columns.insert(columns.end(), {{"F", 'F'}, {"SIGF", 'Q'}}); // F: the French-Wilson amplitude + if (has_anom) + columns.insert(columns.end(), {{"F(+)", 'G'}, {"SIGF(+)", 'L'}, {"F(-)", 'G'}, {"SIGF(-)", 'L'}}); + columns.push_back({"FreeR_flag", 'I'}); + for (const auto &[label, type] : columns) + m.columns.push_back({label, type, type == 'H' ? 0 : 1, {}}); + for (const auto &row : out_rows) { + std::vector v{static_cast(row.h), static_cast(row.k), static_cast(row.l), + row.Imean, row.sImean}; + if (has_anom) + v.insert(v.end(), {row.Ip, row.sIp, row.Im, row.sIm}); + v.insert(v.end(), {row.Fmean, row.sFmean}); + if (has_anom) + v.insert(v.end(), {row.Fp, row.sFp, row.Fm, row.sFm}); + v.push_back(static_cast(row.rfree)); + for (size_t c = 0; c < v.size(); c++) + m.columns[c].values.push_back(v[c]); } return m; } @@ -649,7 +658,7 @@ void WriteMtzReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const OutputProvenance &provenance, const std::string &filename) { - WriteMergedMtz(MergedReflectionsForWriter(reflections, unitCell, experiment, provenance), filename); + WriteMtz(MergedMtzTable(reflections, unitCell, experiment, provenance), filename); } namespace { @@ -906,22 +915,22 @@ std::vector> SumRockingEvents(const std::vector &outcomes, - const UnitCell &unitCell, - const DiffractionExperiment &experiment, - bool sum_partials, - const OutputProvenance &provenance, - size_t nthreads) { - gemmi::Mtz mtz; - mtz.spacegroup = &experiment.GetSpaceGroupOrP1(); - mtz.set_cell_for_all(unitCell); +MtzTable UnmergedMtzTable(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + bool sum_partials, + const OutputProvenance &provenance, + size_t nthreads) { + const gemmi::SpaceGroup *spacegroup = &experiment.GetSpaceGroupOrP1(); + MtzTable mtz; + mtz.space_group_name = spacegroup->xhm(); + mtz.unit_cell = unitCell; mtz.title = "Rugnux unmerged observations"; mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); AddProvenanceHistory(mtz.history, provenance); - mtz.add_base(); // the HKL_base dataset and the H K L columns - gemmi::Mtz::Dataset &ds = mtz.add_dataset("unmerged"); - ds.crystal_name = experiment.GetSampleName(); - ds.wavelength = experiment.GetWavelength_A(); + mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f}, + {"unmerged", experiment.GetSampleName(), "unmerged", experiment.GetWavelength_A()}}; + const int dataset_id = 1; // Every data column belongs to HKL_base and only the batches carry a dataset id, as in a // POINTLESS file. Down to FLAG this is POINTLESS's own column set, plus QE and FLIGHT; the four after it @@ -941,22 +950,11 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, // rises above 1 toward the detector edge where a longer crossing makes the sensor more // efficient. DIALS stores the un-normalised absorbed fraction instead; the two differ by the // per-dataset constant QE(0), which is an overall scale. - mtz.add_column("M/ISYM", 'Y', 0, -1, false); - mtz.add_column("BATCH", 'B', 0, -1, false); - mtz.add_column("I", 'J', 0, -1, false); - mtz.add_column("SIGI", 'Q', 0, -1, false); - mtz.add_column("FRACTIONCALC", 'R', 0, -1, false); - mtz.add_column("XDET", 'R', 0, -1, false); - mtz.add_column("YDET", 'R', 0, -1, false); - mtz.add_column("ROT", 'R', 0, -1, false); - mtz.add_column("LP", 'R', 0, -1, false); - mtz.add_column("QE", 'R', 0, -1, false); - mtz.add_column("FLIGHT", 'R', 0, -1, false); - mtz.add_column("FLAG", 'I', 0, -1, false); - mtz.add_column("DELPHI", 'R', 0, -1, false); - mtz.add_column("ZETA", 'R', 0, -1, false); - mtz.add_column("BGMEAN", 'R', 0, -1, false); - mtz.add_column("BGVAR", 'R', 0, -1, false); + for (const auto &[label, type] : std::vector>{ + {"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"M/ISYM", 'Y'}, {"BATCH", 'B'}, {"I", 'J'}, {"SIGI", 'Q'}, + {"FRACTIONCALC", 'R'}, {"XDET", 'R'}, {"YDET", 'R'}, {"ROT", 'R'}, {"LP", 'R'}, {"QE", 'R'}, + {"FLIGHT", 'R'}, {"FLAG", 'I'}, {"DELPHI", 'R'}, {"ZETA", 'R'}, {"BGMEAN", 'R'}, {"BGVAR", 'R'}}) + mtz.columns.push_back({label, type, 0, {}}); const auto gon = experiment.GetGoniometer(); // A scanning goniometer means rugnux integrated a rocking curve image by image, so each @@ -985,7 +983,7 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, // None of the three is a scale, and a program that reads this file has no way to recompute them. The // partiality is NOT divided out - that IS a scale, FRACTIONCALC carries it, and every program // this file is for wants to handle it its own way. - gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup); + gemmi::UnmergedHklMover hkl_mover(spacegroup); // The batch headers are written for every image the observations span, not only for the images // that produced one. AIMLESS starts a new run wherever the phi series jumps, so an image that // indexed nothing leaves a hole that cuts the sweep in two: measured on a 360 deg sweep with 67 @@ -1000,7 +998,7 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, // roughly halve AIMLESS's . Screw and glide absences ARE written: they are the evidence // the space group was chosen on, so a program reading this file can check that choice, and // removing them turns its test into an assumption. XDS and DIALS draw the line in the same place. - const char centering = mtz.spacegroup ? mtz.spacegroup->hm[0] : 'P'; + const char centering = spacegroup->hm[0]; // On stills only the crystals the merge takes are written: one whose cell is not the run's was indexed // in another lattice (on few spots, typically an axis multiple), so its indices do not name the @@ -1040,13 +1038,14 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, add_row(r); } - // Every row has its own slot, so they are filled on several workers. + // Every row has its own slot in each column, so they are filled on several workers. const size_t w = mtz.columns.size(); - mtz.data.resize(rows.size() * w); + for (auto &c : mtz.columns) + c.values.resize(rows.size()); ParallelChunks(static_cast(rows.size()), ThreadsForWork(rows.size(), nthreads), [&](int lo, int hi) { + std::vector row(w); for (int i = lo; i < hi; ++i) { const Reflection &r = *rows[i]; - float *row = &mtz.data[static_cast(i) * w]; std::array hkl{r.h, r.k, r.l}; const int isym = hkl_mover.move_to_asu(hkl); row[0] = static_cast(hkl[0]); @@ -1073,9 +1072,10 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, row[16] = r.zeta; row[17] = r.bkg; row[18] = r.var_bkg; + for (size_t c = 0; c < w; c++) + mtz.columns[c].values[i] = row[c]; } }); - mtz.nreflections = static_cast(mtz.data.size() / mtz.columns.size()); // The batch header's orientation matrix is the crystal at rotation angle zero - each batch's own // PHISTT is applied on top of it - but the lattice stored with an outcome is the crystal as it @@ -1106,7 +1106,7 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, gemmi::Mtz::Batch batch; batch.title = "Rugnux"; batch.axes.emplace_back("PHI"); - batch.set_dataset_id(ds.id); + batch.set_dataset_id(dataset_id); batch.ints[12] = 1; // ncryst batch.ints[14] = partials ? 1 : 2; // ldtype: oscillation data (2D spots) / area detector (3D) batch.ints[15] = 1; // jsaxs: the goniostat scan axis @@ -1143,38 +1143,41 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, batch.number = number; batch.floats[36] = phi_start_deg(static_cast(number - 1)); // phistt batch.floats[37] = batch.floats[36] + wedge_deg; // phiend - mtz.batches.push_back(batch); + mtz.batches.push_back({batch.number, batch.title, {batch.ints.begin(), batch.ints.end()}, batch.floats, + batch.axes}); } // By H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in. This is Mtz::sort(5) - // a stable sort on those five columns - done on several workers: rows equal in all five keep the // order they were added in, so with the row number as the last key the order is total and the // parallel sort returns the same rows in the same order. - if (!mtz.has_data()) - mtz.sort(5); // no rows: its own "No data." failure, as before + if (rows.empty()) + gemmi::fail("No data."); // as gemmi's own sort has it { - std::vector order(mtz.nreflections); + std::vector order(rows.size()); std::iota(order.begin(), order.end(), 0); ParallelSort(order.begin(), order.end(), nthreads, [&](int i, int j) { - const float *a = &mtz.data[static_cast(i) * w], *b = &mtz.data[static_cast(j) * w]; - for (int n = 0; n < 5; ++n) - if (a[n] != b[n]) - return a[n] < b[n]; + for (int n = 0; n < 5; ++n) { + const float a = mtz.columns[n].values[i], b = mtz.columns[n].values[j]; + if (a != b) + return a < b; + } return i < j; }); - std::vector sorted(mtz.data.size()); - ParallelChunks(static_cast(order.size()), ThreadsForWork(order.size(), nthreads), [&](int lo, int hi) { - for (int i = lo; i < hi; ++i) - std::memcpy(&sorted[static_cast(i) * w], &mtz.data[static_cast(order[i]) * w], - w * sizeof(float)); - }); - mtz.data.swap(sorted); - mtz.sort_order = {{1, 2, 3, 4, 5}}; + for (auto &c : mtz.columns) { + std::vector sorted(order.size()); + ParallelChunks(static_cast(order.size()), ThreadsForWork(order.size(), nthreads), [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) + sorted[i] = c.values[order[i]]; + }); + c.values.swap(sorted); + } + mtz.sort_order = {1, 2, 3, 4, 5}; } return mtz; } -void SetUnmergedMtzMosaicity(gemmi::Mtz &mtz, const std::vector &outcomes) { +void SetUnmergedMtzMosaicity(MtzTable &mtz, const std::vector &outcomes) { // The reflecting range of the same image the orientation matrix is taken from. std::optional mosaicity_deg; for (const auto &outcome : outcomes) { @@ -1194,9 +1197,9 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes const OutputProvenance &provenance, const std::string &filename, size_t nthreads) { - gemmi::Mtz mtz = UnmergedMtz(outcomes, unitCell, experiment, sum_partials, provenance, nthreads); + MtzTable mtz = UnmergedMtzTable(outcomes, unitCell, experiment, sum_partials, provenance, nthreads); SetUnmergedMtzMosaicity(mtz, outcomes); - mtz.write_to_file(filename); + WriteMtz(mtz, filename); } void WriteReflections(const std::vector &reflections, diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 68d2ab778..1d9053fd8 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -9,11 +9,9 @@ #include "../common/Reflection.h" #include "../common/UnitCell.h" #include "../common/DiffractionExperiment.h" -#include "../common/JFJochMessages.h" +#include "../common/MtzTable.h" #include "IntegrationOutcome.h" -#include - struct MergeStatistics; struct TwinningAnalysisResult; @@ -50,13 +48,14 @@ void WriteMmcifReflections(const std::vector &reflections, const std::string &filename, size_t nthreads); -// The merged MTZ as a structured END-message block: the columns, metadata and HISTORY records that -// WriteMtzReflections writes, for jfjoch_writer to write the same file (WriteMergedMtz, -// writer/WriteRugnuxFiles.h). WriteMtzReflections is exactly WriteMergedMtz of this. -MergedReflectionsMessage MergedReflectionsForWriter(const std::vector &reflections, - const UnitCell &unitCell, - const DiffractionExperiment &experiment, - const OutputProvenance &provenance); +// The merged MTZ: H K L IMEAN SIGIMEAN [I(+) SIGI(+) I(-) SIGI(-)] F SIGF [F(+) SIGF(+) F(-) SIGF(-)] +// FreeR_flag, one row per unique reflection in the CCP4 anomalous layout, the bracketed columns only where +// the merge kept a Bijvoet split; with the HISTORY records saying where it came from. WriteMtzReflections +// writes it (WriteMtz, writer/WriteRugnuxFiles.h). +MtzTable MergedMtzTable(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const OutputProvenance &provenance); void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, @@ -98,17 +97,16 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes const std::string &filename, size_t nthreads = 1); -// The same file in two steps, for a caller that builds it while the outcomes' per-image scale fields -// are still being written: UnmergedMtz reads everything but mosaicity_deg, which -// SetUnmergedMtzMosaicity then copies into the batch headers. The Mtz points at the experiment's space -// group, so the experiment has to outlive it. -gemmi::Mtz UnmergedMtz(const std::vector &outcomes, - const UnitCell &unitCell, - const DiffractionExperiment &experiment, - bool sum_partials, - const OutputProvenance &provenance, - size_t nthreads = 1); -void SetUnmergedMtzMosaicity(gemmi::Mtz &mtz, const std::vector &outcomes); +// The same file as its table, in two steps, for a caller that builds it while the outcomes' per-image +// scale fields are still being written: UnmergedMtzTable reads everything but mosaicity_deg, which +// SetUnmergedMtzMosaicity then copies into the batch headers. +MtzTable UnmergedMtzTable(const std::vector &outcomes, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + bool sum_partials, + const OutputProvenance &provenance, + size_t nthreads = 1); +void SetUnmergedMtzMosaicity(MtzTable &mtz, const std::vector &outcomes); void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 86dc97067..b16a5e24a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -12,7 +12,7 @@ #include #include #include -#include +#include #include #include #include @@ -303,10 +303,10 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, std::optional reference_cell, int scaling_iterations, size_t nthreads, Logger &logger, - std::string observation_dump_path) + std::function observation_dump) : x(experiment), partials_out(partial_outcomes), reference_cell(std::move(reference_cell)), nthreads(nthreads == 0 ? AvailableCpus() : nthreads), logger(&logger), - observation_dump_path(std::move(observation_dump_path)) { + observation_dump(std::move(observation_dump)) { const auto s = x.GetScalingSettings(); min_partiality = s.GetMinPartiality(); d_min_limit = s.GetHighResolutionLimit_A(); @@ -401,7 +401,7 @@ void RotationScaleMerge::IngestHost() { // 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_; - resident_ingest = gpu_active_ && observation_dump_path.empty(); + resident_ingest = gpu_active_ && !observation_dump; #endif frame_cell_ok.assign(n_frames, 1); @@ -4403,7 +4403,7 @@ void RotationScaleMerge::Combine() { // Combine one raw-hkl run into fulls, appended to `out`. Independent per run (the events of one hkl // touch no shared state), so runs parallelise cleanly. Returns the number of usable partials seen. // `dump` (serial path only) writes each emitted full for the diagnostic observation dump. - auto process_rawrun = [&](int r, std::vector &out, std::ofstream *dump, int64_t &n_overloaded) -> size_t { + auto process_rawrun = [&](int r, std::vector &out, std::ostream *dump, int64_t &n_overloaded) -> size_t { prefetch_run(r + COMBINE_PREFETCH_RUNS); const int lo = rawrun_start[r], hi = rawrun_start[r] + rawrun_count[r]; std::vector ev; // usable perm-indices of this raw hkl, in image-number order @@ -4545,11 +4545,10 @@ void RotationScaleMerge::Combine() { size_t n_used = 0; std::vector> part; // the fulls of consecutive runs, in run order - if (!observation_dump_path.empty() || nthreads <= 1) { + if (observation_dump || nthreads <= 1) { // Serial (diagnostic dump needs a single writer). - std::ofstream dump; - if (!observation_dump_path.empty()) { - dump.open(observation_dump_path); + std::ostringstream dump; + if (observation_dump) { // Columns are enough to rescale/merge the combined fulls in an external program: peak_frame // is a per-full batch (ZD) and px/py the predicted detector position (XD/YD). dump << "# h k l I sigma d n_frames captured_fraction peak_frame px py\n"; @@ -4557,7 +4556,9 @@ void RotationScaleMerge::Combine() { part.resize(1); part[0].reserve(n_run); for (int r = 0; r < n_run; ++r) - n_used += process_rawrun(r, part[0], dump.is_open() ? &dump : nullptr, fulls_dropped_overloaded); + n_used += process_rawrun(r, part[0], observation_dump ? &dump : nullptr, fulls_dropped_overloaded); + if (observation_dump) + observation_dump(dump.str()); } else { // Parallel over contiguous rawrun chunks; the per-thread fulls are taken in run order below, so // the result is deterministic. @@ -6837,7 +6838,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool // scale-fulls loop below. Only the diagnostic dump (serial, one writer) falls back to CPU. The // decay/absorption correction surfaces run on the downloaded fulls (px/py is carried through the // combine) and their corrected corr is re-uploaded before the resident merge - so they stay on the GPU. - if (gpu_active_ && observation_dump_path.empty()) { + if (gpu_active_ && !observation_dump) { // The smoothed corr is already resident (scaling + smooth-G ran on the device, no round-trip). const int nf = gpu_->Combine(rawrun_group.data(), min_partiality, capture_uncertainty_coeff, min_captured_fraction, max_frame_gap, fulls_dropped_overloaded); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index b1aa6d658..820d3fb7a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -136,7 +136,7 @@ public: int scaling_iterations, size_t nthreads, Logger &logger, - std::string observation_dump_path = {}); + std::function observation_dump = {}); // Copy the per-frame partials into the flat buffers. Call once before the first Run(). void Ingest(); @@ -358,7 +358,8 @@ private: std::optional reference_cell; size_t nthreads; Logger *logger; - std::string observation_dump_path; + // Handed the text of the diagnostic observation dump each time the combine runs (CPU only). + std::function observation_dump; // The group the current Run merges in (see Run); null outside a Run that was given one. const gemmi::SpaceGroup *run_group = nullptr; const gemmi::SpaceGroup &Group() const { return run_group ? *run_group : x.GetSpaceGroupOrP1(); } diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index 0497730c2..6044222da 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -27,8 +27,6 @@ #include // get_size_for_hkl, get_f_phi_on_grid (the maps) #include // SolventMasker #include // Scaling (bulk solvent + anisotropic B) -#include // Ccp4 map I/O -#include // Mtz (map-coefficient output) #include // calculate_data_statistics #include "../../common/CorrelationCoefficient.h" @@ -79,22 +77,9 @@ gemmi::Grid map_from_coefficients(gemmi::AsuData> &co return MapFromFPhi(gemmi::get_f_phi_on_grid(coef, size, true)); } -// A file this run may already have written - the validation in the model's setting writes over the -// first one's maps - is removed before it is written again, not truncated: XFS (and ext4) flush a file -// truncated and rewritten when it is closed, and three 150 MB maps forced out to a disk that way cost -// seconds where a fresh file costs nothing. -void remove_before_rewriting(const std::string &path) { - std::remove(path.c_str()); -} - -// Write a map as CCP4; return its RMS (the sigma the map is read in). -double write_ccp4(const gemmi::Grid &map, const std::string &path) { - gemmi::Ccp4 ccp4; - ccp4.grid = map; - ccp4.update_ccp4_header(2); - remove_before_rewriting(path); - ccp4.write_ccp4_map(path); - return ccp4.hstats.rms; +// The RMS of a map, the sigma it is read in: the one its CCP4 header states (Ccp4::update_ccp4_header). +double map_rms(const gemmi::Grid &map) { + return gemmi::calculate_data_statistics(map.data).rms; } // Cubic, not the default linear, for reading a map at a point. The maps are sampled every d_min/3, @@ -536,7 +521,6 @@ namespace { ModelValidationResult Validate(const std::vector &merged, const UnitCell &cell, const RugnuxModel &model, - const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group, bool probe_indexing_ambiguity, @@ -1711,23 +1695,27 @@ ModelValidationResult Validate(const std::vector &merged, mapfofc.v.push_back({terms[i].hkl, delfwt[i] * ph}); } - // A validation started before the one it follows had decided writes nothing until that one has, - // and nothing at all where it decided otherwise. - if (schedule.write_gate.valid() && !schedule.write_gate.get()) { - result.failure_reason = "superseded before its files were written"; + // A validation started before the one it follows had decided makes no maps until that one has, + // and none at all where it decided otherwise. + if (schedule.maps_gate.valid() && !schedule.maps_gate.get()) { + result.failure_reason = "superseded before its maps were made"; return result; } - // Where the validation in the model's setting, started beside this one, is going to be kept, it - // writes the maps and the map coefficients on the axes the files are written on, and this one's would - // be overwritten unread: they are not made at all. + // Where the validation in the model's setting, started beside this one, is going to be kept, its + // maps and map coefficients are on the axes the files are written on, and this one's would be + // thrown away unread: they are not made at all. const bool maps_superseded = schedule.maps_superseded && schedule.maps_superseded(result); + std::shared_ptr maps; if (!maps_superseded) { - // --- write the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- + maps = std::make_shared(); + // --- make the maps and score the 2mFo-DFc map at atom centres (a real map peaks there) --- // The three maps - this one, the difference map and the anomalous map - share nothing but what - // they read, so the other two are made and written on threads of their own beside this one. + // they read, so the other two are made on threads of their own beside this one. Nothing here + // reads the difference map, so it is made only where it is kept. std::future fofc_map = std::async(std::launch::async, [&] { - write_ccp4(map_from_coefficients(mapfofc, sf_gpu), output_prefix + "_fofc.ccp4"); + if (schedule.keep_map_grids) + maps->map_fofc = map_from_coefficients(mapfofc, sf_gpu); }); // --- anomalous difference map, where the merge kept the Bijvoet split --- @@ -1741,13 +1729,15 @@ ModelValidationResult Validate(const std::vector &merged, result.indexing_op, fmodel); result.anomalous_pairs = static_cast(mapanom.v.size()); if (!mapanom.v.empty()) { - const gemmi::Grid grid = map_from_coefficients(mapanom, sf_gpu); - read_anomalous_map(grid, write_ccp4(grid, output_prefix + "_anom.ccp4"), st, result); + gemmi::Grid grid = map_from_coefficients(mapanom, sf_gpu); + read_anomalous_map(grid, map_rms(grid), st, result); + if (schedule.keep_map_grids) + maps->map_anom = std::move(grid); } }); - const gemmi::Grid grid2fofc = map_from_coefficients(map2fofc, sf_gpu); - const double rms2 = write_ccp4(grid2fofc, output_prefix + "_2fofc.ccp4"); + gemmi::Grid grid2fofc = map_from_coefficients(map2fofc, sf_gpu); + const double rms2 = map_rms(grid2fofc); { double s = 0; int n = 0; for (gemmi::Model &m : st.models) @@ -1756,6 +1746,8 @@ ModelValidationResult Validate(const std::vector &merged, for (gemmi::Atom &a : r.atoms) { s += grid2fofc.interpolate_value(a.pos, MAP_INTERPOLATION_ORDER); ++n; } result.mean_atom_density_sigma = (n > 0 && rms2 > 0) ? (s / n) / rms2 : 0; } + if (schedule.keep_map_grids) + maps->map_2fofc = std::move(grid2fofc); fofc_map.get(); anomalous_map.get(); } @@ -1780,60 +1772,53 @@ ModelValidationResult Validate(const std::vector &merged, : "the model does not fit these data"); if (!maps_superseded) { - // --- MTZ of map coefficients so the maps can be re-opened / rebuilt in Coot etc. --- - try { - gemmi::Mtz mtz(true); - // The group the REFLECTIONS end up in, which is the model's only where its hand was adopted - - // AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes - // identically, so the coefficients are the same numbers either way and only the label moves; - // but the two groups have different screw translations, so a reader that expands symmetry out - // of this file works in the wrong one if the label disagrees with the .mtz beside it. - mtz.spacegroup = result.adopted_model_enantiomorph - ? sg - // No data group means the caller merged in P1, and P1 is what the - // reflection files beside this one carry - not the model's group. - : (data_space_group != nullptr ? data_space_group - : gemmi::find_spacegroup_by_number(1)); - mtz.set_cell_for_all(ucell); - mtz.add_dataset("model_validation"); - mtz.datasets.back().wavelength = wavelength_A; - mtz.add_column("FP", 'F', -1, -1, false); - mtz.add_column("FC", 'F', -1, -1, false); - mtz.add_column("PHIC", 'P', -1, -1, false); - mtz.add_column("FWT", 'F', -1, -1, false); - mtz.add_column("PHWT", 'P', -1, -1, false); - mtz.add_column("DELFWT", 'F', -1, -1, false); - mtz.add_column("PHDELWT", 'P', -1, -1, false); - // The figure of merit the coefficients carry, so the weighting can be read off - and undone - - // from the file rather than having to be taken on trust. - mtz.add_column("FOM", 'W', -1, -1, false); - mtz.add_column("FREE", 'I', -1, -1, false); - std::vector data; - for (size_t i = 0; i < terms.size(); i++) { - // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg - // returns (-180, 180] and a P column is not supposed to. - const double phi_wrapped = terms[i].phi * 180.0 / PI; - const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); - data.insert(data.end(), {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), - static_cast(terms[i].hkl[2]), - static_cast(terms[i].fo), static_cast(terms[i].fc), - phi_deg, - fwt[i], phi_deg, - delfwt[i], phi_deg, - static_cast(sigma_a.weight[i].m), - terms[i].free ? 0.0f : 1.0f}); - } - mtz.nreflections = static_cast(terms.size()); - mtz.data = std::move(data); - remove_before_rewriting(output_prefix + "_maps.mtz"); - mtz.write_to_file(output_prefix + "_maps.mtz"); - } catch (const std::exception &e) { - logger.Warning("Model validation: could not write map MTZ: {}", e.what()); + // --- the map coefficients, so the maps can be re-opened / rebuilt in Coot etc. --- + MtzTable &mtz = maps->coefficients; + // The group the REFLECTIONS end up in, which is the model's only where its hand was adopted - + // AdoptModelFrame decides this the same way a few lines below. An enantiomorphic pair indexes + // identically, so the coefficients are the same numbers either way and only the label moves; + // but the two groups have different screw translations, so a reader that expands symmetry out + // of this file works in the wrong one if the label disagrees with the .mtz beside it. + mtz.space_group_name = result.adopted_model_enantiomorph + ? sg->xhm() + // No data group means the caller merged in P1, and P1 is what the + // reflection files beside this one carry - not the model's group. + : (data_space_group != nullptr ? data_space_group->xhm() + : gemmi::find_spacegroup_by_number(1)->xhm()); + mtz.unit_cell = UnitCell{static_cast(ucell.a), static_cast(ucell.b), + static_cast(ucell.c), static_cast(ucell.alpha), + static_cast(ucell.beta), static_cast(ucell.gamma)}; + mtz.datasets = {{"HKL_base", "HKL_base", "HKL_base", 0.0f}, + {"model_validation", "model_validation", "model_validation", + static_cast(wavelength_A)}}; + // The figure of merit the coefficients carry, so the weighting can be read off - and undone - + // from the file rather than having to be taken on trust. + const std::pair columns[] = {{"H", 'H'}, {"K", 'H'}, {"L", 'H'}, {"FP", 'F'}, + {"FC", 'F'}, {"PHIC", 'P'}, {"FWT", 'F'}, {"PHWT", 'P'}, + {"DELFWT", 'F'}, {"PHDELWT", 'P'}, {"FOM", 'W'}, + {"FREE", 'I'}}; + for (const auto &[label, type] : columns) + mtz.columns.push_back({label, type, type == 'H' ? 0 : 1, {}}); + for (size_t i = 0; i < terms.size(); i++) { + // [0, 360), the convention the MTZ 'P' columns carried before sigma_A weighting: std::arg + // returns (-180, 180] and a P column is not supposed to. + const double phi_wrapped = terms[i].phi * 180.0 / PI; + const auto phi_deg = static_cast(phi_wrapped < 0 ? phi_wrapped + 360.0 : phi_wrapped); + const float row[] = {static_cast(terms[i].hkl[0]), static_cast(terms[i].hkl[1]), + static_cast(terms[i].hkl[2]), + static_cast(terms[i].fo), static_cast(terms[i].fc), + phi_deg, + fwt[i], phi_deg, + delfwt[i], phi_deg, + static_cast(sigma_a.weight[i].m), + terms[i].free ? 0.0f : 1.0f}; + for (size_t c = 0; c < mtz.columns.size(); c++) + mtz.columns[c].values.push_back(row[c]); } } result.ok = true; - result.maps_prefix = output_prefix; + result.maps = std::move(maps); // The placed coordinates, for the caller to write out beside the maps once the frame is settled. result.placed_model = std::make_shared(st); logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) " @@ -1862,12 +1847,7 @@ ModelValidationResult Validate(const std::vector &merged, "data were indexed in the wrong hand, and reindexing would hide which", result.anomalous_deepest_site, result.anomalous_deepest_sigma); if (maps_superseded) - logger.Info("Model validation: no maps made here - the validation in the model's setting writes them"); - else - logger.Info("Model validation: wrote {}_2fofc.ccp4, {}_fofc.ccp4{}, {}_maps.mtz", - output_prefix, output_prefix, - result.anomalous_sites.empty() ? "" : ", " + output_prefix + "_anom.ccp4", - output_prefix); + logger.Info("Model validation: no maps made here - the validation in the model's setting makes them"); if (result.fit_tested && !result.model_fits) logger.Info("Model validation: the model was rejected, so nothing downstream moved - the " "reflection files are byte for byte the ones a run with no model would have " @@ -1917,7 +1897,6 @@ ModelValidationGPUPlan PlanModelValidationGPU(const ModelValidationGPUNeed &need ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, const RugnuxModel &model, - const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group, bool probe_indexing_ambiguity, @@ -1932,7 +1911,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // the validation reports why it did not finish, and a validation that did not finish decides nothing, // so the reflection files are those of a run without a model. try { - return Validate(merged, cell, model, output_prefix, logger, data_space_group, + return Validate(merged, cell, model, logger, data_space_group, probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, true, placed); } catch (const RigidBodyGPUFailure &e) { @@ -1945,7 +1924,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return failed; } #else - return Validate(merged, cell, model, output_prefix, logger, data_space_group, probe_indexing_ambiguity, + return Validate(merged, cell, model, logger, data_space_group, probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false, placed); #endif } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index 594c3fb65..4dac7e8e8 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -11,15 +11,25 @@ #include #include +#include "gemmi/grid.hpp" #include "gemmi/symmetry.hpp" #include "../../common/Reflection.h" // MergedReflection #include "../../common/UnitCell.h" +#include "../../common/MtzTable.h" #include "RugnuxModel.h" class Logger; namespace gemmi { struct Structure; } +// The maps of a validation: their coefficients, as _maps.mtz holds them, and the maps themselves, +// as the .ccp4 files hold them, where the caller asked for them (ModelValidationSchedule::keep_map_grids). +// The anomalous map is there only where the merge kept the Bijvoet split. +struct ModelMaps { + MtzTable coefficients; + std::optional> map_2fofc, map_fofc, map_anom; +}; + // Result of validating merged intensities against an atomic model. struct ModelValidationResult { bool ok = false; @@ -58,7 +68,11 @@ struct ModelValidationResult { double radial_misfit = NAN; int shell_scale_bins = 0; double mean_atom_density_sigma = 0.0; // mean 2mFo-DFc value at atom centres, in sigma - std::string maps_prefix; // where the .ccp4 / _maps.mtz were written + // The maps, where this validation made them; the caller takes them over (null after that). + std::shared_ptr maps; + // Where the caller puts the maps beside the reflection files: "" for the answer's, "_mr" for the + // model's choice. Set by the caller; the report names . + std::string maps_suffix; // The rigid-body placement of the model against these data, where it was committed. The step is // taken only when it lowers R-free, so an unmoved model means the data said no, not that nothing @@ -211,8 +225,8 @@ struct ModelValidationResult { // Given merged intensities and an atomic model, scale the model structure factors (with an // optimized flat bulk-solvent contribution and an overall anisotropic B) to the WORKING observed -// amplitudes, then report R-work / R-free and write 2Fo-Fc and Fo-Fc electron-density maps -// (CCP4) plus an MTZ of map coefficients next to output_prefix. Nothing the reported R-free is +// amplitudes, then report R-work / R-free and make 2Fo-Fc and Fo-Fc electron-density maps and +// their coefficients (ModelValidationResult::maps). Nothing the reported R-free is // computed against was fitted on the free reflections: the scale, like the rigid-body placement // below, sees the working set only and is then applied to all of them. // @@ -266,10 +280,10 @@ struct ModelFrameForecast { struct ModelValidationSchedule { // Called once, before the null, with the forecast above. std::function on_forecast; - // Where valid, waited on before the first file is written: false returns without writing any. - std::shared_future write_gate; + // Where valid, waited on before the maps are made: false returns without making any. + std::shared_future maps_gate; // Asked once the verdict, the setting and the indexing are decided, before any map is made: true - // where another validation will write this one's maps and map coefficients, which are then not made. + // where another validation's maps and map coefficients are the ones kept, so this one's are not made. std::function maps_superseded; // The null that says whether the model fits. False where the verdict is another validation's // (KeepModelVerdict) or was reached elsewhere (a molecular-replacement placement, scored for its maps). @@ -280,6 +294,8 @@ struct ModelValidationSchedule { // GPU memory another validation running beside this one holds there. It limits how many engines // this one takes and never what runs on the GPU, so no number depends on whether the other runs. size_t gpu_beside_bytes = 0; + // Keep the maps themselves (ModelMaps), not only their coefficients: rugnux writes them as .ccp4. + bool keep_map_grids = false; }; // What a model validation needs of its card, part by part: the bytes of one engine of each, how many @@ -330,7 +346,6 @@ ModelValidationGPUPlan PlanModelValidationGPU(const ModelValidationGPUNeed &need ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, const RugnuxModel &model, - const std::string &output_prefix, Logger &logger, const gemmi::SpaceGroup *data_space_group = nullptr, bool probe_indexing_ambiguity = true, diff --git a/receiver/JFJochReceiverRugnux.cpp b/receiver/JFJochReceiverRugnux.cpp index ec3084643..b4e9daf2a 100644 --- a/receiver/JFJochReceiverRugnux.cpp +++ b/receiver/JFJochReceiverRugnux.cpp @@ -246,7 +246,7 @@ void JFJochReceiverRugnux::RugnuxThread() { parallel_detail::WorkerPool::Instance().SetNice(10); try { const auto dataset = reader.GetDataset(); - const RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, scratch_dir + "/p", model); + const RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, model); LiveObserver observer(plots, scan_result, *az_int_mapping, live_status, logger); ProcessResult result; { @@ -288,7 +288,7 @@ void JFJochReceiverRugnux::ProbeThread() { const int64_t frames = std::lround(PROBE_DEG / std::abs(gonio->GetIncrement_deg())); if (2 * frames > static_cast(reader.GetNumberOfImages())) return; // a sweep this short is decided on soon enough by the run itself - RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, "", nullptr); + RugnuxLiveSetup setup = MakeRugnuxLiveSetup(*dataset, rugnux_nthreads, nullptr); setup.config.end_image = static_cast(frames); setup.config.rotation_postrefine_geometry = false; setup.config.run_scaling = false; diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 83444105f..448975f01 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -31,13 +31,18 @@ ADD_LIBRARY(JFJochRugnux STATIC BrokerFeed.h ) -TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter JFJochModelWriter JFJochPreview gemmi fftw3f +TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochPreview gemmi fftw3f httplib::httplib) +# 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) +TARGET_LINK_LIBRARIES(JFJochRugnuxOutput JFJochRugnux JFJochWriter JFJochModelWriter) + # rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the # default), azint, scale (re-scale/merge stored reflections) or calibration (powder-ring geometry). ADD_EXECUTABLE(rugnux rugnux_cli.cpp RugnuxConsole.cpp RugnuxConsole.h) -TARGET_LINK_LIBRARIES(rugnux JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter) +TARGET_LINK_LIBRARIES(rugnux JFJochRugnuxOutput JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS rugnux RUNTIME COMPONENT rugnux) # macOS: sign the installed binary with the Developer ID (hardened runtime, secure timestamp) so the # .tgz can be notarized - see JFJOCH_MACOS_SIGNING_IDENTITY in the top-level CMakeLists.txt. The diff --git a/rugnux/DiagnosticOutput.cpp b/rugnux/DiagnosticOutput.cpp index 96557552d..7083bccaa 100644 --- a/rugnux/DiagnosticOutput.cpp +++ b/rugnux/DiagnosticOutput.cpp @@ -5,7 +5,6 @@ #include #include -#include #include @@ -103,17 +102,13 @@ std::string RenderDetectorDiagnostic(const DiffractionExperiment &experiment, return WriteJPEGToMem(rgb_image, 85); } -void WritePerImagePlot(const std::string &output_prefix, const PerImagePlot &plot, - const std::optional &goniometer, size_t image_count) { - const std::string path = output_prefix + "_plot.txt"; - std::ofstream file(path, std::ofstream::out | std::ofstream::trunc); - if (!file) - throw JFJochException(JFJochExceptionCategory::FileWriteError, - "Cannot open " + path + " for writing"); +std::string PerImagePlotText(const PerImagePlot &plot) { + const std::optional &goniometer = plot.goniometer; + std::string file; // One comment line, its names over their own columns, and nothing else a plotting program has // to be told to skip. The '#' takes the first character of the image-number field. - file << "#" << fmt::format("{:>5} {:>10} {:>12} {:>12} {:>8} {:>12} {:>13} {:>12} {:>8}" + file += "#" + fmt::format("{:>5} {:>10} {:>12} {:>12} {:>8} {:>12} {:>13} {:>12} {:>8}" " {:>7} {:>11}\n", "image", "angle_deg", "bkg", "resolution_A", "spots", "scale_G", "sigma_M_deg", "cc_to_merge", "cc_n", "merged", "disposition"); @@ -125,7 +120,7 @@ void WritePerImagePlot(const std::string &output_prefix, const PerImagePlot &plo || std::any_of(plot.scale_g.begin(), plot.scale_g.end(), [](float g) { return std::isfinite(g); }); - for (size_t i = 0; i < image_count; i++) { + for (size_t i = 0; i < plot.image_count; i++) { const float angle = goniometer ? goniometer->GetAngle_deg(static_cast(i)) + goniometer->GetWedge_deg() / 2.0f @@ -152,11 +147,12 @@ void WritePerImagePlot(const std::string &output_prefix, const PerImagePlot &plo } else if (merge_ran && std::isfinite(scale_g)) disposition = 1.0f; - file << fmt::format("{:>6} {:>10.4f} {:>12.4g} {:>12.4g} {:>8.0f} {:>12.5g} {:>13.4f}" + file += fmt::format("{:>6} {:>10.4f} {:>12.4g} {:>12.4g} {:>8.0f} {:>12.5g} {:>13.4f}" " {:>12.4f} {:>8} {:>7} {:>11}\n", i, angle, Element(plot.bkg_estimate, i), Element(plot.resolution_A, i), Element(plot.spot_count, i), scale_g, Element(plot.mosaicity_deg, i), Element(plot.cc_to_merge, i), Element(plot.cc_n, i), merge_ran ? std::to_string(merged) : "nan", disposition); } + return file; } diff --git a/rugnux/DiagnosticOutput.h b/rugnux/DiagnosticOutput.h index 0566b2da5..04d41d4e2 100644 --- a/rugnux/DiagnosticOutput.h +++ b/rugnux/DiagnosticOutput.h @@ -40,9 +40,10 @@ struct PerImagePlot { std::vector cc_to_merge; std::vector cc_n; std::vector frame_disposition; + std::optional goniometer; // the angle column; none: nan + size_t image_count = 0; // the rows }; -// Writes _plot.txt: one '#' legend line naming the columns, then one whitespace- -// separated row per image, so `plot '' using 1:4` works with nothing in between. -void WritePerImagePlot(const std::string &output_prefix, const PerImagePlot &plot, - const std::optional &goniometer, size_t image_count); +// The text of _plot.txt: one '#' legend line naming the columns, then one whitespace-separated +// row per image, so `plot '' using 1:4` works with nothing in between. +std::string PerImagePlotText(const PerImagePlot &plot); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 9e56e7da6..61d8a5203 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include #include @@ -1459,7 +1458,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, domains ? ", or spots of the further lattice domains found land on those positions" : "", result.doubled_cell ? fmt::format("Both cells were merged; the other one is in {}.mtz - compare them in refinement", - result.doubled_cell->other_prefix) + output_prefix + result.doubled_cell->other_suffix) : fmt::format("The doubled cell {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} was not merged", sc.doubled_cell.a, sc.doubled_cell.b, sc.doubled_cell.c, sc.doubled_cell.alpha, sc.doubled_cell.beta, sc.doubled_cell.gamma))); @@ -1467,7 +1466,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const auto &h = *result.doubled_cell; Add(s, KeyBool("DOUBLED_CELL_ADOPTED", h.adopted)); Add(s, KeyText("DOUBLED_CELL_DECISION", h.decision)); - Add(s, KeyText("DOUBLED_CELL_OTHER_FILE", h.other_prefix + ".mtz")); + Add(s, KeyText("DOUBLED_CELL_OTHER_FILE", output_prefix + h.other_suffix + ".mtz")); Add(s, KeyText("DOUBLED_CELL_OTHER_CELL", CellString(h.other_cell))); Add(s, KeyText("DOUBLED_CELL_OTHER_SPACE_GROUP", h.other_space_group)); Add(s, KeyReal("DOUBLED_CELL_OTHER_D_MIN", h.other_d_min_A, "{:.2f}")); @@ -1932,7 +1931,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (result.twin_subgroup_merge) { const auto &t = *result.twin_subgroup_merge; - Add(s, KeyText("TWIN_SUBGROUP_FILE", t.file_prefix + ".mtz")); + Add(s, KeyText("TWIN_SUBGROUP_FILE", output_prefix + t.file_suffix + ".mtz")); Add(s, KeyText("TWIN_SUBGROUP_SPACE_GROUP", t.space_group)); Add(s, KeyText("TWIN_SUBGROUP_LAW", t.twin_law)); Add(s, KeyReal("TWIN_SUBGROUP_FRACTION", t.twin_fraction, "{:.3f}")); @@ -1952,7 +1951,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " as _pdbx_reflns_twin). Refine against it with that twin law where the model prefers it.\n" " Nothing in this report or in the answer's files depends on it.{}", t.adopted_space_group, t.reason, t.space_group, t.adopted_space_group, t.space_group, - t.space_group, t.file_prefix, t.twin_law, t.twin_fraction, + t.space_group, output_prefix + t.file_suffix, t.twin_law, t.twin_fraction, t.other_subgroups.empty() ? std::string() : "\n Other subgroups the search confirmed, each a twin hypothesis of\n" " its own (re-run with -S to merge in one): see\n" @@ -2441,12 +2440,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, mv.indexing_margin_null_sd), true)); Add(s, KeyReal("MODEL_INDEXING_MARGIN_SIGMA", mv.indexing_margin_sigma, "{:+.2f}", true)); } - if (!mv.maps_prefix.empty()) - Add(s, KeyText("MAPS_PREFIX", mv.maps_prefix)); + if (mv.ok && !output_prefix.empty()) + Add(s, KeyText("MAPS_PREFIX", output_prefix + mv.maps_suffix)); // The model's choice of another group or of the twin hypothesis' subgroup, by either route: files of // their own, never the answer's. - if (!result.model_prior_file_prefix.empty()) { - Add(s, KeyText("MODEL_CHOICE_FILE_PREFIX", result.model_prior_file_prefix)); + if (!result.model_prior_file_suffix.empty()) { + Add(s, KeyText("MODEL_CHOICE_FILE_PREFIX", output_prefix + result.model_prior_file_suffix)); Add(s, KeyText("MODEL_CHOICE_SPACE_GROUP", result.model_prior_file_space_group)); } // The twin hypothesis (section 4): the model scored against that merge too, and which of @@ -2474,7 +2473,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " fraction a scanned on the working set. The maps and the placed model come from {},\n" " the merge the model fits better on R-free ({:.4f} against {:.4f}). This chooses the\n" " maps only: the answer's reflection files and every number above are the answer's.", - t.space_group, t.file_prefix, t.twin_law, + t.space_group, output_prefix + t.file_suffix, t.twin_law, t.maps_from_this_merge ? "the merge in " + t.space_group : "the answer's merge", tv.twin_target ? std::min(tv.r_free, tv.r_free_twin) : tv.r_free, mv.r_free))); } @@ -2704,7 +2703,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n" " chooses is written as {}.mtz/.cif with the placed model.", mp.d_max_search, mp.d_min_search, - result.model_prior_file_prefix.empty() ? std::string("_mr") : result.model_prior_file_prefix))); + result.model_prior_file_suffix.empty() ? std::string("_mr") + : output_prefix + result.model_prior_file_suffix))); Add(s, Prose(fmt::format( " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n" " R-free {:.2f}, null {:.2f}, PAE domains {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, @@ -3136,26 +3136,3 @@ std::string RenderResultReport(const std::string &output_prefix, BuildReportDocument(output_prefix, input_file, experiment, result, provenance), provenance.developer); } - -ReportDocument WriteResultReport(const std::string &output_prefix, - const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - Logger &logger, - const RunProvenance &provenance) { - ReportDocument doc = BuildReportDocument(output_prefix, input_file, experiment, result, provenance); - if (output_prefix.empty()) - return doc; // "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 << RenderReportText(doc, provenance.developer); - } catch (const std::exception &e) { - logger.Warning("Could not write the results report {}: {}", filename, e.what()); - } - return doc; -} diff --git a/rugnux/ResultReport.h b/rugnux/ResultReport.h index d27bd3eb5..bdbbe1b72 100644 --- a/rugnux/ResultReport.h +++ b/rugnux/ResultReport.h @@ -72,14 +72,3 @@ std::string RenderResultReport(const std::string &output_prefix, const DiffractionExperiment &experiment, const ProcessResult &result, const RunProvenance &provenance = {}); - -// Renders and writes _report.txt, and returns the document it rendered (the CLI's -// closing summary reads it). Writes 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. -ReportDocument WriteResultReport(const std::string &output_prefix, - const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - Logger &logger, - const RunProvenance &provenance = {}); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 23c1dfedb..899a7c4de 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -40,7 +38,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" @@ -752,7 +749,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const bool calibration = (config_.mode == ProcessMode::Calibration); const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots; const bool per_image_analysis = full || calibration; - const bool write_files = write_output && !config_.output_prefix.empty(); + // The output pass makes what the run returns to its caller (ProcessResult::outputs). + if (write_output && !geometry_prepass) + result.outputs = std::make_shared(); // Set below the first merge, once the pass can tell that RunAllPasses is going to throw it away and // run again - see there for what that turns off and why it is knowable at that point. bool superseded = false; @@ -762,9 +761,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.BitDepthImage(32).PixelSigned(true); experiment_.Mode(DetectorMode::Standard); experiment_.OverwriteExistingFiles(true); - // Offline processing: the output prefix is a local path the operator chose (rugnux -o or the - // viewer's "next to the input file"), so use the trusted setter that permits an absolute path. - experiment_.FilePrefixTrusted(config_.output_prefix.empty() ? "output" : config_.output_prefix); experiment_.SetFileWriterFormat(FileWriterFormat::NXmxLegacy); experiment_.ImagesPerTrigger(images_to_process); experiment_.NumTriggers(1); @@ -838,7 +834,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b JFJochReceiverPlots plots; plots.Setup(experiment_, mapping); - // Output file (NXmxIntegrated master that links back to the original images). + // The per-image analysis' start message (an NXmxIntegrated master that links back to the original + // images); the image sink names the file. StartMessage start_message; experiment_.FillMessage(start_message); start_message.arm_date = dataset->arm_date; @@ -889,7 +886,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing - // but the process file reads them. Copy them where a writer is actually built rather than here: + // but the process file reads them. Copy them where the start message is handed on rather than here: // the default merging path writes no process file at all. auto fill_writer_maps = [&] { start_message.pixel_mask["default"] = pixel_mask_.GetMask(experiment_); @@ -897,18 +894,18 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b start_message.roi_map = experiment_.ExportROIMap(); }; - std::unique_ptr writer; + // The output pass streams its per-image analysis to the caller's image sink, where there is one. + RugnuxImageSink *const sink = write_output && !geometry_prepass ? config_.image_sink : nullptr; std::unique_ptr writer_queue; // Whether anything is in place to break the indexing ambiguity, for the warning further down: a // reference MTZ, or (serial stills) the reference computed from a model just below. bool ambiguity_reference = !config_.reference_data.empty(); - if (write_files && config_.write_process_h5) { + if (sink) { fill_writer_maps(); - writer = std::make_unique(start_message, /*check_overwrite_at_start=*/true, - /*trusted_path=*/true); + sink->Start(start_message, output_label_); // Deep enough that a worker never waits for the writer in the normal case, shallow enough that // the backlog is bounded by the worker count rather than by the length of the run. - writer_queue = std::make_unique(*writer, 4 * std::max(1, config_.nthreads)); + writer_queue = std::make_unique(*sink, 4 * std::max(1, config_.nthreads)); } const char *mode_name = full ? "full analysis" @@ -1002,26 +999,23 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const int loop_end = geometry_prepass || prepass_frames_only_ ? prepass_end : end_image; // Declared here rather than where they are first filled, so that the phases below can share them - // through PipelineLocals: the bytes the image loop read, the end message (also written to the file), - // and the unmerged MTZ when it is built beside the reflection files (see there) with the - // experiment it is built from - which it points into for its space group, so it lives until the - // file is written. + // through PipelineLocals: the bytes the image loop read, the end message (also handed to the sink), + // and the unmerged MTZ when it is built beside the reflection files (see there). std::atomic total_uncompressed_bytes = 0; EndMessage end_msg; - std::optional unmerged_experiment; - std::future unmerged_mtz; + std::future unmerged_mtz; PipelineLocals p{ .logger = logger, .result = result, .observer = observer, .geometry_prepass = geometry_prepass, .full = full, .calibration = calibration, .calibration_spots = calibration_spots, .per_image_analysis = per_image_analysis, - .write_files = write_files, .start_image = start_image, .end_image = end_image, + .write_output = write_output, .start_image = start_image, .end_image = end_image, .images_to_process = images_to_process, .prepass_images = prepass_images, .loop_end = loop_end, .dataset = dataset, .mapping = mapping, .plots = plots, - .writer = writer, .writer_queue = writer_queue, .indexer_pool = indexer_pool, .indexer = indexer, + .sink = sink, .writer_queue = writer_queue, .indexer_pool = indexer_pool, .indexer = indexer, .indexing_bound_A = indexing_bound_A, .census_future = census_future, .superseded = superseded, .ambiguity_reference = ambiguity_reference, .start_time = start_time, .total_uncompressed_bytes = total_uncompressed_bytes, .end_msg = end_msg, - .unmerged_experiment = unmerged_experiment, .unmerged_mtz = unmerged_mtz}; + .unmerged_mtz = unmerged_mtz}; // First pass of two-pass rotation indexing (full analysis only). if (full && force_rotation_result_.has_value()) { @@ -1098,47 +1092,42 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b return result; } - TimingMark output_mark(logger, "output files"); + TimingMark output_mark(logger, "outputs"); // Unmerged observations. After the merge, because that is where the space group is settled and // the file has to be written in it, but from the integrated observations rather than the merged // ones: the partiality and the per-image scale are left for the reading program, which fits a - // scale model of its own. write_files, not write_output - without an output prefix there is - // nowhere to put it, and the name would come out as a bare "_unmerged.mtz". - if (full && !cancelled_ && write_files && !geometry_prepass && !superseded && result.consensus_cell) { + // scale model of its own. + if (full && !cancelled_ && result.outputs && !superseded && result.consensus_cell) { if (config_.export_unmerged) { - if (observer) observer->OnPhase("Writing unmerged reflections"); - const std::string path = config_.output_prefix + "_unmerged.mtz"; - if (unmerged_mtz.valid()) { - gemmi::Mtz mtz = unmerged_mtz.get(); - SetUnmergedMtzMosaicity(mtz, indexer->GetIntegrationOutcome()); - mtz.write_to_file(path); - } else - WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, - experiment_, true, config_.provenance, path, config_.nthreads); - logger.Info("Unmerged observations written to {}", path); + if (observer) observer->OnPhase("Unmerged reflections"); + MtzTable mtz = unmerged_mtz.valid() + ? unmerged_mtz.get() + : UnmergedMtzTable(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, true, config_.provenance, config_.nthreads); + SetUnmergedMtzMosaicity(mtz, indexer->GetIntegrationOutcome()); + result.outputs->unmerged = std::move(mtz); } if (config_.export_unmerged_partials) { - if (observer) observer->OnPhase("Writing unmerged reflections"); - const std::string path = config_.output_prefix + "_unmerged_partials.mtz"; - WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell, - experiment_, false, config_.provenance, path, config_.nthreads); - logger.Info("Unmerged partials written to {}", path); + if (observer) observer->OnPhase("Unmerged reflections"); + MtzTable mtz = UnmergedMtzTable(indexer->GetIntegrationOutcome(), *result.consensus_cell, + experiment_, false, config_.provenance, config_.nthreads); + SetUnmergedMtzMosaicity(mtz, indexer->GetIntegrationOutcome()); + result.outputs->unmerged_partials = std::move(mtz); } } // When merging is on, the _process.h5 is skipped because the merged reflections are the wanted // output and the per-image file is large. If nothing indexed there are no merged reflections, so - // the run would finish successfully having written no file at all. Write one now: the per-image - // messages went past while there was no writer, but the end message still carries the per-image + // the run would finish successfully having written no file at all. Hand the caller one: the + // per-image messages went past with no sink, but the end message still carries the per-image // scalars (spot counts, background, resolution estimate) and the dataset metadata - which is what // anyone asking "why did nothing index" needs, and it is small precisely because nothing indexed. - if (!writer && full && write_files && !cancelled_ && result.indexing_rate.value_or(0.0f) <= 0.0f) { - logger.Warning("No image indexed, so there are no merged reflections to write - writing " - "{}_process.h5 with the per-image analysis instead", config_.output_prefix); + if (!sink && full && result.outputs && !cancelled_ && result.indexing_rate.value_or(0.0f) <= 0.0f) { + logger.Warning("No image indexed, so there are no merged reflections - the per-image analysis is " + "kept for a _process.h5 instead"); fill_writer_maps(); - writer = std::make_unique(start_message, /*check_overwrite_at_start=*/true, - /*trusted_path=*/true); + result.outputs->process_start = start_message; } // The per-frame sigma_M the merge wrote back onto the outcomes (the _plots.txt column), for the report. @@ -1154,22 +1143,13 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.mosaicity_deg_p90 = mos[mos.size() * 9 / 10]; } } - // The two files a person checks the run by eye with. Not on the geometry pre-pass or a pass that - // is thrown away, which would otherwise leave a stale _01 copy beside the real one. The table is - // written on a stills run too - it is the only per-image scaling table there is, and the columns + // The two files a person checks the run by eye with. Not on a pass that is thrown away. The table + // is made on a stills run too - it is the only per-image scaling table there is, and the columns // a stills run cannot fill come out as nan; the picture is rotation-only at its source, since // only the rotation pre-scan builds the projection it draws. - if (full && write_files && !cancelled_ && !geometry_prepass && !superseded) { - const std::string jpeg = detector_diagnostic_jpeg_.valid() ? detector_diagnostic_jpeg_.get() : std::string(); - if (!jpeg.empty()) { - const std::string path = config_.output_prefix + "_detector.jpg"; - std::ofstream f(path, std::ofstream::out | std::ofstream::trunc | std::ofstream::binary); - if (f) { - f.write(jpeg.data(), static_cast(jpeg.size())); - logger.Info("Detector diagnostic image written to {}", path); - } else - logger.Warning("Cannot open {} for writing", path); - } + if (full && result.outputs && !cancelled_ && !superseded) { + if (detector_diagnostic_jpeg_.valid()) + result.outputs->detector_jpeg = detector_diagnostic_jpeg_.get(); PerImagePlot plot; plots.GetPlotRaw(plot.bkg_estimate, PlotType::BkgEstimate, ""); @@ -1188,17 +1168,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } if (result.has_merge_statistics) plot.frame_disposition = result.merge_statistics.sweep_quality.frame_disposition; - WritePerImagePlot(config_.output_prefix, plot, experiment_.GetGoniometer(), - static_cast(images_to_process)); - logger.Info("Per-image table written to {}_plot.txt", config_.output_prefix); + plot.goniometer = experiment_.GetGoniometer(); + plot.image_count = static_cast(images_to_process); + result.outputs->plot = std::move(plot); } - if (writer) { - writer->WriteHDF5(end_msg); - writer->Finalize(); - result.written_master_path = config_.output_prefix + "_process.h5"; - } + if (sink) + sink->End(end_msg); result.end_message = end_msg; + // A pass that is to be run again makes nothing. + if (superseded) + result.outputs.reset(); output_mark.Stop(); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index c70acc7da..505e637ee 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -44,6 +44,7 @@ #include "../image_analysis/WriteReflections.h" // OutputProvenance #include "../image_analysis/bragg_integration/SpotWidth.h" // FluxCurve #include "../image_analysis/bragg_integration/SpotFootprint.h" // FootprintSpot +#include "DiagnosticOutput.h" // PerImagePlot class JFJochReader; class JFJochHDF5Reader; @@ -51,17 +52,31 @@ struct JFJochReaderRawImage; struct PipelineLocals; // RugnuxPipeline.h, private to the implementation // Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (all of its --mode -// values that read images) and (later) the viewer. The full workflow lives here, not in the CLIs: +// values that read images), the viewer and the broker. The full workflow lives here, not in the CLIs: // setup, an optional two-pass rotation-indexing pre-pass, a parallel per-image loop (std::thread), -// an optional scaling/merging post-pass, and the _process.h5 output. The detector geometry and all -// algorithm settings are configured on the DiffractionExperiment by the caller; ProcessConfig only -// carries run control. Cancellable from any thread (e.g. SIGINT or a GUI button) via Cancel(). +// an optional scaling/merging post-pass. It reads and writes no file: what it makes is returned in +// ProcessResult (RugnuxOutputs) and handed to RugnuxImageSink, and the caller writes it +// (RugnuxOutputFiles.h). The detector geometry and all algorithm settings are configured on the +// DiffractionExperiment by the caller; ProcessConfig only carries run control. Cancellable from any +// thread (e.g. SIGINT or a GUI button) via Cancel(). enum class ProcessMode { AzimuthalIntegration, // preprocess + azimuthal integration only (rugnux --mode=azint) FullAnalysis, // spot finding + indexing + refinement + integration (rugnux --mode=mx) Calibration // detector geometry from powder rings (rugnux --mode=calibration) }; +// Where the per-image analysis of an output pass goes as it is made: rugnux and the viewer put it in +// _process.h5 (ProcessFileSink, RugnuxOutputFiles.h). Called from one thread at a time. +class RugnuxImageSink { +public: + virtual ~RugnuxImageSink() = default; + // A pass starts. `label` tells its files apart from the answer's: empty for every pass but the + // doubled cell's, which is "_doubled" (ProcessResult::DoubledCell). + virtual void Start(const StartMessage &start, const std::string &label) = 0; + virtual void Write(const DataMessage &msg) = 0; + virtual void End(const EndMessage &end) = 0; +}; + struct ProcessConfig { ProcessMode mode = ProcessMode::FullAnalysis; @@ -73,20 +88,23 @@ struct ProcessConfig { // The per-image loop is capped in that case (see RunPipeline); an explicit -N never is. bool nthreads_auto = false; - // Output prefix for the _process.h5 (and scaled reflections). Empty => process without writing. - std::string output_prefix; - // What the reflection files record about where they came from (WriteReflections). + // What the reflection files record about where they came from (the MTZ HISTORY, WriteReflections). OutputProvenance provenance; - // Write the per-image _process.h5. Defaults true; rugnux turns it off when merging - // (the .mtz/.cif is the wanted output and the h5 is large) unless --write-process-h5 is given. - bool write_process_h5 = true; + // Where the per-image analysis of each output pass goes (the _process.h5); none: it is not kept. + // rugnux gives one when not merging, or with --write-process-h5. + RugnuxImageSink *image_sink = nullptr; - // Write the merged reflections (.mtz/.cif). Defaults true (the CLI always writes them when - // merging); the viewer can turn them off independently of the _process.h5 (e.g. a quick run that - // only wants the ISa/R-meas numbers). Only relevant to FullAnalysis with scaling. + // Make the merged reflections' outputs (RugnuxOutputs::merged and the twin hypothesis' and finalist + // ledger's merges). Defaults true (the CLI always writes them when merging); the viewer can turn them + // off (e.g. a quick run that only wants the ISa/R-meas numbers). Only relevant to FullAnalysis with + // scaling. bool write_merged = true; + // Keep the model validation's maps themselves (ModelMaps), not only their coefficients: rugnux writes + // them as .ccp4. The broker never asks for them. + bool keep_map_grids = 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). @@ -191,15 +209,16 @@ struct ProcessConfig { std::string reference_path; std::string reference_column; - // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). - std::string observation_dump_path; + // Diagnostic: if set, handed the unmerged fulls of every -P rot3d combine as text, for comparison + // with XDS (--dump-observations). + std::function observation_dump; - // Write the integrated observations as _unmerged.mtz, for aimless / pointless / careless - // and anything else that scales unmerged data itself - each rocking event summed into one full. - // On by default (--no-export-unmerged turns it off): it is how a run's observations reach the - // rest of the field, and a user who does not know the flag exists should still get the file. - // It is the largest file a run produces, so a batch that only wants the merged numbers turns it - // off. Needs an output prefix. + // Make the integrated observations as _unmerged.mtz (RugnuxOutputs::unmerged), for aimless / + // pointless / careless and anything else that scales unmerged data itself - each rocking event summed + // into one full. On by default (--no-export-unmerged turns it off): it is how a run's observations + // reach the rest of the field, and a user who does not know the flag exists should still get the + // file. It is the largest file a run produces, so a batch that only wants the merged numbers turns + // it off. bool export_unmerged = true; // The same as _unmerged_partials.mtz, with each partial written as its own row for the @@ -322,6 +341,48 @@ struct SearchMergeSignal { int64_t inner_reflections = 0; // ... of this many reflections }; +// The reflection files of one merge, .mtz and .cif (and the answer's .hkl): the MTZ as its +// table, and what the mmCIF and the HKLF-4 file are written from (WriteReflections.h). +struct MergedOutput { + std::string suffix; // "", "_twin_" or "_mr" + MtzTable mtz; + std::vector reflections; + UnitCell cell{}; + DiffractionExperiment experiment; // in the group the files are written in + MergeStatistics statistics; + ErrorModelReport error_model; + TwinningAnalysisResult twinning; // a twin hypothesis' law and fraction, for its mmCIF + std::vector scaled_fulls; // the .hkl's records, where the merge made them (rotation) +}; + +// An atomic model as placed against the data (_model.cif, WritePlacedModel): `cell` and +// `space_group` are those of the reflection files beside it. +struct PlacedModelOutput { + std::shared_ptr model; + UnitCell cell{}; + gemmi::SpaceGroup space_group; +}; + +// What a run makes for its caller to keep - the files rugnux writes, each named here after the file it +// goes to at the caller's prefix. Made by the output pass, the one whose result the run returns. +struct RugnuxOutputs { + std::optional merged; // .mtz/.cif/.hkl + std::optional twin; // the twin hypothesis' merge + std::optional model_prior; // the model's choice, _mr.mtz/.cif + std::shared_ptr maps; // _maps.mtz, _2fofc/_fofc/_anom.ccp4 + std::shared_ptr model_prior_maps; // the same at _mr + std::optional placed_model; // _model.cif/.pdb + std::optional model_prior_placed_model; // _mr_model.cif/.pdb + std::optional unmerged; // _unmerged.mtz + std::optional unmerged_partials; // _unmerged_partials.mtz + std::string detector_jpeg; // _detector.jpg; empty: none drawn + std::optional plot; // _plot.txt + // A run that indexed nothing, with no per-image analysis streamed (ProcessConfig::image_sink): the + // start message of a _process.h5 holding the end message (ProcessResult::end_message), so + // there is a file that says why nothing indexed. + std::optional process_start; +}; + struct ProcessResult { bool cancelled = false; uint64_t images_processed = 0; @@ -336,6 +397,8 @@ struct ProcessResult { // The end message of this pass, as the _process.h5 records it: the lattice, refined geometry, // reindexing and the run means - what the broker puts into its own when it runs rugnux live. EndMessage end_message; + // What this pass made for the caller to keep - the output pass's only; null on every other pass. + std::shared_ptr outputs; // Workers the per-image loop actually ran, which is not config.nthreads once the loop is capped. int image_loop_threads = 0; double frame_rate_hz = 0.0; @@ -392,7 +455,6 @@ struct ProcessResult { // every other class, or where the free refinement is not available. std::optional lattice_angle_departure_deg; MeanProcessingTime mean_processing_time{}; - std::optional written_master_path; std::string merge_statistics_text; // populated when scaling/merging ran // Structured merge statistics (per-shell + overall), populated when scaling/merging ran. ISa is @@ -569,7 +631,7 @@ struct ProcessResult { std::string reason; // why the hypothesis was written std::string adopted_space_group; // the group of the answer, as written std::string space_group; // the subgroup, as written (standard setting) - std::string file_prefix; // _twin_ + std::string file_suffix; // _twin_, or _mr: the files are std::string twin_law; // hkl triplet, on the subgroup file's axes double twin_fraction = NAN; // EstimateTwinLawFraction: a lower bound int twin_fraction_pairs = 0; @@ -591,7 +653,8 @@ struct ProcessResult { struct DoubledCell { bool adopted = false; // the answer is the doubled cell, the other one the probed cell std::string decision; // why, in the words the run reports it in - std::string other_prefix; // _halfcell or _doubled + std::string other_suffix; // _halfcell or _doubled: the other cell's files are + std::shared_ptr other_outputs; UnitCell other_cell{}; std::string other_space_group; double other_d_min_A = NAN, other_cc_half = NAN, other_r_free = NAN; @@ -611,9 +674,9 @@ struct ProcessResult { std::string model_placement; // What molecular replacement found, where it ran; failure_reason says why where it could not. std::optional model_prior; - // The file prefix the model's choice of another group or the twin hypothesis' subgroup was written under - // (_mr), and that group; empty where nothing was. - std::string model_prior_file_prefix; + // The file suffix the model's choice of another group or the twin hypothesis' subgroup is written under + // (_mr), and that group; empty where nothing was. + std::string model_prior_file_suffix; std::string model_prior_file_space_group; // Space group used for the final re-scale/merge and written to the master file: determined by the @@ -883,8 +946,11 @@ class Rugnux { std::optional force_rotation_result_; // The indexing result the last full pass integrated on, which the doubled-cell hypothesis doubles. std::optional last_rotation_result_; - // Keep the written merged reflections on the result (ProcessResult::merged_reflections). + // Keep the merged reflections on the result (ProcessResult::merged_reflections). bool keep_merged_reflections_ = false; + // What the image sink is told the output pass is (RugnuxImageSink::Start): "_doubled" for the doubled + // cell's pass, empty otherwise. + std::string output_label_; // The space group the USER fixed (-S on the CLI, the settings panel in the viewer), captured before // anything can clear it - a pass that determines a group of its own writes it onto the experiment, @@ -1270,7 +1336,8 @@ class Rugnux { RugnuxObserver *observer); // One full run of the pipeline (index + integrate + optional scale/merge/space-group + output). - // write_output == false suppresses all file output (the geometry pre-pass writes nothing). + // write_output == false makes no outputs (ProcessResult::outputs) and streams nothing to the image sink; + // nor does the geometry pre-pass. // geometry_prepass == true stops after post-refining the detector geometry from the first integration // and applies it (+ the refined lattice) to experiment_, so the next pass re-integrates with it. ProcessResult RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass); @@ -1342,7 +1409,6 @@ class Rugnux { }; // What the refined passes start from. struct RefinedPassStart { - std::string base_prefix; std::optional gonio_snapshot; std::array header_geometry{}; RotationScaleOutcome scale_walk; // where the walk has ended before they start @@ -1351,8 +1417,7 @@ class Rugnux { // the rotation-scale walk and it overturned them. // After the refined passes: where the supercell probe's class stands clear, integrate and merge the // doubled cell too, and return whichever of the two is the answer (ProcessResult::doubled_cell). - ProcessResult DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half, - const std::string &base_prefix); + ProcessResult DoubledCellHypothesis(RugnuxObserver *observer, ProcessResult half); std::optional RefinedPasses(RugnuxObserver *observer, ProcessResult pass1, const RefinedPassStart &start); // The rotation-scale walk, its probes on copies of the run, while the refined passes run on a copy diff --git a/rugnux/RugnuxCalibration.cpp b/rugnux/RugnuxCalibration.cpp index 334d5d6e8..11d3e5cb6 100644 --- a/rugnux/RugnuxCalibration.cpp +++ b/rugnux/RugnuxCalibration.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -39,7 +37,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 4dd6d6dc7..9bdb53dd1 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -53,6 +53,8 @@ namespace { std::string RugnuxCommandLine(const ProcessConfig &config, const DiffractionExperiment &experiment, const std::string &input_file, + const std::string &output_prefix, + bool write_process_h5, const std::string &calibrant_name) { std::vector args; const bool azint = (config.mode == ProcessMode::AzimuthalIntegration); @@ -71,8 +73,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, args.push_back(val); }; - if (!config.output_prefix.empty()) - add("-o", config.output_prefix); + if (!output_prefix.empty()) + add("-o", output_prefix); // Every CPU the run may use is what rugnux takes without -N, on whatever machine runs the command. if (!config.nthreads_auto) add("-N", std::to_string(config.nthreads)); @@ -224,7 +226,7 @@ std::string RugnuxCommandLine(const ProcessConfig &config, // When merging, the CLI skips the large _process.h5 unless asked; emit the flag when it is // wanted so a copied command matches the GUI's "Save _process.h5" choice. (write_merged has // no CLI equivalent - the CLI always writes the .mtz/.cif when merging.) - if (config.write_process_h5) + if (write_process_h5) args.emplace_back("--write-process-h5"); // The low-resolution limit is on by default at the same value in both, so emit it only when // it differs - including 0, which is how the CLI spells "no limit". diff --git a/rugnux/RugnuxCommandLine.h b/rugnux/RugnuxCommandLine.h index f387a1686..b9f202709 100644 --- a/rugnux/RugnuxCommandLine.h +++ b/rugnux/RugnuxCommandLine.h @@ -14,7 +14,11 @@ class DiffractionExperiment; // matter for the run, not every obscure flag; geometry is taken from the input file, so geometry // overrides are not emitted. calibrant_name is only used by ProcessMode::Calibration - the config // carries the calibrant's rings, and --calibrant takes the name they were resolved from. +// output_prefix and write_process_h5 are where the run puts its files and whether a merging run writes +// its _process.h5 too. std::string RugnuxCommandLine(const ProcessConfig &config, const DiffractionExperiment &experiment, const std::string &input_file, + const std::string &output_prefix, + bool write_process_h5, const std::string &calibrant_name = ""); diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 2f4840c13..74ba2d3f2 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -39,7 +37,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp index f0eb897c6..aaf31b21c 100644 --- a/rugnux/RugnuxFirstPassAccept.cpp +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -39,7 +37,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp index ccc616dfe..554ca3fb3 100644 --- a/rugnux/RugnuxFirstPassRescues.cpp +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -39,7 +37,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 0510a79e9..8d31c3c00 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -4,7 +4,6 @@ #include #include "Rugnux.h" #include "../image_analysis/structure_refinement/ModelValidation.h" -#include "../writer/WriteModel.h" #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -19,7 +18,6 @@ #include #include #include -#include #include #include #include @@ -39,7 +37,6 @@ #include "../common/CUDAWrapper.h" #include "../common/JFJochException.h" #include "../common/time_utc.h" -#include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" @@ -713,7 +710,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { // it is written in - the pair the volume comparisons in Run() need. result.consensus_centering = rot->search_result.centering; // Log only: what the rest of the sweep says about the decisions taken on its first part. - if (!p.geometry_prepass && p.write_files && experiment_.IsRotationIndexing()) + if (!p.geometry_prepass && p.write_output && experiment_.IsRotationIndexing()) FarEndCheck(p, *rot); } result.consensus_cell = indexer->GetConsensusUnitCell(); diff --git a/rugnux/RugnuxLive.cpp b/rugnux/RugnuxLive.cpp index b05cf60b8..df46c56e6 100644 --- a/rugnux/RugnuxLive.cpp +++ b/rugnux/RugnuxLive.cpp @@ -6,7 +6,7 @@ #include "RugnuxDefaults.h" RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthreads, - const std::string &output_prefix, std::shared_ptr model) { + std::shared_ptr model) { RugnuxLiveSetup s; // The experiment as rugnux_cli.cpp sets it up for a rotation sweep with no options given. @@ -36,8 +36,6 @@ RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthr ProcessConfig &c = s.config; c.mode = ProcessMode::FullAnalysis; c.nthreads = nthreads; - c.output_prefix = output_prefix; - c.write_process_h5 = false; c.export_unmerged = false; c.spot_finding = RugnuxDefaultSpotFinding(true); c.rotation_indexing = true; diff --git a/rugnux/RugnuxLive.h b/rugnux/RugnuxLive.h index c1e0820e8..c24df0e26 100644 --- a/rugnux/RugnuxLive.h +++ b/rugnux/RugnuxLive.h @@ -9,13 +9,13 @@ #include "../reader/JFJochReaderDataset.h" // A rotation collection processed inside the broker while its frames arrive: what a bare `rugnux` on -// the same file would set, with the run made fast_live. The output goes to output_prefix - a place -// private to the broker, from which it builds the END message - and the model, where one is given, -// is placed and validated as `rugnux --model` does. +// the same file would set, with the run made fast_live. What it makes is returned in memory +// (ProcessResult::outputs), for the broker to build the END message from, and the model, where one is +// given, is placed and validated as `rugnux --model` does. struct RugnuxLiveSetup { DiffractionExperiment experiment; ProcessConfig config; }; RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthreads, - const std::string &output_prefix, std::shared_ptr model); + std::shared_ptr model); diff --git a/rugnux/RugnuxOutputFiles.cpp b/rugnux/RugnuxOutputFiles.cpp new file mode 100644 index 000000000..1570f47e0 --- /dev/null +++ b/rugnux/RugnuxOutputFiles.cpp @@ -0,0 +1,176 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RugnuxOutputFiles.h" + +#include +#include + +#include + +#include "../writer/FileWriter.h" +#include "../writer/WriteModel.h" +#include "../writer/WriteRugnuxFiles.h" + +namespace { + +void WriteCcp4(const gemmi::Grid &map, const std::string &path) { + gemmi::Ccp4 ccp4; + ccp4.grid = map; + ccp4.update_ccp4_header(2); + ccp4.write_ccp4_map(path); +} + +// The files of one output pass at `prefix`, all but its _process.h5. +void WriteOutputs(const RugnuxOutputs &o, const std::string &prefix, const OutputProvenance &provenance, + size_t nthreads, Logger &logger) { + if (o.merged) + WriteMergedOutput(*o.merged, prefix, provenance, nthreads, /*hkl=*/true); + if (o.twin) + WriteMergedOutput(*o.twin, prefix, provenance, nthreads, /*hkl=*/false); + if (o.model_prior) + WriteMergedOutput(*o.model_prior, prefix, provenance, nthreads, /*hkl=*/false); + if (o.maps) + WriteModelMaps(*o.maps, prefix, logger); + if (o.model_prior_maps) + WriteModelMaps(*o.model_prior_maps, prefix + "_mr", logger); + if (o.placed_model) + WritePlacedModel(*o.placed_model->model, o.placed_model->cell, o.placed_model->space_group, prefix, logger); + if (o.model_prior_placed_model) + WritePlacedModel(*o.model_prior_placed_model->model, o.model_prior_placed_model->cell, + o.model_prior_placed_model->space_group, prefix + "_mr", logger); + if (o.unmerged) { + WriteMtz(*o.unmerged, prefix + "_unmerged.mtz"); + logger.Info("Unmerged observations written to {}_unmerged.mtz", prefix); + } + if (o.unmerged_partials) { + WriteMtz(*o.unmerged_partials, prefix + "_unmerged_partials.mtz"); + logger.Info("Unmerged partials written to {}_unmerged_partials.mtz", prefix); + } + if (!o.detector_jpeg.empty()) { + const std::string path = prefix + "_detector.jpg"; + std::ofstream f(path, std::ofstream::out | std::ofstream::trunc | std::ofstream::binary); + if (f) { + f.write(o.detector_jpeg.data(), static_cast(o.detector_jpeg.size())); + logger.Info("Detector diagnostic image written to {}", path); + } else + logger.Warning("Cannot open {} for writing", path); + } + if (o.plot) { + WritePerImagePlot(prefix, *o.plot); + logger.Info("Per-image table written to {}_plot.txt", prefix); + } +} + +} // namespace + +ProcessFileSink::ProcessFileSink(std::string prefix) : prefix_(std::move(prefix)) {} + +ProcessFileSink::~ProcessFileSink() = default; + +void ProcessFileSink::Start(const StartMessage &start, const std::string &label) { + writer_.reset(); + StartMessage s = start; + s.file_prefix = prefix_ + label; + writer_ = std::make_unique(s, /*check_overwrite_at_start=*/true, /*trusted_path=*/true); +} + +void ProcessFileSink::Write(const DataMessage &msg) { + writer_->Write(msg); +} + +void ProcessFileSink::End(const EndMessage &end) { + writer_->WriteHDF5(end); + writer_->Finalize(); + writer_.reset(); +} + +void WriteRugnuxOutputs(const ProcessResult &result, const std::string &prefix, const OutputProvenance &provenance, + size_t nthreads, Logger &logger) { + if (result.outputs) { + WriteOutputs(*result.outputs, prefix, provenance, nthreads, logger); + if (result.outputs->process_start) { + StartMessage start = *result.outputs->process_start; + start.file_prefix = prefix; + FileWriter writer(start, /*check_overwrite_at_start=*/true, /*trusted_path=*/true); + writer.WriteHDF5(result.end_message); + writer.Finalize(); + } + } + if (result.doubled_cell) { + const ProcessResult::DoubledCell &h = *result.doubled_cell; + if (h.other_outputs) + WriteOutputs(*h.other_outputs, prefix + h.other_suffix, provenance, nthreads, logger); + // The doubled cell's pass put its per-image analysis under a name of its own; where it is the answer, + // the two _process.h5 files trade names, as the other files of the two passes are named. + std::error_code ec; + if (h.adopted && std::filesystem::exists(prefix + "_doubled_process.h5", ec)) { + if (std::filesystem::exists(prefix + "_process.h5", ec)) + std::filesystem::rename(prefix + "_process.h5", prefix + h.other_suffix + "_process.h5", ec); + std::filesystem::rename(prefix + "_doubled_process.h5", prefix + "_process.h5", ec); + } + } +} + +void WriteMergedOutput(const MergedOutput &merged, const std::string &prefix, const OutputProvenance &provenance, + size_t nthreads, bool hkl) { + // An MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / phenix + // reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 for SHELXL + // refinement and SHELXC / ANODE substructure solution). The .hkl holds the unmerged scaled fulls where the + // merge provides them (rotation), the merged reflections otherwise (stills). + const std::string file = prefix + merged.suffix; + WriteMtz(merged.mtz, file + ".mtz"); + WriteMmcifReflections(merged.reflections, merged.cell, merged.experiment, merged.statistics, merged.error_model, + merged.twinning, provenance, file + ".cif", nthreads); + if (!hkl) + return; + if (!merged.scaled_fulls.empty()) + WriteShelxHklReflections(merged.scaled_fulls, file + ".hkl", nthreads); + else + WriteShelxHklReflections(merged.reflections, merged.experiment, file + ".hkl", nthreads); +} + +void WriteModelMaps(const ModelMaps &maps, const std::string &prefix, Logger &logger) { + WriteMtz(maps.coefficients, prefix + "_maps.mtz"); + if (maps.map_2fofc) + WriteCcp4(*maps.map_2fofc, prefix + "_2fofc.ccp4"); + if (maps.map_fofc) + WriteCcp4(*maps.map_fofc, prefix + "_fofc.ccp4"); + if (maps.map_anom) + WriteCcp4(*maps.map_anom, prefix + "_anom.ccp4"); + logger.Info("Model validation: wrote {}_maps.mtz{}", prefix, + maps.map_2fofc ? fmt::format(", {0}_2fofc.ccp4, {0}_fofc.ccp4{1}", prefix, + maps.map_anom ? ", " + prefix + "_anom.ccp4" : std::string()) + : std::string()); +} + +void WritePerImagePlot(const std::string &prefix, const PerImagePlot &plot) { + const std::string path = prefix + "_plot.txt"; + std::ofstream file(path, std::ofstream::out | std::ofstream::trunc); + if (!file) + throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot open " + path + " for writing"); + file << PerImagePlotText(plot); +} + +ReportDocument WriteResultReport(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + Logger &logger, + const RunProvenance &provenance) { + ReportDocument doc = BuildReportDocument(output_prefix, input_file, experiment, result, provenance); + if (output_prefix.empty()) + return doc; // "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 << RenderReportText(doc, provenance.developer); + } catch (const std::exception &e) { + logger.Warning("Could not write the results report {}: {}", filename, e.what()); + } + return doc; +} diff --git a/rugnux/RugnuxOutputFiles.h b/rugnux/RugnuxOutputFiles.h new file mode 100644 index 000000000..86f425576 --- /dev/null +++ b/rugnux/RugnuxOutputFiles.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 "Rugnux.h" +#include "ResultReport.h" + +class FileWriter; + +// The files of a run, written by the program that ran it - rugnux and the viewer - from what the analysis +// returned: the library itself writes none. + +//