sc_neurocore_engine/neurons/trivial/
quadratic_if.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 — Latham QIF + preserved symmetric SC profile
8
9/// Quadratic Integrate-and-Fire — canonical Type-I excitability.
10/// dv/dt = v² + I, reset at v_peak.
11#[derive(Clone, Debug)]
12pub struct QuadraticIFNeuron {
13    pub v: f64,
14    pub v_reset: f64,
15    pub v_peak: f64,
16    pub dt: f64,
17    pub source_profile: bool,
18}
19
20pub type QuadraticIFCompleteTrace = (Vec<f64>, Vec<u8>, f64);
21
22impl QuadraticIFNeuron {
23    pub fn new(v_reset: f64, v_peak: f64, dt: f64) -> Self {
24        Self {
25            v: v_reset,
26            v_reset,
27            v_peak,
28            dt,
29            source_profile: false,
30        }
31    }
32
33    /// Construct the normalized numerical profile from Latham et al. (2000).
34    pub fn latham_2000() -> Self {
35        Self {
36            v: -1.0,
37            v_reset: -3.0,
38            v_peak: 31.0 / 3.0,
39            dt: 0.05,
40            source_profile: true,
41        }
42    }
43
44    fn valid_numeric_contract(&self) -> bool {
45        self.v.is_finite()
46            && self.v_reset.is_finite()
47            && self.v_peak.is_finite()
48            && self.dt.is_finite()
49            && self.v < self.v_peak
50            && self.v_reset < self.v_peak
51            && self.dt > 0.0
52    }
53
54    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
55        if !self.valid_numeric_contract() || !current.is_finite() {
56            return Err("quadratic-IF state/current violates its finite profile contract");
57        }
58        let (next_v, spiked) = self.exact_candidate(current);
59        if !next_v.is_finite() {
60            return Err("quadratic-IF exact-flow candidate became non-finite");
61        }
62        self.v = next_v;
63        Ok(i32::from(spiked))
64    }
65
66    pub fn step(&mut self, current: f64) -> i32 {
67        self.try_step(current).unwrap_or(0)
68    }
69
70    /// Execute a complete failure-atomic batch against a cloned state.
71    pub fn simulate_complete(
72        &self,
73        n_steps: usize,
74        current: f64,
75    ) -> Result<QuadraticIFCompleteTrace, &'static str> {
76        if !self.valid_numeric_contract() || !current.is_finite() {
77            return Err("invalid quadratic-IF batch contract");
78        }
79        let mut candidate = self.clone();
80        let mut voltage = Vec::with_capacity(n_steps);
81        let mut events = Vec::with_capacity(n_steps);
82        for _ in 0..n_steps {
83            let event = candidate.try_step(current)?;
84            voltage.push(candidate.v);
85            events.push(event as u8);
86        }
87        Ok((voltage, events, candidate.v))
88    }
89
90    pub fn reset(&mut self) {
91        self.v = self.v_reset;
92    }
93
94    fn exact_candidate(&self, current: f64) -> (f64, bool) {
95        if current > 0.0 {
96            let root_i = current.sqrt();
97            let phase = (self.v / root_i).atan();
98            let peak_phase = (self.v_peak / root_i).atan();
99            let next_phase = phase + root_i * self.dt;
100            if next_phase >= peak_phase || next_phase >= std::f64::consts::FRAC_PI_2 {
101                return (self.v_reset, true);
102            }
103            return (root_i * next_phase.tan(), false);
104        }
105        if current == 0.0 {
106            let denominator = 1.0 - self.v * self.dt;
107            if denominator <= 0.0 {
108                return (self.v_reset, true);
109            }
110            let next_v = self.v / denominator;
111            if next_v >= self.v_peak {
112                return (self.v_reset, true);
113            }
114            return (next_v, false);
115        }
116
117        let root_i = (-current).sqrt();
118        if (self.v + root_i).abs() <= 1e-15 {
119            return (self.v, false);
120        }
121        let numerator_ratio = (self.v - root_i) / (self.v + root_i);
122        let evolved_ratio = numerator_ratio * (2.0 * root_i * self.dt).exp();
123        let denominator = 1.0 - evolved_ratio;
124        if (numerator_ratio < 1.0 && evolved_ratio >= 1.0) || denominator.abs() <= 1e-15 {
125            return (self.v_reset, true);
126        }
127        let next_v = root_i * (1.0 + evolved_ratio) / denominator;
128        if next_v >= self.v_peak {
129            (self.v_reset, true)
130        } else {
131            (next_v, false)
132        }
133    }
134}
135
136impl Default for QuadraticIFNeuron {
137    fn default() -> Self {
138        Self::new(-1.0, 1.0, 0.01)
139    }
140}
141
142#[cfg(test)]
143mod tests {
144    use super::*;
145
146    #[test]
147    fn qif_fires_with_positive_input() {
148        let mut n = QuadraticIFNeuron::default();
149        let total: i32 = (0..1000).map(|_| n.step(0.5)).sum();
150        assert!(total > 0);
151    }
152    #[test]
153    fn qif_silent_without_input() {
154        let mut n = QuadraticIFNeuron::default();
155        let t: i32 = (0..1000).map(|_| n.step(0.0)).sum();
156        assert_eq!(t, 0);
157    }
158    #[test]
159    fn qif_reset_clears_state() {
160        let mut n = QuadraticIFNeuron::default();
161        for _ in 0..100 {
162            n.step(0.5);
163        }
164        n.reset();
165        assert!((n.v - n.v_reset).abs() < 1e-10);
166    }
167    #[test]
168    fn qif_bounded() {
169        let mut n = QuadraticIFNeuron::default();
170        for _ in 0..1000 {
171            n.step(10.0);
172        }
173        assert!(n.v.is_finite());
174    }
175    #[test]
176    fn qif_nan_no_panic() {
177        let mut n = QuadraticIFNeuron::default();
178        let before = n.v;
179        assert_eq!(n.step(f64::NAN), 0);
180        assert_eq!(n.v, before);
181    }
182    #[test]
183    fn latham_profile_has_source_boundaries_and_complete_events() {
184        let n = QuadraticIFNeuron::latham_2000();
185        assert_eq!(n.v_reset, -3.0);
186        assert_eq!(n.v_peak, 31.0 / 3.0);
187        assert_eq!(n.dt, 0.05);
188        assert!(n.source_profile);
189        let (trace, events, final_v) = n.simulate_complete(8, 4.0).unwrap();
190        assert_eq!(trace.len(), 8);
191        assert_eq!(events.len(), 8);
192        assert_eq!(trace.last().copied(), Some(final_v));
193    }
194    #[test]
195    fn complete_batch_rejects_without_mutating_source() {
196        let mut n = QuadraticIFNeuron::latham_2000();
197        n.v = -0.25;
198        assert!(n.simulate_complete(2, f64::NAN).is_err());
199        assert_eq!(n.v, -0.25);
200    }
201    #[test]
202    fn qif_nonfinite_increment_preserves_state() {
203        let mut n = QuadraticIFNeuron {
204            v: -0.25,
205            ..Default::default()
206        };
207        let before = n.v;
208        assert_eq!(n.step(-1.0e308), 0);
209        assert_eq!(n.v, before);
210    }
211    #[test]
212    fn qif_matches_exact_positive_current_flow() {
213        let mut n = QuadraticIFNeuron::default();
214        let root_i = 0.5_f64.sqrt();
215        let expected = root_i * ((n.v / root_i).atan() + root_i * n.dt).tan();
216        assert_eq!(n.step(0.5), 0);
217        assert!((n.v - expected).abs() < 1e-12);
218    }
219    #[test]
220    fn qif_preserves_negative_current_fixed_point() {
221        let mut n = QuadraticIFNeuron::default();
222        assert_eq!(n.step(-1.0), 0);
223        assert_eq!(n.v, -1.0);
224    }
225    #[test]
226    fn qif_exact_flow_resets_on_peak_crossing() {
227        let mut n = QuadraticIFNeuron {
228            v: 0.95,
229            dt: 0.5,
230            ..Default::default()
231        };
232        assert_eq!(n.step(1.0), 1);
233        assert_eq!(n.v, n.v_reset);
234    }
235    #[test]
236    fn qif_negative_no_crash() {
237        let mut n = QuadraticIFNeuron::default();
238        for _ in 0..500 {
239            n.step(-5.0);
240        }
241        assert!(n.v.is_finite());
242    }
243}
sc_neurocore_engine/neurons/trivial/quadratic_if.rs

sc_neurocore_engine/neurons/trivial/
quadratic_if.rs
