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

#include "TwinningAnalysis.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <limits>
#include <numbers>
#include <sstream>
#include <unordered_map>
#include <vector>

#include <Eigen/Dense>

namespace {
    int64_t PackHKL(int h, int k, int l) {
        constexpr int64_t bias = 1 << 20;          // indices assumed within +/- 2^20
        return ((h + bias) << 42) | ((k + bias) << 21) | (l + bias);
    }

    bool UsableIntensity(const MergedReflection& r) {
        return std::isfinite(r.I) && std::isfinite(r.d) && r.d > 0.0;
    }

    // Merohedral twinning needs a twin law - a lattice symmetry operation that is not a symmetry of the
    // crystal - which exists only when the Laue class is a proper subgroup of the lattice holohedry.
    // The holohedral high-symmetry Laue classes (4/mmm, 6/mmm, m-3m, and -3m on a rhombohedral lattice)
    // admit no such operation, so twinning is geometrically impossible and the intensity statistics
    // cannot be indicating it. Low-symmetry classes stay eligible because pseudo-merohedral twinning
    // through an accidental metric specialisation cannot be excluded from the symmetry alone.
    // <|L|> of a partial twin with fraction a: I = (1-a) X + a Y for independent Wilson X, Y, and the
    // expectation of |I1 - I2|/(I1 + I2) over two independent such intensities, integrated in closed
    // form (0.500 at a = 0, 0.440 at 0.1, 0.408 at 0.2, 0.375 at 0.5; checked by simulation, and
    // the relation ctruncate tabulates).
    double MeanAbsLForTwinFraction(double a) {
        const double p = 1.0 - a, q = a;
        const double A = 1.0 / q, B = 1.0 / p - 1.0 / q;
        auto F = [&](double t) { return 2.0 / (B * B) * std::log(A + B * t) + (2.0 * A / B + 1.0) / (B * (A + B * t)); };
        const double J = F(1.0) + F(0.0) - 2.0 * F(0.5);
        return ((p * p + q * q) / 2.0 - 2.0 * J) / ((p - q) * (p - q));
    }

    // The inverse, by bisection: the relation is monotonic, and flat near a = 0.5, where <|L|> can no
    // longer tell 0.4 from 0.5.
    double TwinFractionFromMeanAbsL(double mean_abs_l) {
        double lo = 1e-4, hi = 0.4999;
        if (mean_abs_l >= MeanAbsLForTwinFraction(lo))
            return 0.0;
        if (mean_abs_l <= MeanAbsLForTwinFraction(hi))
            return 0.5;
        for (int i = 0; i < 50; ++i) {
            const double mid = 0.5 * (lo + hi);
            (MeanAbsLForTwinFraction(mid) > mean_abs_l ? lo : hi) = mid;
        }
        return 0.5 * (lo + hi);
    }

    bool MerohedralTwinningPossible(const gemmi::SpaceGroup* sg) {
        if (!sg)
            return true;                    // P1 / unknown symmetry: cannot rule twinning out
        switch (sg->laue_class()) {
            case gemmi::Laue::L4mmm:        // 4/mmm - tetragonal holohedry
            case gemmi::Laue::L6mmm:        // 6/mmm - hexagonal holohedry
            case gemmi::Laue::Lm3m:         // m-3m  - cubic holohedry
                return false;
            case gemmi::Laue::L3m:          // -3m is holohedral on a rhombohedral (R) lattice, but a
                return sg->hm[0] != 'R';    //   hexagonal-P 32/3m crystal can still twin towards 6/mmm
            default:
                return true;
        }
    }
}

