writerslogic_warm/core/

entangled.rs

// Copyright 2024-2026 WritersLogic Contributors
// SPDX-License-Identifier: Apache-2.0

use serde::{Deserialize, Serialize};
use std::borrow::Borrow;

/// Default sparsity denominator: rho = 1/256 of dimensions are active.
///
/// At D=16384 this yields ~64 active indices. This is an engineering choice
/// balancing memory (256x smaller than dense) against discrimination; the
/// cited dense ±1 papers (Laiho 2015, Frady 2021) use different representations.
pub const DEFAULT_RHO_DENOM: usize = 256;

/// Number of entanglement groups (hash-seeded, not learned).
///
/// 64 groups at D=16384 gives ~1 active index per group on average.
/// Groups hash into the full `[0, dim)` range (they do NOT partition the
/// dimension space — collisions across groups are possible).
const NUM_ENTANGLEMENTS: usize = 64;

/// Entangled Sparse Hypervector.
///
/// Instead of a dense bitvec (~50% active), stores only sorted active indices.
/// Sparsity rho = 1/256 means ~D/256 active bits out of D total.
/// Indices are grouped into `NUM_ENTANGLEMENTS` entanglement groups,
/// each generated from a seeded hash.
///
/// Benefits:
/// - 256x less memory per vector at D=16384
/// - Hamming via sorted merge (O(k) where k = active count)
#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct EntangledHVec {
    pub(crate) dim: usize,
    /// Sorted active bit indices (rho * dim total).
    pub(crate) indices: Vec<u32>,
}

impl EntangledHVec {
    /// Create a random entangled sparse hypervector.
    /// Uses entanglement groups when dim is large enough, otherwise direct hash generation.
    /// After dedup, if birthday collisions reduced the count below `active_count`,
    /// deterministic backfill slots are generated until the exact count is reached.
    pub fn new_deterministic(dim: usize, seed: u64) -> Self {
        let active_count = (dim / DEFAULT_RHO_DENOM).max(1);
        let mut indices = Vec::with_capacity(active_count);

        let per_group = active_count / NUM_ENTANGLEMENTS;
        if per_group > 0 {
            for e in 0..NUM_ENTANGLEMENTS {
                let group_seed = hash_u64(seed, e as u64);
                for p in 0..per_group {
                    let idx = hash_u64(group_seed, p as u64) % dim as u64;
                    indices.push(idx as u32);
                }
            }
        } else {
            // Too few active indices to give every group one.
            for p in 0..active_count {
                let idx = hash_u64(seed, p as u64) % dim as u64;
                indices.push(idx as u32);
            }
        }

        indices.sort_unstable();
        indices.dedup();

        backfill_collisions(&mut indices, active_count, dim, seed);

        Self { dim, indices }
    }

    pub fn new_with_density(dim: usize, density_denom: usize, seed: u64) -> Self {
        let active_count = (dim / density_denom.max(1)).max(1);
        let mut indices = Vec::with_capacity(active_count);

        for p in 0..active_count {
            let idx = hash_u64(seed, p as u64) % dim as u64;
            indices.push(idx as u32);
        }

        indices.sort_unstable();
        indices.dedup();

        backfill_collisions(&mut indices, active_count, dim, seed);

        Self { dim, indices }
    }

    /// Create directly from a set of active indices.
    /// Indices must be sorted, deduplicated, and within `[0, dim)`.
    pub fn from_indices(indices: Vec<u32>, dim: usize) -> Self {
        debug_assert!(
            indices.windows(2).all(|w| w[0] < w[1]),
            "from_indices: indices must be sorted and deduplicated"
        );
        debug_assert!(
            indices.last().is_none_or(|&last| (last as usize) < dim),
            "from_indices: indices must be within [0, dim)"
        );
        Self { dim, indices }
    }

