lucisearch 0.8.0

//! Distance kernel property tests against a private reference scalar.
//!
//! The reshaped kernels in [`super`] target LLVM auto-vectorization and
//! diverge structurally from the naive "iterator chain + sum" form. The
//! reference scalar kept here is the indisputable correct loop; production
//! kernels must agree with it within ULP tolerance across the input space.
//!
//! Property-test inputs are deterministic (LCG-seeded) so failures
//! reproduce; magnitudes and dimensions cover the SIMD-lane-width
//! boundaries called out in
//! [[optimization-vector-distance-kernel-trait]] §"Property tests".
//!
//! See [[cross-platform-portability]] for the related cross-platform
//! agreement tests (ULP tolerance across writer/reader platforms);
//! those live in `luci-index/tests/cross_platform_portability.rs`
//! since they need the full storage pipeline.

#![cfg(test)]

use super::{DistanceMetric, distance, normalize_in_place};

/// Reference scalar dot product — the indisputably-correct loop.
///
/// Intentionally minimal: no `chunks_exact`, no parallel accumulators,
/// no `target_feature` attributes. The production [`super::dot_product`]
/// is verified against this within ULP tolerance.
fn reference_dot(a: &[f32], b: &[f32]) -> f32 {
    assert_eq!(a.len(), b.len());
    let mut sum = 0.0f32;
    for i in 0..a.len() {
        sum += a[i] * b[i];
    }
    sum
}

/// Reference scalar L2 (squared then sqrt). Same minimal form.
fn reference_l2(a: &[f32], b: &[f32]) -> f32 {
    assert_eq!(a.len(), b.len());
    let mut sum = 0.0f32;
    for i in 0..a.len() {
        let d = a[i] - b[i];
        sum += d * d;
    }
    sum.sqrt()
}

/// Simple LCG for deterministic test inputs.
///
/// `state * 6364136223846793005 + 1442695040888963407` — Knuth's
/// 64-bit constants. The `(value >> 33)` is the high-quality bits.
struct Lcg {
    state: u64,
}

impl Lcg {
    fn new(seed: u64) -> Self {
        Self { state: seed }
    }

    fn next_f32(&mut self, mag: f32) -> f32 {
        self.state = self
            .state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        let bits = (self.state >> 33) as u32;
        let normalized = (bits as f32 / u32::MAX as f32) * 2.0 - 1.0;
        normalized * mag
    }

    fn vec(&mut self, dim: usize, mag: f32) -> Vec<f32> {
        (0..dim).map(|_| self.next_f32(mag)).collect()
    }
}

/// 1-2 ULP equivalent at the magnitudes we test; covers reassociation
/// rounding for vectors of length ≤2048 with magnitudes in [1e-3, 1e3].
const TOL_ABS: f32 = 1e-5;
const TOL_REL: f32 = 1e-4;

fn close(a: f32, b: f32) -> bool {
    let abs_diff = (a - b).abs();
    if abs_diff <= TOL_ABS {
        return true;
    }
    let larger = a.abs().max(b.abs());
    abs_diff <= larger * TOL_REL
}

// --- Property tests over the SIMD-lane-width boundary grid ----------

#[test]
fn dot_matches_reference_on_lcg_inputs() {
    let mut rng = Lcg::new(0xdeadbeef);
    // Cover SIMD-lane boundaries and off-multiples; mid-range
    // magnitudes (the cosine-on-unit-vectors common case) and a
    // wider sweep.
    for &dim in &[
        1_usize, 4, 8, 15, 16, 17, 31, 32, 33, 63, 64, 65, 128, 257, 1024,
    ] {
        for &mag in &[1e-3_f32, 1.0, 10.0, 1e3] {
            let a = rng.vec(dim, mag);
            let b = rng.vec(dim, mag);
            let p = super::dot_product(&a, &b);
            let r = reference_dot(&a, &b);
            assert!(
                close(p, r),
                "dot mismatch at dim={dim} mag={mag}: production={p} reference={r}"
            );
        }
    }
}