TwinningAnalysisResult AnalyzeTwinning(const std::vector<MergedReflection>& merged,
                                       const gemmi::SpaceGroup* space_group, int resolution_shells,
                                       const std::array<double, 3>* tncs_vector, char lattice_centring) {
    TwinningAnalysisResult result;
    if (merged.empty())
        return result;

    // Centric reflections follow different statistics and must be excluded. In P1 (no space group)
    // none are centric.
    const gemmi::GroupOps gops = space_group ? space_group->operations() : gemmi::GroupOps{};
    auto acentric = [&](const MergedReflection& r) {
        return !space_group || !gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}});
    };
    // A merge in P1 of a centred lattice still holds the reflections the centring extinguishes (an
    // R lattice indexed in its hexagonal cell keeps the two thirds with -h+k+l != 3n). They are
    // noise, not intensities, and an L-test partner step that crosses into them pairs a real
    // reflection with an absent one. So they are left out of both statistics.
    gemmi::GroupOps centring;
    centring.sym_ops = {gemmi::Op::identity()};
    centring.cen_ops = gemmi::centring_vectors(lattice_centring);
    auto used = [&](const MergedReflection& r) {
        return UsableIntensity(r) && acentric(r)
               && !centring.is_systematically_absent(gemmi::Op::Miller{{r.h, r.k, r.l}});
    };

    // --- Resolution shells: the normalising mean of both statistics ---
    // Binning by 1/d^2 removes the resolution fall-off, so the second moment is 2.0 (untwinned) or 1.5
    // (perfect twin) regardless of the overall B-factor, and an L-test pair compares two intensities
    // on the same scale. The moment divides by the *square* of the shell-mean intensity, so it is not
    // robust: on weak or mis-integrated data a shell mean can collapse to the noise floor and one
    // outlier reflection then dominates (I/mean)^2 (a single I=158 in a mean~1 shell contributed 78% of
    // a whole dataset's value). To keep these twinning indicators rather than data-quality artefacts -
    // as phenix.xtriage does - we skip noise-only shells (<I/sigma> below 1) and reject Wilson outliers
    // (E^2 above 8, ~exp(-8) upper tail) with one shell-mean re-iteration so the outlier does not
    // corrupt the normalising mean either.
    //
    // The shells are equal in width in 1/d^2. Equal-COUNT shells were measured as the alternative
    // (the first equal-width shell spans a wide range of d, and its fall-off inflates the moment):
    // they move the moment the wrong way on the same data, widening its gap to phenix.xtriage by
    // 0.01-0.04 rather than closing it, so the shells were kept as they are.
    constexpr double min_shell_isig = 1.0;     // shells below this are noise, not signal
    constexpr double wilson_outlier_e2 = 8.0;  // reject improbably strong reflections (P ~ e^-8)
    const int n_shells = std::max(1, resolution_shells);
    double min_s = std::numeric_limits<double>::infinity();
    double max_s = -std::numeric_limits<double>::infinity();
    for (const auto& r : merged) {
        if (!used(r))
            continue;
        const double s = 1.0 / (r.d * r.d);
        min_s = std::min(min_s, s);
        max_s = std::max(max_s, s);
    }
    if (!std::isfinite(min_s) || !(max_s > min_s))
        return result;
    auto shell_of = [&](double d) {
        const double t = (1.0 / (d * d) - min_s) / (max_s - min_s);
        return std::min(n_shells - 1, std::max(0, static_cast<int>(t * n_shells)));
    };
    // Group intensities by shell, and accumulate <I/sigma> to gauge each shell's signal.
    std::vector<std::vector<double>> shell_I(n_shells);
    std::vector<double> shell_isig_sum(n_shells, 0.0);
    std::vector<int> shell_isig_n(n_shells, 0);
    for (const auto& r : merged) {
        if (!used(r))
            continue;
        const int b = shell_of(r.d);
        shell_I[b].push_back(r.I);
        if (std::isfinite(r.sigma) && r.sigma > 0.0) {
            shell_isig_sum[b] += r.I / r.sigma;
            shell_isig_n[b] += 1;
        }
    }
    // The mean intensity of each shell that carries signal; 0 marks a shell that does not.
    std::vector<double> shell_mean(n_shells, 0.0);
    for (int b = 0; b < n_shells; ++b) {
        const auto& intensities = shell_I[b];
        if (intensities.empty() || shell_isig_n[b] == 0
            || shell_isig_sum[b] / shell_isig_n[b] < min_shell_isig)
            continue;
        double sum = 0.0;
        for (double I : intensities)
            sum += I;
        const double mean = sum / intensities.size();
        if (mean <= 0.0)
            continue;
        // Re-fit the mean over the reflections that pass the outlier cut, so the outlier does not
        // inflate the very mean it is measured against.
        sum = 0.0;
        int n_kept = 0;
        for (double I : intensities)
            if (I / mean <= wilson_outlier_e2) {
                sum += I;
                ++n_kept;
            }
        if (n_kept > 0 && sum > 0.0)
            shell_mean[b] = sum / n_kept;
    }

    // --- Second moment <I^2>/<I>^2 of acentric intensities, normalised per resolution shell ---
    double sum_e4 = 0.0;
    int n_moment = 0;
    for (int b = 0; b < n_shells; ++b) {
        if (shell_mean[b] <= 0.0)
            continue;
        for (double I : shell_I[b]) {
            const double e2 = I / shell_mean[b];
            if (e2 > wilson_outlier_e2)
                continue;
            sum_e4 += e2 * e2;
            ++n_moment;
        }
    }
    result.moment_reflections = n_moment;
    if (n_moment > 0)
        result.second_moment = sum_e4 / n_moment;

    // --- L-test ---
    // Following Padilla & Yeates (2003) Acta Cryst. D59, 1124-1130
    // Pair each reflection with a symmetry-independent neighbour two steps away along an axis (the
    // step of 2 keeps the partner local in resolution while avoiding the reflection itself). The
    // merged reflections are unique in the asymmetric unit, so any other merged reflection is
    // genuinely non-equivalent - exactly the pairing the L-test wants. Only acentric reflections
    // with positive intensity enter, which also keeps L = (I1-I2)/(I1+I2) bounded in [-1, 1].
    //
    // The test assumes the two members of a pair have the same EXPECTED intensity, and two index
    // steps are not the same resolution: on a small cell with a large B the pair's expected
    // intensities differ by up to a factor ~3 at 2 A, which pushes <|L|> up by as much as 0.08 (an
    // untwinned crystal read 0.568, a partial twin 0.427 where phenix.xtriage reads 0.486 and 0.360).
    // So each intensity is divided by its shell mean first, as xtriage and ctruncate do, and only the
    // shells the second moment reads enter. The selection stays on the SHELL: cutting individual
    // reflections on their own I/sigma removes the weak tail of every shell and biases <|L|> down
    // (the same untwinned crystal reads 0.423 with I/sigma > 2).
    std::unordered_map<int64_t, double> intensity;
    intensity.reserve(merged.size() * 2);
    for (const auto& r : merged) {
        if (!used(r) || r.I <= 0.0)
            continue;
        const double mean = shell_mean[shell_of(r.d)];
        if (mean > 0.0)
            intensity.emplace(PackHKL(r.h, r.k, r.l), r.I / mean);
    }

    // Step to a nearby, symmetry-independent partner. The axis step of 2 is load-bearing for TWO
    // reasons, and only the first was ever written down.
    //
    // (1) It preserves the reflection condition of P/I/C/F/A/B lattices (parity-based). See below.
    //
    // (2) It preserves the class of a HALF-INTEGER translational pseudo-symmetry. A pseudo-
    //     translation u multiplies the intensity by |1 + exp(2 pi i h.u)|^2; a partner at h + s
    //     carries the same factor exactly when s.u is an integer, which for an axis step of 2 means
    //     every component of u is 0 or 1/2. The commonest tNCS - a second copy at (1/2,1/2,1/2) or
    //     similar - therefore leaves <|L|> untouched, and the L-test works on it by construction
    //     rather than by luck. Measured on a real half-integer case: <|L|> = 0.4755 with these steps,
    //     0.4880 with parity-breaking steps of 1. Do not "shorten" or "lengthen" these steps without
    //     re-measuring that: a step of 6, which would also survive a one-third pseudo-translation,
    //     puts <|L|> above 0.500 on two thirds of a 137-dataset corpus purely because the partner is
    //     too far away in resolution, which is worse than the problem it solves.
    //
    //     A pseudo-translation that is NOT half-integer is a different matter, and which STEP is used
    //     decides it: only the steps with s.u not an integer break the class. A vector along c alone,
    //     u = (0,0,1/3), is paired by (2,0,0) on almost every reflection and is harmless; a vector
    //     with a one-third component along a breaks (2,0,0) and biases <|L|> hard. Measured on
    //     synthetic data, a PERFECT TWIN carrying a one-third pseudo-translation along a reads
    //     <|L|> = 0.454 instead of 0.377 while its second moment rises from 1.50 to 1.97 - both
    //     indicators destroyed and the twin call lost, with nothing in the report to say why.
    //
    //     So when a pseudo-translation is known, the partner steps are restricted to those that
    //     preserve its class. That is a repair, not a compromise: on the same synthetic twin it
    //     returns <|L|> = 0.376 against 0.377 for the tNCS-free control, at the cost of about 0.4%
    //     of the pairs. When no step qualifies - u = (1/3,1/3,1/3), say, where every one of the five
    //     breaks the class - the L-test simply cannot be measured on this crystal, and the verdict
    //     stops reading it in either direction.
    //
    // The axis step violates R-centring (-h+k+l = 0 mod 3,
    // as 2 != 0 mod 3) -> the partner is systematically absent and rhombohedral crystals yield zero
    // pairs. The diagonal (1,1,0)/(1,1,3) steps preserve the mod-3 condition in BOTH obverse and
    // reverse settings; they are tried only when the axis steps find no present partner, so P/I/C/F
    // behaviour is unchanged (first present partner wins).
    const std::vector<std::array<int, 3>> all_steps{
        {2, 0, 0}, {0, 2, 0}, {0, 0, 2}, {1, 1, 0}, {1, 1, 3}};

    // A partner at h + s carries the same pseudo-translation factor as h exactly when s.u is an
    // integer. Keep the steps for which it is, in the original order so the "first present partner
    // wins" behaviour is unchanged among those that remain.
    std::vector<std::array<int, 3>> steps = all_steps;
    if (tncs_vector != nullptr) {
        std::vector<std::array<int, 3>> preserving;
        for (const auto& s : all_steps) {
            const double t = s[0] * (*tncs_vector)[0] + s[1] * (*tncs_vector)[1]
                             + s[2] * (*tncs_vector)[2];
            if (std::fabs(t - std::round(t)) <= 0.06)
                preserving.push_back(s);
        }
        if (preserving.empty())
            result.l_test_contaminated_by_tncs = true;
        else {
            steps = preserving;
            result.l_test_tncs_step_restricted = true;
        }
    }

    auto run_l_test = [&](const std::vector<std::array<int, 3>>& use) {
        double sum_abs_l = 0.0, sum_l2 = 0.0;
        int n_pairs = 0;
        for (const auto& [key, i1] : intensity) {
            const int h = static_cast<int>((key >> 42) & 0x1FFFFF) - (1 << 20);
            const int k = static_cast<int>((key >> 21) & 0x1FFFFF) - (1 << 20);
            const int l = static_cast<int>(key & 0x1FFFFF) - (1 << 20);
            for (const auto& s : use) {
                const auto it = intensity.find(PackHKL(h + s[0], k + s[1], l + s[2]));
                if (it == intensity.end())
                    continue;
                const double lstat = (i1 - it->second) / (i1 + it->second);
                sum_abs_l += std::fabs(lstat);
                sum_l2 += lstat * lstat;
                ++n_pairs;
                break;                                 // one neighbour per reflection
            }
        }
        result.l_test_pairs = n_pairs;
        result.mean_abs_l = n_pairs > 0 ? sum_abs_l / n_pairs : 0.0;
        result.mean_l_squared = n_pairs > 0 ? sum_l2 / n_pairs : 0.0;
    };

    run_l_test(steps);
    // A restricted step set can leave nothing to pair with - an R-centred lattice whose only present
    // partners are the diagonals, and those are the steps the pseudo-translation broke. Fall back to
    // the full set and mark the statistic unreadable rather than reporting <|L|> from a handful of
    // pairs.
    if (result.l_test_tncs_step_restricted && result.l_test_pairs < 100) {
        result.l_test_tncs_step_restricted = false;
        result.l_test_contaminated_by_tncs = true;
        run_l_test(all_steps);
    }

    // Either indicator dropping clearly below its untwinned value is suspicious - in every Laue
    // class, holohedral ones included. A holohedral class admits no twin law WITHIN it, so there a low
    // reading does not mean "twinned in this group"; it means an operator of the group averages
    // intensities that are not equal. Merging I(h) with I(Th) under a false operator T gives
    // (I(h) + I(Th))/2 for every twin fraction - the fraction cancels - so a merge under a false
    // operator has the perfect-twin distribution by construction, whether the crystal is an untwinned
    // pseudo-symmetric one, a partial twin or a perfect twin, and only a genuine operator leaves it
    // at the untwinned 0.5. Noise cannot fake it (with the I > 0 selection it pushes <|L|> towards
    // ~0.43, not 0.375) and a pseudo-translation pushes the second moment UP. Measured over 42
    // holohedral merges of the corpus, the genuine ones read <|L|> 0.436-0.513 and the one
    // over-promotion (an H3 twin merged in R32) 0.374. One more reads low, 0.365: a crystal with a
    // ~490 A axis, whose neighbouring reflections overlap and share counts - averaging of another
    // kind, which narrows the distribution the same way. That is why this is a warning and never a
    // veto on the space group.
    result.merohedral_twinning_possible = MerohedralTwinningPossible(space_group);
    // The two indicators can only move one way under twinning: <|L|> down from 0.500 towards 0.375,
    // the second moment down from 2.0 towards 1.5. A narrow intensity distribution sitting next to an
    // <|L|> at or ABOVE its untwinned value therefore has some other cause, and calling it a twin is
    // the one reading the data rule out. It happens on small-cell rotation data, where a twin fraction
    // of 0.50 was reported for a crystal whose <|L|> was 0.63 - a value a twin cannot produce.
    // A high <|L|> here is NOT evidence of a centrosymmetric structure. Measured across the corpus,
    // the highest <|L|> of all - 0.712, second moment 3.089 - belongs to a small (~10 A) cell in an
    // orthorhombic group whose three 2_1 axes are each proven by absences, no violations against a
    // control class of the same size, so it is chiral. Both numbers sit ABOVE the centric
    // expectations of 0.637 and 3.0, which is the point: no Wilson distribution, centric or acentric,
    // reaches them, so they are not a statement about the structure's symmetry at all. What does put
    // them there is not established here, but it is neither of the obvious guesses - a structure of
    // few atoms gives a second moment of 2 - 1/N, BELOW 2, and a large term common to every structure
    // factor drives it towards 1. Read this statistic as "not a twin", and nothing further.
    //
    // Unless the L-test could not be measured free of a pseudo-translation. Three real crystals
    // carrying a one-third pseudo-translation read <|L|> = 0.527, 0.529 and 0.587, against 0.468-0.491
    // for clean untwinned controls - every one of them trips the >= 0.50 branch and has its twin test
    // vetoed for a reason that has nothing to do with twinning. Where the step restriction above could
    // not repair it, the statistic is dropped from the verdict in BOTH directions: it can no longer
    // say "not a twin", and it can no longer say "a twin" either. The second moment then decides
    // alone, and the text says the L-test could not be read. This REMOVES a veto; it never adds one.
    const bool l_test_usable = result.l_test_pairs > 0 && !result.l_test_contaminated_by_tncs;
    const bool l_test_contradicts_twin = l_test_usable && result.mean_abs_l >= 0.50;
    const bool l_test_low = l_test_usable && result.mean_abs_l < 0.44;
    const bool moment_low = result.moment_reflections > 0 && result.second_moment < 1.85;
    if (result.merohedral_twinning_possible)
        result.twinning_suspected = !l_test_contradicts_twin && (l_test_low || moment_low);
    else {
        // In a holohedral class the alternative is not a partial twin anywhere between 0.5 and 0.375
        // but a false operator, which puts <|L|> at the perfect-twin 0.375 itself - so the bar sits
        // nearer that value. Measured on the corpus (normalised L-test): the weakest genuine
        // holohedral crystal reads 0.436, the over-promoted twin 0.374. The second moment does not
        // enter here: it is the statistic a pseudo-translation inflates (that same twin read 2.08,
        // the untwinned value), and it is the one a pseudo-cubic lattice in a 4/mmm group, whose twin
        // law lies OUTSIDE the Laue class, pulls down (1.66 beside an untwinned <|L|>).
        result.adopted_operators_suspect = l_test_usable && result.mean_abs_l < 0.42;
        result.twinning_suspected = result.adopted_operators_suspect;
    }

    // The twin fraction, from the statistic that carries the verdict. Merged under a false operator
    // the fraction is lost (see above), so none is quoted there. Otherwise the L-test, which is local
    // and so indifferent to anisotropy and to a pseudo-translation it has been repaired for, unless
    // the call rests on the second moment alone. The second moment is not quoted under a detected
    // pseudo-translation, which inflates it: an H3 twin with tNCS read 2.08, the untwinned value,
    // beside a perfect-twin <|L|>.
    if (!result.adopted_operators_suspect) {
        if (l_test_usable && (l_test_low || !moment_low)) {
            result.estimated_twin_fraction = TwinFractionFromMeanAbsL(result.mean_abs_l);
            result.twin_fraction_source = TwinFractionSource::LTest;
        } else if (result.moment_reflections > 0 && tncs_vector == nullptr) {
            // M = 2(1 - a + a^2): a = (1 - sqrt(2M-3))/2.
            const double m = std::clamp(result.second_moment, 1.5, 2.0);
            result.estimated_twin_fraction = (1.0 - std::sqrt(std::max(0.0, 2.0 * m - 3.0))) / 2.0;
            result.twin_fraction_source = TwinFractionSource::SecondMoment;
        }
    }
    return result;
}

