Jungfraujoch/image_analysis/scale_merge/ReindexAmbiguity.cpp

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

#include "ReindexAmbiguity.h"

#include <cmath>
#include <unordered_map>
#include <utility>

#include "gemmi/twin.hpp"

#include "HKLKey.h"

namespace {
    constexpr size_t MIN_REFLECTIONS = 20;

    struct BestReindex {
        gemmi::Op op = gemmi::Op::identity();
        bool is_identity = true;
        double score = 0.0;
        double identity_score = 0.0;
    };

    // Among identity (the baseline, scored `identity_score`) and `ops`, return the operator with the
    // highest `score_op(op)`. A non-finite score never wins (too little overlap to decide). Shared by the
    // stills per-image resolver and the rotation post-merge ChooseReindex.
    BestReindex PickBestReindex(double identity_score, const std::vector<gemmi::Op> &ops,
                                const std::function<double(const gemmi::Op &)> &score_op) {
        BestReindex best;
        best.identity_score = identity_score;
        best.score = identity_score;
        for (const auto &op : ops) {
            const double sc = score_op(op);
            if (std::isfinite(sc) && (!std::isfinite(best.score) || sc > best.score)) {
                best.score = sc;
                best.op = op;
                best.is_identity = false;
            }
        }
        return best;
    }
}

std::vector<gemmi::Op> ReindexAmbiguityOperators(const UnitCell &cell, const gemmi::SpaceGroup &space_group,
                                                 double max_obliquity_deg) {
    return gemmi::find_twin_laws(static_cast<gemmi::UnitCell>(cell), &space_group, max_obliquity_deg,
                                 /*all_ops=*/false);
}

std::vector<MergedReflection> ReindexReflections(const std::vector<MergedReflection> &merged,
                                                 const gemmi::Op &op) {
    std::vector<MergedReflection> out = merged;
    for (auto &r : out) {
        const gemmi::Op::Miller h = op.apply_to_hkl({{static_cast<int>(r.h),
                                                      static_cast<int>(r.k),
                                                      static_cast<int>(r.l)}});
        r.h = h[0];
        r.k = h[1];
        r.l = h[2];
    }
    return out;
}

std::vector<MergedReflection> ReindexMergedIntoAsu(const std::vector<MergedReflection> &merged,
                                                   const gemmi::Op &op, const gemmi::SpaceGroup &space_group,
                                                   bool merge_friedel) {
    // merge_friedel=false so the generator reports which side of the Friedel pair the ASU index was
    // reached from; that sign is what the two hands have to follow.
    const HKLKeyGenerator key_gen(/*merge_friedel=*/false, space_group);

    std::vector<MergedReflection> out = merged;
    for (auto &r : out) {
        const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}});
        const HKLKey key = key_gen(h[0], h[1], h[2]);
        // The hands follow a CHANGE of hand, not the hand the new index happens to land on. An
        // anomalous merge stores the minus mate at -hkl_asu (Merge.cpp), so a row can already be on
        // the minus side before the operator is applied, while I_plus/I_minus are attached in the
        // plus convention on both mates alike (RotationScaleMerge.cpp keys the lookup at plus=true).
        // Since the operator is rotation-type, op(-x) == -op(x), and the minus mate's key.plus is
        // therefore the negation of the plus mate's - testing key.plus alone swaps exactly the rows
        // that must not move, and leaves the ones that must. With merge_friedel every row is stored
        // at the ASU representative, so was_plus is always true and this reduces to the old test.
        const bool was_plus = key_gen(r.h, r.k, r.l).plus;
        if (was_plus != key.plus) {
            std::swap(r.I_plus, r.I_minus);
            std::swap(r.sigma_plus, r.sigma_minus);
            std::swap(r.F_plus, r.F_minus);
            std::swap(r.sigmaF_plus, r.sigmaF_minus);
        }
        const bool at_asu_index = merge_friedel || key.plus;
        r.h = at_asu_index ? key.h : -key.h;
        r.k = at_asu_index ? key.k : -key.k;
        r.l = at_asu_index ? key.l : -key.l;
    }
    return out;
}

