Jungfraujoch/rugnux/RasterReport.cpp

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

#include <cmath>
#include <fstream>
#include <sstream>

#include <nlohmann/json.hpp>
#include <spdlog/fmt/fmt.h>

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

#include "RasterReport.h"

namespace {
    // The version of this file format. Bumped when a key is renamed or removed or a table column
    // moves - a consumer can gate on it.
    constexpr int RASTER_REPORT_VERSION = 3;

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

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

    template <class T> void Key(std::ostream &os, const char *key, const T &value) {
        os << key << "= " << value << "\n";
    }

    // What the per-image scores add up to over the whole raster. Both renderers say the same thing,
    // so both read it from here.
    struct ScanSummary {
        int64_t scored = 0;         // images that came back with a protein score at all
        int64_t above = 0;          // ... of which above the threshold the blob search used
        int64_t ice_above = 0;      // ... whose ICE score is above that same threshold
        float mean_protein = 0;
        float mean_ice = 0;
    };

    ScanSummary Summarize(const ScanResult &scan, float threshold) {
        ScanSummary s;
        float sum_protein = 0, sum_ice = 0;
        for (const auto &e : scan.images) {
            if (!e.protein_score)
                continue;
            s.scored++;
            sum_protein += *e.protein_score;
            if (*e.protein_score > threshold)
                s.above++;
            if (e.ice_score) {
                sum_ice += *e.ice_score;
                if (*e.ice_score > threshold)
                    s.ice_above++;
            }
        }
        if (s.scored > 0) {
            s.mean_protein = sum_protein / static_cast<float>(s.scored);
            s.mean_ice = sum_ice / static_cast<float>(s.scored);
        }
        return s;
    }

    float Fraction(int64_t part, int64_t whole) {
        return whole > 0 ? static_cast<float>(part) / static_cast<float>(whole) : 0.0f;
    }
}