#[test]
fn l2_matches_reference_on_lcg_inputs() {
    let mut rng = Lcg::new(0xc0ffee);
    for &dim in &[
        1_usize, 4, 8, 15, 16, 17, 31, 32, 33, 63, 64, 65, 128, 257, 1024,
    ] {
        for &mag in &[1e-3_f32, 1.0, 10.0, 1e3] {
            let a = rng.vec(dim, mag);
            let b = rng.vec(dim, mag);
            let p = super::l2_distance(&a, &b);
            let r = reference_l2(&a, &b);
            assert!(
                close(p, r),
                "l2 mismatch at dim={dim} mag={mag}: production={p} reference={r}"
            );
        }
    }
}

// --- Edge cases per design doc ----

#[test]
fn dot_identical_unit_vector_is_one() {
    // Pre-normalized cosine path: dot(v, v) = 1.0 for unit v.
    let mut v = vec![1.0f32, 2.0, 3.0, 4.0];
    normalize_in_place(&mut v).unwrap();
    let d = super::dot_product(&v, &v);
    assert!(
        close(d, 1.0),
        "self-dot of unit vector should be 1.0, got {d}"
    );
}

#[test]
fn cosine_orthogonal_distance_is_one() {
    // Standard basis vectors are orthogonal — cosine distance = 1.0.
    let e1 = vec![1.0f32, 0.0, 0.0];
    let e2 = vec![0.0f32, 1.0, 0.0];
    let d = distance(&e1, &e2, DistanceMetric::Cosine);
    assert!(
        close(d, 1.0),
        "orthogonal cosine distance should be 1.0, got {d}"
    );
}

#[test]
fn cosine_identical_distance_is_zero() {
    let v = vec![0.6f32, 0.8, 0.0]; // already unit length
    let d = distance(&v, &v, DistanceMetric::Cosine);
    assert!(
        close(d, 0.0),
        "identical cosine distance should be 0.0, got {d}"
    );
}

#[test]
fn l2_identical_distance_is_zero() {
    let v = vec![1.5f32, -2.7, 3.14, 0.0];
    let d = distance(&v, &v, DistanceMetric::L2);
    assert!(
        close(d, 0.0),
        "identical L2 distance should be 0.0, got {d}"
    );
}

#[test]
fn dot_product_off_boundary_remainder_runs() {
    // dim=17 exercises 4 chunked iterations + 1-element tail. Verifies
    // the tail loop is correct.
    let a: Vec<f32> = (0..17).map(|i| (i as f32) * 0.1).collect();
    let b: Vec<f32> = (0..17).map(|i| (16 - i) as f32 * 0.1).collect();
    let p = super::dot_product(&a, &b);
    let r = reference_dot(&a, &b);
    assert!(
        close(p, r),
        "dim=17 tail mismatch: production={p} reference={r}"
    );
}

#[test]
fn l2_distance_off_boundary_remainder_runs() {
    // Same shape but for L2 — also has a tail loop.
    let a: Vec<f32> = (0..33).map(|i| (i as f32) * 0.05).collect();
    let b: Vec<f32> = (0..33).map(|i| ((33 - i) as f32) * 0.05).collect();
    let p = super::l2_distance(&a, &b);
    let r = reference_l2(&a, &b);
    assert!(
        close(p, r),
        "dim=33 tail mismatch: production={p} reference={r}"
    );
}

#[test]
fn dot_handles_small_magnitudes_without_underflow() {
    // Magnitudes near the f32 subnormal boundary. The accumulator
    // must not underflow to zero for typical vectors with small but
    // non-tiny magnitudes.
    let a = vec![1e-20_f32; 64];
    let b = vec![1e-20_f32; 64];
    let d = super::dot_product(&a, &b);
    // 64 * (1e-20)^2 = 6.4e-39; below f32 normal range. Result is
    // expected to underflow to 0 or subnormal — verify it's finite,
    // not NaN.
    assert!(d.is_finite(), "small-mag dot should be finite, got {d}");
}