sc_neurocore_engine/neurons/simple_spiking/
mckean.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 — McKean Neuron Model
8
9//! McKean piecewise-linear neuron dynamics.
10
11/// McKean 1970 — piecewise-linear FHN variant.
12#[derive(Clone, Debug)]
13pub struct McKeanNeuron {
14    pub v: f64,
15    pub w: f64,
16    pub a: f64,
17    pub epsilon: f64,
18    pub gamma: f64,
19    pub dt: f64,
20    pub v_peak: f64,
21}
22
23impl McKeanNeuron {
24    pub fn new() -> Self {
25        Self {
26            v: 0.0,
27            w: 0.0,
28            a: 0.25,
29            epsilon: 0.01,
30            gamma: 0.5,
31            dt: 0.1,
32            v_peak: 0.8,
33        }
34    }
35    fn f_v(&self, v: f64) -> f64 {
36        let half_a = self.a / 2.0;
37        let mid = (1.0 + self.a) / 2.0;
38        if v < half_a {
39            -v
40        } else if v < mid {
41            v - self.a
42        } else {
43            1.0 - v
44        }
45    }
46    fn valid_numeric_contract(&self) -> bool {
47        self.v.is_finite()
48            && self.w.is_finite()
49            && self.a.is_finite()
50            && self.epsilon.is_finite()
51            && self.gamma.is_finite()
52            && self.dt.is_finite()
53            && self.v_peak.is_finite()
54            && self.a > 0.0
55            && self.a < 1.0
56            && self.epsilon > 0.0
57            && self.gamma > 0.0
58            && self.dt > 0.0
59    }
60    fn derivatives(&self, v: f64, w: f64, current: f64) -> Option<(f64, f64)> {
61        if !(v.is_finite() && w.is_finite() && current.is_finite()) {
62            return None;
63        }
64        let dv = self.f_v(v) - w + current;
65        let dw = self.epsilon * (v - self.gamma * w);
66        if dv.is_finite() && dw.is_finite() {
67            Some((dv, dw))
68        } else {
69            None
70        }
71    }
72    fn rk4_candidate(&self, current: f64) -> Option<(f64, f64)> {
73        let v0 = self.v;
74        let w0 = self.w;
75        let dt = self.dt;
76        let k1 = self.derivatives(v0, w0, current)?;
77        let k2 = self.derivatives(v0 + 0.5 * dt * k1.0, w0 + 0.5 * dt * k1.1, current)?;
78        let k3 = self.derivatives(v0 + 0.5 * dt * k2.0, w0 + 0.5 * dt * k2.1, current)?;
79        let k4 = self.derivatives(v0 + dt * k3.0, w0 + dt * k3.1, current)?;
80        let next_v = v0 + dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0;
81        let next_w = w0 + dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0;
82        if next_v.is_finite() && next_w.is_finite() {
83            Some((next_v, next_w))
84        } else {
85            None
86        }
87    }
88    pub fn step(&mut self, current: f64) -> i32 {
89        if !self.valid_numeric_contract() || !current.is_finite() {
90            return 0;
91        }
92        let v_prev = self.v;
93        let (next_v, next_w) = match self.rk4_candidate(current) {
94            Some(candidate) => candidate,
95            None => return 0,
96        };
97        self.v = next_v;
98        self.w = next_w;
99        if self.v >= self.v_peak && v_prev < self.v_peak {
100            1
101        } else {
102            0
103        }
104    }
105    /// Run `n_steps` RK4 updates under a constant input, returning the `v` trace
106    /// and the upward-`v_peak`-crossing spike count. Reuses `step`, so the trace
107    /// is bit-identical to the per-step path and to the Python reference (the
108    /// piecewise-linear RHS is exact arithmetic — no transcendental functions).
109    /// The final state is left in `self.v` / `self.w`.
110    pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
111        let mut trace = Vec::with_capacity(n_steps);
112        let mut spikes: i64 = 0;
113        for _ in 0..n_steps {
114            let spiked = self.step(current);
115            trace.push(self.v);
116            spikes += spiked as i64;
117        }
118        (trace, spikes)
119    }
120    pub fn reset(&mut self) {
121        self.v = 0.0;
122        self.w = 0.0;
123    }
124}
125impl Default for McKeanNeuron {
126    fn default() -> Self {
127        Self::new()
128    }
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    #[test]
136    fn default_matches_constructor_state() {
137        let default = McKeanNeuron::default();
138        let constructed = McKeanNeuron::new();
139        assert_eq!(default.v, constructed.v);
140    }
141
142    #[test]
143    fn simulate_matches_repeated_step() {
144        let mut simulated = McKeanNeuron::new();
145        let mut repeated = McKeanNeuron::new();
146        let (trace, spikes) = simulated.simulate(2_000, 0.5);
147        let mut expected_trace = Vec::with_capacity(2_000);
148        let mut expected_spikes = 0_i64;
149        for _ in 0..2_000 {
150            if repeated.step(0.5) == 1 {
151                expected_spikes += 1;
152            }
153            expected_trace.push(repeated.v);
154        }
155        assert_eq!(trace, expected_trace);
156        assert_eq!(spikes, expected_spikes);
157    }
158
159    #[test]
160    fn mckean_fires() {
161        let mut n = McKeanNeuron::new();
162        let t: i32 = (0..2000).map(|_| n.step(0.5)).sum();
163        assert!(t > 0);
164    }
165
166    #[test]
167    fn mckean_reset_clears_state() {
168        let mut n = McKeanNeuron::new();
169        for _ in 0..500 {
170            n.step(0.5);
171        }
172        n.reset();
173        assert!((n.v - 0.0).abs() < 1e-10);
174    }
175
176    #[test]
177    fn mckean_bounded() {
178        let mut n = McKeanNeuron::new();
179        for _ in 0..2000 {
180            n.step(50.0);
181        }
182        assert!(n.v.is_finite());
183    }
184
185    #[test]
186    fn mckean_matches_rk4_candidate() {
187        fn f(v: f64, a: f64) -> f64 {
188            let half_a = a / 2.0;
189            let mid = (1.0 + a) / 2.0;
190            if v < half_a {
191                -v
192            } else if v < mid {
193                v - a
194            } else {
195                1.0 - v
196            }
197        }
198        fn rhs(n: &McKeanNeuron, v: f64, w: f64, current: f64) -> (f64, f64) {
199            (f(v, n.a) - w + current, n.epsilon * (v - n.gamma * w))
200        }
201
202        let mut n = McKeanNeuron {
203            v: 0.2,
204            w: -0.1,
205            ..Default::default()
206        };
207        let current = 0.5;
208        let v0 = n.v;
209        let w0 = n.w;
210        let dt = n.dt;
211        let k1 = rhs(&n, v0, w0, current);
212        let k2 = rhs(&n, v0 + 0.5 * dt * k1.0, w0 + 0.5 * dt * k1.1, current);
213        let k3 = rhs(&n, v0 + 0.5 * dt * k2.0, w0 + 0.5 * dt * k2.1, current);
214        let k4 = rhs(&n, v0 + dt * k3.0, w0 + dt * k3.1, current);
215        let expected_v = v0 + dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0;
216        let expected_w = w0 + dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0;
217
218        assert_eq!(n.step(current), 0);
219        assert!((n.v - expected_v).abs() < 1e-14);
220        assert!((n.w - expected_w).abs() < 1e-14);
221    }
222
223    #[test]
224    fn mckean_nan_no_panic() {
225        let mut n = McKeanNeuron::new();
226        let before = (n.v, n.w);
227        assert_eq!(n.step(f64::NAN), 0);
228        assert_eq!((n.v, n.w), before);
229    }
230
231    #[test]
232    fn mckean_overflow_candidate_preserves_state() {
233        let mut n = McKeanNeuron {
234            v: 1.0e308,
235            w: -1.7e308,
236            ..Default::default()
237        };
238        let before = (n.v, n.w);
239        assert_eq!(n.step(1.7e308), 0);
240        assert_eq!((n.v, n.w), before);
241    }
242
243    #[test]
244    fn mckean_negative_no_crash() {
245        let mut n = McKeanNeuron::new();
246        for _ in 0..500 {
247            n.step(-5.0);
248        }
249        assert!(n.v.is_finite());
250    }
251}
sc_neurocore_engine/neurons/simple_spiking/mckean.rs

sc_neurocore_engine/neurons/simple_spiking/
mckean.rs
