sc_neurocore_engine/
wong_wang.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 — Rust N-step simulator for the Wong-Wang 2006 decision unit
8
9//! Batch parity with `WongWangUnit.step` in
10//! `src/sc_neurocore/neurons/models/wong_wang.py` (Wong & Wang 2006,
11//! J. Neurosci. 26:1314–1328).
12//!
13//! Per step:
14//!   1. `i_k = j_n * s_k - j_cross * s_(3-k) + i_0 + stim_k + sigma * xi`
15//!   2. `r_k = phi(i_k)` where
16//!      `phi(i) = (a*i - b) / (1 - exp(-d*(a*i - b)))`
17//!      with singularity guard `|a*i - b| < 1e-6 -> 1/d`.
18//!   3. `s_k += (-s_k/tau_s + (1 - s_k) * gamma * r_k) * dt`
19//!   4. clamp `s_k` into `[0, 1]`.
20//!
21//! The Python primary draws `np.random.randn()` twice per step; the
22//! Rust simulator takes `xi` pre-drawn from the Python RNG so trajectories
23//! are bit-exact for matching seeds. This mirrors the ping.rs +
24//! PINGCircuit pattern: Python owns the RNG, Rust owns the inner loop.
25
26const A: f64 = 270.0;
27const B: f64 = 108.0;
28const D: f64 = 0.154;
29
30#[inline]
31fn phi(i_syn: f64) -> f64 {
32    let x = A * i_syn - B;
33    if x.abs() < 1e-6 {
34        1.0 / D
35    } else {
36        x / (1.0 - (-D * x).exp())
37    }
38}
39
40/// Simulate `n_steps` Wong-Wang iterations; write per-step state +
41/// firing-rate traces. `xi` must be length `2 * n_steps` (two noise
42/// samples per step, consumed in `i1, i2` order).
43#[allow(clippy::too_many_arguments)]
44pub fn simulate(
45    mut s1: f64,
46    mut s2: f64,
47    tau_s: f64,
48    gamma: f64,
49    j_n: f64,
50    j_cross: f64,
51    i_0: f64,
52    sigma: f64,
53    dt: f64,
54    stim1: &[f64],
55    stim2: &[f64],
56    xi: &[f64],
57    s1_out: &mut [f64],
58    s2_out: &mut [f64],
59    r1_out: &mut [f64],
60    r2_out: &mut [f64],
61) -> (f64, f64) {
62    let n = stim1.len();
63    assert_eq!(stim2.len(), n, "stim2 length mismatch");
64    assert_eq!(xi.len(), 2 * n, "xi length must be 2 * n_steps");
65    assert_eq!(s1_out.len(), n, "s1_out length mismatch");
66    assert_eq!(s2_out.len(), n, "s2_out length mismatch");
67    assert_eq!(r1_out.len(), n, "r1_out length mismatch");
68    assert_eq!(r2_out.len(), n, "r2_out length mismatch");
69
70    for t in 0..n {
71        let xi1 = xi[2 * t];
72        let xi2 = xi[2 * t + 1];
73        let i1 = j_n * s1 - j_cross * s2 + i_0 + stim1[t] + sigma * xi1;
74        let i2 = j_n * s2 - j_cross * s1 + i_0 + stim2[t] + sigma * xi2;
75        let r1 = phi(i1);
76        let r2 = phi(i2);
77        s1 += (-s1 / tau_s + (1.0 - s1) * gamma * r1) * dt;
78        s2 += (-s2 / tau_s + (1.0 - s2) * gamma * r2) * dt;
79        s1 = s1.clamp(0.0, 1.0);
80        s2 = s2.clamp(0.0, 1.0);
81        s1_out[t] = s1;
82        s2_out[t] = s2;
83        r1_out[t] = r1;
84        r2_out[t] = r2;
85    }
86    (s1, s2)
87}
88
89#[cfg(test)]
90mod tests {
91    use super::*;
92
93    fn params() -> (f64, f64, f64, f64, f64, f64, f64) {
94        // tau_s, gamma, j_n, j_cross, i_0, sigma, dt
95        (0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001)
96    }
97
98    #[test]
99    fn phi_singularity_guard_returns_finite() {
100        // a*i - b == 0 at i = b/a = 0.4
101        let r = phi(B / A);
102        assert!(r.is_finite());
103        assert!((r - 1.0 / D).abs() < 1e-6);
104    }
105
106    #[test]
107    fn phi_monotone_increasing() {
108        let lo = phi(0.5);
109        let hi = phi(1.0);
110        assert!(hi > lo);
111    }
112
113    #[test]
114    fn zero_noise_zero_stim_converges_to_fixed_point() {
115        // With sigma=0, identical initial conditions, balanced dynamics:
116        // s1 and s2 should stay very close and bounded in [0, 1].
117        let (tau_s, gamma, j_n, j_cross, i_0, _, dt) = params();
118        let n = 10_000;
119        let stim = vec![0.0_f64; n];
120        let xi = vec![0.0_f64; 2 * n];
121        let mut s1o = vec![0.0_f64; n];
122        let mut s2o = vec![0.0_f64; n];
123        let mut r1o = vec![0.0_f64; n];
124        let mut r2o = vec![0.0_f64; n];
125        let (s1_f, s2_f) = simulate(
126            0.1, 0.1, tau_s, gamma, j_n, j_cross, i_0, 0.0, dt, &stim, &stim, &xi, &mut s1o,
127            &mut s2o, &mut r1o, &mut r2o,
128        );
129        assert!(s1_f.is_finite() && s2_f.is_finite());
130        assert!((0.0..=1.0).contains(&s1_f));
131        assert!((0.0..=1.0).contains(&s2_f));
132        assert!(
133            (s1_f - s2_f).abs() < 1e-9,
134            "symmetric init must stay symmetric under zero noise"
135        );
136    }
137
138    #[test]
139    fn biased_stimulus_drives_winner() {
140        // Strong stim1, no stim2, no noise → s1 > s2 within a few time constants.
141        let (tau_s, gamma, j_n, j_cross, i_0, _, dt) = params();
142        let n = 50_000;
143        let stim1 = vec![0.2_f64; n];
144        let stim2 = vec![0.0_f64; n];
145        let xi = vec![0.0_f64; 2 * n];
146        let mut s1o = vec![0.0_f64; n];
147        let mut s2o = vec![0.0_f64; n];
148        let mut r1o = vec![0.0_f64; n];
149        let mut r2o = vec![0.0_f64; n];
150        let (s1_f, s2_f) = simulate(
151            0.1, 0.1, tau_s, gamma, j_n, j_cross, i_0, 0.0, dt, &stim1, &stim2, &xi, &mut s1o,
152            &mut s2o, &mut r1o, &mut r2o,
153        );
154        assert!(s1_f > 0.5, "winner s1 should reach attractor; got {s1_f}");
155        assert!(s2_f < 0.2, "loser s2 should be suppressed; got {s2_f}");
156    }
157
158    #[test]
159    fn output_trace_shape_matches_input() {
160        let n = 128;
161        let stim = vec![0.1_f64; n];
162        let xi = vec![0.0_f64; 2 * n];
163        let mut s1o = vec![0.0_f64; n];
164        let mut s2o = vec![0.0_f64; n];
165        let mut r1o = vec![0.0_f64; n];
166        let mut r2o = vec![0.0_f64; n];
167        simulate(
168            0.1, 0.1, 0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001, &stim, &stim, &xi, &mut s1o,
169            &mut s2o, &mut r1o, &mut r2o,
170        );
171        // Every step wrote; nothing still at sentinel zero (for r_out —
172        // s_out might genuinely be near 0, so check r_out which is phi(I) ≥ some positive
173        // floor when I > 0 on this parameter set).
174        assert!(r1o.iter().all(|&r| r > 0.0));
175        assert!(r2o.iter().all(|&r| r > 0.0));
176    }
177
178    #[test]
179    #[should_panic(expected = "xi length must be 2 * n_steps")]
180    fn mismatched_xi_length_panics() {
181        let n = 10;
182        let stim = vec![0.0_f64; n];
183        let xi = vec![0.0_f64; n]; // wrong — should be 2*n
184        let mut s1o = vec![0.0_f64; n];
185        let mut s2o = vec![0.0_f64; n];
186        let mut r1o = vec![0.0_f64; n];
187        let mut r2o = vec![0.0_f64; n];
188        simulate(
189            0.1, 0.1, 0.1, 0.641, 0.2609, 0.0497, 0.3255, 0.0, 0.001, &stim, &stim, &xi, &mut s1o,
190            &mut s2o, &mut r1o, &mut r2o,
191        );
192    }
193}
sc_neurocore_engine/wong_wang.rs

sc_neurocore_engine/
wong_wang.rs
