Jungfraujoch/tests/ReindexAmbiguityTest.cpp

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

#include <catch2/catch_all.hpp>

#include <cmath>
#include <unordered_map>
#include <unordered_set>
#include <vector>

#include "../image_analysis/scale_merge/CrystalSetting.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
#include "../image_analysis/scale_merge/SublatticeRelabel.h"

namespace {
    const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); }

    UnitCell Tetragonal() { return UnitCell{78.0f, 78.0f, 37.0f, 90.0f, 90.0f, 90.0f}; }

    // A reference set: one distinct intensity per Laue-ASU reflection, so that reflections related by
    // a twin law (which are NOT Laue-equivalent in P4) carry different intensities.
    std::vector<MergedReflection> DistinctReference(int space_group_number) {
        const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(space_group_number));
        std::vector<MergedReflection> v;
        std::unordered_set<uint64_t> seen;
        for (int h = -8; h <= 8; ++h)
            for (int k = -8; k <= 8; ++k)
                for (int l = 0; l <= 8; ++l) {
                    if (h == 0 && k == 0 && l == 0) continue;
                    const uint64_t kk = key(h, k, l).pack();
                    if (!seen.insert(kk).second) continue;   // one row per ASU reflection
                    MergedReflection r;
                    r.h = h; r.k = k; r.l = l;
                    r.I = static_cast<float>(100 + kk % 9973);   // distinct, spread
                    r.sigma = 1.0f;
                    r.d = 60.0f / (1 + h * h + k * k + l * l);
                    v.push_back(r);
                }
        return v;
    }
}

TEST_CASE("Reindex: holohedral crystal has no indexing ambiguity", "[reindex]") {
    // P4(3)2(1)2 (422, holohedral for the tetragonal lattice) -> no twin laws.
    CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(96)).empty());
    // P422 likewise.
    CHECK(ReindexAmbiguityOperators(Tetragonal(), SG(89)).empty());
}

TEST_CASE("Reindex: merohedral crystal exposes the ambiguity operators", "[reindex]") {
    // P4 (point group 4) in a tetragonal lattice (422) -> a non-trivial reindexing coset.
    CHECK_FALSE(ReindexAmbiguityOperators(Tetragonal(), SG(75)).empty());
}

TEST_CASE("Reindex: reference agreement recovers a misindexed dataset", "[reindex]") {
    const int sg = 75;   // P4
    const auto reference = DistinctReference(sg);
    const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg));
    REQUIRE_FALSE(laws.empty());

    // Deliberately mis-index the data by one twin law.
    const auto data = ReindexReflections(reference, laws.front());

    const auto choice = ChooseReindex(
        data, Tetragonal(), SG(sg),
        [&](const std::vector<MergedReflection> &m) { return ReferenceIntensityCC(m, reference, SG(sg)); });

    CHECK_FALSE(choice.is_identity);                 // a reindex was needed
    CHECK(choice.score > 0.99);                      // the winner realigns with the reference
    CHECK(choice.identity_score < 0.9);              // leaving it mis-indexed correlates poorly
    CHECK(choice.n_candidates >= 2);
}

TEST_CASE("Reindex: a correctly indexed dataset keeps identity", "[reindex]") {
    const int sg = 75;
    const auto reference = DistinctReference(sg);
    const auto choice = ChooseReindex(
        reference, Tetragonal(), SG(sg),
        [&](const std::vector<MergedReflection> &m) { return ReferenceIntensityCC(m, reference, SG(sg)); });
    CHECK(choice.is_identity);
    CHECK(choice.score > 0.99);
}

