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

#include <catch2/catch_all.hpp>

#include <cstring>
#include <random>
#include "../image_analysis/bragg_integration/CalcISigma.h"
#include "../image_analysis/scale_merge/ErrorModel.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/scale_merge/SublatticeRelabel.h"
#include "gemmi/reciproc.hpp"

TEST_CASE("HKLKey_NoSG_noMergeFriedel") {
    HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
}

TEST_CASE("HKLKey_NoSG_MergeFriedel") {
    HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
}

TEST_CASE("HKLKey_SG1_MergeFriedel") {
    HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
}

TEST_CASE("HKLKey_SG1_NoMergeFriedel") {
    HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(-1,-2,-3) != hkl_key_gen(1,-2,-3));
}

TEST_CASE("HKLKey_SG96_MergeFriedel") {
    HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(96));
    CHECK(hkl_key_gen(-1, -2, -3) == hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3));
    
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,-1,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2, 1, 3));
}

TEST_CASE("HKLKey_SG96_NoMergeFriedel") {
    HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(96));
    CHECK(hkl_key_gen(-1, -2, -3) != hkl_key_gen(1,2,3));
    CHECK(hkl_key_gen(-1,-2,-3) == hkl_key_gen(-1,-2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2,1,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,-2,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,-1,3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(1,-2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-1,2,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(2,1,-3));
    CHECK(hkl_key_gen(1,2,3) == hkl_key_gen(-2, -1, -3));
    
    CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(-2,-1,3));
    CHECK(hkl_key_gen(1,2,3) != hkl_key_gen(2, 1, 3));
}

TEST_CASE("HKLKey_pack_friedel") {
    HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3).pack() != hkl_key_gen(1,2,3).pack());
    CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack());
    CHECK(hkl_key_gen(-1,-2,-3).pack() != hkl_key_gen(1,-2,-3).pack());
}

TEST_CASE("HKLKey_pack_no_friedel") {
    HKLKeyGenerator hkl_key_gen(true, *gemmi::find_spacegroup_by_number(1));
    CHECK(hkl_key_gen(-1, -2, -3).pack() == hkl_key_gen(1,2,3).pack());
    CHECK(hkl_key_gen(-1,-2,-3).pack() == hkl_key_gen(-1,-2,-3).pack());
    CHECK(hkl_key_gen(-1,-2,-3).pack() != hkl_key_gen(1,-2,-3).pack());
}


TEST_CASE("HKLKey_sys_absence_P212121") {
    HKLKeyGenerator hkl_key_gen(false, *gemmi::find_spacegroup_by_number(19));
    CHECK(hkl_key_gen.IsSystematicallyAbsent(5,0,0));
    CHECK(!hkl_key_gen.IsSystematicallyAbsent(6,0,0));
    CHECK(hkl_key_gen.IsSystematicallyAbsent(0,5,0));
    CHECK(hkl_key_gen.IsSystematicallyAbsent(0,0,5));
    CHECK(!hkl_key_gen.IsSystematicallyAbsent(0,4,0));
    CHECK(!hkl_key_gen.IsSystematicallyAbsent(5,5,5));
}
TEST_CASE("AcceptReflection_ResolutionLimits") {
    Reflection r{};
    r.I = 100.0f;
    r.sigma = 5.0f;
    r.prescaling_corr = 1.0f;
    r.d = 20.0f;

    // No limits: only the finiteness checks apply.
    CHECK(AcceptReflection(r, std::nullopt, std::nullopt));

    // Low-resolution limit rejects anything coarser than the limit, and is exclusive at it.
    CHECK_FALSE(AcceptReflection(r, std::nullopt, std::optional<double>(15.0)));
    CHECK(AcceptReflection(r, std::nullopt, std::optional<double>(20.0)));
    CHECK(AcceptReflection(r, std::nullopt, std::optional<double>(50.0)));

    // High-resolution limit still rejects anything finer, in the same direction as before.
    CHECK_FALSE(AcceptReflection(r, std::optional<double>(25.0), std::nullopt));
    CHECK(AcceptReflection(r, std::optional<double>(2.0), std::optional<double>(50.0)));

    // The plain-double overload treats 0 as "no limit" at both ends.
    CHECK(AcceptReflection(r, 0.0, 0.0));
    CHECK_FALSE(AcceptReflection(r, 0.0, 15.0));
    CHECK(AcceptReflection(r, 2.0, 50.0));

    // A reflection with a saturated pixel is no measurement, whatever its resolution.
    r.overloaded = true;
    CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt));
    CHECK_FALSE(AcceptReflection(r, 0.0, 0.0));
}

// --- Completeness denominator --------------------------------------------------------------------

namespace {
    // A merged set built straight out of the reflections the cell and space group can give, so the
    // test says exactly which of them were measured: every unique reflection between d_min and
    // d_max_measured and none outside. Without Friedel merging an acentric contributes both hands.
    std::vector<MergedReflection> MeasuredBetween(const gemmi::SpaceGroup &sg, const UnitCell &cell,
                                                  double d_min, double d_max_measured,
                                                  bool merge_friedel) {
        const gemmi::UnitCell gemmi_cell = cell;
        const gemmi::GroupOps gops = sg.operations();
        std::vector<MergedReflection> out;
        for (const auto &hkl: gemmi::make_miller_vector(gemmi_cell, &sg, d_min, d_max_measured, true)) {
            MergedReflection r;
            r.h = hkl[0];
            r.k = hkl[1];
            r.l = hkl[2];
            r.d = static_cast<float>(gemmi_cell.calculate_d(hkl));
            r.I = 100.0f;
            r.sigma = 10.0f;
            r.I_half[0] = 100.0f;
            r.I_half[1] = 100.0f;
            out.push_back(r);
            if (!merge_friedel && !gops.is_reflection_centric(hkl)) {
                r.h = -hkl[0];
                r.k = -hkl[1];
                r.l = -hkl[2];
                out.push_back(r);
            }
        }
        return out;
    }

