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)) * gain * 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    fn parameters_are_valid(&self) -> bool {
43        self.theta.is_finite()
44            && self.dt.is_finite()
45            && self.dt > 0.0
46            && self.gain.is_finite()
47            && self.theta_threshold.is_finite()
48    }
49
50    /// Checked update. Rejected input or candidates leave phase unchanged.
51    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
52        if !self.parameters_are_valid() {
53            return Err("invalid Ermentrout-Kopell runtime state");
54        }
55        if !current.is_finite() {
56            return Err("invalid Ermentrout-Kopell current");
57        }
58        let input = self.gain * current;
59        if !input.is_finite() {
60            return Err("invalid Ermentrout-Kopell input drive");
61        }
62        let theta_prev = self.theta;
63        let cos_theta = theta_prev.cos();
64        let d_theta = (1.0 - cos_theta) + (1.0 + cos_theta) * input;
65        let theta_next = theta_prev + self.dt * d_theta;
66        if !d_theta.is_finite() || !theta_next.is_finite() {
67            return Err("invalid Ermentrout-Kopell candidate phase");
68        }
69
70        // Spike detection: crossing pi
71        let fired = if theta_next >= self.theta_threshold && theta_prev < self.theta_threshold {
72            1
73        } else {
74            0
75        };
76
77        self.theta = theta_next.rem_euclid(2.0 * std::f64::consts::PI);
78        Ok(fired)
79    }
80
81    /// Infallible NetworkRunner adapter; invalid updates are event-silent and atomic.
82    pub fn step(&mut self, current: f64) -> i32 {
83        self.try_step(current).unwrap_or(0)
84    }
85
86    /// Run `n_steps` under a constant input, returning the `theta` trace
87    /// (wrapped to `[0, 2*pi)`) and the upward-crossing spike count. Reuses
88    /// `step` so the trace matches the per-step path; on a shared libm it also
89    /// matches the Python reference bit-for-bit (the only transcendental is
90    /// `cos`, and the non-chaotic phase flow does not amplify ULP differences).
91    /// The final state is left in `self.theta`.
92    pub fn simulate(
93        &mut self,
94        n_steps: usize,
95        current: f64,
96    ) -> Result<(Vec<f64>, i64), &'static str> {
97        let mut trace = Vec::with_capacity(n_steps);
98        let mut spikes: i64 = 0;
99        for _ in 0..n_steps {
100            let spiked = self.try_step(current)?;
101            trace.push(self.theta);
102            spikes += spiked as i64;
103        }
104        Ok((trace, spikes))
105    }
106
107    pub fn reset(&mut self) {
108        self.theta = 0.0;
109    }
110}
111
112#[cfg(test)]
113mod tests {
114    use super::*;
115
116    #[test]
117    fn ek_fires_with_input() {
118        let mut n = ErmentroutKopellMapNeuron::new();
119        let t: i32 = (0..5000).map(|_| n.step(0.5)).sum();
120        assert!(t > 0, "EK must fire with input, got {t}");
121    }
122
123    #[test]
124    fn ek_silent_without_input() {
125        // Type I: no firing below threshold (I < 0 is subthreshold for theta model)
126        let mut n = ErmentroutKopellMapNeuron::new();
127        let t: i32 = (0..5000).map(|_| n.step(-0.1)).sum();
128        assert_eq!(t, 0, "EK must be silent with negative input, got {t}");
129    }
130
131    #[test]
132    fn ek_type_i_excitability() {
133        // Type I: arbitrarily low firing rate near threshold
134        let mut n_low = ErmentroutKopellMapNeuron::new();
135        let mut n_high = ErmentroutKopellMapNeuron::new();
136        let spikes_low: i32 = (0..10_000).map(|_| n_low.step(0.01)).sum();
137        let spikes_high: i32 = (0..10_000).map(|_| n_high.step(1.0)).sum();
138        assert!(
139            spikes_high > spikes_low,
140            "Higher input → higher rate: high={spikes_high} vs low={spikes_low}"
141        );
142    }
143
144    #[test]
145    fn ek_theta_wraps() {
146        // Theta should stay in [0, 2*pi)
147        let mut n = ErmentroutKopellMapNeuron::new();
148        for _ in 0..10_000 {
149            n.step(0.5);
150        }
151        let two_pi = 2.0 * std::f64::consts::PI;
152        assert!(
153            n.theta >= 0.0 && n.theta < two_pi,
154            "Theta must wrap to [0, 2pi), theta={}",
155            n.theta
156        );
157    }
158
159    #[test]
160    fn ek_negative_input_no_crash() {
161        let mut n = ErmentroutKopellMapNeuron::new();
162        for _ in 0..10_000 {
163            n.step(-100.0);
164        }
165        assert!(n.theta.is_finite());
166    }
167
168    #[test]
169    fn ek_nan_input_is_event_silent_and_atomic() {
170        let mut n = ErmentroutKopellMapNeuron::new();
171        let before = n.theta;
172        assert!(n.try_step(f64::NAN).is_err());
173        assert_eq!(n.step(f64::NAN), 0);
174        assert_eq!(n.theta, before);
175    }
176
177    #[test]
178    fn ek_extreme_input_bounded() {
179        let mut n = ErmentroutKopellMapNeuron::new();
180        for _ in 0..1000 {
181            n.step(1e6);
182        }
183        assert!(n.theta.is_finite());
184    }
185
186    #[test]
187    fn ek_reset_clears_state() {
188        let mut n = ErmentroutKopellMapNeuron::new();
189        n.dt = 0.05;
190        n.gain = 1.5;
191        n.theta_threshold = 2.75;
192        for _ in 0..100 {
193            n.step(0.5);
194        }
195        n.reset();
196        assert_eq!(n.theta, 0.0);
197        assert_eq!((n.dt, n.gain, n.theta_threshold), (0.05, 1.5, 2.75));
198    }
199
200    #[test]
201    fn ek_uses_true_circular_modulo_for_large_finite_steps() {
202        let mut n = ErmentroutKopellMapNeuron::new();
203        n.try_step(1.0e6).unwrap();
204        assert!((0.0..2.0 * std::f64::consts::PI).contains(&n.theta));
205    }
206
207    #[test]
208    fn ek_checked_simulation_returns_complete_trace() {
209        let mut n = ErmentroutKopellMapNeuron::new();
210        let (trace, events) = n.simulate(2_000, 0.5).unwrap();
211        assert_eq!(trace.len(), 2_000);
212        assert_eq!(events, 45);
213        assert_eq!(trace.last().copied(), Some(n.theta));
214    }
215
216    #[test]
217    fn ek_performance_100k_steps() {
218        let start = std::time::Instant::now();
219        let mut n = ErmentroutKopellMapNeuron::new();
220        for _ in 0..100_000 {
221            std::hint::black_box(n.step(0.5));
222        }
223        let elapsed = start.elapsed();
224        assert!(
225            elapsed.as_millis() < 50,
226            "100k steps must complete in <50ms"
227        );
228    }
229}
sc_neurocore_engine/neurons/ermentrout_kopell_map.rs

sc_neurocore_engine/neurons/
ermentrout_kopell_map.rs
