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

#include "SearchSpaceGroup.h"
#include "../../common/ParallelFor.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <limits>
#include <map>
#include <numeric>
#include <sstream>
#include <string>
#include <tuple>
#include <unordered_map>
#include <vector>

namespace {
    // A merged reflection, folded onto the +/- Friedel-equivalent it represents, used as a
    // hash key to match symmetry-related reflections.
    struct HKLKey {
        int h = 0, k = 0, l = 0;
        bool operator==(const HKLKey& o) const noexcept { return h == o.h && k == o.k && l == o.l; }
    };

    struct HKLKeyHash {
        size_t operator()(const HKLKey& key) const noexcept {
            auto mix = [](uint64_t x) {
                x ^= x >> 33; x *= 0xff51afd7ed558ccdULL;
                x ^= x >> 33; x *= 0xc4ceb9fe1a85ec53ULL;
                x ^= x >> 33; return x;
            };
            return static_cast<size_t>(mix(static_cast<uint64_t>(key.h)) ^
                                       (mix(static_cast<uint64_t>(key.k)) << 1) ^
                                       (mix(static_cast<uint64_t>(key.l)) << 2));
        }
    };

    HKLKey Canonicalize(int h, int k, int l, bool merge_friedel) {
        if (merge_friedel && std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l))
            return {-h, -k, -l};
        return {h, k, l};
    }

    double PearsonCC(const std::vector<double>& x, const std::vector<double>& y) {
        if (x.size() < 2)
            return std::numeric_limits<double>::quiet_NaN();

        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
        for (size_t i = 0; i < x.size(); ++i) {
            sx += x[i]; sy += y[i];
            sxx += x[i] * x[i]; syy += y[i] * y[i]; sxy += x[i] * y[i];
        }
        const double n = static_cast<double>(x.size());
        const double vx = sxx - sx * sx / n;
        const double vy = syy - sy * sy / n;
        if (vx <= 0 || vy <= 0)
            return std::numeric_limits<double>::quiet_NaN();
        return (sxy - sx * sy / n) / std::sqrt(vx * vy);
    }

    // A reflection is extinct from lattice centering alone (independent of any screw/glide) when a
    // centering translation makes its structure factor cancel. Mirrors the centering half of
    // gemmi::GroupOps::is_systematically_absent, so screw absences can be judged separately.
    bool CenteringAbsent(const gemmi::GroupOps& gops, const gemmi::Op::Miller& hkl) {
        for (size_t i = 1; i < gops.cen_ops.size(); ++i) {
            const auto& t = gops.cen_ops[i];
            if ((t[0] * hkl[0] + t[1] * hkl[1] + t[2] * hkl[2]) % gemmi::Op::DEN != 0)
                return true;
        }
        return false;
    }

    // The reciprocal-space ROW a reflection lies on: its direction, reduced by the gcd and
    // sign-canonicalised, so 0,0,l and 0,0,-l are one row and h,0,0 is a different one. Screw
    // absences are judged against the other reflections of their own row (see below).
    using AxialRow = std::array<int, 3>;

    AxialRow RowOf(int h, int k, int l) {
        const int g = std::gcd(std::gcd(std::abs(h), std::abs(k)), std::abs(l));
        if (g > 0) { h /= g; k /= g; l /= g; }
        if (std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l))
            return {-h, -k, -l};
        return {h, k, l};
    }

    // Median of an unordered set (reordered in place); 0 for an empty set.
    double MedianOf(std::vector<double>& v) {
        if (v.empty())
            return 0.0;
        const size_t mid = v.size() / 2;
        std::nth_element(v.begin(), v.begin() + mid, v.end());
        return v[mid];
    }

    // How unlikely a predicted-absent class would be if the condition producing it did not exist,
    // in nats. Used for a SCREW against the rest of its own axial row, and for a CENTERING against
    // the present class.
    //
    // Scoring the absence against its own control class follows the POINTLESS zone test
    // (Evans, Acta Cryst D67, 282-292 (2011), App. A3); the Beta tail here is an analytic null in
    // place of its control transforms.
    //
    // Under "no condition" the absent class and its control are both Wilson-distributed with the SAME
    // mean, so with each absent intensity expressed in units of the control mean, the fraction
    // T = sum_u / (sum_u + n_control) follows Beta(n_absent, n_control) exactly. The control's own
    // strength cancels out of it - which is the property a count does not have, and the reason a
    // uniformly weak class decides nothing here instead of deciding "absent". Returns -log of that
    // Beta lower tail.
    //
    // Only the leading term of the regularized incomplete beta is kept. It is exact as T -> 0, which is
    // where a condition is claimed, and dropping the (1-T)^n_control factor only ever UNDER-states the
    // evidence, which is the safe direction for a test that has to clear a bound.
    // No measurement can place a merged intensity at exactly zero, so sum_u is floored at a
    // thousandth of the control mean per absent reflection. Without it, a zone whose absences all
    // merged non-positive gives sum_u = 0 exactly, T clamps to the epsilon below, and EACH absent
    // reflection is worth ~690 nats. That was harmless while the value only had to clear a bound of
    // 20, but it is now summed across zones and ranks the candidates, and it inverted the ranking
    // outright: a zone of 2 absences that were never measurable scored 1378 nats where a genuine zone
    // of 6 absences at 1% of its row scores 22, so the candidate claiming a screw on an UNMEASURED
    // row beat the one whose rows are actually dead, by 60x.
    //
    // The constant is deliberately an order of magnitude below the precision any real merge reaches -
    // a thousandth of the row mean needs I/sigma ~ 1000 against that row, where ISa tops out near 40 -
    // so it can only ever remove the singularity, never suppress evidence a measurement could have
    // produced. Measured over the realistic range it changes no genuine zone at all (22.0, 34.7 and
    // 65.4 nats to four figures) and takes the unmeasurable ones to 13 and 36.
    constexpr double MIN_U_PER_ABSENT_REFLECTION = 1e-3;


    // The reciprocal-space row {1,0,0} written the way a crystallographer names the zone: h00.
    std::string RowLabel(std::array<int, 3> row) {
        static const char* letter[3] = {"h", "k", "l"};
        // RowOf canonicalises the sign lexicographically, which can leave the leading component
        // negative; hkl and -h-k-l are one row, so name it by the positive one.
        for (int v : row)
            if (v != 0) {
                if (v < 0)
                    for (int& c : row) c = -c;
                break;
            }
        std::string out;
        for (int i = 0; i < 3; ++i)
            out += row[i] == 0 ? "0" : (row[i] == 1 ? letter[i] : std::to_string(row[i]) + letter[i]);
        return out;
    }

    std::string FormatDouble(double v, int decimals) {
        std::ostringstream o;
        o << std::fixed << std::setprecision(decimals) << v;
        return o.str();
    }

    std::array<int, 9> RotKey(const gemmi::Op& op) {
        std::array<int, 9> out{};
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j)
                out[i * 3 + j] = op.rot[i][j];
        return out;
    }

    // Whether a CELL can host a group's rotations. A setting is a statement about DIRECTION -
    // P 1 1 2 puts the 2-fold on c and needs alpha = beta = 90, where P 1 2 1 puts it on b and needs
    // alpha = gamma = 90 - so a candidate offered in a setting the metric does not have would be
    // merged on axes the crystal does not have. Compared on the metric tensor, each element against
    // its own scale, which makes the bound a tolerance on a lattice ANGLE and on an axis-length ratio;
    // gemmi's own is_compatible_with_groupops takes one absolute eps in A^2 instead, which means
    // something different on a 30 A cell and on a 300 A one.
    constexpr double CELL_SETTING_TOLERANCE = 2e-3;   // ~0.11 deg on an angle, ~0.1% on an axis ratio

    bool CellHostsRotations(const gemmi::UnitCell& cell, const gemmi::GroupOps& gops) {
        const auto g = cell.metric_tensor();
        const double G[3][3] = {{g.u11, g.u12, g.u13}, {g.u12, g.u22, g.u23}, {g.u13, g.u23, g.u33}};
        for (const gemmi::Op& op : gops.sym_ops) {
            for (int i = 0; i < 3; ++i)
                for (int j = 0; j < 3; ++j) {
                    double v = 0;
                    for (int p = 0; p < 3; ++p)
                        for (int q = 0; q < 3; ++q)
                            v += static_cast<double>(op.rot[p][i]) * G[p][q] * op.rot[q][j];
                    v /= static_cast<double>(gemmi::Op::DEN) * gemmi::Op::DEN;
                    if (std::fabs(v - G[i][j]) > CELL_SETTING_TOLERANCE * std::sqrt(G[i][i] * G[j][j]))
                        return false;
                }
        }
        return true;
    }

    // How a candidate is NAMED in the report. short_name() is what the reference settings have
    // always been printed as, and it stays that; it cannot express a setting, though - P 1 2 1 and
    // P 1 1 2 are both "P2" - so a non-reference setting is printed as its extended Hermann-Mauguin
    // name, which is the only faithful one.
    std::string SettingName(const gemmi::SpaceGroup& sg) {
        return sg.is_reference_setting() ? sg.short_name() : sg.xhm();
    }

    // The absences a group predicts over a fixed probe of low-index reflections - everything Stage B
    // judges a candidate on. Two settings with the same signature are one hypothesis written twice
    // (an alternative-origin entry, or a screw whose absences the centering already accounts for), so
    // only the first of them is worth scoring. Four indices is enough for every Sohncke setting in
    // the table: taking the probe to six or eight admits exactly the same candidates.
    std::vector<char> AbsenceSignature(const gemmi::SpaceGroup& sg) {
        const gemmi::GroupOps gops = sg.operations();
        std::vector<char> out;
        out.reserve(9 * 9 * 9);
        for (int h = -4; h <= 4; ++h)
            for (int k = -4; k <= 4; ++k)
                for (int l = -4; l <= 4; ++l)
                    out.push_back(gops.is_systematically_absent({{h, k, l}}) ? 1 : 0);
        return out;
    }

    // The rotation part of a space group in the reference setting (identity included), as a
    // sorted list of matrices - the key that groups space groups into a candidate point group.
    // It must be the rotation SET, not gemmi's PointGroup enum: P321 and P312 are both "32" yet
    // have their 2-folds along different directions, and only the matrices tell them apart.
    using RotationSet = std::vector<std::array<int, 9>>;

    RotationSet RotationSetOf(const gemmi::SpaceGroup& sg) {
        RotationSet out;
        for (const auto& op : sg.operations().derive_symmorphic().sym_ops)
            out.push_back(RotKey(op));
        std::sort(out.begin(), out.end());
        return out;
    }

    // Proper rotations of a crystal system's holohedry (the highest lattice symmetry it can host),
    // in the reference setting. Any candidate point group must be a subgroup of this.
    RotationSet HolohedryRotationSet(gemmi::CrystalSystem system) {
        int number = 0;
        switch (system) {
            case gemmi::CrystalSystem::Triclinic:    number = 1;   break;  // P1
            case gemmi::CrystalSystem::Monoclinic:   number = 3;   break;  // P2  (unique axis b)
            case gemmi::CrystalSystem::Orthorhombic: number = 16;  break;  // P222
            case gemmi::CrystalSystem::Tetragonal:   number = 89;  break;  // P422
            case gemmi::CrystalSystem::Trigonal:     number = 155; break;  // R32
            case gemmi::CrystalSystem::Hexagonal:    number = 177; break;  // P622
            case gemmi::CrystalSystem::Cubic:        number = 207; break;  // P432
        }
        const auto* sg = gemmi::find_spacegroup_by_number(number);
        return sg ? RotationSetOf(*sg) : RotationSet{};
    }

    // A candidate point group: its proper rotations (reference setting) and a representative
    // symmorphic space group (used when only the point group is wanted, or for display).
    struct PointGroupInfo {
        RotationSet rotation_set;
        std::vector<gemmi::Op> rotations;       // non-identity proper rotations
        const gemmi::SpaceGroup* representative = nullptr;
        // This rung exists only because no reference setting names it (the a- and c-unique monoclinic
        // 2-folds). It may be ADOPTED like any other, but it must not take part in judging a higher
        // promotion - see the two sites that read it below.
        bool widened = false;
    };

    // Enumerate candidate point groups. When a holohedry is given (from the lattice metric), keep
    // only its subgroups - this both skips operators the lattice forbids and avoids accepting a
    // coincidental higher symmetry; all subgroups down to P1 are still candidates.
    std::vector<PointGroupInfo> EnumeratePointGroups(const std::optional<RotationSet>& holohedry,
                                                     const std::optional<gemmi::UnitCell>& cell,
                                                     bool all_rotation_sets) {
        std::vector<PointGroupInfo> out;
        std::map<RotationSet, size_t> index;

        // Two passes. The reference settings first, exactly as before; then, only if asked and only
        // for rotation sets NO reference setting carries, a setting that names them. The four such
        // sets are the a-unique and c-unique monoclinic 2-folds and the two rhombohedral-axes trigonal
        // groups, and without them a crystal whose only 2-fold lies on a or c has no rung to stand on
        // between P1 and 222, so it falls to P1. Restricting the second pass to sets the first did not
        // reach makes it a pure addition: every point group reachable before is still reached, by the
        // same group, in the same setting.
        for (int pass = 0; pass < 2; ++pass) {
            if (pass == 1 && (!all_rotation_sets || !cell.has_value()))
                break;
            for (const auto& sg : gemmi::spacegroup_tables::main) {
                if (!sg.is_sohncke() || sg.is_reference_setting() != (pass == 0))
                    continue;

                RotationSet rs = RotationSetOf(sg);
                if (holohedry.has_value() &&
                    !std::includes(holohedry->begin(), holohedry->end(), rs.begin(), rs.end()))
                    continue;
                if (pass == 1 &&
                    (index.count(rs) > 0 || !CellHostsRotations(*cell, sg.operations())))
                    continue;

                auto it = index.find(rs);
                size_t pos;
                if (it == index.end()) {
                    PointGroupInfo info;
                    for (const auto& op : sg.operations().derive_symmorphic().sym_ops) {
                        if (op.rot == gemmi::Op::identity().rot)
                            continue;
                        info.rotations.push_back(gemmi::Op{op.rot, {0, 0, 0}, op.notation});
                    }
                    info.rotation_set = rs;
                    info.widened = (pass == 1);
                    pos = out.size();
                    index[rs] = pos;
                    out.push_back(std::move(info));
                } else {
                    pos = it->second;
                }
                // Prefer a symmorphic representative (the plain point-group setting).
                auto& info = out[pos];
                if (info.representative == nullptr ||
                    (!info.representative->is_symmorphic() && sg.is_symmorphic()))
                    info.representative = &sg;
            }
        }
        return out;
    }
}