    MergeStatistics StatsWithLowLimit(const gemmi::SpaceGroup &sg, const std::optional<UnitCell> &cell,
                                      const std::vector<MergedReflection> &merged,
                                      std::optional<double> low_limit, bool merge_friedel) {
        DiffractionExperiment x;
        x.SetSpaceGroup(sg);
        ScalingSettings s = x.GetScalingSettings();
        s.LowResolutionLimit_A(low_limit);
        s.MergeFriedel(merge_friedel);
        x.ImportScalingSettings(s);

        MergeOnTheFly merge(x);
        merge.ReferenceCell(cell);
        return merge.MergeStats(merged, {});
    }

    double Completeness(const MergeStatisticsShell &s) {
        return s.possible_unique_reflections > 0
            ? 100.0 * s.unique_reflections / s.possible_unique_reflections : 0.0;
    }

    const gemmi::SpaceGroup &TestSpaceGroup() { return gemmi::get_spacegroup_by_name("P 1 2 1"); }
    constexpr UnitCell TEST_CELL{40, 45, 50, 90, 100, 90};   // synthetic; coarsest reflection ~49 A
}

// The low-resolution terms a beam stop ate must count as missing: the denominator is the declared
// range, so widening the declared range lowers completeness rather than leaving it alone.
TEST_CASE("MergeStats_CompletenessFallsWhenTheLowBoundCrossesAMaskedRegion") {
    const auto &sg = TestSpaceGroup();
    // Nothing coarser than 20 A was measured - it is all behind the stop.
    const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
    REQUIRE(!merged.empty());

    const auto at_20 = StatsWithLowLimit(sg, TEST_CELL, merged, 20.0, true);
    const auto at_50 = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);

    // Declared exactly where the data stop: everything possible was measured.
    CHECK(Completeness(at_20.overall) > 99.0);

    // Declared out to 50 A: the 20-50 A shell is in the denominator and in nothing else.
    CHECK(at_50.overall.possible_unique_reflections > at_20.overall.possible_unique_reflections);
    CHECK(at_50.overall.unique_reflections == at_20.overall.unique_reflections);
    CHECK(Completeness(at_50.overall) < Completeness(at_20.overall));
    // The innermost shell is where it bites.
    CHECK(Completeness(at_50.shells.front()) < Completeness(at_20.shells.front()));
}

// No low-resolution limit means the whole sphere. The cell has no reflection coarser than 50 A, so
// freeing the 50 A limit must count the same set - the fix does not presuppose either default.
TEST_CASE("MergeStats_CompletenessWithNoDeclaredLowLimit") {
    const auto &sg = TestSpaceGroup();
    const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);

    const auto at_50 = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);
    const auto unlimited = StatsWithLowLimit(sg, TEST_CELL, merged, std::nullopt, true);

    CHECK(unlimited.overall.possible_unique_reflections == at_50.overall.possible_unique_reflections);
    CHECK(Completeness(unlimited.overall) == Catch::Approx(Completeness(at_50.overall)));
    // The shell table stays finite even though the bound is not.
    CHECK(std::isfinite(unlimited.shells.front().d_max));
}

// Counting the two Bijvoet mates of an acentric separately doubles the denominator too, so a fully
// measured anomalous set is 100% complete and not 200%.
TEST_CASE("MergeStats_CompletenessNeverExceeds100") {
    const auto &sg = TestSpaceGroup();
    for (const bool merge_friedel: {true, false}) {
        const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 50.0, merge_friedel);
        const auto stats = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, merge_friedel);
        INFO("merge_friedel = " << merge_friedel);
        CHECK(Completeness(stats.overall) <= 100.0);
        CHECK(Completeness(stats.overall) > 99.0);
        for (const auto &sh: stats.shells)
            CHECK(Completeness(sh) <= 100.0);
    }
}

// The shell grid and the denominator share their bounds, so every possible reflection lands in a
// shell: the shells sum to the overall, and the overall is the sphere the run declared.
TEST_CASE("MergeStats_PossibleSumsOverTheShellsToTheDeclaredSphere") {
    const auto &sg = TestSpaceGroup();
    const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);
    const auto stats = StatsWithLowLimit(sg, TEST_CELL, merged, 50.0, true);

    int sum = 0;
    for (const auto &sh: stats.shells)
        sum += sh.possible_unique_reflections;
    CHECK(sum == stats.overall.possible_unique_reflections);

    // Counted independently over the same declared range - nothing is lost between the two.
    const gemmi::UnitCell gemmi_cell = TEST_CELL;
    const int expected = gemmi::count_reflections(gemmi_cell, &sg, stats.overall.d_min * 0.999, 50.0, true);
    CHECK(stats.overall.possible_unique_reflections == expected);
}

// Without a reference cell there is no set to count against; completeness stays unmeasured rather
// than becoming a number, with or without a declared low limit.
TEST_CASE("MergeStats_NoReferenceCellLeavesCompletenessUnmeasured") {
    const auto &sg = TestSpaceGroup();
    const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);

    for (const std::optional<double> low_limit: {std::optional<double>(50.0), std::optional<double>()}) {
        const auto stats = StatsWithLowLimit(sg, std::nullopt, merged, low_limit, true);
        CHECK(stats.overall.possible_unique_reflections == 0);
        CHECK(stats.overall.unique_reflections > 0);
        CHECK(Completeness(stats.overall) == 0.0);
    }
}

