//! TRACES: FR-CULL-9 | FR-CULL-10
//! Finding the face pairs that could possibly be the same person.
//!
//! [`crate::cluster`] used to begin by computing every pairwise cosine and
//! holding the lot as an `n²` matrix of `f32`. At 1,800 faces that is 13 MB,
//! which is why it survived; at 25,000 it is 2.5 GB, which is why it could not
//! keep surviving.
//!
//! Almost all of that matrix is thrown away unread. Clustering only ever asks
//! whether a pair is *above* the merge threshold, and in a real library the
//! answer is no for well over 99% of pairs — the reference library's 1,813
//! faces produced 7,875 qualifying pairs out of 1.6 million. So this module
//! answers the only question that is actually asked — **which pairs clear the
//! bar** — and returns that sparse list. Memory goes from `O(n²)` to `O(edges)`
//! and the caller never has to hold a matrix at all.
//!
//! # Exact, not approximate
//!
//! The usual way to make this fast is an approximate nearest-neighbour index,
//! which trades recall for speed: it *misses* some true neighbours, and a
//! missed neighbour here is a face that silently never joins its person.
//! Nothing surfaces that — the screen just quietly shows one person as two —
//! so it is a poor trade for a feature whose whole job is to be trusted. This
//! module is exact, and the clusters it produces are identical to those from a
//! full scan.
//!
//! # An exact index was tried, measured, and removed
//!
//! Worth recording so it is not rediscovered as a good idea. The obvious exact
//! index is IVF with a triangle-inequality bound: group the embeddings into
//! cells, and skip a whole cell **pair** when the geometry proves no member of
//! one can reach any member of the other. On the sphere,
//!
//! ```text
//! angle(x, y) >= angle(c_P, c_Q) - radius(P) - radius(Q)
//! ```
//!
//! so a cell pair is impossible when `cos` of that lower bound falls below the
//! threshold. Exact, no recall loss, and it prunes beautifully on synthetic
//! clusters.
//!
//! It prunes **nothing at all** on real face embeddings. Measured over the
//! 1,813-face reference library, at √n = 43 cells:
//!
//! | quantity | measured |
//! |---|---|
//! | median pair angle | 88.5° (cosine 0.026) |
//! | merge threshold | 66.2° (cosine 0.403) |
//! | median cell radius | 80.4° |
//! | median centroid separation | 85.0° |
//! | cell pairs surviving the bound | **946 of 946 — 100%** |
//!
//! The arithmetic is not close. For the bound to exclude a typical cell pair it
//! needs `radius(P) + radius(Q) < 85° - 66° = 19°`, so cells of radius under
//! ~10°. But two photographs of the *same person* sit 36–60° apart, so even a
//! perfect single-identity cell has a radius three times too large. No
//! ball-based partition of this space can have cells tight enough for the
//! inequality to bite — 512-d embeddings are near-orthogonal, and that is the
//! curse of dimensionality doing exactly what it says.
//!
//! So the scan stayed exhaustive, and the effort went where it actually pays:
//! not materialising the matrix, an unrolled dot product, and spreading the
//! blocks across cores. That is `O(n²)` time and `O(edges)` memory, which for
//! this problem is the honest answer.

use crate::calibrate::Calibration;

/// Rows of the similarity triangle handed to one thread at a time.
///
/// Small enough that the tail of the triangle divides evenly across cores —
/// row `i` does `n - i` comparisons, so equal *row counts* are very unequal
/// work — and large enough that the per-block overhead disappears.
const BLOCK: usize = 64;

/// Columns compared against one row block before moving on.
///
/// The other half of the tiling: 64 embeddings of 512 floats is 128 KB, which
/// sits in L2 beside the row block instead of being re-read from memory for
/// every row. See [`scan_rows`] for what it was worth.
const COLUMN_TILE: usize = 64;

/// Below this many faces, do the whole thing on the calling thread.
///
/// Spawning threads for a set this small costs more than the scan.
const THREADS_ABOVE: usize = 2048;