    /// Create from a dense float vector via sparse ternary sketching.
    /// Selects top rho*dim indices by Achlioptas ternary projection magnitude.
    pub fn from_dense(dense: &[f32], target_dim: usize) -> Self {
        let active_count = (target_dim / DEFAULT_RHO_DENOM).max(1);
        if dense.is_empty() {
            return Self {
                dim: target_dim,
                indices: Vec::new(),
            };
        }

        // Achlioptas (2003) sparse ternary projection: P(+1) = P(-1) = 1/6, P(0) = 2/3,
        // which preserves JL distances with 3x fewer multiplications. Each ternary value is
        // `hash_u64(output dimension, input position)`, so no RNG is seeded per dimension.
        //
        // The hash pass is branchless: surviving terms are compacted into
        // `terms` (hash outcomes are unpredictable, so branching mispredicts),
        // then summed in ascending-j order, which keeps the f64 accumulation
        // sequence identical to adding each surviving term in index order.
        let vals: Vec<f64> = dense.iter().map(|&v| f64::from(v)).collect();
        let mut terms = vec![0.0f64; vals.len()];
        let mut projections: Vec<(u32, f64)> = (0..target_dim.div_ceil(2))
            .map(|i| {
                let mut kept = 0usize;
                for (j, &val) in vals.iter().enumerate() {
                    let r = hash_u64(i as u64, j as u64) % 6;
                    terms[kept] = if r == 5 { -val } else { val };
                    kept += usize::from((r == 0) | (r == 5));
                }
                let mut dot_product = 0.0f64;
                for &t in &terms[..kept] {
                    dot_product += t;
                }
                let index = i * 2 + usize::from(dot_product < 0.0);
                (
                    index.min(target_dim.saturating_sub(1)) as u32,
                    dot_product.abs(),
                )
            })
            .collect();

        if projections.len() > active_count {
            projections.select_nth_unstable_by(active_count - 1, |a, b| {
                b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)
            });
            projections.truncate(active_count);
        }

        let mut indices: Vec<u32> = projections.into_iter().map(|(idx, _)| idx).collect();
        indices.sort_unstable();

