sc_neurocore_engine/neuron/
exp_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 — Exponential integrate-and-fire neuron
8
9/// Exponential IF (no adaptation). Fourcaud-Trocmé et al. 2003.
10#[derive(Clone, Debug)]
11pub struct ExpIfNeuron {
12    pub v: f64,
13    pub v_rest: f64,
14    pub v_reset: f64,
15    pub v_threshold: f64,
16    pub v_rh: f64,
17    pub delta_t: f64,
18    pub tau: f64,
19    pub dt: f64,
20    pub refractory_period: f64,
21    pub refractory_remaining: f64,
22    /// False preserves the historical SC RK4 recurrence; true selects source RK2.
23    pub source_profile: bool,
24    pub inv_delta_t: f64,
25    pub dt_div_tau: f64,
26}
27
28/// Explicit rejection classes for checked ExpIF execution.
29#[derive(Clone, Copy, Debug, PartialEq, Eq)]
30pub enum ExpIfError {
31    InvalidInput,
32    InvalidState,
33    NonFiniteUpdate,
34}
35
36/// Aligned voltage, refractory-state, and event traces from one complete batch.
37pub type ExpIfCompleteTrace = (Vec<f64>, Vec<f64>, Vec<u8>);
38
39impl Default for ExpIfNeuron {
40    fn default() -> Self {
41        Self::new()
42    }
43}
44
45impl ExpIfNeuron {
46    pub fn new() -> Self {
47        Self {
48            v: -65.0,
49            v_rest: -65.0,
50            v_reset: -68.0,
51            v_threshold: 30.0,
52            v_rh: -59.9,
53            delta_t: 3.48,
54            tau: 10.0,
55            dt: 0.02,
56            refractory_period: 0.0,
57            refractory_remaining: 0.0,
58            source_profile: false,
59            inv_delta_t: 1.0 / 3.48,
60            dt_div_tau: 0.02 / 10.0,
61        }
62    }
63
64    /// Construct the fitted Fourcaud-Trocmé protocol's deterministic RK2 lane.
65    pub fn fourcaud_trocme_2003() -> Self {
66        Self {
67            v_threshold: -30.0,
68            dt: 0.01,
69            refractory_period: 1.7,
70            source_profile: true,
71            dt_div_tau: 0.01 / 10.0,
72            ..Self::new()
73        }
74    }
75
76    /// Preserve the historical scalar ABI while failing closed on rejection.
77    pub fn step(&mut self, current: f64) -> i32 {
78        self.try_step(current).unwrap_or(0)
79    }
80
81    /// Advance one checked update without mutating the receiver on rejection.
82    pub fn try_step(&mut self, current: f64) -> Result<i32, ExpIfError> {
83        if !self.v.is_finite()
84            || !current.is_finite()
85            || !self.v_rest.is_finite()
86            || !self.v_reset.is_finite()
87            || !self.v_threshold.is_finite()
88            || !self.v_rh.is_finite()
89            || !self.delta_t.is_finite()
90            || !self.tau.is_finite()
91            || !self.dt.is_finite()
92            || !self.refractory_period.is_finite()
93            || !self.refractory_remaining.is_finite()
94            || self.delta_t <= 0.0
95            || self.tau <= 0.0
96            || self.dt <= 0.0
97            || self.refractory_period < 0.0
98            || self.refractory_remaining < 0.0
99            || self.refractory_remaining > self.refractory_period
100            || self.v_threshold <= self.v_rh
101            || self.v >= self.v_threshold
102            || self.v_rest >= self.v_threshold
103            || self.v_reset >= self.v_threshold
104            || (self.source_profile
105                && (self.v_rest != -65.0
106                    || self.v_reset != -68.0
107                    || self.v_threshold != -30.0
108                    || self.v_rh != -59.9
109                    || self.delta_t != 3.48
110                    || self.tau != 10.0
111                    || self.dt >= 0.02
112                    || self.refractory_period != 1.7))
113        {
114            return Err(if !current.is_finite() {
115                ExpIfError::InvalidInput
116            } else {
117                ExpIfError::InvalidState
118            });
119        }
120
121        if self.refractory_remaining > 0.0 {
122            self.refractory_remaining = (self.refractory_remaining - self.dt).max(0.0);
123            self.v = self.v_reset;
124            return Ok(0);
125        }
126
127        let inv_delta_t = 1.0 / self.delta_t;
128        let k1 = self.rhs(self.v, current, inv_delta_t);
129        let predictor = self.v + self.dt * k1;
130        let k2 = self.rhs(
131            if self.source_profile {
132                predictor
133            } else {
134                self.v + 0.5 * self.dt * k1
135            },
136            current,
137            inv_delta_t,
138        );
139        let (k3, k4, next_v) = if self.source_profile {
140            (0.0, 0.0, self.v + 0.5 * self.dt * (k1 + k2))
141        } else {
142            let k3 = self.rhs(self.v + 0.5 * self.dt * k2, current, inv_delta_t);
143            let k4 = self.rhs(self.v + self.dt * k3, current, inv_delta_t);
144            (
145                k3,
146                k4,
147                self.v + (self.dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4),
148            )
149        };
150        if !k1.is_finite()
151            || !k2.is_finite()
152            || !k3.is_finite()
153            || !k4.is_finite()
154            || !next_v.is_finite()
155        {
156            return Err(ExpIfError::NonFiniteUpdate);
157        }
158
159        self.inv_delta_t = inv_delta_t;
160        self.dt_div_tau = self.dt / self.tau;
161        if next_v >= self.v_threshold {
162            self.v = self.v_reset;
163            self.refractory_remaining = self.refractory_period;
164            Ok(1)
165        } else {
166            self.v = next_v;
167            Ok(0)
168        }
169    }
170
171    /// Run a checked, failure-atomic batch and return aligned state/event rows.
172    pub fn simulate_complete(
173        &mut self,
174        n_steps: usize,
175        current: f64,
176    ) -> Result<ExpIfCompleteTrace, ExpIfError> {
177        let mut candidate = self.clone();
178        let mut voltage = Vec::with_capacity(n_steps);
179        let mut refractory = Vec::with_capacity(n_steps);
180        let mut events = Vec::with_capacity(n_steps);
181        for _ in 0..n_steps {
182            let event = candidate.try_step(current)?;
183            voltage.push(candidate.v);
184            refractory.push(candidate.refractory_remaining);
185            events.push(event as u8);
186        }
187        *self = candidate;
188        Ok((voltage, refractory, events))
189    }
190
191    pub fn reset(&mut self) {
192        self.v = self.v_rest;
193        self.refractory_remaining = 0.0;
194    }
195
196    fn rhs(&self, v: f64, current: f64, inv_delta_t: f64) -> f64 {
197        if !v.is_finite() {
198            return f64::NAN;
199        }
200        let bounded_v = v.min(self.v_threshold);
201        let exp_arg = (bounded_v - self.v_rh) * inv_delta_t;
202        let exp_term = self.delta_t * exp_arg.exp();
203        (-(bounded_v - self.v_rest) + exp_term + current) / self.tau
204    }
205}
206
207#[cfg(test)]
208#[path = "exp_if_tests.rs"]
209mod tests;
sc_neurocore_engine/neuron/exp_if.rs

sc_neurocore_engine/neuron/
exp_if.rs