/// One face pair that clears the merge threshold.
///
/// `i < j` always, and the probability is carried because the caller would
/// otherwise recompute the sigmoid it took a dot product to reach.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Pair {
    pub i: usize,
    pub j: usize,
    pub probability: f32,
}

/// What the caller has to tell us about each face.
///
/// Deliberately not [`crate::cluster::Candidate`]: this module has no business
/// knowing what a person or a photograph is, and taking the three arrays it
/// actually reads keeps it testable on bare vectors.
pub struct Faces<'a> {
    /// Every embedding end to end, L2-normalised, [`Faces::dim`] floats each.
    ///
    /// **Flat, not a slice of vectors**, and the difference is measurable: a
    /// `&[Vec<f32>]` is one heap allocation per face and the scan chases a
    /// pointer per row, which defeats both the prefetcher and the tiling this
    /// module does to stay in cache. One buffer is also the layout a GPU would
    /// want, which is where this is eventually going.
    pub embeddings: &'a [f32],
    /// Floats per embedding — [`crate::EMBEDDING_DIM`] in practice, a parameter
    /// so the tests can work in 64 dimensions.
    pub dim: usize,
    /// Source pixels across the aligned crop, for the calibration's size term.
    pub crop_px: &'a [f32],
    /// Which photograph each face came from. Two faces in one frame are not
    /// the same person, so those pairs are never returned (docs/dev/faces.md §9).
    pub images: &'a [u64],
    /// Which faces may be compared *against* — the gallery
    /// ([`crate::embedding::MIN_GALLERY_QUALITY`]).
    ///
    /// A pair needs at least one gallery side: a probe measured against a
    /// reference is a comparison, two short vectors measured against each
    /// other is noise agreeing with noise, and those pairs are never returned.
    /// Filtered here rather than by the caller for the same reason
    /// co-occurrence is: what this module leaves out of the list stays out of
    /// the graph, the components and the merge order, so nothing downstream
    /// has to remember the rule.
    pub gallery: &'a [bool],
}

impl Faces<'_> {
    /// How many faces there are.
    pub fn len(&self) -> usize {
        if self.dim == 0 {
            0
        } else {
            self.embeddings.len() / self.dim
        }
    }

    pub fn is_empty(&self) -> bool {
        self.len() == 0
    }

    #[inline(always)]
    fn row(&self, i: usize) -> &[f32] {
        &self.embeddings[i * self.dim..(i + 1) * self.dim]
    }
}

/// Every pair whose calibrated probability reaches `min_probability`.
///
/// Excludes pairs from the same photograph, which the clusterer would refuse
/// anyway — dropping them here keeps them out of the graph the caller builds
/// and out of the connected components it derives from it.
///
/// Ordered by `(i, j)`, which is what the caller's determinism rests on.
pub fn above_threshold(faces: &Faces, cal: &Calibration, min_probability: f32) -> Vec<Pair> {
    let n = faces.len();
    if n < 2 {
        return Vec::new();
    }

    // The loosest cosine that could clear the bar for *any* pair in the set.
    // Cheaper than the sigmoid by far, and it rejects almost everything.
    let scan = Scan {
        faces,
        cal,
        min_probability,
        tau: loosest_cosine(faces.crop_px, cal, min_probability),
        dot: fastest_dot(),
    };

    let blocks: Vec<(usize, usize)> = (0..n)
        .step_by(BLOCK)
        .map(|start| (start, (start + BLOCK).min(n)))
        .collect();

    if n < THREADS_ABOVE {
        let mut out = Vec::new();
        for &(from, to) in &blocks {
            scan_rows(&scan, from, to, &mut out);
        }
        return out;
    }

    // One worker per core bar one. This runs on a background thread behind a
    // button the user pressed, and NFR-ARCH-2 puts it behind the UI: taking
    // every core would stall the window it is reporting progress to.
    let workers = std::thread::available_parallelism()
        .map(|p| p.get().saturating_sub(1).max(1))
        .unwrap_or(1)
        .min(blocks.len());

    let next = std::sync::atomic::AtomicUsize::new(0);
    let mut parts: Vec<Vec<Vec<Pair>>> = std::thread::scope(|scope| {
        let handles: Vec<_> = (0..workers)
            .map(|_| {
                let next = &next;
                let blocks = &blocks;
                scope.spawn(move || {
                    // Results stay tagged with their block index, so the order
                    // of the output does not depend on which thread got there
                    // first. Determinism is a promise this module keeps.
                    let mut mine: Vec<(usize, Vec<Pair>)> = Vec::new();
                    loop {
                        let b = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
                        let Some(&(from, to)) = blocks.get(b) else {
                            break;
                        };
                        let mut out = Vec::new();
                        scan_rows(&scan, from, to, &mut out);
                        mine.push((b, out));
                    }
                    mine
                })
            })
            .collect();

        let mut slots: Vec<Vec<Vec<Pair>>> = vec![Vec::new(); blocks.len()];
        for h in handles {
            for (b, pairs) in h.join().unwrap_or_default() {
                slots[b].push(pairs);
            }
        }
        slots
    });

    let mut out = Vec::new();
    for slot in &mut parts {
        for pairs in slot.drain(..) {
            out.extend(pairs);
        }
    }
    out
}

