Skip to main content

sc_neurocore_engine/neurons/trivial/
theta.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 — Theta Neuron
8
9/// Theta neuron — Ermentrout & Kopell canonical form.
10/// dθ/dt = (1 - cosθ) + (1 + cosθ)·I, spike at θ crossing π.
11#[derive(Clone, Debug)]
12pub struct ThetaNeuron {
13    pub theta: f64,
14    pub dt: f64,
15}
16
17pub type ThetaCompleteTrace = (Vec<f64>, Vec<u8>, f64);
18
19impl ThetaNeuron {
20    pub fn new(dt: f64) -> Self {
21        Self { theta: 0.0, dt }
22    }
23
24    fn wrap_phase(theta: f64) -> f64 {
25        (theta + std::f64::consts::PI).rem_euclid(2.0 * std::f64::consts::PI) - std::f64::consts::PI
26    }
27
28    fn valid(&self) -> bool {
29        self.theta.is_finite() && self.dt.is_finite() && self.dt > 0.0
30    }
31
32    fn event_packet_representable(&self, current: f64) -> bool {
33        current <= 0.0 || current.sqrt() * self.dt <= std::f64::consts::PI
34    }
35
36    fn exact_candidate(&self, current: f64) -> (f64, bool) {
37        let y = (self.theta / 2.0).tan();
38        if current > 0.0 {
39            let root_i = current.sqrt();
40            let phase = (y / root_i).atan();
41            let next_phase = phase + root_i * self.dt;
42            if next_phase.cos().abs() <= 1.0e-15 {
43                return (
44                    -std::f64::consts::PI,
45                    next_phase >= std::f64::consts::FRAC_PI_2,
46                );
47            }
48            return (
49                Self::wrap_phase(2.0 * (root_i * next_phase.tan()).atan()),
50                next_phase >= std::f64::consts::FRAC_PI_2,
51            );
52        }
53        if current == 0.0 {
54            let denominator = 1.0 - y * self.dt;
55            if denominator.abs() <= 1.0e-15 {
56                return (-std::f64::consts::PI, true);
57            }
58            return (
59                Self::wrap_phase(2.0 * (y / denominator).atan()),
60                denominator <= 0.0,
61            );
62        }
63
64        let root_i = (-current).sqrt();
65        if (y + root_i).abs() <= 1.0e-15 {
66            return (self.theta, false);
67        }
68        let ratio = (y - root_i) / (y + root_i);
69        let evolved = ratio * (2.0 * root_i * self.dt).exp();
70        let denominator = 1.0 - evolved;
71        let spiked = (ratio < 1.0 && evolved >= 1.0) || denominator.abs() <= 1.0e-15;
72        if spiked && denominator.abs() <= 1.0e-15 {
73            return (-std::f64::consts::PI, true);
74        }
75        (
76            Self::wrap_phase(2.0 * (root_i * (1.0 + evolved) / denominator).atan()),
77            spiked,
78        )
79    }
80
81    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
82        if !current.is_finite() || !self.valid() {
83            return Err("theta state/current must be finite with positive dt");
84        }
85        if !self.event_packet_representable(current) {
86            return Err("theta step can contain more than one source event");
87        }
88        let (next_theta, spiked) = self.exact_candidate(current);
89        if !next_theta.is_finite() {
90            return Err("theta exact-flow candidate became non-finite");
91        }
92        self.theta = Self::wrap_phase(next_theta);
93        Ok(i32::from(spiked))
94    }
95
96    pub fn step(&mut self, current: f64) -> i32 {
97        self.try_step(current).unwrap_or(0)
98    }
99
100    /// Execute a checked complete batch without mutating the source state.
101    pub fn simulate_complete(
102        &self,
103        n_steps: usize,
104        current: f64,
105    ) -> Result<ThetaCompleteTrace, &'static str> {
106        if !self.valid() || !current.is_finite() || !self.event_packet_representable(current) {
107            return Err("invalid theta batch contract");
108        }
109        let mut candidate = self.clone();
110        let mut phase = Vec::with_capacity(n_steps);
111        let mut events = Vec::with_capacity(n_steps);
112        for _ in 0..n_steps {
113            let event = candidate.try_step(current)?;
114            phase.push(candidate.theta);
115            events.push(event as u8);
116        }
117        Ok((phase, events, candidate.theta))
118    }
119
120    pub fn reset(&mut self) {
121        self.theta = 0.0;
122    }
123}
124
125impl Default for ThetaNeuron {
126    fn default() -> Self {
127        Self::new(0.01)
128    }
129}
130
131#[cfg(test)]
132mod tests {
133    use super::*;
134
135    #[test]
136    fn theta_fires() {
137        let mut n = ThetaNeuron::default();
138        let total: i32 = (0..1000).map(|_| n.step(0.5)).sum();
139        assert!(total > 0);
140    }
141    #[test]
142    fn theta_silent_without_input() {
143        let mut n = ThetaNeuron::default();
144        let t: i32 = (0..1000).map(|_| n.step(0.0)).sum();
145        assert_eq!(t, 0);
146    }
147    #[test]
148    fn theta_reset_clears_state() {
149        let mut n = ThetaNeuron::default();
150        for _ in 0..100 {
151            n.step(0.5);
152        }
153        n.reset();
154        assert!((n.theta - 0.0).abs() < 1e-10);
155    }
156    #[test]
157    fn theta_bounded() {
158        let mut n = ThetaNeuron::default();
159        for _ in 0..1000 {
160            n.step(10.0);
161        }
162        assert!(n.theta.is_finite());
163    }
164    #[test]
165    fn theta_nan_no_panic() {
166        ThetaNeuron::default().step(f64::NAN);
167    }
168    #[test]
169    fn theta_exact_positive_flow() {
170        let mut n = ThetaNeuron {
171            theta: 1.0,
172            dt: 0.2,
173        };
174        let root_i = 2.0_f64.sqrt();
175        let phase = ((n.theta / 2.0).tan() / root_i).atan();
176        let expected =
177            ThetaNeuron::wrap_phase(2.0 * (root_i * (phase + root_i * n.dt).tan()).atan());
178        let spike = n.step(2.0);
179        assert_eq!(spike, 0);
180        assert!((n.theta - expected).abs() < 1.0e-12);
181    }
182    #[test]
183    fn theta_exact_flow_reports_within_step_crossing() {
184        let mut n = ThetaNeuron {
185            theta: 2.5,
186            dt: 1.0,
187        };
188        assert_eq!(n.step(1.0), 1);
189        assert!(n.theta >= -std::f64::consts::PI && n.theta <= std::f64::consts::PI);
190    }
191    #[test]
192    fn theta_stable_fixed_point_preserved() {
193        let mut n = ThetaNeuron {
194            theta: -std::f64::consts::FRAC_PI_2,
195            dt: 100.0,
196        };
197        assert_eq!(n.step(-1.0), 0);
198        assert!((n.theta + std::f64::consts::FRAC_PI_2).abs() < 1.0e-12);
199    }
200    #[test]
201    fn theta_non_finite_exact_candidate_preserves_state() {
202        let mut n = ThetaNeuron {
203            theta: 0.25,
204            dt: 1.0e308,
205        };
206        let before = n.theta;
207        assert_eq!(n.step(-1.0e308), 0);
208        assert_eq!(n.theta, before);
209    }
210    #[test]
211    fn theta_complete_batch_is_aligned_and_failure_atomic() {
212        let n = ThetaNeuron {
213            theta: 0.37,
214            dt: 0.037,
215        };
216        let (phase, events, final_theta) = n.simulate_complete(400, 2.2).unwrap();
217        assert_eq!(phase.len(), 400);
218        assert_eq!(events.len(), 400);
219        assert_eq!(phase.last().copied(), Some(final_theta));
220        assert_eq!(n.theta, 0.37);
221        assert!(n.simulate_complete(1, f64::NAN).is_err());
222        assert_eq!(n.theta, 0.37);
223    }
224    #[test]
225    fn theta_rejects_multi_event_step_without_mutation() {
226        let mut n = ThetaNeuron {
227            theta: 0.25,
228            dt: 1.0,
229        };
230        let before = n.theta;
231        assert!(n.try_step(16.0).is_err());
232        assert_eq!(n.theta, before);
233    }
234}