sc_neurocore_engine/
sobol.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SC-NeuroCore — Sobol quasi-random bitstream generator
7
8//! Sobol quasi-random bitstream generator.
9//!
10//! 1D Sobol (Van der Corput base-2) with Gray-code indexing.
11//! LDS convergence: O(1/N) vs O(1/√N) for Bernoulli.
12
13/// 1D Sobol quasi-random engine using Gray-code traversal.
14#[derive(Clone, Debug)]
15pub struct SobolEngine {
16    state: u32,
17    index: u32,
18    direction: [u32; 32],
19}
20
21impl SobolEngine {
22    /// Initialize with seed (XOR-scramble on direction numbers).
23    pub fn new(seed: u32) -> Self {
24        let mut direction = [0u32; 32];
25        for (i, d) in direction.iter_mut().enumerate() {
26            *d = 1u32 << (31 - i);
27        }
28        Self {
29            state: seed,
30            index: 0,
31            direction,
32        }
33    }
34
35    /// Draw next quasi-random sample in [0, 1).
36    pub fn sample(&mut self) -> f64 {
37        let c = self.index.trailing_ones() as usize;
38        self.state ^= self.direction[c.min(31)];
39        self.index += 1;
40        self.state as f64 / (1u64 << 32) as f64
41    }
42
43    pub fn reset(&mut self, seed: u32) {
44        self.state = seed;
45        self.index = 0;
46    }
47}
48
49/// Generate a Sobol bitstream: threshold comparison of quasi-random samples.
50pub fn generate_sobol_bitstream(p: f64, length: usize, seed: u32) -> Vec<u8> {
51    assert!((0.0..=1.0).contains(&p), "p must be in [0,1]");
52    let mut engine = SobolEngine::new(seed);
53    (0..length)
54        .map(|_| if engine.sample() < p { 1u8 } else { 0u8 })
55        .collect()
56}
57
58/// Pack Sobol bitstream into u64 words.
59pub fn generate_sobol_packed(p: f64, length: usize, seed: u32) -> Vec<u64> {
60    let bits = generate_sobol_bitstream(p, length, seed);
61    crate::bitstream::pack_fast(&bits).data
62}
63
64#[cfg(test)]
65mod tests {
66    use super::*;
67
68    #[test]
69    fn deterministic_output() {
70        let a = generate_sobol_bitstream(0.5, 128, 0);
71        let b = generate_sobol_bitstream(0.5, 128, 0);
72        assert_eq!(a, b);
73    }
74
75    #[test]
76    fn p_zero_all_zeros() {
77        let bits = generate_sobol_bitstream(0.0, 256, 42);
78        assert!(bits.iter().all(|&b| b == 0));
79    }
80
81    #[test]
82    fn p_one_all_ones() {
83        let bits = generate_sobol_bitstream(1.0, 256, 42);
84        assert!(bits.iter().all(|&b| b == 1));
85    }
86
87    #[test]
88    fn proportion_close_to_p() {
89        let p = 0.3;
90        let bits = generate_sobol_bitstream(p, 4096, 0);
91        let ones = bits.iter().filter(|&&b| b == 1).count();
92        let ratio = ones as f64 / 4096.0;
93        assert!(
94            (ratio - p).abs() < 0.02,
95            "Sobol ratio {ratio} should be close to {p}"
96        );
97    }
98
99    #[test]
100    fn lower_discrepancy_than_bernoulli() {
101        // Sobol at p=0.5 should have exactly 50% ones for power-of-2 lengths
102        let bits = generate_sobol_bitstream(0.5, 1024, 0);
103        let ones = bits.iter().filter(|&&b| b == 1).count();
104        let error = (ones as f64 / 1024.0 - 0.5).abs();
105        assert!(error < 0.05, "Sobol discrepancy {error} too high");
106    }
107
108    #[test]
109    fn packed_roundtrip() {
110        let bits = generate_sobol_bitstream(0.7, 200, 99);
111        let packed = generate_sobol_packed(0.7, 200, 99);
112        let unpacked = crate::bitstream::unpack(&crate::bitstream::BitStreamTensor {
113            data: packed,
114            length: 200,
115        });
116        assert_eq!(bits, unpacked);
117    }
118}
sc_neurocore_engine/sobol.rs

sc_neurocore_engine/
sobol.rs