std::string RenderRasterReport(const std::string &input_file,
                               const GridScanSettings &grid,
                               const ScanResult &scan,
                               const GridScanResult &crystals,
                               const RasterSettings &settings,
                               const RunProvenance &provenance) {
    std::ostringstream os;
    const ScanSummary summary = Summarize(scan, settings.analysis.GetProteinScoreThreshold());

    os << BANNER << "\n"
       << "  RUGNUX RASTER REPORT\n"
       << BANNER << "\n\n"
       << " What this grid scan found: where the crystals are, how large they are and how well they\n"
       << " diffract. The `KEY= value` lines and the table below are a stable interface - a script\n"
       << " greps them, and RASTER_REPORT_VERSION says when that interface last changed. The same\n"
       << " content, typed, is in the _raster.json beside this file.\n\n";

    Key(os, "RASTER_REPORT_VERSION", RASTER_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);
    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));

    // ------------------------------------------------------------------------------ 1. GRID
    Section(os, "1. GRID");
    os << " The raster as it was collected. Positions everywhere below are in the DISPLAY grid -\n"
       << " column 0 the lowest x, row 0 the lowest y - whichever way the stage moved, so a snake\n"
       << " scan and a raster scan of the same area name the same cell the same way.\n\n";
    Key(os, "GRID_SIZE", fmt::format("{} {}", grid.GetGridSizeX_step(), grid.GetGridSizeY_step()));
    Key(os, "GRID_STEP_UM", fmt::format("{:.2f} {:.2f}", std::fabs(grid.GetGridStepX_um()),
                                        std::fabs(grid.GetGridStepY_um())));
    Key(os, "GRID_EXTENT_UM", fmt::format("{:.1f} {:.1f}", grid.GetGridSizeX_um(), grid.GetGridSizeY_um()));
    Key(os, "GRID_SNAKE", grid.IsSnakeScan() ? "TRUE" : "FALSE");
    Key(os, "GRID_VERTICAL", grid.IsVerticalScan() ? "TRUE" : "FALSE");
    Key(os, "GRID_N_ELEM", grid.GetNElem());
    Key(os, "IMAGES_PROCESSED", scan.images.size());
    Key(os, "BEAM_SIZE_UM", fmt::format("{:.2f} {:.2f}", settings.beam_size_x_um, settings.beam_size_y_um));
    Key(os, "BEAM_SIZE_SOURCE", settings.beam_size_source);

    // -------------------------------------------------------------------- 2. PER-IMAGE SCORES
    Section(os, "2. PER-IMAGE SCORES");
    os << " Every image carries a protein score and an ice score, each between 0 and 1. A cell counts\n"
       << " as diffracting above PROTEIN_SCORE_THRESHOLD, and a patch holding one grows out to\n"
       << " GROW_SCORE_THRESHOLD, so a crystal is not broken in two by a single cell that fell just\n"
       << " under the higher bar; a patch that never reaches the higher bar is not a patch at all.\n"
       << " The crystals below are the grown patches of at least MIN_BLOB_CELLS cells - or smaller\n"
       << " patches whose best cell reaches DECISIVE_PROTEIN_SCORE, so one cell is enough where that\n"
       << " cell is clearly protein.\n"
       << " ICE_ABOVE_THRESHOLD counts the cells whose ICE score passes the protein threshold - the one\n"
       << " number the run already has, rather than a second one invented for ice - so it says how much\n"
       << " of the loop is ice rather than sample.\n\n";
    Key(os, "PROTEIN_SCORE_THRESHOLD", fmt::format("{:.2f}", settings.analysis.GetProteinScoreThreshold()));
    Key(os, "GROW_SCORE_THRESHOLD", fmt::format("{:.2f}", settings.analysis.GetGrowScoreThreshold()));
    Key(os, "MIN_BLOB_CELLS", settings.analysis.GetMinBlobCells());
    Key(os, "DECISIVE_PROTEIN_SCORE", fmt::format("{:.2f}", settings.analysis.GetDecisiveSingleCellScore()));
    Key(os, "MAX_CRYSTALS", settings.analysis.GetMaxCrystals()
            ? std::to_string(*settings.analysis.GetMaxCrystals()) : std::string("none"));
    Key(os, "IMAGES_SCORED", summary.scored);
    Key(os, "IMAGES_ABOVE_THRESHOLD", summary.above);
    Key(os, "FRACTION_ABOVE_THRESHOLD", fmt::format("{:.4f}", Fraction(summary.above, summary.scored)));
    Key(os, "ICE_ABOVE_THRESHOLD", summary.ice_above);
    Key(os, "ICE_FRACTION", fmt::format("{:.4f}", Fraction(summary.ice_above, summary.scored)));
    Key(os, "MEAN_PROTEIN_SCORE", fmt::format("{:.4f}", summary.mean_protein));
    Key(os, "MEAN_ICE_SCORE", fmt::format("{:.4f}", summary.mean_ice));

    // ------------------------------------------------------------------------- 3. CRYSTALS
    Section(os, "3. CRYSTALS");
    os << " One row per crystal, best first. NX/NY and X_UM/Y_UM are the centre, the micrometres\n"
       << " measured from the centre of cell (0,0) along the grid axes; IMAGE is the nearest image\n"
       << " that was actually collected, so a DAQ can drive back to it. MAJOR_UM x MINOR_UM is the\n"
       << " extent along the crystal's own axes and ANGLE_DEG points along MAJOR_UM, counter-clockwise\n"
       << " from +x and folded into [0,180) because an axis has no sign. The sizes are MEASURED and\n"
       << " still contain the beam (BEAM_SIZE_UM above says which one) - nothing here deconvolves it.\n"
       << " RES_A is the robust best resolution inside the patch, blank where none was measured. SCORE\n"
       << " is the patch mean and ranks the rows; PEAK is its best cell, which is what admits a patch\n"
       << " with fewer than MIN_BLOB_CELLS cells in it. N_CELLS and the extents describe the GROWN\n"
       << " patch, so they reach the GROW_SCORE_THRESHOLD contour and not the seed one.\n\n";
    Key(os, "CRYSTAL_COUNT", crystals.crystals.size());
    os << "\n"
       << "        NX      NY      X_UM      Y_UM      IMAGE  MAJOR_UM  MINOR_UM  ANGLE_DEG   SCORE    PEAK  ICE_SCORE   RES_A  N_CELLS\n"
       << "  --------  ------  --------  --------  ---------  --------  --------  ---------  ------  ------  ---------  ------  -------\n";
    for (const auto &c : crystals.crystals) {
        os << fmt::format("  {:8.2f}  {:6.2f}  {:8.1f}  {:8.1f}  {:9d}  {:8.1f}  {:8.1f}  {:9.1f}  {:6.3f}  {:6.3f}  {:9.3f}  {:>6}  {:7d}\n",
                          c.nx, c.ny, c.x_um, c.y_um, c.image_number, c.major_um, c.minor_um,
                          c.angle_deg, c.score, c.peak_score, c.ice_score,
                          std::isfinite(c.res_A) ? fmt::format("{:.2f}", c.res_A) : std::string("-"),
                          c.n_images);
    }
    if (crystals.crystals.empty())
        os << " (none)\n";
    else
        os << "  --------  ------  --------  --------  ---------  --------  --------  ---------  ------  ------  ---------  ------  -------\n";

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