std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) {
    std::ostringstream os;
    os << std::fixed << std::setprecision(3);
    os << "Twinning analysis\n";
    if (result.l_test_pairs > 0) {
        os << "  L-test (Padilla-Yeates): <|L|> = " << result.mean_abs_l
           << ", <L^2> = " << result.mean_l_squared
           << "  [untwinned 0.500 / 0.333, perfect twin 0.375 / 0.200; " << result.l_test_pairs
           << " pairs]\n";
        if (result.l_test_contaminated_by_tncs)
            os << "     NOT READ: this crystal has a translational pseudo-symmetry, and no partner\n"
               << "     reflection available to this test shares the class that pseudo-translation\n"
               << "     defines, so <|L|> is biased upwards by it. The number above is not a statement\n"
               << "     about twinning in either direction; the verdict below rests on the second\n"
               << "     moment alone.\n";
        else if (result.l_test_tncs_step_restricted)
            os << "     Measured against partner reflections chosen to share the class of this\n"
               << "     crystal's pseudo-translation, so the pseudo-symmetry does not bias it.\n";
    }
    if (result.moment_reflections > 0)
        os << "  Second moment <I^2>/<I>^2 = " << result.second_moment
           << "  [untwinned 2.00, perfect twin 1.50]\n";
    const bool l_test_read = result.l_test_pairs > 0 && !result.l_test_contaminated_by_tncs;
    if (result.adopted_operators_suspect)
        os << "  => An adopted operator relates UNEQUAL intensities. The Laue class is holohedral, so no\n"
           << "     twin law exists within it - yet the intensities read like a twin's, and merging under\n"
           << "     a false operator produces exactly that distribution, whatever the twin fraction.\n"
           << (result.laue_class_was_chosen_by_promotion
                       ? "     This class was chosen by the space-group search: the promotion is suspect.\n"
                       : "     The space group in use is higher than these intensities support.\n")
           << "     Either the crystal has lower symmetry, or it is a twin whose law the point group has\n"
           << "     absorbed; either way it has to be refined in a subgroup. The twin fraction cannot be\n"
           << "     measured on this merge.\n";
    else if (result.twinning_suspected) {
        os << "  => Twinning suspected";
        if (result.twin_fraction_source == TwinFractionSource::LTest)
            os << " (estimated twin fraction ~" << result.estimated_twin_fraction << ", from the L-test)";
        else if (result.twin_fraction_source == TwinFractionSource::SecondMoment)
            os << " (estimated twin fraction ~" << result.estimated_twin_fraction
               << ", from the second moment)";
        else
            os << " on the second moment; no twin fraction is quoted, because this crystal's\n"
               << "     pseudo-translation inflates the second moment and <|L|> does not agree";
        os << ".\n";
        if (l_test_read && result.mean_abs_l < 0.44 && result.moment_reflections > 0
            && result.second_moment >= 1.85)
            os << "     The second moment does not agree (" << result.second_moment << ")"
               << (result.l_test_tncs_step_restricted
                           ? " - this crystal's translational pseudo-symmetry is a known reason for it\n"
                             "     to run high, so the L-test carries the verdict here.\n"
                           : ".\n");
        os << "     Statistics flag the presence of twinning, not the twin law; the law itself is a\n"
           << "     question for a dedicated twin-law analysis.\n";
    } else if (!result.merohedral_twinning_possible && result.laue_class_was_chosen_by_promotion)
        os << "  => No twinning indicated. The Laue class is holohedral, so no twin law exists within it,\n"
           << "     and it was chosen by the space-group search - but its operators read as genuine: a\n"
           << "     merge under a false operator has the perfect-twin distribution, not the untwinned one\n"
           << "     measured here. An operator that is only APPROXIMATE (non-crystallographic symmetry\n"
           << "     along it) also passes; the subgroup statistics reported by the search speak to that.\n";
    else if (!result.merohedral_twinning_possible)
        os << "  => No twinning: the Laue class is holohedral, so no merohedral twin law exists, and no\n"
           << "     operator of it averages unequal intensities.\n";
    else if (result.l_test_contaminated_by_tncs)
        os << "  => No twinning indicated by the second moment. The L-test, which is normally the\n"
           << "     stronger of the two, could not be read on this crystal (see above), so this is a\n"
           << "     weaker statement than usual - where the space-group search reported subgroup\n"
           << "     statistics, those are the stronger evidence.\n";
    else
        os << "  => No twinning indicated.\n";
    return os.str();
}

