sc_neurocore_engine/analysis/
rate.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 — Instantaneous and population firing rate estimation
8
9/// Build a 1-D kernel for rate estimation.
10fn build_kernel(kernel: &str, sigma_steps: usize) -> Vec<f64> {
11    match kernel {
12        "gaussian" => {
13            let hw = 3 * sigma_steps;
14            let mut k: Vec<f64> = (0..2 * hw + 1)
15                .map(|i| {
16                    let x = i as f64 - hw as f64;
17                    (-0.5 * (x / sigma_steps as f64).powi(2)).exp()
18                })
19                .collect();
20            let s: f64 = k.iter().sum();
21            k.iter_mut().for_each(|v| *v /= s);
22            k
23        }
24        "exponential" => {
25            let hw = 5 * sigma_steps;
26            let mut k: Vec<f64> = (0..hw)
27                .map(|i| (-(i as f64) / sigma_steps as f64).exp())
28                .collect();
29            let s: f64 = k.iter().sum();
30            k.iter_mut().for_each(|v| *v /= s);
31            k
32        }
33        "rectangular" => {
34            let hw = sigma_steps;
35            let len = 2 * hw + 1;
36            vec![1.0 / len as f64; len]
37        }
38        _ => vec![1.0],
39    }
40}
41
42/// Convolve signal with kernel (mode="same").
43fn convolve_same(signal: &[f64], kernel: &[f64]) -> Vec<f64> {
44    let n = signal.len();
45    let kn = kernel.len();
46    let pad = kn / 2;
47    let mut out = vec![0.0; n];
48    for i in 0..n {
49        let mut acc = 0.0;
50        for j in 0..kn {
51            let si = i as isize + j as isize - pad as isize;
52            if si >= 0 && (si as usize) < n {
53                acc += signal[si as usize] * kernel[j];
54            }
55        }
56        out[i] = acc;
57    }
58    out
59}
60
61/// Instantaneous firing rate via kernel convolution (Hz).
62pub fn instantaneous_rate(binary_train: &[f64], dt: f64, kernel: &str, sigma_ms: f64) -> Vec<f64> {
63    let sigma_steps = (sigma_ms / (dt * 1000.0)).round().max(1.0) as usize;
64    let mut k = build_kernel(kernel, sigma_steps);
65    // Normalise to produce Hz: divide by dt
66    let s: f64 = k.iter().sum();
67    if s > 0.0 {
68        let norm = s * dt;
69        k.iter_mut().for_each(|v| *v /= norm);
70    }
71    // Re-normalise because build_kernel already normalised by sum
72    // The Python does: k /= k.sum() * dt, applied to un-normalised kernel.
73    // So we need to undo the build_kernel normalisation and redo correctly.
74    let mut k_raw = build_kernel(kernel, sigma_steps);
75    let s_raw: f64 = k_raw.iter().sum();
76    let norm = s_raw * dt;
77    k_raw.iter_mut().for_each(|v| *v /= norm);
78
79    convolve_same(binary_train, &k_raw)
80}
81
82/// Population-level instantaneous firing rate (Hz).
83pub fn population_rate(trains: &[Vec<f64>], dt: f64, sigma_ms: f64) -> Vec<f64> {
84    if trains.is_empty() {
85        return vec![];
86    }
87    let min_len = trains.iter().map(|t| t.len()).min().unwrap_or(0);
88    if min_len == 0 {
89        return vec![];
90    }
91    let mut total = vec![0.0; min_len];
92    for t in trains {
93        for (i, &v) in t.iter().take(min_len).enumerate() {
94            total[i] += v;
95        }
96    }
97    instantaneous_rate(&total, dt, "gaussian", sigma_ms)
98}
99
100/// Peri-stimulus time histogram.
101/// Returns (rates_hz, bin_centers_ms).
102pub fn psth(trials: &[Vec<f64>], bin_ms: f64, dt: f64) -> (Vec<f64>, Vec<f64>) {
103    if trials.is_empty() {
104        return (vec![], vec![]);
105    }
106    let max_len = trials.iter().map(|t| t.len()).max().unwrap_or(0);
107    let bin_steps = (bin_ms / (dt * 1000.0)).round().max(1.0) as usize;
108    let n_bins = max_len / bin_steps;
109    if n_bins == 0 {
110        return (vec![], vec![]);
111    }
112    let mut counts = vec![0.0; n_bins];
113    for trial in trials {
114        let usable = trial.len().min(n_bins * bin_steps);
115        let n_full = usable / bin_steps;
116        for b in 0..n_full {
117            let s: f64 = trial[b * bin_steps..(b + 1) * bin_steps].iter().sum();
118            counts[b] += s;
119        }
120    }
121    let n_trials = trials.len() as f64;
122    let bin_sec = bin_ms / 1000.0;
123    let rates: Vec<f64> = counts.iter().map(|&c| c / (n_trials * bin_sec)).collect();
124    let centers: Vec<f64> = (0..n_bins).map(|i| (i as f64 + 0.5) * bin_ms).collect();
125    (rates, centers)
126}
127
128#[cfg(test)]
129mod tests {
130    use super::*;
131
132    #[test]
133    fn test_instantaneous_rate_gaussian() {
134        let mut train = vec![0.0; 100];
135        train[50] = 1.0;
136        let rate = instantaneous_rate(&train, 0.001, "gaussian", 5.0);
137        assert_eq!(rate.len(), 100);
138        // Peak should be near index 50
139        let peak_idx = rate
140            .iter()
141            .enumerate()
142            .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
143            .unwrap()
144            .0;
145        assert!((peak_idx as i32 - 50).abs() <= 1);
146        assert!(rate[50] > 0.0);
147    }
148
149    #[test]
150    fn test_instantaneous_rate_rectangular() {
151        let train = vec![1.0; 10];
152        let rate = instantaneous_rate(&train, 0.001, "rectangular", 1.0);
153        assert_eq!(rate.len(), 10);
154        // Middle values should be ~1000 Hz (1 spike per ms)
155        assert!(rate[5] > 500.0);
156    }
157
158    #[test]
159    fn test_population_rate() {
160        let t1 = vec![1.0, 0.0, 1.0, 0.0, 1.0];
161        let t2 = vec![0.0, 1.0, 0.0, 1.0, 0.0];
162        let rate = population_rate(&[t1, t2], 0.001, 1.0);
163        assert_eq!(rate.len(), 5);
164        // Combined: all ones → high rate
165        assert!(rate[2] > 0.0);
166    }
167
168    #[test]
169    fn test_psth() {
170        let trial = vec![0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0];
171        let (rates, centers) = psth(&[trial.clone(), trial], 5.0, 0.001);
172        assert_eq!(rates.len(), 2);
173        assert_eq!(centers.len(), 2);
174        assert!((centers[0] - 2.5).abs() < 0.01);
175        assert!(rates[0] > 0.0);
176    }
177
178    #[test]
179    fn test_psth_empty() {
180        let (r, c) = psth(&[], 10.0, 0.001);
181        assert!(r.is_empty());
182        assert!(c.is_empty());
183    }
184}
sc_neurocore_engine/analysis/rate.rs

sc_neurocore_engine/analysis/
rate.rs
