From ce12d697518de200b79e795dc7ad572c9b1df627 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 16 Jul 2026 10:37:56 +0200 Subject: [PATCH] rugnux: write standard CCP4 anomalous MTZ (IMEAN + I(+)/I(-)) In anomalous mode the merge keeps the two Friedel mates as separate rows, which WriteMtzReflections emitted verbatim - two IMEAN rows per Bijvoet pair that downstream tools had to re-collapse. Pair the mates into one row per reflection with the standard CCP4 anomalous layout: IMEAN + I(+)/I(-) and the matching F/F(+)/F(-) amplitudes. A single HKLKeyGenerator(merge_friedel=false) yields both the shared ASU group key and which mate a row is (.plus); IMEAN/F are the inverse-variance Friedel mean, centrics/unpaired keep one hand with the other missing. aimless/mtz2sca/ANODE now read the file directly. The non-anomalous MTZ path is unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/CHANGELOG.md | 1 + image_analysis/WriteReflections.cpp | 121 +++++++++++++++++++++++----- 2 files changed, 102 insertions(+), 20 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 472cee57..1f71828d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -20,6 +20,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * jfjoch_writer: Remove the CBF and TIFF image writers - only NXmx HDF5 is written now (all three HDF5 layouts remain: legacy, VDS, integrated). * Reader: Treat a negative `total_flux` in a stored dataset as unknown/absent rather than a real (invalid) flux. * Packaging: Build the self-contained Linux viewer against a static libdbus with glib disabled, so the `.tgz` no longer drags a tail of system `.so` dependencies; add parallel image-build and in-container viewer-verification scripts. +* rugnux: Write anomalous data as a standard CCP4 anomalous MTZ - one row per reflection with `IMEAN`, `I(+)`/`I(-)` (and the matching `F`/`F(+)`/`F(-)` amplitudes) - instead of two `IMEAN` rows per Bijvoet pair, so `aimless`/`mtz2sca`/ANODE read it directly. The non-anomalous MTZ output is unchanged. ### 1.0.0-rc.159 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index b01bc0a7..e617d674 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -3,9 +3,12 @@ #include "WriteReflections.h" #include "scale_merge/Merge.h" +#include "scale_merge/HKLKey.h" #include "scale_merge/TwinningAnalysis.h" #include +#include +#include #include #include #include @@ -220,36 +223,114 @@ void WriteMtzReflections(const std::vector &reflections, ds.crystal_name = experiment.GetSampleName(); ds.wavelength = experiment.GetWavelength_A(); - int dataset_id = ds.id; + const int dataset_id = ds.id; + + // In anomalous mode the merge keeps the two Friedel mates as separate rows (I+ under the ASU + // representative hkl, I- under -hkl). Emitting those verbatim gives a file with two rows per + // reflection that downstream tools have to re-collapse. Instead pair the mates into one row per + // reflection with the standard CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for + // the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly. + if (experiment.GetScalingSettings().GetMergeFriedel()) { + mtz.add_column("H", 'H', dataset_id, -1, false); + mtz.add_column("K", 'H', dataset_id, -1, false); + mtz.add_column("L", 'H', dataset_id, -1, false); + mtz.add_column("IMEAN", 'J', dataset_id, -1, false); + mtz.add_column("SIGIMEAN", 'Q', 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); + mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false); + + mtz.nreflections = static_cast(reflections.size()); + mtz.data.reserve(reflections.size() * 8); + for (const auto& r : reflections) { + mtz.data.push_back(static_cast(r.h)); + mtz.data.push_back(static_cast(r.k)); + mtz.data.push_back(static_cast(r.l)); + mtz.data.push_back(r.I); + mtz.data.push_back(r.sigma); + mtz.data.push_back(r.F); + mtz.data.push_back(r.sigmaF); + mtz.data.push_back(r.rfree_flag ? 1.0f : 0.0f); + } + mtz.write_to_file(filename); + return; + } + + // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative. A single + // generator gives both the group key (its hkl, identical for +hkl and -hkl) and which mate this + // row is (.plus). + const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1)); + + struct AnomRow { + int h = 0, k = 0, l = 0; + float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; + float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN; + int rfree = 0; + }; + std::map, AnomRow> rows; + for (const auto& r : reflections) { + const HKLKey key = key_gen(r); + AnomRow& row = rows[{key.h, key.k, key.l}]; + row.h = key.h; row.k = key.k; row.l = key.l; + row.rfree = r.rfree_flag ? 1 : 0; + if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; } + else { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; } + } + + // Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured). + const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) { + const bool ok_a = std::isfinite(a) && sa > 0.0f; + const bool ok_b = std::isfinite(b) && sb > 0.0f; + if (ok_a && ok_b) { + const double wa = 1.0 / (static_cast(sa) * sa); + const double wb = 1.0 / (static_cast(sb) * sb); + val = static_cast((wa * a + wb * b) / (wa + wb)); + sig = static_cast(1.0 / std::sqrt(wa + wb)); + } else if (ok_a) { val = a; sig = sa; } + else if (ok_b) { val = b; sig = sb; } + else { val = NAN; sig = NAN; } + }; - // Add columns mtz.add_column("H", 'H', dataset_id, -1, false); mtz.add_column("K", 'H', dataset_id, -1, false); mtz.add_column("L", 'H', dataset_id, -1, false); mtz.add_column("IMEAN", 'J', dataset_id, -1, false); mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false); - mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude + 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 (mean) mtz.add_column("SIGF", 'Q', 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("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); - // Number of rows - mtz.nreflections = reflections.size(); - - // Allocate data table - mtz.data.reserve(reflections.size() * 8); - - for (const auto& r : reflections) { - mtz.data.push_back(static_cast(r.h)); - mtz.data.push_back(static_cast(r.k)); - mtz.data.push_back(static_cast(r.l)); - mtz.data.push_back(r.I); - mtz.data.push_back(r.sigma); - mtz.data.push_back(r.F); - mtz.data.push_back(r.sigmaF); - mtz.data.push_back(r.rfree_flag ? 1 : 0); + mtz.nreflections = static_cast(rows.size()); + mtz.data.reserve(rows.size() * 16); + for (const auto& [hkl, row] : rows) { + float i_mean, sig_i_mean, f_mean, sig_f_mean; + combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean); + combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean); + 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(i_mean); + mtz.data.push_back(sig_i_mean); + 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(f_mean); + mtz.data.push_back(sig_f_mean); + 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)); } - - // Write MTZ mtz.write_to_file(filename); }