std::string TwinningVerdictLine(const TwinningAnalysisResult& result) {
    std::ostringstream os;
    os << std::fixed << std::setprecision(3);
    const bool l_test_read = result.l_test_pairs > 0 && !result.l_test_contaminated_by_tncs;
    if (result.adopted_operators_suspect)
        os << "SYMMETRY SUSPECT (<|L|> " << result.mean_abs_l << ", second moment "
           << std::setprecision(2) << result.second_moment
           << ": an adopted operator averages unequal intensities)";
    else if (result.twinning_suspected) {
        os << "INDICATED (<|L|> " << result.mean_abs_l;
        if (result.twin_fraction_source != TwinFractionSource::None)
            os << std::setprecision(2) << ", twin fraction ~" << result.estimated_twin_fraction << " from the "
               << (result.twin_fraction_source == TwinFractionSource::LTest ? "L-test" : "second moment");
        os << ")";
    } else if (l_test_read)
        os << "no indication (<|L|> " << result.mean_abs_l << ")";
    else
        os << "no indication (second moment " << std::setprecision(2) << result.second_moment
           << "; <|L|> not readable)";
    return os.str();
}

namespace {
    // Measurement noise inflates <|E^2-1|> towards centric in every class (by +0.44 at a shell
    // <I/sigma> of 2.5 and ~0 at 16, measured), so only shells above this are read. At 5 the acentric
    // control of four reference merges reads 0.65-0.74.
    constexpr double ZONE_MIN_SHELL_ISIG = 5.0;
    // Each class is normalised against its own mean in resolution bins of this many reflections, so a
    // zone that is a plane through reciprocal space is not judged on another direction's fall-off.
    constexpr size_t ZONE_BIN = 100;
    // A reflection this far above its bin mean is not a Wilson draw (centric P ~ 1e-5): one such
    // reflection read 2.33 for a whole zone.
    constexpr double ZONE_OUTLIER_E2 = 20.0;
    // What CentricOverAcentric averages over an error-free acentric population: with x = E^2 ~ Exp(1),
    // <-ln(2 pi x)/2 + x/2> = (1 + gamma_E - ln 2 pi)/2 = -0.130 nats (a centric one reads +0.216).
    // Errors move it towards zero, so an error-free baseline takes the larger excess off a noisy
    // control - the conservative side.
    constexpr double ACENTRIC_EVIDENCE_PER_REFLECTION = -0.1303;