/// Everything [`scan_rows`] needs that does not change between blocks.
///
/// A struct rather than eight arguments, and the grouping is real: these five
/// are fixed for a whole scan and only the row range moves.
#[derive(Clone, Copy)]
struct Scan<'a> {
    faces: &'a Faces<'a>,
    cal: &'a Calibration,
    min_probability: f32,
    tau: f32,
    dot: DotFn,
}

/// Compare rows `from..to` against everything after them.
///
/// The upper triangle, split by rows, and **tiled on the column side too**.
/// Walking `j` from `i + 1` to `n` for one row at a time streams the whole
/// embedding array past the core once per row — 336 GB of traffic for an
/// 18,000-face library — where a column tile small enough to sit in L2 is read
/// once per *tile* of rows. Measured on that library, the tiling alone took the
/// scan from 4.64 s to 2.81 s before any change to the kernel.
///
/// Pairs come out ordered by `(i, j)`: the tiles are walked in ascending order
/// but the row loop is inside them, so the block's own output is sorted before
/// it is returned. Blocks are concatenated in row order, so the whole list is
/// ordered — which is the promise the caller's determinism rests on.
fn scan_rows(scan: &Scan, from: usize, to: usize, out: &mut Vec<Pair>) {
    let Scan {
        faces,
        cal,
        min_probability,
        tau,
        dot,
    } = *scan;
    let n = faces.len();
    for tile in (from..n).step_by(COLUMN_TILE) {
        let tile_end = (tile + COLUMN_TILE).min(n);
        for i in from..to {
            // The diagonal: nothing at or before `i` is this row's business.
            let start = tile.max(i + 1);
            if start >= tile_end {
                continue;
            }
            let a = faces.row(i);
            let crop_a = faces.crop_px[i];
            let image_a = faces.images[i];
            let gallery_a = faces.gallery[i];
            for j in start..tile_end {
                if image_a == faces.images[j] || !(gallery_a || faces.gallery[j]) {
                    continue;
                }
                let cos = dot(a, faces.row(j));
                // The cheap rejection, and it takes well over 99% of pairs.
                if cos < tau {
                    continue;
                }
                let probability = cal.probability(cos, crop_a.min(faces.crop_px[j]), 0.0);
                if probability >= min_probability {
                    out.push(Pair { i, j, probability });
                }
            }
        }
    }
    out.sort_unstable_by_key(|p| (p.i, p.j));
}