TEST_CASE("Reindex into the ASU: a twin law permutes the reflections without losing any", "[reindex]") {
    const int sg = 75;   // P4
    const auto reference = DistinctReference(sg);
    const auto laws = ReindexAmbiguityOperators(Tetragonal(), SG(sg));
    REQUIRE_FALSE(laws.empty());

    // Mis-index by a twin law, then reindex back into the ASU: the labels must land where an export
    // needs them, and the reflection the label carries must be the one the reference has there.
    const auto misindexed = ReindexReflections(reference, laws.front());
    const auto fixed = ReindexMergedIntoAsu(misindexed, laws.front(), SG(sg), /*merge_friedel=*/true);
    REQUIRE(fixed.size() == reference.size());

    const HKLKeyGenerator key(true, *gemmi::find_spacegroup_by_number(sg));
    std::unordered_map<uint64_t, float> ref_by_key;
    for (const auto &r : reference)
        ref_by_key[key(r).pack()] = r.I;

    std::unordered_set<uint64_t> seen;
    for (const auto &r : fixed) {
        const HKLKey k = key(r);
        CHECK(k.h == r.h);                       // already at its own ASU index
        CHECK(k.k == r.k);
        CHECK(k.l == r.l);
        CHECK(seen.insert(k.pack()).second);     // no two reflections land on one label
        const auto it = ref_by_key.find(k.pack());
        REQUIRE(it != ref_by_key.end());
        CHECK(it->second == r.I);
    }
}

TEST_CASE("Reindex into the ASU: the change of hand swaps the Bijvoet halves", "[reindex]") {
    // The change of hand is the inversion, so it leaves the Laue-ASU label alone and moves the
    // anomalous signal instead - which is the whole of what adopting a model's hand does to intensities.
    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(96);   // P4(3)2(1)2
    REQUIRE(sg != nullptr);
    const HKLKeyGenerator key(true, *sg);
    const HKLKey asu = key(3, 1, 2);

    MergedReflection r;
    r.h = asu.h; r.k = asu.k; r.l = asu.l;
    r.I = 100.0f; r.sigma = 1.0f; r.d = 10.0f;
    r.I_plus = 110.0f; r.sigma_plus = 2.0f; r.I_minus = 90.0f; r.sigma_minus = 3.0f;
    r.F_plus = 10.5f; r.F_minus = 9.5f;

    const auto merged = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/true);
    REQUIRE(merged.size() == 1);
    CHECK(merged[0].h == asu.h);
    CHECK(merged[0].k == asu.k);
    CHECK(merged[0].l == asu.l);
    CHECK(merged[0].I == 100.0f);
    CHECK(merged[0].I_plus == 90.0f);
    CHECK(merged[0].I_minus == 110.0f);
    CHECK(merged[0].sigma_plus == 3.0f);
    CHECK(merged[0].sigma_minus == 2.0f);
    CHECK(merged[0].F_plus == 9.5f);
    CHECK(merged[0].F_minus == 10.5f);

    // With the mates kept apart the merge stores the minus hand at -hkl, so the row moves there.
    const auto anom = ReindexMergedIntoAsu({r}, sg->change_of_hand_op(), SG(96), /*merge_friedel=*/false);
    REQUIRE(anom.size() == 1);
    CHECK(anom[0].h == -asu.h);
    CHECK(anom[0].k == -asu.k);
    CHECK(anom[0].l == -asu.l);
    CHECK(anom[0].I == 100.0f);
}