double AbsenceEvidence(double sum_u, int n_absent, int n_control) {
    if (n_absent <= 0 || n_control <= 0)
        return 0.0;
    const double a = n_absent, b = n_control;
    const double u = std::max(sum_u, a * MIN_U_PER_ABSENT_REFLECTION);
    const double T = u / (u + b);
    return -(a * std::log(T) + std::lgamma(a + b) - std::lgamma(a + 1) - std::lgamma(b));
}

// The same for a SCREW zone, with the control COUNT taken out of it - the b -> infinity limit of
// AbsenceEvidence, i.e. -log P(Gamma(n_absent, 1) <= sum_u) to leading order.
//
// A screw's control class is the COMPLEMENT of its absent class on one axial row, so the two move
// together: a candidate that predicts more of the row absent leaves fewer reflections to be judged
// against. AbsenceEvidence grows with the control count, so that candidate is charged for the very
// reflections it correctly called extinct, and a group whose absent class is a strict SUPERSET of
// another's, with the extra reflections equally dead, could score LOWER. Measured on a tetragonal
// 4_1/4_3 crystal: 29 dead 00l against 8 control read 47.7 nats where a subset of 19 of them against
// 18 control read 55.5 - the wrong order, from the control count alone.
//
// The count belongs in a p-value for ONE candidate and not in a ranking of several: each candidate's
// tail is computed against its own null, and -log p from different nulls is not one scale. What is
// left is the likelihood ratio of the absent class - dead against Wilson at the row's own mean -
// which is a sum over reflections and so comparable across candidates. Asymptotically it is
// n_absent * (log(1/ubar) - 1), linear in the number of absences at fixed deadness, so an equally
// dead superset can no longer score lower. It is also the same number wherever the control class
// dwarfs the absent one, and it keeps the property the whole test is built on: a class as strong as
// its control reads about -n_absent, so a uniformly weak axial row still decides nothing.
//
// NOT used for the CENTERING class (AbsenceEvidence above), whose control is the whole present
// population rather than the complement of the claim on one row.
double ScrewZoneEvidence(double sum_u, int n_absent) {
    if (n_absent <= 0)
        return 0.0;
    const double a = n_absent;
    const double u = std::max(sum_u, a * MIN_U_PER_ABSENT_REFLECTION);
    return -a * std::log(u) + std::lgamma(a + 1);
}

