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
17impl ThetaNeuron {
18    pub fn new(dt: f64) -> Self {
19        Self { theta: 0.0, dt }
20    }
21
22    fn wrap_phase(theta: f64) -> f64 {
23        (theta + std::f64::consts::PI).rem_euclid(2.0 * std::f64::consts::PI) - std::f64::consts::PI
24    }
25
26    fn valid(&self) -> bool {
27        self.theta.is_finite() && self.dt.is_finite() && self.dt > 0.0
28    }
29
30    fn exact_candidate(&self, current: f64) -> (f64, bool) {
31        let y = (self.theta / 2.0).tan();
32        if current > 0.0 {
33            let root_i = current.sqrt();
34            let phase = (y / root_i).atan();
35            let next_phase = phase + root_i * self.dt;
36            if next_phase.cos().abs() <= 1.0e-15 {
37                return (
38                    -std::f64::consts::PI,
39                    next_phase >= std::f64::consts::FRAC_PI_2,
40                );
41            }
42            return (
43                Self::wrap_phase(2.0 * (root_i * next_phase.tan()).atan()),
44                next_phase >= std::f64::consts::FRAC_PI_2,
45            );
46        }
47        if current == 0.0 {
48            let denominator = 1.0 - y * self.dt;
49            if denominator.abs() <= 1.0e-15 {
50                return (-std::f64::consts::PI, true);
51            }
52            return (
53                Self::wrap_phase(2.0 * (y / denominator).atan()),
54                denominator <= 0.0,
55            );
56        }
57
58        let root_i = (-current).sqrt();
59        if (y + root_i).abs() <= 1.0e-15 {
60            return (self.theta, false);
61        }
62        let ratio = (y - root_i) / (y + root_i);
63        let evolved = ratio * (2.0 * root_i * self.dt).exp();
64        let denominator = 1.0 - evolved;
65        let spiked = (ratio < 1.0 && evolved >= 1.0) || denominator.abs() <= 1.0e-15;
66        if spiked && denominator.abs() <= 1.0e-15 {
67            return (-std::f64::consts::PI, true);
68        }
69        (
70            Self::wrap_phase(2.0 * (root_i * (1.0 + evolved) / denominator).atan()),
71            spiked,
72        )
73    }
74
75    pub fn step(&mut self, current: f64) -> i32 {
76        if !current.is_finite() || !self.valid() {
77            return 0;
78        }
79        let (next_theta, spiked) = self.exact_candidate(current);
80        if !next_theta.is_finite() {
81            return 0;
82        }
83        self.theta = Self::wrap_phase(next_theta);
84        if spiked {
85            1
86        } else {
87            0
88        }
89    }
90
91    pub fn reset(&mut self) {
92        self.theta = 0.0;
93    }
94}
95
96impl Default for ThetaNeuron {
97    fn default() -> Self {
98        Self::new(0.01)
99    }
100}
101
102#[cfg(test)]
103mod tests {
104    use super::*;
105
106    #[test]
107    fn theta_fires() {
108        let mut n = ThetaNeuron::default();
109        let total: i32 = (0..1000).map(|_| n.step(0.5)).sum();
110        assert!(total > 0);
111    }
112    #[test]
113    fn theta_silent_without_input() {
114        let mut n = ThetaNeuron::default();
115        let t: i32 = (0..1000).map(|_| n.step(0.0)).sum();
116        assert_eq!(t, 0);
117    }
118    #[test]
119    fn theta_reset_clears_state() {
120        let mut n = ThetaNeuron::default();
121        for _ in 0..100 {
122            n.step(0.5);
123        }
124        n.reset();
125        assert!((n.theta - 0.0).abs() < 1e-10);
126    }
127    #[test]
128    fn theta_bounded() {
129        let mut n = ThetaNeuron::default();
130        for _ in 0..1000 {
131            n.step(10.0);
132        }
133        assert!(n.theta.is_finite());
134    }
135    #[test]
136    fn theta_nan_no_panic() {
137        ThetaNeuron::default().step(f64::NAN);
138    }
139    #[test]
140    fn theta_exact_positive_flow() {
141        let mut n = ThetaNeuron {
142            theta: 1.0,
143            dt: 0.2,
144        };
145        let root_i = 2.0_f64.sqrt();
146        let phase = ((n.theta / 2.0).tan() / root_i).atan();
147        let expected =
148            ThetaNeuron::wrap_phase(2.0 * (root_i * (phase + root_i * n.dt).tan()).atan());
149        let spike = n.step(2.0);
150        assert_eq!(spike, 0);
151        assert!((n.theta - expected).abs() < 1.0e-12);
152    }
153    #[test]
154    fn theta_exact_flow_reports_within_step_crossing() {
155        let mut n = ThetaNeuron {
156            theta: 2.5,
157            dt: 1.0,
158        };
159        assert_eq!(n.step(1.0), 1);
160        assert!(n.theta >= -std::f64::consts::PI && n.theta <= std::f64::consts::PI);
161    }
162    #[test]
163    fn theta_stable_fixed_point_preserved() {
164        let mut n = ThetaNeuron {
165            theta: -std::f64::consts::FRAC_PI_2,
166            dt: 100.0,
167        };
168        assert_eq!(n.step(-1.0), 0);
169        assert!((n.theta + std::f64::consts::FRAC_PI_2).abs() < 1.0e-12);
170    }
171    #[test]
172    fn theta_non_finite_exact_candidate_preserves_state() {
173        let mut n = ThetaNeuron {
174            theta: 0.25,
175            dt: 1.0e308,
176        };
177        let before = n.theta;
178        assert_eq!(n.step(-1.0e308), 0);
179        assert_eq!(n.theta, before);
180    }
181}