// The reference-range table (--report-resolution) is binned over the range it is GIVEN, whether or
// not the data reach it: its last shell ends at the declared d_min, and past the run's own limit
// nothing is read from anywhere - the reflections the run did not keep count as missing.
TEST_CASE("MergeStats_ReferenceRangeIsBinnedOverTheDeclaredRange") {
    const auto &sg = TestSpaceGroup();
    const auto merged = MeasuredBetween(sg, TEST_CELL, 2.0, 20.0, true);   // what the run kept
    DiffractionExperiment x;
    x.SetSpaceGroup(sg);
    MergeOnTheFly merge(x);
    merge.ReferenceCell(TEST_CELL);
    const auto own = merge.MergeStats(merged, {});

    // A range coarser than the run's own: a complete subset, on a grid ending at the declared bound.
    const auto within = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{3.0, 50.0});
    CHECK(within.shells.back().d_min == Catch::Approx(3.0));
    CHECK(within.shells.front().d_max == Catch::Approx(50.0));
    int coarser_than_3 = 0;
    for (const auto &m: merged)
        if (m.d > 3.0f) ++coarser_than_3;
    CHECK(within.overall.unique_reflections == coarser_than_3);
    CHECK(within.overall.unique_reflections < own.overall.unique_reflections);
    CHECK(Completeness(within.overall) > 99.0);

    // A range finer than the run kept: the grid still ends at 1.5 A, the outer shell is empty and its
    // possible reflections are in the denominator, and the measured range says where the data stop.
    const auto beyond = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.5, 50.0});
    CHECK(beyond.shells.back().d_min == Catch::Approx(1.5));
    CHECK(beyond.overall.unique_reflections == own.overall.unique_reflections);
    CHECK(beyond.overall.possible_unique_reflections > own.overall.possible_unique_reflections);
    CHECK(Completeness(beyond.overall) < Completeness(own.overall));
    CHECK(beyond.shells.back().unique_reflections == 0);
    CHECK(beyond.shells.back().possible_unique_reflections > 0);
    CHECK(beyond.overall.d_min == Catch::Approx(own.overall.d_min));

    // A range the data never reach at all is still a grid, not an error.
    const auto empty = merge.MergeStats(merged, {}, {}, std::nullopt, ReportResolutionRange{1.0, 1.5});
    CHECK(empty.overall.unique_reflections == 0);
    CHECK(empty.overall.possible_unique_reflections > 0);
    CHECK(empty.overall.d_min == 0.0f);
}

// ---------------------------------------------------------------- the automatic resolution cutoff
namespace {
    // Half-set pairs spread uniformly in s = 1/d^2 over [s_from, s_to), either correlated with each
    // other (signal) or drawn independently (noise, CC1/2 ~ 0), so a whole CC1/2 curve can be built
    // band by band.
    void AddBand(std::vector<MergedReflection> &v, std::mt19937 &rng,
                 double s_from, double s_to, int n, bool correlated, double sigma = 1.0) {
        std::normal_distribution<double> g(0.0, sigma);
        for (int j = 0; j < n; ++j) {
            MergedReflection m;
            m.h = static_cast<int32_t>(v.size());   // every reflection its own index
            m.d = static_cast<float>(1.0 / std::sqrt(s_from + (j + 0.5) * (s_to - s_from) / n));
            const double a = g(rng), b = g(rng);
            m.I_half[0] = static_cast<float>(a);
            m.I_half[1] = static_cast<float>(correlated ? a : b);
            // The pair's weight in a CC1/2: the precision it would have had at the typical frame scale
            // over the precision it has - 1/sigma^2 for a band scaled up by sigma from dead frames.
            m.cc_weight = static_cast<float>(1.0 / (sigma * sigma));
            v.push_back(m);
        }
    }
}

// A clean fall-off: CC1/2 crosses the target where the signal stops, and the cut is written one
// shell past it.
TEST_CASE("ResolutionCutoff_CleanFallOff") {
    Logger logger("test");
    std::mt19937 rng(12345);
    std::vector<MergedReflection> merged;
    AddBand(merged, rng, 0.01, 0.25, 480, true);    // signal to 1/sqrt(0.25) = 2.00 A
    AddBand(merged, rng, 0.25, 0.51, 520, false);   // noise beyond it

    const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
    REQUIRE(rc.d_fit);
    CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
    REQUIRE(rc.d_cut);
    CHECK(*rc.d_cut < *rc.d_fit);                   // the deliberate one-shell extension
    CHECK(*rc.d_cut == Catch::Approx(1.92).margin(0.15));
}

// Past the fall-off, a few single reflections read hundreds of times the shell's noise in both halves -
// an artefact the half-sets reproduce. Each carries nearly all of its shell's covariance, so a Pearson
// CC1/2 reads those shells as signal and the cut runs on through them; one reflection per shell must not
// decide where the signal stops.
TEST_CASE("ResolutionCutoff_SingleReflectionsDoNotHoldTheCutOpen") {
    Logger logger("test");
    std::mt19937 rng(12345);
    std::vector<MergedReflection> merged;
    AddBand(merged, rng, 0.01, 0.25, 480, true);
    AddBand(merged, rng, 0.25, 0.51, 520, false);
    for (const double s : {0.26, 0.28, 0.30, 0.32}) {
        MergedReflection m;
        m.h = static_cast<int32_t>(merged.size());
        m.d = static_cast<float>(1.0 / std::sqrt(s));
        m.I_half[0] = m.I_half[1] = 300.0f;
        merged.push_back(m);
    }

    const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
    REQUIRE(rc.d_fit);
    CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));
}