SearchSpaceGroupResult SearchSpaceGroup(
    const std::vector<MergedReflection>& merged,
    const SearchSpaceGroupOptions& opt) {

    SearchSpaceGroupResult result;
    if (merged.empty())
        return result;

    const size_t n = merged.size();

    // Flatten the reflections and mark which ones each stage may use. The correlation stage drops
    // weak reflections; the absence stage must keep them - that is where the screw-axis signal is.
    std::vector<int> H(n), K(n), L(n);
    std::vector<double> I(n), Sigma(n), IoverSigma(n);
    std::vector<HKLKey> key(n);
    std::vector<char> pass_absence(n, 0), pass_cc(n, 0);

    for (size_t i = 0; i < n; ++i) {
        const auto& r = merged[i];
        H[i] = r.h; K[i] = r.k; L[i] = r.l;
        I[i] = r.I;
        Sigma[i] = std::isfinite(r.sigma) && r.sigma > 0 ? r.sigma : 0.0;
        key[i] = Canonicalize(r.h, r.k, r.l, opt.merge_friedel);

        const bool finite = std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0 &&
                            std::isfinite(r.d) && r.d > 0;
        const bool in_range = finite && (opt.d_min_limit_A <= 0 || r.d >= opt.d_min_limit_A);
        IoverSigma[i] = finite ? r.I / r.sigma : 0.0;
        pass_absence[i] = in_range;
    }

    // present_i_over_sigma is a cut on the reflection's own significance, and the merged I/sigma is not
    // that: it carries the error model's (b*I)^2 term, so it saturates at ISa = 1/b for a reflection
    // measured once and stops rising with the intensity above that knee. Convert the cut to the
    // quantity the merge exports, once, here - see SearchSpaceGroupOptions::merge_isa for the
    // derivation and for what it is worth. An unknown ISa (or b = 0, no systematic term) leaves the cut
    // exactly where the caller set it.
    double present_cut = opt.present_i_over_sigma;
    if (opt.merge_isa > 0.0) {
        const double r = opt.present_i_over_sigma / opt.merge_isa;
        present_cut = opt.present_i_over_sigma / std::sqrt(1.0 + r * r);
    }

    // The correlation stage uses only genuinely-present reflections. Near-zero (systematically
    // absent) reflections would otherwise form a second cluster at the origin and fake a high
    // correlation for false operators - fatal on centered lattices, where half the reflections
    // are extinct.
    for (size_t i = 0; i < n; ++i)
        pass_cc[i] = pass_absence[i] && IoverSigma[i] >= present_cut &&
                     (opt.min_i_over_sigma <= 0 || IoverSigma[i] >= opt.min_i_over_sigma);

    // Resolution-normalised intensity E^2 = I / <I>(shell), from equal-count resolution shells over a
    // given subset of the merge. Which subset matters: E^2 is only free of the resolution fall-off on
    // the population it was normalised over, so a caller has to normalise over the reflections it will
    // actually use (see Ecc below).
    auto shell_normalised = [&](const std::vector<char>& subset) {
        std::vector<double> E(n, 0.0);
        std::vector<size_t> order;
        order.reserve(n);
        for (size_t i = 0; i < n; ++i)
            if (subset[i])
                order.push_back(i);
        std::sort(order.begin(), order.end(),
                  [&](size_t a, size_t b) { return merged[a].d > merged[b].d; });   // low res -> high res
        // A hundred reflections a shell, at most twenty-five shells - and integer division, so
        // below two hundred there is a single shell and E^2 is I over one global mean. That
        // leaves the operator correlation exactly where it was, Pearson being invariant to a
        // common scale, but not the E^2 overlap cap below: it then cuts on I against the whole
        // merge instead of against each reflection's own shell.
        const int bins = std::clamp(static_cast<int>(order.size() / 100), 1, 25);
        const size_t per = (order.size() + bins - 1) / std::max(1, bins);
        for (size_t b = 0; b * per < order.size(); ++b) {
            const size_t lo = b * per, hi = std::min(order.size(), lo + per);
            double sum = 0.0;
            for (size_t j = lo; j < hi; ++j)
                sum += I[order[j]];
            const double mean = (hi > lo) ? sum / static_cast<double>(hi - lo) : 0.0;
            for (size_t j = lo; j < hi; ++j)
                E[order[j]] = mean > 0.0 ? I[order[j]] / mean : 0.0;
        }
        return E;
    };

    // Over the reflections the absence test uses, so it can judge "present" by intensity magnitude
    // rather than by a possibly under-estimated sigma (see present_e_squared).
    const std::vector<double> Esq = shell_normalised(pass_absence);

    // Overlap guard (Stage A / correlation only): drop the extreme resolution-normalised-E tail, which
    // on a two-lattice crystal is the one-sided overlap contamination that poisons the operator CC.
    // See SearchSpaceGroupOptions::max_e_squared_for_cc. Absences (pass_absence) keep the full range.
    if (opt.max_e_squared_for_cc > 0.0)
        for (size_t i = 0; i < n; ++i)
            if (pass_cc[i] && Esq[i] > opt.max_e_squared_for_cc)
                pass_cc[i] = false;

    // The operator correlation is scored on E^2, not on raw I. Both members of a symmetry pair sit at
    // the same |s|, so the resolution fall-off is variance shared perfectly between the two arms of
    // every pair: a Pearson CC on raw I measures the fall-off as well as the symmetry, and the fall-off
    // lifts a FALSE operator's CC as much as a true one's. Measured over the rotation battery with
    // shell-matched random pairing as the null for a metrically-allowed false operator, that raw-I noise
    // floor has a median of 0.31 and reaches 0.53 on one crystal - above the old bound of 0.5 outright -
    // and it varies more between crystals (spread 0.46) than the whole true/false gap is wide (0.38), so
    // an absolute bound on it was a different test on every crystal. It also moves with the search
    // resolution cut, by a median 0.09 and up to 0.23, which is what made that cut decide symmetries.
    // Normalised, that floor has a median of 0.015 and a maximum of 0.06, and moves by a median 0.03.
    // Following POINTLESS, which likewise scores each symmetry element on normalised intensities
    // (Evans, Acta Cryst. D62, 72-82 (2006)).
    //
    // Normalised over pass_cc - the reflections the correlation actually pairs - and NOT over the
    // pass_absence set Esq uses. pass_cc keeps only the stronger reflections and the fraction it keeps
    // itself falls with resolution, so an E^2 built on pass_absence still carries a resolution trend
    // inside the pass_cc subset; reusing Esq here makes the coupling WORSE than raw I (floor movement
    // 0.14). Esq stays as it is: the absence tests need their own set, and the E^2 cap above helps
    // DEFINE pass_cc, so normalising that over pass_cc would be circular.
    //
    // The price is a new dependence in place of the old one: the CC now moves with whatever defines
    // pass_cc, i.e. with present_i_over_sigma and the merge_isa conversion above. That conversion is
    // what keeps the dependence harmless - it holds the cut at one counting significance on every
    // crystal, so the population this is normalised over means the same thing on all of them.
    const std::vector<double> Ecc = shell_normalised(pass_cc);

    std::unordered_map<HKLKey, int, HKLKeyHash> key_to_index;
    key_to_index.reserve(n * 2);
    for (size_t i = 0; i < n; ++i)
        if (pass_absence[i])
            key_to_index.emplace(key[i], static_cast<int>(i));

    // --- Stage A: score each distinct rotation operator once ---
    // `visited` and `epoch` are scratch, taken as arguments rather than captured so that several
    // operators can be scored at once - each worker below keeps its own pair.
    auto score_operator = [&](const gemmi::Op& op, std::vector<uint32_t>& visited,
                              uint32_t& epoch) -> SpaceGroupOperatorScore {
        ++epoch;
        std::vector<double> x, y;      // raw merged I of each pair, for the H statistic
        std::vector<double> ex, ey;    // the same pairs as E^2, for the correlation
        for (size_t i = 0; i < n; ++i) {
            if (!pass_cc[i] || visited[i] == epoch)
                continue;
            const auto m2 = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}});
            const HKLKey k2 = Canonicalize(m2[0], m2[1], m2[2], opt.merge_friedel);
            if (k2 == key[i])
                continue;                       // reflection lies on this rotation axis
            const auto it = key_to_index.find(k2);
            if (it == key_to_index.end())
                continue;
            const int j = it->second;
            if (!pass_cc[j])
                continue;
            x.push_back(I[i]);
            y.push_back(I[j]);
            ex.push_back(Ecc[i]);
            ey.push_back(Ecc[j]);
            visited[i] = epoch;
            visited[j] = epoch;
        }
        SpaceGroupOperatorScore s;
        s.op_triplet_hkl = op.as_hkl().triplet('h');
        s.n_pairs = static_cast<int>(x.size());
        s.cc = PearsonCC(ex, ey);
        // Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio).
        // On RAW I, deliberately - which is why the pairs are collected twice. The shell divisor cancels
        // in |I1-I2|/(I1+I2) exactly for a pair whose members share a shell, but not for one that
        // straddles a shell boundary: measured over the battery's operators, normalising moves H by a
        // median 0.08% but by 3.7% at p95 and 12% at worst, and max_operator_h_ratio has been decided on
        // a margin of 1.8%. H is calibrated on raw I and stays there.
        std::vector<double> hv;
        hv.reserve(x.size());
        for (size_t p = 0; p < x.size(); ++p) {
            const double denom = x[p] + y[p];
            if (denom > 0.0)
                hv.push_back(std::fabs(x[p] - y[p]) / denom);
        }
        if (!hv.empty()) {
            const size_t mid = hv.size() / 2;
            std::nth_element(hv.begin(), hv.begin() + mid, hv.end());
            s.h_stat = hv[mid];
        }
        s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) &&
                    s.cc >= opt.min_operator_cc;
        return s;
    };

    std::map<std::array<int, 9>, SpaceGroupOperatorScore> op_cache;
    std::vector<uint32_t> visited(n, 0);
    uint32_t epoch = 0;
    auto operator_score = [&](const gemmi::Op& op) -> const SpaceGroupOperatorScore& {
        const auto rk = RotKey(op);
        auto it = op_cache.find(rk);
        if (it != op_cache.end())
            return it->second;
        return op_cache.emplace(rk, score_operator(op, visited, epoch)).first->second;
    };

    // Conjugate rotations (symmetry-equivalent within the point group) relate symmetry-equivalent
    // reflection sets, so on real data their CCs cluster; a noisy crystal can push one class member
    // below min_operator_cc while the class is unmistakably present (this was first seen on a cubic
    // crystal whose three 3-folds spread over 0.13 of CC with only the weakest below the bound).
    // Judge each conjugacy class by its mean CC, not its weakest
    // member, so a genuine high-symmetry point group is not lost to one marginal operator. chi2_under
    // (below) remains the safety net against a truly false promotion. Returns {all classes present,
    // worst class-mean CC}.
    auto point_group_present = [&](const std::vector<gemmi::Op>& rots) -> std::pair<bool, double> {
        const size_t m = rots.size();
        std::vector<int> cls(m, -1);
        int n_cls = 0;
        for (size_t i = 0; i < m; ++i) {
            if (cls[i] >= 0)
                continue;
            cls[i] = n_cls;
            for (size_t j = i + 1; j < m; ++j)
                if (cls[j] < 0)
                    for (const auto& p : rots)
                        if ((p * rots[i] * p.inverse()).rot == rots[j].rot) {
                            cls[j] = n_cls;
                            break;
                        }
            ++n_cls;
        }
        bool ok = true;
        double worst_mean = 1.0;
        for (int c = 0; c < n_cls; ++c) {
            // Average over the class members that actually have enough pairs to score; a single
            // low-multiplicity / degenerate (NaN) operator in an otherwise strong class is skipped,
            // not allowed to veto the class. The class must still have at least one scored member.
            double sum_cc = 0.0;
            int n_valid = 0;
            for (size_t i = 0; i < m; ++i)
                if (cls[i] == c) {
                    const auto& s = operator_score(rots[i]);
                    if (s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc)) {
                        sum_cc += s.cc;
                        ++n_valid;
                    }
                }
            const double mean_cc = n_valid > 0 ? sum_cc / n_valid : 0.0;
            worst_mean = std::min(worst_mean, mean_cc);
            if (n_valid == 0 || mean_cc < opt.min_operator_cc)
                ok = false;
        }
        return {ok, worst_mean};
    };

    std::optional<RotationSet> holohedry;
    if (opt.lattice_system.has_value())
        holohedry = HolohedryRotationSet(opt.lattice_system.value());
    const auto point_groups =
        EnumeratePointGroups(holohedry, opt.cell, opt.enumerate_all_rotation_sets);

    // Every operator the search can ask about comes from this list, and scoring one is a pass over the
    // whole merge with a hash lookup per reflection - the most expensive thing in here. They do not
    // depend on each other, so score the distinct ones now and let the search below read the cache.
    // One `visited` per worker, not per operator: it is as long as the merge, so allocating it per
    // operator would cost more than the scoring.
    {
        std::vector<gemmi::Op> distinct;
        std::vector<std::array<int, 9>> keys;
        for (const auto& pg : point_groups)
            for (const auto& rot : pg.rotations) {
                const auto rk = RotKey(rot);
                if (std::find(keys.begin(), keys.end(), rk) == keys.end()) {
                    keys.push_back(rk);
                    distinct.push_back(rot);
                }
            }
        std::vector<SpaceGroupOperatorScore> scored(distinct.size());
        ParallelChunks(static_cast<int>(distinct.size()),
                       std::min(opt.nthreads, distinct.size()), [&](int lo, int hi) {
            std::vector<uint32_t> scratch(n, 0);
            uint32_t ep = 0;
            for (int i = lo; i < hi; ++i)
                scored[i] = score_operator(distinct[i], scratch, ep);
        });
        for (size_t i = 0; i < distinct.size(); ++i)
            op_cache.emplace(keys[i], scored[i]);
    }

    // Mapping every observation onto its symmetry representative under a candidate's rotations - one
    // apply_to_hkl + Canonicalize per observation per operator - is the expensive half of BOTH
    // quantities below, and both need exactly the same mapping. Build it once per point group.
    struct Acc { double sw = 0.0, swI = 0.0; int n = 0; };
    struct Orbits {
        std::vector<Acc> acc;          // one inverse-variance accumulator per orbit, in first-seen order
        std::vector<int32_t> orbit;    // observation -> its orbit, -1 for one neither quantity below uses
    };
    auto build_orbits = [&](const std::vector<gemmi::Op>& rotations) -> Orbits {
        Orbits orb;
        orb.orbit.assign(n, -1);
        // The representative is interned to a dense index right here, so the two quantities below
        // index an array rather than hashing a key per observation - each of them is a pass over the
        // whole merge, and the lookup was the larger half of both. Same orbits, same order, same sums.
        std::unordered_map<HKLKey, int32_t, HKLKeyHash> orbit_id;
        for (size_t i = 0; i < n; ++i) {
            if (!pass_cc[i] || !(Sigma[i] > 0.0))
                continue;
            HKLKey best = key[i];
            for (const auto& op : rotations) {
                const auto m = op.apply_to_hkl(gemmi::Op::Miller{{H[i], K[i], L[i]}});
                const HKLKey k2 = Canonicalize(m[0], m[1], m[2], opt.merge_friedel);
                if (std::make_tuple(k2.h, k2.k, k2.l) < std::make_tuple(best.h, best.k, best.l))
                    best = k2;
            }
            const auto [it, fresh] = orbit_id.emplace(best, static_cast<int32_t>(orb.acc.size()));
            if (fresh)
                orb.acc.emplace_back();
            orb.orbit[i] = it->second;
            Acc& g = orb.acc[it->second];
            const double w = 1.0 / (Sigma[i] * Sigma[i]);
            g.sw += w; g.swI += w * I[i]; g.n += 1;
        }
        return orb;
    };

    // Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry
    // equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces
    // non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is
    // more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well
    // - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over
    // the present (pass_cc) reflections.
    auto chi2_under = [&](const Orbits& orb) -> double {
        double chi2 = 0.0;
        long dof = 0;
        for (size_t i = 0; i < n; ++i) {
            if (orb.orbit[i] < 0)
                continue;
            const Acc& g = orb.acc[orb.orbit[i]];
            if (g.n < 2)
                continue;
            const double mean = g.swI / g.sw, dev = I[i] - mean;
            chi2 += dev * dev / (Sigma[i] * Sigma[i]);
        }
        for (const Acc& g : orb.acc)
            if (g.n >= 2)
                dof += g.n - 1;
        return dof > 0 ? chi2 / static_cast<double>(dof) : std::numeric_limits<double>::quiet_NaN();
    };

    // Extra intensity-proportional systematic error a point group's merge has to invoke to reconcile
    // its symmetry equivalents: the smallest b for which sigma^2 + (b I)^2 brings the merged reduced
    // chi^2 down to 1. A genuine symmetry needs almost none - its equivalents already agree within
    // their errors, so the extra scatter is random and multiplicity absorbs it. A twin or pseudo-
    // symmetry forces non-equivalent reflections together, and that disagreement scales with I, so b
    // has to grow to swallow it (mirroring the merge error model's b / ISa collapse). This isolates
    // the systematic part of the scatter, which the fixed-sigma chi^2 ratio cannot: a genuine but
    // imperfectly-scaled high-symmetry merge and a twin can share a chi^2 ratio (~2) yet differ
    // sharply here (a genuine cubic step b x1.04 vs a merohedral twin b x1.6).
    // This `b` is a fraction of I fitted with the sigma^2 coefficient held at 1 - it is NOT the merge
    // error model's b, and NOT XDS's b either. The gate constants below are calibrated in this
    // convention; converting them to any other silently squares the ratios and makes the absolute
    // floor a-dependent, on a quantity that has no a. Leave it alone.
    auto merge_systematic_b = [&](const Orbits& orb) -> double {
        std::vector<std::array<double, 3>> obs;   // I, sigma, deviation-from-orbit-mean
        for (size_t i = 0; i < n; ++i) {
            if (orb.orbit[i] < 0)
                continue;
            const Acc& g = orb.acc[orb.orbit[i]];
            if (g.n < 2)
                continue;
            obs.push_back({I[i], Sigma[i], I[i] - g.swI / g.sw});
        }
        if (obs.size() < 20)
            return 0.0;
        auto reduced_chi2 = [&](double b) {
            double s = 0.0;
            for (const auto& o : obs)
                s += o[2] * o[2] / (o[1] * o[1] + (b * o[0]) * (b * o[0]));
            return s / static_cast<double>(obs.size());
        };
        if (reduced_chi2(0.0) <= 1.0)
            return 0.0;
        double lo = 0.0, hi = 2.0;   // b is a fraction of I; 2.0 = 200% is far past any real error model
        for (int it = 0; it < 40; ++it) {
            const double mid = 0.5 * (lo + hi);
            (reduced_chi2(mid) > 1.0 ? lo : hi) = mid;
        }
        return 0.5 * (lo + hi);
    };

    // Operator-CC-confirmed candidates, each with its merge chi^2 and systematic-error b; chi2_ref =
    // the most consistent.
    struct PGCand { const PointGroupInfo* pg; int order; double min_class_cc; double chi2; double b_extra;
                    // Filled by the selection loop below and carried so the adopted candidate's H
                    // ratio can be reported whether or not the bound had anything to say about it.
                    double h_ratio = std::numeric_limits<double>::quiet_NaN(); };
    int refused_order = 0;
    std::string refused_pg_hm, refused_why;
    std::vector<PGCand> pg_cands;
    double chi2_ref = std::numeric_limits<double>::infinity();
    // Which point groups their own operators confirm. Serial: the operators were all scored in
    // parallel above, so this is a handful of cache lookups per group.
    std::vector<const PointGroupInfo*> confirmed;
    std::vector<double> confirmed_cc;
    for (const auto& pg : point_groups) {
        const auto [present, min_class_cc] = point_group_present(pg.rotations);
        if (!present)
            continue;
        confirmed.push_back(&pg);
        confirmed_cc.push_back(min_class_cc);
    }
    // Merging under a candidate costs three passes over the whole merge, and there are a handful of
    // candidates, each reading nothing but the shared reflection arrays - so give each one a thread.
    // Every candidate's sums are still formed on one thread, over the same reflections in the same
    // order, so the numbers this hands back do not depend on the split.
    std::vector<double> cand_chi2(confirmed.size(), std::numeric_limits<double>::quiet_NaN());
    std::vector<double> cand_b(confirmed.size(), 0.0);
    ParallelFor(static_cast<int>(confirmed.size()),
                std::min<size_t>(opt.nthreads, confirmed.size()), [&](int i) {
        if (confirmed[i]->rotations.empty())
            return;
        const Orbits orb = build_orbits(confirmed[i]->rotations);
        cand_chi2[i] = chi2_under(orb);
        cand_b[i] = merge_systematic_b(orb);
    });
    for (size_t i = 0; i < confirmed.size(); ++i) {
        const PointGroupInfo& pg = *confirmed[i];
        pg_cands.push_back({&pg, static_cast<int>(pg.rotations.size()) + 1, confirmed_cc[i],
                            cand_chi2[i], cand_b[i]});
        // A rung the widened enumeration added does not set the reference chi^2. Otherwise adding it
        // would make every HIGHER promotion harder to reach - chi2_ref is a minimum - and the search
        // could answer LOWER because a candidate was offered, which is the opposite of what offering
        // it is for. Measured before this line existed: an F-cubic crystal's 432 and an orthorhombic
        // crystal's 222 were both refused once the a- and c-unique 2-folds joined the list.
        if (!pg.rotations.empty() && std::isfinite(cand_chi2[i]) && !pg.widened)
            chi2_ref = std::min(chi2_ref, cand_chi2[i]);
    }

    // Choose the largest point group that is both operator-confirmed AND self-consistent (its merge
    // chi^2 is not inflated past the miscalibration-widened bound below; ties -> higher min class CC).
    // Identity (no operators) is always consistent, so it stays the P1 fallback.
    const PointGroupInfo* best_pg = nullptr;
    int best_pg_order = 0;
    double best_pg_min_cc = -2.0;
    for (auto& c : pg_cands) {
        // A genuine symmetry operator merges equivalent reflections, so it barely changes the reduced
        // chi^2 relative to the best subgroup - across the whole rotation-test battery every correct
        // point group stays within ~1.7x, even on weak or badly-integrated data (a cubic F432 chi2_ref
        // 8.3 -> 1.15; a tetragonal P41212 -> 1.71). A twin law or pseudo-symmetry forces non-equivalent
        // reflections together, so its ratio is markedly higher (a merohedral twin 2-fold: R3 3.02 ->
        // R32 6.07, ratio 2.01). max_merge_chi2_ratio sits between the two. (An earlier log10(chi2_ref) widening
        // compensated for an under-calibrated error model that inflated real-symmetry ratios with data
        // weakness; the variance-floor fix removed that inflation, and the widening now only let the
        // twin through, so it is gone.)
        bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) ||
                          !std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio;

        // Systematic-error test vs the largest confirmed subgroup (by rotation-set inclusion): merging
        // under a genuine operator gains multiplicity without intensity-proportional disagreement, so the
        // merge error model's b barely moves; a merohedral twin forces non-equivalent reflections together
        // and b balloons. It both RESCUES a genuine step whose chi^2 drifts just past the ratio bound
        // (imperfectly scaled data) and VETOES a twin whose chi^2 now looks self-consistent but whose b
        // balloons - the chi^2 ratio alone no longer separates them.
        double parent_b = -1.0;
        // Every confirmed subgroup of the largest order below this candidate. There can be more than one
        // - 422 has both 4 and 222 - and on a twinned crystal the rival is not a harmless alternative: a
        // P4 crystal twinned by 2[100] has its two twin 2-folds confirmed, so 222 is CC-confirmed too and
        // CONTAINS the twin laws. Normalising the H test against it hides the twin among the promotion's
        // own real operators (measured on the synthetic grid: ratio 8.19 against the true parent 4, 0.78
        // against the rival 222). Which one is the true parent is exactly what is unknown here, so the
        // promotion must answer to all of them.
        std::vector<const PointGroupInfo *> parents;
        if (!c.pg->rotations.empty()) {
            int parent_order = 0;
            for (const auto& s : pg_cands)
                // Same reason as chi2_ref above: a rung only the widened enumeration offers is not a
                // parent. The promotion answers to the most damning of its parents, so admitting two
                // more order-2 subgroups of 222 makes the 222 step strictly harder than it was before
                // the rung was offered at all.
                if (!s.pg->widened && s.order < c.order && s.order >= parent_order
                    && std::includes(c.pg->rotation_set.begin(), c.pg->rotation_set.end(),
                                     s.pg->rotation_set.begin(), s.pg->rotation_set.end())) {
                    if (s.order > parent_order) {
                        parent_order = s.order;
                        parent_b = s.b_extra;
                        parents.clear();
                    } else {
                        // Tied parent. The b tests answer to the most damning of them, exactly as the H
                        // test below does: the SMALLEST parent b is the one that makes the veto easiest to
                        // trip and the rescue hardest to pass, and a rival subgroup that already contains
                        // the twin laws has its own b ballooned - taking it would hide the twin.
                        parent_b = std::min(parent_b, s.b_extra);
                    }
                    parents.push_back(s.pg);
                }
        }

        // Sigma-free twin test: compare the disagreement H of the operators this promotion ADDS with the
        // disagreement of the parent's own operators, measured on the same reflections. A real operator
        // relates equal intensities and matches the parent; a twin law relates different ones and reads
        // systematically higher. Skipped when either side has too few pairs to mean anything, and when
        // there is no parent group to normalise against (the first step out of P1). Where several parents
        // tie (see above), the promotion is judged on the most damning of them.
        double h_ratio = std::numeric_limits<double>::quiet_NaN();
        double h_added = 0.0;                 // the added operators' own H, for the twin fraction below
        for (const auto *parent : parents) {
            double h_new = 0.0, h_par = 0.0;
            int n_new = 0, n_par = 0, pairs_new = 0, pairs_par = 0;
            for (const auto &rot : c.pg->rotations) {
                if (rot.rot == gemmi::Op::identity().rot)
                    continue;
                const auto &os = operator_score(rot);
                if (os.n_pairs < opt.min_pairs_per_operator)
                    continue;
                const bool in_parent = std::binary_search(parent->rotation_set.begin(),
                                                          parent->rotation_set.end(), RotKey(rot));
                if (in_parent) { h_par += os.h_stat; ++n_par; pairs_par += os.n_pairs; }
                else           { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; }
            }
            if (n_new > 0 && n_par > 0 && pairs_new >= opt.min_pairs_for_h
                && pairs_par >= opt.min_pairs_for_h && h_par > 0.0) {
                const double r = (h_new / n_new) / (h_par / n_par);
                if (!std::isfinite(h_ratio) || r > h_ratio) {
                    h_ratio = r;
                    h_added = h_new / n_new;
                }
            }
        }
        c.h_ratio = h_ratio;
        // The chi^2 ratio is only trustworthy when the error model is calibrated. When even the best
        // subgroup's reduced chi^2 (chi2_ref) is far above 1 - weak, low-resolution data whose merged
        // sigmas are badly under-estimated - the ratio grows with point-group order for genuine high
        // symmetry too and wrongly rejects it (a true weak F432 reaches ratio ~14). The systematic-b test
        // re-fits its own error, so it stays valid under a broken sigma model: a genuine step's b barely
        // moves (b-ratio ~1) while a twin's balloons. So once chi2_ref shows the error model is unreliable,
        // a promotion is rescued on the b-test alone (subject to the balloon veto below); otherwise the
        // rescue is confined to the narrow chi^2 band just past the ratio bound.
        const bool miscalibrated = std::isfinite(chi2_ref) && chi2_ref > opt.chi2_ref_reliable;
        if (!consistent && parent_b > 1e-4 && c.b_extra <= parent_b * opt.max_systematic_b_ratio
            && (miscalibrated || (std::isfinite(c.chi2) && std::isfinite(chi2_ref)
                                  && c.chi2 <= chi2_ref * opt.max_merge_chi2_rescue)))
            consistent = true;
        // Veto a chi^2-passing promotion whose b clearly ballooned (above the largest genuine step, below a
        // twin); a genuine but imperfectly-scaled high-symmetry merge stays under the bound and is untouched.
        // The parent b is floored (min_systematic_b_for_veto) so a near-zero parent on excellent data cannot
        // fabricate a huge ratio out of a still-tiny absolute b (a genuine 422 at b=0.05 over a 222 parent at
        // b=0.008 is not a twin - a real twin drives b to ~0.19 regardless).
        //
        // The veto also answers to the H test: it does not fire on a promotion H has CONFIRMED. The two
        // can disagree because they measure different things - b is fitted against the merged sigmas, so
        // it moves with the measurement, while H compares intensities with intensities - and where they
        // do, H is the one that was measured to separate a twin from a real symmetry (see
        // max_operator_h_ratio and max_systematic_b_veto). A twin fails both.
        const bool h_confirmed = std::isfinite(h_ratio) && h_ratio <= opt.max_operator_h_ratio;
        bool b_refused = false;
        if (consistent && !h_confirmed && parent_b > 1e-4
            && c.b_extra > std::max(parent_b, opt.min_systematic_b_for_veto) * opt.max_systematic_b_veto) {
            consistent = false;
            b_refused = true;
        }

        // The H test is a necessary condition for promotion where it can be computed: it is the only
        // statistic measured to separate genuine symmetry from a merohedral twin across data amounts.
        const bool h_refused = std::isfinite(h_ratio) && h_ratio > opt.max_operator_h_ratio;
        if (h_refused)
            consistent = false;

        if (!consistent) {
            // Record the highest-order refusal so the caller can say WHY it is processing lower.
            if (c.order > refused_order && c.pg->representative) {
                refused_order = c.order;
                refused_pg_hm = c.pg->representative->point_group_hm();
                if (h_refused)
                    // Twinning by a fraction a scales every twin-related difference by (1-2a), and
                    // |I1-I2|/(I1+I2) is uniform on [0,1] for untwinned Wilson intensities, so the
                    // added operator's median H reads (1-2a)/2 and implies a = 0.5 - H. Quoted because
                    // it is the number a user acts on. It is a LOWER bound: measurement error only
                    // adds to H, and so only subtracts from a.
                    refused_why = "operator disagreement H is " + FormatDouble(h_ratio, 2)
                                + "x the parent's (bound " + FormatDouble(opt.max_operator_h_ratio, 2)
                                + ") - the added operator relates unequal intensities, as a twin law of "
                                  "fraction " + FormatDouble(std::max(0.0, 0.5 - h_added), 2)
                                + " or more would";
                else if (b_refused)
                    // Against the FLOORED parent b, which is what the test compared - quoting the raw
                    // parent b reads as a far bigger balloon than the one that actually tripped the veto.
                    refused_why = "its merge's systematic error b is " + FormatDouble(c.b_extra, 3)
                                + " against the subgroup's "
                                + FormatDouble(std::max(parent_b, opt.min_systematic_b_for_veto), 3)
                                + " (bound " + FormatDouble(opt.max_systematic_b_veto, 2)
                                + "x) - the added operator forces unequal intensities together";
                else if (std::isfinite(c.chi2) && std::isfinite(chi2_ref))
                    refused_why = "merge chi^2 is " + FormatDouble(c.chi2 / chi2_ref, 2)
                                + "x the subgroup's (bound " + FormatDouble(opt.max_merge_chi2_ratio, 2) + ")";
                else
                    refused_why = "the merge under it is not self-consistent";
            }
            continue;
        }
        if (c.order > best_pg_order || (c.order == best_pg_order && c.min_class_cc > best_pg_min_cc)) {
            best_pg = c.pg;
            best_pg_order = c.order;
            best_pg_min_cc = c.min_class_cc;
        }
    }

    for (const auto& [rk, s] : op_cache)
        result.operator_scores.push_back(s);
    std::sort(result.operator_scores.begin(), result.operator_scores.end(),
              [](const auto& a, const auto& b) { return a.cc > b.cc; });

    // A caller that already decided the point group elsewhere overrides the choice here, keeping the
    // operator scores and the refusal report Stage A just produced. Nothing else is bypassed: Stage B
    // below judges the absences on THIS merge's reflections, which is the whole point of pinning.
    if (opt.fixed_point_group.has_value()) {
        const RotationSet want = RotationSetOf(*opt.fixed_point_group);
        for (const auto& pg : point_groups)
            if (pg.rotation_set == want) {
                best_pg = &pg;
                best_pg_order = static_cast<int>(pg.rotations.size()) + 1;
                break;
            }
    }

    if (best_pg == nullptr)                          // should not happen (C1 always qualifies)
        return result;

    if (best_pg->representative) {
        result.point_group_hm = best_pg->representative->point_group_hm();
        result.point_group_representative = *best_pg->representative;
    }
    result.point_group_order = best_pg_order;

    // The H ratio of the promotion that was ADOPTED, reported whether the bound had anything to say
    // about it or not. Read after the choice is final, so a fixed_point_group override reports the
    // ratio of the group it forced rather than of the one Stage A would have taken.
    for (const auto& c : pg_cands)
        if (c.pg == best_pg)
            result.h_ratio = c.h_ratio;
    result.h_ratio_bound = opt.max_operator_h_ratio;

    // Only report a refusal that is actually ABOVE what was adopted.
    if (refused_order > best_pg_order) {
        result.refused_point_group_hm = refused_pg_hm;
        result.refused_reason = refused_why;
    }

    // --- Stage B: pick the space group within the point group ---
    // Without screw/centering determination, return the symmorphic representative.
    if (!opt.determine_space_group || best_pg->rotations.empty()) {
        if (best_pg->representative)
            result.best_space_group = *best_pg->representative;
        return result;
    }

    // The candidate space groups of the chosen point group. Scoring one is a pass over the whole
    // merge with three absence tests per reflection, and there are up to a dozen of them; they read
    // nothing but the shared reflection arrays, so give each one a thread. Every candidate's own pass
    // is unchanged and they are appended in table order, so the ranking below sees what it saw before.
    std::vector<const gemmi::SpaceGroup*> sg_cands;
    for (const auto& sg : gemmi::spacegroup_tables::main)
        if (sg.is_sohncke() && sg.is_reference_setting() && RotationSetOf(sg) == best_pg->rotation_set)
            sg_cands.push_back(&sg);
    // The non-reference SETTINGS of the same point group, when asked for - and unconditionally when
    // the point group is one only a non-reference setting carries (Stage A's second pass), since
    // otherwise there is no candidate at all and the group would be lost after being found. These are
    // alternative namings at the same order, so nothing here can promote the point group; what they
    // add is a screw or a centering on the axis the data show it on. Two refusals bound them: a
    // setting whose axes the cell does not have is not offered, and one predicting exactly the
    // absences a candidate already offered predicts is the same hypothesis under another name.
    if (opt.cell.has_value() &&
        (opt.enumerate_all_settings || (opt.enumerate_all_rotation_sets && sg_cands.empty()))) {
        std::vector<std::vector<char>> signatures;
        for (const auto* c : sg_cands)
            signatures.push_back(AbsenceSignature(*c));
        for (const auto& sg : gemmi::spacegroup_tables::main) {
            if (!sg.is_sohncke() || sg.is_reference_setting() ||
                RotationSetOf(sg) != best_pg->rotation_set ||
                !CellHostsRotations(*opt.cell, sg.operations()))
                continue;
            auto sig = AbsenceSignature(sg);
            if (std::find(signatures.begin(), signatures.end(), sig) != signatures.end())
                continue;
            signatures.push_back(std::move(sig));
            sg_cands.push_back(&sg);
        }
    }
    std::vector<SpaceGroupCandidateScore> sg_scored(sg_cands.size());
    ParallelFor(static_cast<int>(sg_cands.size()),
                std::min<size_t>(opt.nthreads, sg_cands.size()), [&](int ci) {
        const gemmi::SpaceGroup& sg = *sg_cands[ci];
        const gemmi::GroupOps gops = sg.operations();
        SpaceGroupCandidateScore s{.space_group = sg};
        double absent_sum = 0, present_sum = 0;
        int present_n = 0;

        // Judge centering and screw/glide absences on separate reflection sets. Lumping them lets
        // a large, correct centering-absent set hide a few strong screw violations and over-claim
        // screw axes (e.g. I4_132 on I432 data).
        int centering_absent = 0, centering_violations = 0;
        double centering_absent_sum = 0;
        // The same two classes as E^2, which is what AbsenceEvidence is stated in: the centring-absent
        // class and the present class that is its control.
        double centering_absent_esq = 0, present_esq_sum = 0;
        int present_strong = 0;
        // A screw axis extinguishes only the reflections that lie ON it, so its absent class and the
        // rest of that same axial row are collected apart from the general reflections and judged
        // against each other, ROW BY ROW (see screw_e_squared below).
        struct ScrewAbsent { AxialRow row; double e_squared; double i_over_sigma; };
        std::vector<ScrewAbsent> screw_absent_refl;
        std::map<AxialRow, std::vector<double>> row_present_esq;

        for (size_t i = 0; i < n; ++i) {
            if (!pass_absence[i])
                continue;
            const gemmi::Op::Miller hkl{{H[i], K[i], L[i]}};
            // Present := statistically significant AND intensity-significant. The E^2 gate keeps a
            // weak axial reflection with an under-estimated sigma (fake high I/sigma) from faking a
            // screw-axis violation; it only relaxes "present", so it cannot over-call a screw whose
            // predicted-absent class carries real intensity.
            //
            // present_cut, not the fixed cut: on a merge weak enough that nothing clears the fixed cut,
            // screw_violations is identically zero, so every screw axis passes unchallenged, and
            // present_strong is zero, so the centering rescue below switches itself off on exactly the
            // weak data it exists for. present_cut is the same cut converted to the counting scale the
            // two stages share (see merge_isa); on a healthy merge it is the fixed cut to within 1%.
            const bool present = IoverSigma[i] > present_cut &&
                                 (opt.present_e_squared <= 0.0 || Esq[i] > opt.present_e_squared);

            if (CenteringAbsent(gops, hkl)) {
                s.absent_observed += 1;
                absent_sum += IoverSigma[i];
                centering_absent += 1;
                centering_absent_sum += IoverSigma[i];
                centering_absent_esq += std::max(0.0, Esq[i]);
                if (present) { s.absent_violations += 1; centering_violations += 1; }
            } else if (gops.is_systematically_absent(hkl)) {
                s.absent_observed += 1;
                absent_sum += IoverSigma[i];
                screw_absent_refl.push_back({RowOf(H[i], K[i], L[i]), Esq[i], IoverSigma[i]});
            } else {
                present_n += 1;
                present_sum += IoverSigma[i];
                present_esq_sum += std::max(0.0, Esq[i]);
                if (present) present_strong += 1;
                // A non-identity rotation of the group maps this reflection to itself, i.e. it lies on
                // a rotation axis - the control class for that axis's screw absences.
                if (gops.epsilon_factor_without_centering(hkl) > 1)
                    row_present_esq[RowOf(H[i], K[i], L[i])].push_back(Esq[i]);
            }
        }

        // "Too strong to be systematically absent" is judged, for a SCREW, against the axial row the
        // screw constrains rather than against the shell mean over all reflections. An axial row can be
        // far stronger than an average reflection, and <I>(shell) falls with resolution while a
        // systematically-absent reflection keeps a small non-decaying residual (background / profile
        // leakage) - so at high resolution the plain E^2 cut turns those residuals into violations even
        // though the reflections next to them in the same row are tens of times stronger. That cost a
        // real P4_1 2_1 2 crystal its 4_1: 18 of its 47 absent 00l crossed the cut, all beyond 3.7 A,
        // at 1-2% of the l=4n reflections beside them. Scaling by the row's own median E^2 removes the
        // resolution dependence; the scale is floored at 1 so a row weaker than average keeps the plain
        // cut, which makes this a rescue only - a screw can be recovered by it, never lost.
        //
        // Row by row, not pooled: a 4_1 along c and a 2_1 along a are separate conditions with separate
        // control rows, and on the same crystal one row can be 20x an average reflection while another
        // is half of one. Pooling them lets the weak rows set the threshold for the strong one and the
        // rescue never fires (that very crystal pooled to a row median of 0.7 and stayed at P42_12).
        // A row needs more than a reflection or two behind its median before that median may set the
        // threshold. The scale only ever RAISES the cut, so a row whose control class holds one strong
        // reflection would license a screw claim the row does not support - a genuine 4_2 whose 00l are
        // observed only at l=4n, strongly, would read its l=4n+2 reflections as absent and be reported as
        // 4_1/4_3. Below this count the row falls back to the plain cut, i.e. no rescue.
        constexpr size_t MIN_ROW_CONTROL_REFLECTIONS = 3;
        // That count guards the row MEDIAN, which sets a violation THRESHOLD. The row MEAN is a
        // different thing - it is the scale the zone's evidence is stated in, and a scale needs a
        // number to divide by, not a middle. Holding both to the same count made a candidate FORFEIT
        // a zone by being right: a screw predicts more of its own row absent, which leaves fewer
        // reflections in its control, and below three the zone was declared undetermined and its
        // evidence discarded. Measured on a tetragonal 4_1/4_3 wedge: thirteen 00l reflections at
        // 0.1% of the two l = 4n beside them scored ZERO, while the nine of them a 4_2 also predicts
        // absent, judged against a control that still held four dead l = 4n+2, scored 38.7 - so the
        // group that explains the row lost to one that explains part of it.
        constexpr size_t MIN_ROW_CONTROL_FOR_SCALE = 2;
        std::map<AxialRow, double> row_median, row_mean;
        for (auto& [row, esq] : row_present_esq) {
            if (esq.size() >= MIN_ROW_CONTROL_FOR_SCALE)
                row_mean[row] = std::accumulate(esq.begin(), esq.end(), 0.0) / esq.size();
            if (esq.size() >= MIN_ROW_CONTROL_REFLECTIONS)
                row_median[row] = MedianOf(esq);
        }

        const int screw_absent = static_cast<int>(screw_absent_refl.size());
        int screw_violations = 0;
        // The absent intensities in units of their own row's control mean, kept ZONE BY ZONE - one
        // axial row is one screw condition, and the rows are measured to wildly different depths. A
        // row with no control class of its own cannot say whether its absences are weak or its whole
        // row is, so it contributes no evidence either way rather than being judged against the shell
        // mean; it is carried as an undetermined zone so the report can say so.
        struct ZoneSums { double sum_u = 0; int n_absent = 0; int n_control = 0; };
        std::map<AxialRow, ZoneSums> zones;
        for (const auto& a : screw_absent_refl) {
            const auto it = row_median.find(a.row);
            const double row_scale = it == row_median.end() ? 1.0 : std::max(1.0, it->second);
            // present_cut for the same reason as the general "present" above - with the fixed cut this
            // test cannot fire at all on a low-ISa merge. The row-relative E^2 gate is the other half of
            // the test and is unchanged, so an absence still has to look strong against its OWN axial row
            // to count as a violation.
            if (a.i_over_sigma > present_cut &&
                (opt.present_e_squared <= 0.0 || a.e_squared > opt.present_e_squared * row_scale))
                ++screw_violations;
            auto& z = zones[a.row];
            z.n_absent += 1;
            const auto mit = row_mean.find(a.row);
            if (mit != row_mean.end() && mit->second > 0.0) {
                z.sum_u += std::max(0.0, a.e_squared) / mit->second;
                z.n_control = static_cast<int>(row_present_esq.at(a.row).size());
            }
        }
        // Each zone is scored on its OWN evidence, and the strongest one speaks for the group. Pooling
        // the rows into one verdict lets a row the sweep barely sampled outweigh one that is decisively
        // dead - measured: a crystal whose h zone is confirmed on its own row and whose k zone is noise
        // pooled to a refusal and came out in the symmorphic group, with the confirmed screw lost and
        // nothing said about it. The pooled value is kept as a floor so several jointly-decisive zones
        // still count together, which is the case pooling gets right.
        double screw_sum_u = 0.0;
        int screw_scored = 0;
        double best_zone_evidence = 0.0;
        for (const auto& [row, z] : zones) {
            const double ev = z.n_control > 0 ? ScrewZoneEvidence(z.sum_u, z.n_absent) : 0.0;
            s.screw_zones.push_back({row, z.n_absent, z.n_control, ev});
            if (z.n_control > 0) {
                screw_sum_u += z.sum_u;
                screw_scored += z.n_absent;
                best_zone_evidence = std::max(best_zone_evidence, ev);
            }
        }
        s.absent_violations += screw_violations;
        s.centering_absent = centering_absent;
        // A centred group whose centering-absent class is EMPTY was not tested, it was skipped: the
        // reflections that would decide it are not in this merge (see centering_untested).
        s.centering_untested = sg.centring_type() != 'P' && centering_absent == 0;
        s.screw_absent = screw_absent;
        s.screw_absence_evidence =
            std::max(best_zone_evidence, ScrewZoneEvidence(screw_sum_u, screw_scored));

        if (s.absent_observed > 0)
            s.absent_mean_i_over_sigma = absent_sum / s.absent_observed;
        if (present_n > 0)
            s.present_mean_i_over_sigma = present_sum / present_n;

        // The centering-absent class gets the same likelihood a screw zone gets, with the present class
        // as its control: how unlikely that class would be if the centering did not exist. A COUNT of
        // net absences cannot separate a true centering from a super-centering that swallows it, because
        // the super-centering predicts every true absence PLUS a block of present reflections and can
        // still net the larger count - measured, an F222 candidate with a third of its "absent" class as
        // strong as the present one out-counted the C222_1 whose class was at 0.4% of it. Under the
        // likelihood those two are orders of magnitude apart, on the same numbers, with no bound to tune.
        const double present_esq_mean = present_n > 0 ? present_esq_sum / present_n : 0.0;
        if (centering_absent > 0 && present_esq_mean > 0.0)
            s.centering_absence_evidence =
                AbsenceEvidence(centering_absent_esq / present_esq_mean, centering_absent, present_n);

        // Centering is judged by class STRENGTH, not a per-reflection violation count. A real centering
        // cancels structure factors, so its absent class is systematically weak - its mean signed
        // I/sigma sits well below the present class - regardless of noise or obverse/reverse twinning;
        // a false centering leaves the "absent" class as strong as the present one (mean ratio ~1). The
        // count-of-strong-violations gate is brittle on noisy/twinned data, where enough genuinely-absent
        // reflections randomly clear I/sigma>3 to trip the 10% bound though the class is 3-4x weaker (a
        // true R3 at 13.5% violations, absent 1.7 vs present 6.0). The mean is well-determined here
        // because a centering-absent class holds a third-to-half of all reflections. Screws keep the
        // count gate: their predicted-absent class is a handful of axial reflections, too few to average
        // - what they get instead is a row-relative threshold for counting a violation at all.
        const double present_mean = present_n > 0 ? present_sum / present_n : 0.0;
        const double centering_absent_mean =
            centering_absent > 0 ? centering_absent_sum / centering_absent : 0.0;
        // The centering-absent class proves itself weak in either of two floor-independent ways; a
        // FALSE centering (absent as strong as present) fails both:
        //  (1) mean signed I/sigma well below the present class, OR
        //  (2) its strong-reflection RATE well below the present class's own strong rate.
        // (2) is needed because weak / low-energy data carry a positive intensity floor (background /
        // profile leakage) that lifts <I/s>abs to ~1.5-2.3 even for genuinely extinct reflections; when
        // the present class is itself weak (small present_mean) that additive floor inflates the mean
        // ratio past the bound and hides a real centering - e.g. an I-centred cubic crystal at low
        // energy, whose true I-centering sat at ratio ~0.57. Normalising the violation count by the
        // present class's own strong rate cancels the shared floor and stays reliable on weak data
        // (both rates shrink together).
        const double present_strong_rate =
            present_n > 0 ? static_cast<double>(present_strong) / present_n : 0.0;
        const double centering_violation_rate =
            centering_absent > 0 ? static_cast<double>(centering_violations) / centering_absent : 0.0;
        // The mean-ratio test only means anything while the present class carries signal: with a
        // present mean at or below zero the bound is non-positive and the comparison turns on the sign
        // of the absent mean rather than on its size, accepting or rejecting a centering by accident.
        // Leave that case to the rate test below, which counts violations and cannot change sign.
        const bool centering_ok = centering_absent == 0
            || (present_n > 0 && present_mean > 0.0
                && centering_absent_mean <= opt.max_absent_present_ratio * present_mean)
            || (present_strong_rate > 0.0
                && centering_violation_rate <= opt.max_absent_present_ratio * present_strong_rate);
        const bool screw_ok = screw_absent == 0 ||
            screw_violations <= opt.max_absent_violation_fraction * screw_absent;
        s.consistent = centering_ok && screw_ok;
        sg_scored[ci] = std::move(s);
    });
    for (auto& s : sg_scored)
        result.candidates.push_back(std::move(s));

    // A candidate is eligible when its absences are confirmed and there are enough of them to
    // trust (the symmorphic group, with no absences, is always eligible as the fallback). Rank
    // eligible candidates by the EVIDENCE their predicted absences carry, centering and screws
    // together - both are -log Beta tails of an absent class against its own control, so they are in
    // the same units (nats) and add as independent evidence does. A count, even netted against the
    // violations, is the wrong scale: a false super-centering predicts every true absence PLUS a block
    // of present reflections, so it can net the larger count while a third of its class is as strong as
    // the present one. The likelihood collapses for exactly that class and grows for a genuine one, so
    // a real F222 still beats C222_1 on F-centred data and no bound has to be tuned to tell them apart.
    // The screw part is the SUM over the zones, not the group's single gating number: each axial row
    // is a separate condition tested on its own reflections, so the log-likelihoods add. Summing is
    // what makes an extra condition pay its own way - a group claiming one screw more gains that
    // zone's evidence when the row really is dead and gains little when the row is not, where a count
    // could only ever go up. (It gains LITTLE, not nothing: a zone whose absences were never
    // measurable still scores a bounded positive value, see MIN_U_PER_ABSENT_REFLECTION. Requiring a
    // minimum number of absences per zone before it may contribute, as the centring class does
    // through min_absent_observed, would close that too - it needs its own battery.) Pooling the rows into one Beta instead would make three genuine screws
    // read as weaker evidence than two whenever the third row is measured less deeply, which is a
    // property of the pooling, not of the crystal.
    auto absence_evidence = [](const SpaceGroupCandidateScore& s) {
        double total = s.centering_absence_evidence;
        for (const auto& z : s.screw_zones)
            if (z.n_control > 0)
                total += z.evidence;
        return total;
    };
    // A candidate also needs enough EVIDENCE behind its absences before they may be claimed, and the two
    // kinds of absence need different measures of it. A centering class is a third to a half of every
    // reflection in the data set, so a count is a fair measure and min_absent_observed is never the
    // binding constraint. A screw class is a handful of axial reflections BY CONSTRUCTION - one row of
    // reciprocal space, often lying near the spindle where a rotation sweep records least - and there a
    // count measures the sweep's geometry, not the evidence: six axial reflections measured at zero
    // against a row that averages 1.4x the shell mean settle the question, while twenty uniformly weak
    // ones settle nothing. So the screw class is judged by its zones' AbsenceEvidence instead.
    // A candidate in a non-reference setting must also have had its CENTERING tested here. The
    // reference-setting path can adopt an untested centering because something else backs it - the
    // caller's centred-lattice re-test reindexes the data into the metric candidate's conventional
    // cell and only then commits - and a non-reference setting has no such backing: adopting one on
    // an absence class this merge does not contain would name a lattice the data never showed, which
    // is worse than the group it displaces. Refused, not warned about.
    auto eligible = [&](const SpaceGroupCandidateScore& s) {
        return s.consistent
            && !(s.centering_untested && !s.space_group.is_reference_setting())
            && (s.centering_absent == 0 || s.centering_absent >= opt.min_absent_observed)
            && (s.screw_absent == 0 || s.screw_absence_evidence >= opt.min_screw_absence_evidence);
    };
    std::sort(result.candidates.begin(), result.candidates.end(),
              [&](const SpaceGroupCandidateScore& a, const SpaceGroupCandidateScore& b) {
                  if (eligible(a) != eligible(b))
                      return eligible(a);
                  if (absence_evidence(a) != absence_evidence(b))
                      return absence_evidence(a) > absence_evidence(b);
                  if (a.absent_violations != b.absent_violations)
                      return a.absent_violations < b.absent_violations;   // prefer the honest, less over-claiming group
                  // Genuinely indistinguishable (e.g. I23 vs I2_13, or an enantiomorphic pair): lower
                  // space-group number is the representative.
                  return a.space_group.number < b.space_group.number;
              });

    if (!result.candidates.empty() && eligible(result.candidates.front())) {
        // Alternatives are only the candidates with the SAME absence signature - identical absent AND
        // violation counts - as the winner: the enantiomorphic / origin-ambiguous partners the data
        // truly cannot separate. A super-centering that nets the same count but over-claims differs in
        // its violation count and is therefore not reported as an equal alternative.
        const int sel_absent = result.candidates.front().absent_observed;
        const int sel_violations = result.candidates.front().absent_violations;
        for (auto& s : result.candidates) {
            if (!eligible(s) || s.absent_observed != sel_absent || s.absent_violations != sel_violations)
                continue;
            s.selected = true;
            if (!result.best_space_group.has_value())
                result.best_space_group = s.space_group;   // representative (lowest number)
            else
                result.alternatives.push_back(s.space_group);
        }
    }

    return result;
}

