// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <algorithm>
#include <fstream>
#include <sstream>

#include <spdlog/fmt/fmt.h>

#include "../common/JFJochMath.h"   // PI - M_PI is not portable to MSVC
#include "../common/GitInfo.h"
#include "../common/time_utc.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"

#include "ResultReport.h"

namespace {
    // The version of this file format. Bumped when a key is renamed or removed, a table column moves,
    // or a reason code changes meaning - a consumer can gate on it.
    constexpr int REPORT_VERSION = 5;

    const char *BANNER = " ******************************************************************************";

    void Section(std::ostream &os, const std::string &title) {
        os << "\n" << BANNER << "\n  " << title << "\n" << BANNER << "\n\n";
    }

    // Every number a consumer might want is written as one of these, so it is one grep away.
    template <class T> void Key(std::ostream &os, const char *key, const T &value) {
        os << key << "= " << value << "\n";
    }

    std::string CellString(const UnitCell &c) {
        return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c,
                           c.alpha, c.beta, c.gamma);
    }

    // Which reciprocal axis a principal anisotropy direction lies along. The tensor is fitted on
    // s = frac.mat * (h,k,l), so in that Cartesian frame a*, b*, c* are the rows of frac.mat and the
    // eigenvector is named by whichever it makes the smallest angle with. The label is exact in every
    // Laue class the tensor has a free direction in except triclinic; the cosine says how well it fits.
    std::string ReciprocalAxisLabel(const UnitCell &cell, const double v[3]) {
        static const char *NAME[3] = {"a*", "b*", "c*"};
        const gemmi::UnitCell gc = cell;
        int best = 0;
        double best_cos = -1.0;
        for (int i = 0; i < 3; ++i) {
            const gemmi::Vec3 axis = gc.frac.mat.left_multiply(
                    gemmi::Vec3(i == 0 ? 1.0 : 0.0, i == 1 ? 1.0 : 0.0, i == 2 ? 1.0 : 0.0));
            const double len = axis.length();
            if (!(len > 0.0))
                continue;
            const double c = std::fabs((axis.x * v[0] + axis.y * v[1] + axis.z * v[2]) / len);
            if (c > best_cos) {
                best_cos = c;
                best = i;
            }
        }
        return fmt::format("{} (cos {:.2f})", NAME[best], best_cos);
    }
}