// A sweep with a long stretch where the crystal barely diffracted: the reflections measured only there
// are noise scaled up by 1/G, sigmas with them, and there are many of them at every resolution. Counted
// by the information they carry they must not hide where the well-measured reflections stop.
TEST_CASE("ResolutionCutoff_ScaledUpNoiseDoesNotHideTheFallOff") {
    Logger logger("test");
    std::mt19937 rng(12345);
    std::vector<MergedReflection> merged;
    AddBand(merged, rng, 0.01, 0.25, 480, true);
    AddBand(merged, rng, 0.25, 0.51, 520, false);
    AddBand(merged, rng, 0.01, 0.51, 300, false, 50.0);

    const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
    REQUIRE(rc.d_fit);
    CHECK(*rc.d_fit == Catch::Approx(2.0).margin(0.15));

    // Counted as equals, the same reflections hide it: the curve reads noise from the first bin on.
    for (auto &m : merged) m.cc_weight = 1.0f;
    const auto unweighted = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
    CHECK((!unweighted.d_fit || *unweighted.d_fit > 3.0));
}

// The frame factor of the CC1/2 weight is 1 on a sweep whose frames all sit at one scale, however that
// scale is spread over the frames, and 1 for any frame brighter than typical.
TEST_CASE("CCHalfFrameFactors_UniformScaleIsOne") {
    const auto f = CCHalfFrameFactors({1.7, 1.7, 1.7, 1.7}, {100, 3, 250, 0});
    for (double x : f) CHECK(x == 1.0);
    const auto g = CCHalfFrameFactors({1.0, 2.0}, {100, 100});
    CHECK(g[1] == 1.0);
    CHECK(g[0] == Catch::Approx(9.0 / 5.0 * 9.0 / 5.0));   // G_ref = (1 + 8) / (1 + 4)
}

// Most of the sweep at 2% of the good frames' scale: the typical scale is still the good frames', not
// the run median, and an observation from a dead frame carries 2500x the variance of a good one.
TEST_CASE("CCHalfFrameFactors_DeadStretchDoesNotSetTheTypicalScale") {
    std::vector<double> scale;
    std::vector<int64_t> n;
    for (int f = 0; f < 100; ++f) {
        scale.push_back(f < 60 ? 0.02 : 1.0);
        n.push_back(50);
    }
    scale.push_back(NAN);   // a frame without a scale counts as G = 1
    n.push_back(50);
    const auto factor = CCHalfFrameFactors(scale, n);
    CHECK(factor[0] == Catch::Approx(2500.0).epsilon(0.01));
    CHECK(factor[99] == Catch::Approx(1.0).epsilon(1e-3));
    CHECK(factor[100] == Catch::Approx(1.0).epsilon(1e-3));
}

// Equal weights give the ordinary R_meas terms, whatever their common value.
TEST_CASE("WeightedRmeas_EqualWeightsAreTheOrdinaryRmeas") {
    // Three observations 9, 10, 11 of a reflection with <I> = 10: sum|dev| = 2, sum I = 30.
    double num, den;
    REQUIRE(WeightedRmeasTerms(0.5 * 2.0, 0.5 * 30.0, 0.5 * 3, 0.25 * 3, num, den));
    CHECK(num == Catch::Approx(std::sqrt(3.0 / 2.0) * 2.0));
    CHECK(den == Catch::Approx(30.0));
}

// An observation that carries almost no information counts for almost nothing: the terms tend to those
// of the reflection without it, and a reflection left with one effective observation has none.
TEST_CASE("WeightedRmeas_NegligibleWeightDropsOut") {
    // 9 and 11 at weight 1, <I> = 10; a third at 40 with weight 1e-6.
    const double v = 1e-6;
    double num, den;
    REQUIRE(WeightedRmeasTerms(2.0 + v * 30.0, 20.0 + v * 40.0, 2.0 + v, 2.0 + v * v, num, den));
    CHECK(num == Catch::Approx(std::sqrt(2.0) * 2.0).epsilon(1e-4));
    CHECK(den == Catch::Approx(20.0).epsilon(1e-4));
    CHECK_FALSE(WeightedRmeasTerms(0.0, 10.0, 1.0, 1.0, num, den));
}

// A fall-off region a logistic cannot follow: CC1/2 drops through the target and comes straight back
// up. The fitted crossing is then an extrapolation far past the bins it was made over, and reading
// the cut off it writes the data deep into the noise; the crossing the bins themselves show is where
// the signal stopped, and that is what must be used.
TEST_CASE("ResolutionCutoff_RaggedFallOffIsReadOffTheBins") {
    Logger logger("test");
    std::mt19937 rng(12345);
    std::vector<MergedReflection> merged;
    AddBand(merged, rng, 0.01, 0.13, 240, true);    // signal to 1/sqrt(0.13) = 2.77 A
    AddBand(merged, rng, 0.13, 0.17, 80, false);    // a hole below the target
    AddBand(merged, rng, 0.17, 0.25, 160, true);    // correlated again - not a fall-off
    AddBand(merged, rng, 0.25, 0.51, 520, false);

    const auto rc = ComputeCCHalfLogisticCutoff(merged, 0.30, logger);
    REQUIRE(rc.d_fit);
    CHECK(*rc.d_fit == Catch::Approx(2.77).margin(0.20));
    REQUIRE(rc.d_cut);
    CHECK(*rc.d_cut > 2.30);                        // coarser than the band that correlates again
}