    // ln p_centric(E^2) - ln p_acentric(E^2) for a measured E^2 with error s: both Wilson densities
    // convolved with the measurement error. Without the convolution the centric density, infinite at
    // zero, reads every noisy weak reflection as decisive; flooring E^2 at s instead turns the centric
    // excess of weak reflections into moderate values and read a genuine centric zone of a 1.2 A
    // crystal as acentric. Integrated over u = sqrt(E^2_true), where neither density is singular,
    // across +-6 s of the measurement (trapezoid).
    double CentricOverAcentric(double e2, double s) {
        constexpr int n = 200;
        const double u_lo = std::sqrt(std::max(0.0, e2 - 6.0 * s));
        const double u_hi = std::sqrt(std::max(0.0, e2 + 6.0 * s));
        if (!(u_hi > u_lo))
            return 0.0;
        const double du = (u_hi - u_lo) / n;
        double p_centric = 0.0, p_acentric = 0.0;
        for (int i = 0; i <= n; ++i) {
            const double u = u_lo + i * du;
            const double w = (i == 0 || i == n ? 0.5 : 1.0) * std::exp(-0.5 * std::pow((e2 - u * u) / s, 2));
            p_centric += w * std::sqrt(2.0 / std::numbers::pi) * std::exp(-0.5 * u * u);
            p_acentric += w * 2.0 * u * std::exp(-u * u);
        }
        return p_centric > 0.0 && p_acentric > 0.0 ? std::log(p_centric / p_acentric) : 0.0;
    }