        Self {
            dim: target_dim,
            indices,
        }
    }

    /// Encode a scalar into a sparse hypervector via permuted sliding window.
    ///
    /// The normalized value selects a contiguous window of `active_count`
    /// logical slots in `[0, target_dim - active_count]`. Each logical slot
    /// is then scattered across the index space by multiplying with a prime
    /// constant modulo `target_dim` (golden-ratio permutation). This gives:
    /// - Nearby values share most window slots → high Jaccard
    /// - Distant values share few slots → low Jaccard
    /// - Monotonic: similarity decreases with scalar distance
    pub fn from_scalar(value: f64, min_val: f64, max_val: f64, target_dim: usize) -> Self {
        let range = max_val - min_val;
        let normalized = if range.abs() < f64::EPSILON {
            0.5
        } else {
            ((value - min_val) / range).clamp(0.0, 1.0)
        };
        let active_count = (target_dim / DEFAULT_RHO_DENOM).max(1);

        let max_offset = (target_dim - active_count) as f64;
        let start_offset = (normalized * max_offset).round() as usize;

        // Scatters adjacent window slots across the index space; nearby values still share
        // most slots, so similarity falls monotonically with distance. A constant divisible
        // by 5 would confine output to a fifth of the dimensions when `dim % 5 == 0`.
        const PERM_PRIME: u64 = 0x9E3779B97F4A7C55;

        let mut indices = Vec::with_capacity(active_count);
        for i in 0..active_count {
            let logical_idx = (start_offset + i) as u64;
            let idx = (logical_idx.wrapping_mul(PERM_PRIME) % target_dim as u64) as u32;
            indices.push(idx);
        }

        indices.sort_unstable();
        // The wrapping product is not guaranteed injective modulo `target_dim`.
        indices.dedup();

        Self {
            dim: target_dim,
            indices,
        }
    }

    /// Active bit indices (sorted).
    pub fn indices(&self) -> &[u32] {
        &self.indices
    }

    /// Convert sorted indices to deltas (gaps between consecutive indices).
    /// First element is the absolute value.
    pub fn to_deltas(&self) -> Vec<u32> {
        if self.indices.is_empty() {
            return Vec::new();
        }
        let mut deltas = Vec::with_capacity(self.indices.len());
        deltas.push(self.indices[0]);
        for i in 1..self.indices.len() {
            deltas.push(self.indices[i] - self.indices[i - 1]);
        }
        deltas
    }

    /// Restore absolute indices from deltas.
    pub fn from_deltas(deltas: &[u32], dim: usize) -> Self {
        if deltas.is_empty() {
            return Self {
                dim,
                indices: Vec::new(),
            };
        }
        let mut indices = Vec::with_capacity(deltas.len());
        let mut current: u32 = 0;
        for &d in deltas {
            current = match current.checked_add(d) {
                Some(v) if (v as usize) < dim => v,
                _ => break,
            };
            indices.push(current);
        }
        Self { dim, indices }
    }

    /// Hamming distance via sorted-merge intersection count on sparse index sets.
    pub fn hamming(&self, other: &Self) -> u32 {
        let intersection =
            super::intersection::sparse_intersection_count(&self.indices, &other.indices);
        (self.indices.len() + other.indices.len() - 2 * intersection) as u32
    }

    /// Jaccard similarity: |A∩B| / |A∪B|.
    /// Proper metric for sparse index sets (ρ=1/256).
    /// The old 1−hamming/dim collapsed to ~0.992 for all sparse pairs.
    pub fn similarity(&self, other: &Self) -> f64 {
        let intersection =
            super::intersection::sparse_intersection_count(&self.indices, &other.indices);
        let union = self.indices.len() + other.indices.len() - intersection;
        if union == 0 {
            return 1.0;
        }
        intersection as f64 / union as f64
    }

    /// Entangled bind: symmetric XOR of index sets.
    /// For sparse vectors, XOR of active sets = symmetric difference.
    pub fn bind(&self, other: &Self) -> Self {
        let indices = sorted_symmetric_difference(&self.indices, &other.indices);
        Self {
            dim: self.dim,
            indices,
        }
    }

    /// Permute by shifting all indices by `shifts` modulo dim.
    pub fn permute(&self, shifts: usize) -> Self {
        if self.dim == 0 {
            return self.clone();
        }
        let dim = self.dim as u32;
        let shift = (shifts % self.dim) as u32;
        let mut indices: Vec<u32> = self
            .indices
            .iter()
            .map(|&idx| (idx + shift) % dim)
            .collect();
        indices.sort_unstable();
        Self {
            dim: self.dim,
            indices,
        }
    }

    /// Nonlinear index permutation via keyed hash.
    ///
    /// Maps each active index through hash(idx, seed) % dim, producing a
    /// completely scrambled view. Unlike circular shift (permute), this
    /// destroys spatial correlation between indices, making containment
    /// values independent across different seeds.
    ///
    /// Nearly density-preserving: hash collisions reduce count by ~0.1%
    /// at density 1/256.
    pub fn hash_permute(&self, seed: u64) -> Self {
        if self.dim == 0 {
            return self.clone();
        }
        let d = self.dim as u64;
        let mut indices: Vec<u32> = self
            .indices
            .iter()
            .map(|&idx| (hash_u64(idx as u64, seed) % d) as u32)
            .collect();
        indices.sort_unstable();
        indices.dedup();
        Self {
            dim: self.dim,
            indices,
        }
    }

    /// Bundle multiple sparse vectors with threshold at rho.
    /// Counts per-index frequency, keeps only those above threshold.
    pub fn bundle<V: Borrow<Self>>(vectors: &[V]) -> Self {
        if vectors.is_empty() {
            return Self {
                dim: 0,
                indices: Vec::new(),
            };
        }
        let dim = vectors[0].borrow().dim;
        let n = vectors.len();
        let threshold = (n as u32).div_ceil(2);
        let target_count = (dim / DEFAULT_RHO_DENOM).max(1);

        if let Some(selected) = count_majority(vectors, dim, threshold) {
            return Self::from_selected(selected, dim, target_count);
        }

        let mut all_indices: Vec<u32> = vectors
            .iter()
            .flat_map(|v| v.borrow().indices.iter().copied())
            .collect();
        all_indices.sort_unstable();

        if all_indices.is_empty() {
            return Self {
                dim,
                indices: Vec::new(),
            };
        }

        // Run-length count on sorted indices, filtering by threshold while counting.
        let mut selected: Vec<(u32, u32)> = Vec::new();

        let mut current = all_indices[0];
        let mut count: u32 = 1;
        for &idx in &all_indices[1..] {
            if idx == current {
                count += 1;
            } else {
                if count >= threshold {
                    selected.push((current, count));
                }
                current = idx;
                count = 1;
            }
        }
        if count >= threshold {
            selected.push((current, count));
        }
        Self::from_selected(selected, dim, target_count)
    }

    /// The bundle of the `(index, count)` pairs that reached the majority threshold, in
    /// ascending index order: the `target_count` most frequent, ties resolved by the position
    /// `select_nth_unstable_by` leaves them in, so the order of `selected` must not change.
    pub(crate) fn from_selected(
        mut selected: Vec<(u32, u32)>,
        dim: usize,
        target_count: usize,
    ) -> Self {
        if selected.len() > target_count {
            selected.select_nth_unstable_by(target_count - 1, |a, b| b.1.cmp(&a.1));
            selected.truncate(target_count);
        }

        let mut indices: Vec<u32> = selected.into_iter().map(|(idx, _)| idx).collect();
        indices.sort_unstable();

        Self { dim, indices }
    }

    /// Bloom-filter bundling: set UNION of all input vectors.
    ///
    /// Instead of majority vote (which fails at K≥3 for sparse vectors with
    /// density 1/256), takes the union of all active indices. The bundle grows
    /// denser with more items. Use `containment_similarity` to query.
    ///
    /// Validated by HyperCam (arxiv 2501.10547): Bloom filter bundling achieves
    /// 2 orders of magnitude reduction in encoding operations for sparse HDC.
    pub fn bundle_bloom<V: Borrow<Self>>(vectors: &[V]) -> Self {
        if vectors.is_empty() {
            return Self {
                dim: 0,
                indices: Vec::new(),
            };
        }
        let dim = vectors[0].borrow().dim;

        let mut all_indices: Vec<u32> = vectors
            .iter()
            .flat_map(|v| v.borrow().indices.iter().copied())
            .collect();
        all_indices.sort_unstable();
        all_indices.dedup();

        Self {
            dim,
            indices: all_indices,
        }
    }

    /// Containment similarity: |A∩B| / |A|.
    ///
    /// Measures what fraction of `self`'s active indices appear in `other`.
    /// The correct metric for Bloom-filter bundles: if A was bundled into B,
    /// containment(A, B) ≈ 1.0 regardless of how many other items are in B.
    pub fn containment_similarity(&self, other: &Self) -> f64 {
        if self.indices.is_empty() {
            return 1.0;
        }
        let intersection =
            super::intersection::sparse_intersection_count(&self.indices, &other.indices);
        intersection as f64 / self.indices.len() as f64
    }

    /// Density-corrected containment: (raw - density) / (1 - density).
    ///
    /// Subtracts the expected noise floor (bundle density) and normalizes.
    /// Members: raw=1.0, corrected=1.0. Non-members: raw≈density, corrected≈0.
    /// This is the correct metric for large Bloom bundles where raw containment
    /// has a high floor due to bundle density.
    pub fn corrected_containment(&self, bloom_bundle: &Self) -> f64 {
        let raw = self.containment_similarity(bloom_bundle);
        let density = bloom_bundle.indices.len() as f64 / bloom_bundle.dim.max(1) as f64;
        if density >= 1.0 - 1e-15 {
            return 0.0;
        }
        ((raw - density) / (1.0 - density)).max(0.0)
    }

    /// Non-commutative binding: permute(self, 1) XOR other.
    ///
    /// Captures order: bind_ordered(A, B) ≠ bind_ordered(B, A).
    /// ~2x cost of plain XOR (one permute + one XOR). Inverse: unbind
    /// B first (XOR), then inverse-permute to recover A.
    pub fn bind_ordered(&self, other: &Self) -> Self {
        self.permute(1).bind(other)
    }

    /// Inverse of bind_ordered: given C = bind_ordered(A, B), recover A.
    /// A = permute_inv(C XOR B, 1)
    pub fn unbind_ordered_left(&self, other: &Self) -> Self {
        let xored = self.bind(other);
        xored.permute(self.dim - 1)
    }

    /// Fixed-domain, contribution-clipped histogram followed by noisy top-k.
    /// The dimension must be public; prefer WarmCore::bundle for empty datasets.
    pub fn bundle_dp<V: Borrow<Self>>(vectors: &[V], epsilon: f64) -> Self {
        let dim = vectors.first().map_or(0, |v| v.borrow().dim);
        Self::bundle_dp_in_space(vectors, epsilon, dim)
    }

    pub(crate) fn dp_scale(dim: usize, epsilon: f64) -> Option<f64> {
        let scale = (dim / DEFAULT_RHO_DENOM).max(1) as f64 / epsilon;
        // The inverse-CDF sampler's 53-bit uniform draw has |ln(tail)| < 64.
        (epsilon.is_finite() && epsilon > 0.0 && scale > 0.0 && scale <= f64::MAX / 64.0)
            .then_some(scale)
    }

    pub(crate) fn bundle_dp_in_space<V: Borrow<Self>>(
        vectors: &[V],
        epsilon: f64,
        dim: usize,
    ) -> Self {
        assert!(
            epsilon.is_finite() && epsilon > 0.0,
            "epsilon must be finite and positive"
        );
        assert!(
            vectors.iter().all(|v| v.borrow().dim == dim),
            "bundle dimensions must match"
        );
        if dim == 0 {
            return Self::from_indices(Vec::new(), 0);
        }
        let scale =
            Self::dp_scale(dim, epsilon).expect("epsilon is outside the supported numeric range");
        let contribution_limit = (dim / DEFAULT_RHO_DENOM).max(1);
        let mut counts = vec![0.0; dim];
        for vector in vectors {
            for &index in vector.borrow().indices.iter().take(contribution_limit) {
                counts[index as usize] += 1.0;
            }
        }
        use rand::SeedableRng;
        let mut rng = rand::rngs::StdRng::try_from_rng(&mut rand::rngs::SysRng)
            .expect("operating system randomness unavailable");
        let mut scores: Vec<(u32, f64)> = counts
            .into_iter()
            .enumerate()
            .map(|(i, count)| (i as u32, count + laplace_sample(&mut rng, scale)))
            .collect();
        scores.select_nth_unstable_by(contribution_limit - 1, |a, b| {
            b.1.total_cmp(&a.1).then(a.0.cmp(&b.0))
        });
        scores.truncate(contribution_limit);
        let mut indices: Vec<u32> = scores.into_iter().map(|(index, _)| index).collect();
        indices.sort_unstable();
        Self::from_indices(indices, dim)
    }
}