namespace {
    // Samples of var = a*s2 + b^2*I2 over four decades of counting I/sigma, with counting variances
    // spread over a factor of 100 at every intensity, and a fraction of gross outliers.
    std::vector<ErrorModelSample> SyntheticErrorModelSamples(double a, double b, double max_snr,
                                                             double outlier_fraction) {
        std::mt19937 rng(7);
        std::uniform_real_distribution<double> u(0.0, 1.0);
        std::normal_distribution<double> z(0.0, 1.0);
        std::vector<ErrorModelSample> out;
        for (int i = 0; i < 200000; ++i) {
            const double s2 = std::pow(10.0, 2.0 * u(rng));
            const double snr = max_snr * std::pow(10.0, -4.0 * u(rng));
            const double I2 = snr * snr * s2;
            const double e = z(rng);
            double dev2 = (a * s2 + b * b * I2) * e * e;
            if (u(rng) < outlier_fraction) dev2 *= 400.0;
            out.push_back({s2, I2, dev2, 2.0f});
        }
        return out;
    }
}

// The fit recovers a and b from observations whose counting variances differ widely inside every
// intensity bin (where a ratio of bin medians does not), with 0.3% gross outliers in the pool.
TEST_CASE("ErrorModel_RecoversAAndBThroughOutliers") {
    const auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003);
    std::vector<ErrorModelBinned> scratch;
    const auto fit = FitErrorModel(pool, scratch, 4);
    REQUIRE(fit.active);
    CHECK(fit.b_measured);
    CHECK(fit.b_resolved);
    CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03));
    CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05));
}

// A few observations whose counting variance is hugely overstated - their scatter is ordinary - do not
// set the fit: each sample counts by its deviation relative to its own variance.
TEST_CASE("ErrorModel_OverstatedCountingVarianceDoesNotSetTheFit") {
    auto pool = SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.0);
    for (size_t i = 0; i < pool.size(); i += 200)
        pool[i].s2 *= 1e4;
    std::vector<ErrorModelBinned> scratch;
    const auto fit = FitErrorModel(pool, scratch, 4);
    REQUIRE(fit.active);
    CHECK(fit.a == Catch::Approx(1.3).epsilon(0.03));
    CHECK(1.0 / std::sqrt(fit.b2) == Catch::Approx(1.0 / 0.03).epsilon(0.05));
}

// The same pool on one thread and on many gives the same numbers.
TEST_CASE("ErrorModel_SameOnAnyNumberOfThreads") {
    const auto pool = SyntheticErrorModelSamples(1.1, 0.05, 100.0, 0.001);
    std::vector<ErrorModelBinned> scratch;
    const auto one = FitErrorModel(pool, scratch, 1);
    const auto many = FitErrorModel(pool, scratch, 16);
    CHECK(one.a == many.a);
    CHECK(one.b2 == many.b2);
}

// Data that never reach where b could be seen: a is fitted alone and b held at 0.
TEST_CASE("ErrorModel_BNotMeasuredOnWeakData") {
    const auto pool = SyntheticErrorModelSamples(0.9, 0.05, 1.5, 0.0);
    std::vector<ErrorModelBinned> scratch;
    const auto fit = FitErrorModel(pool, scratch, 4);
    REQUIRE(fit.active);
    CHECK(!fit.b_measured);
    CHECK(fit.b2 == 0.0);
    CHECK(fit.a == Catch::Approx(0.9).epsilon(0.03));
}

#ifdef JFJOCH_USE_CUDA
#include "../common/CUDAWrapper.h"
#include "../image_analysis/scale_merge/ErrorModelGPU.h"

// The device fit against the host one. Same bins and same medians; only the order each bin's two sums
// are added in differs, so a and b^2 agree to rounding and every flag agrees exactly. The same pool
// twice on the device gives the same bits.
TEST_CASE("ErrorModel_DeviceFitIsTheHostFit") {
    if (get_gpu_count() == 0) {
        WARN("No CUDA GPU present. Skipping ErrorModel_DeviceFitIsTheHostFit");
        return;
    }
    for (const auto &pool : {SyntheticErrorModelSamples(1.3, 0.03, 300.0, 0.003),
                             SyntheticErrorModelSamples(0.9, 0.05, 1.5, 0.0)}) {
        std::vector<ErrorModelBinned> scratch;
        const auto host = FitErrorModel(pool, scratch, 4);
        const auto dev = FitErrorModelGPU(pool);
        const auto again = FitErrorModelGPU(pool);
        CHECK(dev.active == host.active);
        CHECK(dev.b_measured == host.b_measured);
        CHECK(dev.b_resolved == host.b_resolved);
        CHECK(dev.a == Catch::Approx(host.a).epsilon(1e-12));
        CHECK(dev.b2 == Catch::Approx(host.b2).epsilon(1e-12).margin(1e-300));
        CHECK(again.a == dev.a);
        CHECK(again.b2 == dev.b2);
    }
}
#endif