std::string RenderResultReport(const std::string &output_prefix,
                               const std::string &input_file,
                               const DiffractionExperiment &experiment,
                               const ProcessResult &result,
                               const RunProvenance &provenance) {
    std::ostringstream os;
    const bool rotation = experiment.IsRotationIndexing();
    const bool merged = result.has_merge_statistics;
    std::vector<std::string> warnings = result.warnings;

    os << BANNER << "\n"
       << "  RUGNUX PROCESSING REPORT\n"
       << BANNER << "\n\n"
       << " What this run determined, written next to its other output. The `KEY= value` lines and\n"
       << " the tables below are a stable interface - a script greps them, and REPORT_VERSION says\n"
       << " when that interface last changed. Rates and per-image progress are not here; they are\n"
       << " on stdout.\n\n";

    Key(os, "REPORT_VERSION", REPORT_VERSION);
    Key(os, "RUGNUX_VERSION", jfjoch_version());
    if (!jfjoch_git_sha1().empty())
        Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date());
    Key(os, "DATE", time_UTC(std::chrono::system_clock::now()));
    Key(os, "INPUT_FILE", input_file);
    Key(os, "OUTPUT_PREFIX", output_prefix);
    // How the result was produced, what it cost and what it ran on, so the report stands on its own
    // once the shell history it came from is gone. Absent rather than zero where the caller does not
    // know them - the library and the viewer have no command line and no invocation to time.
    if (!provenance.command_line.empty())
        Key(os, "COMMAND_LINE", provenance.command_line);
    if (provenance.wall_time_s > 0.0)
        Key(os, "WALL_TIME", fmt::format("{:.2f}", provenance.wall_time_s));
    if (provenance.gpu_count >= 0) {
        Key(os, "GPU_COUNT", provenance.gpu_count);
        // GPU_COUNT= 0 with no GPU= line is the CPU-only case, and saying so is the point: whether
        // the GPUs were there is the first question about how long the run took.
        if (!provenance.gpu_description.empty())
            Key(os, "GPU", provenance.gpu_description);
    }

    // ---------------------------------------------------------------- 1. DATA SET
    Section(os, "1. DATA SET");
    Key(os, "EXPERIMENT_TYPE", rotation ? "ROTATION" : "STILLS");
    Key(os, "IMAGES_PROCESSED", result.images_processed);
    Key(os, "WAVELENGTH", fmt::format("{:.5f}", experiment.GetWavelength_A()));
    if (const auto gonio = experiment.GetGoniometer()) {
        Key(os, "OSCILLATION_RANGE", fmt::format("{:.4f}", gonio->GetIncrement_deg()));
        Key(os, "STARTING_ANGLE", fmt::format("{:.3f}", gonio->GetStart_deg()));
        const auto ax = gonio->GetAxis();
        Key(os, "ROTATION_AXIS", fmt::format("{:.6f} {:.6f} {:.6f}", ax.x, ax.y, ax.z));
    }
    Key(os, "DETECTOR_DISTANCE", fmt::format("{:.3f}", result.used_distance_mm));
    Key(os, "BEAM_CENTRE", fmt::format("{:.2f} {:.2f}", result.used_beam_x_pxl, result.used_beam_y_pxl));
    // The tilt this run integrated at, and the two points it separates. BEAM_CENTRE is the PONI - the
    // foot of the perpendicular from the sample - so on a tilted detector it is NOT where the direct
    // beam lands, and the two were previously indistinguishable in this report because only one of
    // them was printed. Degrees here; the JSON below carries radians, as the API spells it.
    {
        Key(os, "DETECTOR_TILT", fmt::format("{:.4f} {:.4f} {:.4f}", result.used_detector_tilt_deg[0],
                                             result.used_detector_tilt_deg[1],
                                             result.used_detector_tilt_deg[2]));
        Key(os, "DIRECT_BEAM", fmt::format("{:.2f} {:.2f}", result.used_direct_beam_x_pxl,
                                           result.used_direct_beam_y_pxl));
        if (result.refined_detector_tilt_deg)
            Key(os, "REFINED_DETECTOR_TILT",
                fmt::format("{:.4f} {:.4f}", (*result.refined_detector_tilt_deg)[0],
                            (*result.refined_detector_tilt_deg)[1]));
    }
    os << "\n"
       << " The distance and beam centre above are the ones this result was integrated at, which on\n"
       << " a rotation run is the post-refined geometry rather than the values in the input file.\n"
       << " DETECTOR_TILT is rot1/rot2/rot3 in degrees, as the run INTEGRATED at them. Rotation\n"
       << " indexing does refine rot1/rot2, but the result is never written back onto the geometry,\n"
       << " so what it measured is reported separately as REFINED_DETECTOR_TILT.\n"
       << " BEAM_CENTRE is the PONI and DIRECT_BEAM is where the beam actually lands; they differ by\n"
       << " distance*tan(tilt)/pixel and are identical only when the tilt is zero. Quote whichever the\n"
       << " program you are feeding expects, and check which one it means.\n";
    if (result.refined_detector_tilt_deg)
        os << " REFINED_DETECTOR_TILT is rot1/rot2 in degrees as rotation indexing MEASURED them, and is\n"
           << " NOT what this run integrated at - nothing writes a refined tilt back onto the geometry.\n"
           << " One crystal does not measure a detector tilt. The tilt is aliased with the beam centre,\n"
           << " so what a single sweep determines is the DIRECT_BEAM above, not the split between the\n"
           << " two, and this fit stays close to the tilt it was started from - it moves a fraction of\n"
           << " the way to the true value, and the fraction varies from crystal to crystal. Comparing\n"
           << " one run against a powder calibration will therefore disagree with it by far more than\n"
           << " either number's uncertainty. What is worth comparing is the MEDIAN of this value over\n"
           << " several crystals collected on the same detector, which does track the calibration well\n"
           << " enough to show up a placeholder or a stale tilt in the file - it does not replace the\n"
           << " calibration. Do not feed a single run's value back into the instrument.\n"
           << " rot3 is omitted because a rotation about the beam is an exact null of this experiment\n"
           << " and the fit cannot move it.\n";

    // The same geometry once more, as the object jfjoch_broker takes it in: the four required
    // properties of dataset_settings in broker/jfjoch_api.yaml, spelled the way the API spells them.
    // A run that refined the geometry is usually the best measurement of it anyone has, and without
    // this the only way back into the instrument is to read two numbers off this report by eye and
    // retype them. One line, valid JSON, so a script can lift it with a grep and POST it.
    os << "\n";
    {
        // The rotations belong here whenever they are not zero: dataset_settings carries
        // poni_rot1/2/3_rad, and a block that omits them describes a FLAT detector - a different
        // geometry from the one this run used, silently, on every tilted setup. Omitted when zero
        // because the API's own default is 0.0, so the shorter block means the same thing.
        constexpr double RAD = PI / 180.0;
        const auto &t = result.used_detector_tilt_deg;
        std::string rot;
        if (t[0] != 0.0 || t[1] != 0.0 || t[2] != 0.0)
            rot = fmt::format(R"(, "poni_rot1_rad": {:.6f}, "poni_rot2_rad": {:.6f}, )"
                              R"("poni_rot3_rad": {:.6f})",
                              t[0] * RAD, t[1] * RAD, t[2] * RAD);
        Key(os, "JFJOCH_DATASET_SETTINGS",
            fmt::format(R"({{"beam_x_pxl": {:.2f}, "beam_y_pxl": {:.2f}, "detector_distance_mm": {:.3f}, )"
                        R"("incident_energy_keV": {:.4f}{}}})",
                        result.used_beam_x_pxl, result.used_beam_y_pxl, result.used_distance_mm,
                        experiment.GetDatasetSettings().GetPhotonEnergy_keV(), rot));
    }
    os << "\n"
       << " The geometry above as jfjoch_broker's dataset_settings, to carry a refined beam centre and\n"
       << " distance back to the instrument for the next collection. The PONI rotations ride with it\n"
       << " when they are non-zero, because without them the block describes a flat detector.\n";

    if (result.spot_resolution_estimate_A.has_value()) {
        os << "\n";
        Key(os, "SPOT_RESOLUTION_ESTIMATE", fmt::format("{:.2f}", *result.spot_resolution_estimate_A));
        os << "\n"
           << " How far the merged data are expected to reach, read off the found spots alone - no\n"
           << " lattice, no integration, no merge. It is a prediction, good to about 0.2 A on the\n"
           << " rotation data it was calibrated on, and it is not what the run achieved: compare it\n"
           << " with INCLUDE_RESOLUTION_RANGE in section 5. It is not limited to what this detector\n"
           << " records: where it reads finer than the high-resolution end of that range, the crystal\n"
           << " diffracts past the corner and the run is detector-limited.\n";
    }

    if (result.pass_count > 1) {
        os << "\n";
        Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count));
        Key(os, "PASS_DECISION", result.pass_decision);
        os << "\n"
           << " A rotation run integrates twice: once at the geometry in the input file, then again at\n"
           << " the post-refined geometry. Every number in this report describes the pass named above,\n"
           << " whose files are " << output_prefix << ".*; the header-geometry pass is not written.\n";
    }

    // ---------------------------------------------------------------- 2. INDEXING
    Section(os, "2. INDEXING");
    if (result.indexing_rate.has_value())
        Key(os, "INDEXING_RATE", fmt::format("{:.4f}", result.indexing_rate.value()));
    Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE"));
    if (result.consensus_cell.has_value())
        Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell));
    if (result.space_group.has_value()) {
        Key(os, "SPACE_GROUP_NUMBER", result.space_group->number);
        // The number alone does not name the setting, and the search can now return one that is not
        // gemmi's reference setting - P 1 1 2_1 and P 1 2_1 1 are both number 4, on different axes.
        Key(os, "SPACE_GROUP_NAME", result.space_group->xhm());
    }
    if (result.indexing_rate.value_or(0.0f) <= 0.0f)
        warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset");

    // ---------------------------------------------- 3. GEOMETRY POST-REFINEMENT
    if (result.post_refine.has_value()) {
        const auto &pr = *result.post_refine;
        Section(os, "3. GEOMETRY POST-REFINEMENT");
        os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n"
           << " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n"
           << " step is committed only if it improves a held-out residual.\n\n";
        Key(os, "POSTREFINE_EVENTS_USED", pr.events_used);
        Key(os, "POSTREFINE_OBS_USED", pr.obs_used);
        Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE");
        Key(os, "POSTREFINE_DETECTOR_COMMITTED", pr.detector_refined ? "TRUE" : "FALSE");
        Key(os, "POSTREFINE_DISTANCE", fmt::format("{:.3f} -> {:.3f}", pr.distance_before_mm,
                                                   pr.distance_after_mm));
        Key(os, "POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}",
                                                      pr.beam_x_before_px, pr.beam_y_before_px,
                                                      pr.beam_x_after_px, pr.beam_y_after_px));
        Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale));
        Key(os, "GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect ? "TRUE" : "FALSE");
        os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n"
           << " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It drives\n"
           << " the second integration pass only when SUSPECT is TRUE - both cross-validated and outside\n"
           << " the tolerance - since a stage that is in fact well calibrated must be left alone. A\n"
           << " manual --rotation-scale replaces it and is applied to both passes.\n";
        if (pr.rotation_scale_suspect)
            warnings.emplace_back(fmt::format(
                    "The goniometer turned by a factor {:.5f} of the angles stored in the file - the "
                    "stage rotation looks mis-calibrated by {:+.2f}%. The correction was applied to this "
                    "run, but the fault is in the hardware and should be fixed there",
                    pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0)));
    }

    // ---------------------------------------------- 4. SPACE GROUP DETERMINATION
    Section(os, "4. SPACE GROUP DETERMINATION");

    // What the data could NOT decide, as keys rather than only as prose. SPACE_GROUP_NAME above is a
    // scalar and reads like a determination; where several groups predict the same absences it is one
    // of them, chosen by convention. A script that greps only the name records a coin-flip as an
    // answer - measured, refining against the deposited model in the wrong enantiomorph gives
    // R = 0.549 - so the ambiguity has to survive the same grep.
    {
        std::string alts;
        if (result.space_group_search.has_value() && result.space_group.has_value())
            for (const auto &alt : result.space_group_search->alternatives)
                if (alt.number != result.space_group->number)
                    alts += (alts.empty() ? "" : " | ") + alt.xhm();
        Key(os, "SPACE_GROUP_ALTERNATIVES", alts.empty() ? "NONE" : alts);

        // The hand is a separate question from the alternatives list, and it is decidable from the
        // 22 groups' numbers alone, so it is answered even where the search did not run. Merged
        // intensities never decide it: an enantiomorphic pair has the same absences and the same
        // Laue class, and only a model, a substructure or an anomalous signal names the hand.
        const char *enantiomorph = "NOT_APPLICABLE";
        if (result.space_group.has_value() && result.space_group->is_enantiomorphic()) {
            if (result.model_validation.has_value() && result.model_validation->ok)
                enantiomorph = "DETERMINED_FROM_MODEL";
            else if (!result.space_group_search.has_value())
                enantiomorph = "GIVEN";       // -S, or a reference MTZ: the user's assertion
            else
                enantiomorph = "UNDETERMINED";
        }
        Key(os, "SPACE_GROUP_ENANTIOMORPH", enantiomorph);

        // A higher point group whose operators the intensities confirmed and whose promotion was
        // refused. It exists only as prose in the search text below, so nothing can act on it - and
        // a refusal is exactly the case where a user might want to try the higher group as well.
        const bool refused = result.space_group_search.has_value()
                             && !result.space_group_search->refused_point_group_hm.empty();
        Key(os, "SPACE_GROUP_REFUSED_POINT_GROUP",
            refused ? result.space_group_search->refused_point_group_hm : std::string("NONE"));
        if (refused)
            Key(os, "SPACE_GROUP_REFUSED_REASON", result.space_group_search->refused_reason);
        os << "\n SPACE_GROUP_ALTERNATIVES names every group these data cannot separate from the one\n"
           << " adopted (enantiomorphic partners, origin-ambiguous pairs, groups a gap in the data\n"
           << " leaves untested); NONE means the absences single the answer out. SPACE_GROUP_ENANTIOMORPH\n"
           << " is UNDETERMINED whenever the group is one of the 22 that come in enantiomorphic pairs and\n"
           << " nothing outside the merged intensities named the hand - the intensities cannot, so this\n"
           << " is the normal outcome, not a failure. SPACE_GROUP_REFUSED_POINT_GROUP is a higher point\n"
           << " group the operator correlations supported and the consistency tests would not take.\n";
    }

    if (result.space_group_search.has_value()) {
        Key(os, "SPACE_GROUP_SEARCH", "DE_NOVO");
        os << "\n" << SearchSpaceGroupResultToText(*result.space_group_search) << "\n";
        // A centering the data could not test must not read like one they confirmed. The group may
        // still be right - the lattice metric says so - but nothing in these intensities backs it,
        // and that belongs beside the warnings rather than in a table column alone.
        const auto &search = *result.space_group_search;
        if (search.best_space_group.has_value())
            for (const auto &c : search.candidates)
                if (c.space_group.number == search.best_space_group->number && c.centering_untested)
                    warnings.emplace_back(fmt::format(
                            "The {} centering of {} was NOT confirmed from these data: the crystal was "
                            "indexed and integrated on the primitive sub-cell, so the reflections a "
                            "{}-centred lattice extinguishes are not in this merge at all. It comes "
                            "from the lattice metric. The point group is confirmed from the "
                            "intensities; the centering is not",
                            search.best_space_group->centring_type(),
                            search.best_space_group->short_name(),
                            search.best_space_group->centring_type()));
    } else if (result.space_group.has_value()) {
        Key(os, "SPACE_GROUP_SEARCH", "FIXED");
        os << "\n The space group was given, not determined here.\n";
    } else {
        Key(os, "SPACE_GROUP_SEARCH", "NONE");
        os << "\n No space group was determined.\n";
    }

    // ---------------------------------------------------- 5. SCALING AND MERGING
    Section(os, "5. SCALING AND MERGING");
    if (!merged) {
        Key(os, "MERGE", "NOT_PERFORMED");
        os << "\n No scaling or merging was performed on this run, so there are no merging statistics, no\n"
           << " error model, and no sweep-quality diagnosis below. The integrated reflections are in\n"
           << " " << output_prefix << "_process.h5.\n";
    } else {
        const auto &o = result.merge_statistics.overall;
        Key(os, "MERGE", "PERFORMED");
        Key(os, "INCLUDE_RESOLUTION_RANGE", fmt::format("{:.3f} {:.3f}", o.d_max, o.d_min));
        if (result.resolution_fit_A)
            Key(os, "FITTED_RESOLUTION", fmt::format("{:.2f}", *result.resolution_fit_A));
        Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE");
        Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections);
        Key(os, "TOTAL_OBSERVATIONS", o.total_observations);
        Key(os, "COMPLETENESS", o.possible_unique_reflections > 0
                ? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections)
                : std::string("nan"));
        Key(os, "MULTIPLICITY", o.unique_reflections > 0
                ? fmt::format("{:.2f}", static_cast<double>(o.total_observations) / o.unique_reflections)
                : std::string("nan"));
        Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma));
        Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas));
        Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half));
        Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous));
        Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b));
        // The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP.
        Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a));
        Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b));
        Key(os, "ISA", fmt::format("{:.2f}", result.error_model_isa));
        if (result.error_model_isa_asymptotic > 0.0)
            Key(os, "ISA_ASYMPTOTIC", fmt::format("{:.2f}", result.error_model_isa_asymptotic));
        Key(os, "REFERENCE_DATA_USED", result.has_reference ? "TRUE" : "FALSE");
        // The shell table straight off the statistics rather than result.merge_statistics_text: that
        // string also carries the twinning analysis and the advisories, which have sections of their own.
        os << "\n";
        if (result.resolution_fit_A)
            os << fmt::format(
                    " INCLUDE_RESOLUTION_RANGE is the range the reflections were WRITTEN to;"
                    " FITTED_RESOLUTION\n is where the CC1/2 fall-off crosses {:.2f}, and is the number to"
                    " quote. The data are kept one\n shell past it on purpose: a shell that is included can"
                    " still be downweighted or dropped by\n refinement, while one that was truncated cannot"
                    " be put back.\n\n",
                    experiment.GetScalingSettings().GetResolutionCCTarget());
        os << " ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n"
           << " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n"
           << " is the strong-reflection tier only.\n\n"
           << result.merge_statistics;
    }

    // --------------------------------------------------------------- 6. TWINNING
    if (merged && result.twinning.l_test_pairs > 0) {
        Section(os, "6. TWINNING");
        Key(os, "TWINNING_SUSPECTED", result.twinning.twinning_suspected ? "TRUE" : "FALSE");
        Key(os, "L_TEST_MEAN_ABS_L", fmt::format("{:.4f}", result.twinning.mean_abs_l));
        Key(os, "L_TEST_MEAN_L_SQUARED", fmt::format("{:.4f}", result.twinning.mean_l_squared));
        Key(os, "SECOND_MOMENT_I", fmt::format("{:.4f}", result.twinning.second_moment));
        // The same two statistics measured before the group was chosen. Printed next to the ones above
        // rather than instead of them, because they answer different questions and can disagree: these
        // are the only ones a promotion cannot have contaminated.
        if (result.pre_promotion_twinning) {
            Key(os, "L_TEST_MEAN_ABS_L_BEFORE_SEARCH",
                fmt::format("{:.4f}", result.pre_promotion_twinning->mean_abs_l));
            Key(os, "SECOND_MOMENT_I_BEFORE_SEARCH",
                fmt::format("{:.4f}", result.pre_promotion_twinning->second_moment));
        }
        Key(os, "ESTIMATED_TWIN_FRACTION", fmt::format("{:.3f}", result.twinning.estimated_twin_fraction));
        os << "\n" << TwinningAnalysisToText(result.twinning) << "\n";
        if (result.pre_promotion_twinning)
            os << " The _BEFORE_SEARCH pair was measured on the merge the space-group search was given,\n"
               << " before any point group was adopted. That ordering is the whole point: the statistics\n"
               << " above are computed in the Laue class this run ADOPTED, so if the search promoted the\n"
               << " point group they can only report that no twin law exists inside the class it chose -\n"
               << " and a twin is precisely what would have caused that promotion. Where the two pairs\n"
               << " disagree, believe the _BEFORE_SEARCH one about whether the crystal is twinned.\n";
        if (result.twinning.twinning_suspected)
            warnings.emplace_back(fmt::format(
                    "Twinning is indicated (<|L|> = {:.3f}, <I^2>/<I>^2 = {:.3f}, estimated twin "
                    "fraction {:.2f}) - refine against the merged data with care",
                    result.twinning.mean_abs_l, result.twinning.second_moment,
                    result.twinning.estimated_twin_fraction));
    }

    // ------------------------------------------------------- 7. RADIATION DAMAGE
    if (!result.radiation_damage_text.empty()) {
        Section(os, "7. RADIATION DAMAGE");
        // A number, or a word saying why there is none: NOT_A_TREND where the per-batch curve was measured
        // but no straight line describes it (damage is progressive, so that curve is not dose), NOT_MEASURED
        // where the monitor could not run at all.
        const double db = result.merge_statistics.radiation_damage_delta_b;
        Key(os, "RADIATION_DAMAGE_RELATIVE_B",
            std::isfinite(db) ? fmt::format("{:.2f}", db)
            : result.merge_statistics.radiation_damage_b_batch.empty() ? std::string("NOT_MEASURED")
                                                                      : std::string("NOT_A_TREND"));
        os << "\n" << result.radiation_damage_text << "\n";
    }

    // ------------------------------------------------------ 8. SWEEP QUALITY
    const auto &sq = result.merge_statistics.sweep_quality;
    Section(os, "8. SWEEP QUALITY");
    os << " Stretches of the sweep over which the crystal delivered much less than the rest of the run.\n"
       << " REASON comes from a closed vocabulary, listed below so a consumer can tell an unknown code\n"
       << " from a missing one. SEVERITY is the fraction of the run's typical diffracting power missing\n"
       << " over the range (0 = as good as the run, 1 = nothing at all); SCALE and CC are the range's\n"
       << " mean per-image scale and CC-to-merge relative to the run median; INDEXED is the fraction of\n"
       << " the range's frames that were scaled at all. Nothing is excluded on the strength of this.\n\n";
    Key(os, "SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED");
    Key(os, "SWEEP_QUALITY_COUNT", sq.ranges.size());
    {
        std::string codes;
        for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
            codes += (codes.empty() ? "" : " ")
                     + std::string(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
        Key(os, "SWEEP_QUALITY_REASONS", codes);
    }
    if (sq.measured) {
        Key(os, "SWEEP_ROTATION", fmt::format("{:.1f}", sq.sweep_deg));
        Key(os, "FLUX_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.flux_peak_to_trough));
        Key(os, "SCALE_MODULATION_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.modulation_peak_to_trough));
    }
    os << "\n"
       << "  FIRST_IMAGE   LAST_IMAGE   N_IMAGES  ROTATION  REASON                SEVERITY   SCALE      CC   INDEXED\n"
       << "  -----------  -----------  ---------  --------  --------------------  --------  ------  ------  --------\n";
    for (const auto &r : sq.ranges) {
        os << fmt::format("  {:11d}  {:11d}  {:9d}  {:8.1f}  {:<20}  {:8.2f}  {:6.2f}  {:6.2f}  {:8.2f}\n",
                          r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg,
                          SweepQualityReasonCode(r.reason), r.severity, r.mean_relative_scale,
                          r.mean_relative_cc, r.indexed_fraction);
        warnings.push_back(fmt::format(
                "Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)",
                r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg,
                r.mean_relative_scale, r.mean_relative_cc, 100.0 * r.indexed_fraction));
    }
    os << "  -----------  -----------  ---------  --------  --------------------  --------  ------  ------  --------\n";

    // ---------------------------------------------------------- 9. DIFFRACTION ANISOTROPY
    const auto &an = result.merge_statistics.anisotropy;
    if (merged && an.n_reflections > 0) {
        Section(os, "9. DIFFRACTION ANISOTROPY");
        os << " How much the fall-off depends on direction, and whether that is established above this\n"
           << " data set's own systematic error. Nothing here corrects an intensity or removes a\n"
           << " reflection: the merged data and the written files do not depend on direction at all.\n"
           << " ANISOTROPY_DELTA_B is the range of the principal components of the anisotropy tensor,\n"
           << " on the ordinary crystallographic B scale (the same scale as phenix.xtriage's B_cart and\n"
           << " ctruncate's anisotropic B), fitted on intensities with nothing dropped;\n"
           << " ANISOTROPY_SIGNIFICANCE gates ANISOTROPY_DELTA_B_LINEAR, the part of it that follows\n"
           << " exp(-1/2 s^T B s), which is not the same number. A 1 in ANISOTROPY_D_MIN_CENSORED marks\n"
           << " a direction whose limit is the edge of the measured data rather than the crystal's own.\n\n";
        Key(os, "ANISOTROPY_VERDICT", AnisotropyVerdictCode(an.verdict));
        Key(os, "ANISOTROPY_FREE_DIRECTIONS", an.n_free_parameters);
        Key(os, "ANISOTROPY_DELTA_B", fmt::format("{:.2f}", an.delta_b));
        Key(os, "ANISOTROPY_DELTA_B_LINEAR", fmt::format("{:.2f}", an.delta_b_linear));
        Key(os, "ANISOTROPY_PRINCIPAL_B", fmt::format("{:.2f} {:.2f} {:.2f}",
                                                      an.eigenvalue[0] - an.eigenvalue[2],
                                                      an.eigenvalue[1] - an.eigenvalue[2], 0.0));
        Key(os, "ANISOTROPY_FOLD_WEAKENING", fmt::format("{:.1f}", an.fold_weakening));
        Key(os, "ANISOTROPY_D_MIN_PRINCIPAL", fmt::format("{:.2f} {:.2f} {:.2f}", an.d_min_axis[0],
                                                          an.d_min_axis[1], an.d_min_axis[2]));
        Key(os, "ANISOTROPY_D_MIN_CENSORED", fmt::format("{} {} {}", an.d_min_censored[0] ? 1 : 0,
                                                          an.d_min_censored[1] ? 1 : 0,
                                                          an.d_min_censored[2] ? 1 : 0));
        Key(os, "ANISOTROPY_D_MIN_SPREAD", fmt::format("{:.2f}", an.d_min_spread));
        // The finest of the three principal limits, on its own line so it can be grepped. The cut this
        // run applied is isotropic, so on an anisotropic crystal the outer shells are complete in count
        // and empty in signal along the weak directions; this says how far the crystal actually reaches
        // where it reaches furthest.
        Key(os, "ANISOTROPY_D_MIN_BEST",
            fmt::format("{:.2f}", *std::min_element(an.d_min_axis, an.d_min_axis + 3)));
        Key(os, "ANISOTROPY_SHAPE", AnisotropyShapeCode(an.shape));
        Key(os, "ANISOTROPY_SHAPE_INTERCEPT", fmt::format("{:.3f}", an.shape_intercept));
        Key(os, "ANISOTROPY_SHAPE_INTERCEPT_Z", fmt::format("{:.1f}", an.shape_intercept_z));
        Key(os, "ANISOTROPY_SHAPE_SLOPE", fmt::format("{:.2f}", an.shape_slope));
        Key(os, "ANISOTROPY_SHAPE_RESIDUAL", fmt::format("{:.1f}", an.shape_residual));
        Key(os, "ANISOTROPY_N_OBSERVATIONS", an.n_observations);
        Key(os, "ANISOTROPY_SIGMA_SYSTEMATIC", fmt::format("{:.3f}", an.sigma_systematic));
        Key(os, "ANISOTROPY_FORBIDDEN_Z", fmt::format("{:.1f}", an.forbidden_z));
        Key(os, "ANISOTROPY_FLOOR", fmt::format("{:.3f}", an.floor));
        Key(os, "ANISOTROPY_SIGNIFICANCE", fmt::format("{:.2f}", an.significance));
        Key(os, "ANISOTROPY_DETECTION_LIMIT", fmt::format("{:.2f}", an.detection_limit));
        os << "\n" << AnisotropyToText(an) << "\n";
        // Three different resolutions appear in this report and they answer three different questions.
        // Say so here rather than let a reader assume one of them is "the" resolution.
        if (result.resolution_fit_A)
            os << fmt::format(
                    " Three resolutions, three questions. INCLUDE_RESOLUTION_RANGE ({:.2f} A) is what was\n"
                    " WRITTEN, deliberately one shell past the fit. FITTED_RESOLUTION ({:.2f} A) is where\n"
                    " the isotropic CC1/2 fall-off crosses its target - the overall measurement.\n"
                    " ANISOTROPY_D_MIN_BEST ({:.2f} A) is how far the crystal reaches along its strongest\n"
                    " direction. It is read off the fitted tensor, not off the cut, so it can land either\n"
                    " side of the other two - and where ANISOTROPY_D_MIN_CENSORED flags that direction it\n"
                    " is the edge of the measured data rather than the crystal's limit. None of the three\n"
                    " is wrong; quoting one without saying which it is, is.\n\n",
                    result.merge_statistics.overall.d_min, *result.resolution_fit_A,
                    *std::min_element(an.d_min_axis, an.d_min_axis + 3));
        // CC1/2 that falls and then climbs again is not a fall-off, so the single number the cut is read
        // off does not describe these data. Reported as an observation, with no advice to re-cut: the
        // generous cut is deliberate and refinement can downweight what it does not want.
        {
            // A climb only counts while the shell still carries signal: past CC1/2 ~ 0.3 the data are
            // beyond any cut this run would take and the number is oscillating in noise, which says
            // nothing about the fall-off. 0.05 is well above that surviving noise.
            constexpr double CC_HALF_RISE = 0.05, CC_HALF_ALIVE = 0.30;
            const auto &sh = result.merge_statistics.shells;
            size_t fell = 0, rose_after = 0;
            for (size_t i = 1; i < sh.size(); ++i) {
                if (sh[i].cc_half < sh[i - 1].cc_half) fell = i;
                else if (fell > 0 && sh[i].cc_half > sh[i - 1].cc_half + CC_HALF_RISE
                         && sh[i].cc_half > CC_HALF_ALIVE) rose_after = i;
            }
            if (rose_after > 0)
                warnings.emplace_back(fmt::format(
                        "CC1/2 is not monotone with resolution (it climbs again at {:.2f}-{:.2f} A) - the "
                        "fall-off the resolution cut is read off does not describe these data",
                        sh[rose_after].d_max, sh[rose_after].d_min));
        }
        if (an.verdict == AnisotropyVerdict::Detected && an.d_min_spread > 0.5) {
            // Say WHICH direction each limit belongs to. Without it the warning states that the crystal
            // is anisotropic and leaves the reader no way to act on it; the eigenvectors are measured
            // here and reach the mmCIF, so the name costs nothing.
            const double *worst = std::max_element(an.d_min_axis, an.d_min_axis + 3);
            const double *best = std::min_element(an.d_min_axis, an.d_min_axis + 3);
            const auto along = [&](const double *it) {
                const int n = static_cast<int>(it - an.d_min_axis);
                return result.consensus_cell
                       ? ReciprocalAxisLabel(*result.consensus_cell, an.eigenvector[n])
                       : fmt::format("principal direction {}", n + 1);
            };
            warnings.emplace_back(fmt::format(
                    "Diffraction is anisotropic (deltaB {:.1f} A^2; the diffraction limit is {:.2f} A "
                    "along {} and {:.2f} A along {}) - refinement and map interpretation should allow "
                    "for it; no intensity has been corrected for it here",
                    an.delta_b, *worst, along(worst), *best, along(best)));
        }
    }

    // ------------------------------------------------------- 10. MODEL VALIDATION
    if (result.model_validation.has_value()) {
        const auto &mv = *result.model_validation;
        Section(os, "10. MODEL VALIDATION");
        Key(os, "MODEL_FILE", mv.model_path);
        if (!mv.ok) {
            Key(os, "MODEL_VALIDATION", "NOT_PERFORMED");
            Key(os, "MODEL_VALIDATION_REASON", mv.failure_reason);
            os << "\n A model was given but could not be used, so there are no R-factors and no maps.\n"
               << " Everything else in this report is unaffected: the merge does not depend on the model.\n";
        } else {
            Key(os, "MODEL_SPACE_GROUP_NUMBER", mv.model_space_group_number);
            Key(os, "R_WORK", fmt::format("{:.4f}", mv.r_work));
            Key(os, "R_FREE", fmt::format("{:.4f}", mv.r_free));
            Key(os, "R_WORK_REFLECTIONS", mv.n_work);
            Key(os, "R_FREE_REFLECTIONS", mv.n_free);
            Key(os, "BULK_SOLVENT_K_SOL", fmt::format("{:.3f}", mv.k_sol));
            Key(os, "BULK_SOLVENT_B_SOL", fmt::format("{:.1f}", mv.b_sol));
            Key(os, "SCALE_OVERALL", fmt::format("{:.4f}", mv.k_overall));
            Key(os, "MEAN_ATOM_DENSITY_SIGMA", fmt::format("{:.2f}", mv.mean_atom_density_sigma));
            // The anomalous scatterers the data themselves found, named by the model's atoms.
            if (!mv.anomalous_sites.empty()) {
                Key(os, "ANOMALOUS_BIJVOET_PAIRS", mv.anomalous_pairs);
                for (size_t i = 0; i < mv.anomalous_sites.size(); i++)
                    // Two digits so the ten keys are the same width and the values line up.
                    Key(os, fmt::format("ANOMALOUS_SITE_{:02}", i + 1).c_str(),
                        fmt::format("{:<18} {:6.2f} sigma", mv.anomalous_sites[i].label,
                                    mv.anomalous_sites[i].sigma));
            }
            // What was applied to the written reflections, so a reader can tell whether the file is
            // in the indexing it was merged in or in the model's. The enantiomorph is a label only.
            Key(os, "MODEL_ENANTIOMORPH_ADOPTED", mv.adopted_model_enantiomorph ? "TRUE" : "FALSE");
            Key(os, "MODEL_INDEXING_OPERATOR", mv.indexing_op.triplet());
            if (!mv.maps_prefix.empty())
                Key(os, "MAPS_PREFIX", mv.maps_prefix);
            os << "\n R-free here measures the merged intensities against an external structure, which is\n"
               << " what CC1/2 and R_meas cannot do - they only measure the data against themselves. The\n"
               << " model is not refined: it is scaled to the data with a flat bulk solvent and an overall\n"
               << " anisotropic B, so these R-factors are higher than a refined structure's and are a\n"
               << " data-quality reading, not a refinement result.\n";
            if (!mv.anomalous_sites.empty())
                os << "\n The anomalous sites are the highest peaks of the anomalous difference map -\n"
                   << " F(+)-F(-) on the model phase turned back by 90 degrees - read at the model's own\n"
                   << " atom centres, so each one is named rather than left as a coordinate. A dataset with\n"
                   << " no anomalous signal still lists ten sites: it is their height, a few sigma at most,\n"
                   << " that says so. The map itself is written as _anom.ccp4, where a scatterer the model\n"
                   << " does not contain would show up as a peak on nothing.\n";
            if (!(mv.indexing_op == gemmi::Op::identity()))
                os << "\n The written reflections were reindexed into the model's frame - the operator above\n"
                   << " says how - so the reflection files, the R-factors and the maps all describe one\n"
                   << " indexing.\n";
            if (mv.adopted_model_enantiomorph)
                os << "\n The written reflections carry the model's enantiomorph as their space group. That is a\n"
                   << " change of label and nothing else: the two groups of an enantiomorphic pair have the same\n"
                   << " rotation operations, so no reflection moved. Reindexing by the change-of-hand operator\n"
                   << " would have swapped I(+) with I(-) - flipping the anomalous differences, not correcting\n"
                   << " them - on the strength of a label the space-group search reports as undetermined.\n";
            if (mv.anomalous_hands_disagree)
                warnings.emplace_back(fmt::format(
                        "The anomalous density at the model's atoms is inverted ({} reads {:.1f} sigma, deeper "
                        "than the highest peak) - the data and the model are in opposite hands. Either the model "
                        "is the wrong enantiomorph for this crystal or the data were indexed in the wrong hand; "
                        "the reflections have not been reindexed, which would have hidden which of the two it is",
                        mv.anomalous_deepest_site, mv.anomalous_deepest_sigma));
        }
    }

    // --------------------------------------------------------------- 11. WARNINGS
    if (result.cancelled)
        warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report "
                                          "describes an incomplete run", result.images_processed));
    Section(os, "11. WARNINGS");
    os << " Everything that needs a person's attention, one line each, marked so a script can find\n"
       << " them with a single grep for \"WARNING:\".\n\n";
    Key(os, "WARNING_COUNT", warnings.size());
    os << "\n";
    for (const auto &w : warnings)
        os << "WARNING: " << w << "\n";
    if (warnings.empty())
        os << " (none)\n";

    os << "\n" << BANNER << "\n  END OF REPORT\n" << BANNER << "\n";
    return os.str();
}

void WriteResultReport(const std::string &output_prefix,
                       const std::string &input_file,
                       const DiffractionExperiment &experiment,
                       const ProcessResult &result,
                       Logger &logger,
                       const RunProvenance &provenance) {
    if (output_prefix.empty())
        return;   // "compute the statistics, persist nothing"

    const std::string filename = output_prefix + "_report.txt";
    // The report is unconditional, so it must never be the reason a run fails: a run that produced a
    // good .mtz must survive an unwritable path or a full disk. Report the failure and carry on.
    try {
        std::ofstream file(filename);
        file.exceptions(std::ios::failbit | std::ios::badbit);
        file << RenderResultReport(output_prefix, input_file, experiment, result, provenance);
    } catch (const std::exception &e) {
        logger.Warning("Could not write the results report {}: {}", filename, e.what());
    }
}