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. + +//