Skip to main content

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
8//! Source-bound space-clamped McKean Heaviside system.
9
10const STATE_BOUND: f64 = 1.0e6;
11
12/// Complete state and configuration for the McKean/Tonnelier source equations.
13#[derive(Clone, Debug)]
14pub struct McKeanNeuron {
15    pub v: f64,
16    pub w: f64,
17    pub a: f64,
18    pub lambda: f64,
19    pub mu: f64,
20    pub b: f64,
21    pub dt: f64,
22}
23
24impl McKeanNeuron {
25    /// Construct the normalized source profile with `H(0)=1`.
26    pub fn new() -> Self {
27        Self {
28            v: 0.0,
29            w: 0.0,
30            a: 0.25,
31            lambda: 1.0,
32            mu: 1.0,
33            b: 0.01,
34            dt: 0.1,
35        }
36    }
37    fn rhs(&self, v: f64, w: f64, current: f64) -> (f64, f64) {
38        let h = if v >= self.a { 1.0 } else { 0.0 };
39        (-self.lambda * v + self.mu * h - w + current, self.b * v)
40    }
41    fn candidate(&self, current: f64) -> (f64, f64) {
42        let dt = self.dt;
43        let k1 = self.rhs(self.v, self.w, current);
44        let k2 = self.rhs(self.v + dt * k1.0 / 2.0, self.w + dt * k1.1 / 2.0, current);
45        let k3 = self.rhs(self.v + dt * k2.0 / 2.0, self.w + dt * k2.1 / 2.0, current);
46        let k4 = self.rhs(self.v + dt * k3.0, self.w + dt * k3.1, current);
47        (
48            self.v + dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0,
49            self.w + dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0,
50        )
51    }
52    /// Validate the source state, parameter inequalities, and RK4 envelope.
53    pub fn valid(&self) -> bool {
54        [
55            self.v,
56            self.w,
57            self.a,
58            self.lambda,
59            self.mu,
60            self.b,
61            self.dt,
62        ]
63        .into_iter()
64        .all(f64::is_finite)
65            && self.v.abs() <= STATE_BOUND
66            && self.w.abs() <= STATE_BOUND
67            && self.a > 0.0
68            && self.lambda > 0.0
69            && self.mu > self.lambda * self.a
70            && self.b > 0.0
71            && self.dt > 0.0
72            && self.dt <= 1.0
73    }
74    /// Advance atomically and report invalid transitions as an error.
75    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
76        if !self.valid() || !current.is_finite() {
77            return Err("invalid McKean state, configuration, or current");
78        }
79        let previous = self.v;
80        let (v, w) = self.candidate(current);
81        if !(v.is_finite() && w.is_finite() && v.abs() <= STATE_BOUND && w.abs() <= STATE_BOUND) {
82            return Err("McKean RK4 candidate outside safety envelope");
83        }
84        let event = i32::from(previous < self.a && v >= self.a);
85        self.v = v;
86        self.w = w;
87        Ok(event)
88    }
89    pub fn step(&mut self, current: f64) -> i32 {
90        self.try_step(current).unwrap_or(-1)
91    }
92    /// Restore the source equilibrium state without changing configuration.
93    pub fn reset(&mut self) {
94        self.v = 0.0;
95        self.w = 0.0;
96    }
97}
98impl Default for McKeanNeuron {
99    fn default() -> Self {
100        Self::new()
101    }
102}
103
104#[cfg(test)]
105mod tests {
106    use super::*;
107    #[test]
108    fn source_transition_is_atomic_and_uses_switching_event() {
109        let mut n = McKeanNeuron::new();
110        assert_eq!(n.step(3.0), 1);
111        let before = (n.v, n.w);
112        assert_eq!(n.step(f64::NAN), -1);
113        assert_eq!((n.v, n.w), before);
114    }
115}