    TwinImmuneZone ReadZone(std::vector<const MergedReflection*> refl, const std::array<double, 3>* tncs) {
        TwinImmuneZone z;
        // A pseudo-translation u multiplies the intensity by 1 + f cos(2 pi h.u): normalising within
        // each sign of the cosine keeps the two populations from reading as one broad distribution.
        std::vector<const MergedReflection*> classes[2];
        for (const auto* r : refl) {
            const double c = tncs == nullptr ? 1.0
                    : std::cos(2.0 * std::numbers::pi * (r->h * (*tncs)[0] + r->k * (*tncs)[1] + r->l * (*tncs)[2]));
            classes[c >= 0.0 ? 0 : 1].push_back(r);
        }
        double sum = 0.0, sum2 = 0.0;
        for (auto& cls : classes) {
            std::sort(cls.begin(), cls.end(), [](const auto* a, const auto* b) { return a->d > b->d; });
            const size_t n_bins = std::max<size_t>(1, cls.size() / ZONE_BIN);
            for (size_t b = 0; b < n_bins; ++b) {
                const size_t first = b * cls.size() / n_bins, last = (b + 1) * cls.size() / n_bins;
                double mean = 0.0;
                for (size_t i = first; i < last; ++i)
                    mean += cls[i]->I;
                mean /= static_cast<double>(last - first);
                if (mean <= 0.0)
                    continue;
                double kept = 0.0;
                int n_kept = 0;
                for (size_t i = first; i < last; ++i)
                    if (cls[i]->I / mean <= ZONE_OUTLIER_E2) {
                        kept += cls[i]->I;
                        ++n_kept;
                    }
                if (n_kept == 0 || kept <= 0.0)
                    continue;
                mean = kept / n_kept;
                for (size_t i = first; i < last; ++i) {
                    const double e2 = cls[i]->I / mean;
                    if (e2 > ZONE_OUTLIER_E2)
                        continue;
                    const double v = std::fabs(e2 - 1.0);
                    sum += v;
                    sum2 += v * v;
                    ++z.n;
                    z.evidence_nats += CentricOverAcentric(e2, cls[i]->sigma / mean);
                }
            }
        }
        if (z.n > 0) {
            z.mean_abs_e2_minus_1 = sum / z.n;
            const double var = std::max(0.0, sum2 / z.n - z.mean_abs_e2_minus_1 * z.mean_abs_e2_minus_1);
            z.standard_error = std::sqrt(var / z.n);
        }
        return z;
    }

    // The control's excess over the acentric expectation is the normalisation's (see
    // TwinImmuneZoneResult::control_excess_per_reflection); every zone carries it per reflection.
    void Calibrate(TwinImmuneZoneResult& result) {
        if (result.control.n > 0)
            result.control_excess_per_reflection = std::max(
                    0.0, result.control.evidence_nats / result.control.n - ACENTRIC_EVIDENCE_PER_REFLECTION);
        result.control.calibrated_evidence_nats =
                result.control.evidence_nats - result.control.n * result.control_excess_per_reflection;
        for (auto& z : result.zones)
            z.calibrated_evidence_nats = z.evidence_nats - z.n * result.control_excess_per_reflection;
    }
}

