From aca27be862ca50dcb08ceddf44005ae5f98f15d7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 17 Oct 2026 12:19:51 +0200 Subject: [PATCH] jfjoch_writer writes rugnux's .mtz, _maps.mtz and _model.cif from the END blocks New gemmi-only library JFJochRugnuxFiles (writer/), linked by JFJochWriter, so the writer needs no analysis library and builds with JFJOCH_WRITER_ONLY: - WriteMergedMtz / WriteMapCoefficientsMtz / PlacedModelStructure make the files from the structured blocks; WriteRugnuxResultFiles writes whichever blocks END carries beside the master, each under a temporary name and renamed into place, an existing file replaced only when the run's overwrite is set. A file that fails costs only itself (logged), not the run. - FillRugnuxResultBlocks reads the blocks back from the files rugnux wrote at a prefix (p.mtz, p_maps.mtz, p_model.cif), for the broker to fill END. WriteMtzReflections now goes through WriteMergedMtz (MergedReflectionsForWriter builds the block); its output is byte-identical to before, checked by md5 on a synthetic set in both the Friedel-merged and the anomalous layout. StreamWriter writes the files after the master and before the END ACK, so the ACK means they are in place; END receipt and each rename are logged at INFO with the wall-clock time, for last-frame-to-files latency. The broker's direct HDF5 pusher writes them too. Test: rugnux's three files -> blocks -> CBOR -> writer reproduces all of them (and _model.pdb) byte for byte; a ~25k-reflection END is 2.7 MB. Co-Authored-By: Claude Opus 5.5 --- image_analysis/CMakeLists.txt | 2 +- image_analysis/WriteReflections.cpp | 119 ++++++--------- image_analysis/WriteReflections.h | 9 ++ image_pusher/HDF5FilePusher.cpp | 1 + tests/CMakeLists.txt | 1 + tests/RugnuxResultFilesTest.cpp | 219 ++++++++++++++++++++++++++++ writer/CMakeLists.txt | 18 ++- writer/FileWriter.cpp | 6 + writer/FileWriter.h | 4 + writer/RugnuxResultBlocks.cpp | 140 ++++++++++++++++++ writer/RugnuxResultBlocks.h | 22 +++ writer/StreamWriter.cpp | 7 +- writer/WriteRugnuxFiles.cpp | 210 ++++++++++++++++++++++++++ writer/WriteRugnuxFiles.h | 33 +++++ 14 files changed, 706 insertions(+), 85 deletions(-) create mode 100644 tests/RugnuxResultFilesTest.cpp create mode 100644 writer/RugnuxResultBlocks.cpp create mode 100644 writer/RugnuxResultBlocks.h create mode 100644 writer/WriteRugnuxFiles.cpp create mode 100644 writer/WriteRugnuxFiles.h diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 098307e75..916ebb0e2 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -131,4 +131,4 @@ ADD_SUBDIRECTORY(azint) ADD_SUBDIRECTORY(roi) ADD_SUBDIRECTORY(structure_refinement) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochROIIntegration JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochStructureRefinement gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochROIIntegration JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochStructureRefinement JFJochRugnuxFiles gemmi) diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index faf3b5075..8cf24c2a5 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -7,6 +7,7 @@ #include "scale_merge/TwinningAnalysis.h" #include "bragg_integration/SystematicAbsence.h" #include "../common/ParallelFor.h" +#include "../writer/WriteRugnuxFiles.h" #include #include @@ -121,16 +122,16 @@ std::string GitDescription() { // (mtzformat.html), so a long value runs on over as many records as it needs, up to that limit, each // continuation starting with '+'. Readers strip the blanks at both ends of a record, so a record is // never cut just after a space: the space starts the next one, behind its '+'. -void AddProvenanceHistory(gemmi::Mtz &mtz, const OutputProvenance &provenance) { - const auto add = [&mtz](const std::string &label, const std::string &value) { +void AddProvenanceHistory(std::vector &history, const OutputProvenance &provenance) { + const auto add = [&history](const std::string &label, const std::string &value) { if (value.empty()) return; std::string rest = label + value; - while (!rest.empty() && mtz.history.size() < 30) { + while (!rest.empty() && history.size() < 30) { size_t n = std::min(80, rest.size()); while (n > 1 && n < rest.size() && rest[n - 1] == ' ') --n; - mtz.history.push_back(rest.substr(0, n)); + history.push_back(rest.substr(0, n)); rest = n < rest.size() ? "+" + rest.substr(n) : std::string(); } }; @@ -606,81 +607,49 @@ 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; + m.unit_cell = unitCell; + m.space_group_name = experiment.GetSpaceGroupOrP1().xhm(); + m.wavelength_A = experiment.GetWavelength_A(); + m.crystal_name = experiment.GetSampleName(); + m.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); + AddProvenanceHistory(m.history, provenance); + + 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)); + } + return m; +} + void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const OutputProvenance &provenance, const std::string &filename) { - gemmi::Mtz mtz; - - // Optional but recommended metadata - mtz.spacegroup = &experiment.GetSpaceGroupOrP1(); - mtz.set_cell_for_all(unitCell); - - // Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention). - mtz.title = "Rugnux merged reflections"; - mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); - AddProvenanceHistory(mtz, provenance); - - // Add dataset. 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. - mtz.add_base(); // the HKL_base dataset and the H K L columns - gemmi::Mtz::Dataset& ds = mtz.add_dataset("native"); - ds.crystal_name = experiment.GetSampleName(); - ds.wavelength = experiment.GetWavelength_A(); - - const int dataset_id = ds.id; - - // 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); - - mtz.add_column("IMEAN", 'J', dataset_id, -1, false); - mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false); - if (has_anom) { - mtz.add_column("I(+)", 'K', dataset_id, -1, false); - mtz.add_column("SIGI(+)", 'M', dataset_id, -1, false); - mtz.add_column("I(-)", 'K', dataset_id, -1, false); - mtz.add_column("SIGI(-)", 'M', dataset_id, -1, false); - } - mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude - mtz.add_column("SIGF", 'Q', dataset_id, -1, false); - if (has_anom) { - mtz.add_column("F(+)", 'G', dataset_id, -1, false); - mtz.add_column("SIGF(+)", 'L', dataset_id, -1, false); - mtz.add_column("F(-)", 'G', dataset_id, -1, false); - mtz.add_column("SIGF(-)", 'L', dataset_id, -1, false); - } - mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false); - - mtz.nreflections = static_cast(out_rows.size()); - mtz.data.reserve(out_rows.size() * 16); - for (const auto& row : out_rows) { - mtz.data.push_back(static_cast(row.h)); - mtz.data.push_back(static_cast(row.k)); - mtz.data.push_back(static_cast(row.l)); - mtz.data.push_back(row.Imean); - mtz.data.push_back(row.sImean); - if (has_anom) { - mtz.data.push_back(row.Ip); - mtz.data.push_back(row.sIp); - mtz.data.push_back(row.Im); - mtz.data.push_back(row.sIm); - } - mtz.data.push_back(row.Fmean); - mtz.data.push_back(row.sFmean); - if (has_anom) { - mtz.data.push_back(row.Fp); - mtz.data.push_back(row.sFp); - mtz.data.push_back(row.Fm); - mtz.data.push_back(row.sFm); - } - mtz.data.push_back(static_cast(row.rfree)); - } - mtz.write_to_file(filename); + WriteMergedMtz(MergedReflectionsForWriter(reflections, unitCell, experiment, provenance), filename); } namespace { @@ -948,7 +917,7 @@ gemmi::Mtz UnmergedMtz(const std::vector &outcomes, mtz.set_cell_for_all(unitCell); mtz.title = "Rugnux unmerged observations"; mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction"); - AddProvenanceHistory(mtz, provenance); + 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(); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 559b3e662..68d2ab778 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -9,6 +9,7 @@ #include "../common/Reflection.h" #include "../common/UnitCell.h" #include "../common/DiffractionExperiment.h" +#include "../common/JFJochMessages.h" #include "IntegrationOutcome.h" #include @@ -49,6 +50,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); + void WriteMtzReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, diff --git a/image_pusher/HDF5FilePusher.cpp b/image_pusher/HDF5FilePusher.cpp index bec85f54f..ac43e1661 100644 --- a/image_pusher/HDF5FilePusher.cpp +++ b/image_pusher/HDF5FilePusher.cpp @@ -66,6 +66,7 @@ bool HDF5FilePusher::EndDataCollection(const EndMessage &message) { if (!writer) throw JFJochException(JFJochExceptionCategory::WrongDAQState, "Image pusher not ready for writing images"); writer->WriteHDF5(message); + writer->WriteResultFiles(message, logger); writer->Finalize(); } catch (...) { writer.reset(); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c6dd46391..5f5c4b207 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -172,6 +172,7 @@ ADD_EXECUTABLE(jfjoch_test LoadReferenceMtzFreeFlagsTest.cpp LoadReferenceSfCifTest.cpp WriteReflectionsTest.cpp + RugnuxResultFilesTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/RugnuxResultFilesTest.cpp b/tests/RugnuxResultFilesTest.cpp new file mode 100644 index 000000000..d0910ab98 --- /dev/null +++ b/tests/RugnuxResultFilesTest.cpp @@ -0,0 +1,219 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include + +#include +#include +#include // count_atom_sites + +#include "../common/Definitions.h" +#include "../common/DiffractionExperiment.h" +#include "../common/Logger.h" +#include "../image_analysis/WriteReflections.h" +#include "../frame_serialize/CBORStream2Serializer.h" +#include "../frame_serialize/CBORStream2Deserializer.h" +#include "../writer/WriteModel.h" +#include "../writer/WriteRugnuxFiles.h" +#include "../writer/RugnuxResultBlocks.h" +#include "SyntheticMergedReflections.h" + +namespace { + constexpr UnitCell TEST_CELL{47.0f, 47.0f, 63.0f, 90.0f, 90.0f, 90.0f}; + + std::string ReadFile(const std::string &path) { + std::ifstream f(path, std::ios::binary); + return {std::istreambuf_iterator(f), std::istreambuf_iterator()}; + } + + std::vector Labels(const gemmi::Mtz &mtz) { + std::vector v; + for (const auto &c : mtz.columns) + v.push_back(c.label); + return v; + } + + // A small model with what a placed model carries: two chains, a ligand, a water, an alternate + // conformation, an anisotropic atom, an insertion code and a charge. + gemmi::Structure TestModel() { + gemmi::Structure st; + st.name = "test"; + st.cell = TEST_CELL; + st.spacegroup_hm = "P 43 21 2"; + gemmi::Model &m = st.models.emplace_back(1); + double x = 1.0; + auto add_atom = [&x](gemmi::Residue &r, const std::string &name, const char *el, char altloc = 0) { + gemmi::Atom &a = r.atoms.emplace_back(); + a.name = name; + a.element = gemmi::Element(el); + a.altloc = altloc; + a.pos = gemmi::Position(x, 2.0 * x + 0.123, 3.0 * x - 0.456); + a.b_iso = 15.0f + static_cast(x); + a.occ = altloc ? 0.5f : 1.0f; + x += 1.37; + return &a; + }; + for (const std::string chain : {"A", "B"}) { + gemmi::Chain &ch = m.chains.emplace_back(chain); + for (int i = 1; i <= 3; i++) { + gemmi::Residue &r = ch.residues.emplace_back(gemmi::ResidueId{gemmi::SeqId(i, i == 3 ? 'A' : ' '), {}, "ALA"}); + r.het_flag = 'A'; + add_atom(r, "N", "N"); + add_atom(r, "CA", "C"); + add_atom(r, "C", "C"); + add_atom(r, "O", "O"); + add_atom(r, "CB", "C", 'A'); + add_atom(r, "CB", "C", 'B'); + } + } + gemmi::Residue &zn = m.chains[0].residues.emplace_back(gemmi::ResidueId{gemmi::SeqId(101, ' '), {}, "ZN"}); + zn.het_flag = 'H'; + gemmi::Atom *a = add_atom(zn, "ZN", "Zn"); + a->charge = 2; + a->aniso = {0.2f, 0.25f, 0.3f, 0.01f, -0.02f, 0.03f}; + gemmi::Residue &w = m.chains[1].residues.emplace_back(gemmi::ResidueId{gemmi::SeqId(201, ' '), {}, "HOH"}); + w.het_flag = 'H'; + add_atom(w, "O", "O"); + st.setup_cell_images(); + return st; + } +} + +TEST_CASE("RugnuxResultFiles: rugnux's files -> END blocks -> CBOR -> writer, byte for byte", "[rugnux_result_files]") { + Logger logger("RugnuxResultFiles"); + const std::filesystem::path dir = std::filesystem::temp_directory_path() / "rugnux_result_files_test"; + std::filesystem::remove_all(dir); + std::filesystem::create_directories(dir); + const std::string original = (dir / "rugnux").string(); + const std::string written = (dir / "writer").string(); + + // About the size of a standard-crystal run (~25k reflections), with an anomalous split. + jfjoch_test::SyntheticMergeParams params; + params.true_space_group = "P 43 21 2"; + params.twin_supergroup = "P 43 21 2"; + params.d_min_A = 2.25; + auto reflections = jfjoch_test::GenerateSyntheticMerged(params); + for (size_t i = 0; i < reflections.size(); i++) { + auto &r = reflections[i]; + r.I_plus = r.I * 1.01f; r.sigma_plus = r.sigma * 1.4f; + r.F_plus = std::sqrt(std::max(r.I, 1.0f)); r.sigmaF_plus = 0.1f; + if (i % 7 != 0) { + r.I_minus = r.I * 0.99f; r.sigma_minus = r.sigma * 1.4f; + r.F_minus = std::sqrt(std::max(r.I, 1.0f)); r.sigmaF_minus = 0.12f; + } + r.rfree_flag = (i % 20 == 0); + } + DiffractionExperiment experiment; + experiment.IncidentEnergy_keV(12.7f); + experiment.SpaceGroupNumber(96); + experiment.SetUnitCell(TEST_CELL); + OutputProvenance provenance; + provenance.command_line = "rugnux -N 8 input_master.h5"; + provenance.input_file = "/some/where/input_master.h5"; + WriteMtzReflections(reflections, TEST_CELL, experiment, provenance, original + ".mtz"); + + MapCoefficientsMessage maps; + maps.unit_cell = TEST_CELL; + maps.space_group_name = "P 43 21 2"; + maps.wavelength_A = experiment.GetWavelength_A(); + for (size_t i = 0; i < reflections.size(); i++) { + const float phi = static_cast((i * 37) % 360); + maps.h.push_back(reflections[i].h); + maps.k.push_back(reflections[i].k); + maps.l.push_back(reflections[i].l); + maps.fp.push_back(reflections[i].F); + maps.fc.push_back(reflections[i].F * 0.9f); + maps.phic.push_back(phi); + maps.fwt.push_back(reflections[i].F * 1.1f); + maps.phwt.push_back(phi); + maps.delfwt.push_back(reflections[i].F * 0.1f); + maps.phdelwt.push_back(phi); + maps.fom.push_back(0.8f); + maps.free_flag.push_back(reflections[i].rfree_flag ? 0 : 1); + } + WriteMapCoefficientsMtz(maps, original + "_maps.mtz"); + + WritePlacedModel(TestModel(), TEST_CELL, experiment.GetSpaceGroupOrP1(), original, logger); + + EndMessage end{}; + end.max_image_number = 100; + FillRugnuxResultBlocks(original, end); + REQUIRE(end.merged); + REQUIRE(end.map_coefficients); + REQUIRE(end.placed_model); + CHECK(end.merged->h.size() == reflections.size()); + CHECK(end.merged->i_plus.size() == reflections.size()); + CHECK(end.placed_model->atom_name.size() == 2 * 3 * 6 + 2); + + std::vector buffer(MESSAGE_SIZE_FOR_START_END); + CBORStream2Serializer serializer(buffer.data(), buffer.size()); + REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(end)); + std::cout << "END message with " << reflections.size() << " merged reflections, as many map coefficients and " + << end.placed_model->atom_name.size() << " atoms: " << serializer.GetBufferSize() << " bytes" << std::endl; + auto deserialized = CBORStream2Deserialize(buffer.data(), serializer.GetBufferSize()); + REQUIRE(deserialized); + REQUIRE(deserialized->end_message); + const EndMessage &out = *deserialized->end_message; + REQUIRE(out.merged); + REQUIRE(out.map_coefficients); + REQUIRE(out.placed_model); + CHECK(out.merged->history == end.merged->history); + CHECK(out.placed_model->x == end.placed_model->x); + + WriteRugnuxResultFiles(out, written, false, logger); + + CHECK(ReadFile(written + ".mtz") == ReadFile(original + ".mtz")); + CHECK(ReadFile(written + "_maps.mtz") == ReadFile(original + "_maps.mtz")); + CHECK(ReadFile(written + "_model.cif") == ReadFile(original + "_model.cif")); + CHECK(ReadFile(written + "_model.pdb") == ReadFile(original + "_model.pdb")); + + const gemmi::Mtz mtz = gemmi::read_mtz_file(written + ".mtz"); + CHECK(Labels(mtz) == std::vector{"H", "K", "L", "IMEAN", "SIGIMEAN", "I(+)", "SIGI(+)", "I(-)", + "SIGI(-)", "F", "SIGF", "F(+)", "SIGF(+)", "F(-)", "SIGF(-)", + "FreeR_flag"}); + CHECK(mtz.nreflections == static_cast(reflections.size())); + CHECK(mtz.datasets.back().wavelength == Catch::Approx(experiment.GetWavelength_A()).epsilon(1e-5)); + const gemmi::Mtz maps_mtz = gemmi::read_mtz_file(written + "_maps.mtz"); + CHECK(Labels(maps_mtz) == std::vector{"H", "K", "L", "FP", "FC", "PHIC", "FWT", "PHWT", "DELFWT", + "PHDELWT", "FOM", "FREE"}); + CHECK(maps_mtz.get_column_with_label("K")[5] == static_cast(maps.k[5])); + CHECK(maps_mtz.get_column_with_label("PHWT")[5] == maps.phwt[5]); + const gemmi::Structure st = gemmi::read_structure_gz(written + "_model.cif"); + CHECK(gemmi::count_atom_sites(st) == 2 * 3 * 6 + 2); + + std::filesystem::remove_all(dir); +} + +TEST_CASE("RugnuxResultFiles: an existing file is kept when overwrite is off", "[rugnux_result_files]") { + Logger logger("RugnuxResultFiles"); + const std::filesystem::path dir = std::filesystem::temp_directory_path() / "rugnux_result_files_overwrite"; + std::filesystem::remove_all(dir); + std::filesystem::create_directories(dir); + const std::string prefix = (dir / "run").string(); + + EndMessage end{}; + end.merged = MergedReflectionsMessage{}; + end.merged->unit_cell = TEST_CELL; + end.merged->space_group_name = "P 43 21 2"; + end.merged->h = {1}; end.merged->k = {2}; end.merged->l = {3}; + end.merged->i_mean = {10}; end.merged->sigma_i_mean = {1}; + end.merged->f = {3}; end.merged->sigma_f = {0.2}; + end.merged->free_flag = {1}; + + { std::ofstream f(prefix + ".mtz"); f << "keep"; } + WriteRugnuxResultFiles(end, prefix, false, logger); + CHECK(ReadFile(prefix + ".mtz") == "keep"); + CHECK(!std::filesystem::exists(prefix + ".mtz.tmp")); + + WriteRugnuxResultFiles(end, prefix, true, logger); + const gemmi::Mtz mtz = gemmi::read_mtz_file(prefix + ".mtz"); + CHECK(mtz.nreflections == 1); + CHECK(Labels(mtz) == std::vector{"H", "K", "L", "IMEAN", "SIGIMEAN", "F", "SIGF", "FreeR_flag"}); + + std::filesystem::remove_all(dir); +} diff --git a/writer/CMakeLists.txt b/writer/CMakeLists.txt index 2f6f43738..c96b830b7 100644 --- a/writer/CMakeLists.txt +++ b/writer/CMakeLists.txt @@ -31,14 +31,18 @@ ADD_LIBRARY(JFJochWriter STATIC HDF5DataFilePluginPerformance.h ) -TARGET_LINK_LIBRARIES(JFJochWriter JFJochPreview JFJochLogger JFJochHDF5Wrappers CBORStream2FrameSerialize) +# The placed atomic model (mmCIF/PDB) that rugnux writes beside its merged output. +ADD_LIBRARY(JFJochModelWriter STATIC WriteModel.cpp WriteModel.h) +TARGET_LINK_LIBRARIES(JFJochModelWriter JFJochLogger gemmi) -# The placed atomic model (mmCIF/PDB) that rugnux writes beside its merged output. Its own target, -# because it needs gemmi, which JFJochWriter does not link and a writer-only build does not want. -IF (NOT JFJOCH_WRITER_ONLY) - ADD_LIBRARY(JFJochModelWriter STATIC WriteModel.cpp WriteModel.h) - TARGET_LINK_LIBRARIES(JFJochModelWriter JFJochLogger gemmi) -ENDIF() +# rugnux's merged MTZ, map-coefficient MTZ and placed model, to and from the structured END-message +# blocks, so that jfjoch_writer writes them with gemmi alone and no analysis library. +ADD_LIBRARY(JFJochRugnuxFiles STATIC WriteRugnuxFiles.cpp WriteRugnuxFiles.h + RugnuxResultBlocks.cpp RugnuxResultBlocks.h) +TARGET_LINK_LIBRARIES(JFJochRugnuxFiles JFJochModelWriter JFJochLogger gemmi) + +TARGET_LINK_LIBRARIES(JFJochWriter JFJochPreview JFJochLogger JFJochHDF5Wrappers CBORStream2FrameSerialize + JFJochRugnuxFiles) # The streaming writer and the jfjoch_writer service belong to the Linux server stack, not the # portable viewer subset: JFJochStreamWriter pulls in JFJochImagePuller (not built in a viewer-only diff --git a/writer/FileWriter.cpp b/writer/FileWriter.cpp index 5d6b42432..92d9a3804 100644 --- a/writer/FileWriter.cpp +++ b/writer/FileWriter.cpp @@ -9,6 +9,7 @@ #include "../common/Logger.h" #include "../common/JFJochException.h" #include "JFJochDecompress.h" +#include "WriteRugnuxFiles.h" FileWriter::FileWriter(const StartMessage &request, bool check_overwrite_at_start, bool trusted_path) : start_message(request) { @@ -339,3 +340,8 @@ void FileWriter::WriteHDF5(const EndMessage &msg) { std::rethrow_exception(first_exception); } } + +void FileWriter::WriteResultFiles(const EndMessage &msg, Logger &logger) { + if (master_file) + WriteRugnuxResultFiles(msg, start_message.file_prefix, start_message.overwrite.value_or(false), logger); +} diff --git a/writer/FileWriter.h b/writer/FileWriter.h index c4ba5494e..df96a3e1c 100644 --- a/writer/FileWriter.h +++ b/writer/FileWriter.h @@ -9,6 +9,7 @@ #include "../common/JFJochMessages.h" #include "../common/ZMQWrappers.h" #include "HDF5NXmx.h" +#include "../common/Logger.h" #include #include #include @@ -59,6 +60,9 @@ public: void WriteHDF5(const DataMessage& msg); void WriteHDF5(const CompressedImage& msg); void WriteHDF5(const EndMessage& msg); + // The merged result files the END message carries (WriteRugnuxResultFiles), beside the master + // file, by the writer that owns it. Call after WriteHDF5(EndMessage). + void WriteResultFiles(const EndMessage& msg, Logger &logger); std::vector Finalize(); void SetupFinalizedFileSocket(const std::string &addr); // Called for every data file once it is closed and renamed into place. diff --git a/writer/RugnuxResultBlocks.cpp b/writer/RugnuxResultBlocks.cpp new file mode 100644 index 000000000..c7ac1d621 --- /dev/null +++ b/writer/RugnuxResultBlocks.cpp @@ -0,0 +1,140 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RugnuxResultBlocks.h" + +#include + +#include +#include + +namespace { + +UnitCell ToUnitCell(const gemmi::UnitCell &c) { + return UnitCell{static_cast(c.a), static_cast(c.b), static_cast(c.c), + static_cast(c.alpha), static_cast(c.beta), static_cast(c.gamma)}; +} + +// The column with this label, empty where the file has none. +std::vector Column(const gemmi::Mtz &mtz, const std::string &label) { + std::vector v; + if (const gemmi::Mtz::Column *col = mtz.column_with_label(label)) + for (int i = 0; i < col->size(); i++) + v.push_back((*col)[i]); + return v; +} + +std::vector IntColumn(const gemmi::Mtz &mtz, const std::string &label) { + std::vector v; + for (float x : Column(mtz, label)) + v.push_back(static_cast(x)); + return v; +} + +std::vector FlagColumn(const gemmi::Mtz &mtz, const std::string &label) { + std::vector v; + for (float x : Column(mtz, label)) + v.push_back(static_cast(x)); + return v; +} + +std::string SpaceGroupName(const gemmi::Mtz &mtz) { + return mtz.spacegroup != nullptr ? mtz.spacegroup->xhm() : std::string(); +} + +} // namespace + +MergedReflectionsMessage MergedReflectionsFromMtz(const std::string &path) { + const gemmi::Mtz mtz = gemmi::read_mtz_file(path); + MergedReflectionsMessage m; + m.unit_cell = ToUnitCell(mtz.cell); + m.space_group_name = SpaceGroupName(mtz); + m.wavelength_A = static_cast(mtz.datasets.back().wavelength); + m.crystal_name = mtz.datasets.back().crystal_name; + m.history = mtz.history; + m.h = IntColumn(mtz, "H"); + m.k = IntColumn(mtz, "K"); + m.l = IntColumn(mtz, "L"); + m.i_mean = Column(mtz, "IMEAN"); + m.sigma_i_mean = Column(mtz, "SIGIMEAN"); + m.i_plus = Column(mtz, "I(+)"); + m.sigma_i_plus = Column(mtz, "SIGI(+)"); + m.i_minus = Column(mtz, "I(-)"); + m.sigma_i_minus = Column(mtz, "SIGI(-)"); + m.f = Column(mtz, "F"); + m.sigma_f = Column(mtz, "SIGF"); + m.f_plus = Column(mtz, "F(+)"); + m.sigma_f_plus = Column(mtz, "SIGF(+)"); + m.f_minus = Column(mtz, "F(-)"); + m.sigma_f_minus = Column(mtz, "SIGF(-)"); + m.free_flag = FlagColumn(mtz, "FreeR_flag"); + return m; +} + +MapCoefficientsMessage MapCoefficientsFromMtz(const std::string &path) { + const gemmi::Mtz mtz = gemmi::read_mtz_file(path); + MapCoefficientsMessage c; + c.unit_cell = ToUnitCell(mtz.cell); + c.space_group_name = SpaceGroupName(mtz); + c.wavelength_A = static_cast(mtz.datasets.back().wavelength); + c.h = IntColumn(mtz, "H"); + c.k = IntColumn(mtz, "K"); + c.l = IntColumn(mtz, "L"); + c.fp = Column(mtz, "FP"); + c.fc = Column(mtz, "FC"); + c.phic = Column(mtz, "PHIC"); + c.fwt = Column(mtz, "FWT"); + c.phwt = Column(mtz, "PHWT"); + c.delfwt = Column(mtz, "DELFWT"); + c.phdelwt = Column(mtz, "PHDELWT"); + c.fom = Column(mtz, "FOM"); + c.free_flag = FlagColumn(mtz, "FREE"); + return c; +} + +PlacedModelMessage PlacedModelFromStructure(const gemmi::Structure &st, const std::string &data_space_group_name) { + PlacedModelMessage p; + p.name = st.name; + p.unit_cell = ToUnitCell(st.cell); + p.space_group_name = data_space_group_name; + p.model_space_group_name = st.spacegroup_hm; + bool has_aniso = false; + for (const gemmi::Model &m : st.models) + for (const gemmi::Chain &ch : m.chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) { + p.model_number.push_back(m.num); + p.chain.push_back(ch.name); + p.residue_name.push_back(r.name); + p.seq_num.push_back(*r.seqid.num); + p.icode.push_back(static_cast(r.seqid.icode)); + p.het_flag.push_back(static_cast(r.het_flag)); + p.atom_name.push_back(a.name); + p.element.push_back(a.element.name()); + p.altloc.push_back(static_cast(a.altloc)); + p.charge.push_back(a.charge); + p.x.push_back(a.pos.x); + p.y.push_back(a.pos.y); + p.z.push_back(a.pos.z); + p.b_iso.push_back(a.b_iso); + p.occupancy.push_back(a.occ); + p.aniso.insert(p.aniso.end(), {a.aniso.u11, a.aniso.u22, a.aniso.u33, + a.aniso.u12, a.aniso.u13, a.aniso.u23}); + has_aniso = has_aniso || a.aniso.nonzero(); + } + if (!has_aniso) + p.aniso.clear(); + return p; +} + +void FillRugnuxResultBlocks(const std::string &prefix, EndMessage &end) { + if (std::filesystem::exists(prefix + ".mtz")) + end.merged = MergedReflectionsFromMtz(prefix + ".mtz"); + if (std::filesystem::exists(prefix + "_maps.mtz")) + end.map_coefficients = MapCoefficientsFromMtz(prefix + "_maps.mtz"); + if (std::filesystem::exists(prefix + "_model.cif")) { + const gemmi::Structure st = gemmi::read_structure_gz(prefix + "_model.cif"); + end.placed_model = PlacedModelFromStructure(st, end.merged ? end.merged->space_group_name + : st.spacegroup_hm); + } +} diff --git a/writer/RugnuxResultBlocks.h b/writer/RugnuxResultBlocks.h new file mode 100644 index 000000000..f0198e46f --- /dev/null +++ b/writer/RugnuxResultBlocks.h @@ -0,0 +1,22 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "gemmi/model.hpp" + +#include "../common/JFJochMessages.h" + +// The END-message blocks (common/JFJochMessages.h) read back from the files rugnux wrote at +// `prefix`: .mtz into `merged`, _maps.mtz into `map_coefficients` and +// _model.cif into `placed_model`, each only where its file exists. The writer makes the same +// files again from them (WriteRugnuxResultFiles, writer/WriteRugnuxFiles.h). +void FillRugnuxResultBlocks(const std::string &prefix, EndMessage &end); + +MergedReflectionsMessage MergedReflectionsFromMtz(const std::string &path); +MapCoefficientsMessage MapCoefficientsFromMtz(const std::string &path); +// data_space_group_name: the group of the reflection files beside the model, which is what +// WritePlacedModel is to be given; the model's own group is taken from `st`. +PlacedModelMessage PlacedModelFromStructure(const gemmi::Structure &st, const std::string &data_space_group_name); diff --git a/writer/StreamWriter.cpp b/writer/StreamWriter.cpp index efac1d632..616af4955 100644 --- a/writer/StreamWriter.cpp +++ b/writer/StreamWriter.cpp @@ -4,6 +4,7 @@ #include "../common/JFJochException.h" #include "StreamWriter.h" +#include #include #include "FileWriter.h" @@ -203,14 +204,16 @@ void StreamWriter::ProcessEndMessage() { if (state == StreamWriterState::Idle || state == StreamWriterState::Finalized) return; - if (verbose) - logger.Info("Received end message"); + logger.Info("Received end message at wall clock {:.3f} s", + std::chrono::duration(std::chrono::system_clock::now().time_since_epoch()).count()); if (state != StreamWriterState::Error) { try { if ((image_puller_output.cbor->end_message->max_image_number == 0) && (max_image_number > 0)) image_puller_output.cbor->end_message->max_image_number = max_image_number; file_writer->WriteHDF5(*image_puller_output.cbor->end_message); + // Before the END ACK below, so that the ACK says these files are in place too. + file_writer->WriteResultFiles(*image_puller_output.cbor->end_message, logger); } catch (const std::exception &e) { // std::exception for the same reason as in ProcessDataImage. Here the cost of a leak is the // highest: the END ACK and the notification below are never sent, the state never reaches diff --git a/writer/WriteRugnuxFiles.cpp b/writer/WriteRugnuxFiles.cpp new file mode 100644 index 000000000..2f9449485 --- /dev/null +++ b/writer/WriteRugnuxFiles.cpp @@ -0,0 +1,210 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "WriteRugnuxFiles.h" + +#include +#include + +#include + +#include "WriteModel.h" +#include "../common/Logger.h" + +void WriteMergedMtz(const MergedReflectionsMessage &merged, const std::string &filename) { + gemmi::Mtz mtz; + + // Optional but recommended metadata + mtz.spacegroup = gemmi::find_spacegroup_by_name(merged.space_group_name); + mtz.set_cell_for_all(merged.unit_cell); + + // Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention). + mtz.title = "Rugnux merged reflections"; + mtz.history = merged.history; + + // Add dataset. 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. + mtz.add_base(); // the HKL_base dataset and the H K L columns + gemmi::Mtz::Dataset& ds = mtz.add_dataset("native"); + ds.crystal_name = merged.crystal_name; + ds.wavelength = merged.wavelength_A; + + const int dataset_id = ds.id; + + // 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. + const bool has_anom = !merged.i_plus.empty(); + + mtz.add_column("IMEAN", 'J', dataset_id, -1, false); + mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false); + if (has_anom) { + mtz.add_column("I(+)", 'K', dataset_id, -1, false); + mtz.add_column("SIGI(+)", 'M', dataset_id, -1, false); + mtz.add_column("I(-)", 'K', dataset_id, -1, false); + mtz.add_column("SIGI(-)", 'M', dataset_id, -1, false); + } + mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude + mtz.add_column("SIGF", 'Q', dataset_id, -1, false); + if (has_anom) { + mtz.add_column("F(+)", 'G', dataset_id, -1, false); + mtz.add_column("SIGF(+)", 'L', dataset_id, -1, false); + mtz.add_column("F(-)", 'G', dataset_id, -1, false); + mtz.add_column("SIGF(-)", 'L', dataset_id, -1, false); + } + mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false); + + const size_t n = merged.h.size(); + mtz.nreflections = static_cast(n); + mtz.data.reserve(n * 16); + for (size_t i = 0; i < n; i++) { + mtz.data.push_back(static_cast(merged.h[i])); + mtz.data.push_back(static_cast(merged.k[i])); + mtz.data.push_back(static_cast(merged.l[i])); + mtz.data.push_back(merged.i_mean[i]); + mtz.data.push_back(merged.sigma_i_mean[i]); + if (has_anom) { + mtz.data.push_back(merged.i_plus[i]); + mtz.data.push_back(merged.sigma_i_plus[i]); + mtz.data.push_back(merged.i_minus[i]); + mtz.data.push_back(merged.sigma_i_minus[i]); + } + mtz.data.push_back(merged.f[i]); + mtz.data.push_back(merged.sigma_f[i]); + if (has_anom) { + mtz.data.push_back(merged.f_plus[i]); + mtz.data.push_back(merged.sigma_f_plus[i]); + mtz.data.push_back(merged.f_minus[i]); + mtz.data.push_back(merged.sigma_f_minus[i]); + } + mtz.data.push_back(static_cast(merged.free_flag[i])); + } + mtz.write_to_file(filename); +} + +void WriteMapCoefficientsMtz(const MapCoefficientsMessage &coefficients, const std::string &filename) { + gemmi::Mtz mtz(true); + mtz.spacegroup = gemmi::find_spacegroup_by_name(coefficients.space_group_name); + mtz.set_cell_for_all(coefficients.unit_cell); + mtz.add_dataset("model_validation"); + mtz.datasets.back().wavelength = coefficients.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); + const size_t n = coefficients.h.size(); + mtz.data.reserve(n * 12); + for (size_t i = 0; i < n; i++) + mtz.data.insert(mtz.data.end(), {static_cast(coefficients.h[i]), + static_cast(coefficients.k[i]), + static_cast(coefficients.l[i]), + coefficients.fp[i], coefficients.fc[i], coefficients.phic[i], + coefficients.fwt[i], coefficients.phwt[i], + coefficients.delfwt[i], coefficients.phdelwt[i], + coefficients.fom[i], + static_cast(coefficients.free_flag[i])}); + mtz.nreflections = static_cast(n); + mtz.write_to_file(filename); +} + +gemmi::Structure PlacedModelStructure(const PlacedModelMessage &model) { + gemmi::Structure st; + st.name = model.name; + st.cell = model.unit_cell; + st.spacegroup_hm = model.model_space_group_name; + const bool has_aniso = !model.aniso.empty(); + gemmi::Model *m = nullptr; + gemmi::Chain *ch = nullptr; + gemmi::Residue *res = nullptr; + for (size_t i = 0; i < model.atom_name.size(); i++) { + if (m == nullptr || m->num != model.model_number[i]) { + m = &st.models.emplace_back(model.model_number[i]); + ch = nullptr; + } + if (ch == nullptr || ch->name != model.chain[i]) { + ch = &m->chains.emplace_back(model.chain[i]); + res = nullptr; + } + const gemmi::SeqId seqid(model.seq_num[i], static_cast(model.icode[i])); + if (res == nullptr || res->name != model.residue_name[i] || res->seqid != seqid + || res->het_flag != static_cast(model.het_flag[i])) { + res = &ch->residues.emplace_back(gemmi::ResidueId{seqid, {}, model.residue_name[i]}); + res->het_flag = static_cast(model.het_flag[i]); + } + gemmi::Atom &a = res->atoms.emplace_back(); + a.name = model.atom_name[i]; + a.element = gemmi::Element(model.element[i]); + a.altloc = static_cast(model.altloc[i]); + a.charge = static_cast(model.charge[i]); + a.pos = gemmi::Position(model.x[i], model.y[i], model.z[i]); + a.b_iso = model.b_iso[i]; + a.occ = model.occupancy[i]; + if (has_aniso) + a.aniso = {model.aniso[6 * i], model.aniso[6 * i + 1], model.aniso[6 * i + 2], + model.aniso[6 * i + 3], model.aniso[6 * i + 4], model.aniso[6 * i + 5]}; + } + st.setup_cell_images(); + return st; +} + +namespace { + +double WallClockSeconds() { + return std::chrono::duration(std::chrono::system_clock::now().time_since_epoch()).count(); +} + +// Rename tmp onto final, unless final is there already and may not be overwritten. +void RenameIntoPlace(const std::string &tmp, const std::string &final_name, bool overwrite, Logger &logger) { + if (!overwrite && std::filesystem::exists(final_name)) { + std::filesystem::remove(tmp); + logger.Error("{} already exists and overwrite is off - not written", final_name); + return; + } + std::filesystem::rename(tmp, final_name); + logger.Info("{} written, renamed into place at wall clock {:.3f} s", final_name, WallClockSeconds()); +} + +} // namespace + +void WriteRugnuxResultFiles(const EndMessage &end, const std::string &file_prefix, bool overwrite, + Logger &logger) { + // A file that cannot be written costs only itself: the master is already in place. + try { + if (end.merged) { + WriteMergedMtz(*end.merged, file_prefix + ".mtz.tmp"); + RenameIntoPlace(file_prefix + ".mtz.tmp", file_prefix + ".mtz", overwrite, logger); + } + } catch (const std::exception &e) { + logger.Error("Cannot write {}.mtz: {}", file_prefix, e.what()); + } + try { + if (end.map_coefficients) { + WriteMapCoefficientsMtz(*end.map_coefficients, file_prefix + "_maps.mtz.tmp"); + RenameIntoPlace(file_prefix + "_maps.mtz.tmp", file_prefix + "_maps.mtz", overwrite, logger); + } + } catch (const std::exception &e) { + logger.Error("Cannot write {}_maps.mtz: {}", file_prefix, e.what()); + } + try { + if (end.placed_model) { + const PlacedModelMessage &model = *end.placed_model; + // WritePlacedModel appends _model.cif and _model.pdb to the prefix it is given. + const std::string tmp_prefix = file_prefix + ".tmp"; + WritePlacedModel(PlacedModelStructure(model), model.unit_cell, + *gemmi::find_spacegroup_by_name(model.space_group_name), tmp_prefix, logger); + for (const std::string suffix : {"_model.cif", "_model.pdb"}) + if (std::filesystem::exists(tmp_prefix + suffix)) + RenameIntoPlace(tmp_prefix + suffix, file_prefix + suffix, overwrite, logger); + } + } catch (const std::exception &e) { + logger.Error("Cannot write {}_model.cif: {}", file_prefix, e.what()); + } +} diff --git a/writer/WriteRugnuxFiles.h b/writer/WriteRugnuxFiles.h new file mode 100644 index 000000000..f2fc6fe5f --- /dev/null +++ b/writer/WriteRugnuxFiles.h @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "gemmi/model.hpp" + +#include "../common/JFJochMessages.h" + +class Logger; + +// The files rugnux writes beside its merged output, made from the structured END-message blocks +// (common/JFJochMessages.h), so that jfjoch_writer can write them with no analysis library. rugnux +// writes its own .mtz and _maps.mtz through the same two functions. + +// .mtz: H K L IMEAN SIGIMEAN [I(+) SIGI(+) I(-) SIGI(-)] F SIGF [F(+) SIGF(+) F(-) SIGF(-)] +// FreeR_flag, the bracketed columns only where the block has them. +void WriteMergedMtz(const MergedReflectionsMessage &merged, const std::string &filename); + +// _maps.mtz: H K L FP FC PHIC FWT PHWT DELFWT PHDELWT FOM FREE. +void WriteMapCoefficientsMtz(const MapCoefficientsMessage &coefficients, const std::string &filename); + +// The model the block describes, for WritePlacedModel (writer/WriteModel.h). +gemmi::Structure PlacedModelStructure(const PlacedModelMessage &model); + +// The three files above, from whichever blocks `end` carries, at .mtz, _maps.mtz and +// _model.cif (WritePlacedModel adds _model.pdb where a PDB can hold the cell). Each is written under a +// temporary name and renamed into place. An existing file is replaced only with `overwrite`; otherwise +// it is kept and the new one is not written. +void WriteRugnuxResultFiles(const EndMessage &end, const std::string &file_prefix, bool overwrite, + Logger &logger);