Jungfraujoch/image_analysis/scale_merge/AnisotropyAnalysis.cpp

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

#include "AnisotropyAnalysis.h"

#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <numeric>
#include <tuple>
#include <sstream>

#include <Eigen/Dense>

#include "HKLKey.h"

#include "gemmi/eig3.hpp"
#include "gemmi/scaling.hpp"

namespace {

    using Tensor = gemmi::SMat33<double>;

    // ---------------------------------------------------------------- tuned constants
    // The binning is the one the shape test was validated on: 12 equal-count resolution shells x 60
    // directions on a hemisphere, a cell entering the fit once it holds 12 independent reflections.
    // Shell count, direction count and the cell-population floor were each varied by a factor of two
    // either way without moving a verdict, so none of the three is sharp.
    constexpr int N_SHELLS = 12;
    constexpr int N_DIRECTIONS = 60;
    constexpr int MIN_CELL_REFLECTIONS = 12;
    constexpr int MIN_CELLS = 40;               // fewer than this and no tensor is fitted at all
    constexpr int MIN_SHAPE_SHELLS = 5;         // fewer than this and the s^2 signature is not fitted
    constexpr int MIN_CELLS_PER_SHAPE_SHELL = 8;

    // The whole diagnostic is refused below this: at <I/sigma> ~ 1 neither the tensor nor its shape is
    // interpretable, and one battery dataset at 0.65 produced numbers that mean nothing.
    constexpr double MIN_MEAN_I_OVER_SIGMA = 1.0;

    // Below ~90 deg of OBSERVED rotation a lab-fixed systematic reaches a second and a third tensor
    // direction (canonical correlation 0.74-0.85 against 0.05 at a full sweep), the anisotropy such a
    // systematic manufactures grows 5-9x, and the power to detect a real 10 A^2 falls to <= 0.25.
    constexpr double MIN_ROTATION_DEG = 90.0;

    // Symmetry-expanded reflections are capped by striding the unique list: the estimator is a fit of a
    // few tens of cell means and does not improve past this, while the memory does grow.
    constexpr size_t MAX_ENTRIES = 4000000;

    // The unmerged arm strides whole unique reflections down to this many. The systematic-error floor is
    // set by the systematic, not by counting: over a 12x range in the number of unique reflections and a
    // 10x range in <I/sigma> the smallest anisotropy the gate can establish moves by less than 2x, so
    // more than this buys nothing and costs time and memory.
    constexpr int MAX_CLUSTERS = 20000;

    // The most partials ScaledObservations assembles on rotation data (see there). In partials because
    // partials are what its sort costs, and generous because how many unique reflections that is
    // depends on how many frames a rocking curve spans and how many mates the Laue class gives it,
    // which together vary by an order of magnitude: over the rotation battery it leaves between 25
    // thousand and 183 thousand of them, against the MAX_CLUSTERS the floor above actually reads.
    constexpr size_t MAX_SAMPLE_PARTS = 8u << 20;

    // splitmix64's finalizer. Sampling on this rather than on the key itself matters: a packed hkl key
    // has h in its low bits, so any direct test on it selects a set of lattice planes, not a sample.
    uint64_t MixKey(uint64_t x) {
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
        return x ^ (x >> 31);
    }

    // The cone half-angle of the directional diffraction limits, and the <I/sigma(I)> the limit is read
    // at. Both are AIMLESS's shipped defaults (ANALYSIS CONE 20, ISIGMINIMUM 1.5 -> here 2.0, the value
    // this project already uses for a resolution statement).
    constexpr double CONE_HALF_ANGLE_DEG = 20.0;
    constexpr double CONE_I_OVER_SIGMA = 2.0;

    // The s^2 signature. FLAT is declared at c0 > +0.70 with c0/sigma(c0) > 6: over a 38-crystal battery
    // that fires on six of the seven crystals whose deficit does not follow a Debye-Waller B, with no
    // false positive among the other 24, and the false-positive count stays 0 anywhere between z > 5 and
    // z > 8. A pure-B simulation on the same reflection lists returns |c0| <= 0.18, so the machinery
    // cannot invent an intercept of this size. The curvature arm uses the same z for symmetry; it has no
    // false-positive calibration of its own, and it is reported rather than acted on.
    constexpr double FLAT_INTERCEPT = 0.70;
    constexpr double SHAPE_Z = 6.0;
    // CONVEX asks for the s^4 term to carry more than half the deficit at the resolution limit, not
    // merely to be significant: a mild positive curvature is common and is not a different mechanism.
    constexpr double CONVEX_SHARE = 0.5;

    // Gate bands. The denominator is this dataset's own systematic error scale, not a counting-statistics
    // error bar: against known ground truth a counting-statistics gate calls "real" in 15-58% of clean
    // isotropic datasets, because the symmetry-forbidden directions are determined from contrasts BETWEEN
    // symmetry mates of one reflection (which share their true |F|^2, so Wilson scatter cancels) while
    // the symmetry-allowed ones can only be determined BETWEEN different reflections (which do not).
    // 3.5 is a 2-5% false-positive rate pooled over six artefact classes and 6-14% in the worst of them;
    // 5.0 is 0-0.5% and 0-3%.
    constexpr double GATE_MARGINAL = 2.0;
    constexpr double GATE_ESTABLISHED = 3.5;
    constexpr double GATE_STRONG = 5.0;

    // The forbidden subspace of a real dataset is dominated by systematics, not by counting noise, so
    // this ratio is above 1 on essentially every dataset: measured over the battery its median is 6.5
    // and only two of 22 exceed 20, so "it is above its counting null" would fire on almost all of
    // them. This is the upper tail. Together with p == 1 and a gate that established nothing it picks
    // out two datasets in the battery - one of them the crystal whose space-group determination is
    // itself known to be unstable - where the unconditional p == 1 caution fired on 17 of 38.
    constexpr double FORBIDDEN_Z_SUSPECT = 20.0;

    // E[dB] / sigma_per_component for a zero-mean Gaussian tensor confined to a subspace: the eigenvalues
    // of a noisily estimated tensor repel, so dB has a strictly positive expectation even when the tensor
    // is exactly zero. Per-subspace Monte-Carlo medians, indexed by the dimension of the subspace;
    // K_ALLOWED is for the symmetry-allowed subspaces, K_FORBIDDEN for the symmetry-forbidden ones (the
    // two differ at the same dimension because they are different subspaces of the deviatoric space).
    constexpr double K_ALLOWED[6] = {0.0, 0.824, 1.583, 2.075, 0.0, 2.872};
    constexpr double K_FORBIDDEN[6] = {0.0, 0.0, 1.66, 2.13, 2.53, 2.87};

    // ---------------------------------------------------------------- tensor algebra
    double FrobDot(const Tensor &a, const Tensor &b) {
        return a.u11 * b.u11 + a.u22 * b.u22 + a.u33 * b.u33
               + 2.0 * (a.u12 * b.u12 + a.u13 * b.u13 + a.u23 * b.u23);
    }

    Tensor Scale(const Tensor &a, double s) {
        return {a.u11 * s, a.u22 * s, a.u33 * s, a.u12 * s, a.u13 * s, a.u23 * s};
    }

    Tensor Deviatoric(const Tensor &a) { return a.added_kI(-a.trace() / 3.0); }

    // s^T E s from the cell-averaged outer product q = <(sx^2, sy^2, sz^2, sx sy, sx sz, sy sz)>.
    double QuadForm(const Tensor &e, const std::array<double, 6> &q) {
        return e.u11 * q[0] + e.u22 * q[1] + e.u33 * q[2]
               + 2.0 * (e.u12 * q[3] + e.u13 * q[4] + e.u23 * q[5]);
    }

    // Gram-Schmidt in the Frobenius inner product; directions that the earlier ones already span drop out.
    std::vector<Tensor> Orthonormalise(const std::vector<Tensor> &raw) {
        std::vector<Tensor> out;
        for (Tensor t : raw) {
            for (const auto &o : out)
                t = t - Scale(o, FrobDot(t, o));
            const double n = std::sqrt(FrobDot(t, t));
            if (n > 1e-9)
                out.push_back(Scale(t, 1.0 / n));
        }
        return out;
    }

    // The cell with its lengths and angles snapped to what the crystal system requires.
    gemmi::UnitCell MetricIdealCell(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) {
        if (!sg)
            return cell;
        double a = cell.a, b = cell.b, c = cell.c;
        double alpha = cell.alpha, beta = cell.beta, gamma = cell.gamma;
        switch (sg->crystal_system()) {
            case gemmi::CrystalSystem::Triclinic:
                break;
            case gemmi::CrystalSystem::Monoclinic:
                if (sg->monoclinic_unique_axis() == 'a') { beta = 90.0; gamma = 90.0; }
                else if (sg->monoclinic_unique_axis() == 'c') { alpha = 90.0; beta = 90.0; }
                else { alpha = 90.0; gamma = 90.0; }
                break;
            case gemmi::CrystalSystem::Orthorhombic:
                alpha = beta = gamma = 90.0;
                break;
            case gemmi::CrystalSystem::Tetragonal:
                a = b = 0.5 * (a + b);
                alpha = beta = gamma = 90.0;
                break;
            case gemmi::CrystalSystem::Trigonal:
            case gemmi::CrystalSystem::Hexagonal:
                if (sg->ext == 'R') {
                    a = b = c = (a + b + c) / 3.0;
                    alpha = beta = gamma = (alpha + beta + gamma) / 3.0;
                } else {
                    a = b = 0.5 * (a + b);
                    alpha = beta = 90.0;
                    gamma = 120.0;
                }
                break;
            case gemmi::CrystalSystem::Cubic:
                a = b = c = (a + b + c) / 3.0;
                alpha = beta = gamma = 90.0;
                break;
        }
        gemmi::UnitCell out;
        out.set(a, b, c, alpha, beta, gamma);
        return out;
    }

