sc_neurocore_engine/neurons/trivial/
energy_lif.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
8//! Source-faithful Fardet-Levina energy-based LIF.
9
10const V_MIN: f64 = -200.0;
11const V_MAX: f64 = 100.0;
12const ENERGY_MAX: f64 = 5.0;
13
14/// Complete state and configuration for the authors' Brian RK4 profile.
15#[derive(Clone, Debug)]
16pub struct EnergyLIFNeuron {
17    pub v: f64,
18    pub epsilon: f64,
19    pub capacitance: f64,
20    pub g_leak: f64,
21    pub e_0: f64,
22    pub e_u: f64,
23    pub e_d: f64,
24    pub e_f: f64,
25    pub v_threshold: f64,
26    pub v_reset: f64,
27    pub alpha: f64,
28    pub epsilon_0: f64,
29    pub epsilon_c: f64,
30    pub delta: f64,
31    pub tau_e: f64,
32    pub dt: f64,
33}
34
35impl EnergyLIFNeuron {
36    /// Construct the pinned author-Brian state and parameter profile.
37    pub fn new() -> Self {
38        Self {
39            v: -61.0,
40            epsilon: 0.32,
41            capacitance: 100.0,
42            g_leak: 9.0,
43            e_0: -62.5,
44            e_u: -58.5,
45            e_d: -40.0,
46            e_f: -62.0,
47            v_threshold: -59.0,
48            v_reset: -62.0,
49            alpha: 1.0,
50            epsilon_0: 0.5,
51            epsilon_c: 0.18,
52            delta: 0.01,
53            tau_e: 200.0,
54            dt: 0.1,
55        }
56    }
57
58    fn rhs(&self, v: f64, epsilon: f64, current: f64) -> (f64, f64) {
59        let leak = self.e_0 + (self.e_u - self.e_0) * (1.0 - epsilon / self.epsilon_0);
60        let dv = (self.g_leak * (leak - v) + current) / self.capacitance;
61        let production = (1.0 - epsilon / (self.alpha * self.epsilon_0)).powi(3);
62        let cost = (v - self.e_f) / (self.e_d - self.e_f);
63        (dv, (production - cost) / self.tau_e)
64    }
65
66    fn candidate(&self, current: f64) -> (f64, f64) {
67        let dt = self.dt;
68        let k1 = self.rhs(self.v, self.epsilon, current);
69        let k2 = self.rhs(
70            self.v + dt * k1.0 / 2.0,
71            self.epsilon + dt * k1.1 / 2.0,
72            current,
73        );
74        let k3 = self.rhs(
75            self.v + dt * k2.0 / 2.0,
76            self.epsilon + dt * k2.1 / 2.0,
77            current,
78        );
79        let k4 = self.rhs(self.v + dt * k3.0, self.epsilon + dt * k3.1, current);
80        (
81            self.v + dt * (k1.0 + 2.0 * k2.0 + 2.0 * k3.0 + k4.0) / 6.0,
82            self.epsilon + dt * (k1.1 + 2.0 * k2.1 + 2.0 * k3.1 + k4.1) / 6.0,
83        )
84    }
85
86    /// Validate the complete source state and integration envelope.
87    pub fn valid(&self) -> bool {
88        [
89            self.v,
90            self.e_0,
91            self.e_u,
92            self.e_d,
93            self.e_f,
94            self.v_threshold,
95            self.v_reset,
96        ]
97        .into_iter()
98        .all(f64::is_finite)
99            && (V_MIN..=V_MAX).contains(&self.v)
100            && (V_MIN..=V_MAX).contains(&self.v_reset)
101            && self.epsilon.is_finite()
102            && (0.0..=ENERGY_MAX).contains(&self.epsilon)
103            && [self.epsilon_0, self.epsilon_c, self.delta]
104                .into_iter()
105                .all(|x| x.is_finite() && x >= 0.0)
106            && [
107                self.capacitance,
108                self.g_leak,
109                self.alpha,
110                self.tau_e,
111                self.dt,
112            ]
113            .into_iter()
114            .all(|x| x.is_finite() && x > 0.0)
115            && self.e_d != self.e_f
116            && self.v_threshold > self.v_reset
117            && self.dt <= 1.0
118            && self.dt <= self.tau_e
119    }
120
121    /// Advance atomically and report invalid transitions as an error.
122    pub fn try_step(&mut self, current: f64) -> Result<i32, &'static str> {
123        if !self.valid() || !current.is_finite() {
124            return Err("invalid EnergyLIF state, configuration, or current");
125        }
126        let (v, epsilon) = self.candidate(current);
127        if !(v.is_finite()
128            && (V_MIN..=V_MAX).contains(&v)
129            && epsilon.is_finite()
130            && (0.0..=ENERGY_MAX).contains(&epsilon))
131        {
132            return Err("EnergyLIF RK4 candidate outside safety envelope");
133        }
134        if v > self.v_threshold && epsilon > self.epsilon_c {
135            let after = epsilon - self.delta;
136            if !(0.0..=ENERGY_MAX).contains(&after) {
137                return Err("EnergyLIF post-spike energy outside safety envelope");
138            }
139            self.v = self.v_reset;
140            self.epsilon = after;
141            return Ok(1);
142        }
143        self.v = v;
144        self.epsilon = epsilon;
145        Ok(0)
146    }
147
148    pub fn step(&mut self, current: f64) -> i32 {
149        self.try_step(current).unwrap_or(-1)
150    }
151
152    /// Restore the source equilibrium-oriented reset state.
153    pub fn reset(&mut self) {
154        self.v = self.e_0;
155        self.epsilon = self.alpha * self.epsilon_0;
156    }
157}
158
159impl Default for EnergyLIFNeuron {
160    fn default() -> Self {
161        Self::new()
162    }
163}
164
165#[cfg(test)]
166mod tests {
167    use super::*;
168    #[test]
169    fn source_transition_is_atomic() {
170        let mut n = EnergyLIFNeuron::new();
171        assert_eq!(n.step(80.0), 0);
172        let before = (n.v, n.epsilon);
173        assert_eq!(n.step(f64::NAN), -1);
174        assert_eq!((n.v, n.epsilon), before);
175    }
176}
sc_neurocore_engine/neurons/trivial/energy_lif.rs

sc_neurocore_engine/neurons/trivial/
energy_lif.rs
