Jungfraujoch/image_analysis/WriteReflections.cpp

// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only

#include "WriteReflections.h"
#include "scale_merge/Merge.h"
#include "scale_merge/HKLKey.h"
#include "scale_merge/TwinningAnalysis.h"
#include "../common/ParallelFor.h"

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <map>
#include <tuple>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <ctime>
#include <chrono>

#include <gemmi/mtz.hpp>

#include "../common/GitInfo.h"

namespace {

/// Current date in ISO-8601 (YYYY-MM-DD) for the _audit block.
std::string CurrentDateISO() {
    auto now = std::chrono::system_clock::now();
    auto t   = std::chrono::system_clock::to_time_t(now);
    std::tm tm{};
#ifdef _WIN32
    gmtime_s(&tm, &t);
#else
    gmtime_r(&t, &tm);
#endif
    char buf[32];
    std::strftime(buf, sizeof(buf), "%Y-%m-%d", &tm);
    return buf;
}

/// Format a double with given decimal places; returns "?" for non-finite.
/// snprintf rather than an ostringstream: the reflection loop below calls this twelve times per
/// reflection, and building a stream (and its locale) per call dominated the time spent writing a
/// merged file. Same digits - both go through the C locale's %.*f.
std::string Fmt(double val, int decimals = 4) {
    if (!std::isfinite(val))
        return "?";
    char buf[512];
    const int n = std::snprintf(buf, sizeof(buf), "%.*f", decimals, val);
    return std::string(buf, n);
}

/// Quote a CIF string value; returns "?" for empty.
std::string CifStr(const std::string& s) {
    if (s.empty())
        return "?";
    // If it contains spaces or special chars, single-quote it
    if (s.find(' ') != std::string::npos ||
        s.find('\'') != std::string::npos ||
        s.find('#') != std::string::npos)
        return "'" + 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<MergedOutRow> BuildMergedRows(const std::vector<MergedReflection> &reflections,
                                          const DiffractionExperiment &experiment,
                                          bool &has_anom) {
    std::vector<MergedOutRow> 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<std::tuple<int, int, int>, 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<double>(sa) * sa);
                const double wb = 1.0 / (static_cast<double>(sb) * sb);
                val = static_cast<float>((wa * a + wb * b) / (wa + wb));
                sig = static_cast<float>(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


void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
                           const UnitCell &unitCell,
                           const DiffractionExperiment &experiment,
                           const MergeStatistics &statistics,
                           const ErrorModelReport &error_model,
                           const TwinningAnalysisResult &twinning,
                           const std::string &filename,
                           size_t nthreads) {

    std::ofstream out(filename);
    if (!out)
        throw std::runtime_error("WriteMmcifReflections: cannot open " + filename);

    out << std::fixed;

    // ---------- data block ----------
    out << "data_sample" << "\n";
    out << "#\n";

    // ---------- _audit ----------
    out << "_audit.revision_id     1\n";
    out << "_audit.creation_date   " << CurrentDateISO() << "\n";
    out << "_audit.update_record   'Initial release'\n";
    out << "#\n";

    // ---------- _software ----------
    out << "_software.name         'Rugnux'\n";
    out << "_software.version      " << CifStr(jfjoch_version()) << "\n";
    out << "_software.classification 'data reduction'\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";

    auto *sg = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1));
    if (sg == nullptr)
        throw std::runtime_error("WriteMmcifReflections: invalid space group number");

    // ---------- _symmetry ----------
    out << "_symmetry.space_group_name_H-M  " << CifStr(sg->hm) << "\n";
    out << "_symmetry.Int_Tables_number      " << sg->number << "\n";
    out << "#\n";