    // The symmetry-allowed anisotropy directions, deviatoric, in the Cartesian frame.
    //
    // gemmi::adp_symmetry_constraints gives the standard ADP symmetry constraints in the crystal-axis
    // (B*) parameterisation - 6 / 4 / 3 / 2 / 2 / 1 vectors for triclinic / monoclinic / orthorhombic /
    // tetragonal / trigonal-hexagonal / cubic. Transforming each to Cartesian is a bijection on symmetric
    // matrices, and the isotropic direction is always allowed, so removing the trace leaves exactly
    // 5 / 3 / 2 / 1 / 1 / 0 free deviatoric parameters.
    //
    // Following Sheriff & Hendrickson (1987) Acta Cryst. A43, 118-121
    std::vector<Tensor> AllowedDeviatoricBasis(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg) {
        std::vector<Tensor> raw;
        // The constraints are a statement about the IDEAL metric, and rugnux writes the unconstrained
        // refined cell: a beta of 89.97 deg in a C2 crystal leaves the isotropic tensor outside the span
        // of the monoclinic constraints, the trace is then not fully removed, and the count comes out
        // one too high. A cell whose angle is a hundredth of a degree off does not create a fourth
        // anisotropy direction, so the basis is built on the idealised metric.
        const gemmi::UnitCell ideal = MetricIdealCell(cell, sg);
        for (const gemmi::Vec6 &v : gemmi::adp_symmetry_constraints(sg)) {
            const Tensor b_star{v[0], v[1], v[2], v[3], v[4], v[5]};
            raw.push_back(Deviatoric(b_star.transformed_by<double>(ideal.orth.mat)));
        }
        return Orthonormalise(raw);
    }

    std::vector<Tensor> FullDeviatoricBasis() {
        const std::vector<Tensor> raw{{1, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0},
                                      {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}};
        std::vector<Tensor> dev;
        dev.reserve(raw.size());
        for (const auto &t : raw)
            dev.push_back(Deviatoric(t));
        return Orthonormalise(dev);
    }

    // What is left of the deviatoric space once the symmetry-allowed directions are taken out. In those
    // directions the true tensor is exactly zero by symmetry, whatever the crystal is - so whatever is
    // measured there is this dataset's own systematic error, measured on this dataset.
    std::vector<Tensor> ForbiddenBasis(const std::vector<Tensor> &full, const std::vector<Tensor> &allowed) {
        std::vector<Tensor> raw;
        for (Tensor t : full) {
            for (const auto &a : allowed)
                t = t - Scale(a, FrobDot(t, a));
            raw.push_back(t);
        }
        return Orthonormalise(raw);
    }

    // ---------------------------------------------------------------- binning
    // A Fibonacci spiral on the hemisphere: near-uniform, and no pole or seam where a lattice direction
    // could pile up.
    std::vector<gemmi::Vec3> HemisphereDirections(int n) {
        std::vector<gemmi::Vec3> d;
        d.reserve(n);
        const double golden = gemmi::pi() * (1.0 + std::sqrt(5.0));
        for (int i = 0; i < n; ++i) {
            const double z = (i + 0.5) / n;
            const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
            const double phi = golden * (i + 0.5);
            d.emplace_back(r * std::cos(phi), r * std::sin(phi), z);
        }
        return d;
    }

    struct Cells {
        std::vector<double> mu, se, s2;
        std::vector<std::array<double, 6>> q;
        std::vector<int> shell, n;
        // The cluster structure, kept so the covariance can be made robust to it: one entry per
        // measurement that went into a cell, grouped by the unique reflection it belongs to.
        std::vector<int32_t> entry_cell, entry_cluster;
        std::vector<float> entry_value;
        int n_shells = 0;
    };

    struct Usable {
        std::array<int, 3> hkl;
        double value;      // epsilon-corrected merged intensity, negatives kept
        double s2;
        int shell = 0;
    };