TEST_CASE("Reindex into the ASU: both mates of an anomalous pair follow the same hand", "[reindex]") {
    // With the mates kept apart, the merge stores the plus mate at +hkl_asu and the minus mate at
    // -hkl_asu, while I_plus/I_minus are attached in the plus convention on BOTH rows alike. An
    // alternative-indexing operator is rotation-type, so op(-x) == -op(x): whatever it does to one
    // mate's hand it does to the other's, and the two rows must still agree afterwards about which
    // intensity is which. Deciding the swap from the hand the new index lands on instead of from the
    // CHANGE of hand moves exactly one of the two, and the pair ends up self-contradictory.
    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(75);   // P4, Laue 4/m
    REQUIRE(sg != nullptr);
    const HKLKeyGenerator key_gen(/*merge_friedel=*/false, *sg);
    const HKLKey asu = key_gen(3, 1, 2);
    REQUIRE(asu.plus);

    // The alternative indexing of a P4 lattice: 4/mmm holds it, 4/m does not.
    const gemmi::Op op = gemmi::parse_triplet("k,h,-l");

    MergedReflection plus;
    plus.h = asu.h; plus.k = asu.k; plus.l = asu.l;
    plus.I = 100.0f; plus.sigma = 1.0f; plus.d = 10.0f;
    plus.I_plus = 110.0f; plus.sigma_plus = 2.0f; plus.I_minus = 90.0f; plus.sigma_minus = 3.0f;
    plus.F_plus = 10.5f; plus.F_minus = 9.5f;

    MergedReflection minus = plus;                   // same anomalous split, stored at -hkl_asu
    minus.h = -asu.h; minus.k = -asu.k; minus.l = -asu.l;

    const auto out = ReindexMergedIntoAsu({plus, minus}, op, SG(75), /*merge_friedel=*/false);
    REQUIRE(out.size() == 2);

    // The two rows are still one pair: same Laue label up to the Friedel sign, opposite hands.
    CHECK(out[0].h == -out[1].h);
    CHECK(out[0].k == -out[1].k);
    CHECK(out[0].l == -out[1].l);

    // ... and they still tell the same story about which hand carries which intensity.
    CHECK(out[0].I_plus == out[1].I_plus);
    CHECK(out[0].I_minus == out[1].I_minus);
    CHECK(out[0].sigma_plus == out[1].sigma_plus);
    CHECK(out[0].sigma_minus == out[1].sigma_minus);
    CHECK(out[0].F_plus == out[1].F_plus);
    CHECK(out[0].F_minus == out[1].F_minus);

    // And the swap did have to happen: "k,h,-l" takes (3,1,2) to (1,3,-2), whose Laue class has its
    // ASU representative at l > 0, so the row that was the plus mate of its pair is the minus mate of
    // the new one. The intensity now standing at the plus index is the one that was at -hkl before.
    CHECK(out[0].I_plus == 90.0f);
    CHECK(out[0].I_minus == 110.0f);
}

TEST_CASE("Reindex: every ambiguity operator is a proper rotation", "[reindex]") {
    // A relabelling by a lattice rotation keeps the hand; one with determinant -1 would swap every
    // Bijvoet pair. The candidates come from the lattice symmetry without the inversion, so none may
    // be improper, in any crystal system that has an ambiguity.
    const std::vector<std::pair<UnitCell, int>> cases{
        {Tetragonal(), 75},                                           // P4
        {UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 143},   // P3
        {UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 149},   // P312
        {UnitCell{60.0f, 60.0f, 90.0f, 90.0f, 90.0f, 120.0f}, 168},   // P6
        {UnitCell{77.0f, 77.0f, 77.0f, 90.0f, 90.0f, 90.0f}, 197},    // I23
        {UnitCell{77.0f, 77.0f, 77.0f, 90.0f, 90.0f, 90.0f}, 199},    // I2(1)3
    };
    for (const auto &[cell, sg] : cases) {
        const auto ops = ReindexAmbiguityOperators(cell, SG(sg));
        CHECK_FALSE(ops.empty());
        for (const auto &op : ops)
            CHECK(op.det_rot() == gemmi::Op::DEN * gemmi::Op::DEN * gemmi::Op::DEN);
    }
}

TEST_CASE("Reindex: a near-symmetry law moves the cell with the indices", "[reindex]") {
    // Pseudo-merohedry: b within 2% of c in a primitive orthorhombic cell, and beta a degree off 90 in
    // a monoclinic one. Each admits an indexing law that is not a symmetry of the cell, so relabelling
    // the reflections by it permutes the axes (or flips beta), and the cell they are written with has
    // to follow - every reflection then keeps its own d-spacing.
    const std::vector<std::pair<UnitCell, int>> cases{
        {UnitCell{50.0f, 60.0f, 61.2f, 90.0f, 90.0f, 90.0f}, 16},    // P222, b ~ c
        {UnitCell{70.0f, 79.0f, 83.0f, 90.0f, 91.1f, 90.0f}, 3},     // P2, beta ~ 90
    };
    for (const auto &[cell, sg] : cases) {
        const auto laws = ReindexAmbiguityOperators(cell, SG(sg));
        REQUIRE_FALSE(laws.empty());
        bool moved = false;
        for (const gemmi::Op &law : laws) {
            const UnitCell written = CellInBasis(cell, HklOperator(law));
            moved = moved || std::abs(written.b - cell.b) > 0.5f || std::abs(written.beta - cell.beta) > 0.5f;
            const gemmi::UnitCell before = cell, after = written;
            for (const gemmi::Op::Miller hkl : {gemmi::Op::Miller{1, 2, 3}, gemmi::Op::Miller{-4, 0, 7},
                                                gemmi::Op::Miller{5, -3, 1}})
                CHECK(after.calculate_d(law.apply_to_hkl(hkl)) == Catch::Approx(before.calculate_d(hkl)).epsilon(1e-5));
        }
        CHECK(moved);
    }
}

