Jungfraujoch/image_analysis/scale_merge/ScaleAndMerge.h

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

#pragma once

#include <vector>
#include <cstdint>
#include <optional>

#include "../../common/Reflection.h"
#include "gemmi/symmetry.hpp"

struct ScaleMergeOptions {
    int max_num_iterations = 100;
    double max_solver_time_s = 1.0;

    double image_number_rounding = 1.0;
    double min_sigma = 1e-3;

    // Symmetry canonicalization of HKL prior to merging/scaling.
    // If not set, the routine uses raw HKL as-is.
    std::optional<gemmi::SpaceGroup> space_group;

    // If true, treat Friedel mates as equivalent (merge anomalous pairs).
    // If false, keep them separate by including a sign flag in the HKL key.
    bool merge_friedel = true;

    // --- Kabsch(XDS)-style partiality model ---
    // Rotation range (wedge) used in partiality calculation.
    // Set to 0 to disable partiality correction.
    std::optional<double> wedge_deg;

    // --- Mosaicity (user input in degrees; internally converted to radians) ---
    double mosaicity_init_deg = 0.17;
    double mosaicity_min_deg = 1e-3;
    double mosaicity_max_deg = 2.0;
    std::vector<double> mosaicity_init_deg_vec;

    // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) ---
    bool regularize_scale_to_one = true;
    double scale_regularization_sigma = 0.05;

    double min_partiality_for_merge = 0.2;
    bool smoothen_g = true;
    bool smoothen_mos = true;

    double d_min_limit_A = 0.0;

    int64_t image_cluster = 1;

    bool refine_wedge = false;

    bool selection = true;
    double filter_sigma_cutoff = 1.0;
    int filter_min_per_bin = 10000; // needs optimization
    int filter_max_per_bin = 80000;
    int filter_n_resolution_bins = 20;

    enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
};

struct MergedReflection {
    int h;
    int k;
    int l;
    double I;
    double sigma;
    double d = 0.0;
};

/// Per-resolution-shell merging statistics
struct MergeStatisticsShell {
    float d_min = 0.0f;          ///< High-resolution limit of this shell (Å)
    float d_max = 0.0f;          ///< Low-resolution limit of this shell (Å)
    float mean_one_over_d2 = 0;  ///< Mean 1/d² in shell
    int total_observations = 0;  ///< Total number of (corrected) observations
    int unique_reflections = 0;  ///< Number of unique HKLs with observations
    double rmeas = 0.0;          ///< Redundancy-independent merging R-factor
    double mean_i_over_sigma = 0.0; ///< Mean I/σ(I) of the merged reflections
    double completeness = 0.0;   ///< Fraction of possible reflections observed (0 if unknown)
    int possible_reflections = 0;///< Theoretical number of reflections in this shell (0 if unknown)
};

/// Overall + per-shell merging statistics
struct MergeStatistics {
    std::vector<MergeStatisticsShell> shells;
    MergeStatisticsShell overall;  ///< Statistics over all shells combined
};

struct ScaleMergeResult {
    std::vector<MergedReflection> merged;
    std::vector<float> image_scale_g;
    std::vector<float> mosaicity_deg;

    /// Per-shell and overall merging statistics (populated after merging)
    MergeStatistics statistics;
};

ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection>>& observations,
                                              const ScaleMergeOptions& opt = {});