    // ---------- _diffrn_source / _diffrn_detector ----------
    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_radiation_wavelength.wavelength  " << Fmt(experiment.GetWavelength_A(), 5) << "\n";
    out << "_diffrn_detector.detector  " << CifStr(experiment.GetDetectorDescription()) << "\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.
    out << "_diffrn_detector.jfjoch_distance_mm       " << Fmt(experiment.GetDetectorDistance_mm(), 4) << "\n";
    out << "_diffrn_detector.jfjoch_beam_center_x_pxl " << Fmt(experiment.GetBeamX_pxl(), 2) << "\n";
    out << "_diffrn_detector.jfjoch_beam_center_y_pxl " << Fmt(experiment.GetBeamY_pxl(), 2) << "\n";
    out << "#\n";

    // ---------- merging statistics (_reflns overall + _reflns_shell loop) ----------
    // cc_half and r_meas are stored as fractions (0-1), which is the mmCIF convention. ISa (the
    // Diederichs asymptotic I/sigma, 1/b of the a*sigma^2 + (b*I)^2 error model) and the twinning
    // indicators below have no standard mmCIF item. They are written under the "jfjoch" reserved
    // prefix (_reflns.jfjoch_*), the IUCr-sanctioned local-data-name extension for private items -
    // NOT the "pdbx_" prefix, which is owned by the wwPDB PDBx/mmCIF dictionary and must not label
    // items that dictionary does not define. (The other pdbx_ items here are genuine PDBx items.)
    const auto mult = [](const MergeStatisticsShell &s) {
        return s.unique_reflections > 0 ? static_cast<double>(s.total_observations) / s.unique_reflections : 0.0; };
    const auto compl_pct = [](const MergeStatisticsShell &s) {
        return s.possible_unique_reflections > 0
               ? 100.0 * static_cast<double>(s.unique_reflections) / s.possible_unique_reflections : 0.0; };
    if (!statistics.shells.empty()) {
        const auto &ov = statistics.overall;
        // 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.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";
        out << "_reflns.pdbx_number_measured_all  " << ov.total_observations << "\n";
        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_Rrim_I_all           " << Fmt(ov.r_meas, 4) << "\n";
        out << "_reflns.pdbx_CC_half              " << Fmt(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|>/<sigma(dano)>\n";
        // ISa in XDS's sense: the whole-range 1/sqrt(a*b) of the error model below, so this item can
        // be read straight against a CORRECT.LP. The strong-reflection asymptote - a tier XDS does not
        // have, and always the more optimistic of the two - is written separately rather than here.
        out << "_reflns.jfjoch_diffrn_ISa             " << CifStr(error_model.isa)
            << "   # 1/sqrt(a*b), the XDS convention\n";
        if (!error_model.isa_asymptotic.empty())
            out << "_reflns.jfjoch_diffrn_ISa_asymptotic  " << CifStr(error_model.isa_asymptotic)
                << "   # strong-reflection asymptote (Diederichs); rotation path only\n";
        if (!error_model.a.empty())
            out << "_reflns.jfjoch_error_model_a          " << CifStr(error_model.a)
                << "   # sigma^2 = a*(sigma0^2 + b*I^2), XDS convention\n";
        if (!error_model.b.empty())
            out << "_reflns.jfjoch_error_model_b          " << CifStr(error_model.b) << "\n";
        // Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's
        // "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded.
        if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0)
            out << "_reflns.B_iso_Wilson_estimate     " << Fmt(statistics.wilson_b, 2) << "\n";
        // Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above).
        if (twinning.l_test_pairs > 0) {
            out << "_reflns.jfjoch_L_test_mean_abs_L      " << Fmt(twinning.mean_abs_l, 3)
                << "   # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n";
            out << "_reflns.jfjoch_L_test_mean_L_squared  " << Fmt(twinning.mean_l_squared, 3)
                << "   # <L^2> (untwinned 0.333, perfect twin 0.200)\n";
        }
        if (twinning.moment_reflections > 0)
            out << "_reflns.jfjoch_second_moment_I        " << Fmt(twinning.second_moment, 3)
                << "   # <I^2>/<I>^2 (untwinned 2.00, perfect twin 1.50)\n";
        // Radiation-damage monitor (rotation): the relative Debye-Waller B change from the first to the last
        // frame (A^2). A large magnitude flags a dose-dependent resolution-scale change = radiation damage;
        // positive is the typical direction (high-resolution intensity fades with dose). No standard mmCIF item.
        if (std::isfinite(statistics.radiation_damage_delta_b))
            out << "_reflns.jfjoch_radiation_damage_relative_B  " << Fmt(statistics.radiation_damage_delta_b, 2)
                << "   # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
        out << "#\n";

        // Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
        // per rotation-range batch, measured before any correction. rotation_start_deg = id * batch_deg.
        if (!statistics.radiation_damage_b_batch.empty()) {
            out << "loop_\n";
            out << "_jfjoch_radiation_damage_batch.id\n";
            out << "_jfjoch_radiation_damage_batch.rotation_start_deg\n";
            out << "_jfjoch_radiation_damage_batch.relative_B\n";
            for (size_t i = 0; i < statistics.radiation_damage_b_batch.size(); ++i)
                out << "  " << (i + 1) << "  "
                    << Fmt(static_cast<double>(i) * statistics.radiation_damage_batch_deg, 1) << "  "
                    << Fmt(statistics.radiation_damage_b_batch[i], 2) << "\n";
            out << "#\n";
        }

        out << "loop_\n";
        out << "_reflns_shell.d_res_high\n";
        out << "_reflns_shell.d_res_low\n";
        out << "_reflns_shell.number_measured_obs\n";
        out << "_reflns_shell.number_unique_obs\n";
        out << "_reflns_shell.pdbx_redundancy\n";
        out << "_reflns_shell.percent_possible_obs\n";
        out << "_reflns_shell.meanI_over_sigI_obs\n";
        out << "_reflns_shell.pdbx_Rrim_I_all\n";
        out << "_reflns_shell.pdbx_CC_half\n";
        if (has_anom)
            out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n";
        for (const auto &s : statistics.shells) {
            if (s.unique_reflections == 0)
                continue;
            out << 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);
            if (has_anom)
                out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3);
            out << "\n";
        }
        out << "#\n";
    }

    // ---------- _refln loop ----------
    out << "loop_\n";
    out << "_refln.index_h\n";
    out << "_refln.index_k\n";
    out << "_refln.index_l\n";
    out << "_refln.intensity_meas\n";
    out << "_refln.intensity_sigma\n";
    out << "_refln.pdbx_I_plus\n";
    out << "_refln.pdbx_I_plus_sigma\n";
    out << "_refln.pdbx_I_minus\n";
    out << "_refln.pdbx_I_minus_sigma\n";
    out << "_refln.F_meas_au\n";
    out << "_refln.F_meas_sigma_au\n";
    out << "_refln.pdbx_F_plus\n";
    out << "_refln.pdbx_F_plus_sigma\n";
    out << "_refln.pdbx_F_minus\n";
    out << "_refln.pdbx_F_minus_sigma\n";
    out << "_refln.status_free\n";
    out << "_refln.status\n";

    // One row per unique reflection, twelve formatted floats each - tens of megabytes on a crowded
    // crystal, and the largest single-threaded stretch left in a run. Nothing about a row depends on
    // any other, so each worker formats its own block into its own string and the blocks go to the
    // file in order. The columns are written exactly as the stream wrote them: the same "%.4f" (or
    // "?" where the value is not finite), right-aligned in the same width.
    {
        const size_t nrow = reflections.size();
        const size_t nw = std::max<size_t>(nthreads, 1);
        const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
        std::vector<std::string> 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) * 208);
                const auto column = [&s](const std::string &v, size_t width) {
                    if (v.size() < width) s.append(width - v.size(), ' ');
                    s.append(v);
                    s.push_back(' ');
                };
                for (size_t i = lo; i < hi; ++i) {
                    const auto &r = reflections[i];
                    column(std::to_string(r.h), 5);
                    column(std::to_string(r.k), 5);
                    column(std::to_string(r.l), 5);
                    column(Fmt(r.I, 4), 14);
                    column(Fmt(r.sigma, 4), 14);
                    column(Fmt(r.I_plus, 4), 14);
                    column(Fmt(r.sigma_plus, 4), 14);
                    column(Fmt(r.I_minus, 4), 14);
                    column(Fmt(r.sigma_minus, 4), 14);
                    column(Fmt(r.F, 4), 14);
                    column(Fmt(r.sigmaF, 4), 14);
                    column(Fmt(r.F_plus, 4), 14);
                    column(Fmt(r.sigmaF_plus, 4), 14);
                    column(Fmt(r.F_minus, 4), 14);
                    column(Fmt(r.sigmaF_minus, 4), 14);
                    s.push_back(r.rfree_flag ? '1' : '0');
                    s.append(" o\n");   // 'o' = observed
                }
            }
        });
        for (const std::string &s : block)
            out.write(s.data(), static_cast<std::streamsize>(s.size()));
    }

    out << "#\n";
    out << "# End of reflections\n";
    out.close();
}