    std::vector<Usable> SelectReflections(const std::vector<MergedReflection> &merged,
                                          const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg,
                                          int n_shells) {
        const gemmi::GroupOps gops = sg ? sg->operations() : gemmi::GroupOps{};
        std::vector<Usable> u;
        u.reserve(merged.size());
        for (const auto &r : merged) {
            if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f))
                continue;
            const std::array<int, 3> hkl{r.h, r.k, r.l};
            const double eps = sg ? std::max(1, gops.epsilon_factor(hkl)) : 1;
            u.push_back({hkl, r.I / eps, 1.0 / (static_cast<double>(r.d) * r.d), 0});
        }
        if (u.empty())
            return u;
        // Equal-count shells in s^2. Every symmetry copy of a reflection has the same |s|, so equal-count
        // over the expanded set and over the unique set are the same edges.
        std::vector<int> order(u.size());
        std::iota(order.begin(), order.end(), 0);
        std::sort(order.begin(), order.end(), [&](int a, int b) { return u[a].s2 < u[b].s2; });
        for (size_t i = 0; i < order.size(); ++i)
            u[order[i]].shell = std::min(n_shells - 1, static_cast<int>(i * n_shells / order.size()));
        return u;
    }

    // Keep the cells that hold enough measurements, renumber them, and carry the entry list over.
    void FinishCells(Cells &c, const std::vector<double> &sum, const std::vector<double> &sumsq,
                     const std::vector<double> &s2sum, const std::vector<std::array<double, 6>> &qsum,
                     const std::vector<int> &count, const std::vector<int32_t> &raw_cell,
                     const std::vector<int32_t> &raw_cluster, const std::vector<float> &raw_value) {
        std::vector<int32_t> remap(count.size(), -1);
        for (size_t cid = 0; cid < count.size(); ++cid) {
            if (count[cid] < MIN_CELL_REFLECTIONS)
                continue;
            const double n = count[cid];
            const double mean = sum[cid] / n;
            const double var = std::max(0.0, (sumsq[cid] - n * mean * mean) / (n - 1.0));
            remap[cid] = static_cast<int32_t>(c.mu.size());
            c.mu.push_back(mean);
            c.se.push_back(std::sqrt(var / n));
            c.s2.push_back(s2sum[cid] / n);
            c.q.push_back({qsum[cid][0] / n, qsum[cid][1] / n, qsum[cid][2] / n,
                           qsum[cid][3] / n, qsum[cid][4] / n, qsum[cid][5] / n});
            c.shell.push_back(static_cast<int>(cid) / N_DIRECTIONS);
            c.n.push_back(count[cid]);
        }
        // A cell whose measurements happen to agree exactly has no error bar of its own; give it the
        // typical one rather than an infinite weight.
        std::vector<double> positive;
        for (double s : c.se)
            if (s > 0.0)
                positive.push_back(s);
        if (!positive.empty()) {
            std::nth_element(positive.begin(), positive.begin() + positive.size() / 2, positive.end());
            const double median = positive[positive.size() / 2];
            for (double &s : c.se)
                if (!(s > 0.0))
                    s = median;
        }
        for (size_t e = 0; e < raw_cell.size(); ++e)
            if (remap[raw_cell[e]] >= 0) {
                c.entry_cell.push_back(remap[raw_cell[e]]);
                c.entry_cluster.push_back(raw_cluster[e]);
                c.entry_value.push_back(raw_value[e]);
            }
    }

    Cells BuildCells(const std::vector<Usable> &refl, const gemmi::UnitCell &cell,
                     const gemmi::SpaceGroup *sg, int n_shells) {
        Cells c;
        c.n_shells = n_shells;
        const std::vector<gemmi::Op> ops = sg ? sg->operations().sym_ops : std::vector<gemmi::Op>{gemmi::Op::identity()};
        const auto dirs = HemisphereDirections(N_DIRECTIONS);
        // Friedel needs no separate operation: -s folds onto the same hemisphere direction as s.
        const size_t stride = std::max<size_t>(1, (refl.size() * ops.size() + MAX_ENTRIES - 1) / MAX_ENTRIES);

        const int n_cells_total = n_shells * N_DIRECTIONS;
        std::vector<double> sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0);
        std::vector<std::array<double, 6>> qsum(n_cells_total, {0, 0, 0, 0, 0, 0});
        std::vector<int> count(n_cells_total, 0);
        std::vector<int32_t> raw_cell, raw_cluster;
        std::vector<float> raw_value;

        std::vector<std::pair<int, gemmi::Vec3>> hits;      // (cell, s) of this reflection's copies
        for (size_t i = 0; i < refl.size(); i += stride) {
            const auto &r = refl[i];
            hits.clear();
            for (const auto &op : ops) {
                const gemmi::Op::Miller m = op.apply_to_hkl(r.hkl);
                gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2]));
                if (s.z < 0)
                    s = gemmi::Vec3(-s.x, -s.y, -s.z);
                const double len = s.length();
                if (!(len > 0))
                    continue;
                const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len);
                int best = 0;
                double best_dot = -1.0;
                for (int j = 0; j < N_DIRECTIONS; ++j) {
                    const double dot = std::fabs(unit.dot(dirs[j]));
                    if (dot > best_dot) { best_dot = dot; best = j; }
                }
                hits.emplace_back(r.shell * N_DIRECTIONS + best, s);
            }
            // A reflection contributes to a cell ONCE however many of its symmetry copies land there, so
            // a cell's standard error counts independent measurements.
            std::sort(hits.begin(), hits.end(), [](const auto &a, const auto &b) { return a.first < b.first; });
            const int32_t cluster = static_cast<int32_t>(i);
            int previous = -1;
            for (const auto &[cid, s] : hits) {
                if (cid == previous)
                    continue;
                previous = cid;
                sum[cid] += r.value;
                sumsq[cid] += r.value * r.value;
                s2sum[cid] += r.s2;
                qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z;
                qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z;
                count[cid] += 1;
                raw_cell.push_back(static_cast<int32_t>(cid));
                raw_cluster.push_back(cluster);
                raw_value.push_back(static_cast<float>(r.value));
            }
        }

        FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value);
        return c;
    }

    // ---------------------------------------------------------------- the tensor fit
    struct TensorFit {
        bool ok = false;
        Tensor b{0, 0, 0, 0, 0, 0};
        Eigen::MatrixXd cov;      // nb x nb, cluster-robust
        double chi2red = NAN;
    };

    // mu_cell = K_shell * exp(-1/2 s^T B s), fitted on the INTENSITY scale by weighted Gauss-Newton with
    // weights 1/SE^2, one free K per resolution shell profiled alongside the tensor. Because K is free per
    // shell, every isotropic feature - the Wilson curve, an ice ring, a noise floor, a scaling error - is
    // absorbed exactly, and only the l=2 angular part drives the tensor. Nothing is censored: a cell whose
    // mean is zero or negative is kept, which is what keeps the fit honest in a direction that has died.
    //
    // The covariance is a cluster-robust sandwich, clustering by unique reflection. That is not a detail:
    // symmetry mates of one reflection share their true |F|^2, so a symmetry-forbidden tensor direction -
    // determined from contrasts WITHIN a cluster - carries none of the Wilson scatter that a
    // symmetry-allowed direction carries, and a covariance that assumed independent cells would get the
    // ratio of the two backwards by a factor of several.
    //
    // Following Popov & Bourenkov (2003) Acta Cryst. D59, 1145-1153
    TensorFit FitTensor(const Cells &c, const std::vector<Tensor> &basis, bool want_covariance) {
        TensorFit fit;
        const int nb = static_cast<int>(basis.size());
        const int nc = static_cast<int>(c.mu.size());
        if (nb == 0 || nc < MIN_CELLS)
            return fit;

        // Shells that carry a positive mean intensity: ln K is undefined for one that does not, and a
        // shell with no signal at all cannot say anything about direction either.
        std::vector<double> shell_num(c.n_shells, 0.0), shell_den(c.n_shells, 0.0);
        for (int i = 0; i < nc; ++i) {
            const double w = 1.0 / (c.se[i] * c.se[i]);
            shell_num[c.shell[i]] += w * c.mu[i];
            shell_den[c.shell[i]] += w;
        }
        std::vector<int> slot(c.n_shells, -1);
        int ns = 0;
        std::vector<double> ln_k;
        for (int s = 0; s < c.n_shells; ++s)
            if (shell_den[s] > 0.0 && shell_num[s] / shell_den[s] > 0.0) {
                slot[s] = ns++;
                ln_k.push_back(std::log(shell_num[s] / shell_den[s]));
            }
        std::vector<int> keep;
        for (int i = 0; i < nc; ++i)
            if (slot[c.shell[i]] >= 0)
                keep.push_back(i);
        if (static_cast<int>(keep.size()) < MIN_CELLS || ns == 0)
            return fit;

        const int n = static_cast<int>(keep.size());
        const int np = ns + nb;
        Eigen::MatrixXd G(n, nb);
        Eigen::VectorXd mu(n), se(n);
        std::vector<int> row_slot(n);
        for (int i = 0; i < n; ++i) {
            const int ci = keep[i];
            for (int k = 0; k < nb; ++k)
                G(i, k) = -0.5 * QuadForm(basis[k], c.q[ci]);
            mu(i) = c.mu[ci];
            se(i) = c.se[ci];
            row_slot[i] = slot[c.shell[ci]];
        }

        Eigen::VectorXd beta = Eigen::VectorXd::Zero(np);
        for (int s = 0; s < ns; ++s)
            beta(s) = ln_k[s];

        auto predict = [&](const Eigen::VectorXd &b, Eigen::VectorXd &pred) {
            for (int i = 0; i < n; ++i) {
                double e = b(row_slot[i]);
                for (int k = 0; k < nb; ++k)
                    e += G(i, k) * b(ns + k);
                pred(i) = std::exp(std::clamp(e, -60.0, 60.0));
            }
        };
        auto chi2_of = [&](const Eigen::VectorXd &pred) {
            double s = 0.0;
            for (int i = 0; i < n; ++i) {
                const double r = (mu(i) - pred(i)) / se(i);
                s += r * r;
            }
            return s;
        };

        Eigen::VectorXd pred(n), resid(n);
        predict(beta, pred);
        double chi2 = chi2_of(pred);
        Eigen::MatrixXd J(n, np);
        for (int iter = 0; iter < 60; ++iter) {
            J.setZero();
            for (int i = 0; i < n; ++i) {
                const double f = -pred(i) / se(i);
                J(i, row_slot[i]) = f;
                for (int k = 0; k < nb; ++k)
                    J(i, ns + k) = f * G(i, k);
                resid(i) = (mu(i) - pred(i)) / se(i);
            }
            const Eigen::MatrixXd A = J.transpose() * J;
            const Eigen::VectorXd grad = J.transpose() * resid;
            const Eigen::VectorXd step =
                    A.ldlt().solve(-grad);                    // Gauss-Newton on chi2/2
            if (!step.allFinite())
                break;
            // Halve the step until chi2 falls; a fit that cannot improve at all is converged.
            double t = 1.0;
            bool improved = false;
            Eigen::VectorXd trial(np), trial_pred(n);
            for (int back = 0; back < 20; ++back) {
                trial = beta + t * step;
                predict(trial, trial_pred);
                const double trial_chi2 = chi2_of(trial_pred);
                if (trial_chi2 < chi2) {
                    const bool converged = chi2 - trial_chi2 < 1e-10 * std::max(1.0, chi2);
                    beta = trial;
                    pred = trial_pred;
                    chi2 = trial_chi2;
                    improved = !converged;
                    break;
                }
                t *= 0.5;
            }
            if (!improved)
                break;
        }

        Tensor b{0, 0, 0, 0, 0, 0};
        for (int k = 0; k < nb; ++k)
            b = b + Scale(basis[k], beta(ns + k));
        fit.b = b;
        fit.chi2red = chi2 / std::max(1, n - np);
        fit.ok = true;
        if (!want_covariance)
            return fit;

        // Sandwich: A^-1 B A^-1 with A the Gauss-Newton normal matrix and B the sum of outer products of
        // per-cluster scores. A cell mean is the mean of its reflections, so a reflection's share of the
        // cell's score is its own share of the cell's residual.
        for (int i = 0; i < n; ++i) {
            const double f = -pred(i) / se(i);
            J(i, row_slot[i]) = f;
            for (int k = 0; k < nb; ++k)
                J(i, ns + k) = f * G(i, k);
        }
        const Eigen::MatrixXd A = J.transpose() * J;
        std::vector<int32_t> cell_row(c.mu.size(), -1);
        for (int i = 0; i < n; ++i)
            cell_row[keep[i]] = i;
        Eigen::MatrixXd B = Eigen::MatrixXd::Zero(np, np);
        Eigen::VectorXd score = Eigen::VectorXd::Zero(np);
        int32_t current = -1;
        auto flush = [&]() {
            if (current >= 0)
                B.noalias() += score * score.transpose();
            score.setZero();
        };
        for (size_t e = 0; e < c.entry_cell.size(); ++e) {
            if (c.entry_cluster[e] != current) {
                flush();
                current = c.entry_cluster[e];
            }
            const int i = cell_row[c.entry_cell[e]];
            if (i < 0)
                continue;
            const int ci = keep[i];
            const double w = (c.entry_value[e] - pred(i)) * pred(i)
                             / (c.n[ci] * c.se[ci] * c.se[ci]);
            score(row_slot[i]) += w;
            for (int k = 0; k < nb; ++k)
                score(ns + k) += w * G(i, k);
        }
        flush();
        const Eigen::MatrixXd ainv = A.ldlt().solve(Eigen::MatrixXd::Identity(np, np));
        const Eigen::MatrixXd full = ainv * B * ainv;
        fit.cov = full.bottomRightCorner(nb, nb);
        return fit;
    }

    struct Eigen3 {
        double value[3];
        double vec[3][3];
    };

    Eigen3 Decompose(const Tensor &b) {
        Eigen3 out{};
        double d[3];
        const gemmi::Mat33 v = gemmi::eigen_decomposition(b, d);
        int idx[3] = {0, 1, 2};
        std::sort(idx, idx + 3, [&](int a, int c) { return d[a] > d[c]; });
        for (int n = 0; n < 3; ++n) {
            out.value[n] = d[idx[n]];
            for (int j = 0; j < 3; ++j)
                out.vec[n][j] = v[j][idx[n]];    // eigenvectors are the columns
        }
        return out;
    }

    double DeltaB(const Tensor &b) {
        double d[3];
        gemmi::eigen_decomposition(b, d);
        return *std::max_element(d, d + 3) - *std::min_element(d, d + 3);
    }

    // The dB that counting statistics alone would produce in a subspace: E[dB] = k * sigma, with sigma the
    // per-component standard error and k the constant for that subspace's dimension.
    double CountingNull(const Eigen::MatrixXd &cov, double k) {
        if (cov.rows() == 0)
            return 0.0;
        return k * std::sqrt(std::max(0.0, cov.trace() / cov.rows()));
    }

    // ---------------------------------------------------------------- the s^2 signature
    struct ShapeFit {
        // The regressions below all ran. Distinct from shells: there can be plenty of shells and
        // still no fit, and a ShapeFit that returned early carries NaN coefficients with flat and
        // convex both false - which reads as a clean "Linear" to anyone testing those two alone.
        bool ok = false;
        int shells = 0;
        double c0 = NAN, c0_err = NAN, c1 = NAN, slope_through_origin = NAN;
        double curvature = NAN, curvature_err = NAN, curvature_share = NAN;
        double residual = NAN;
        bool flat = false, convex = false;
    };

    // One shell's l=2 amplitude: mu_cell = K exp(A u_cell), u = -1/2 (s^T bhat s)/|s|^2, K profiled out,
    // fitted on the intensity scale so no cell is dropped. A genuine tensor gives A = |B_dev| s^2, i.e. a
    // straight line through the origin; the fit makes no assumption at all about how A depends on s.
    bool FitShellAmplitude(const std::vector<double> &u, const std::vector<double> &mu,
                           const std::vector<double> &se, double &a_out, double &a_err) {
        auto chi2 = [&](double a) {
            double num = 0.0, den = 0.0;
            std::vector<double> g(u.size());
            for (size_t i = 0; i < u.size(); ++i) {
                g[i] = std::exp(std::clamp(a * u[i], -60.0, 60.0));
                const double w = 1.0 / (se[i] * se[i]);
                num += mu[i] * g[i] * w;
                den += g[i] * g[i] * w;
            }
            const double k = den > 0.0 ? num / den : 0.0;
            double s = 0.0;
            for (size_t i = 0; i < u.size(); ++i) {
                const double r = (mu[i] - k * g[i]) / se[i];
                s += r * r;
            }
            return s;
        };
        // A coarse scan brackets the minimum, golden section refines it. A is |B_dev| s^2, so a few
        // hundred covers any crystal; the search is on a smooth one-parameter curve.
        double lo = -300.0, hi = 300.0, best = 0.0, best_chi2 = chi2(0.0);
        for (int i = 0; i <= 120; ++i) {
            const double a = lo + (hi - lo) * i / 120.0;
            const double v = chi2(a);
            if (v < best_chi2) { best_chi2 = v; best = a; }
        }
        lo = best - 5.0;
        hi = best + 5.0;
        constexpr double INV_PHI = 0.6180339887498949;
        double x1 = hi - INV_PHI * (hi - lo), x2 = lo + INV_PHI * (hi - lo);
        double f1 = chi2(x1), f2 = chi2(x2);
        for (int i = 0; i < 60; ++i) {
            if (f1 < f2) { hi = x2; x2 = x1; f2 = f1; x1 = hi - INV_PHI * (hi - lo); f1 = chi2(x1); }
            else { lo = x1; x1 = x2; f1 = f2; x2 = lo + INV_PHI * (hi - lo); f2 = chi2(x2); }
        }
        a_out = 0.5 * (lo + hi);
        // One-parameter Gauss-Newton error bar, scaled by the residual so a shell the model does not
        // describe reports a large one.
        double num = 0.0, den = 0.0;
        std::vector<double> g(u.size());
        for (size_t i = 0; i < u.size(); ++i) {
            g[i] = std::exp(std::clamp(a_out * u[i], -60.0, 60.0));
            const double w = 1.0 / (se[i] * se[i]);
            num += mu[i] * g[i] * w;
            den += g[i] * g[i] * w;
        }
        const double k = den > 0.0 ? num / den : 0.0;
        double jtj = 0.0, ss = 0.0;
        for (size_t i = 0; i < u.size(); ++i) {
            const double dr = k * g[i] * u[i] / se[i];
            jtj += dr * dr;
            const double r = (mu[i] - k * g[i]) / se[i];
            ss += r * r;
        }
        if (!(jtj > 0.0) || u.size() < 3)
            return false;
        a_err = std::sqrt(ss / (u.size() - 2) / jtj);
        return std::isfinite(a_out) && std::isfinite(a_err) && a_err > 0.0;
    }

    // Weighted least squares of `y` on the columns of `x`; returns the parameters and the covariance.
    bool Wls(const Eigen::MatrixXd &x, const Eigen::VectorXd &y, const Eigen::VectorXd &w,
             Eigen::VectorXd &p, Eigen::MatrixXd &cov, double &chi2) {
        const Eigen::MatrixXd xtw = x.transpose() * w.asDiagonal();
        const Eigen::MatrixXd normal = xtw * x;
        const Eigen::LDLT<Eigen::MatrixXd> ldlt(normal);
        if (ldlt.info() != Eigen::Success)
            return false;
        p = ldlt.solve(xtw * y);
        const Eigen::VectorXd r = y - x * p;
        chi2 = r.transpose() * w.asDiagonal() * r;
        cov = ldlt.solve(Eigen::MatrixXd::Identity(x.cols(), x.cols()));
        return p.allFinite();
    }

    ShapeFit FitShape(const Cells &c, const Tensor &tensor) {
        ShapeFit out;
        const double norm = std::sqrt(FrobDot(tensor, tensor));
        if (!(norm > 0.0))
            return out;
        const Tensor bhat = Scale(tensor, 1.0 / norm);

        std::vector<double> xs, as, errs;
        for (int s = 0; s < c.n_shells; ++s) {
            std::vector<double> u, mu, se;
            double s2sum = 0.0;
            for (size_t i = 0; i < c.mu.size(); ++i) {
                if (c.shell[i] != s)
                    continue;
                u.push_back(-0.5 * QuadForm(bhat, c.q[i]) / std::max(c.s2[i], 1e-12));
                mu.push_back(c.mu[i]);
                se.push_back(c.se[i]);
                s2sum += c.s2[i];
            }
            if (static_cast<int>(u.size()) < MIN_CELLS_PER_SHAPE_SHELL)
                continue;
            const auto [lo, hi] = std::minmax_element(u.begin(), u.end());
            if (*hi - *lo < 1e-6)
                continue;
            double a = 0.0, err = 0.0;
            if (!FitShellAmplitude(u, mu, se, a, err))
                continue;
            xs.push_back(s2sum / u.size());
            as.push_back(a);
            errs.push_back(err);
        }
        out.shells = static_cast<int>(xs.size());
        if (out.shells < MIN_SHAPE_SHELLS)
            return out;

        const int n = out.shells;
        Eigen::VectorXd y(n), w(n);
        for (int i = 0; i < n; ++i) {
            y(i) = as[i];
            w(i) = 1.0 / (errs[i] * errs[i]);
        }
        Eigen::VectorXd p;
        Eigen::MatrixXd cov;
        double chi2 = 0.0;
        Eigen::MatrixXd x1(n, 1), x2(n, 2), xq(n, 2);
        for (int i = 0; i < n; ++i) {
            x1(i, 0) = xs[i];
            x2(i, 0) = 1.0;
            x2(i, 1) = xs[i];
            xq(i, 0) = xs[i];
            xq(i, 1) = xs[i] * xs[i];
        }
        if (!Wls(x1, y, w, p, cov, chi2))
            return out;
        out.slope_through_origin = p(0);
        out.residual = chi2 / std::max(1, n - 1);
        if (!Wls(x2, y, w, p, cov, chi2))
            return out;
        out.c0 = p(0);
        out.c0_err = std::sqrt(std::max(0.0, cov(0, 0)));
        out.c1 = p(1);
        if (Wls(xq, y, w, p, cov, chi2)) {
            out.curvature = p(1);
            out.curvature_err = std::sqrt(std::max(0.0, cov(1, 1)));
            // How much of the deficit at the resolution limit the s^4 term carries.
            const double s2max = *std::max_element(xs.begin(), xs.end());
            const double total = p(0) * s2max + p(1) * s2max * s2max;
            out.curvature_share = total != 0.0 ? p(1) * s2max * s2max / total : NAN;
        }
        out.flat = out.c0 > FLAT_INTERCEPT && out.c0_err > 0.0 && out.c0 / out.c0_err > SHAPE_Z;
        out.convex = !out.flat && out.curvature > 0.0 && out.curvature_err > 0.0
                     && out.curvature / out.curvature_err > SHAPE_Z
                     && out.curvature_share > CONVEX_SHARE;
        out.ok = true;
        return out;
    }

    // ---------------------------------------------------------------- directional limits
    // <I/sigma(I)> in a cone about each principal direction, read where it falls through 2.0. This uses no
    // model of the fall-off at all, which is why it and the tensor are reported side by side.
    //
    // Following Evans & Murshudov (2013) Acta Cryst. D69, 1204-1214
    void ConeLimits(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
                    const gemmi::SpaceGroup *sg, const Eigen3 &axes, int n_shells,
                    double (&d_min)[3], bool (&censored)[3]) {
        const std::vector<gemmi::Op> ops = sg ? sg->operations().sym_ops
                                              : std::vector<gemmi::Op>{gemmi::Op::identity()};
        std::vector<double> s2;
        s2.reserve(merged.size());
        for (const auto &r : merged)
            if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f)
                s2.push_back(1.0 / (static_cast<double>(r.d) * r.d));
        if (s2.size() < static_cast<size_t>(n_shells * MIN_CELL_REFLECTIONS))
            return;
        std::vector<double> sorted = s2;
        std::sort(sorted.begin(), sorted.end());
        std::vector<double> edge(n_shells + 1);
        for (int i = 0; i <= n_shells; ++i)
            edge[i] = sorted[std::min(sorted.size() - 1, sorted.size() * i / n_shells)];

        const double cos_cone = std::cos(CONE_HALF_ANGLE_DEG * gemmi::pi() / 180.0);
        std::vector<std::vector<double>> isig_sum(3, std::vector<double>(n_shells, 0.0));
        std::vector<std::vector<int>> isig_n(3, std::vector<int>(n_shells, 0));
        size_t iu = 0;
        for (const auto &r : merged) {
            if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || !(r.d > 0.0f))
                continue;
            const double this_s2 = s2[iu++];
            int shell = 0;
            while (shell + 1 < n_shells && this_s2 > edge[shell + 1])
                ++shell;
            bool in_cone[3] = {false, false, false};
            for (const auto &op : ops) {
                const gemmi::Op::Miller m = op.apply_to_hkl({{r.h, r.k, r.l}});
                const gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(m[0], m[1], m[2]));
                const double len = s.length();
                if (!(len > 0))
                    continue;
                for (int n = 0; n < 3; ++n) {
                    const double dot = (s.x * axes.vec[n][0] + s.y * axes.vec[n][1] + s.z * axes.vec[n][2]) / len;
                    if (std::fabs(dot) >= cos_cone)
                        in_cone[n] = true;
                }
            }
            for (int n = 0; n < 3; ++n)
                if (in_cone[n]) {
                    isig_sum[n][shell] += r.I / r.sigma;
                    isig_n[n][shell] += 1;
                }
        }
        for (int n = 0; n < 3; ++n) {
            double last_s2 = NAN, last_isig = NAN, limit_s2 = NAN;
            bool crossed = false;
            for (int s = 0; s < n_shells; ++s) {
                if (isig_n[n][s] < MIN_CELL_REFLECTIONS)
                    continue;
                const double mid = 0.5 * (edge[s] + edge[s + 1]);
                const double isig = isig_sum[n][s] / isig_n[n][s];
                if (isig >= CONE_I_OVER_SIGMA) {
                    last_s2 = mid;
                    last_isig = isig;
                    limit_s2 = edge[s + 1];
                } else if (std::isfinite(last_s2)) {
                    // Linear interpolation of <I/sigma> in s^2 between the last shell above the threshold
                    // and the first below it.
                    const double f = (last_isig - CONE_I_OVER_SIGMA) / (last_isig - isig);
                    limit_s2 = last_s2 + f * (mid - last_s2);
                    crossed = true;
                    break;
                } else {
                    break;
                }
            }
            if (std::isfinite(limit_s2) && limit_s2 > 0.0) {
                d_min[n] = 1.0 / std::sqrt(limit_s2);
                // <I/sigma> never fell through the threshold: the outermost shell that carries this cone
                // is still above it, so the value is the edge of the measured data and not this crystal's
                // own limit along that direction.
                censored[n] = !crossed;
            }
        }
    }

    const char *Band(double g) {
        if (g < GATE_MARGINAL) return "not established";
        if (g < GATE_ESTABLISHED) return "marginal";
        if (g < GATE_STRONG) return "established";
        return "strong";
    }

    // ---------------------------------------------------------------- the systematic-error floor
    // The same cells, built from UNMERGED observations at the Miller index each was measured at, with
    // the unique reflection they reduce to as the cluster. Nothing is symmetry-expanded and nothing is
    // de-duplicated: two observations of the same reflection recorded on different frames are two
    // measurements, and the difference between them is exactly the signal this arm is after.
    Cells BuildObservationCells(const std::vector<AnisotropyObservation> &obs, const gemmi::UnitCell &cell,
                                const gemmi::SpaceGroup *sg, int n_shells) {
        Cells c;
        c.n_shells = n_shells;
        const HKLKeyGenerator key(/*merge_friedel=*/true, sg ? *sg : *gemmi::find_spacegroup_by_number(1));

        // Order the observations by the unique reflection they belong to, so the cluster-robust sum
        // below can walk them one cluster at a time; then stride whole clusters if there are more than
        // the fit needs. The estimate is flat in the number of unique reflections over a factor of 12,
        // so striding costs nothing and bounds the work.
        std::vector<std::pair<uint64_t, int32_t>> order;
        order.reserve(obs.size());
        for (size_t i = 0; i < obs.size(); ++i) {
            const auto &o = obs[i];
            if (!std::isfinite(o.I) || !std::isfinite(o.sigma) || o.sigma <= 0.0f || !(o.d > 0.0f))
                continue;
            order.emplace_back(key(o.h, o.k, o.l).pack(), static_cast<int32_t>(i));
        }
        if (order.empty())
            return c;
        std::sort(order.begin(), order.end());
        int32_t n_clusters = 1;
        for (size_t i = 1; i < order.size(); ++i)
            if (order[i].first != order[i - 1].first)
                ++n_clusters;
        const int cluster_stride = std::max<int>(1, n_clusters / MAX_CLUSTERS);

        const auto dirs = HemisphereDirections(N_DIRECTIONS);
        const int n_cells_total = n_shells * N_DIRECTIONS;
        std::vector<double> sum(n_cells_total, 0.0), sumsq(n_cells_total, 0.0), s2sum(n_cells_total, 0.0);
        std::vector<std::array<double, 6>> qsum(n_cells_total, {0, 0, 0, 0, 0, 0});
        std::vector<int> count(n_cells_total, 0);
        std::vector<int32_t> raw_cell, raw_cluster;
        std::vector<float> raw_value;

        std::vector<int32_t> kept;
        std::vector<int32_t> cluster_of;
        uint64_t previous_key = 0;
        int32_t seen = -1;
        for (const auto &[k, idx] : order) {
            if (seen < 0 || k != previous_key) {
                previous_key = k;
                ++seen;
            }
            if (seen % cluster_stride != 0)
                continue;
            kept.push_back(idx);
            cluster_of.push_back(seen / cluster_stride);
        }
        if (kept.size() < static_cast<size_t>(MIN_CELLS * MIN_CELL_REFLECTIONS))
            return c;

        // Equal-count shells over the observations that are kept.
        std::vector<double> s2(kept.size());
        for (size_t i = 0; i < kept.size(); ++i) {
            const auto &o = obs[kept[i]];
            s2[i] = 1.0 / (static_cast<double>(o.d) * o.d);
        }
        std::vector<int32_t> by_s2(kept.size());
        std::iota(by_s2.begin(), by_s2.end(), 0);
        std::sort(by_s2.begin(), by_s2.end(), [&](int32_t a, int32_t b) { return s2[a] < s2[b]; });
        std::vector<int> shell(kept.size());
        for (size_t i = 0; i < by_s2.size(); ++i)
            shell[by_s2[i]] = std::min(n_shells - 1, static_cast<int>(i * n_shells / by_s2.size()));

        raw_cell.reserve(kept.size());
        raw_value.reserve(kept.size());
        raw_cluster.reserve(kept.size());
        for (size_t i = 0; i < kept.size(); ++i) {
            const auto &o = obs[kept[i]];
            gemmi::Vec3 s = cell.frac.mat.left_multiply(gemmi::Vec3(o.h, o.k, o.l));
            if (s.z < 0)
                s = gemmi::Vec3(-s.x, -s.y, -s.z);
            const double len = s.length();
            if (!(len > 0))
                continue;
            const gemmi::Vec3 unit(s.x / len, s.y / len, s.z / len);
            int best = 0;
            double best_dot = -1.0;
            for (int j = 0; j < N_DIRECTIONS; ++j) {
                const double dot = std::fabs(unit.dot(dirs[j]));
                if (dot > best_dot) { best_dot = dot; best = j; }
            }
            const int cid = shell[i] * N_DIRECTIONS + best;
            sum[cid] += o.I;
            sumsq[cid] += static_cast<double>(o.I) * o.I;
            s2sum[cid] += s2[i];
            qsum[cid][0] += s.x * s.x; qsum[cid][1] += s.y * s.y; qsum[cid][2] += s.z * s.z;
            qsum[cid][3] += s.x * s.y; qsum[cid][4] += s.x * s.z; qsum[cid][5] += s.y * s.z;
            count[cid] += 1;
            raw_cell.push_back(cid);
            raw_cluster.push_back(cluster_of[i]);
            raw_value.push_back(o.I);
        }
        FinishCells(c, sum, sumsq, s2sum, qsum, count, raw_cell, raw_cluster, raw_value);
        return c;
    }

    // The systematic error scale this dataset carries, and the part of it that is merely counting
    // noise, both measured in the tensor directions the Laue class forbids.
    struct SystematicFloor {
        bool ok = false;
        double sigma_excess = NAN;     // per tensor component, A^2
        double forbidden_z = NAN;      // that scatter over its own counting null
        int n_observations = 0;
    };

    SystematicFloor MeasureSystematicFloor(const std::vector<AnisotropyObservation> &obs,
                                           const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg,
                                           const std::vector<Tensor> &full,
                                           const std::vector<Tensor> &forbidden) {
        SystematicFloor out;
        const int q = static_cast<int>(forbidden.size());
        if (q == 0 || obs.empty())
            return out;
        const Cells cells = BuildObservationCells(obs, cell, sg, N_SHELLS);
        if (cells.mu.empty())
            return out;
        const TensorFit fit = FitTensor(cells, full, true);
        if (!fit.ok || fit.cov.rows() != static_cast<int>(full.size()))
            return out;

        Eigen::MatrixXd basis_map(static_cast<int>(full.size()), q);
        for (int j = 0; j < q; ++j)
            for (size_t k = 0; k < full.size(); ++k)
                basis_map(static_cast<int>(k), j) = FrobDot(full[k], forbidden[j]);
        Tensor forb{0, 0, 0, 0, 0, 0};
        Eigen::VectorXd theta(static_cast<int>(full.size()));
        for (size_t k = 0; k < full.size(); ++k)
            theta(static_cast<int>(k)) = FrobDot(fit.b, full[k]);
        const Eigen::VectorXd theta_forb = basis_map.transpose() * theta;
        for (int j = 0; j < q; ++j)
            forb = forb + Scale(forbidden[j], theta_forb(j));
        const Eigen::MatrixXd cov_forb = basis_map.transpose() * fit.cov * basis_map;

        const double k_q = K_FORBIDDEN[q];
        if (!(k_q > 0.0))
            return out;
        const double sigma_sys = DeltaB(forb) / k_q;
        const double sigma_stat = CountingNull(cov_forb, k_q) / k_q;
        // Take the counting part out: what is left is the systematic alone.
        out.sigma_excess = std::sqrt(std::max(0.0, sigma_sys * sigma_sys - sigma_stat * sigma_stat));
        if (sigma_stat > 0.0)
            out.forbidden_z = sigma_sys / sigma_stat;
        out.n_observations = static_cast<int>(cells.entry_cell.size());
        out.ok = true;
        return out;
    }
}

