sc_neurocore_engine/neurons/
ermentrout_kopell_map.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 — Ermentrout-Kopell theta map neuron
8
9//! Ermentrout-Kopell theta map neuron.
10
11/// Ermentrout-Kopell canonical Type I — theta neuron in map form.
12///
13/// The canonical model for Type I (saddle-node) excitability.
14/// theta(n+1) = theta(n) + dt * (1 - cos(theta)) + (1 + cos(theta)) * I
15/// Spike when theta crosses pi.
16///
17/// Ermentrout & Kopell, SIAM J Appl Math 46:233, 1986.
18#[derive(Clone, Debug)]
19pub struct ErmentroutKopellMapNeuron {
20    pub theta: f64, // Phase variable [0, 2*pi)
21    pub dt: f64,
22    pub gain: f64,
23    pub theta_threshold: f64,
24}
25
26impl Default for ErmentroutKopellMapNeuron {
27    fn default() -> Self {
28        Self::new()
29    }
30}
31
32impl ErmentroutKopellMapNeuron {
33    pub fn new() -> Self {
34        Self {
35            theta: 0.0,
36            dt: 0.1, // Discrete step size
37            gain: 1.0,
38            theta_threshold: std::f64::consts::PI,
39        }
40    }
41
42    pub fn step(&mut self, current: f64) -> i32 {
43        let input = self.gain * current;
44        let theta_prev = self.theta;
45
46        let d_theta = (1.0 - self.theta.cos()) + (1.0 + self.theta.cos()) * input;
47        self.theta += self.dt * d_theta;
48
49        // Spike detection: crossing pi
50        let fired = if self.theta >= self.theta_threshold && theta_prev < self.theta_threshold {
51            1
52        } else {
53            0
54        };
55
56        // Wrap theta to [0, 2*pi)
57        let two_pi = 2.0 * std::f64::consts::PI;
58        if self.theta >= two_pi {
59            self.theta -= two_pi;
60        }
61        if self.theta < 0.0 {
62            self.theta += two_pi;
63        }
64
65        if !self.theta.is_finite() {
66            self.theta = 0.0;
67        }
68
69        fired
70    }
71
72    /// Run `n_steps` under a constant input, returning the `theta` trace
73    /// (wrapped to `[0, 2*pi)`) and the upward-crossing spike count. Reuses
74    /// `step` so the trace matches the per-step path; on a shared libm it also
75    /// matches the Python reference bit-for-bit (the only transcendental is
76    /// `cos`, and the non-chaotic phase flow does not amplify ULP differences).
77    /// The final state is left in `self.theta`.
78    pub fn simulate(&mut self, n_steps: usize, current: f64) -> (Vec<f64>, i64) {
79        let mut trace = Vec::with_capacity(n_steps);
80        let mut spikes: i64 = 0;
81        for _ in 0..n_steps {
82            let spiked = self.step(current);
83            trace.push(self.theta);
84            spikes += spiked as i64;
85        }
86        (trace, spikes)
87    }
88
89    pub fn reset(&mut self) {
90        *self = Self::new();
91    }
92}
93
94#[cfg(test)]
95mod tests {
96    use super::*;
97
98    #[test]
99    fn ek_fires_with_input() {
100        let mut n = ErmentroutKopellMapNeuron::new();
101        let t: i32 = (0..5000).map(|_| n.step(0.5)).sum();
102        assert!(t > 0, "EK must fire with input, got {t}");
103    }
104
105    #[test]
106    fn ek_silent_without_input() {
107        // Type I: no firing below threshold (I < 0 is subthreshold for theta model)
108        let mut n = ErmentroutKopellMapNeuron::new();
109        let t: i32 = (0..5000).map(|_| n.step(-0.1)).sum();
110        assert_eq!(t, 0, "EK must be silent with negative input, got {t}");
111    }
112
113    #[test]
114    fn ek_type_i_excitability() {
115        // Type I: arbitrarily low firing rate near threshold
116        let mut n_low = ErmentroutKopellMapNeuron::new();
117        let mut n_high = ErmentroutKopellMapNeuron::new();
118        let spikes_low: i32 = (0..10_000).map(|_| n_low.step(0.01)).sum();
119        let spikes_high: i32 = (0..10_000).map(|_| n_high.step(1.0)).sum();
120        assert!(
121            spikes_high > spikes_low,
122            "Higher input → higher rate: high={spikes_high} vs low={spikes_low}"
123        );
124    }
125
126    #[test]
127    fn ek_theta_wraps() {
128        // Theta should stay in [0, 2*pi)
129        let mut n = ErmentroutKopellMapNeuron::new();
130        for _ in 0..10_000 {
131            n.step(0.5);
132        }
133        let two_pi = 2.0 * std::f64::consts::PI;
134        assert!(
135            n.theta >= 0.0 && n.theta < two_pi,
136            "Theta must wrap to [0, 2pi), theta={}",
137            n.theta
138        );
139    }
140
141    #[test]
142    fn ek_negative_input_no_crash() {
143        let mut n = ErmentroutKopellMapNeuron::new();
144        for _ in 0..10_000 {
145            n.step(-100.0);
146        }
147        assert!(n.theta.is_finite());
148    }
149
150    #[test]
151    fn ek_nan_input_stays_finite() {
152        let mut n = ErmentroutKopellMapNeuron::new();
153        n.step(f64::NAN);
154        assert!(n.theta.is_finite());
155    }
156
157    #[test]
158    fn ek_extreme_input_bounded() {
159        let mut n = ErmentroutKopellMapNeuron::new();
160        for _ in 0..1000 {
161            n.step(1e6);
162        }
163        assert!(n.theta.is_finite());
164    }
165
166    #[test]
167    fn ek_reset_clears_state() {
168        let mut n = ErmentroutKopellMapNeuron::new();
169        for _ in 0..100 {
170            n.step(0.5);
171        }
172        n.reset();
173        assert_eq!(n.theta, 0.0);
174    }
175
176    #[test]
177    fn ek_performance_100k_steps() {
178        let start = std::time::Instant::now();
179        let mut n = ErmentroutKopellMapNeuron::new();
180        for _ in 0..100_000 {
181            std::hint::black_box(n.step(0.5));
182        }
183        let elapsed = start.elapsed();
184        assert!(
185            elapsed.as_millis() < 50,
186            "100k steps must complete in <50ms"
187        );
188    }
189}
sc_neurocore_engine/neurons/ermentrout_kopell_map.rs

sc_neurocore_engine/neurons/
ermentrout_kopell_map.rs