namespace {
    // The rotations of a group, as gemmi Ops without translations (its symmorphic representative's
    // proper operations).
    std::vector<gemmi::Op> ProperRotations(const gemmi::GroupOps& gops) {
        std::vector<gemmi::Op> rotations;
        for (const auto& op : gops.derive_symmorphic().sym_ops)
            if (op.det_rot() > 0)
                rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, 'x'});
        return rotations;
    }

    // An isotropic normalisation leaves an anisotropic crystal's E^2 a mixture of scales - at 2.5 A a
    // 67 A^2 anisotropy spreads one shell's expected intensity over a factor of 15 between its
    // directions - and every class then reads centric: the acentric control of such a crystal read
    // 0.99 against its 0.74. Nor can the control calibrate that away, because a zone is a plane
    // through reciprocal space and carries a different mix of directions from the control's sphere
    // (a synthetic P4 twin with 60 A^2 along c*: control 0.77, zone +49 nats after the calibration).
    // So the anisotropy is taken out first: ln I = c + s^T Q s on the acentric reflections of each
    // shell that read above 2 sigma, by least squares, and I and sigma divided by exp(s^T Q_dev s),
    // Q_dev the deviatoric part - the isotropic part is the shell mean, normalised per bin later.
    // The B convention is I ~ exp(-B |s|^2 / 2) with |s| = 1/d, so delta B = 2 (Q_max - Q_min).
    void RemoveAnisotropy(std::vector<MergedReflection>& strong, const std::vector<int>& shell,
                          const gemmi::UnitCell& cell, const gemmi::GroupOps& gops, TwinImmuneZoneResult& result) {
        auto s_of = [&](const MergedReflection& r) {
            return cell.frac.mat.left_multiply(gemmi::Vec3(r.h, r.k, r.l));
        };
        auto monomials = [](const gemmi::Vec3& s) {
            return Eigen::Matrix<double, 7, 1>{1.0, s.x * s.x, s.y * s.y, s.z * s.z,
                                               2.0 * s.x * s.y, 2.0 * s.x * s.z, 2.0 * s.y * s.z};
        };
        std::vector<size_t> fit;
        std::vector<double> shell_sum(*std::max_element(shell.begin(), shell.end()) + 1, 0.0);
        std::vector<int> shell_n(shell_sum.size(), 0);
        for (size_t i = 0; i < strong.size(); ++i) {
            const auto& r = strong[i];
            if (r.I > 2.0 * r.sigma && !gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}})) {
                fit.push_back(i);
                shell_sum[shell[i]] += r.I;
                shell_n[shell[i]] += 1;
            }
        }
        if (fit.size() < ZONE_BIN)          // seven parameters want more than a bin's worth of points
            return;
        Eigen::Matrix<double, 7, 7> ata = Eigen::Matrix<double, 7, 7>::Zero();
        Eigen::Matrix<double, 7, 1> aty = Eigen::Matrix<double, 7, 1>::Zero();
        for (const size_t i : fit) {
            const auto x = monomials(s_of(strong[i]));
            ata += x * x.transpose();
            aty += x * std::log(strong[i].I / (shell_sum[shell[i]] / shell_n[shell[i]]));
        }
        const Eigen::Matrix<double, 7, 1> c = ata.ldlt().solve(aty);
        Eigen::Matrix3d q;
        q << c[1], c[4], c[5], c[4], c[2], c[6], c[5], c[6], c[3];
        q -= q.trace() / 3.0 * Eigen::Matrix3d::Identity();
        const Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(q);
        result.anisotropy_delta_b_A2 = 2.0 * (eig.eigenvalues().maxCoeff() - eig.eigenvalues().minCoeff());
        for (auto& r : strong) {
            const gemmi::Vec3 s = s_of(r);
            const Eigen::Vector3d v(s.x, s.y, s.z);
            const double f = std::exp(v.dot(q * v));
            r.I = static_cast<float>(r.I / f);
            r.sigma = static_cast<float>(r.sigma / f);
        }
    }

    // The reflections a zone test can be read on: present in the group, general position, in the
    // shells with signal enough (ZONE_MIN_SHELL_ISIG), copied with the anisotropy taken out of I and
    // sigma (RemoveAnisotropy). Fills the resolution range read.
    std::vector<MergedReflection> StrongReflections(const std::vector<MergedReflection>& p1_merged,
                                                    const gemmi::UnitCell& cell, const gemmi::GroupOps& gops,
                                                    TwinImmuneZoneResult& result) {
        std::vector<const MergedReflection*> usable;
        for (const auto& r : p1_merged) {
            const gemmi::Op::Miller h{{r.h, r.k, r.l}};
            if (UsableIntensity(r) && std::isfinite(r.sigma) && r.sigma > 0.0
                && !gops.is_systematically_absent(h) && gops.epsilon_factor_without_centering(h) == 1)
                usable.push_back(&r);
        }
        // The shells with signal enough to read (ZONE_MIN_SHELL_ISIG), 20 equal in 1/d^2.
        constexpr int n_shells = 20;
        double min_s = std::numeric_limits<double>::infinity(), max_s = 0.0;
        for (const auto* r : usable) {
            min_s = std::min(min_s, 1.0 / (r->d * r->d));
            max_s = std::max(max_s, 1.0 / (r->d * r->d));
        }
        std::vector<MergedReflection> strong;
        if (!(max_s > min_s))
            return strong;
        auto shell_of = [&](const MergedReflection* r) {
            const double t = (1.0 / (r->d * r->d) - min_s) / (max_s - min_s);
            return std::min(n_shells - 1, std::max(0, static_cast<int>(t * n_shells)));
        };
        std::vector<double> isig_sum(n_shells, 0.0);
        std::vector<int> isig_n(n_shells, 0);
        for (const auto* r : usable) {
            isig_sum[shell_of(r)] += r->I / r->sigma;
            isig_n[shell_of(r)] += 1;
        }
        std::vector<int> shell;
        for (const auto* r : usable) {
            const int b = shell_of(r);
            if (isig_n[b] > 0 && isig_sum[b] / isig_n[b] >= ZONE_MIN_SHELL_ISIG) {
                strong.push_back(*r);
                shell.push_back(b);
            }
        }
        result.d_max_A = 0.0;
        result.d_min_A = std::numeric_limits<double>::infinity();
        for (const auto& r : strong) {
            result.d_max_A = std::max<double>(result.d_max_A, r.d);
            result.d_min_A = std::min<double>(result.d_min_A, r.d);
        }
        if (!strong.empty())
            RemoveAnisotropy(strong, shell, cell, gops, result);
        return strong;
    }

    // The reflections acentric in the group: the control every zone is read beside.
    TwinImmuneZone ReadControl(const std::vector<MergedReflection>& strong, const gemmi::GroupOps& gops,
                               const std::array<double, 3>* tncs_vector) {
        std::vector<const MergedReflection*> control;
        for (const auto& r : strong)
            if (!gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}}))
                control.push_back(&r);
        TwinImmuneZone z = ReadZone(control, tncs_vector);
        z.operators = "acentric";
        return z;
    }

    // The zone of one subgroup hypothesis H (given by its rotations): the reflections centric in the
    // group but not in H - the self-mates of the coset the group adds over H.
    TwinImmuneZone ReadSubgroupZone(const std::vector<MergedReflection>& strong,
                                    const gemmi::GroupOps& gops, const std::vector<gemmi::Op>& rotations,
                                    const std::vector<gemmi::Op>& h_ops,
                                    const std::array<double, 3>* tncs_vector) {
        auto contains = [&](const std::vector<gemmi::Op>& set, const gemmi::Op& op) {
            return std::any_of(set.begin(), set.end(), [&](const gemmi::Op& o) { return o.rot == op.rot; });
        };
        std::vector<const MergedReflection*> zone;
        for (const auto& r : strong) {
            const gemmi::Op::Miller h{{r.h, r.k, r.l}};
            const gemmi::Op::Miller minus{{-gemmi::Op::DEN * r.h, -gemmi::Op::DEN * r.k, -gemmi::Op::DEN * r.l}};
            const bool centric_in_h = std::any_of(h_ops.begin(), h_ops.end(), [&](const gemmi::Op& op) {
                return op.apply_to_hkl_without_division(h) == minus;
            });
            if (gops.is_reflection_centric(h) && !centric_in_h)
                zone.push_back(&r);
        }
        TwinImmuneZone z = ReadZone(zone, tncs_vector);
        for (const auto& op : rotations)
            if (!contains(h_ops, op))
                z.operators += (z.operators.empty() ? "" : " ") + op.as_hkl().triplet();
        return z;
    }
}

TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector<MergedReflection>& p1_merged,
                                           const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group,
                                           const gemmi::SpaceGroup& subgroup,
                                           const std::array<double, 3>* tncs_vector) {
    TwinImmuneZoneResult result;
    result.tncs_normalised = tncs_vector != nullptr;
    const gemmi::GroupOps gops = group.operations();
    const std::vector<gemmi::Op> rotations = ProperRotations(gops);
    const std::vector<gemmi::Op> h_ops = ProperRotations(subgroup.operations());
    const auto strong = StrongReflections(p1_merged, cell, gops, result);
    if (strong.empty())
        return result;
    result.control = ReadControl(strong, gops, tncs_vector);
    result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector));
    Calibrate(result);
    return result;
}

TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector<MergedReflection>& p1_merged,
                                            const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group,
                                            const std::array<double, 3>* tncs_vector) {
    TwinImmuneZoneResult result;
    result.tncs_normalised = tncs_vector != nullptr;
    const gemmi::GroupOps gops = group.operations();

    // The index-2 subgroups H of the group's rotations: each is a hypothesis "the crystal is H, and the
    // rest of the group (G minus H, one coset) is a twin law or a pseudo-symmetry". H holds every square
    // of G, so it is the subgroup the squares generate, or that plus one more rotation; for a proper
    // point group that finds them all (G over its squares is at most a four-group).
    auto same = [](const gemmi::Op& a, const gemmi::Op& b) { return a.rot == b.rot; };
    auto contains = [&](const std::vector<gemmi::Op>& set, const gemmi::Op& op) {
        return std::any_of(set.begin(), set.end(), [&](const gemmi::Op& o) { return same(o, op); });
    };
    auto closure = [&](const std::vector<gemmi::Op>& generators) {
        std::vector<gemmi::Op> ops;
        for (const auto& g : generators)
            if (!contains(ops, g))
                ops.push_back(g);
        for (size_t i = 0; i < ops.size(); ++i)
            for (size_t j = 0; j <= i; ++j)
                for (const gemmi::Op& p : {ops[i] * ops[j], ops[j] * ops[i]})
                    if (!contains(ops, p))
                        ops.push_back(p);
        return ops;
    };
    const std::vector<gemmi::Op> rotations = ProperRotations(gops);
    std::vector<gemmi::Op> squares{gemmi::Op::identity()};
    for (const auto& g : rotations)
        squares.push_back(g * g);
    const auto square_group = closure(squares);
    std::vector<std::vector<gemmi::Op>> subgroups;
    auto add_subgroup = [&](const std::vector<gemmi::Op>& h) {
        if (2 * h.size() != rotations.size())
            return;
        for (const auto& known : subgroups)
            if (known.size() == h.size()
                && std::all_of(h.begin(), h.end(), [&](const gemmi::Op& o) { return contains(known, o); }))
                return;
        subgroups.push_back(h);
    };
    add_subgroup(square_group);
    for (const auto& g : rotations) {
        auto gen = square_group;
        gen.push_back(g);
        add_subgroup(closure(gen));
    }
    if (subgroups.empty())
        return result;

    const auto strong = StrongReflections(p1_merged, cell, gops, result);
    if (strong.empty())
        return result;
    result.control = ReadControl(strong, gops, tncs_vector);

    // Per hypothesis H, the zone: centric in the group, acentric in H - the reflections some added
    // operator sends to -h, i.e. the self-mates of that coset read as a twin law.
    for (const auto& h_ops : subgroups)
        result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector));
    Calibrate(result);
    return result;
}

std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result) {
    std::ostringstream os;
    os << std::fixed;
    os << "Twin-immune zones of the operators the adopted group adds over each index-2 subgroup (P1 merge, " << std::setprecision(2) << result.d_max_A
       << "-" << result.d_min_A << " A, shells with <I/sigma> >= " << std::setprecision(0)
       << ZONE_MIN_SHELL_ISIG << ", a deviatoric quadratic form of delta B " << std::setprecision(1)
       << result.anisotropy_delta_b_A2 << " A^2 fitted over these shells and taken out"
       << (result.tncs_normalised ? ", normalised per pseudo-translation class" : "")
       << ")\n"
       << "  Reflections centric in the group but not in the subgroup are their own twin mates: centric if\n"
       << "  the added operators are real, acentric if they are a twin law or a pseudo-symmetry. <|E^2-1|>\n"
       << "  reads 0.968 centric, 0.736 acentric; evidence is the centric/acentric log-likelihood ratio\n"
       << "  (> 0: the added operators are real). The acentric control reads -0.130 nats per reflection\n"
       << "  when the normalisation is right; what it reads above that is the normalisation's (" << std::showpos
       << std::setprecision(3) << result.control_excess_per_reflection << std::noshowpos
       << " here), every\n  zone carries the same per reflection, and the calibrated evidence has it taken off.\n";
    auto row = [&](const TwinImmuneZone& z) {
        os << "  " << std::left << std::setw(34) << z.operators << std::right << " n " << std::setw(6) << z.n
           << "  <|E^2-1|> " << std::setprecision(3) << z.mean_abs_e2_minus_1 << " +- " << z.standard_error
           << "  evidence " << std::showpos << std::setprecision(1) << z.evidence_nats << " nats, calibrated "
           << z.calibrated_evidence_nats << std::noshowpos << "\n";
    };
    for (const auto& z : result.zones)
        row(z);
    row(result.control);
    return os.str();
}