diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index dd68afaee..bcb4dcd1e 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -368,6 +368,59 @@ 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); +// model_space_group_name is the group the coordinates were built in. +struct PlacedModelMessage { + std::string name; + UnitCell unit_cell; + std::string space_group_name; + std::string model_space_group_name; + std::vector model_number; + std::vector chain; + std::vector residue_name; + std::vector seq_num; + std::vector icode; // ' ' where there is none + std::vector het_flag; // 'A', 'H', or 0 where unknown + std::vector atom_name; + std::vector element; + std::vector altloc; // 0 where there is none + std::vector charge; + std::vector x, y, z; + std::vector b_iso, occupancy; + std::vector aniso; // 6 per atom (u11 u22 u33 u12 u13 u23), empty when none has one +}; + struct EndMessage { // The sample transformation chain, base first (see DetectorTransformation). Optional: when it is // empty the writer builds the same chain itself from the start message, so a producer that does @@ -481,6 +534,12 @@ struct EndMessage { // this says what was DONE about it. Offline (rugnux) only; the broker does not merge. std::vector frame_disposition; std::vector frame_disposition_codes; + + // 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 placed_model; }; struct MetadataMessage { diff --git a/docs/CBOR.md b/docs/CBOR.md index 9bcb8c441..b3e2e4b33 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -342,6 +342,9 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | pixel_sum | Array(int64) | Per-image sum of all valid pixels, excluding error/saturated pixels | | | image_scale_mosaicity | Array(float) | Scaling result: Image scale mosaicity \[deg\] | | | image_scale_cc | Array(float) | Scaling result: Image scale CC | | +| merged | object (optional) | Merged reflections of an in-broker analysis run, written by the writer as `.mtz`. `unit_cell`, `space_group_name` (gemmi extended H-M), `wavelength`, `crystal_name`, `history` (Array(string), MTZ HISTORY records), then one typed array per MTZ column: `h`, `k`, `l` (int32), `i_mean`, `sigma_i_mean`, `i_plus`, `sigma_i_plus`, `i_minus`, `sigma_i_minus`, `f`, `sigma_f`, `f_plus`, `sigma_f_plus`, `f_minus`, `sigma_f_minus` (float; the eight anomalous arrays empty where there is no Bijvoet split), `free_flag` (uint8, 0 = test set) | | +| map_coefficients | object (optional) | Map coefficients, written as `_maps.mtz`. `unit_cell`, `space_group_name`, `wavelength`, then `h`, `k`, `l` (int32), `fp`, `fc`, `phic`, `fwt`, `phwt`, `delfwt`, `phdelwt` (phases in degrees), `fom` (float), `free_flag` (uint8, 0 = test set) | | +| placed_model | object (optional) | Atomic model placed against the data, written as `_model.cif` (and `_model.pdb`). `name`, `unit_cell`, `space_group_name` (the data's), `model_space_group_name` (the coordinates'), then one array per atom: `model_number`, `seq_num`, `charge` (int32), `chain`, `residue_name`, `atom_name`, `element` (Array(string)), `icode`, `het_flag`, `altloc` (uint8), `x`, `y`, `z` (double, Angstrom), `b_iso`, `occupancy` (float), `aniso` (float, 6 per atom: U11 U22 U33 U12 U13 U23; empty where no atom has one) | | End-message vector fields are optional. When present, they provide master-file summary data so readers can inspect scan-level and per-image analysis results without opening every linked data file. Missing optional per-image values are encoded by the producer as zero unless otherwise noted. diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index a1f5ea29d..15d07c0e3 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -320,6 +320,19 @@ namespace { memcpy(v.data(), ptr, len); } + void GetCBORDoubleArray(CborValue &value, std::vector &v) { + if (GetCBORTag(value) != TagDoubleLE) + throw JFJochException(JFJochExceptionCategory::CBORError, "Incorrect array type tag"); + + auto [ptr, len] = GetCBORByteString(value); + + if (len % sizeof(double)) + throw JFJochException(JFJochExceptionCategory::CBORError, "Size mismatch"); + + v.resize(len / sizeof(double)); + memcpy(v.data(), ptr, len); + } + void GetCBORStringArray(CborValue &value, std::vector &v) { if (!cbor_value_is_array(&value)) throw JFJochException(JFJochExceptionCategory::CBORError, "Array expected"); @@ -1490,6 +1503,105 @@ namespace { cborErr(cbor_value_leave_container(&value, &map_value)); } + // A key this build does not know: a tagged item takes two advances (see ProcessEndMessageElement). + void SkipCBORValue(CborValue &value) { + if (cbor_value_is_tag(&value)) + cbor_value_advance(&value); + cbor_value_advance(&value); + } + + MergedReflectionsMessage ProcessMergedElement(CborValue &value) { + MergedReflectionsMessage m; + 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); + else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return m; + } + + MapCoefficientsMessage ProcessMapCoefficientsElement(CborValue &value) { + MapCoefficientsMessage m; + 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); + else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return m; + } + + PlacedModelMessage ProcessPlacedModelElement(CborValue &value) { + PlacedModelMessage m; + 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 == "name") m.name = GetCBORString(map_value); + else 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 == "model_space_group_name") m.model_space_group_name = GetCBORString(map_value); + else if (key == "model_number") GetCBORInt32Array(map_value, m.model_number); + else if (key == "chain") GetCBORStringArray(map_value, m.chain); + else if (key == "residue_name") GetCBORStringArray(map_value, m.residue_name); + else if (key == "seq_num") GetCBORInt32Array(map_value, m.seq_num); + else if (key == "icode") GetCBORUInt8Array(map_value, m.icode); + else if (key == "het_flag") GetCBORUInt8Array(map_value, m.het_flag); + else if (key == "atom_name") GetCBORStringArray(map_value, m.atom_name); + else if (key == "element") GetCBORStringArray(map_value, m.element); + else if (key == "altloc") GetCBORUInt8Array(map_value, m.altloc); + else if (key == "charge") GetCBORInt32Array(map_value, m.charge); + else if (key == "x") GetCBORDoubleArray(map_value, m.x); + else if (key == "y") GetCBORDoubleArray(map_value, m.y); + else if (key == "z") GetCBORDoubleArray(map_value, m.z); + else if (key == "b_iso") GetCBORFloatArray(map_value, m.b_iso); + else if (key == "occupancy") GetCBORFloatArray(map_value, m.occupancy); + else if (key == "aniso") GetCBORFloatArray(map_value, m.aniso); + else SkipCBORValue(map_value); + } + cborErr(cbor_value_leave_container(&value, &map_value)); + return m; + } + bool ProcessEndMessageElement(EndMessage &message, CborValue &value) { if (cbor_value_at_end(&value)) return false; @@ -1586,6 +1698,12 @@ namespace { message.space_group_number = GetCBORUInt(value); else if (key == "pixel_sum") GetCBORInt64Array(value, message.pixel_sum); + else if (key == "merged") + message.merged = ProcessMergedElement(value); + else if (key == "map_coefficients") + message.map_coefficients = ProcessMapCoefficientsElement(value); + else if (key == "placed_model") + message.placed_model = ProcessPlacedModelElement(value); else if (key == "rotation_lattice") { std::vector tmp; GetCBORFloatArray(value, tmp); diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 385723e54..82ee03e51 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -515,6 +515,89 @@ inline void CBOR_ENC(CborEncoder &encoder, const char* key, const UnitCell &val) cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); } +inline void CBOR_ENC(CborEncoder &encoder, const char* key, const std::vector& v) { + cborErr(cbor_encode_text_stringz(&encoder, key)); + cborErr(cbor_encode_tag(&encoder, TagDoubleLE)); + cborErr(cbor_encode_byte_string(&encoder, reinterpret_cast(v.data()), v.size() * sizeof(double))); +} + +inline void CBOR_ENC(CborEncoder &encoder, const char* key, const MergedReflectionsMessage &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, "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)); +} + +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_encoder_close_container(&encoder, &mapEncoder)); +} + +inline void CBOR_ENC(CborEncoder &encoder, const char* key, const PlacedModelMessage &m) { + CborEncoder mapEncoder; + cborErr(cbor_encode_text_stringz(&encoder, key)); + cborErr(cbor_encoder_create_map(&encoder, &mapEncoder, CborIndefiniteLength)); + CBOR_ENC(mapEncoder, "name", m.name); + CBOR_ENC(mapEncoder, "unit_cell", m.unit_cell); + CBOR_ENC(mapEncoder, "space_group_name", m.space_group_name); + CBOR_ENC(mapEncoder, "model_space_group_name", m.model_space_group_name); + CBOR_ENC(mapEncoder, "model_number", m.model_number); + CBOR_ENC(mapEncoder, "chain", m.chain); + CBOR_ENC(mapEncoder, "residue_name", m.residue_name); + CBOR_ENC(mapEncoder, "seq_num", m.seq_num); + CBOR_ENC(mapEncoder, "icode", m.icode); + CBOR_ENC(mapEncoder, "het_flag", m.het_flag); + CBOR_ENC(mapEncoder, "atom_name", m.atom_name); + CBOR_ENC(mapEncoder, "element", m.element); + CBOR_ENC(mapEncoder, "altloc", m.altloc); + CBOR_ENC(mapEncoder, "charge", m.charge); + CBOR_ENC(mapEncoder, "x", m.x); + CBOR_ENC(mapEncoder, "y", m.y); + CBOR_ENC(mapEncoder, "z", m.z); + CBOR_ENC(mapEncoder, "b_iso", m.b_iso); + CBOR_ENC(mapEncoder, "occupancy", m.occupancy); + CBOR_ENC(mapEncoder, "aniso", m.aniso); + cborErr(cbor_encoder_close_container(&encoder, &mapEncoder)); +} + // Add encoder for LatticeMessage inline void CBOR_ENC(CborEncoder &encoder, const char* key, const LatticeMessage &val) { CborEncoder mapEncoder; @@ -857,6 +940,9 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { // should take; the number is beside it for a reader written before the name existed. CBOR_ENC(mapEncoder, "space_group_name", message.space_group_name); CBOR_ENC(mapEncoder, "space_group_number", message.space_group_number); + CBOR_ENC(mapEncoder, "merged", message.merged); + CBOR_ENC(mapEncoder, "map_coefficients", message.map_coefficients); + CBOR_ENC(mapEncoder, "placed_model", message.placed_model); cborErr(cbor_encoder_close_container(&encoder, &mapEncoder));