std::string RenderRasterJson(const std::string &input_file,
                             const GridScanSettings &grid,
                             const ScanResult &scan,
                             const GridScanResult &crystals,
                             const RasterSettings &settings,
                             const RunProvenance &provenance) {
    const ScanSummary summary = Summarize(scan, settings.analysis.GetProteinScoreThreshold());

    nlohmann::json out;
    out["raster_report_version"] = RASTER_REPORT_VERSION;
    out["rugnux_version"] = jfjoch_version();
    if (!jfjoch_git_sha1().empty())
        out["rugnux_git"] = jfjoch_git_sha1().substr(0, 6);
    out["date"] = time_UTC(std::chrono::system_clock::now());
    out["input_file"] = input_file;
    if (!provenance.command_line.empty())
        out["command_line"] = provenance.command_line;
    if (provenance.wall_time_s > 0.0)
        out["wall_time_s"] = provenance.wall_time_s;

    nlohmann::json g;
    g["size_x"] = grid.GetGridSizeX_step();
    g["size_y"] = grid.GetGridSizeY_step();
    g["step_x_um"] = std::fabs(grid.GetGridStepX_um());
    g["step_y_um"] = std::fabs(grid.GetGridStepY_um());
    g["extent_x_um"] = grid.GetGridSizeX_um();
    g["extent_y_um"] = grid.GetGridSizeY_um();
    g["snake"] = grid.IsSnakeScan();
    g["vertical"] = grid.IsVerticalScan();
    g["n_elem"] = grid.GetNElem();
    out["grid"] = g;

    nlohmann::json s;
    s["protein_score_threshold"] = settings.analysis.GetProteinScoreThreshold();
    s["grow_score_threshold"] = settings.analysis.GetGrowScoreThreshold();
    s["min_blob_cells"] = settings.analysis.GetMinBlobCells();
    s["decisive_protein_score"] = settings.analysis.GetDecisiveSingleCellScore();
    if (const auto cap = settings.analysis.GetMaxCrystals())
        s["max_crystals"] = *cap;
    else
        s["max_crystals"] = nullptr;
    s["beam_size_x_um"] = settings.beam_size_x_um;
    s["beam_size_y_um"] = settings.beam_size_y_um;
    s["beam_size_source"] = settings.beam_size_source;
    out["settings"] = s;

    nlohmann::json sc;
    sc["images_processed"] = scan.images.size();
    sc["images_scored"] = summary.scored;
    sc["images_above_threshold"] = summary.above;
    sc["fraction_above_threshold"] = Fraction(summary.above, summary.scored);
    sc["ice_above_threshold"] = summary.ice_above;
    sc["ice_fraction"] = Fraction(summary.ice_above, summary.scored);
    sc["mean_protein_score"] = summary.mean_protein;
    sc["mean_ice_score"] = summary.mean_ice;
    out["scores"] = sc;

    out["crystal_count"] = crystals.crystals.size();
    nlohmann::json list = nlohmann::json::array();
    for (const auto &c : crystals.crystals) {
        nlohmann::json j;
        j["nx"] = c.nx;
        j["ny"] = c.ny;
        j["x_um"] = c.x_um;
        j["y_um"] = c.y_um;
        j["image_number"] = c.image_number;
        j["major_um"] = c.major_um;
        j["minor_um"] = c.minor_um;
        j["angle_deg"] = c.angle_deg;
        j["score"] = c.score;
        j["peak_score"] = c.peak_score;
        j["ice_score"] = c.ice_score;
        // JSON has no NaN, so a resolution that was never measured is null rather than a number.
        if (std::isfinite(c.res_A))
            j["res_A"] = c.res_A;
        else
            j["res_A"] = nullptr;
        j["n_images"] = c.n_images;
        list.push_back(j);
    }
    out["crystals"] = list;

    return out.dump(1) + "\n";
}

void WriteRasterReport(const std::string &output_prefix,
                       const std::string &input_file,
                       const GridScanSettings &grid,
                       const ScanResult &scan,
                       const GridScanResult &crystals,
                       const RasterSettings &settings,
                       Logger &logger,
                       const RunProvenance &provenance) {
    if (output_prefix.empty())
        return;   // "analyse the raster, persist nothing"

    // As for the results report: a run that found the crystals must not fail because a side file
    // could not be written. Report the failure and carry on.
    const std::string txt_path = output_prefix + "_raster_report.txt";
    const std::string json_path = output_prefix + "_raster.json";
    try {
        std::ofstream txt(txt_path);
        txt.exceptions(std::ios::failbit | std::ios::badbit);
        txt << RenderRasterReport(input_file, grid, scan, crystals, settings, provenance);
    } catch (const std::exception &e) {
        logger.Warning("Could not write the raster report {}: {}", txt_path, e.what());
    }
    try {
        std::ofstream json(json_path);
        json.exceptions(std::ios::failbit | std::ios::badbit);
        json << RenderRasterJson(input_file, grid, scan, crystals, settings, provenance);
    } catch (const std::exception &e) {
        logger.Warning("Could not write the raster JSON {}: {}", json_path, e.what());
    }
}