sc_neurocore_engine/analysis/
statistics.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later
2// Commercial license available
3// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
4// © Code 2020–2026 Miroslav Šotek. All rights reserved.
5// ORCID: 0009-0009-3560-0851
6// Contact: www.anulum.li | protoscience@anulum.li
7// SC-NeuroCore — Bootstrap and permutation significance testing
8//
9// NOTE: significance_bootstrap takes a function pointer, so it cannot
10// be exposed via PyO3 directly. It is Rust-only, callable from other
11// engine code and tests.
12
13/// Bootstrap significance test for a pairwise statistic.
14///
15/// `statistic_func(a, b) -> f64` is any function computing a scalar
16/// between two spike trains.  Returns `(observed_value, p_value)`.
17pub fn significance_bootstrap<F>(
18    statistic_func: F,
19    train_a: &[f64],
20    train_b: &[f64],
21    n_surrogates: usize,
22    seed: u64,
23) -> (f64, f64)
24where
25    F: Fn(&[f64], &[f64]) -> f64,
26{
27    let observed = statistic_func(train_a, train_b);
28
29    let mut combined: Vec<f64> = Vec::with_capacity(train_a.len() + train_b.len());
30    combined.extend_from_slice(train_a);
31    combined.extend_from_slice(train_b);
32
33    let n_a = train_a.len();
34    let n_total = combined.len();
35
36    // Simple xoshiro256-based PRNG for permutations (Fisher-Yates)
37    let mut rng = SimpleRng::new(seed);
38    let mut count_extreme = 0u64;
39    let mut perm = combined.clone();
40
41    for _ in 0..n_surrogates {
42        perm.copy_from_slice(&combined);
43        // Fisher-Yates shuffle
44        for i in (1..n_total).rev() {
45            let j = rng.next_u64() as usize % (i + 1);
46            perm.swap(i, j);
47        }
48        let surr_val = statistic_func(&perm[..n_a], &perm[n_a..]);
49        if surr_val.abs() >= observed.abs() {
50            count_extreme += 1;
51        }
52    }
53
54    let p_value = (count_extreme + 1) as f64 / (n_surrogates + 1) as f64;
55    (observed, p_value)
56}
57
58/// Minimal PRNG (splitmix64) for deterministic shuffling.
59struct SimpleRng {
60    state: u64,
61}
62
63impl SimpleRng {
64    fn new(seed: u64) -> Self {
65        Self { state: seed }
66    }
67
68    fn next_u64(&mut self) -> u64 {
69        self.state = self.state.wrapping_add(0x9e3779b97f4a7c15);
70        let mut z = self.state;
71        z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
72        z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
73        z ^ (z >> 31)
74    }
75}
76
77#[cfg(test)]
78mod tests {
79    use super::*;
80
81    fn dummy_statistic(a: &[f64], b: &[f64]) -> f64 {
82        let mean_a: f64 = if a.is_empty() {
83            0.0
84        } else {
85            a.iter().sum::<f64>() / a.len() as f64
86        };
87        let mean_b: f64 = if b.is_empty() {
88            0.0
89        } else {
90            b.iter().sum::<f64>() / b.len() as f64
91        };
92        mean_a - mean_b
93    }
94
95    #[test]
96    fn test_bootstrap_identical() {
97        let a = vec![1.0, 1.0, 1.0, 1.0, 1.0];
98        let b = vec![1.0, 1.0, 1.0, 1.0, 1.0];
99        let (obs, p) = significance_bootstrap(dummy_statistic, &a, &b, 200, 42);
100        assert!((obs).abs() < 1e-12);
101        // p-value should be close to 1 (no difference)
102        assert!(p > 0.5);
103    }
104
105    #[test]
106    fn test_bootstrap_different() {
107        let a = vec![10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0];
108        let b = vec![0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
109        let (obs, p) = significance_bootstrap(dummy_statistic, &a, &b, 500, 42);
110        assert!((obs - 10.0).abs() < 1e-12);
111        // Should be significant
112        assert!(
113            p < 0.05,
114            "p-value {p} not significant for clearly different groups"
115        );
116    }
117
118    #[test]
119    fn test_bootstrap_p_range() {
120        let a = vec![1.0, 2.0, 3.0];
121        let b = vec![4.0, 5.0, 6.0];
122        let (_, p) = significance_bootstrap(dummy_statistic, &a, &b, 100, 42);
123        assert!(p > 0.0 && p <= 1.0);
124    }
125
126    #[test]
127    fn test_bootstrap_deterministic() {
128        let a = vec![1.0, 2.0, 3.0, 4.0];
129        let b = vec![5.0, 6.0, 7.0, 8.0];
130        let (obs1, p1) = significance_bootstrap(dummy_statistic, &a, &b, 100, 42);
131        let (obs2, p2) = significance_bootstrap(dummy_statistic, &a, &b, 100, 42);
132        assert!((obs1 - obs2).abs() < 1e-12);
133        assert!((p1 - p2).abs() < 1e-12);
134    }
135
136    #[test]
137    fn test_bootstrap_single_element() {
138        let a = vec![5.0];
139        let b = vec![0.0];
140        let (obs, p) = significance_bootstrap(dummy_statistic, &a, &b, 50, 42);
141        assert!((obs - 5.0).abs() < 1e-12);
142        assert!(p > 0.0);
143    }
144
145    #[test]
146    fn test_simple_rng_deterministic() {
147        let mut rng1 = SimpleRng::new(123);
148        let mut rng2 = SimpleRng::new(123);
149        for _ in 0..100 {
150            assert_eq!(rng1.next_u64(), rng2.next_u64());
151        }
152    }
153}
sc_neurocore_engine/analysis/statistics.rs

sc_neurocore_engine/analysis/
statistics.rs