/// The lowest cosine that could yield `min_probability` for any pair in the set.
///
/// The calibration is `sigmoid(a·cos + b + w_size·log2(crop))`, so for a fixed
/// size term the cosine boundary is exact. The size term is *not* fixed — it
/// varies per pair with the smaller of the two faces — so the safe bound uses
/// whichever face size pushes the boundary lowest: the largest face when
/// `w_size` is positive, the smallest when it is negative.
///
/// Returns [`f32::NEG_INFINITY`] — a filter that rejects nothing — where no
/// boundary exists: a non-positive steepness, for which probability does not
/// increase with cosine, or a threshold at the ends of the sigmoid. Those are
/// degenerate calibrations rather than impossible ones, and the right response
/// is to stop pruning, not to guess.
fn loosest_cosine(crop_px: &[f32], cal: &Calibration, min_probability: f32) -> f32 {
    if cal.a <= 0.0 || !(min_probability > 0.0 && min_probability < 1.0) {
        return f32::NEG_INFINITY;
    }

    let (mut lo, mut hi) = (f32::INFINITY, 0.0_f32);
    for &c in crop_px {
        let c = c.max(1.0);
        lo = lo.min(c);
        hi = hi.max(c);
    }
    if !lo.is_finite() {
        return f32::NEG_INFINITY;
    }

    let extreme = if cal.w_size >= 0.0 { hi } else { lo };
    let tau = cal.boundary_at(min_probability, extreme, 0.0);
    if tau.is_nan() {
        return f32::NEG_INFINITY;
    }
    // Cosines never exceed 1, so a boundary above it legitimately rejects
    // everything. Clamped rather than left free so the comparison stays cheap.
    tau.min(1.0)
}

/// Cosine of two L2-normalised embeddings.
///
/// Eight accumulators rather than one. Floating-point addition is not
/// associative, so the compiler may not re-associate a single running total and
/// the loop serialises on the adder's latency; eight independent chains give it
/// something to pipeline and vectorise. The order is fixed and identical on
/// every run, which is what the caller's determinism needs — it is a different
/// order from the naive sum, not a variable one.
/// A dot product over two equal-length, L2-normalised rows.
///
/// Chosen once per scan rather than per pair — see [`fastest_dot`].
pub(crate) type DotFn = fn(&[f32], &[f32]) -> f32;

/// The widest dot product this machine can actually run.
///
/// # Why this is worth unsafe code
///
/// The scan is the arithmetic floor of the whole subsystem and it was running
/// at **0.7 flops per cycle**. The workspace builds for baseline `x86-64`,
/// which is SSE2 and no FMA, and the portable loop below was not being
/// vectorised into even that. Measured over a real 18,143-face library, on
/// twenty cores:
///
/// | kernel | scan | GFLOP/s |
/// |---|---|---|
/// | portable, untiled (what this replaced) | 4.64 s | 36 |
/// | portable, tiled | 2.81 s | 60 |
/// | **AVX2 + FMA, tiled** | **0.86 s** | **195 |
///
/// Identical pair lists — 1,531,969 — all three ways.
///
/// **Runtime detection on x86-64, unconditional on aarch64.** Advanced SIMD is
/// in the aarch64 baseline, so every Android device that runs this has NEON and
/// there is nothing to detect; on x86-64 AVX2 is not baseline and a binary that
/// assumed it would not start on older hardware.
///
/// The three kernels sum in different orders, so a cosine may differ in its
/// last bit between them. That only matters for a pair sitting exactly on the
/// threshold, and the portable kernel already sums in eight accumulators rather
/// than one, so the module was never bit-comparable with a naive sum.
pub(crate) fn fastest_dot() -> DotFn {
    #[cfg(target_arch = "x86_64")]
    {
        if std::arch::is_x86_feature_detected!("avx2") && std::arch::is_x86_feature_detected!("fma")
        {
            return dot_avx2;
        }
    }
    #[cfg(target_arch = "aarch64")]
    {
        return dot_neon;
    }
    #[allow(unreachable_code)]
    dot
}