std::vector<AnisotropyObservation> ScaledObservations(const std::vector<IntegrationOutcome> &outcomes,
                                                     bool rotation, const gemmi::SpaceGroup *space_group,
                                                     double min_partiality) {
    // Per-image scale, indexed the way the outcomes are.
    std::vector<double> g(outcomes.size(), 0.0);
    for (size_t i = 0; i < outcomes.size(); ++i)
        if (outcomes[i].image_scale_g.has_value() && std::isfinite(*outcomes[i].image_scale_g)
            && *outcomes[i].image_scale_g > 0.0f)
            g[i] = *outcomes[i].image_scale_g;

    std::vector<AnisotropyObservation> out;
    // The same acceptance the merge itself uses (Merge.cpp): everything in an IntegrationOutcome came
    // out of the integrator, and Reflection::observed is not set by the _process.h5 reader, so testing
    // it would silently empty this on the --mode scale path.
    auto usable = [](const Reflection &r) {
        return std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && r.d > 0.0f
               && std::isfinite(r.rlp) && r.rlp > 0.0f && std::isfinite(r.partiality);
    };

    if (!rotation) {
        // A still is a whole measurement of its reflection, and consecutive stills are different
        // crystals, so there is nothing to assemble.
        for (size_t i = 0; i < outcomes.size(); ++i) {
            if (!(g[i] > 0.0))
                continue;
            for (const auto &r : outcomes[i].reflections) {
                if (!usable(r) || !(r.partiality >= min_partiality))
                    continue;
                const double corr = r.rlp / (r.partiality * g[i]);
                out.push_back({r.h, r.k, r.l, static_cast<float>(r.I * corr),
                               static_cast<float>(r.sigma * corr), r.d});
            }
        }
        return out;
    }

    // A rotation reflection is integrated image by image, so it arrives as a run of partials over
    // consecutive frames. Assemble each run into one full the way the 3D combine and the unmerged
    // export do - same raw hkl, frames no further apart than the combine's own gap, parts added with
    // their variances in quadrature - because a single partial divided by its own partiality carries
    // the rocking-curve model's error as well as its intensity, and that error is a function of the
    // reflection's direction relative to the spindle, which is precisely the direction this
    // diagnostic is measuring.
    constexpr float MAX_FRAME_GAP = 2.0f;
    // The sort key travels with the part instead of being read back through the pointer, the way the
    // merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
    // and an indirect compare is a cache miss on every one of them. The keys are the same values in
    // the same order, so introsort makes the same comparisons and the same swaps and leaves the same
    // order - which matters, because two parts can genuinely share (h,k,l) and image_number and the
    // event sums below are floating point.
    struct Part {
        int32_t h, k, l;
        float image_number;
        const Reflection *r;
        size_t image;
    };
    // Assemble a SAMPLE of the unique reflections, not all of them. Everything below exists to feed
    // MeasureSystematicFloor, which strides whole unique reflections down to MAX_CLUSTERS and on a
    // large sweep keeps about one in seventy of what it is handed; taking the sample here instead of
    // there removes the sort rather than shrinking it. Two things decide what the sample may be:
    //   * whole unique reflections. The floor is measured from the contrast BETWEEN a reflection's
    //     symmetry mates, so a reflection is taken with all of its mates or not at all - which is why
    //     this needs the space group the floor will be measured in, and samples on the ASU key.
    //   * coverage in direction and in resolution, which is what the fit's 12 shells x 60 directions
    //     read. A hash of the ASU key is a uniform sample of unique reflections and carries both; a
    //     stride over the parts, or over hkl order, carries neither.
    // How many unique reflections the partials cover is not known here, and does not need to be: one
    // part in 2^n is one unique reflection in 2^n whatever their number, so the rate comes from the
    // partial count, which is.
    size_t n_parts = 0;
    for (size_t i = 0; i < outcomes.size(); ++i)
        if (g[i] > 0.0)
            n_parts += outcomes[i].reflections.size();
    const HKLKeyGenerator asu_key(/*merge_friedel=*/true,
                                  space_group ? *space_group : *gemmi::find_spacegroup_by_number(1));
    const uint64_t sample_mask =
            std::bit_ceil((n_parts + MAX_SAMPLE_PARTS - 1) / MAX_SAMPLE_PARTS) - 1;
    std::vector<Part> parts;
    // Reserve what the sample holds rather than letting the vector grow: the growth chain copies a few
    // hundred megabytes on a large sweep.
    parts.reserve(std::min(n_parts, MAX_SAMPLE_PARTS));
    for (size_t i = 0; i < outcomes.size(); ++i) {
        if (!(g[i] > 0.0))
            continue;
        for (const auto &r : outcomes[i].reflections) {
            if (!usable(r))
                continue;
            if ((MixKey(asu_key(r.h, r.k, r.l).pack()) & sample_mask) != 0)
                continue;
            parts.push_back({r.h, r.k, r.l, r.image_number, &r, i});
        }
    }
    std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) {
        return std::tie(a.h, a.k, a.l, a.image_number) < std::tie(b.h, b.k, b.l, b.image_number);
    });
    for (size_t i = 0; i < parts.size();) {
        size_t j = i + 1;
        while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k
               && parts[j].l == parts[i].l
               && parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP)
            ++j;
        double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, p_g = 0.0;
        for (size_t m = i; m < j; ++m) {
            const Reflection &r = *parts[m].r;
            sum_p += r.partiality;
            sum_I += static_cast<double>(r.I) * r.rlp;
            sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
            p_g += r.partiality * g[parts[m].image];
        }
        const Reflection &first = *parts[i].r;
        i = j;
        // The scale is the rocking curve's own weighted mean over the frames the event spans; the
        // partiality divisor only corrects an event the sweep cut short, since a complete one sums
        // to 1 by construction.
        if (!(sum_p >= min_partiality) || !(p_g > 0.0))
            continue;
        const double scale = p_g / sum_p;               // partiality-weighted mean G
        const double corr = 1.0 / (sum_p * scale);
        out.push_back({first.h, first.k, first.l, static_cast<float>(sum_I * corr),
                       static_cast<float>(std::sqrt(sum_var) * corr), first.d});
    }
    return out;
}