// The dataset-wide Wilson B: a pure exp(-B / 2d^2) fall-off, read back over the fit window. A crystal
// diffracting to just under 4 A has almost nothing between its d_min and the 4 A convention, so the
// window is widened to the 1/d^2 width a 3 A dataset gets; and a high B is a measurement, not garbage.
TEST_CASE("Global Wilson B over a widened window, with no upper cap", "[wilson]") {
    const auto merge = [](double b, double d_min) {
        std::vector<MergedReflection> out;
        for (int i = 0; i < 20000; ++i) {
            // Uniform in 1/d^3 from 20 A to d_min, as reflections are.
            const double s3 = (i + 0.5) / 20000.0 * (1.0 / std::pow(d_min, 3) - 1.0 / 8000.0) + 1.0 / 8000.0;
            MergedReflection r;
            r.h = i; r.k = 0; r.l = 0;
            r.d = static_cast<float>(1.0 / std::cbrt(s3));
            r.I = static_cast<float>(1e6 * std::exp(-b / (2.0 * r.d * r.d)));
            r.sigma = 0.01f * r.I;
            out.push_back(r);
        }
        return out;
    };
    const GlobalWilsonB high_res = CalcGlobalWilsonB(merge(40.0, 2.0));
    CHECK(high_res.b == Catch::Approx(40.0).epsilon(0.02));
    CHECK(high_res.d_low_A == Catch::Approx(4.0));

    const GlobalWilsonB low_res = CalcGlobalWilsonB(merge(300.0, 3.9));
    CHECK(low_res.b == Catch::Approx(300.0).epsilon(0.02));
    CHECK(low_res.n_shells >= 3);
    CHECK(1.0 / (3.9 * 3.9) - 1.0 / (low_res.d_low_A * low_res.d_low_A)
          == Catch::Approx(1.0 / 9.0 - 1.0 / 16.0).epsilon(0.01));
}

// The parallel AddImages is the serial AddImage loop, bit for bit: the per-observation work is spread
// over the workers, but every sum and every half-set choice is still made in image order. More images
// than one of its blocks, so the block boundary is crossed.
TEST_CASE("MergeOnTheFly_AddImagesIsTheSerialMerge") {
    DiffractionExperiment x;
    x.SetSpaceGroup(TestSpaceGroup());

    std::mt19937 rng(7);
    std::uniform_int_distribution<int> index(-6, 6);
    std::uniform_real_distribution<float> unit(0.0f, 1.0f);
    std::vector<IntegrationOutcome> outcomes(5000);
    for (auto &o: outcomes) {
        for (int j = 0; j < 20; ++j) {
            Reflection r{};
            r.h = index(rng);
            r.k = index(rng);
            r.l = index(rng);
            if (r.h == 0 && r.k == 0 && r.l == 0)
                continue;
            r.d = 20.0f / std::sqrt(static_cast<float>(r.h * r.h + r.k * r.k + r.l * r.l));
            r.I = 100.0f * unit(rng) - 5.0f;
            r.sigma = 1.0f + 5.0f * unit(rng);
            r.var_bkg = r.sigma * r.sigma * unit(rng);
            r.prescaling_corr = 1.0f;
            r.partiality = 0.5f + 0.5f * unit(rng);
            r.image_scale_corr = 0.5f + unit(rng);
            r.on_ice_ring = unit(rng) < 0.05f;
            o.reflections.push_back(r);
        }
    }

    MergeOnTheFly serial(x);
    serial.RefineErrorModel(outcomes);
    for (size_t i = 0; i < outcomes.size(); ++i)
        serial.AddImage(outcomes[i], static_cast<int64_t>(i));
    const auto serial_merged = serial.ExportReflections();

    MergeOnTheFly parallel(x);
    parallel.Threads(8);
    parallel.RefineErrorModel(outcomes);
    parallel.AddImages(outcomes);
    const auto parallel_merged = parallel.ExportReflections();

    REQUIRE(serial.ErrorModelActive());
    CHECK(serial.ErrorModelA() == parallel.ErrorModelA());
    CHECK(serial.ErrorModelB() == parallel.ErrorModelB());
    REQUIRE(serial_merged.size() == parallel_merged.size());
    for (size_t i = 0; i < serial_merged.size(); ++i) {
        CHECK(serial_merged[i].h == parallel_merged[i].h);
        CHECK(serial_merged[i].k == parallel_merged[i].k);
        CHECK(serial_merged[i].l == parallel_merged[i].l);
        CHECK(serial_merged[i].I == parallel_merged[i].I);
        CHECK(serial_merged[i].sigma == parallel_merged[i].sigma);
        for (int half = 0; half < 2; ++half)   // NaN where a half is empty, on both sides
            CHECK(std::memcmp(&serial_merged[i].I_half[half], &parallel_merged[i].I_half[half], sizeof(float)) == 0);
    }

    const auto serial_stats = serial.MergeStats(serial_merged, outcomes);
    const auto parallel_stats = parallel.MergeStats(parallel_merged, outcomes);
    CHECK(serial_stats.overall.total_observations == parallel_stats.overall.total_observations);
    CHECK(serial_stats.overall.cc_half == parallel_stats.overall.cc_half);
    CHECK(serial_stats.overall.r_meas == parallel_stats.overall.r_meas);
}