/// Eight-wide fused multiply-add, on the half of desktops that have it.
#[cfg(target_arch = "x86_64")]
fn dot_avx2(a: &[f32], b: &[f32]) -> f32 {
    // SAFETY: `fastest_dot` is the only thing that hands this out, and only
    // after `is_x86_feature_detected!` has said both features are present.
    unsafe { dot_avx2_inner(a, b) }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2", enable = "fma")]
unsafe fn dot_avx2_inner(a: &[f32], b: &[f32]) -> f32 {
    use std::arch::x86_64::*;

    let n = a.len().min(b.len());
    let (mut acc0, mut acc1) = (_mm256_setzero_ps(), _mm256_setzero_ps());
    let mut k = 0;
    // Two accumulators, because one FMA cannot start until the previous one
    // retires and the unit is pipelined several deep.
    while k + 16 <= n {
        // SAFETY: `k + 16 <= n`, and `n` is within both slices.
        acc0 = _mm256_fmadd_ps(
            _mm256_loadu_ps(a.as_ptr().add(k)),
            _mm256_loadu_ps(b.as_ptr().add(k)),
            acc0,
        );
        acc1 = _mm256_fmadd_ps(
            _mm256_loadu_ps(a.as_ptr().add(k + 8)),
            _mm256_loadu_ps(b.as_ptr().add(k + 8)),
            acc1,
        );
        k += 16;
    }

    let mut lanes = [0.0_f32; 8];
    _mm256_storeu_ps(lanes.as_mut_ptr(), _mm256_add_ps(acc0, acc1));
    let mut total = ((lanes[0] + lanes[1]) + (lanes[2] + lanes[3]))
        + ((lanes[4] + lanes[5]) + (lanes[6] + lanes[7]));
    for k in k..n {
        total += a[k] * b[k];
    }
    total
}

/// The same, four-wide, for the phone and the tablet.
///
/// No feature detection and no `target_feature`: Advanced SIMD is mandatory in
/// the aarch64 baseline, so this compiles for every Android target the app
/// builds for. The explicit `vfmaq` matters — LLVM will not fuse a multiply and
/// an add on its own without fast-math, which is most of the win.
#[cfg(target_arch = "aarch64")]
fn dot_neon(a: &[f32], b: &[f32]) -> f32 {
    use std::arch::aarch64::*;

    let n = a.len().min(b.len());
    // SAFETY: every load below is bounded by `k + 16 <= n`, and `n` is within
    // both slices. NEON needs no feature detection on aarch64.
    unsafe {
        let mut acc = [vdupq_n_f32(0.0); 4];
        let mut k = 0;
        while k + 16 <= n {
            for (l, slot) in acc.iter_mut().enumerate() {
                *slot = vfmaq_f32(
                    *slot,
                    vld1q_f32(a.as_ptr().add(k + l * 4)),
                    vld1q_f32(b.as_ptr().add(k + l * 4)),
                );
            }
            k += 16;
        }

        let mut total = vaddvq_f32(vaddq_f32(
            vaddq_f32(acc[0], acc[1]),
            vaddq_f32(acc[2], acc[3]),
        ));
        for k in k..n {
            total += a[k] * b[k];
        }
        total
    }
}

/// The portable fallback, and the definition the others have to agree with.
///
/// Eight accumulators so the adds are independent; that is as much as can be
/// asked of a loop that has to compile for anything.
pub(crate) fn dot(a: &[f32], b: &[f32]) -> f32 {
    const LANES: usize = 8;
    let mut acc = [0.0_f32; LANES];
    let n = a.len().min(b.len());
    let chunks = n / LANES;

    for c in 0..chunks {
        let base = c * LANES;
        for (l, slot) in acc.iter_mut().enumerate() {
            *slot += a[base + l] * b[base + l];
        }
    }
    let mut total =
        ((acc[0] + acc[1]) + (acc[2] + acc[3])) + ((acc[4] + acc[5]) + (acc[6] + acc[7]));
    for k in chunks * LANES..n {
        total += a[k] * b[k];
    }
    total
}

#[cfg(test)]
mod tests {
    use super::*;

    const DIM: usize = 64;

    fn cal() -> Calibration {
        Calibration {
            a: 30.0,
            b: -30.0 * 0.35,
            w_size: 0.0,
            valid: true,
            positive_pairs: 1000,
            negative_pairs: 10_000,
        }
    }

    /// A unit vector in a reproducible pseudo-random direction.
    ///
    /// Hashed from the seed rather than drawn from an RNG, so a failure is
    /// reproducible from the test alone.
    fn vector(seed: u64) -> Vec<f32> {
        let mut s = seed.wrapping_mul(0x9E37_79B9_7F4A_7C15) | 1;
        let mut v = Vec::with_capacity(DIM);
        for _ in 0..DIM {
            s ^= s << 13;
            s ^= s >> 7;
            s ^= s << 17;
            v.push(((s >> 11) as f64 / (1u64 << 53) as f64) as f32 - 0.5);
        }
        normalise(v)
    }

    fn normalise(mut v: Vec<f32>) -> Vec<f32> {
        let n = v.iter().map(|x| x * x).sum::<f32>().sqrt();
        for x in &mut v {
            *x /= n;
        }
        v
    }

    /// A vector a known cosine away from `base`.
    fn near(base: &[f32], other: &[f32], cosine: f32) -> Vec<f32> {
        let d = dot(base, other);
        let mut perp: Vec<f32> = other.iter().zip(base).map(|(o, b)| o - d * b).collect();
        let n = perp.iter().map(|x| x * x).sum::<f32>().sqrt();
        for x in &mut perp {
            *x /= n;
        }
        let s = (1.0 - cosine * cosine).max(0.0).sqrt();
        normalise(
            base.iter()
                .zip(&perp)
                .map(|(b, p)| cosine * b + s * p)
                .collect(),
        )
    }

    struct Set {
        embeddings: Vec<f32>,
        crop_px: Vec<f32>,
        images: Vec<u64>,
        gallery: Vec<bool>,
    }

    impl Set {
        fn faces(&self) -> Faces<'_> {
            Faces {
                embeddings: &self.embeddings,
                dim: DIM,
                crop_px: &self.crop_px,
                images: &self.images,
                gallery: &self.gallery,
            }
        }

        fn len(&self) -> usize {
            self.embeddings.len() / DIM
        }
    }

    /// `groups` identities, `per` faces each, every face in its own photograph.
    fn population(groups: usize, per: usize, tightness: f32) -> Set {
        let mut embeddings = Vec::new();
        let mut images = Vec::new();
        let mut image = 0u64;
        for g in 0..groups {
            let base = vector(g as u64 + 1);
            let off = vector(g as u64 + 9_999);
            for m in 0..per {
                embeddings.push(if m == 0 {
                    base.clone()
                } else {
                    near(&base, &off, tightness)
                });
                images.push(image);
                image += 1;
            }
        }
        let crop_px = vec![150.0; embeddings.len()];
        let gallery = vec![true; embeddings.len()];
        Set {
            embeddings: embeddings.concat(),
            crop_px,
            images,
            gallery,
        }
    }

    /// The unpruned, unthreaded, unblocked definition of the answer.
    fn reference(faces: &Faces, cal: &Calibration, min_probability: f32) -> Vec<Pair> {
        let n = faces.len();
        let mut out = Vec::new();
        for i in 0..n {
            for j in i + 1..n {
                if faces.images[i] == faces.images[j] || !(faces.gallery[i] || faces.gallery[j]) {
                    continue;
                }
                let cos: f32 = faces
                    .row(i)
                    .iter()
                    .zip(faces.row(j))
                    .map(|(x, y)| x * y)
                    .sum();
                let p = cal.probability(cos, faces.crop_px[i].min(faces.crop_px[j]), 0.0);
                if p >= min_probability {
                    out.push(Pair {
                        i,
                        j,
                        probability: p,
                    });
                }
            }
        }
        out
    }

    fn same_pairs(a: &[Pair], b: &[Pair]) -> bool {
        a.len() == b.len() && a.iter().zip(b).all(|(x, y)| x.i == y.i && x.j == y.j)
    }

    /// The guard on the SIMD kernels, and the only check the aarch64 one gets
    /// on a machine that is not aarch64: whatever [`fastest_dot`] picked has to
    /// agree with the portable definition. A wrong lane index or a mishandled
    /// tail would show up here as a wildly different number, not a rounding
    /// difference.
    #[test]
    fn the_fastest_kernel_agrees_with_the_portable_one() {
        let fast = fastest_dot();
        for seed in 0..64u64 {
            let a = vector(seed);
            let b = vector(seed + 1_000);
            let (want, got) = (dot(&a, &b), fast(&a, &b));
            assert!(
                (want - got).abs() < 1e-5,
                "kernel disagreed on seed {seed}: {want} vs {got}"
            );
        }

        // A length that is not a multiple of the widest step, so the tail is
        // exercised rather than assumed away.
        let a: Vec<f32> = (0..37).map(|k| k as f32 * 0.01).collect();
        let b: Vec<f32> = (0..37).map(|k| 1.0 - k as f32 * 0.02).collect();
        assert!((dot(&a, &b) - fast(&a, &b)).abs() < 1e-5, "tail mishandled");
    }

    #[test]
    fn nothing_to_pair_is_no_pairs() {
        let s = population(1, 1, 0.9);
        assert!(above_threshold(&s.faces(), &cal(), 0.9).is_empty());
    }

    #[test]
    fn the_blocked_scan_finds_exactly_what_the_reference_does() {
        let s = population(60, 8, 0.97);
        let f = s.faces();
        let got = above_threshold(&f, &cal(), 0.9);
        let want = reference(&f, &cal(), 0.9);
        assert!(!want.is_empty(), "the reference found nothing to check");
        assert!(same_pairs(&got, &want), "{} vs {}", got.len(), want.len());
    }

    /// Past `THREADS_ABOVE` the work is split across cores and stitched back
    /// together, and the stitching is where an order bug would live.
    #[test]
    fn the_threaded_scan_finds_exactly_what_the_reference_does() {
        let s = population(300, 8, 0.97);
        assert!(
            s.len() > THREADS_ABOVE,
            "population is below the threading cutoff"
        );
        let f = s.faces();
        let got = above_threshold(&f, &cal(), 0.9);
        let want = reference(&f, &cal(), 0.9);
        assert!(!want.is_empty());
        assert!(same_pairs(&got, &want), "{} vs {}", got.len(), want.len());
    }

    /// The size term moves the cosine boundary per pair, so the pre-filter has
    /// to be built from the most permissive size in the set or it will drop a
    /// pair that would have qualified.
    #[test]
    fn a_size_weighted_calibration_still_matches_the_reference() {
        let mut s = population(60, 8, 0.97);
        for (i, c) in s.crop_px.iter_mut().enumerate() {
            *c = 40.0 + (i % 17) as f32 * 30.0;
        }
        let sized = Calibration {
            w_size: 0.5,
            b: -30.0 * 0.35 - 0.5 * 7.0,
            ..cal()
        };
        let f = s.faces();
        assert!(same_pairs(
            &above_threshold(&f, &sized, 0.9),
            &reference(&f, &sized, 0.9)
        ));
    }

    /// A negative size weight flips which extreme is permissive. Cheap to get
    /// wrong and silent when it is, so it gets its own case.
    #[test]
    fn a_negative_size_weight_prunes_from_the_other_end() {
        let mut s = population(60, 8, 0.97);
        for (i, c) in s.crop_px.iter_mut().enumerate() {
            *c = 40.0 + (i % 17) as f32 * 30.0;
        }
        let sized = Calibration {
            w_size: -0.5,
            b: -30.0 * 0.35 + 0.5 * 7.0,
            ..cal()
        };
        let f = s.faces();
        assert!(same_pairs(
            &above_threshold(&f, &sized, 0.9),
            &reference(&f, &sized, 0.9)
        ));
    }

    /// A degenerate calibration has no cosine boundary to prune against, and
    /// must stop pruning rather than prune on a bound that does not hold.
    #[test]
    fn a_flat_calibration_prunes_nothing_and_still_agrees() {
        let flat = Calibration {
            a: 0.0,
            b: 4.0,
            ..cal()
        };
        assert_eq!(loosest_cosine(&[150.0], &flat, 0.9), f32::NEG_INFINITY);
        let s = population(20, 6, 0.97);
        let f = s.faces();
        assert!(same_pairs(
            &above_threshold(&f, &flat, 0.9),
            &reference(&f, &flat, 0.9)
        ));
    }

    #[test]
    fn two_faces_in_one_photograph_are_never_paired() {
        let mut s = population(1, 2, 1.0);
        s.images = vec![7, 7];
        assert!(above_threshold(&s.faces(), &cal(), 0.9).is_empty());
    }

    /// A probe against a reference is a comparison; two probes against each
    /// other is not. The rule lives here so that nothing downstream sees the
    /// pair at all.
    #[test]
    fn two_faces_outside_the_gallery_are_never_paired() {
        let mut s = population(1, 3, 1.0);
        s.gallery = vec![false, false, true];
        let pairs = above_threshold(&s.faces(), &cal(), 0.9);
        assert!(
            !pairs.iter().any(|p| p.i == 0 && p.j == 1),
            "two probes were paired with each other"
        );
        // Each probe is still measured against the one reference.
        assert!(pairs.iter().any(|p| p.i == 0 && p.j == 2));
        assert!(pairs.iter().any(|p| p.i == 1 && p.j == 2));
    }

    #[test]
    fn the_gallery_rule_matches_the_reference_at_scale() {
        let mut s = population(60, 8, 0.97);
        for (i, g) in s.gallery.iter_mut().enumerate() {
            *g = i % 3 != 0;
        }
        let f = s.faces();
        let got = above_threshold(&f, &cal(), 0.9);
        let want = reference(&f, &cal(), 0.9);
        assert!(same_pairs(&got, &want), "{} vs {}", got.len(), want.len());
    }

    #[test]
    fn pairs_come_back_in_index_order() {
        let s = population(300, 8, 0.97);
        let pairs = above_threshold(&s.faces(), &cal(), 0.9);
        assert!(pairs
            .windows(2)
            .all(|w| (w[0].i, w[0].j) < (w[1].i, w[1].j)));
        assert!(pairs.iter().all(|p| p.i < p.j));
    }

    /// Threads must not make the answer depend on which one finished first.
    #[test]
    fn the_same_input_yields_the_same_pairs() {
        let s = population(300, 8, 0.97);
        assert_eq!(
            above_threshold(&s.faces(), &cal(), 0.9),
            above_threshold(&s.faces(), &cal(), 0.9)
        );
    }

    /// The unrolled dot has to agree with the obvious one, tail included — the
    /// lengths here are deliberately not multiples of the lane count.
    #[test]
    fn the_unrolled_dot_matches_the_naive_one() {
        for len in [1usize, 7, 8, 9, 63, 64, 65, 512] {
            let a: Vec<f32> = (0..len).map(|i| (i as f32 * 0.37).sin()).collect();
            let b: Vec<f32> = (0..len).map(|i| (i as f32 * 0.11).cos()).collect();
            let naive: f32 = a.iter().zip(&b).map(|(x, y)| x * y).sum();
            assert!(
                (dot(&a, &b) - naive).abs() < 1e-4,
                "len {len}: {} vs {naive}",
                dot(&a, &b)
            );
        }
    }

    /// The pre-filter is the whole speed story, so it is worth asserting it
    /// actually rejects the bulk of the population rather than trusting it to.
    #[test]
    fn the_threshold_filter_rejects_almost_everything() {
        let s = population(60, 8, 0.97);
        let n = s.embeddings.len();
        let total = n * (n - 1) / 2;
        let kept = above_threshold(&s.faces(), &cal(), 0.9).len();
        assert!(
            kept * 20 < total,
            "kept {kept} of {total} pairs, which is not sparse"
        );
    }
}