const char *AnisotropyShapeCode(AnisotropyShape shape) {
    switch (shape) {
        case AnisotropyShape::Linear: return "LINEAR";
        case AnisotropyShape::Flat: return "FLAT";
        case AnisotropyShape::Convex: return "CONVEX";
        default: return "UNDETERMINED";
    }
}

const char *AnisotropyVerdictCode(AnisotropyVerdict verdict) {
    switch (verdict) {
        case AnisotropyVerdict::NotDetected: return "NOT_DETECTED";
        case AnisotropyVerdict::Detected: return "DETECTED";
        default: return "CANNOT_DETERMINE";
    }
}

AnisotropyResult AnalyzeAnisotropy(const std::vector<MergedReflection> &merged,
                                   const std::vector<AnisotropyObservation> &unmerged,
                                   const gemmi::UnitCell &cell,
                                   const gemmi::SpaceGroup *space_group,
                                   const AnisotropyRunInfo &run) {
    AnisotropyResult result;
    if (merged.empty() || !cell.is_crystal())
        return result;

    const auto allowed = AllowedDeviatoricBasis(cell, space_group);
    const auto full = FullDeviatoricBasis();
    const auto forbidden = ForbiddenBasis(full, allowed);
    const int p = static_cast<int>(allowed.size());
    const int q = static_cast<int>(forbidden.size());
    result.n_free_parameters = p;

    if (p == 0) {
        // Cubic. The anisotropy tensor is forced isotropic by symmetry, so its deviatoric part is exactly
        // zero - there is nothing to fit and nothing to gate.
        result.n_reflections = static_cast<int>(merged.size());
        result.delta_b = 0.0;
        result.delta_b_linear = 0.0;
        result.delta_b_flat = 0.0;
        result.eigenvalue[0] = result.eigenvalue[1] = result.eigenvalue[2] = 0.0;
        result.verdict = AnisotropyVerdict::NotDetected;
        result.band = "no anisotropy is possible in this Laue class";
        return result;
    }

    double isig_sum = 0.0;
    int isig_n = 0;
    for (const auto &r : merged)
        if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
            isig_sum += r.I / r.sigma;
            ++isig_n;
        }
    const double mean_isig = isig_n > 0 ? isig_sum / isig_n : 0.0;

    const auto refl = SelectReflections(merged, cell, space_group, N_SHELLS);
    if (refl.size() < static_cast<size_t>(MIN_CELLS * MIN_CELL_REFLECTIONS))
        return result;

    const Cells cells = BuildCells(refl, cell, space_group, N_SHELLS);
    const TensorFit constrained = FitTensor(cells, allowed, true);
    if (!constrained.ok)
        return result;
    result.n_reflections = static_cast<int>(refl.size());
    result.n_cells = static_cast<int>(cells.mu.size());

    const Eigen3 axes = Decompose(constrained.b);
    result.delta_b = axes.value[0] - axes.value[2];
    for (int n = 0; n < 3; ++n) {
        result.eigenvalue[n] = axes.value[n];
        for (int j = 0; j < 3; ++j)
            result.eigenvector[n][j] = axes.vec[n][j];
    }
    double d_min = 0.0;
    for (const auto &r : merged)
        if (r.d > 0.0f && (d_min == 0.0 || r.d < d_min))
            d_min = r.d;
    if (d_min > 0.0)
        result.fold_weakening = std::exp(result.delta_b / (2.0 * d_min * d_min));
    ConeLimits(merged, cell, space_group, axes, N_SHELLS, result.d_min_axis, result.d_min_censored);
    {
        double lo = INFINITY, hi = -INFINITY;
        for (double d : result.d_min_axis)
            if (std::isfinite(d)) { lo = std::min(lo, d); hi = std::max(hi, d); }
        if (std::isfinite(lo) && std::isfinite(hi))
            result.d_min_spread = hi - lo;
    }

    // --- the resolution signature, and whether the binning decides it ---
    const ShapeFit shape = FitShape(cells, constrained.b);
    // A fit that did not run was fitted over no shells. Reported that way rather than as the shells
    // it had available, so the renderer and the cautions below do not print NaN coefficients beside a
    // shell count that suggests they mean something.
    result.shape_shells = shape.ok ? shape.shells : 0;
    result.shape_intercept = shape.c0;
    result.shape_intercept_z = shape.c0_err > 0.0 ? shape.c0 / shape.c0_err : NAN;
    result.shape_slope = shape.c1;
    result.shape_curvature_z = shape.curvature_err > 0.0 ? shape.curvature / shape.curvature_err : NAN;
    result.shape_curvature_share = shape.curvature_share;
    result.shape_residual = shape.residual;
    if (shape.ok && shape.shells >= MIN_SHAPE_SHELLS) {
        result.shape = shape.flat ? AnisotropyShape::Flat
                     : shape.convex ? AnisotropyShape::Convex : AnisotropyShape::Linear;
        if (std::isfinite(shape.slope_through_origin) && shape.slope_through_origin != 0.0) {
            result.delta_b_linear = result.delta_b * shape.c1 / shape.slope_through_origin;
            result.delta_b_flat = result.delta_b - result.delta_b_linear;
        }
        // Rebin at 8 and at 16 shells and refit the tensor each time, so the direction the amplitude is
        // measured along is re-derived too. A verdict that moves is not a measurement.
        result.shape_stable = true;
        for (int alt : {8, 16}) {
            const auto alt_refl = SelectReflections(merged, cell, space_group, alt);
            const Cells alt_cells = BuildCells(alt_refl, cell, space_group, alt);
            const TensorFit alt_fit = FitTensor(alt_cells, allowed, false);
            if (!alt_fit.ok) {
                result.shape_stable = false;
                break;
            }
            const ShapeFit alt_shape = FitShape(alt_cells, alt_fit.b);
            if (!alt_shape.ok || alt_shape.shells < MIN_SHAPE_SHELLS || alt_shape.flat != shape.flat) {
                result.shape_stable = false;
                break;
            }
        }
        if (!result.shape_stable)
            result.shape = AnisotropyShape::Undetermined;
    }
    // The gate is applied to the part of delta_b that behaves as a Debye-Waller B, when that could be
    // separated; to delta_b itself when it could not.
    // A FLAT crystal can fit a negative linear slope; a negative anisotropy is not a smaller one, it is
    // no established Debye-Waller component at all.
    const double magnitude = std::max(0.0, std::isfinite(result.delta_b_linear) ? result.delta_b_linear
                                                                                : result.delta_b);

    // --- the gate ---
    auto refuse = [&](std::string why) {
        result.verdict = AnisotropyVerdict::CannotDetermine;
        result.refusal = std::move(why);
    };
    if (mean_isig < MIN_MEAN_I_OVER_SIGMA) {
        refuse("the merged data are at the noise floor (<I/sigma(I)> below 1), where neither the tensor "
               "nor its resolution signature means anything");
    } else if (q == 0) {
        // Triclinic. Every other Laue class leaves directions in which the true tensor is zero by
        // symmetry, and measuring those is how this dataset's systematic error is established. Triclinic
        // leaves none, and three candidate substitutes (a degree-4 spherical-harmonic block, a
        // low/high-resolution split, a floor borrowed from similar datasets) were built and measured
        // against known truth: the first two are blind to real anisotropy, as a denominator must be, but
        // neither delivers a usable scale, and the borrowed floor spans three orders of magnitude.
        refuse("the Laue class is triclinic, which leaves no symmetry-forbidden tensor direction, so this "
               "dataset carries no measurement of its own systematic error");
    } else if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < MIN_ROTATION_DEG) {
        refuse("the observed rotation range is too short for the tensor to be separated from a lab-fixed "
               "systematic");
    } else if (std::isfinite(run.observed_rotation_deg) && !run.dose_term_in_scale_model) {
        // A dose ramp is not a significance failure - it is a confident false detection. A simulated
        // isotropic crystal given a 25 A^2 relative-B ramp over 180 deg returns dB ~ 9 A^2 with a
        // significance no test can see through, because the ramp really is a smooth quadratic function of
        // direction once the sweep geometry is folded in.
        refuse("the scale model carried no dose term, and an uncorrected dose ramp manufactures anisotropy "
               "that no significance test can see through");
    } else if (unmerged.empty()) {
        // Merged intensities have exact Laue symmetry by construction, so the symmetry-forbidden tensor
        // directions - the only place a dataset can measure its own systematic error - are identically
        // zero there whatever the crystal carries. Gating on what is left would be gating on counting
        // statistics, which against known ground truth calls "real" in 15-58% of clean isotropic datasets.
        refuse("the systematic-error scale can only be measured on unmerged observations, and none were "
               "available for this merge");
    } else {
        const SystematicFloor sys = MeasureSystematicFloor(unmerged, cell, space_group, full, forbidden);
        result.n_observations = sys.n_observations;
        result.forbidden_z = sys.forbidden_z;
        if (!sys.ok) {
            refuse("the symmetry-forbidden tensor directions of the unmerged observations could not be "
                   "measured, so this dataset's systematic-error scale is unknown");
        } else {
            const double k_p = K_ALLOWED[p];
            result.sigma_systematic = sys.sigma_excess;
            const double allowed_null = CountingNull(constrained.cov, k_p);
            result.floor = std::sqrt(k_p * result.sigma_systematic * k_p * result.sigma_systematic
                                     + allowed_null * allowed_null);
            if (result.floor > 0.0) {
                result.significance = magnitude / result.floor;
                result.detection_limit = GATE_ESTABLISHED * result.floor;
                result.band = Band(result.significance);
                result.verdict = result.significance > GATE_ESTABLISHED ? AnisotropyVerdict::Detected
                                                                        : AnisotropyVerdict::NotDetected;
            } else {
                refuse("the systematic-error floor could not be measured");
            }
        }
    }

    // --- cautions: measured things that bias the verdict, in either direction ---
    if (result.shape == AnisotropyShape::Flat)
        result.cautions.emplace_back(
                "the directional deficit does not follow exp(-1/2 s^T B s), so the fitted deltaB is a fit "
                "of the wrong functional form and may be an under-estimate");
    if (result.shape == AnisotropyShape::Undetermined && result.shape_shells >= MIN_SHAPE_SHELLS)
        result.cautions.emplace_back(
                "the resolution signature changed when the shells were rebinned, so it is reported as "
                "undetermined rather than as a measurement");
    if (std::isfinite(run.observed_rotation_deg) && run.observed_rotation_deg < 180.0)
        result.cautions.emplace_back(
                "the observed rotation range is under 180 deg, which leaves a second tensor direction "
                "reachable by a lab-fixed systematic and reduces the power of the test");
    if (std::isfinite(run.radiation_damage_relative_b) && std::fabs(run.radiation_damage_relative_b) > 10.0)
        result.cautions.emplace_back(
                "the run carries substantial radiation damage; the scale model removes its average "
                "monotone part, but a non-monotone dose ramp needs --relative-b to be taken out too");
    // R2's Laue-class caution, fired on its measured signature rather than on the Laue class alone.
    // Too high a symmetry assignment pushes real anisotropy into the directions the systematic error is
    // measured in, so sigma_sys is inflated by the signal itself and the gate collapses (measured: a
    // 42x inflation of sigma_sys, a gate of 106 falling to 3.1). The signature is a forbidden subspace
    // far above its own counting noise TOGETHER with a gate that has established nothing. p == 1 on its
    // own is not the signature - it is true of every tetragonal, trigonal and hexagonal dataset, 17 of
    // 38 in the battery, most of them correctly assigned.
    if (p == 1 && result.verdict == AnisotropyVerdict::NotDetected
        && std::isfinite(result.forbidden_z) && result.forbidden_z > FORBIDDEN_Z_SUSPECT)
        result.cautions.emplace_back(
                "the space group is worth checking before this verdict is believed: the tensor directions "
                "this Laue class forbids carry far more scatter than their own counting noise, yet nothing "
                "was established in the single direction it allows - which is what a too-high symmetry "
                "assignment looks like, real anisotropy being measured as systematic error");
    return result;
}