namespace {
    // A P2(1) cell with a near-symmetry: a+c is within 1% of a, so the lattice is metrically close to
    // a C-centred orthorhombic one (obliquity ~0.7 deg) without being one.
    UnitCell PseudoMonoclinic() { return UnitCell{64.0f, 29.0f, 36.0f, 90.0f, 107.0f, 90.0f}; }

    Coord Recip(const CrystalLattice &latt, int h, int k, int l) {
        return latt.Astar() * static_cast<float>(h) + latt.Bstar() * static_cast<float>(k)
             + latt.Cstar() * static_cast<float>(l);
    }
}

TEST_CASE("Reindex: an ambiguity operator of a pseudo-symmetric monoclinic cell keeps every reflection's resolution",
          "[reindex]") {
    const UnitCell cell = PseudoMonoclinic();
    const auto laws = ReindexAmbiguityOperators(cell, SG(4));
    REQUIRE_FALSE(laws.empty());
    const CrystalLattice latt(cell);
    for (const auto &op : laws)
        for (int h = -6; h <= 6; ++h)
            for (int l = -6; l <= 6; ++l) {
                if (h == 0 && l == 0) continue;
                // As the resolver applies it: an operator that is a near-symmetry of the metric moves an
                // index to one of nearly the same length, never to a finer one.
                const auto m = op.apply_to_hkl({{h, 3, l}});
                const float d0 = 1.0f / Recip(latt, h, 3, l).Length();
                const float d1 = 1.0f / Recip(latt, m[0], m[1], m[2]).Length();
                CHECK(std::fabs(d1 - d0) / d0 < 0.03f);
            }
}

TEST_CASE("Reindex: the reindexed lattice predicts every spot where it was", "[reindex]") {
    const UnitCell cell = PseudoMonoclinic();
    const auto laws = ReindexAmbiguityOperators(cell, SG(4));
    REQUIRE_FALSE(laws.empty());
    const CrystalLattice latt(cell);
    for (const auto &op : laws) {
        const CrystalLattice moved = ReindexLattice(latt, op);
        for (int h = -4; h <= 4; ++h)
            for (int k = -4; k <= 4; ++k)
                for (int l = -4; l <= 4; ++l) {
                    const auto m = op.apply_to_hkl({{h, k, l}});
                    const Coord before = Recip(latt, h, k, l);
                    const Coord after = Recip(moved, m[0], m[1], m[2]);
                    CHECK((after - before).Length() < 1e-5f);
                }
    }
}

TEST_CASE("Sublattice: a frame on the given cell in another setting is described on the given cell's axes",
          "[reindex]") {
    // The given C222(1)-style cell, and the same lattice with b and c swapped (and a flipped to keep the
    // basis right-handed): not close axis by axis, but the same lattice.
    const UnitCell given{40.0f, 180.0f, 143.0f, 90.0f, 90.0f, 90.0f};
    const UnitCell frame{40.2f, 143.4f, 180.5f, 90.0f, 90.0f, 90.0f};
    IndexingSettings settings;
    REQUIRE_FALSE(frame.is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
    const auto on_given = CellOnReference(frame, given, settings);
    REQUIRE(on_given.has_value());
    CHECK(on_given->is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
    CHECK(on_given->b == Catch::Approx(180.5f).epsilon(0.001));
    CHECK(on_given->c == Catch::Approx(143.4f).epsilon(0.001));

    // An unrelated cell is not put on it.
    CHECK_FALSE(CellOnReference(UnitCell{55.0f, 66.0f, 77.0f, 90.0f, 90.0f, 90.0f}, given, settings).has_value());
}