Skip to main content

sc_neurocore_engine/neurons/hardware/
dpi_neuron.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 — DPI Neuron Circuit Emulator
8
9/// DPI current-mode adaptive integrate-and-fire circuit.
10///
11/// Implements Indiveri, Stefanini & Chicca (2010), Eqs. (2)–(3), with a
12/// simultaneous explicit-Euler macro-step and a spike-driven refractory pulse.
13#[derive(Clone, Debug)]
14pub struct DPINeuron {
15    pub i_mem: f64,
16    pub i_ahp: f64,
17    pub refractory_time: f64,
18    pub i_threshold: f64,
19    pub i_reset: f64,
20    pub i_rest: f64,
21    pub i_tau: f64,
22    pub i_g: f64,
23    pub i_tau_ahp: f64,
24    pub i_ga: f64,
25    pub i_spike: f64,
26    pub i_0: f64,
27    pub kappa: f64,
28    pub alpha: f64,
29    pub tau: f64,
30    pub tau_ahp: f64,
31    pub refractory_period: f64,
32    pub dt: f64,
33}
34
35impl DPINeuron {
36    pub fn new() -> Self {
37        Self {
38            i_mem: 0.01,
39            i_ahp: 0.01,
40            refractory_time: 0.0,
41            i_threshold: 1.0,
42            i_reset: 0.01,
43            i_rest: 0.1,
44            i_tau: 1.0,
45            i_g: 1.0,
46            i_tau_ahp: 0.1,
47            i_ga: 1.0,
48            i_spike: 5.0,
49            i_0: 0.01,
50            kappa: 0.7,
51            alpha: 10.0,
52            tau: 20.0,
53            tau_ahp: 100.0,
54            refractory_period: 2.0,
55            dt: 0.1,
56        }
57    }
58
59    fn valid(&self) -> bool {
60        self.i_mem.is_finite()
61            && self.i_mem > 0.0
62            && self.i_ahp.is_finite()
63            && self.i_ahp >= 0.0
64            && self.refractory_time.is_finite()
65            && self.refractory_time >= 0.0
66            && self.i_threshold.is_finite()
67            && self.i_threshold > 0.0
68            && self.i_reset.is_finite()
69            && self.i_reset > 0.0
70            && self.i_reset < self.i_threshold
71            && self.i_rest.is_finite()
72            && self.i_rest >= 0.0
73            && self.i_tau.is_finite()
74            && self.i_tau > 0.0
75            && self.i_g.is_finite()
76            && self.i_g > 0.0
77            && self.i_tau_ahp.is_finite()
78            && self.i_tau_ahp > 0.0
79            && self.i_ga.is_finite()
80            && self.i_ga > 0.0
81            && self.i_spike.is_finite()
82            && self.i_spike > 0.0
83            && self.i_0.is_finite()
84            && self.i_0 > 0.0
85            && self.kappa.is_finite()
86            && self.kappa > 0.0
87            && self.alpha.is_finite()
88            && self.alpha > 0.0
89            && self.tau.is_finite()
90            && self.tau > 0.0
91            && self.tau_ahp.is_finite()
92            && self.tau_ahp > 0.0
93            && self.refractory_period.is_finite()
94            && self.refractory_period > 0.0
95            && self.dt.is_finite()
96            && self.dt > 0.0
97            && self.refractory_period >= self.dt
98    }
99
100    fn sigmoid(value: f64) -> f64 {
101        if value >= 0.0 {
102            1.0 / (1.0 + (-value).exp())
103        } else {
104            let exponential = value.exp();
105            exponential / (1.0 + exponential)
106        }
107    }
108
109    fn feedback_current(&self) -> f64 {
110        let log_current = (self.i_0.ln() + self.kappa * self.i_mem.ln()) / (self.kappa + 1.0);
111        log_current.exp() * Self::sigmoid(self.alpha * (self.i_mem - self.i_threshold))
112    }
113
114    pub fn step(&mut self, current: f64) -> i32 {
115        if !current.is_finite() || !self.valid() {
116            return 0;
117        }
118        let total_input = self.i_rest + current;
119        if !total_input.is_finite() || total_input < 0.0 {
120            return 0;
121        }
122
123        let spike_active = self.refractory_time > 0.0;
124        let spike_current = if spike_active { self.i_spike } else { 0.0 };
125        let d_i_ahp = self.i_ahp / (self.tau_ahp * self.i_tau_ahp)
126            * (spike_current / (1.0 + self.i_ahp / self.i_ga) - self.i_tau_ahp);
127        let next_i_ahp = self.i_ahp + self.dt * d_i_ahp;
128
129        let (next_i_mem, next_refractory, spiked) = if spike_active {
130            (
131                self.i_reset,
132                (self.refractory_time - self.dt).max(0.0),
133                false,
134            )
135        } else {
136            let i_fb = self.feedback_current();
137            let d_i_mem = self.i_mem / (self.tau * self.i_tau)
138                * (total_input / (1.0 + self.i_mem / self.i_g) - self.i_tau + i_fb - self.i_ahp);
139            let candidate = self.i_mem + self.dt * d_i_mem;
140            if !candidate.is_finite() || candidate <= 0.0 {
141                return 0;
142            }
143            if candidate >= self.i_threshold {
144                (self.i_reset, self.refractory_period, true)
145            } else {
146                (candidate, 0.0, false)
147            }
148        };
149
150        if !next_i_mem.is_finite()
151            || !next_i_ahp.is_finite()
152            || !next_refractory.is_finite()
153            || next_i_mem <= 0.0
154            || next_i_ahp < 0.0
155            || next_refractory < 0.0
156        {
157            return 0;
158        }
159
160        self.i_mem = next_i_mem;
161        self.i_ahp = next_i_ahp;
162        self.refractory_time = next_refractory;
163        i32::from(spiked)
164    }
165
166    pub fn reset(&mut self) {
167        self.i_mem = self.i_reset;
168        self.i_ahp = self.i_0;
169        self.refractory_time = 0.0;
170    }
171}
172impl Default for DPINeuron {
173    fn default() -> Self {
174        Self::new()
175    }
176}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181
182    #[test]
183    fn dpi_fires() {
184        let mut n = DPINeuron::new();
185        let t: i32 = (0..20_000).map(|_| n.step(5.0)).sum();
186        assert_eq!(t, 44);
187    }
188    #[test]
189    fn dpi_silent() {
190        let mut n = DPINeuron::new();
191        let t: i32 = (0..100).map(|_| n.step(0.0)).sum();
192        assert_eq!(t, 0);
193    }
194    #[test]
195    fn dpi_reset() {
196        let mut n = DPINeuron::new();
197        for _ in 0..500 {
198            n.step(5.0);
199        }
200        n.reset();
201        assert_eq!(n.i_mem, n.i_reset);
202        assert_eq!(n.i_ahp, n.i_0);
203        assert_eq!(n.refractory_time, 0.0);
204    }
205    #[test]
206    fn dpi_nan_no_panic() {
207        let mut n = DPINeuron::new();
208        let old = (n.i_mem, n.i_ahp, n.refractory_time);
209        assert_eq!(n.step(f64::NAN), 0);
210        assert_eq!((n.i_mem, n.i_ahp, n.refractory_time), old);
211    }
212    #[test]
213    fn dpi_coupled_euler_anchor() {
214        let mut n = DPINeuron::new();
215        assert_eq!(n.step(5.0), 0);
216        assert!((n.i_mem - 0.010201975272610835).abs() < 1.0e-17);
217        assert!((n.i_ahp - 0.00999).abs() < 1.0e-17);
218    }
219    #[test]
220    fn dpi_refractory_pulse_drives_adaptation() {
221        let mut n = DPINeuron {
222            refractory_time: 2.0,
223            ..DPINeuron::new()
224        };
225        assert_eq!(n.step(0.0), 0);
226        assert_eq!(n.i_mem, n.i_reset);
227        assert!(n.i_ahp > 0.01);
228        assert_eq!(n.refractory_time, 1.9);
229    }
230}