std::string AnisotropyToText(const AnisotropyResult &result) {
    std::ostringstream os;
    if (result.n_reflections == 0)
        return os.str();
    os << std::fixed;
    os << "Diffraction anisotropy\n";
    if (result.n_free_parameters == 0) {
        os << "  The Laue class forces the anisotropy tensor to be isotropic: deltaB = 0 exactly, with no\n"
           << "  free parameter to fit. This is symmetry, not a measurement.\n";
        return os.str();
    }
    os << std::setprecision(2);
    os << "  Anisotropic deltaB (range of principal components) = " << result.delta_b << " A^2"
       << "  [" << result.n_free_parameters << " free direction"
       << (result.n_free_parameters == 1 ? "" : "s") << " in this Laue class]\n";
    os << "  Principal components (A^2, relative):";
    for (double v : result.eigenvalue)
        os << " " << (v - result.eigenvalue[2]);
    os << "\n";
    if (std::isfinite(result.fold_weakening))
        os << "  Strongest / weakest direction at the resolution limit: " << result.fold_weakening << "x\n";
    // Every axis, not just the outer two: ConeLimits fills each independently, so a crystal whose
    // first and third directions fall in cones too sparse to measure - but whose second is well
    // covered - would otherwise suppress the whole line, d_min_spread included.
    if (std::isfinite(result.d_min_axis[0]) || std::isfinite(result.d_min_axis[1])
        || std::isfinite(result.d_min_axis[2])) {
        bool any_censored = false;
        os << "  d_min along the principal directions (<I/sigma(I)> = " << CONE_I_OVER_SIGMA
           << " in a " << CONE_HALF_ANGLE_DEG << " deg cone):";
        for (int n = 0; n < 3; ++n) {
            if (!std::isfinite(result.d_min_axis[n])) {
                os << " -";
            } else if (result.d_min_censored[n]) {
                os << " <" << result.d_min_axis[n];
                any_censored = true;
            } else {
                os << " " << result.d_min_axis[n];
            }
        }
        os << " A\n";
        if (any_censored)
            os << "     \"<\" marks a direction in which <I/sigma(I)> had not fallen to "
               << CONE_I_OVER_SIGMA << " where the measured data\n"
               << "     end: that figure is the edge of the data, not this crystal's limit along that\n"
               << "     direction, and the spread between the directions is therefore a lower bound.\n";
    }
    if (result.shape_shells >= MIN_SHAPE_SHELLS) {
        os << "  Resolution signature of the deficit over " << result.shape_shells << " shells: "
           << AnisotropyShapeCode(result.shape) << "  (intercept " << std::showpos << result.shape_intercept
           << std::noshowpos << ", z = " << result.shape_intercept_z << "; slope " << result.shape_slope
           << " A^2; through-origin chi2/dof " << result.shape_residual << ")\n";
        if (result.shape == AnisotropyShape::Flat)
            os << "     The deficit does not follow exp(-1/2 s^T B s): " << result.delta_b_flat
               << " A^2 of the deltaB above does not behave as a Debye-Waller B. The fall-off is being\n"
               << "     described by the wrong functional form, so the number may be an UNDER-estimate,\n"
               << "     not an over-estimate.\n";
        else if (result.shape == AnisotropyShape::Convex)
            os << "     The deficit grows faster than s^2, which a Debye-Waller B cannot do.\n";
    }
    if (result.verdict == AnisotropyVerdict::CannotDetermine) {
        os << "  => CANNOT DETERMINE: " << result.refusal << ".\n";
    } else {
        // The magnitude judged is the Debye-Waller part of the deltaB, which is what the floor is a
        // floor on; on a crystal whose deficit is not a Debye-Waller fall-off the two differ. Which of
        // the two numbers is which is said here rather than left to the reader: deltaB_linear exceeds
        // the deltaB on roughly a third of a real battery, and a verdict line that quoted the larger
        // one without saying so read as a contradiction of the headline above.
        const bool linear = std::isfinite(result.delta_b_linear);
        const double judged = std::max(0.0, linear ? result.delta_b_linear : result.delta_b);
        os << std::setprecision(2)
           << (result.verdict == AnisotropyVerdict::Detected ? "  => DETECTED (" : "  => NOT DETECTED (")
           << result.band << "): ratio " << result.significance << " = "
           << (linear ? "deltaB_linear " : "deltaB ") << judged
           << " A^2 / this dataset's own systematic-error floor of " << result.floor << " A^2.\n";
        if (linear) {
            os << "     The gate is applied to deltaB_linear, the part of the deltaB that follows\n"
               << "     exp(-1/2 s^T B s)";
            if (result.delta_b_linear > result.delta_b)
                os << " - it exceeds the deltaB here because the fitted line passes below\n"
                   << "     the origin";
            else if (result.delta_b_linear < 0.0)
                os << " - the fitted linear slope is negative here, which is no established\n"
                   << "     Debye-Waller component at all rather than a small one, so it enters as zero";
            os << ". The anisotropy to quote is the deltaB of " << result.delta_b << " A^2 above.\n";
        }
    }
    if (std::isfinite(result.detection_limit))
        os << "     Below about " << result.detection_limit << " A^2 nothing could be established on this\n"
           << "     dataset. That limit is set by the systematic error, not by counting, so it does NOT\n"
           << "     improve with more reflections or a longer exposure.\n";
    for (const auto &c : result.cautions)
        os << "     Note: " << c << "\n";
    os << "  Reported and not corrected: no intensity is changed and no reflection is removed.\n";
    return os.str();
}