holographic_memory/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 {
            // Standard entangled generation with groups
            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 {
            // Fallback for small dims: generate active_count indices directly
            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();

        // Birthday collision backfill: if collisions reduced count, fill remaining
        // slots deterministically until we reach exactly active_count.
        let max_backfill_attempts = active_count * 10;
        let mut backfill_counter = 0u64;
        while indices.len() < active_count {
            if backfill_counter as usize >= max_backfill_attempts {
                break;
            }
            let idx = (hash_u64(seed.wrapping_add(0xDEAD), backfill_counter) % dim as u64) as u32;
            backfill_counter += 1;
            if indices.binary_search(&idx).is_err() {
                let pos = indices.partition_point(|&x| x < idx);
                indices.insert(pos, idx);
            }
        }

        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();

        let mut backfill_counter = 0u64;
        while indices.len() < active_count {
            if backfill_counter as usize >= active_count * 10 {
                break;
            }
            let idx = (hash_u64(seed.wrapping_add(0xDEAD), backfill_counter) % dim as u64) as u32;
            backfill_counter += 1;
            if indices.binary_search(&idx).is_err() {
                let pos = indices.partition_point(|&x| x < idx);
                indices.insert(pos, idx);
            }
        }

        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(),
            };
        }

        // Compute random projection dot products.
        // Achlioptas (2003) sparse ternary projection:
        // P(+1) = P(-1) = 1/6, P(0) = 2/3.
        // Preserves JL distances with 3x fewer multiplications.
        // Uses FxHash per (dim, input_pos) pair for deterministic ternary values
        // instead of StdRng per dimension (avoids target_dim RNG inits).
        //
        // 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();

        // Select top active_count by magnitude
        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);

        // Map normalized value to a starting offset in the index space.
        // We use (target_dim - active_count) possible start positions.
        let max_offset = (target_dim - active_count) as f64;
        let start_offset = (normalized * max_offset).round() as usize;

        // Deterministic permutation for global spread with local overlap monotonicity.
        // Must be prime so gcd(PERM_PRIME, dim) = 1 for all dim, ensuring a
        // full-period bijection. (The old golden-ratio constant was divisible by 5,
        // restricting output to 20% of dimensions for 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;
            // Deterministic shuffle: (logical_idx * prime) % target_dim
            let idx = (logical_idx.wrapping_mul(PERM_PRIME) % target_dim as u64) as u32;
            indices.push(idx);
        }

        indices.sort_unstable();
        // PERM_PRIME is prime ⇒ gcd(PERM_PRIME, target_dim)=1 ⇒ mapping is injective.
        // Dedup is a defensive no-op; no backfill needed (unlike random() where
        // hash collisions across entanglement groups are possible).
        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 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 threshold = (n as u32).div_ceil(2);
        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));
        }

        let target_count = (dim / DEFAULT_RHO_DENOM).max(1);

        // If too many pass threshold, keep top target_count by count
        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 HmsCore::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 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 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;
            }
        }
        // IMPORTANT: add/remove adjacency has L1 sensitivity contribution_limit.
        // Every public-domain bin receives noise, including unobserved bins.
        let scale = contribution_limit as f64 / epsilon;
        assert!(
            scale.is_finite() && scale > 0.0,
            "epsilon is outside the supported numeric range"
        );
        let mut rng = rand::rng();
        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)
    }
}

/// 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()
}

/// Splitmix64-based mixer with full avalanche. Two inputs are combined via
/// wrapping arithmetic then passed through the splitmix64 finalizer.
/// Replaces FxHasher which has poor avalanche for close seeds.
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 test_entangled_self_similarity() {
        let v = EntangledHVec::new_deterministic(16384, 42);
        assert!(
            (v.similarity(&v) - 1.0).abs() < 0.0001,
            "Self-similarity should be 1.0"
        );
    }

    #[test]
    fn test_entangled_random_pair_distance() {
        // Two random sparse vectors should have near-zero Jaccard overlap
        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);
        // Create 5 copies of base + 2 random noise
        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);
        // Bundled should be most similar to base (majority)
        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);
        // Should have ~dim/256 active indices
        let expected = 10000 / 256;
        assert_eq!(e.indices.len(), expected);
    }

    #[test]
    fn test_from_scalar() {
        let e = EntangledHVec::from_scalar(0.5, 0.0, 1.0, 10000);
        assert_eq!(e.dim, 10000);
        assert!(!e.indices.is_empty());
    }

    #[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]
    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);
        // Run multiple low-epsilon bundles and check they differ from deterministic
        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_containment_self() {
        let v = EntangledHVec::new_deterministic(16384, 42);
        assert!(
            (v.containment_similarity(&v) - 1.0).abs() < 0.0001,
            "Self-containment should be 1.0"
        );
    }

    #[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);

        let mut prev_sim = 1.0;
        // Check points moving away from 0.5
        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,
        }
    }

    #[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);
    }
}