/// Above this dimension [`EntangledHVec::bundle`] sorts the concatenated indices instead of
/// counting them in a per-thread histogram of `dim` slots (at most 256 KiB).
const BUNDLE_COUNTING_MAX_DIM: usize = 1 << 16;

thread_local! {
    /// [`count_majority`]'s histogram: all-zero between calls.
    static BUNDLE_COUNTS: std::cell::RefCell<Vec<u32>> = const { std::cell::RefCell::new(Vec::new()) };
}

/// The `(index, count)` pairs of `vectors` whose count reaches `threshold`, in ascending index
/// order (what the run-length count over the sorted multiset yields), counted in a histogram.
/// `None` when `dim` is above [`BUNDLE_COUNTING_MAX_DIM`] or an index lies outside `[0, dim)`;
/// the caller then sorts.
fn count_majority<V: Borrow<EntangledHVec>>(
    vectors: &[V],
    dim: usize,
    threshold: u32,
) -> Option<Vec<(u32, u32)>> {
    if dim > BUNDLE_COUNTING_MAX_DIM {
        return None;
    }
    let mut total = 0usize;
    for v in vectors {
        let indices = &v.borrow().indices;
        // Checked in full rather than by the last index: deserialized vectors need not be
        // sorted, and an out-of-range index must fall back rather than panic.
        if indices.iter().any(|&idx| idx as usize >= dim) {
            return None;
        }
        total += indices.len();
    }
    BUNDLE_COUNTS.with(|counts| {
        let mut counts = counts.borrow_mut();
        if counts.len() < dim {
            counts.resize(dim, 0);
        }
        let counts = &mut counts[..dim];
        debug_assert!(counts.iter().all(|&c| c == 0));
        for v in vectors {
            for &idx in &v.borrow().indices {
                counts[idx as usize] += 1;
            }
        }
        let mut selected = Vec::new();
        if total < dim {
            // Visit only the counted slots, zeroing each so its index is taken once, then
            // order the few that reached the threshold.
            for v in vectors {
                for &idx in &v.borrow().indices {
                    let n = std::mem::take(&mut counts[idx as usize]);
                    if n >= threshold {
                        selected.push((idx, n));
                    }
                }
            }
            selected.sort_unstable_by_key(|&(idx, _)| idx);
        } else {
            for (idx, slot) in counts.iter_mut().enumerate() {
                let n = std::mem::take(slot);
                if n >= threshold {
                    selected.push((idx as u32, n));
                }
            }
        }
        Some(selected)
    })
}

