sc_neurocore_engine/neurons/simple_spiking/
alpha.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 — Dual alpha-synapse leaky integrate-and-fire engine
8
9//! Exact piecewise-constant-input flow for the dual alpha-synapse LIF: a leaky
10//! membrane relaxation driven by two five-state alpha-filter cascades, with
11//! the exact alpha-current convolution including the equal-time-constant limit.
12
13use std::error::Error;
14use std::fmt;
15
16/// Typed caller-contract and numerical failures.
17#[derive(Clone, Copy, Debug, Eq, PartialEq)]
18pub enum AlphaError {
19    /// A state or parameter is not finite.
20    NonFiniteConfiguration,
21    /// A magnitude or ordering constraint is violated.
22    InvalidScale,
23    /// One external-current value is not finite.
24    NonFiniteInput,
25    /// The exact-flow candidate is not finite.
26    NonFiniteCandidate,
27}
28
29impl fmt::Display for AlphaError {
30    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
31        let message = match self {
32            Self::NonFiniteConfiguration => "alpha state and parameters must be finite",
33            Self::InvalidScale => {
34                "alpha tau_v/tau_exc/tau_inh/dt must be positive and v_threshold must exceed v_rest"
35            }
36            Self::NonFiniteInput => "alpha input must contain only finite values",
37            Self::NonFiniteCandidate => "alpha exact-flow candidate must remain finite",
38        };
39        formatter.write_str(message)
40    }
41}
42
43impl Error for AlphaError {}
44
45/// Dual excitatory/inhibitory alpha-synapse leaky integrate-and-fire neuron.
46#[derive(Clone, Debug)]
47pub struct AlphaNeuron {
48    /// Membrane potential.
49    pub v: f64,
50    /// Excitatory alpha-rise state.
51    pub a_exc: f64,
52    /// Excitatory synaptic current.
53    pub i_exc: f64,
54    /// Inhibitory alpha-rise state.
55    pub a_inh: f64,
56    /// Inhibitory synaptic current.
57    pub i_inh: f64,
58    /// Leak reversal potential, also the somatic spike reset.
59    pub v_rest: f64,
60    /// Spike threshold; must exceed `v_rest`.
61    pub v_threshold: f64,
62    /// Positive membrane time constant.
63    pub tau_v: f64,
64    /// Positive excitatory alpha time constant.
65    pub tau_exc: f64,
66    /// Positive inhibitory alpha time constant.
67    pub tau_inh: f64,
68    /// Positive piecewise-constant-input sampling interval.
69    pub dt: f64,
70}
71
72impl AlphaNeuron {
73    /// Construct the catalogue model-family defaults.
74    pub fn new() -> Self {
75        Self {
76            v: 0.0,
77            a_exc: 0.0,
78            i_exc: 0.0,
79            a_inh: 0.0,
80            i_inh: 0.0,
81            v_rest: 0.0,
82            v_threshold: 1.0,
83            tau_v: 20.0,
84            tau_exc: 5.0,
85            tau_inh: 10.0,
86            dt: 1.0,
87        }
88    }
89
90    /// Construct and validate a complete numerical configuration.
91    #[allow(clippy::too_many_arguments)]
92    pub fn with_parameters(
93        v: f64,
94        a_exc: f64,
95        i_exc: f64,
96        a_inh: f64,
97        i_inh: f64,
98        v_rest: f64,
99        v_threshold: f64,
100        tau_v: f64,
101        tau_exc: f64,
102        tau_inh: f64,
103        dt: f64,
104    ) -> Result<Self, AlphaError> {
105        let neuron = Self {
106            v,
107            a_exc,
108            i_exc,
109            a_inh,
110            i_inh,
111            v_rest,
112            v_threshold,
113            tau_v,
114            tau_exc,
115            tau_inh,
116            dt,
117        };
118        neuron.validate()?;
119        Ok(neuron)
120    }
121
122    fn validate(&self) -> Result<(), AlphaError> {
123        if ![
124            self.v,
125            self.a_exc,
126            self.i_exc,
127            self.a_inh,
128            self.i_inh,
129            self.v_rest,
130            self.v_threshold,
131            self.tau_v,
132            self.tau_exc,
133            self.tau_inh,
134            self.dt,
135        ]
136        .into_iter()
137        .all(f64::is_finite)
138        {
139            return Err(AlphaError::NonFiniteConfiguration);
140        }
141        if self.tau_v <= 0.0
142            || self.tau_exc <= 0.0
143            || self.tau_inh <= 0.0
144            || self.dt <= 0.0
145            || self.v_threshold <= self.v_rest
146        {
147            return Err(AlphaError::InvalidScale);
148        }
149        Ok(())
150    }
151
152    fn filter_candidates(
153        rise_state: f64,
154        current_state: f64,
155        drive: f64,
156        tau: f64,
157        dt: f64,
158        decay: f64,
159    ) -> Result<(f64, f64), AlphaError> {
160        let steady_state = tau * drive;
161        let rise_delta = rise_state - steady_state;
162        let current_delta = current_state - steady_state;
163        let rise_next = steady_state + rise_delta * decay;
164        let current_next = steady_state + decay * (current_delta + rise_delta * dt / tau);
165        if !rise_next.is_finite() || !current_next.is_finite() {
166            return Err(AlphaError::NonFiniteCandidate);
167        }
168        Ok((rise_next, current_next))
169    }
170
171    fn drive_contribution(
172        current_delta: f64,
173        rise_delta: f64,
174        tau_drive: f64,
175        dt: f64,
176        rates: [f64; 2],
177        decays: [f64; 2],
178    ) -> Result<f64, AlphaError> {
179        let [rate_v, rate_drive] = rates;
180        let [decay_v, decay_drive] = decays;
181        let contribution = if (rate_v - rate_drive).abs() <= 1.0e-14 {
182            rate_v * decay_v * (current_delta * dt + rise_delta * dt * dt / (2.0 * tau_drive))
183        } else {
184            let rate_delta = rate_v - rate_drive;
185            let first_order = current_delta * (decay_drive - decay_v) / rate_delta;
186            let second_order = rise_delta / tau_drive
187                * (decay_drive * (rate_delta * dt - 1.0) + decay_v)
188                / (rate_delta * rate_delta);
189            rate_v * (first_order + second_order)
190        };
191        if !contribution.is_finite() {
192            return Err(AlphaError::NonFiniteCandidate);
193        }
194        Ok(contribution)
195    }
196
197    /// Advance one exact-flow interval with explicit error reporting.
198    #[inline]
199    pub fn try_step(&mut self, exc_current: f64, inh_current: f64) -> Result<i32, AlphaError> {
200        self.validate()?;
201        if !exc_current.is_finite() || !inh_current.is_finite() {
202            return Err(AlphaError::NonFiniteInput);
203        }
204        let rate_v = 1.0 / self.tau_v;
205        let rate_exc = 1.0 / self.tau_exc;
206        let rate_inh = 1.0 / self.tau_inh;
207        let decay_v = (-self.dt / self.tau_v).exp();
208        let decay_exc = (-self.dt / self.tau_exc).exp();
209        let decay_inh = (-self.dt / self.tau_inh).exp();
210        let (a_exc_next, i_exc_next) = Self::filter_candidates(
211            self.a_exc,
212            self.i_exc,
213            exc_current,
214            self.tau_exc,
215            self.dt,
216            decay_exc,
217        )?;
218        let (a_inh_next, i_inh_next) = Self::filter_candidates(
219            self.a_inh,
220            self.i_inh,
221            inh_current,
222            self.tau_inh,
223            self.dt,
224            decay_inh,
225        )?;
226        let exc_steady = self.tau_exc * exc_current;
227        let inh_steady = self.tau_inh * inh_current;
228        let v_steady = self.v_rest + exc_steady - inh_steady;
229        let v_next = v_steady
230            + (self.v - v_steady) * decay_v
231            + Self::drive_contribution(
232                self.i_exc - exc_steady,
233                self.a_exc - exc_steady,
234                self.tau_exc,
235                self.dt,
236                [rate_v, rate_exc],
237                [decay_v, decay_exc],
238            )?
239            - Self::drive_contribution(
240                self.i_inh - inh_steady,
241                self.a_inh - inh_steady,
242                self.tau_inh,
243                self.dt,
244                [rate_v, rate_inh],
245                [decay_v, decay_inh],
246            )?;
247        if !v_next.is_finite() {
248            return Err(AlphaError::NonFiniteCandidate);
249        }
250        self.a_exc = a_exc_next;
251        self.i_exc = i_exc_next;
252        self.a_inh = a_inh_next;
253        self.i_inh = i_inh_next;
254        if v_next >= self.v_threshold {
255            self.v = self.v_rest;
256            return Ok(1);
257        }
258        self.v = v_next;
259        Ok(0)
260    }
261
262    /// Preserve the legacy scalar API while failing closed on invalid input.
263    #[inline]
264    pub fn step(&mut self, exc_current: f64, inh_current: f64) -> i32 {
265        self.try_step(exc_current, inh_current).unwrap_or(0)
266    }
267
268    /// Restore the documented rest state while preserving configuration.
269    pub fn reset(&mut self) {
270        self.v = self.v_rest;
271        self.a_exc = 0.0;
272        self.i_exc = 0.0;
273        self.a_inh = 0.0;
274        self.i_inh = 0.0;
275    }
276}
277
278impl Default for AlphaNeuron {
279    fn default() -> Self {
280        Self::new()
281    }
282}
283
284/// Complete post-update trace and final-state receipt.
285pub struct AlphaTrace {
286    /// Membrane-potential trace.
287    pub v: Vec<f64>,
288    /// Excitatory rise-state trace.
289    pub a_exc: Vec<f64>,
290    /// Excitatory current trace.
291    pub i_exc: Vec<f64>,
292    /// Inhibitory rise-state trace.
293    pub a_inh: Vec<f64>,
294    /// Inhibitory current trace.
295    pub i_inh: Vec<f64>,
296    /// Binary spike-event trace represented as floating point for ABI parity.
297    pub spikes: Vec<f64>,
298    /// Final ``[v, a_exc, i_exc, a_inh, i_inh]`` state, including the initial
299    /// state for an empty batch.
300    pub final_state: [f64; 5],
301    /// Number of candidate-crossing spikes.
302    pub spike_count: usize,
303}
304
305/// Simulate a complete caller-owned piecewise-constant drive batch atomically.
306#[allow(clippy::too_many_arguments)]
307pub fn simulate(
308    v: f64,
309    a_exc: f64,
310    i_exc: f64,
311    a_inh: f64,
312    i_inh: f64,
313    v_rest: f64,
314    v_threshold: f64,
315    tau_v: f64,
316    tau_exc: f64,
317    tau_inh: f64,
318    dt: f64,
319    exc_current: &[f64],
320    inh_current: &[f64],
321) -> Result<AlphaTrace, AlphaError> {
322    if exc_current.len() != inh_current.len() {
323        return Err(AlphaError::NonFiniteInput);
324    }
325    let mut neuron = AlphaNeuron::with_parameters(
326        v,
327        a_exc,
328        i_exc,
329        a_inh,
330        i_inh,
331        v_rest,
332        v_threshold,
333        tau_v,
334        tau_exc,
335        tau_inh,
336        dt,
337    )?;
338    if !exc_current.iter().all(|value| value.is_finite())
339        || !inh_current.iter().all(|value| value.is_finite())
340    {
341        return Err(AlphaError::NonFiniteInput);
342    }
343    let steps = exc_current.len();
344    let mut v_trace = Vec::with_capacity(steps);
345    let mut a_exc_trace = Vec::with_capacity(steps);
346    let mut i_exc_trace = Vec::with_capacity(steps);
347    let mut a_inh_trace = Vec::with_capacity(steps);
348    let mut i_inh_trace = Vec::with_capacity(steps);
349    let mut spikes = Vec::with_capacity(steps);
350    let mut spike_count = 0usize;
351    for index in 0..steps {
352        let spike = neuron.try_step(exc_current[index], inh_current[index])?;
353        spike_count += spike as usize;
354        v_trace.push(neuron.v);
355        a_exc_trace.push(neuron.a_exc);
356        i_exc_trace.push(neuron.i_exc);
357        a_inh_trace.push(neuron.a_inh);
358        i_inh_trace.push(neuron.i_inh);
359        spikes.push(f64::from(spike));
360    }
361    Ok(AlphaTrace {
362        v: v_trace,
363        a_exc: a_exc_trace,
364        i_exc: i_exc_trace,
365        a_inh: a_inh_trace,
366        i_inh: i_inh_trace,
367        spikes,
368        final_state: [
369            neuron.v,
370            neuron.a_exc,
371            neuron.i_exc,
372            neuron.a_inh,
373            neuron.i_inh,
374        ],
375        spike_count,
376    })
377}
378
379#[cfg(test)]
380mod tests {
381    use super::*;
382
383    #[test]
384    fn defaults_match_catalogue_model_family() {
385        let neuron = AlphaNeuron::new();
386        assert_eq!(
387            (
388                neuron.v,
389                neuron.a_exc,
390                neuron.i_exc,
391                neuron.a_inh,
392                neuron.i_inh,
393                neuron.v_rest,
394                neuron.v_threshold,
395                neuron.tau_v,
396                neuron.tau_exc,
397                neuron.tau_inh,
398                neuron.dt,
399            ),
400            (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 20.0, 5.0, 10.0, 1.0)
401        );
402    }
403
404    #[test]
405    fn filter_matches_exact_alpha_cascade() {
406        let decay = (-0.5_f64 / 5.0).exp();
407        let (rise_next, current_next) =
408            AlphaNeuron::filter_candidates(0.25, 0.1, 2.0, 5.0, 0.5, decay).unwrap();
409        let steady = 5.0 * 2.0;
410        let expected_rise = steady + (0.25 - steady) * decay;
411        let expected_current = steady + decay * ((0.1 - steady) + (0.25 - steady) * 0.5 / 5.0);
412        assert!((rise_next - expected_rise).abs() < 1.0e-12);
413        assert!((current_next - expected_current).abs() < 1.0e-12);
414    }
415
416    #[test]
417    fn drive_contribution_handles_equal_time_constants() {
418        let rate = 1.0 / 20.0;
419        let decay = (-0.5_f64 / 20.0).exp();
420        let exact =
421            AlphaNeuron::drive_contribution(0.3, 0.2, 20.0, 0.5, [rate, rate], [decay, decay])
422                .unwrap();
423        let expected = rate * decay * (0.3 * 0.5 + 0.2 * 0.5 * 0.5 / (2.0 * 20.0));
424        assert!((exact - expected).abs() < 1.0e-12);
425    }
426
427    #[test]
428    fn spike_resets_only_the_membrane() {
429        let mut neuron = AlphaNeuron {
430            v: 0.9,
431            a_exc: 0.4,
432            i_exc: 0.6,
433            a_inh: 0.2,
434            i_inh: 0.1,
435            v_threshold: 0.5,
436            ..AlphaNeuron::new()
437        };
438        let (a_exc_before, i_exc_before, a_inh_before, i_inh_before) =
439            (neuron.a_exc, neuron.i_exc, neuron.a_inh, neuron.i_inh);
440        assert_eq!(neuron.try_step(0.0, 0.0), Ok(1));
441        assert_eq!(neuron.v, 0.0);
442        let decay_exc = (-1.0_f64 / 5.0).exp();
443        let decay_inh = (-1.0_f64 / 10.0).exp();
444        assert!((neuron.a_exc - a_exc_before * decay_exc).abs() < 1.0e-12);
445        assert!(
446            (neuron.i_exc - decay_exc * (i_exc_before + a_exc_before * 1.0 / 5.0)).abs() < 1.0e-12
447        );
448        assert!((neuron.a_inh - a_inh_before * decay_inh).abs() < 1.0e-12);
449        assert!(
450            (neuron.i_inh - decay_inh * (i_inh_before + a_inh_before * 1.0 / 10.0)).abs() < 1.0e-12
451        );
452    }
453
454    #[test]
455    fn invalid_step_is_atomic() {
456        let mut neuron = AlphaNeuron::new();
457        neuron.v = 0.25;
458        neuron.a_exc = 0.5;
459        let before = (
460            neuron.v,
461            neuron.a_exc,
462            neuron.i_exc,
463            neuron.a_inh,
464            neuron.i_inh,
465        );
466        assert_eq!(
467            neuron.try_step(f64::NAN, 0.0),
468            Err(AlphaError::NonFiniteInput)
469        );
470        assert_eq!(
471            (
472                neuron.v,
473                neuron.a_exc,
474                neuron.i_exc,
475                neuron.a_inh,
476                neuron.i_inh
477            ),
478            before
479        );
480    }
481
482    #[test]
483    fn invalid_configuration_is_rejected() {
484        assert!(matches!(
485            AlphaNeuron::with_parameters(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.5, 20.0, 5.0, 10.0, 1.0),
486            Err(AlphaError::InvalidScale)
487        ));
488        assert!(matches!(
489            AlphaNeuron::with_parameters(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 5.0, 10.0, 1.0),
490            Err(AlphaError::InvalidScale)
491        ));
492    }
493
494    #[test]
495    fn batch_matches_scalar_and_empty_preserves_initial_state() {
496        let empty = simulate(
497            0.1,
498            0.2,
499            0.3,
500            0.4,
501            0.5,
502            0.0,
503            1.0,
504            20.0,
505            5.0,
506            10.0,
507            1.0,
508            &[],
509            &[],
510        )
511        .unwrap();
512        assert!(empty.v.is_empty() && empty.spikes.is_empty());
513        assert_eq!(empty.final_state, [0.1, 0.2, 0.3, 0.4, 0.5]);
514
515        let exc = [0.5, 0.6, 0.7, 0.8];
516        let inh = [0.1, 0.2, 0.1, 0.2];
517        let batch = simulate(
518            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 20.0, 5.0, 10.0, 1.0, &exc, &inh,
519        )
520        .unwrap();
521        let mut scalar = AlphaNeuron::new();
522        let mut count = 0;
523        for index in 0..exc.len() {
524            count += scalar.try_step(exc[index], inh[index]).unwrap() as usize;
525        }
526        assert_eq!(
527            batch.final_state,
528            [
529                scalar.v,
530                scalar.a_exc,
531                scalar.i_exc,
532                scalar.a_inh,
533                scalar.i_inh
534            ]
535        );
536        assert_eq!(batch.spike_count, count);
537    }
538
539    #[test]
540    fn reset_restores_documented_rest_state_not_configuration() {
541        let mut neuron = AlphaNeuron::new();
542        neuron.v = 0.4;
543        neuron.a_exc = 0.3;
544        neuron.i_inh = 0.2;
545        neuron.reset();
546        assert_eq!(
547            (
548                neuron.v,
549                neuron.a_exc,
550                neuron.i_exc,
551                neuron.a_inh,
552                neuron.i_inh
553            ),
554            (0.0, 0.0, 0.0, 0.0, 0.0)
555        );
556        assert_eq!(
557            (
558                neuron.v_threshold,
559                neuron.tau_v,
560                neuron.tau_exc,
561                neuron.tau_inh,
562                neuron.dt
563            ),
564            (1.0, 20.0, 5.0, 10.0, 1.0)
565        );
566    }
567}
sc_neurocore_engine/neurons/simple_spiking/alpha.rs

sc_neurocore_engine/neurons/simple_spiking/
alpha.rs