ReindexChoice ChooseReindex(const std::vector<MergedReflection> &merged,
                            const UnitCell &cell, const gemmi::SpaceGroup &space_group,
                            const std::function<double(const std::vector<MergedReflection> &)> &score,
                            double max_obliquity_deg) {
    const auto ops = ReindexAmbiguityOperators(cell, space_group, max_obliquity_deg);
    const BestReindex best = PickBestReindex(
        score(merged), ops,
        [&](const gemmi::Op &op) { return score(ReindexReflections(merged, op)); });

    ReindexChoice choice;
    choice.op = best.op;
    choice.is_identity = best.is_identity;
    choice.score = best.score;
    choice.identity_score = best.identity_score;
    choice.n_candidates = 1 + static_cast<int>(ops.size());
    return choice;
}

double ReferenceIntensityCC(const std::vector<MergedReflection> &merged,
                            const std::vector<MergedReflection> &reference,
                            const gemmi::SpaceGroup &space_group) {
    const HKLKeyGenerator key(/*merge_friedel=*/true, space_group);

    std::unordered_map<uint64_t, double> ref;
    ref.reserve(reference.size());
    for (const auto &r : reference)
        if (std::isfinite(r.I))
            ref[key(r).pack()] = r.I;

    double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
    int n = 0;
    for (const auto &m : merged) {
        if (!std::isfinite(m.I))
            continue;
        const auto it = ref.find(key(m).pack());
        if (it == ref.end())
            continue;
        const double x = m.I, y = it->second;
        sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y;
        ++n;
    }
    if (n < 10)
        return 0.0;
    const double cov = n * sxy - sx * sy;
    const double vx = n * sxx - sx * sx;
    const double vy = n * syy - sy * sy;
    return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : 0.0;
}

ReindexAmbiguityResolver::ReindexAmbiguityResolver(const DiffractionExperiment &x,
                                                   const std::vector<MergedReflection> &reference)
    : s(x.GetScalingSettings()),
      hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupOrP1()) {
    for (const auto &r : reference)
        reference_data[hkl_key_generator(r)] = r.I;
    if (x.GetUnitCell().has_value() && x.GetGemmiSpaceGroup().has_value())
        ops = ReindexAmbiguityOperators(*x.GetUnitCell(), x.GetSpaceGroupOrP1());
}

bool ReindexAmbiguityResolver::Accept(const Reflection &r) const {
    if (r.on_ice_ring)   // ice-contaminated intensity would bias the correlation; keep it out
        return false;
    return AcceptReflection(r, s.GetHighResolutionLimit_A(), s.GetLowResolutionLimit_A());
}

double ReindexAmbiguityResolver::ReferenceCC(const std::vector<Reflection> &reflections,
                                             const gemmi::Op &op) const {
    double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
    size_t n = 0;
    for (const auto &r : reflections) {
        if (!Accept(r) || r.partiality < s.GetMinPartiality())
            continue;
        if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f)
            continue;
        const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}});
        const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2]));
        if (it == reference_data.end())
            continue;
        const double x = static_cast<double>(r.I) * r.rlp / r.partiality;
        const double y = it->second;
        if (!std::isfinite(x) || !std::isfinite(y))
            continue;
        sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y;
        ++n;
    }
    if (n < MIN_REFLECTIONS)
        return NAN;
    const double nd = static_cast<double>(n);
    const double cov = sxy - sx * sy / nd;
    const double vx = sxx - sx * sx / nd;
    const double vy = syy - sy * sy / nd;
    return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN;
}

// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this
// image alone, the reindexing whose intensities correlate best with the external reference and apply it.
void ReindexAmbiguityResolver::Resolve(std::vector<Reflection> &reflections) const {
    if (ops.empty())
        return;

    const BestReindex best = PickBestReindex(
        ReferenceCC(reflections, gemmi::Op::identity()), ops,
        [&](const gemmi::Op &op) { return ReferenceCC(reflections, op); });
    if (best.is_identity)
        return;

    for (auto &r : reflections) {
        const gemmi::Op::Miller h = best.op.apply_to_hkl({{r.h, r.k, r.l}});
        r.h = h[0]; r.k = h[1]; r.l = h[2];
    }
}