// Stills partiality post-refinement on noiseless synthetic crystals: each crystal's true orientation is
// tilted from the stored one by up to 0.05 deg, its intensities are the true partials of a fixed set of
// full intensities, and its stored partialities those of the stored orientation. The refinement has to
// move the partialities towards the true ones. (The Ceres solve this replaced gave the same numbers: error
// 0.0559 -> 0.0302, mean tilt 0.0206 deg.)
TEST_CASE("StillsPartialityRefine_RecoversTheTilt") {
    constexpr float CELL = 50.0f, WAVELENGTH = 1.0f, D_MIN = 2.0f;
    constexpr double GAMMA0 = 0.001, ACCEPT = 1.7 * GAMMA0;   // RLP width; |dist_ewald| accepted at prediction
    DiffractionExperiment x;
    x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
    DiffractionGeometry geom;
    geom.Wavelength_A(WAVELENGTH);
    const Coord S0 = geom.GetScatteringVector();

    std::mt19937 rng(11);
    std::normal_distribution<double> gauss(0.0, 1.0);
    std::uniform_real_distribution<double> unit(0.0, 1.0);
    auto full_I = [](int h, int k, int l) {   // a fixed, Friedel-symmetric pseudo-random full intensity
        const uint64_t key = HKLKeyGenerator(true, gemmi::get_spacegroup_by_name("P 1"))(h, k, l).pack();
        return 10.0 + 990.0 * static_cast<double>((key * 0x9e3779b97f4a7c15ULL) >> 40) / static_cast<double>(1ULL << 24);
    };
    auto rotate = [](const double R[9], const Coord &v) {
        return Coord(static_cast<float>(R[0] * v.x + R[1] * v.y + R[2] * v.z),
                     static_cast<float>(R[3] * v.x + R[4] * v.y + R[5] * v.z),
                     static_cast<float>(R[6] * v.x + R[7] * v.y + R[8] * v.z));
    };
    auto rotation = [](double ax, double ay, double az, double angle, double R[9]) {   // Rodrigues
        const double n = std::sqrt(ax * ax + ay * ay + az * az);
        ax /= n; ay /= n; az /= n;
        const double c = std::cos(angle), s = std::sin(angle), t = 1.0 - c;
        const double M[9] = {c + t * ax * ax, t * ax * ay - s * az, t * ax * az + s * ay,
                             t * ax * ay + s * az, c + t * ay * ay, t * ay * az - s * ax,
                             t * ax * az - s * ay, t * ay * az + s * ax, c + t * az * az};
        std::copy(M, M + 9, R);
    };
    auto partiality = [&](const Coord &q) {
        const double de = (q + S0).Length() - 1.0 / WAVELENGTH;
        return std::exp(-0.5 * de * de / (GAMMA0 * GAMMA0));
    };

    std::vector<IntegrationOutcome> outcomes(150);
    std::vector<std::vector<double>> true_p(outcomes.size());
    const int hmax = static_cast<int>(CELL / D_MIN);
    for (size_t c = 0; c < outcomes.size(); ++c) {
        double R[9], T[9];
        rotation(gauss(rng), gauss(rng), gauss(rng), 2.0 * M_PI * unit(rng), R);
        const double tilt = (0.02 + 0.03 * unit(rng)) * M_PI / 180.0, phi = 2.0 * M_PI * unit(rng);
        rotation(std::cos(phi), std::sin(phi), 0.0, tilt, T);   // about an axis normal to the beam
        auto &o = outcomes[c];
        o.geom = geom;
        o.latt = CrystalLattice(rotate(R, Coord(CELL, 0, 0)), rotate(R, Coord(0, CELL, 0)),
                                rotate(R, Coord(0, 0, CELL)));
        const double G = 0.5 + unit(rng);
        for (int h = -hmax; h <= hmax; ++h)
            for (int k = -hmax; k <= hmax; ++k)
                for (int l = -hmax; l <= hmax; ++l) {
                    const Coord q = o.latt.Astar() * h + o.latt.Bstar() * k + o.latt.Cstar() * l;
                    const float dstar = q.Length();
                    if (dstar == 0.0f || dstar > 1.0f / D_MIN)
                        continue;
                    if (std::fabs((q + S0).Length() - 1.0 / WAVELENGTH) > ACCEPT)
                        continue;
                    const double p_true = partiality(rotate(T, q));
                    Reflection r{};
                    r.h = h; r.k = k; r.l = l;
                    r.d = 1.0f / dstar;
                    r.I = static_cast<float>(G * p_true * full_I(h, k, l));
                    r.sigma = static_cast<float>(1.0 + 0.01 * r.I);
                    r.prescaling_corr = 1.0f;
                    r.partiality = static_cast<float>(partiality(q));
                    r.image_scale_corr = 1.0f / r.partiality;
                    o.reflections.push_back(r);
                    true_p[c].push_back(p_true);
                }
        REQUIRE(o.reflections.size() >= 40);
    }

    auto mean_abs_error = [&]() {
        double sum = 0.0;
        size_t n = 0;
        for (size_t c = 0; c < outcomes.size(); ++c)
            for (size_t i = 0; i < outcomes[c].reflections.size(); ++i, ++n)
                sum += std::fabs(outcomes[c].reflections[i].partiality - true_p[c][i]);
        return sum / static_cast<double>(n);
    };
    const double before = mean_abs_error();
    ScaleOnTheFly(x, MergeAll(x, outcomes, 4, std::nullopt)).Scale(outcomes, 4);
    const double mean_tilt = StillsPartialityRefine(x).Run(outcomes, 4, std::nullopt);
    const double after = mean_abs_error();
    CHECK(mean_tilt > 0.015);
    CHECK(mean_tilt < 0.06);
    CHECK(after < 0.6 * before);   // measured 0.54: the prior holds the larger tilts back
}

namespace {
    // The basis vectors a, b, c of a cell, turned by a fixed arbitrary rotation, as a frame's lattice.
    std::array<Coord, 3> TurnedAxes(const UnitCell &cell) {
        const CrystalLattice l(cell);
        const double c = std::cos(0.7), s = std::sin(0.7);   // about (1, 2, 3)/sqrt(14)
        const double n[3] = {1.0 / std::sqrt(14.0), 2.0 / std::sqrt(14.0), 3.0 / std::sqrt(14.0)};
        auto turn = [&](const Coord &v) {
            const double dot = n[0] * v.x + n[1] * v.y + n[2] * v.z;
            const double x[3] = {n[1] * v.z - n[2] * v.y, n[2] * v.x - n[0] * v.z, n[0] * v.y - n[1] * v.x};
            return Coord(static_cast<float>(v.x * c + x[0] * s + n[0] * dot * (1 - c)),
                         static_cast<float>(v.y * c + x[1] * s + n[1] * dot * (1 - c)),
                         static_cast<float>(v.z * c + x[2] * s + n[2] * dot * (1 - c)));
        };
        return {turn(l.Vec0()), turn(l.Vec1()), turn(l.Vec2())};
    }

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