void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
                         const UnitCell &unitCell,
                         const DiffractionExperiment &experiment,
                         const std::string &filename) {
    gemmi::Mtz mtz;

    // Optional but recommended metadata
    mtz.spacegroup = gemmi::find_spacegroup_by_number(
        experiment.GetSpaceGroupNumber().value_or(1));
    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");

    // Add dataset
    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<MergedOutRow> 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);
    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);
    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<int>(out_rows.size());
    mtz.data.reserve(out_rows.size() * (has_anom ? 16 : 8));
    for (const auto& row : out_rows) {
        mtz.data.push_back(static_cast<float>(row.h));
        mtz.data.push_back(static_cast<float>(row.k));
        mtz.data.push_back(static_cast<float>(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<float>(row.rfree));
    }
    mtz.write_to_file(filename);
}

void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
                              const DiffractionExperiment &experiment,
                              const std::string &filename,
                              size_t nthreads) {
    bool has_anom = true;
    const std::vector<MergedOutRow> 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);
    // Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker
    // blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in
    // the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write.
    const auto column = [](std::string &s, const std::string &v, size_t width) {
        if (v.size() < width) s.append(width - v.size(), ' ');
        s.append(v);
    };
    const auto num2 = [](double v) {
        char b[64];
        const int n = std::snprintf(b, sizeof b, "%.2f", v);
        return std::string(b, static_cast<size_t>(std::clamp(n, 0, static_cast<int>(sizeof b) - 1)));
    };
    const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) {
        column(s, std::to_string(h), 4);
        column(s, std::to_string(k), 4);
        column(s, std::to_string(l), 4);
        column(s, num2(scale * I), 8);
        column(s, num2(scale * sigma), 8);
        s.push_back('\n');
    };
    {
        const size_t nrow = rows.size();
        const size_t nw = std::max<size_t>(nthreads, 1);
        const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
        std::vector<std::string> 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) * 2 * 29);
                for (size_t i = lo; i < hi; ++i) {
                    const auto &r = rows[i];
                    const bool plus = usable(r.Ip, r.sIp);
                    const bool minus = usable(r.Im, r.sIm);
                    if (plus)  emit(s, r.h, r.k, r.l, r.Ip, r.sIp);
                    if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm);
                    if (!plus && !minus && usable(r.Imean, r.sImean))
                        emit(s, r.h, r.k, r.l, r.Imean, r.sImean);
                }
            }
        });
        for (const std::string &s : block)
            out.write(s.data(), static_cast<std::streamsize>(s.size()));
    }
    std::string tail;
    emit(tail, 0, 0, 0, 0.0f, 0.0f);   // HKLF-4 end-of-data marker
    out.write(tail.data(), static_cast<std::streamsize>(tail.size()));
    out.close();
}

void WriteReflections(const std::vector<MergedReflection> &reflections,
                      const UnitCell &unitCell,
                      const DiffractionExperiment &experiment,
                      const MergeStatistics &statistics,
                      const ErrorModelReport &error_model,
                      const TwinningAnalysisResult &twinning,
                      const std::string &filename,
                      size_t nthreads) {
    // 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, error_model, twinning,
                          filename + ".cif", nthreads);
    WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
}