diff --git a/common/Reflection.h b/common/Reflection.h index a6b5ad9c5..cc19f303f 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -68,13 +68,16 @@ struct Reflection { // measurement, every correction and the per-frame scale applied, its sigma under the error model the // merge weighted it with - but NOT averaged with its symmetry mates. h k l are the index it was // measured at, not its ASU representative. Written as the SHELX HKLF 4 file, so that the program -// reading it sees the equivalents and computes Rint itself. +// reading it sees the equivalents and computes Rint itself, and as the scaled unmerged block of the +// mmCIF. struct ScaledFull { int32_t h = 0; int32_t k = 0; int32_t l = 0; float I = NAN; float sigma = NAN; + float scale = NAN; // the factor the scaling applied: I is the integrated intensity times this + float image_number = NAN; // the frame the full peaks on (0-based, fractional) }; struct MergedReflection { diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index ca82c3cbc..c2de28fed 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -253,7 +253,10 @@ reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol. 4, 269-275 [doi:10.1038/nsb0497-269](https://doi.org/10.1038/nsb0497-269); P. A. Karplus and K. Diederichs, "Linking crystallographic model and data quality" (2012), Science 336, 1030-1033 -[doi:10.1126/science.1218231](https://doi.org/10.1126/science.1218231); K. Diederichs and +[doi:10.1126/science.1218231](https://doi.org/10.1126/science.1218231); M. S. Weiss, "Global +indicators of X-ray data quality" (2001), J. Appl. Cryst. 34, 130-135 +[doi:10.1107/S0021889800018227](https://doi.org/10.1107/S0021889800018227) (R_pim beside R_meas and +R_merge); K. Diederichs and P. A. Karplus, "Better models by discarding data?" (2013), Acta Cryst. D69, 1215-1222 [doi:10.1107/S0907444913001121](https://doi.org/10.1107/S0907444913001121). diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index bfb7de4f6..b5b39ccce 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -76,6 +76,7 @@ The methods draw on, and in places reimplement, solutions from: - P. R. Evans & G. N. Murshudov, "How good are my data and what is the resolution?", *Acta Cryst.* **D69** (2013), 1204-1214 (AIMLESS: the anisotropic deltaB as the range of the principal components, and diffraction limits from a cone about each principal direction). - G. Assmann, W. Brehm & K. Diederichs, "Identification of rogue datasets in serial crystallography", *J. Appl. Cryst.* **49** (2016), 1021-1028, and G. M. Assmann, M. Wang & K. Diederichs, *Acta Cryst.* **D76** (2020), 636-652 (XDSCC12: sigma-tau CC1/2, delta-CC1/2, the Fisher transformation and the rejection discipline the frame disposition follows). - K. Diederichs & P. A. Karplus, *Nat. Struct. Biol.* **4** (1997), 269-275, and P. A. Karplus & K. Diederichs, *Science* **336** (2012), 1030-1033 (R_meas / R_pim, CC1/2 and CC\*). +- M. S. Weiss, "Global indicators of X-ray data quality", *J. Appl. Cryst.* **34** (2001), 130-135 (R_pim). - IUCr Commission on Crystallographic Nomenclature, "Statistical descriptors in crystallography", *Acta Cryst.* **A45** (1989), 63-75, and *Acta Cryst.* **A51** (1995), 565-569 (uncertainty conventions). (list is not exhaustive; the full citations, with DOIs, are in [ACKNOWLEDGEMENT.md](ACKNOWLEDGEMENT.md)) diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 2a42120c3..6abcea11e 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -243,6 +243,32 @@ describe the same frame. An older file that was affected by this cannot be repai not recoverable after the fact); such a file now stops with a message naming both cells and the exact `-S`/`-C` override to merge it in its own setting, instead of failing inside the merge. +## The merged mmCIF (`.cif`) + +The `.cif` is written for deposition and follows the PDBx/mmCIF dictionary. It holds two data +blocks: + +- **`data_sample`** — the merged reflections (`_refln`: I, F and the Bijvoet mates, with the free set + in `_refln.status`) and what describes them: the cell and space group; `_software` with Rugnux + named twice, once for `data reduction` and once for `data scaling`; the experiment as the input + file states it (`_diffrn.ambient_temp` from the sample temperature, the wavelength, the source's + facility and beamline and its general class, the detector's make and model in + `_diffrn_detector.type` and its class in `_diffrn_detector.detector`, and the collection date from + the NXmx start time); and the merging statistics overall (`_reflns`) and per shell + (`_reflns_shell`): completeness, multiplicity, , R_merge, R_meas, R_pim, CC1/2 and CC\*, + and SigAno where an anomalous split was made. A quantity the input does not state is left out + rather than guessed. +- **`data_sample_unmerged`** (rotation) — the scaled **unmerged** observations the merge averaged, + as `_diffrn_refln`: the index each was measured at, the frame it peaks on (`pdbx_image_id`, + numbered from 1 as the unmerged MTZ's `BATCH`), the scaled intensity and its σ under the error + model, and the scale factor applied (`pdbx_scale_value`). These are the same records as the `.hkl`. + The wwPDB accepts scaled unmerged data beside the merged set. The `_unmerged.mtz` is a different + thing and stays unscaled, for the programs that fit a scale model of their own. + +Items PDBx does not define — ISa, the error model, the twinning and anisotropy diagnostics, the +refined distance and beam centre — are written under the `jfjoch_` local prefix, which a dictionary +check reports as unknown and a deposition ignores. + ## Command-line options General: diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ffaacd498..8f553e584 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -60,6 +60,13 @@ scaled-up counting error gives it, and is in `R_MEAS_WEIGHTED=` at that weight t the number with its noise; a large gap between the two says weak frames are kept at low weight. `MULTIPLICITY=` stays a count of observations. +**`R_MERGE=`** and **`R_PIM=`** are the other two merging R-factors over the same observations as +`R_MEAS=`: R_merge without the multiplicity correction, and the precision-indicating R_pim (Weiss, +2001), which falls as the multiplicity rises. **`CC_STAR=`** is CC\* = sqrt(2 CC1/2 / (1 + CC1/2)) +(Karplus & Diederichs, 2012), the estimated correlation of the merged data with the true signal; it is +absent where CC1/2 is not positive. All three are also columns of the shell tables (`R_merge`, +`R_pim`, `CC*`) and are written per shell into the merged mmCIF. + **`OBSERVATIONS_REJECTED=`** counts the observations the merge's outlier tests removed, and **`OBSERVATIONS_REJECTED_WILSON=`** the part of them the Wilson test took (each listed in the developer report). **`OBSERVATIONS_REJECTED_OVERLOAD=`** is separate from both: the rotation measurements left diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index dd31bae2b..da0fbe3dd 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -68,6 +69,66 @@ std::string CifStr(const std::string& s) { return s; } +std::string Upper(std::string s) { + for (char &c : s) + c = static_cast(std::toupper(static_cast(c))); + return s; +} + +// The general class of the detector, _diffrn_detector.detector in the PDBx enumeration, where its make +// and model name it; empty where they do not. +std::string DetectorClass(const std::string &description) { + const std::string d = Upper(description); + if (d.find("CCD") != std::string::npos) + return "CCD"; + if (d.find("IMAGE PLATE") != std::string::npos) + return "IMAGE PLATE"; + for (const char *pixel : {"EIGER", "PILATUS", "JUNGFRAU"}) + if (d.find(pixel) != std::string::npos) + return "PIXEL"; + return ""; +} + +// The general class of the source, _diffrn_source.source in the PDBx enumeration, from the NXmx +// /entry/source/type the file states; empty where it states none or one PDBx has no class for. +std::string SourceClass(const std::string &nxmx_type) { + const std::string t = Upper(nxmx_type); + if (t.find("SYNCHROTRON") != std::string::npos) return "SYNCHROTRON"; + if (t.find("FREE-ELECTRON") != std::string::npos || t.find("FREE ELECTRON") != std::string::npos) + return "FREE ELECTRON LASER"; + if (t.find("ROTATING ANODE") != std::string::npos) return "ROTATING ANODE"; + if (t.find("FIXED TUBE") != std::string::npos) return "SEALED TUBE"; + if (t.find("METAL JET") != std::string::npos) return "LIQUID ANODE"; + return ""; +} + +// The date part of an ISO-8601 time stamp (YYYY-MM-DD...), which is all the PDBx collection date +// takes; empty for anything else. +std::string IsoDate(const std::string ×tamp) { + if (timestamp.size() < 10 || timestamp[4] != '-' || timestamp[7] != '-') + return ""; + return timestamp.substr(0, 10); +} + +// _entry, _cell and _symmetry: what every data block of the file restates, so each one stands alone. +void WriteBlockHeader(std::ostream &out, const std::string &entry_id, const UnitCell &unitCell, + const gemmi::SpaceGroup &sg) { + out << "_entry.id " << entry_id << "\n"; + out << "#\n"; + out << "_cell.entry_id " << entry_id << "\n"; + out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n"; + out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n"; + out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n"; + out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n"; + out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n"; + out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n"; + out << "#\n"; + out << "_symmetry.entry_id " << entry_id << "\n"; + out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n"; + out << "_symmetry.Int_Tables_number " << sg.number << "\n"; + out << "#\n"; +} + // One output row per reflection in the standard CCP4 anomalous layout: the merged mean (IMEAN / F) // plus the two Bijvoet mates (I(+)/I(-), F(+)/F(-)). The merge keeps the two mates as separate rows // (I+ under the ASU representative hkl, I- under -hkl); this collapses them into one row so the MTZ @@ -162,6 +223,8 @@ void WriteMmcifReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const std::vector &scaled_fulls, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads) { @@ -171,9 +234,16 @@ void WriteMmcifReflections(const std::vector &reflections, out << std::fixed; + // The keys PDBx asks for are written with the one value this file has for each: one entry, one + // crystal, one diffraction experiment (diffrn_id 1, which _pdbx_reflns_twin below also points at) + // and one wavelength. The merged data come first, as the wwPDB expects of a structure-factor file. + const std::string entry_id = "sample"; + const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); + // ---------- data block ---------- out << "data_sample" << "\n"; out << "#\n"; + WriteBlockHeader(out, entry_id, unitCell, sg); // ---------- _audit ---------- out << "_audit.revision_id 1\n"; @@ -182,35 +252,54 @@ void WriteMmcifReflections(const std::vector &reflections, out << "#\n"; // ---------- _software ---------- - out << "_software.name 'Rugnux'\n"; - out << "_software.version " << CifStr(jfjoch_version()) << "\n"; - out << "_software.classification 'data reduction'\n"; + // Rugnux integrated the images (data reduction) and scaled and merged them (data scaling), the + // two rows the wwPDB asks a data-processing program to be named under. + out << "loop_\n"; + out << "_software.pdbx_ordinal\n"; + out << "_software.name\n"; + out << "_software.version\n"; + out << "_software.classification\n"; + out << "1 Rugnux " << CifStr(jfjoch_version()) << " 'data reduction'\n"; + out << "2 Rugnux " << CifStr(jfjoch_version()) << " 'data scaling'\n"; out << "#\n"; - // ---------- _cell ---------- - out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n"; - out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n"; - out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n"; - out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n"; - out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n"; - out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n"; - - const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); - - // ---------- _symmetry ---------- - out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n"; - out << "_symmetry.Int_Tables_number " << sg.number << "\n"; + // ---------- _exptl_crystal / _diffrn ---------- + out << "_exptl_crystal.id 1\n"; + out << "#\n"; + out << "_diffrn.id 1\n"; + out << "_diffrn.crystal_id 1\n"; + if (const auto t = experiment.GetSampleTemperature_K()) + out << "_diffrn.ambient_temp " << Fmt(*t, 1) << "\n"; out << "#\n"; - // ---------- _diffrn_source / _diffrn_detector ---------- + // ---------- _diffrn_radiation / _diffrn_source / _diffrn_detector ---------- + out << "_diffrn_radiation.diffrn_id 1\n"; + out << "_diffrn_radiation.wavelength_id 1\n"; + out << "_diffrn_radiation.pdbx_diffrn_protocol 'SINGLE WAVELENGTH'\n"; + out << "#\n"; + out << "_diffrn_radiation_wavelength.id 1\n"; + out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; + out << "#\n"; + + // The source and the detector as the input file names them: the facility and the beamline (NXmx + // /entry/source/name and /entry/instrument/name) and the detector's make and model. The general + // classes PDBx enumerates are written only where those names settle them. + out << "_diffrn_source.diffrn_id 1\n"; + if (const std::string source = SourceClass(experiment.GetSourceType()); !source.empty()) + out << "_diffrn_source.source " << CifStr(source) << "\n"; if (!experiment.GetSourceName().empty()) out << "_diffrn_source.pdbx_synchrotron_site " << CifStr(experiment.GetSourceName()) << "\n"; - if (!experiment.GetInstrumentName().empty()) out << "_diffrn_source.pdbx_synchrotron_beamline " << CifStr(experiment.GetInstrumentName()) << "\n"; + out << "_diffrn_source.pdbx_wavelength_list " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; + out << "#\n"; - out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; - out << "_diffrn_detector.detector " << CifStr(experiment.GetDetectorDescription()) << "\n"; + out << "_diffrn_detector.diffrn_id 1\n"; + if (const std::string cls = DetectorClass(experiment.GetDetectorDescription()); !cls.empty()) + out << "_diffrn_detector.detector " << CifStr(cls) << "\n"; + out << "_diffrn_detector.type " << CifStr(experiment.GetDetectorDescription()) << "\n"; + if (const std::string date = IsoDate(provenance.collection_date); !date.empty()) + out << "_diffrn_detector.pdbx_collection_date " << date << "\n"; // Detector geometry actually used for integration (refined, when geometry refinement ran - the rotation // two-pass or the stills global refinement update it on experiment_ before the written pass). jfjoch_ // local-data-name items: feedback of the distance / beam centre the data was reduced with. @@ -236,6 +325,9 @@ void WriteMmcifReflections(const std::vector &reflections, // Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a // non-anomalous merge keeps its previous stats block / shell-loop columns unchanged. const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous); + out << "_reflns.entry_id " << entry_id << "\n"; + out << "_reflns.pdbx_ordinal 1\n"; + out << "_reflns.pdbx_diffrn_id 1\n"; out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n"; out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n"; out << "_reflns.number_obs " << ov.unique_reflections << "\n"; @@ -243,8 +335,11 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_redundancy " << Fmt(mult(ov), 2) << "\n"; out << "_reflns.percent_possible_obs " << Fmt(compl_pct(ov), 1) << "\n"; out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n"; + out << "_reflns.pdbx_Rmerge_I_obs " << Fmt(ov.r_merge, 4) << "\n"; out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; + out << "_reflns.pdbx_Rpim_I_all " << Fmt(ov.r_pim, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; + out << "_reflns.pdbx_CC_star " << Fmt(CCStar(ov.cc_half), 4) << "\n"; if (has_anom) out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3) << " # SigAno = <|dano|>/\n"; @@ -360,6 +455,8 @@ void WriteMmcifReflections(const std::vector &reflections, } out << "loop_\n"; + out << "_reflns_shell.pdbx_ordinal\n"; + out << "_reflns_shell.pdbx_diffrn_id\n"; out << "_reflns_shell.d_res_high\n"; out << "_reflns_shell.d_res_low\n"; out << "_reflns_shell.number_measured_obs\n"; @@ -367,17 +464,23 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns_shell.pdbx_redundancy\n"; out << "_reflns_shell.percent_possible_obs\n"; out << "_reflns_shell.meanI_over_sigI_obs\n"; + out << "_reflns_shell.Rmerge_I_obs\n"; out << "_reflns_shell.pdbx_Rrim_I_all\n"; + out << "_reflns_shell.pdbx_Rpim_I_all\n"; out << "_reflns_shell.pdbx_CC_half\n"; + out << "_reflns_shell.pdbx_CC_star\n"; if (has_anom) out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n"; + int ordinal = 0; for (const auto &s : statistics.shells) { if (s.unique_reflections == 0) continue; - out << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " " + out << ++ordinal << " 1 " + << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " " << s.total_observations << " " << s.unique_reflections << " " << Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " " - << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4); + << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_merge, 4) << " " << Fmt(s.r_meas, 4) << " " + << Fmt(s.r_pim, 4) << " " << Fmt(s.cc_half, 4) << " " << Fmt(CCStar(s.cc_half), 4); if (has_anom) out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3); out << "\n"; @@ -454,6 +557,67 @@ void WriteMmcifReflections(const std::vector &reflections, out << "#\n"; out << "# End of reflections\n"; + + // ---------- the scaled unmerged observations, as a second data block ---------- + // The observations the merge above averaged - each full with the scale and the error-model sigma + // the merge gave it, at the index it was measured at - which is what the wwPDB accepts as scaled + // unmerged data. The same records as the .hkl: an observation without a usable sigma is in + // neither. pdbx_image_id is the frame the full peaks on, numbered from 1 as the unmerged MTZ's + // BATCH is. (The unmerged MTZ itself stays unscaled, for the programs that fit their own scales.) + if (!scaled_fulls.empty()) { + out << "data_sample_unmerged\n"; + out << "#\n"; + WriteBlockHeader(out, entry_id, unitCell, sg); + out << "_exptl_crystal.id 1\n"; + out << "#\n"; + out << "_diffrn.id 1\n"; + out << "_diffrn.crystal_id 1\n"; + out << "#\n"; + out << "_diffrn_radiation_wavelength.id 1\n"; + out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; + out << "#\n"; + out << "loop_\n"; + out << "_diffrn_refln.diffrn_id\n"; + out << "_diffrn_refln.id\n"; + out << "_diffrn_refln.pdbx_image_id\n"; + out << "_diffrn_refln.index_h\n"; + out << "_diffrn_refln.index_k\n"; + out << "_diffrn_refln.index_l\n"; + out << "_diffrn_refln.intensity_net\n"; + out << "_diffrn_refln.intensity_sigma\n"; + out << "_diffrn_refln.pdbx_scale_value\n"; + + // Numbered in order over the usable fulls, so the id of each is known before the rows are + // formatted on several workers. + std::vector id(scaled_fulls.size() + 1, 0); + for (size_t i = 0; i < scaled_fulls.size(); ++i) { + const auto &f = scaled_fulls[i]; + id[i + 1] = id[i] + ((std::isfinite(f.I) && std::isfinite(f.sigma) && f.sigma > 0.0f) ? 1 : 0); + } + const size_t nrow = scaled_fulls.size(); + const size_t nw = std::max(nthreads, 1); + const int nch = static_cast(ThreadsForWork(nrow, nw, 4096)); + std::vector block(nch); + ParallelChunks(nch, nw, [&](int tlo, int thi) { + for (int t = tlo; t < thi; ++t) { + const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch; + std::string &s = block[t]; + s.reserve((hi - lo) * 64); + for (size_t i = lo; i < hi; ++i) { + if (id[i + 1] == id[i]) + continue; + const auto &f = scaled_fulls[i]; + const long frame = std::isfinite(f.image_number) ? 1 + std::lround(f.image_number) : 0; + s += "1 " + std::to_string(id[i + 1]) + " " + (frame > 0 ? std::to_string(frame) : "?") + " " + + std::to_string(f.h) + " " + std::to_string(f.k) + " " + std::to_string(f.l) + " " + + Fmt(f.I, 3) + " " + Fmt(f.sigma, 3) + " " + Fmt(f.scale, 5) + "\n"; + } + } + }); + for (const std::string &s : block) + out.write(s.data(), static_cast(s.size())); + out << "#\n"; + } out.close(); } @@ -1072,6 +1236,7 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads, const std::vector *scaled_fulls) { @@ -1081,6 +1246,7 @@ void WriteReflections(const std::vector &reflections, // fulls where the merge provides them (rotation), the merged reflections otherwise (stills). WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, + scaled_fulls ? *scaled_fulls : std::vector{}, provenance, filename + ".cif", nthreads); if (scaled_fulls && !scaled_fulls->empty()) WriteShelxHklReflections(*scaled_fulls, filename + ".hkl", nthreads); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 697cbe87c..bdac24f65 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -27,13 +27,24 @@ struct ErrorModelReport { std::string b; }; -// nthreads: workers for the per-reflection row formatting, which is the bulk of the file. +// What the reflection files record about where they came from, beyond the experiment itself. +// Empty strings are left out. +struct OutputProvenance { + std::string collection_date; // as the input states it (NXmx /entry/start_time), ISO-8601 +}; + +// The merged reflections, their statistics and the experiment, as one PDBx/mmCIF data block; and, when +// scaled_fulls holds any, the same observations unmerged and scaled (_diffrn_refln) as a second block - +// the scaled unmerged data the wwPDB accepts beside the merged set. nthreads: workers for the +// per-reflection row formatting, which is the bulk of the file. void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const std::vector &scaled_fulls, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads); @@ -92,6 +103,7 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads, const std::vector *scaled_fulls = nullptr); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 4c63aa15c..af996e7d6 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -669,15 +669,23 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m } // R_meas per shell: sum over reflections of sqrt(n/(n-1)) * sum_i|I_i - |, over sum of I_i. + // R_merge and R_pim beside it, with the factor 1 and sqrt(1/(n-1)) - R_pim following Weiss (2001) + // J. Appl. Cryst. 34, 130-135. std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); - double rmeas_num_all = 0.0, rmeas_den_all = 0.0; + std::vector rmerge_num(n_shells, 0.0), rpim_num(n_shells, 0.0); + double rmeas_num_all = 0.0, rmeas_den_all = 0.0, rmerge_num_all = 0.0, rpim_num_all = 0.0; for (const auto &[key, ra]: rmeas_obs) { if (ra.n < 2 || ra.shell < 0 || ra.shell >= n_shells) continue; const double factor = std::sqrt(static_cast(ra.n) / (ra.n - 1)); + const double pim_factor = std::sqrt(1.0 / (ra.n - 1)); rmeas_num[ra.shell] += factor * ra.sum_abs_dev; + rmerge_num[ra.shell] += ra.sum_abs_dev; + rpim_num[ra.shell] += pim_factor * ra.sum_abs_dev; rmeas_den[ra.shell] += ra.sum_I; rmeas_num_all += factor * ra.sum_abs_dev; + rmerge_num_all += ra.sum_abs_dev; + rpim_num_all += pim_factor * ra.sum_abs_dev; rmeas_den_all += ra.sum_I; } @@ -701,6 +709,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m ss.cc_half = sa.cc_half.GetCC(); ss.cc_ref = sa.cc_ref.GetCC(); ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; + ss.r_merge = rmeas_den[s] > 0.0 ? rmerge_num[s] / rmeas_den[s] : NAN; + ss.r_pim = rmeas_den[s] > 0.0 ? rpim_num[s] / rmeas_den[s] : NAN; } auto &overall = out.overall; @@ -728,6 +738,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m overall.cc_half = cc_half_overall.GetCC(); overall.cc_ref = cc_ref_overall.GetCC(); overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; + overall.r_merge = rmeas_den_all > 0.0 ? rmerge_num_all / rmeas_den_all : NAN; + overall.r_pim = rmeas_den_all > 0.0 ? rpim_num_all / rmeas_den_all : NAN; return out; } @@ -798,15 +810,18 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { double multiplicity = in.unique_reflections > 0 ? static_cast(in.total_observations) / in.unique_reflections : NAN; - output << fmt::format("{:8d} {:8d} {:8d} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + output << fmt::format("{:8d} {:8d} {:8d} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", in.total_observations, in.unique_reflections, in.possible_unique_reflections, Measured(completeness, "%"), Measured(multiplicity), Measured(in.mean_i_over_sigma), + Measured(in.r_merge*100.0, "%"), Measured(in.r_meas*100.0, "%"), + Measured(in.r_pim*100.0, "%"), Measured(in.cc_half*100.0, "%"), + Measured(CCStar(in.cc_half)*100.0, "%"), Measured(in.cc_ref*100.0, "%"), Measured(in.cc_anom*100.0, "%"), Measured(in.abs_diff_over_sigma_anomalous, "", 2)); @@ -814,14 +829,14 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { } std::string ShellTableHeader() { - return fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", - "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref", - "CCanom", "SigAno"); + return fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_merge", "R_meas", "R_pim", + "CC1/2", "CC*", "CCref", "CCanom", "SigAno"); } std::string ShellTableRule() { - return fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "", "", ""); + return fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", + "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""); } std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 6119ab318..9ab944876 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -39,6 +39,12 @@ struct MergeStatisticsShell { // error model (WeightedRmeasTerms); with equal sigmas it is r_meas. A frame kept at a small weight // then no longer sets the number with its noise. NaN where it is not computed (stills, per hand). double r_meas_weighted = NAN; + // The other two merging R-factors, over the same observations and means as r_meas and differing + // from it only in the per-reflection factor on sum_i|I_i-|: 1 for R_merge, sqrt(1/(n-1)) for + // the precision-indicating R_pim (Weiss (2001) J. Appl. Cryst. 34, 130-135). Reported for + // deposition, which asks for all three. + double r_merge = NAN; + double r_pim = NAN; // Anomalous signal-to-noise (XDS "SigAno" / mmCIF pdbx_absDiff_over_sigma_anomalous): // <|I(+)-I(-)|> / over acentric reflections measured in both hands. NaN when @@ -207,6 +213,12 @@ struct MergeStatistics { }; +// CC* = sqrt(2 CC1/2 / (1 + CC1/2)), the correlation of the merged data with the noise-free signal +// (Karplus & Diederichs (2012) Science 336, 1030-1033). NaN where CC1/2 is not positive. +inline double CCStar(double cc_half) { + return cc_half > 0.0 ? std::sqrt(2.0 * cc_half / (1.0 + cc_half)) : NAN; +} + std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in); std::ostream &operator<<(std::ostream &output, const MergeStatistics &in); // The header and the rule of the table operator<< prints, for a table with the same columns rendered diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 57e35e070..8e6c9c39b 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -5409,7 +5409,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (int j = lo; j < hi; ++j) { const Obs &o = fulls[kept[j]]; const float I_corr = o.I * o.corr; - result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr)}; + result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr), + o.corr, o.image_number}; } }); } @@ -5639,6 +5640,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool double rmeas_num_all = 0.0, rmeas_den_all = 0.0; std::vector urmeas_num(n_shells, 0.0), urmeas_den(n_shells, 0.0); double urmeas_num_all = 0.0, urmeas_den_all = 0.0; + std::vector rmerge_num(n_shells, 0.0), rpim_num(n_shells, 0.0); + double rmerge_num_all = 0.0, rpim_num_all = 0.0; // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last // full pass over the observations there is - so, like the correction surfaces, it is not made on a // merge whose numbers nobody reads. @@ -5733,6 +5736,11 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); urmeas_num[r.shell] += factor * r.sum_abs_dev; urmeas_den[r.shell] += r.sum_I; urmeas_num_all += factor * r.sum_abs_dev; urmeas_den_all += r.sum_I; + // R_merge and R_pim on the same sums, unweighted as r_meas is: factor 1, and + // sqrt(1/(n-1)) following Weiss (2001) J. Appl. Cryst. 34, 130-135. + const double pim_factor = std::sqrt(1.0 / (r.n - 1)); + rmerge_num[r.shell] += r.sum_abs_dev; rmerge_num_all += r.sum_abs_dev; + rpim_num[r.shell] += pim_factor * r.sum_abs_dev; rpim_num_all += pim_factor * r.sum_abs_dev; } } @@ -6003,6 +6011,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool ss.cc_ref = NAN; ss.r_meas = urmeas_den[s] > 0.0 ? urmeas_num[s] / urmeas_den[s] : NAN; ss.r_meas_weighted = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; + ss.r_merge = urmeas_den[s] > 0.0 ? rmerge_num[s] / urmeas_den[s] : NAN; + ss.r_pim = urmeas_den[s] > 0.0 ? rpim_num[s] / urmeas_den[s] : NAN; ss.abs_diff_over_sigma_anomalous = sig_den[s] > 0.0 ? sig_num[s] / sig_den[s] : NAN; ss.cc_anom = cc_anom_shell[s].GetCC(); } @@ -6019,6 +6029,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool overall.cc_ref = NAN; overall.r_meas = urmeas_den_all > 0.0 ? urmeas_num_all / urmeas_den_all : NAN; overall.r_meas_weighted = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; + overall.r_merge = urmeas_den_all > 0.0 ? rmerge_num_all / urmeas_den_all : NAN; + overall.r_pim = urmeas_den_all > 0.0 ? rpim_num_all / urmeas_den_all : NAN; overall.abs_diff_over_sigma_anomalous = sig_den_all > 0.0 ? sig_num_all / sig_den_all : NAN; overall.cc_anom = cc_anom_overall.GetCC(); @@ -6039,6 +6051,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool hs.cc_half = ha[sh].cc_half.GetCC(); hs.r_meas = ha[sh].rmeas_den > 0.0 ? ha[sh].rmeas_num / ha[sh].rmeas_den : NAN; hs.r_meas_weighted = NAN; + hs.r_merge = NAN; // not formed per hand + hs.r_pim = NAN; out.per_hand.overall.total_observations += ha[sh].total_obs; out.per_hand.overall.unique_reflections += ha[sh].unique; out.per_hand.overall.possible_unique_reflections += ha[sh].possible; diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 0444a668f..12ba9a265 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -855,6 +855,7 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen InstrumentMetadata metadata; metadata.InstrumentName(master_file->GetString("/entry/instrument/name")); metadata.SourceName(master_file->GetString("/entry/source/name")); + metadata.SourceType(master_file->GetString("/entry/source/type")); dataset->experiment.ImportInstrumentMetadata(metadata); // The rotation axis is whatever the file calls it. The name is free-form throughout the API, diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 7a752dae1..2ac62f675 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -995,8 +995,14 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("R_MEAS", o.r_meas, "{:.4f}")); if (std::isfinite(o.r_meas_weighted)) Add(s, KeyReal("R_MEAS_WEIGHTED", o.r_meas_weighted, "{:.4f}")); + if (std::isfinite(o.r_merge)) + Add(s, KeyReal("R_MERGE", o.r_merge, "{:.4f}")); + if (std::isfinite(o.r_pim)) + Add(s, KeyReal("R_PIM", o.r_pim, "{:.4f}")); if (std::isfinite(o.cc_half)) Add(s, KeyReal("CC_HALF", o.cc_half, "{:.4f}")); + if (std::isfinite(CCStar(o.cc_half))) + Add(s, KeyReal("CC_STAR", CCStar(o.cc_half), "{:.4f}")); if (std::isfinite(o.abs_diff_over_sigma_anomalous)) Add(s, KeyReal("SIGANO", o.abs_diff_over_sigma_anomalous, "{:.3f}")); if (std::isfinite(o.cc_anom)) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4df80b6b8..51b284640 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -40,6 +40,7 @@ #include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" // BeamCenterEstimate #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult #include "../image_analysis/geom_refinement/PowderCalibration.h" // CalibrationMethod, CalibrationResult +#include "../image_analysis/WriteReflections.h" // OutputProvenance class JFJochReader; class JFJochHDF5Reader; @@ -70,6 +71,8 @@ struct ProcessConfig { // 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). + 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. diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 823d15df6..e00a0f95a 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3477,7 +3477,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { f.h = h[0]; f.k = h[1]; f.l = h[2]; } WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, - em_report, result.twinning, config_.output_prefix, + em_report, result.twinning, config_.provenance, config_.output_prefix, static_cast(std::max(1, config_.nthreads)), &sm.scaled_fulls); writing_mark.Stop(); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9e74ef787..003f59bde 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1664,6 +1664,9 @@ static int RunRugnux(int argc, char **argv) { } logger.Info("Loaded dataset from {}", input_file); + // Recorded in the reflection files. + OutputProvenance output_provenance; + output_provenance.collection_date = dataset->arm_date; { // The header as the file states it, before anything is measured: what a question about a run // is answered from first. @@ -2229,7 +2232,8 @@ static int RunRugnux(int argc, char **argv) { error_model_a > 0 ? fmt::format("{:.3f}", error_model_a) : std::string(), error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, - em_report, twinning, output_prefix, static_cast(nthreads), &scaled_fulls); + em_report, twinning, output_provenance, output_prefix, static_cast(nthreads), + &scaled_fulls); } // --mode scale re-merges stored reflections, so it determines a merging result and gets the @@ -2986,6 +2990,7 @@ static int RunRugnux(int argc, char **argv) { // Rugnux; the experiment above carries all algorithm settings). ProcessConfig config; config.mode = ProcessMode::FullAnalysis; + config.provenance = output_provenance; config.start_image = start_image; config.end_image = end_image; config.stride = image_stride; diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp index 7f3e070ab..9c95b7c07 100644 --- a/tests/WriteReflectionsTest.cpp +++ b/tests/WriteReflectionsTest.cpp @@ -9,10 +9,14 @@ #include #include +#include +#include #include "../common/DiffractionExperiment.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/IntegrationOutcome.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" #include #include "SyntheticMergedReflections.h" @@ -217,3 +221,85 @@ TEST_CASE("SHELX .hkl of unmerged fulls: fixed 3I4,2F8.2, scaled to fit, 0 0 0 t CHECK(lines[2] == " 1 2 3 -2.50 1.00"); CHECK(lines[3] == " 0 0 0 0.00 0.00"); } + +TEST_CASE("CC* from CC1/2", "[write_reflections][portable]") { + CHECK(CCStar(1.0) == Catch::Approx(1.0)); + CHECK(CCStar(0.5) == Catch::Approx(std::sqrt(2.0 / 3.0))); + CHECK(std::isnan(CCStar(0.0))); + CHECK(std::isnan(CCStar(-0.1))); + CHECK(std::isnan(CCStar(NAN))); +} + +TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged block", + "[write_reflections][portable]") { + jfjoch_test::SyntheticMergeParams params; + params.true_space_group = "P 43 21 2"; + params.twin_supergroup = "P 43 21 2"; + params.d_min_A = 5.0; + const auto reflections = jfjoch_test::GenerateSyntheticMerged(params); + REQUIRE(!reflections.empty()); + + auto experiment = TestExperiment(); + experiment.SampleTemperature_K(100.0f); + + MergeStatistics statistics; + MergeStatisticsShell shell; + shell.d_min = 5.0f; shell.d_max = 30.0f; + shell.total_observations = 400; shell.unique_reflections = 100; shell.possible_unique_reflections = 110; + shell.mean_i_over_sigma = 10.0; + shell.cc_half = 0.98; + shell.r_merge = 0.05; shell.r_meas = 0.06; shell.r_pim = 0.03; + statistics.shells = {shell}; + statistics.overall = shell; + + // Three fulls, one without a usable sigma - that one is left out, as from the .hkl. + const std::vector fulls{ + {1, 2, 3, 200.0f, 20.0f, 1.25f, 0.4f}, + {-1, -2, -3, 100.0f, 12.0f, 0.8f, 9.6f}, + {4, 5, 6, 10.0f, 0.0f, 1.0f, 3.0f}, + }; + OutputProvenance provenance; + provenance.collection_date = "2026-07-13T13:48:52.413Z"; + + const auto path = (std::filesystem::temp_directory_path() / "rugnux_merged_test.cif").string(); + WriteMmcifReflections(reflections, TETRAGONAL_CELL, experiment, statistics, ErrorModelReport{}, + TwinningAnalysisResult{}, fulls, provenance, path, 2); + gemmi::cif::Document doc = gemmi::read_cif_gz(path); + std::filesystem::remove(path); + + REQUIRE(doc.blocks.size() == 2); + gemmi::cif::Block &merged = doc.blocks[0]; + CHECK(merged.name == "sample"); + CHECK(*merged.find_value("_reflns.pdbx_ordinal") == "1"); + CHECK(*merged.find_value("_reflns.pdbx_diffrn_id") == "1"); + CHECK(*merged.find_value("_reflns.pdbx_Rmerge_I_obs") == "0.0500"); + CHECK(*merged.find_value("_reflns.pdbx_Rpim_I_all") == "0.0300"); + CHECK(gemmi::cif::as_number(*merged.find_value("_reflns.pdbx_CC_star")) + == Catch::Approx(std::sqrt(2 * 0.98 / 1.98)).margin(1e-4)); + CHECK(*merged.find_value("_diffrn.ambient_temp") == "100.0"); + CHECK(*merged.find_value("_diffrn_detector.pdbx_collection_date") == "2026-07-13"); + + gemmi::cif::Table software = merged.find("_software.", {"pdbx_ordinal", "classification"}); + REQUIRE(software.length() == 2); + CHECK(gemmi::cif::as_string(software[0][1]) == "data reduction"); + CHECK(gemmi::cif::as_string(software[1][1]) == "data scaling"); + + gemmi::cif::Table shells = merged.find("_reflns_shell.", {"pdbx_ordinal", "Rmerge_I_obs", + "pdbx_Rpim_I_all", "pdbx_CC_star"}); + REQUIRE(shells.length() == 1); + CHECK(shells[0][1] == "0.0500"); + CHECK(shells[0][2] == "0.0300"); + + gemmi::cif::Block &unmerged = doc.blocks[1]; + CHECK(unmerged.name == "sample_unmerged"); + gemmi::cif::Table obs = unmerged.find("_diffrn_refln.", {"id", "pdbx_image_id", "index_h", + "intensity_net", "intensity_sigma", + "pdbx_scale_value"}); + REQUIRE(obs.length() == 2); + CHECK(obs[0][0] == "1"); + CHECK(obs[0][1] == "1"); // peaks on frame 0.4 -> image 1 + CHECK(obs[1][1] == "11"); // 9.6 -> image 11 + CHECK(obs[1][2] == "-1"); + CHECK(gemmi::cif::as_number(obs[0][3]) == Catch::Approx(200.0)); + CHECK(gemmi::cif::as_number(obs[1][5]) == Catch::Approx(0.8)); +}