/// Sample from Laplace(0, scale) using inverse CDF.
fn laplace_sample(rng: &mut impl rand::RngExt, scale: f64) -> f64 {
    let u = loop {
        let sample = rng.random::<f64>() - 0.5;
        if sample.abs() < 0.5 {
            break sample;
        }
    };
    -scale * u.signum() * (1.0 - 2.0 * u.abs()).ln()
}

/// Fills slots lost to birthday collisions with deterministic indices until `indices` holds
/// `active_count` values, giving up after `10 * active_count` probes. `indices` must be sorted
/// and deduplicated, and stays so.
fn backfill_collisions(indices: &mut Vec<u32>, active_count: usize, dim: usize, seed: u64) {
    let mut counter = 0u64;
    while indices.len() < active_count && (counter as usize) < active_count * 10 {
        let idx = (hash_u64(seed.wrapping_add(0xDEAD), counter) % dim as u64) as u32;
        counter += 1;
        if let Err(pos) = indices.binary_search(&idx) {
            indices.insert(pos, idx);
        }
    }
}

/// Splitmix64-based mixer with full avalanche. Two inputs are combined via
/// wrapping arithmetic then passed through the splitmix64 finalizer.
/// Close seeds give unrelated outputs (Fx hashing's avalanche is too weak for that).
pub fn hash_u64(a: u64, b: u64) -> u64 {
    let mut z = a.wrapping_mul(0x9E3779B97F4A7C15).wrapping_add(b);
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
    z ^ (z >> 31)
}

