diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index d7c49114..953c5034 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -4,7 +4,7 @@ data-processing half of the system (see [Naming](NAMING.md) for where the name comes from). It takes an existing HDF5 dataset, runs the full analysis pipeline — spot finding, indexing, geometry refinement, Bragg integration and (optionally) scaling and merging — and writes the -results to a `_process.h5` file, plus reflection files (`.mtz`/`.cif`) when merging is +results to a `_process.h5` file, plus reflection files (`.mtz`/`.cif`/`.hkl`) when merging is requested. It runs the *same* analysis code as the online and interactive tools, just driven from the @@ -20,7 +20,7 @@ command line over a file rather than a live detector stream. | --- | --- | --- | --- | | [`jfjoch_broker`](JFJOCH_BROKER.md) | Online, real-time streaming analysis on FPGA + GPU | HTTP/REST + ZeroMQ | Live results and statistics, images streamed to [`jfjoch_writer`](JFJOCH_WRITER.md) | | [`jfjoch_viewer`](JFJOCH_VIEWER.md) | Interactive, on-screen exploration | Qt desktop application | Displayed on screen (results not saved to disk) | -| **`rugnux`** | **Offline batch processing of a stored dataset** | **Command-line interface** | **`_process.h5`, and `.mtz`/`.cif` when merging** | +| **`rugnux`** | **Offline batch processing of a stored dataset** | **Command-line interface** | **`_process.h5`, and `.mtz`/`.cif`/`.hkl` when merging** | Use `rugnux` to re-analyse data after acquisition, to experiment with processing parameters, or to produce merged intensities for downstream structure solution. @@ -44,11 +44,21 @@ the first pass. integration, azimuthal integration, per-image statistics). See [HDF5 / NeXus data format](HDF5.md) for the layout. Written by default only when **not** merging (i.e. under `--no-merge`); add `--write-process-h5` to also write it when merging. -- Merging is **on by default** (`--no-merge` disables it). The merged reflections are written as - **both** `.mtz` and `.cif` (mmCIF) — each has its uses downstream (MTZ for the - CCP4 / phenix tools, mmCIF for deposition). Both carry the **refined unit cell** (from rotation - indexing) and the **space group determined from systematic absences** (constrained to the indexed - lattice symmetry). No-reference scaling additionally emits per-iteration `_iterN_scale.dat`. +- Merging is **on by default** (`--no-merge` disables it). The merged reflections are written in + **three** formats — each has its uses downstream: + - `.mtz` — CCP4 MTZ (`IMEAN`/`I(+)`/`I(-)`, French–Wilson `F`, `FreeR_flag`) for the CCP4 / + phenix reflection tools. + - `.cif` — mmCIF, for deposition and as the self-describing native format (also carries the + merging statistics, ISa, twinning and radiation-damage indicators). + - `.hkl` — SHELX **HKLF 4** text (`h k l I σ(I)`, fixed `3I4,2F8.2`), the direct input for + **SHELXC / ANODE / SHELXD**. Bijvoet mates are written separately (`I(+)` at `+hkl`, `I(-)` at + `-hkl`) so the anomalous signal is preserved; intensities are put on a common scale so the largest + value fits the fixed-width field (the absolute scale is irrelevant to SHELXC/ANODE), and the file + ends with the `0 0 0` terminator record. + + All three carry the **refined unit cell** (from rotation indexing) and the **space group determined + from systematic absences** (constrained to the indexed lattice symmetry). No-reference scaling + additionally emits per-iteration `_iterN_scale.dat`. Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the console. By default the written resolution is trimmed automatically where CC1/2 falls off diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 18d3d6cd..610fdb3a 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -59,6 +59,84 @@ std::string CifStr(const std::string& s) { return s; } +// 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 +// and SHELX writers share one row list. has_anom is set false when NO reflection carries an anomalous +// split (e.g. the stills path) - callers then omit the +/- columns. +struct MergedOutRow { + int h = 0, k = 0, l = 0; + float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; + float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN; + int rfree = 0; +}; + +std::vector BuildMergedRows(const std::vector &reflections, + const DiffractionExperiment &experiment, + bool &has_anom) { + std::vector out_rows; + has_anom = true; + + if (experiment.GetScalingSettings().GetMergeFriedel()) { + // Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim + // from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak + // anomalous signal is preserved without reprocessing. A reflection with only one mate, or a + // centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous + // split (e.g. the stills path, which does not compute one) the anomalous columns are omitted. + has_anom = std::any_of(reflections.begin(), reflections.end(), + [](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); }); + out_rows.reserve(reflections.size()); + for (const auto& r : reflections) + out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus, + r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus, + r.sigmaF_minus, r.rfree_flag ? 1 : 0}); + } else { + // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form + // IMEAN / F as their inverse-variance Friedel mean. 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; } + }; + + out_rows.reserve(rows.size()); + 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); + out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im, + row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree}); + } + } + return out_rows; +} + } // namespace @@ -286,77 +364,10 @@ void WriteMtzReflections(const std::vector &reflections, const int dataset_id = ds.id; - // The merge keeps the two Friedel mates as separate rows (I+ under the ASU representative hkl, I- - // under -hkl). Both output modes collapse them 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. The modes differ only in how each row is - // formed, so they build one common row list and share the column/row emission below. - struct OutRow { - int h = 0, k = 0, l = 0; - float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN; - float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN; - int rfree = 0; - }; - std::vector out_rows; + // 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; - - if (experiment.GetScalingSettings().GetMergeFriedel()) { - // Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim - // from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak - // anomalous signal is preserved without reprocessing. A reflection with only one mate, or a - // centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous - // split (e.g. the stills path, which does not compute one) the anomalous columns are omitted. - has_anom = std::any_of(reflections.begin(), reflections.end(), - [](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); }); - out_rows.reserve(reflections.size()); - for (const auto& r : reflections) - out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus, - r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus, - r.sigmaF_minus, r.rfree_flag ? 1 : 0}); - } else { - // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form - // IMEAN / F as their inverse-variance Friedel mean. 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; } - }; - - out_rows.reserve(rows.size()); - 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); - out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im, - row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree}); - } - } + const std::vector out_rows = BuildMergedRows(reflections, experiment, has_anom); mtz.add_column("H", 'H', dataset_id, -1, false); mtz.add_column("K", 'H', dataset_id, -1, false); @@ -406,6 +417,49 @@ void WriteMtzReflections(const std::vector &reflections, mtz.write_to_file(filename); } +void WriteShelxHklReflections(const std::vector &reflections, + const DiffractionExperiment &experiment, + const std::string &filename) { + bool has_anom = true; + const std::vector rows = BuildMergedRows(reflections, experiment, has_anom); + + // SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as + // h k l I sigma(I). The Bijvoet mates are written separately - I(+) at +hkl, I(-) at -hkl - so the + // anomalous differences survive; a reflection with no anomalous split is written once as its mean. + // Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale + // is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the + // anomalous signal is untouched. The file ends with a 0 0 0 terminator record. + const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; }; + + double max_abs = 0.0; + for (const auto& r : rows) { + if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)}); + if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)}); + if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean)) + max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)}); + } + const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0; + + std::ofstream out(filename); + if (!out) + throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename); + out << std::fixed << std::setprecision(2); + const auto emit = [&out, scale](int h, int k, int l, float I, float sigma) { + out << std::setw(4) << h << std::setw(4) << k << std::setw(4) << l + << std::setw(8) << scale * I << std::setw(8) << scale * sigma << "\n"; + }; + for (const auto& r : rows) { + const bool plus = usable(r.Ip, r.sIp); + const bool minus = usable(r.Im, r.sIm); + if (plus) emit(r.h, r.k, r.l, r.Ip, r.sIp); + if (minus) emit(-r.h, -r.k, -r.l, r.Im, r.sIm); + if (!plus && !minus && usable(r.Imean, r.sImean)) + emit(r.h, r.k, r.l, r.Imean, r.sImean); + } + emit(0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker + out.close(); +} + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, @@ -413,8 +467,10 @@ void WriteReflections(const std::vector &reflections, const std::string &isa, const TwinningAnalysisResult &twinning, const std::string &filename) { - // Always write both an MTZ and an mmCIF - each has its uses downstream (MTZ for the CCP4 / - // phenix reflection tools, mmCIF for deposition and as the self-describing native format). + // Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / + // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as + // the SHELXC / ANODE substructure-solution input). WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, twinning, filename + ".cif"); + WriteShelxHklReflections(reflections, experiment, filename + ".hkl"); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 5cdc68a2..c6537e81 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -26,6 +26,11 @@ void WriteMtzReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const std::string &filename); +// SHELX HKLF-4 text file (h k l I sigma(I), Bijvoet mates separate) for SHELXC / ANODE. +void WriteShelxHklReflections(const std::vector &reflections, + const DiffractionExperiment &experiment, + const std::string &filename); + void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment,