sc_neurocore_engine/
wilson_cowan.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 Wilson-Cowan 1972 E/I rate model
8
9//! Batch parity with `WilsonCowanUnit.step` in
10//! `src/sc_neurocore/neurons/models/wilson_cowan.py` (Wilson & Cowan
11//! 1972, Biophys. J. 12:1–24).
12//!
13//! Per step:
14//!   s_e = sigmoid(w_ee · E − w_ei · I + ext)
15//!   s_i = sigmoid(w_ie · E − w_ii · I)
16//!   E += (−E + s_e) · dt / τ_e
17//!   I += (−I + s_i) · dt / τ_i
18//!
19//! where `sigmoid(x) = 1 / (1 + exp(−a·(x − θ)))`.
20//!
21//! The model is deterministic (no noise), so bit-exact parity with the
22//! Python primary requires only matching arithmetic — no pre-drawn
23//! RNG buffer needed.
24
25#[inline]
26fn sigmoid(a: f64, theta: f64, x: f64) -> f64 {
27    // Published Wilson-Cowan 1972 two-term form:
28    //   S(x) = 1/(1+exp(-a(x-θ))) − 1/(1+exp(aθ))
29    // Range is [-β, 1-β] where β = 1/(1+exp(aθ)).
30    let baseline = 1.0 / (1.0 + (a * theta).exp());
31    1.0 / (1.0 + (-a * (x - theta)).exp()) - baseline
32}
33
34/// Simulate `ext_input.len()` Wilson-Cowan iterations, writing per-step
35/// `E` and `I` traces into caller-allocated buffers. Returns final
36/// `(E, I)` for convenience.
37#[allow(clippy::too_many_arguments)]
38pub fn simulate(
39    mut e: f64,
40    mut i: f64,
41    w_ee: f64,
42    w_ei: f64,
43    w_ie: f64,
44    w_ii: f64,
45    tau_e: f64,
46    tau_i: f64,
47    a: f64,
48    theta: f64,
49    dt: f64,
50    ext_input: &[f64],
51    e_out: &mut [f64],
52    i_out: &mut [f64],
53) -> (f64, f64) {
54    let n = ext_input.len();
55    assert_eq!(e_out.len(), n, "e_out length mismatch");
56    assert_eq!(i_out.len(), n, "i_out length mismatch");
57
58    for t in 0..n {
59        let s_e = sigmoid(a, theta, w_ee * e - w_ei * i + ext_input[t]);
60        let s_i = sigmoid(a, theta, w_ie * e - w_ii * i);
61        e += (-e + s_e) / tau_e * dt;
62        i += (-i + s_i) / tau_i * dt;
63        e_out[t] = e;
64        i_out[t] = i;
65    }
66    (e, i)
67}
68
69#[cfg(test)]
70mod tests {
71    use super::*;
72
73    fn defaults() -> (f64, f64, f64, f64, f64, f64, f64, f64, f64) {
74        // w_ee, w_ei, w_ie, w_ii, tau_e, tau_i, a, theta, dt
75        (10.0, 6.0, 10.0, 1.0, 1.0, 2.0, 1.2, 4.0, 0.1)
76    }
77
78    #[test]
79    fn sigmoid_monotone_increasing() {
80        let (_, _, _, _, _, _, a, theta, _) = defaults();
81        let lo = sigmoid(a, theta, 0.0);
82        let mid = sigmoid(a, theta, 4.0);
83        let hi = sigmoid(a, theta, 10.0);
84        assert!(lo < mid && mid < hi);
85    }
86
87    #[test]
88    fn sigmoid_at_zero_is_zero() {
89        // Two-term form zeroes the baseline so S(0) = 0 exactly.
90        let (_, _, _, _, _, _, a, theta, _) = defaults();
91        assert!(sigmoid(a, theta, 0.0).abs() < 1e-12);
92    }
93
94    #[test]
95    fn sigmoid_at_theta_equals_half_minus_baseline() {
96        let (_, _, _, _, _, _, a, theta, _) = defaults();
97        let baseline = 1.0 / (1.0 + (a * theta).exp());
98        let r = sigmoid(a, theta, theta);
99        assert!((r - (0.5 - baseline)).abs() < 1e-12);
100    }
101
102    #[test]
103    fn sigmoid_asymptotes_respect_baseline() {
104        // As x → +∞, S(x) → 1 − baseline. As x → −∞, S(x) → −baseline.
105        let (_, _, _, _, _, _, a, theta, _) = defaults();
106        let baseline = 1.0 / (1.0 + (a * theta).exp());
107        assert!((sigmoid(a, theta, 1e6) - (1.0 - baseline)).abs() < 1e-50);
108        assert!((sigmoid(a, theta, -1e6) - (-baseline)).abs() < 1e-50);
109    }
110
111    #[test]
112    fn quiescent_converges() {
113        let (w_ee, w_ei, w_ie, w_ii, tau_e, tau_i, a, theta, dt) = defaults();
114        let n = 20_000;
115        let ext = vec![0.0_f64; n];
116        let mut e_out = vec![0.0_f64; n];
117        let mut i_out = vec![0.0_f64; n];
118        let (e_f, i_f) = simulate(
119            0.1, 0.05, w_ee, w_ei, w_ie, w_ii, tau_e, tau_i, a, theta, dt, &ext, &mut e_out,
120            &mut i_out,
121        );
122        assert!(e_f.is_finite() && i_f.is_finite());
123        assert!(e_f < 0.2, "quiescent E must stay low, got {e_f}");
124        assert!(i_f < 0.2, "quiescent I must stay low, got {i_f}");
125    }
126
127    #[test]
128    fn strong_drive_elevates_activity() {
129        let (w_ee, w_ei, w_ie, w_ii, tau_e, tau_i, a, theta, dt) = defaults();
130        let n = 10_000;
131        let ext = vec![10.0_f64; n];
132        let mut e_out = vec![0.0_f64; n];
133        let mut i_out = vec![0.0_f64; n];
134        let (e_f, _) = simulate(
135            0.1, 0.05, w_ee, w_ei, w_ie, w_ii, tau_e, tau_i, a, theta, dt, &ext, &mut e_out,
136            &mut i_out,
137        );
138        assert!(e_f > 0.3, "strong external drive must elevate E, got {e_f}");
139    }
140
141    #[test]
142    fn output_trace_shape_matches_input() {
143        let n = 64;
144        let ext = vec![1.0_f64; n];
145        let mut e_out = vec![f64::NAN; n];
146        let mut i_out = vec![f64::NAN; n];
147        simulate(
148            0.1, 0.05, 10.0, 6.0, 10.0, 1.0, 1.0, 2.0, 1.2, 4.0, 0.1, &ext, &mut e_out, &mut i_out,
149        );
150        assert!(e_out.iter().all(|v| v.is_finite()));
151        assert!(i_out.iter().all(|v| v.is_finite()));
152    }
153
154    #[test]
155    #[should_panic(expected = "e_out length mismatch")]
156    fn mismatched_e_out_panics() {
157        let n = 10;
158        let ext = vec![0.0_f64; n];
159        let mut e_out = vec![0.0_f64; n + 1];
160        let mut i_out = vec![0.0_f64; n];
161        simulate(
162            0.1, 0.05, 10.0, 6.0, 10.0, 1.0, 1.0, 2.0, 1.2, 4.0, 0.1, &ext, &mut e_out, &mut i_out,
163        );
164    }
165}
sc_neurocore_engine/wilson_cowan.rs

sc_neurocore_engine/
wilson_cowan.rs
