sc_neurocore_engine/
ei_network.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 — Floating-point E-I balanced network for Studio
8
9//! Fused E-I balanced LIF network simulation in f64 arithmetic.
10//!
11//! Runs a complete excitatory-inhibitory network with CSR weight matrix,
12//! Poisson external drive, and spike-scatter coupling in a single Rust
13//! call — no per-step Python overhead.
14
15use rand::{RngExt, SeedableRng};
16use rand_xoshiro::Xoshiro256PlusPlus;
17
18/// Result of an E-I network simulation.
19pub struct EIResult {
20    pub spike_times: Vec<f64>,
21    pub spike_neurons: Vec<u32>,
22    pub n_exc: u32,
23    pub n_inh: u32,
24    pub rate_time: Vec<f64>,
25    pub exc_rates: Vec<f64>,
26    pub inh_rates: Vec<f64>,
27    pub mean_exc_rate: f64,
28    pub mean_inh_rate: f64,
29}
30
31/// Run a complete E-I balanced LIF network simulation.
32///
33/// All computation happens in Rust — connectivity build, Poisson input,
34/// Euler integration, spike detection, rate binning.
35#[allow(clippy::too_many_arguments)]
36pub fn simulate_ei(
37    n_exc: usize,
38    n_inh: usize,
39    w_ee: f64,
40    w_ei: f64,
41    w_ie: f64,
42    w_ii: f64,
43    p_conn: f64,
44    ext_rate: f64,
45    duration: f64,
46    dt: f64,
47    seed: u64,
48) -> EIResult {
49    let n = n_exc + n_inh;
50    let n_steps = ((duration / dt) as usize).min(50_000);
51    let mut rng = Xoshiro256PlusPlus::seed_from_u64(seed);
52
53    // LIF params
54    let tau_m = 20.0_f64;
55    let v_rest = -65.0_f64;
56    let v_threshold = -50.0_f64;
57    let v_reset = -65.0_f64;
58    let tau_ref = 2.0_f64;
59
60    // Build CSR weight matrix
61    let mut row_offsets = Vec::with_capacity(n + 1);
62    let mut col_indices = Vec::new();
63    let mut values = Vec::new();
64    row_offsets.push(0usize);
65
66    for i in 0..n {
67        let i_exc = i < n_exc;
68        for j in 0..n {
69            if i == j {
70                continue;
71            }
72            if rng.random::<f64>() >= p_conn {
73                continue;
74            }
75            let j_exc = j < n_exc;
76            let w = match (i_exc, j_exc) {
77                (true, true) => w_ee,
78                (true, false) => w_ie,
79                (false, true) => -w_ei,
80                (false, false) => -w_ii,
81            };
82            col_indices.push(j);
83            values.push(w);
84        }
85        row_offsets.push(col_indices.len());
86    }
87
88    // State arrays
89    let mut v = vec![v_rest; n];
90    let mut refractory = vec![0.0_f64; n];
91    let mut prev_spiked = vec![false; n];
92
93    // Recording
94    let mut spike_times: Vec<f64> = Vec::new();
95    let mut spike_neurons: Vec<u32> = Vec::new();
96    let bin_size = (n_steps / 100).max(1);
97    let n_bins = n_steps / bin_size;
98    let mut exc_rates = vec![0.0_f64; n_bins];
99    let mut inh_rates = vec![0.0_f64; n_bins];
100    let mut exc_bin = 0u32;
101    let mut inh_bin = 0u32;
102
103    let ext_lambda = (ext_rate * dt / 1000.0).max(0.0);
104    let exp_neg_lambda = (-ext_lambda).exp();
105
106    for t in 0..n_steps {
107        // Decay refractory
108        for r in refractory.iter_mut() {
109            *r = (*r - dt).max(0.0);
110        }
111
112        // Synaptic input from previous spikes (CSR scatter)
113        let mut syn = vec![0.0_f64; n];
114        for i in 0..n {
115            if !prev_spiked[i] {
116                continue;
117            }
118            let start = row_offsets[i];
119            let end = row_offsets[i + 1];
120            for k in start..end {
121                syn[col_indices[k]] += values[k];
122            }
123        }
124
125        // External Poisson + Euler step
126        for i in 0..n {
127            if refractory[i] > 0.0 {
128                continue;
129            }
130            // Knuth Poisson sampling (fast for small lambda)
131            let mut k = 0u32;
132            let mut p = 1.0_f64;
133            loop {
134                p *= rng.random::<f64>();
135                if p <= exp_neg_lambda {
136                    break;
137                }
138                k += 1;
139            }
140            let ext = k as f64 * 5.0;
141            let dv = (-(v[i] - v_rest) / tau_m + ext + syn[i]) * dt;
142            v[i] += dv;
143        }
144
145        // Spike detection
146        prev_spiked.fill(false);
147        for i in 0..n {
148            if refractory[i] > 0.0 {
149                continue;
150            }
151            if v[i] >= v_threshold {
152                v[i] = v_reset;
153                refractory[i] = tau_ref;
154                prev_spiked[i] = true;
155                spike_times.push(t as f64 * dt);
156                spike_neurons.push(i as u32);
157                if i < n_exc {
158                    exc_bin += 1;
159                } else {
160                    inh_bin += 1;
161                }
162            }
163        }
164
165        // Rate binning
166        if (t + 1) % bin_size == 0 {
167            let bi = t / bin_size;
168            if bi < n_bins {
169                let bin_t = bin_size as f64 * dt / 1000.0;
170                exc_rates[bi] = exc_bin as f64 / n_exc.max(1) as f64 / bin_t.max(0.001);
171                inh_rates[bi] = inh_bin as f64 / n_inh.max(1) as f64 / bin_t.max(0.001);
172            }
173            exc_bin = 0;
174            inh_bin = 0;
175        }
176    }
177
178    let rate_time: Vec<f64> = (0..n_bins)
179        .map(|i| i as f64 * bin_size as f64 * dt)
180        .collect();
181
182    let mean_exc = if exc_rates.iter().any(|&r| r > 0.0) {
183        let pos: Vec<f64> = exc_rates.iter().copied().filter(|&r| r > 0.0).collect();
184        pos.iter().sum::<f64>() / pos.len() as f64
185    } else {
186        0.0
187    };
188    let mean_inh = if inh_rates.iter().any(|&r| r > 0.0) {
189        let pos: Vec<f64> = inh_rates.iter().copied().filter(|&r| r > 0.0).collect();
190        pos.iter().sum::<f64>() / pos.len() as f64
191    } else {
192        0.0
193    };
194
195    EIResult {
196        spike_times,
197        spike_neurons,
198        n_exc: n_exc as u32,
199        n_inh: n_inh as u32,
200        rate_time,
201        exc_rates,
202        inh_rates,
203        mean_exc_rate: (mean_exc * 10.0).round() / 10.0,
204        mean_inh_rate: (mean_inh * 10.0).round() / 10.0,
205    }
206}
207
208#[cfg(test)]
209mod tests {
210    use super::*;
211
212    #[test]
213    fn ei_network_runs_without_panic() {
214        let r = simulate_ei(20, 5, 0.1, 0.4, 0.1, 0.4, 0.2, 10.0, 50.0, 0.1, 42);
215        assert_eq!(r.n_exc, 20);
216        assert_eq!(r.n_inh, 5);
217        assert!(!r.rate_time.is_empty());
218        assert_eq!(r.exc_rates.len(), r.rate_time.len());
219    }
220
221    #[test]
222    fn ei_network_with_high_drive_produces_spikes() {
223        // ext_rate=5000 Hz → lambda=0.5 per step → frequent Poisson events
224        let r = simulate_ei(40, 10, 0.1, 0.4, 0.1, 0.4, 0.2, 5000.0, 100.0, 0.1, 42);
225        assert!(
226            !r.spike_times.is_empty(),
227            "high drive should produce spikes"
228        );
229    }
230}
sc_neurocore_engine/ei_network.rs

sc_neurocore_engine/
ei_network.rs