    IntegrationOutcome FrameWithAllIndices(const CrystalLattice &latt, int hmax) {
        IntegrationOutcome o;
        o.latt = latt;
        for (int h = -hmax; h <= hmax; ++h)
            for (int k = -hmax; k <= hmax; ++k)
                for (int l = -hmax; l <= hmax; ++l) {
                    Reflection r{};
                    r.h = h; r.k = k; r.l = l;
                    r.predicted_x = static_cast<float>(o.reflections.size());   // identifies the observation
                    r.I = 100.0f;
                    r.sigma = 10.0f;
                    const float q = ScatteringVector(latt, r).Length();
                    r.d = q > 0.0f ? 1.0f / q : 100.0f;
                    o.reflections.push_back(r);
                }
        return o;
    }
}

// A frame indexed with its c axis doubled: its reciprocal lattice holds every node of the true one plus
// the planes in between. Relabelling keeps exactly the true nodes, each at the same scattering vector,
// and gives the frame the true cell.
TEST_CASE("SublatticeRelabel_DoubledAxisGoesBackToTheTrueCell") {
    const UnitCell reference{60, 60, 25, 90, 90, 90};   // synthetic tetragonal
    const auto ax = TurnedAxes(reference);
    const CrystalLattice super(ax[0], ax[1], ax[2] * 2.0f);
    // A frame on the true cell beside it sets how fine the relabelled reflections may go.
    std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(super, 4),
                                             FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
    const auto before = outcomes[0];

    IndexingSettings settings;
    const auto result = RelabelOnSublattice(outcomes, reference, settings);
    CHECK(result.frames == 1);
    CHECK(result.reflections_kept == 9 * 9 * 5);       // l even: -4, -2, 0, 2, 4
    CHECK(result.reflections_dropped == 9 * 9 * 4);
    CHECK(outcomes[0].latt.GetUnitCell().is_close(reference, 0.001, 0.01));

    for (const auto &r: outcomes[0].reflections) {
        const auto it = std::find_if(before.reflections.begin(), before.reflections.end(),
                                     [&](const Reflection &b) { return b.predicted_x == r.predicted_x; });
        REQUIRE(it != before.reflections.end());
        CHECK(it->l % 2 == 0);
        const Coord q_before = ScatteringVector(before.latt, *it);
        const Coord q_after = ScatteringVector(outcomes[0].latt, r);
        CHECK((q_after - q_before).Length() < 1e-5f);
    }
}

// The true lattice in another setting (axes permuted) is the same lattice: every reflection is kept and
// only re-indexed. A frame whose lattice does not contain the reference is left alone.
TEST_CASE("SublatticeRelabel_OtherSettingAndUnrelatedLattice") {
    const UnitCell reference{40, 50, 60, 90, 90, 90};   // synthetic orthorhombic
    const auto ax = TurnedAxes(reference);
    std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(CrystalLattice(ax[1], ax[2], ax[0]), 3),
                                             FrameWithAllIndices(CrystalLattice(UnitCell{37, 53, 71, 90, 90, 90}), 3),
                                             FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
    const auto permuted = outcomes[0];
    const auto unrelated = outcomes[1];

    IndexingSettings settings;
    const auto result = RelabelOnSublattice(outcomes, reference, settings);
    CHECK(result.frames == 1);
    CHECK(result.reflections_dropped == 0);
    CHECK(outcomes[0].latt.GetUnitCell().is_close(reference, 0.001, 0.01));
    for (size_t i = 0; i < outcomes[0].reflections.size(); ++i) {
        const Coord q_before = ScatteringVector(permuted.latt, permuted.reflections[i]);
        const Coord q_after = ScatteringVector(outcomes[0].latt, outcomes[0].reflections[i]);
        CHECK((q_after - q_before).Length() < 1e-5f);
    }
    CHECK(outcomes[1].reflections.size() == unrelated.reflections.size());
    CHECK(outcomes[1].reflections[5].h == unrelated.reflections[5].h);
}

// A doubled axis that makes the supercell metric cubic fits the reference along any of the three axes;
// only the true one has signal at its nodes, and that is the one taken.
TEST_CASE("SublatticeRelabel_SignalPicksTheAxisThatWasDoubled") {
    const UnitCell reference{60, 60, 30, 90, 90, 90};   // synthetic tetragonal; doubled c gives a cube
    const auto ax = TurnedAxes(reference);
    std::vector<IntegrationOutcome> outcomes{FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2] * 2.0f), 4),
                                             FrameWithAllIndices(CrystalLattice(ax[0], ax[1], ax[2]), 4)};
    for (auto &r: outcomes[0].reflections)
        r.I = r.l % 2 == 0 ? 100.0f : 0.0f;   // the nodes between the true ones are empty
    const auto before = outcomes[0];

    IndexingSettings settings;
    const auto result = RelabelOnSublattice(outcomes, reference, settings);
    REQUIRE(result.frames == 1);
    for (const auto &r: outcomes[0].reflections) {
        const auto &b = before.reflections[static_cast<size_t>(r.predicted_x)];
        CHECK(b.l % 2 == 0);
        CHECK((ScatteringVector(outcomes[0].latt, r) - ScatteringVector(before.latt, b)).Length() < 1e-5f);
    }
}