/// Symmetric difference of two sorted slices (XOR for index sets).
pub(crate) fn sorted_symmetric_difference(a: &[u32], b: &[u32]) -> Vec<u32> {
    let mut result = Vec::with_capacity(a.len() + b.len());
    let mut i = 0;
    let mut j = 0;
    while i < a.len() && j < b.len() {
        match a[i].cmp(&b[j]) {
            std::cmp::Ordering::Less => {
                result.push(a[i]);
                i += 1;
            }
            std::cmp::Ordering::Greater => {
                result.push(b[j]);
                j += 1;
            }
            std::cmp::Ordering::Equal => {
                // In both → cancel (XOR)
                i += 1;
                j += 1;
            }
        }
    }
    result.extend_from_slice(&a[i..]);
    result.extend_from_slice(&b[j..]);
    result
}

/// Count of symmetric difference (without collecting).
pub(crate) fn sorted_symmetric_difference_count(a: &[u32], b: &[u32]) -> usize {
    let mut count = 0;
    let mut i = 0;
    let mut j = 0;
    while i < a.len() && j < b.len() {
        match a[i].cmp(&b[j]) {
            std::cmp::Ordering::Less => {
                count += 1;
                i += 1;
            }
            std::cmp::Ordering::Greater => {
                count += 1;
                j += 1;
            }
            std::cmp::Ordering::Equal => {
                i += 1;
                j += 1;
            }
        }
    }
    count += a.len() - i;
    count += b.len() - j;
    count
}

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

    #[test]
    fn self_similarity_and_containment_are_one() {
        let v = EntangledHVec::new_deterministic(16384, 42);
        assert!(
            (v.similarity(&v) - 1.0).abs() < 0.0001,
            "Self-similarity should be 1.0"
        );
        assert!(
            (v.containment_similarity(&v) - 1.0).abs() < 0.0001,
            "Self-containment should be 1.0"
        );
    }

    #[test]
    fn test_entangled_random_pair_distance() {
        let a = EntangledHVec::new_deterministic(16384, 1);
        let b = EntangledHVec::new_deterministic(16384, 2);
        let sim = a.similarity(&b);
        // Jaccard for random sparse pairs: ~64*64/16384 intersect / ~128 union ≈ 0.002
        assert!(
            sim < 0.05,
            "Random sparse pair Jaccard {:.4} should be near 0",
            sim
        );
    }

    #[test]
    fn test_entangled_bind_involution() {
        let a = EntangledHVec::new_deterministic(16384, 1);
        let b = EntangledHVec::new_deterministic(16384, 2);
        let ab = a.bind(&b);
        let recovered = ab.bind(&b);
        // XOR is involutory: (A⊕B)⊕B = A
        assert_eq!(recovered.indices, a.indices);
    }

    #[test]
    fn test_entangled_bundle_majority() {
        let dim = 16384;
        let base = EntangledHVec::new_deterministic(dim, 100);
        let mut vecs = vec![base.clone(); 5];
        vecs.push(EntangledHVec::new_deterministic(dim, 200));
        vecs.push(EntangledHVec::new_deterministic(dim, 300));
        let bundled = EntangledHVec::bundle(&vecs);
        let sim_base = bundled.similarity(&base);
        let sim_random = bundled.similarity(&EntangledHVec::new_deterministic(dim, 400));
        assert!(
            sim_base > sim_random,
            "Bundle should be closer to majority element"
        );
    }

    #[test]
    fn test_from_dense_produces_sparse() {
        let dense: Vec<f32> = (0..128).map(|i| (i as f32 - 64.0) / 64.0).collect();
        let e = EntangledHVec::from_dense(&dense, 10000);
        assert_eq!(e.dim, 10000);
        let expected = 10000 / 256;
        assert_eq!(e.indices.len(), expected);
    }

    #[test]
    fn test_bundle_dp_produces_valid_output() {
        let dim = 16384;
        let base = EntangledHVec::new_deterministic(dim, 100);
        let vecs = vec![base.clone(); 5];
        let bundled = EntangledHVec::bundle_dp(&vecs, 1.0);
        assert_eq!(bundled.dim, dim);
        assert!(!bundled.indices.is_empty());
        assert!(bundled.indices.windows(2).all(|w| w[0] < w[1]));
    }

    #[test]
    #[should_panic(expected = "epsilon is outside the supported numeric range")]
    fn dp_rejects_finite_scales_whose_inverse_cdf_can_overflow() {
        let vectors = [EntangledHVec::from_indices(vec![1], 256)];
        EntangledHVec::bundle_dp(&vectors, 1e-308);
    }

    #[test]
    fn test_bundle_dp_high_epsilon_matches_deterministic() {
        let dim = 16384;
        let base = EntangledHVec::new_deterministic(dim, 100);
        let vecs = vec![base.clone(); 7];
        let det = EntangledHVec::bundle(&vecs);
        // Very high epsilon = almost no noise
        let dp = EntangledHVec::bundle_dp(&vecs, 1000.0);
        let sim = det.similarity(&dp);
        assert!(
            sim > 0.8,
            "High-epsilon DP bundle should closely match deterministic: {:.4}",
            sim
        );
    }

    #[test]
    fn test_bundle_dp_low_epsilon_adds_noise() {
        let dim = 16384;
        let base = EntangledHVec::new_deterministic(dim, 100);
        let noise = EntangledHVec::new_deterministic(dim, 200);
        let vecs = vec![base.clone(), base.clone(), base.clone(), noise];
        let det = EntangledHVec::bundle(&vecs);
        let mut diffs = 0;
        for _ in 0..5 {
            let dp = EntangledHVec::bundle_dp(&vecs, 0.1);
            if dp.indices != det.indices {
                diffs += 1;
            }
        }
        assert!(diffs > 0, "Low-epsilon DP should produce varied outputs");
    }

    #[test]
    fn test_bloom_bundle_capacity() {
        let dim = 16384;
        let vecs: Vec<EntangledHVec> = (0..200)
            .map(|i| EntangledHVec::new_deterministic(dim, i * 100))
            .collect();
        let bundle = EntangledHVec::bundle_bloom(&vecs);
        assert_eq!(bundle.dim, dim);
        for (i, v) in vecs.iter().enumerate() {
            let sim = v.containment_similarity(&bundle);
            assert!(
                sim > 0.9,
                "Item {} containment in bundle of 200 should be >0.9, got {:.4}",
                i,
                sim
            );
        }
        let noise = EntangledHVec::new_deterministic(dim, 999999);
        let noise_sim = noise.containment_similarity(&bundle);
        assert!(
            noise_sim < 0.9,
            "Non-member containment should be low, got {:.4}",
            noise_sim
        );
    }

    #[test]
    fn test_bind_ordered_non_commutative() {
        let a = EntangledHVec::new_deterministic(16384, 1);
        let b = EntangledHVec::new_deterministic(16384, 2);
        let ab = a.bind_ordered(&b);
        let ba = b.bind_ordered(&a);
        assert_ne!(
            ab.indices, ba.indices,
            "bind_ordered must be non-commutative"
        );
    }

    #[test]
    fn test_bind_ordered_invertible() {
        let a = EntangledHVec::new_deterministic(16384, 1);
        let b = EntangledHVec::new_deterministic(16384, 2);
        let c = a.bind_ordered(&b);
        let recovered = c.unbind_ordered_left(&b);
        assert_eq!(
            recovered.indices, a.indices,
            "unbind_ordered_left should recover A"
        );
    }

    #[test]
    fn test_from_scalar_locality() {
        let dim = 16384;
        let base_val = 0.5;
        let v_base = EntangledHVec::from_scalar(base_val, 0.0, 1.0, dim);
        assert_eq!(v_base.dim, dim);
        assert!(!v_base.indices.is_empty());

        let mut prev_sim = 1.0;
        for i in 1..20 {
            let offset = i as f64 * 0.02;
            let v_curr = EntangledHVec::from_scalar(base_val + offset, 0.0, 1.0, dim);
            let sim = v_base.similarity(&v_curr);

            assert!(
                sim <= prev_sim,
                "Monotonicity violation at offset {}: sim={:.4} was prev={:.4}",
                offset,
                sim,
                prev_sim
            );
            prev_sim = sim;
        }
    }

    /// Verbatim pre-optimisation `from_dense` (encoder "signed-projection-v2").
    fn from_dense_reference(dense: &[f32], target_dim: usize) -> EntangledHVec {
        let active_count = (target_dim / DEFAULT_RHO_DENOM).max(1);
        if dense.is_empty() {
            return EntangledHVec {
                dim: target_dim,
                indices: Vec::new(),
            };
        }
        let mut projections: Vec<(u32, f64)> = (0..target_dim.div_ceil(2))
            .map(|i| {
                let mut dot_product = 0.0f64;
                for (j, &val) in dense.iter().enumerate() {
                    let r = hash_u64(i as u64, j as u64) % 6;
                    if r == 0 {
                        dot_product += val as f64;
                    } else if r == 5 {
                        dot_product -= val as f64;
                    }
                }
                let index = i * 2 + usize::from(dot_product < 0.0);
                (
                    index.min(target_dim.saturating_sub(1)) as u32,
                    dot_product.abs(),
                )
            })
            .collect();
        if projections.len() > active_count {
            projections.select_nth_unstable_by(active_count - 1, |a, b| {
                b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)
            });
            projections.truncate(active_count);
        }
        let mut indices: Vec<u32> = projections.into_iter().map(|(idx, _)| idx).collect();
        indices.sort_unstable();
        EntangledHVec {
            dim: target_dim,
            indices,
        }
    }

    /// `EntangledHVec::bundle` as it was before the histogram: sort the multiset, count runs.
    fn bundle_reference(vectors: &[EntangledHVec]) -> Vec<u32> {
        let dim = vectors[0].dim;
        let threshold = (vectors.len() as u32).div_ceil(2);
        let mut all: Vec<u32> = vectors
            .iter()
            .flat_map(|v| v.indices.iter().copied())
            .collect();
        all.sort_unstable();
        let mut selected: Vec<(u32, u32)> = Vec::new();
        for run in all.chunk_by(|a, b| a == b) {
            if run.len() as u32 >= threshold {
                selected.push((run[0], run.len() as u32));
            }
        }
        let target = (dim / DEFAULT_RHO_DENOM).max(1);
        if selected.len() > target {
            selected.select_nth_unstable_by(target - 1, |a, b| b.1.cmp(&a.1));
            selected.truncate(target);
        }
        let mut indices: Vec<u32> = selected.into_iter().map(|(i, _)| i).collect();
        indices.sort_unstable();
        indices
    }

    #[test]
    fn bundle_matches_sorting_reference() {
        // 65535..=65537 straddle BUNDLE_COUNTING_MAX_DIM; small dims with many vectors make
        // the multiset outgrow dim (dense scan) and force ties at the target count.
        for dim in [1usize, 2, 255, 256, 257, 1024, 16384, 65535, 65536, 65537] {
            for n in [1usize, 2, 3, 4, 7, 10, 64] {
                for overlap in [0u64, 3] {
                    let vectors: Vec<EntangledHVec> = (0..n as u64)
                        .map(|i| {
                            let seed = if overlap > 0 { i % overlap } else { i };
                            EntangledHVec::new_with_density(dim, 16, seed * 7919 + dim as u64)
                        })
                        .collect();
                    assert_eq!(
                        EntangledHVec::bundle(&vectors).indices,
                        bundle_reference(&vectors),
                        "dim {dim} n {n} overlap {overlap}"
                    );
                }
            }
        }
        // An index outside [0, dim) falls back to sorting instead of panicking.
        let bad = vec![
            EntangledHVec {
                dim: 8,
                indices: vec![1, 9],
            },
            EntangledHVec {
                dim: 8,
                indices: vec![1, 9],
            },
        ];
        assert_eq!(EntangledHVec::bundle(&bad).indices, bundle_reference(&bad));
    }

    #[test]
    fn from_dense_matches_reference() {
        let mut state = 0x1234_5678_9ABC_DEF0u64;
        let mut next = move || {
            state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
            hash_u64(state, 0)
        };
        let unit = |x: u64| (x >> 11) as f64 / (1u64 << 53) as f64;
        let mut cases = 0usize;
        for emb in [64usize, 384, 768, 1536] {
            for (d, n) in [(4096usize, 100usize), (16384, 30), (65536, 3)] {
                for k in 0..n {
                    let v: Vec<f32> = (0..emb)
                        .map(|_| {
                            let r = next();
                            match k % 6 {
                                // continuous
                                0 | 1 => (unit(r) * 2.0 - 1.0) as f32,
                                // all zero
                                2 if k < 6 => 0.0,
                                // heavy ties: few distinct values
                                2 | 3 => (r % 3) as f32 - 1.0,
                                // constant
                                4 => 0.5,
                                // extremes and subnormals, finite only
                                _ => match r % 5 {
                                    0 => f32::MAX,
                                    1 => f32::MIN,
                                    2 => f32::MIN_POSITIVE / 4.0,
                                    3 => -0.0,
                                    _ => 0.0,
                                },
                            }
                        })
                        .collect();
                    let got = EntangledHVec::from_dense(&v, d);
                    let want = from_dense_reference(&v, d);
                    assert_eq!(got.dim, want.dim);
                    assert_eq!(got.indices, want.indices, "emb={emb} D={d} case={k}");
                    cases += 1;
                }
            }
        }
        assert!(cases >= 500);
    }
}