std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result,
                                         size_t max_candidates_to_print) {
    std::ostringstream os;

    if (!result.refused_point_group_hm.empty())
        os << "Higher symmetry " << result.refused_point_group_hm << " was confirmed by the operator "
              "correlations but REFUSED: " << result.refused_reason << ".\n"
              "  Processing in the lower symmetry, which is the recoverable direction - if this is a "
              "twin, merging in the higher group would average non-equivalent reflections together and "
              "hide the twin law.\n";
    os << "Point group: " << (result.point_group_hm.empty() ? "?" : result.point_group_hm)
       << " (from intensity correlations)\n";

    os << "  " << std::setw(14) << std::left << "operator" << std::right
       << std::setw(9) << "CC" << std::setw(10) << "pairs" << std::setw(9) << "symm"
       << std::setw(9) << "H" << "\n";
    for (const auto& s : result.operator_scores) {
        os << "  " << std::setw(14) << std::left << s.op_triplet_hkl << std::right
           << std::setw(9) << std::fixed << std::setprecision(3) << s.cc
           << std::setw(10) << s.n_pairs
           << std::setw(9) << (s.present ? "yes" : "no")
           << std::setw(9) << std::fixed << std::setprecision(3) << s.h_stat << "\n";
    }
    os << "  H = median |I1-I2|/(I1+I2) over the operator's pairs - the disagreement it implies, with\n"
          "  no sigma in it. The promotion gate is the RATIO of the mean H over the operators a\n"
          "  promotion adds to the mean over the parent group's own, which is what separates a real\n"
          "  symmetry (ratio near 1) from a merohedral twin law.\n";
    if (std::isfinite(result.h_ratio))
        os << "  H ratio " << FormatDouble(result.h_ratio, 2) << " for the adopted point group (bound "
           << FormatDouble(result.h_ratio_bound, 2) << ").\n";
    else
        os << "  H ratio not available (no parent group to normalise against, or too few pairs).\n";

    os << "\nSpace-group candidates\n";
    os << "  " << std::setw(10) << std::left << "SG" << std::right
       << std::setw(9) << "absent" << std::setw(7) << "viol"
       << std::setw(11) << "<I/s>abs" << std::setw(11) << "<I/s>pres"
       << std::setw(9) << "screw" << std::setw(11) << "screw evid"
       << std::setw(11) << "cent evid"
       << std::setw(6) << "OK" << std::setw(11) << "centering" << "\n";

    const size_t count = std::min(max_candidates_to_print, result.candidates.size());
    for (size_t i = 0; i < count; ++i) {
        const auto& c = result.candidates[i];
        os << (c.selected ? "* " : "  ")
           << std::setw(10) << std::left << SettingName(c.space_group) << std::right
           << std::setw(9) << c.absent_observed << std::setw(7) << c.absent_violations
           << std::setw(11) << std::fixed << std::setprecision(2) << c.absent_mean_i_over_sigma
           << std::setw(11) << std::fixed << std::setprecision(2) << c.present_mean_i_over_sigma
           << std::setw(9) << c.screw_absent
           << std::setw(11) << std::fixed << std::setprecision(1) << c.screw_absence_evidence
           << std::setw(11) << std::fixed << std::setprecision(1) << c.centering_absence_evidence
           << std::setw(6) << (c.consistent ? "yes" : "no")
           << std::setw(11) << (c.centering_untested ? "UNTESTED"
                                                     : (c.centering_absent > 0 ? "tested" : "-"))
           << "\n";
    }
    os << "  absent/viol = reflections the group predicts absent, and how many are nonetheless present.\n"
          "  screw = how many of those lie on an axial row, i.e. are extinguished by a screw rather than\n"
          "  by the centering; screw evid = how much likelier that class is if the screw exists than if it\n"
          "  does not, judged against the rest of its own axial row (in nats - a real screw reads tens to\n"
          "  hundreds, a false one at or below zero). cent evid is how unlikely the centering-absent class\n"
          "  would be without the centering, judged against the present class; the two together rank the\n"
          "  candidates. The <I/s> columns say little about screws, because the merged sigma shrinks with\n"
          "  I on absent and present alike.\n"
          "  centering = whether the group's centering was decided here at all. UNTESTED means this\n"
          "  merge holds none of the reflections that centering extinguishes - the data are indexed on\n"
          "  the primitive sub-cell, so those reflections were never predicted or integrated - and the\n"
          "  candidate scores zero absences for want of evidence, not because the centering is real.\n";

    // Per-zone screw verdicts for the group that was chosen. One axial row is one screw condition, and
    // the rows are measured to very different depths, so a row this merge cannot judge is reported as
    // undetermined rather than being folded into a single yes/no for the whole group.
    for (const auto& c : result.candidates) {
        if (!c.selected || c.screw_zones.empty())
            continue;
        os << "Screw conditions of " << SettingName(c.space_group) << ", zone by zone:\n";
        for (const auto& z : c.screw_zones) {
            os << "  " << RowLabel(z.row) << "  " << std::setw(3) << z.n_absent << " absent  ";
            if (z.n_control == 0)
                os << "UNDETERMINED - this row carries no control reflections here\n";
            else
                os << std::setw(4) << z.n_control << " control  " << std::fixed << std::setprecision(1)
                   << std::setw(8) << z.evidence << " nats\n";
        }
        break;
    }

    if (result.best_space_group.has_value()) {
        os << "Best space group: " << SettingName(*result.best_space_group);
        for (const auto& alt : result.alternatives)
            os << " or " << SettingName(alt);
        if (!result.alternatives.empty())
            os << " (indistinguishable from these data)";
        os << "\n";
        // A group with a different CENTERING is a different lattice, not just a different group, and
        // a run reports one cell - the chosen group's. "C2 or P21 or P2" printed beside a single
        // C-centred cell is not something a user can act on: P2 and P21 live on the primitive
        // sub-cell, with their own cell constants and their own Miller indices. Name the ones the
        // reported cell does not describe, and say how far off it is for each.
        const char sel_centring = result.best_space_group->centring_type();
        bool said_setting = false;
        for (const auto& alt : result.alternatives) {
            if (alt.centring_type() == sel_centring)
                continue;
            if (!said_setting) {
                os << "  These are NOT all the same lattice, and one cell is reported for this run -\n"
                      "  " << SettingName(*result.best_space_group) << "'s (centering " << sel_centring
                   << "). These need a cell of their own:\n";
                said_setting = true;
            }
            os << "    " << std::setw(10) << std::left << SettingName(alt) << std::right
               << "centering " << alt.centring_type() << ", cell volume "
               << std::fixed << std::setprecision(2)
               << static_cast<double>(gemmi::centring_vectors(alt.centring_type()).size())
                  / static_cast<double>(gemmi::centring_vectors(sel_centring).size())
               << "x the reported one\n";
        }
        if (said_setting)
            os << "  Adopting one of those means reindexing to its cell; the cell constants reported\n"
                  "  and the Miller indices written are the chosen group's alone.\n";
    } else {
        os << "Best space group: none determined\n";
    }

    return os.str();
}